Genome data for P. trichocarpa were downloaded from Phytozome v12.1 (https://phytozome.jgi.doe.gov/pz/portal.html) (Goodstein et al., 2012). The typical MFS protein domains (PF07690, PF16983, PF05631, PF07672, PF05977, PF06779, PF12832, and PF13347) were obtained from Pfam (http://pfam.xfam.org/) (Finn et al., 2014). Scanning the poplar genome for potential PtrMFS proteins was conducted with HMMER3 (Jaina et al., 2013). Verification was performed using the SMART database (http://smart.embl-heidelberg.de/) (Ivica and Peer, 2018) and Pfam database to remove proteins without MFS domains.

MFS protein sequences of P. trichocarpa and A. thaliana were downloaded from Phytozome. We used MEGA v5.1 (Tamura et al., 2011) with the Maximum-Likelihood (ML) method to construct a phylogenetic tree of MFS proteins using the JTT (protein mutation data matrix) + G (site-specific variations in mutation rate) + F (mutation frequency data) model.

Physical and chemical properties of PtrMFS proteins were predicted with ProtParam (https://web.expasy.org/protparam/) (Gasteiger et al., 2003), including the number of amino acids, molecular weight, theoretical isoelectric point (pI), aliphatic index, grand average of hydrophilicity (GRAVY), chemical formulas, total number of atoms, and instability index.

To analyze gene structures of PtrMFSs, we downloaded DNA and coding sequences of PtrMFS from Phytozome database. DNA and coding sequences for each PtrMFS gene were aligned to obtain gene structures. TBtools (Chen et al., 2018) was used to visualize gene structures of PtrMFSs, and MEME (http://meme-suite.org/) (Bailey et al., 2006) was used to identify conserved motifs in PtrMFS proteins. Annotations of identified motifs were obtained from InterProScan (https://www.ebi.ac.uk/interpro/search/sequence/) (Quevillon et al., 2005).

Secondary structures of PtrMFS proteins were predicted by SOPMA (https://npsa-prabi.ibcp.fr/cgi-bin/npsa_automat.pl?page=npsa_sopma.html) (Geourjon and Delé age, 1995), and tertiary structures were constructed by SWISS-MODEL (https://swissmodel.expasy.org/) (Torsten et al., 2003).

Topological heterogeneity models of PtrMFS proteins were predicted with Protter (http://wlab.ethz.ch/protter/start/) (Ulrich et al., 2014). Subcellular localization was predicted using WoLF PSORT (https://wolfpsort.hgc.jp) (Paul et al., 2007).

For each PtrMFS gene, sequences starting at 2,000 bp upstream of the start codon were downloaded from Phytozome. Cis-acting elements were extracted using PlantCARE (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/) (Rombauts et al., 1999). TBtools was used to visualize cis-acting elements.

PtrMFS genes were mapped to the genome of P. trichocarpa, and the chromosomal distribution of PtrMFSs in poplar was visualized with TBtools. TBtools and MCScanX (Wang et al., 2012) were used to analyze tandem duplication events in the PtrMFS gene family, and Dual Synteny Plotter (Tang et al., 2008) was used to analyze segmental duplication events and collinearity between PtrMFSs and homologous gene pairs from other species (A. thaliana, Eucalyptus grandis, Oryza sativa, and Solanum lycopersicum). TBtools was used to visualize the results. The ratio of non-synonymous substitution to synonymous substitution (Ka/Ks) of duplicate gene pairs was determined with KaKs_Calculator (Zhang et al., 2006).

WT Pdpap seedlings used in this research were cultured on 0.5X Murashige and Skoog (MS) medium supplemented with 0.01 mg/ml 1-naphthaleneacetic acid (NAA). Two-month-old WT Pdpap plants at the same growth stage was selected. Fifty milliliters of F. oxysporum at a spore density of 1 × 10⁵/ml was poured onto Pdpap roots (Kang et al., 2018). Infection times were 6, 12, 24, or 48 h. WT Pdpap treated with 50 ml ddH₂O was used as control. Every treatment consisted of 3 duplicated samples. The materials involved in this research and their preparation methods have been reported in detail in our previous research (Diao et al., 2021).

RNA-Seq was used to explore gene expression patterns of the MFS gene family from PdPap under F. oxysporum stress. The data were described in detail in our previous research (Diao et al., 2021). Differentially expressed genes (DEGs) were identified with DESeq2 using log₂ (fold change) ≥ and adjusted p-value (padj) ≤ 0.05 as the criteria. Expression patterns of PdPapMFS genes under different treatments of F. oxysporum at 6, 12, 24, and 48 h were visualized with heatmaps. We identified upregulated and downregulated DEGs in response to F. oxysporum and displayed the data by Venn diagrams. To validate the RNA-Seq data, we further analyzed expression levels of DEGs under F. oxysporum stress by qRT-PCR. The qRT-PCR was performed on a Stratagene Mx3000P real-time PCR system (Agilent Technologies, Santa Clara, CA, United States) using the 2 × SYBR Green qPCR Master Mix kit (Bimake, Shanghai, China). Reaction systems are shown in Supplementary Table S1 and primer sequences designed by Primer5.0 (Zhai et al., 2008) are shown in Supplementary Table S2. qRT-PCR amplification conditions were as follows: initial denaturation at 94°C for 30 s; 44 cycles of 94°C for 12 s, 58°C for 30 s, and 72°C for 45 s, then 79°C for 1 s. The reaction specificity was determined by performing a melting-curve analysis from 55°C to 99°C, with fluorescence readings taken every 0.5°C for 1 s. The amplification curve was obtained after analyzing the raw data, and the cycle threshold (Ct) was set with a fluorescence threshold of 0.01 (Wang et al., 2012). Relative expression level of target genes was calculated by the 2^−△△Ct method (Kenneth and Thomas, 2002). Three duplicates were set for each gene. PdpapActin and PdPapEF1-α were used as the internal control genes (Huang et al., 2008).

We analyzed co-expression-based gene networks using STRING (Damian et al., 2011), and visualized the results with Cytoscape (Ideker, 2011). Co-expression analysis was performed on the 41 PdPapMFS genes identified from RNA-Seq data analysis as described above. Genes were annotated with gene ontology (GO)-based functions using agriGO v2.0 (http://systemsbiology.cau.edu.cn/agriGOv2/index.php) (Du et al., 2010).

Data were analyzed with the Statistical Software Package for Social Science (SPSS) version 17.0 (Kumar, 2014). Using Student’s t-test to compare the data, p < 0.05 was considered statistically significant (Choi et al., 2012). Significant differences (p < 0.05) are indicated by different lowercase letters in figures.

We identified 41 MFS genes in P. trichocarpa (each named PtrMFS with a number based on position in the poplar genome; Supplementary Table S3). To determine evolutionary relationships of genes in this family, we constructed a ML phylogenetic tree (Figure 1) using protein sequences from poplar and Arabidopsis. As shown in Figure 1, we classified PtrMFSs into six groups of various sizes. Group MFS_4 is the largest with nine genes, while group MFS_5 is the smallest with five genes.

The physicochemical properties of the identified PtrMFSs vary significantly (Supplementary Table S3). The average length of PtrMFS proteins is 498 amino acids (range = 372–698 aa) and the average molecular weight is 54.27 kDa (range = 40.20–78.23 kDa). The theoretical pI of PtrMFSs is in the range of 5.54–10.07 and the aliphatic index is in the range of 85.96–118.46. The hydrophilicity value of PtrMFSs ranges from 0.183 to 0.758. All of the proteins are predicted to be hydrophilic. PtrMFSs consist of five elements: C, H, N, O, and S. The total number of atoms in each protein range from 5,724 to 11,110. The instability indexes of PtrMFSs range from 24.80 to 52.06; an index >40 indicates an unstable protein, while <40 indicates a stable protein.

PtrMFS genes within the same groups share similar structures with respect to introns and exons (Figure 2). Four members in group MFS_2 contain 14 exons. Three members in group MFS_3 contain 2 exons and only 1 intron. Three members in group MFS_6 contain 2 exons and another three members contain 3 exons. Members with closer relationships share more similar gene structures and exon lengths.

PtrMFS proteins with similar evolutionary relationships have the same or similar conserved domains. We identified 20 conserved motifs using MEME, and motif annotations were predicted with Pfam and InterProScan (Supplementary Table S4). We found that motifs 2, 3, 4, 5, 6, 7, 8, 12, 13, 14, 17, and 19 were annotated as MFS transporters. Motifs 2, 3, 11, and 16 were high affinity nitrate transporter-related; motif 8 was tetracycline resistance domain. Motif 13 was annotated as protein zinc induced facilitator. Motif 15 was annotated as anion transporter and solute carrier family. Motif 1 was annotated as LytB protein and motif 20 was annotated as SPX domain. Motifs 9, 10, and 18 were unknown and could not be annotated. Results show that motifs 13 and 20 are present only in all members of group MFS_1 (Figure 3). Motif 8 is present only in all members of group MFS_2. Motif 6, 7, 12, 17, and 18 are present only in members of group MFS_3. Motif 4, 5, 14, 15, and 19 are present only in all members of group MFS_4. Motif 1, 2, 5, 7, 9, 11, and 16 only occurred in members of group MFS_6. These results demonstrate that there are many group-specific motifs, which may be correlated to specific biological functions.

In the secondary structure of PtrMFS proteins (Supplementary Table S5), an average of 234.59 aa belonged to alpha helixes, accounting for 46.93% of total secondary structures. On average, a further 74.73 aa (15.08%) belonged to extended strands, 23.29 aa (4.71%) belonged to beta turns, and 165.37 aa (33.28%) belonged to random coils.

Tertiary structures were generated primarily using 6s4m.1.A and 6e9c.1.A as templates (Supplementary Table S6). They have functions of major facilitator superfamily domain-containing proteins. Sequence identity of the 41 PtrMFS proteins identified in this study and their corresponding templates is 20.80% (11.02–37.71%) on average.

Topological heterogeneity models of PtrMFSs proteins (Supplementary Table S7) showed that all members of the PtrMFS family have transmembrane helical segments. PtrMFS19, PtrMFS23, and PtrMFS26 have no N-glycosylation sites. PtrMFS8, PtrMFS14, PtrMFS18, PtrMFS19, and PtrMFS28 have signal peptides.

Subcellular localization predictions (Supplementary Table S8) showed that a majority of PtrMFSs are located at the plasma membrane. PtrMFS26 is predicted to be located at the endoplasmic reticulum, PtrMFS31 at the vacuole, and PtrMFS4, PtrMFS5, PtrMFS13, PtrMFS22, PtrMFS29, PtrMFS39, and PtrMFS41 at the chloroplast.

We used PlantCARE (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/) to predict cis-acting elements within 2,000 bp upstream of PtrMFS genes (Supplementary Table S9; Figure 4). Many predicted cis-acting elements identified are involved in hormone responses, such as response to auxin, response to gibberellin, response to salicylic acid, and response to abscisic acid. Many elements are also predicted to be involved in defense and stress responses, meristem expression, root specific expression, flavonoid biosynthetic gene regulation, and wounding response.

The location of PtrMFS genes within chromosomes is variable (Figure 5); they are unevenly distributed in 16 of the 19 poplar chromosomes. The number of genes located in each chromosome is irrelevant to the size of chromosomes. Chromosome 9 contains the most PtrMFS genes, with 8. Only one PtrMFS gene each was located on chromosome 4, chromosome 5, chromosome 10, and chromosome 12, and no PtrMFS genes were found on chromosome 11, chromosome 13, chromosome 17, or chromosome 19.

We analyzed within-genome duplication events of PtrMFS genes using MCScanX. PtrMFS19 and PtrMFS20 (on chromosome 7), PtrMFS24 and PtrMFS25 (on chromosome 9), PtrMFS36 and PtrMFS37 (on chromosome 15) were categorized as derived from tandem duplication events. Fourteen genes exhibit 8 pairs of segmental duplication events (Figure 6; Supplementary Table S10). They are unevenly distributed in eight of 19 chromosomes. Segmental duplication events may be the important driving force for diversity of PtrMFS genes. The ratio of Ka/Ks is an important indicator of selective pressure in evolution, with a Ka/Ks < 1 indicating negative selection. Ka/Ks values of PtrMFS duplicate genes range from 0.16 to 0.37, with an average value of 0.24, suggesting that PtrMFS genes have been subject to purifying selection during evolution.

To explore the DNA sequence similarity between PtrMFS genes and related genes from other representative species, we constructed collinearity relationship maps of Populus trichocarpa with three dicotyledons (E. grandis, S. lycopersicum, and A. thaliana) and one monocotyledon (O. sativa) (Figure 7). We identified 32 repetitive events in E. grandis, 26 in S. lycopersicum, 23 in A. thaliana, and seven in O. sativa (Supplementary Table S11). Collinearity blocks were mainly distributed in the first 10 chromosomes of P. trichocarpa. There was greater collinearity between PtrMFS genes and those in dicotyledons than those in the monocotyledon.

We analyzed expression patterns of PdPapMFS genes in response to F. oxysporum with different lengths of time post-inoculation using RNA-Seq data (Supplementary Table S12). Statistical results indicated that PdPapMFSs expression are all in different fold-changes (Figure 8). The relative gene expression values of PdPapMFS are shown in Supplementary Table S13.

In the samples collected 6 h after inoculation with F. oxysporum, we identified 40 DEGs, including 13 upregulated genes, and 27 downregulated genes. PdPapMFS38 was the most upregulated (3.14X) and PdPapMFS25 was the most downregulated (−2.07X). In the 12 h treatment group, 41 DEGs were found, including 14 upregulated genes and 27 downregulated genes. PdPapMFS38 was the most upregulated (3.22X) and PdPapMFS17 was the most downregulated (−1.43X). In the 24 h treatment group, 40 DEGs were found, including 14 upregulated genes, and 26 downregulated genes. PdPapMFS38 was the most upregulated (2.93X) and PdPapMFS17 was the most downregulated (−1.16X). In the 48 h treatment group, 40 DEGs were found, including 23 upregulated genes, and 17 downregulated genes. PdPapMFS38 was the most upregulated (2.32X) and PdPapMFS17 was the most downregulated (−1.14X). DEG expression in the four treatment groups was visualized with heatmaps (Figure 9) and Venn diagrams (Figure 10).

To validate the RNA-Seq data, we performed qRT-PCR to analyze expression levels of putative DEGs in response to F. oxysporum stress. Results of RNA-Seq and qRT-PCR were generally consistent, with a few exceptions such as PdPapMFS8, PdPapMFS16, PdPapMFS34, which may be caused by experimental errors of RNA-Seq or qRT-PCR (Figure 11).

We constructed a co-expression network using RNA-Seq data of the 41 PtrMFS genes (Figure 12). A large proportion of the genes are shared in the networks, such as PtrMFS32 (LOC7458066), PtrMFS38 (LOC7455882), PtrMFS33 (LOC7496936), PtrMFS6 (LOC7470780), PtrMFS5 (LOC7478542), PtrMFS39 (LOC7466045), PtrMFS3 (LOC18094568), PtrMFS13 (LOC7453664), PtrMFS25 (LOC7463307). We have observed that the shared gene expression patterns in the co-expression network are not exactly the same, which indicated that MFS genes under the stress of F. oxysporum may have complex regulatory characteristics.

Using agriGO, we performed gene set enrichment analysis with PdPapMFS genes (Figure 13). Based on the annotations, these genes can participate in a variety of biological processes (Figure 14A; Supplementary Table S14), such as in process terms “transmembrane transport,” “establishment of localization,” and “cellular process.” Many PdPapMFS genes are enriched in various molecular functions (Figure 14B; Supplementary Table S14), such as in function terms “transporter activity,” and “transmembrane transporter activity.” Regarding cellular localization (Figure 14C; Supplementary Table S14), PdPapMFS genes are enriched in localization terms “integral to membrane,” “intrinsic to membrane,” and “cell part.” It is worth mentioning that, many MFS genes are enriched in “transmembrane transport,” “transport,” “transporter activity,” “transmembrane transporter activity,” “membrane,” “integral to membrane,” and “intrinsic to membrane.” That also proved our speculation that MFS genes are related to the transmembrane transport function of substances.