Genome-Wide Identification of Metal Tolerance Protein Genes in Populus trichocarpa and Their Roles in Response to Various Heavy Metal Stresses

Metal tolerance proteins (MTPs) are plant divalent cation transporters that play important roles in plant metal tolerance and homeostasis. Poplar is an ideal candidate for the phytoremediation of heavy metals because of its numerous beneficial attributes. However, the definitive phylogeny and heavy metal transport mechanisms of the MTP family in poplar remain unknown. Here, 22 MTP genes in P. trichocarpa were identified and classified into three major clusters and seven groups according to phylogenetic relationships. An evolutionary analysis suggested that PtrMTP genes had undergone gene expansion through tandem or segmental duplication events. Moreover, all PtrMTPs were predicted to localize in the vacuole and/or cell membrane, and contained typical structural features of the MTP family, cation efflux domain. The temporal and spatial expression pattern analysis results indicated the involvement of PtrMTP genes in poplar developmental control. Under heavy metal stress, most of PtrMTP genes were induced by at least two metal ions in roots, stems or leaves. In addition, PtrMTP8.1, PtrMTP9 and PtrMTP10.4 displayed the ability of Mn transport in yeast cells, and PtrMTP6 could transport Co, Fe and Mn. These findings will provide an important foundation to elucidate the biological functions of PtrMTP genes, and especially their role in regulating heavy metal tolerance in poplar.


Introduction
Heavy metal pollution is becoming a more and more serious problem globally. It poses a grave risk for food and environmental safety, as well as human health [1]. For example, Cd and Ni are carcinogenic, teratogenic and mutagenic, while Cu, Hg and Pb can cause neurological diseases [2]. Some plants have hyperaccumulation and/or hypertolerance characteristics for specific heavy metals, and are excellent candidates for phytoremediation to reduce environmental heavy metal pollutant levels. These plants must possess a distinct system that contributes to the high capacity of heavy metal enrichment and tolerance, including metal uptake, efflux, chelation, translocation, intracellular sequestration and storage; among these, the metal transporters play a crucial role [3,4].
Cation diffusion facilitators (CDFs) are divalent cation (Zn 2+ , Co 2+ , Fe 2+ , Cd 2+ , Ni 2+ and Mn 2+ ) transporters which have important functions in maintenance of metal homeostasis. Since the identification of Cupriavidus metallidurans in 1995, increasing numbers of CDF genes have been cloned, functionally investigated and classified into three clusters: Zn-CDF, Fe/Zn-CDF and Mn-CDF [5,6].  To further analyze the phylogenetic relationship of PtrMTP proteins to their counterparts from er plants, a total of 118 MTP sequences from eight plant species were used to construct ylogentic tree. The 22 PtrMTP proteins could be categorized to seven groups (1, 5, 6, 7, 8, 9 an ). Among them, group 9 is the largest one, containing 7 PtrMTP members; group 1 includes rMTP members, while groups 5, 6, 7 and 12 contain only one PtrMTP each. Amazingly, group ntains 5 tandem repeat PtrMTP members from PtrMTP8.2 to PtrMTP8.6 and a single PtrMTP8 igures 2). Further, the seven PtrMTP groups are classified into Zn-CDFs, Zn/Fe-CDFs an n-CDFs clusters. To further analyze the phylogenetic relationship of PtrMTP proteins to their counterparts from other plants, a total of 118 MTP sequences from eight plant species were used to construct a phylogentic tree. The 22 PtrMTP proteins could be categorized to seven groups (1, 5, 6, 7, 8, 9 and 12). Among them, group 9 is the largest one, containing 7 PtrMTP members; group 1 includes 5 PtrMTP members, while groups 5, 6, 7 and 12 contain only one PtrMTP each. Amazingly, group 8 contains 5 tandem repeat PtrMTP members from PtrMTP8.2 to PtrMTP8.6 and a single PtrMTP8.1 ( Figure 2). Further, the seven PtrMTP groups are classified into Zn-CDFs, Zn/Fe-CDFs and Mn-CDFs clusters.

Figure 2.
Phylogenetic relationships of MTP proteins in P. trichocarpae and other plant species. One hundred and eighteen MTP proteins are clustered into three major substrate-specific groups and seven primary groups, which are highlight in different colors. The different symbols represent the MTP proteins of different species as follows.

Structure and Characteristics Analysis of PtrMTP Genes
The genome annotation files of poplar were applied to the TBtools software for an exon-intron organizations analysis of PtrMTP genes. As shown in Figure 3a-b, although the number of introns in the PtrMTP genes of the three clusters ranged from 0 to 12, most of the related members that clustered closely shared similar introns in terms of number and phase. Of the three clusters, Zn-CDFs contained the smallest number of introns (group 1 contained only one, and group 12 contained none), except for group 5, which harbored 9 introns; Mn-CDFs contained 3-6 introns (group 8 contained 3, 5 or 6, and group 9 all contained 5), whereas Zn/Fe-CDFs contained the largest number of introns (group 6 contained 10, group 7 contained 12) (Figure 3a-b). Additionally, all the PtrMTP genes contained phase 0 and phase 2 introns, while only the members of group 5-7 contained phase 1 intron (Figure 3a-b).
Next, the physicochemical parameters of the 22 PtrMTPs were analyzed further. The length of the coding sequence (CDS) of PtrMTP genes ranged from 498 bp (PtrMTP8.6) to 2610 bp (PtrMTP12), One hundred and eighteen MTP proteins are clustered into three major substrate-specific groups and seven primary groups, which are highlight in different colors. The different symbols represent the MTP proteins of different species as follows. Solid triangles: Arabidopsis thaliana; hollow triangles: Brachypodium diastychon; reverse hollow triangles: Zea mays; solid diamonds: Cucumis sativus; hollow diamonds: Vitis vinifera; solid circles: Populus trichocarpae; hollow circles: Nicotiana tabacum; solid squares: Sorghum bicolor; hollow squares: Oryza sativa.

Structure and Characteristics Analysis of PtrMTP Genes
The genome annotation files of poplar were applied to the TBtools software for an exon-intron organizations analysis of PtrMTP genes. As shown in Figure 3a-b, although the number of introns in the PtrMTP genes of the three clusters ranged from 0 to 12, most of the related members that clustered closely shared similar introns in terms of number and phase. Of the three clusters, Zn-CDFs contained the smallest number of introns (group 1 contained only one, and group 12 contained none), except for group 5, which harbored 9 introns; Mn-CDFs contained 3-6 introns (group 8 contained 3, 5 or 6, and group 9 all contained 5), whereas Zn/Fe-CDFs contained the largest number of introns (group 6 contained 10, group 7 contained 12) (Figure 3a-b). Additionally, all the PtrMTP genes contained phase 0 and phase 2 introns, while only the members of group 5-7 contained phase 1 intron (Figure 3a   Next, the physicochemical parameters of the 22 PtrMTPs were analyzed further. The length of the coding sequence (CDS) of PtrMTP genes ranged from 498 bp (PtrMTP8.6) to 2610 bp (PtrMTP12), with 165 to 869 amino acids, as well as a relative molecular weight (MW) ranging from 18.375 to 97.498 KDa ( Table 1). The total average of hydropathicity (GRAVY) of the PtrMTP proteins ranged from -0.235 (PtrMTP10.1) to 0.329 (PtrMTP4) ( Table 1). Moreover, PtrMTP7 has the highest isoelectric point (pI), i.e., 7.24, whereas those of other PtrMTPs were below 7 (Table 1). Furthermore, all PtrMTP proteins were expected to localize in the vacuole, and notably, some (PtrMTP9, PtrMTP10.2, PtrMTP10.3 and PtrMTP10.4) were also localized in the cytoplasm membrane (Table 1). In addition, most PtrMTP proteins contained 4-6 typical TMDs, whereas PtrMTP11.2 had only three, and PtrMTP12 harbored 12; none was found in PtrMTP6 and PtrMTP8.6 proteins ( Table 1).

Chromosomal Localization and Gene Duplication Analysis of PtrMTP Genes
To explore the physical locations of the PtrMTP genes, genome annotation files were downloaded from the phytozome12 database and analyzed using the TBtools software. The results showed that the 21 out of the 22 PtrMTP genes were located on nine poplar chromosomes with an uneven distribution pattern ( Figure 4). Most PtrMTP genes were assigned to chromosomes 01 and 10, which contained 7 and 6 PtrMTP genes, respectively. Interestingly, some PtrMTP genes were closely located to one another in a chromosome, such as five PtrMTP genes (PtrMTP8. To better understand the selection type of these duplication genes, the ratios of the number of nonsynonymous substitutions per nonsynonymous site (Ka) to the number of synonymous substitutions per synonymous site (Ks) of the 10 gene pairs mentioned above were further calculated. As shown in Table 2, the Ka/Ks ratios of all duplicated pairs of the PtrMTP gene were less than 1, which suggested that all these duplication events were under negative selection, based on the summaries from Hurst [21].

Conserved Motif and Domain Architectures Analysis of PtrMTP Proteins
Our study found that PtrMTP proteins contained a total of twelve conserved motifs (Figure 3c), among which only three encode functional domains according to the annotation from the Pfam or InterProScan tools (Figure 3c and Supplementary Table S4). Motifs 1 and 7 encode Cation efflux, while motif 2 encodes ZT_dimer. Motif 6 was widely shared by all PtrMTPs, except for PtrMTP5, PtrMTP7 and PtrMTP8.6 ( Figure 3c). Motifs 7, 11 and 12 were mainly distributed in the Zn-CDFs cluster, while motifs 1, 2, 3, 4, 5, 8, 9 and 10 were specifically distributed in the Mn-CDFs cluster. It was also found that the members of the same cluster or group had similar motif types and distributions (Figure 3c). Of the three clusters, Zn/Fe-CDFs contained the smallest number of motifs (group 6 only contained two, and group 7 contained none), Zn-CDFs contained 2-6 motifs (group 1 contained 4 or 5, group 5 contained 2, and group 12 contained 6), whereas Mn-CDFs had the largest number and the most similar types (group 8 contained 8 or 9, and group 9 contained 9 or 10), except for PtrMTP8.6, which contained only three (Figure 3c).
A conserved domain analysis showed that all the PtrMTP proteins contained the cation efflux domain ( Figure 5), a typical feature of the MTP protein [5], whereas the members of groups 6, 8 (except for PtrMTP8.6) and 9 possessed a ZT dimer, an important zinc transporter dimerization domain.

Cis-Acting Elements in the Promoter Regions of PtrMTP Genes
A total of 1271 cis-acting regulatory elements were identified, which were classified into nine major classes, i.e., 917 elements for gene transcription, 52 elements for abiotic stress, 1 element for biotic stress, 8 elements for tissue expression, 5 elements for secondary metabolism, 80 elements for phytohormonal response, 168 elements for light response, 5 elements for circadian control and 35 elements for site binding (Table 3 and Supplementary Table S5).
Among these, gene transcription elements including 366 CAAT elements and 551 TATA-box elements, which were the most abundance elements, and light responsiveness elements, such as ACE, ATCT-motif, Box 4 and CATT-motif, were commonly present in all PtrMTPs. Responsive elements of various phytohormones, such as ABRE, P-box, GARE-motif, TATC-box and SARE, were found in all PtrMTP genes, except for the PtrMTP1.1 and PtrMTP1.2 genes. Abiotic stress elements including LTR, MBS, TC-rich repeat, WUN-motif, ARE and GC-motif were distributed in the promoters of all PtrMTP genes except for PtrMTP1. 1. In comparison, the AT-rich element involved in biotic stress responsive was detected only in the promoter of the PtrMTP3.2 gene. Additionally, tissue expression elements including CAT-box and GCN4_motif were present in promoters of PtrMTP1. 2

Cis-Acting Elements in the Promoter Regions of PtrMTP Genes
A total of 1271 cis-acting regulatory elements were identified, which were classified into nine major classes, i.e., 917 elements for gene transcription, 52 elements for abiotic stress, 1 element for biotic stress, 8 elements for tissue expression, 5 elements for secondary metabolism, 80 elements for phytohormonal response, 168 elements for light response, 5 elements for circadian control and 35 elements for site binding (Table 3 and Supplementary Table S5).
Among these, gene transcription elements including 366 CAAT elements and 551 TATA-box elements, which were the most abundance elements, and light responsiveness elements, such as ACE, ATCT-motif, Box 4 and CATT-motif, were commonly present in all PtrMTPs. Responsive elements of various phytohormones, such as ABRE, P-box, GARE-motif, TATC-box and SARE, were found in all PtrMTP genes, except for the PtrMTP1.1 and PtrMTP1.2 genes. Abiotic stress elements including LTR, MBS, TC-rich repeat, WUN-motif, ARE and GC-motif were distributed in the promoters of all PtrMTP genes except for PtrMTP1.1. In comparison, the AT-rich element involved in biotic stress responsive was detected only in the promoter of the PtrMTP3.2 gene. Additionally, tissue expression elements including CAT-box and GCN4_motif were present in promoters of PtrMTP1.2, PtrMTP3.2, PtrMTP4, PtrMTP6, PtrMTP9, PtrMTP10.1, PtrMTP10.3 and PtrMTP10.4 genes. Moreover, secondary metabolism elements were only detected in the promoters of PtrMTP1.1, PtrMTP3.2, PtrMTP6 and PtrMTP12, including MBSI involved in flavonoid metabolism and O2-site zein involved in zein metabolism. Notably, site-binding elements were found in all PtrMTP genes, except for PtrMTP1.1, PtrMTP7, PtrMTP8.5, PtrMTP6, PtrMTP10.4 and PtrMTP12, whereas circadian control elements were only present in promoters of PtrMTP8.4, PtrMTP8.5 and PtrMTP8.6 ( Table 3 and Supplementary  Table S5). These results indicated a diverse and complicated control of PtrMTP gene expression at the transcriptional level.

The Temporal and Spatial Expression Patterns of PtrMTP Genes
The tissue expression patterns of PtrMTPs were investigated by using transcriptome data. As shown in Figure 6, all 22 PtrMTP genes were expressed in the 12 tested tissues (log2(FPKM+1) > 0), except for PtrMTP8.6 (which had weak expression only in late dormant bud, root and male catkin) and

Expression Profiles of PtrMTPs under Different Heavy Metal Treatments
To gain more insight into the gene expression regulatory mechanism of PtrMTPs, four-week-old tested tube plantlets of P. trichocarpa were subjected to seven different metal treatments. The relative expression levels of PtrMTPs in roots, stems and leaves were investigated.
Under normal conditions, the expression levels of PtrMTP4, PtrMTP8.3, PtrMTP8.4, PtrMTP8.5, PtrMTP10.2 and PtrMTP10.4 were higher in roots, whereas those of the PtrMTP1.1, PtrMTP7, PtrMTP9, PtrMTP10.1, PtrMTP10.3, PtrMTP11.1 and PtrMTP12 genes displayed higher expression levels in stems, and PtrMTP3.1, PtrMTP3.2, PtrMTP11.2 genes displayed higher expression levels in leaves. However, the PtrMTP1.2, PtrMTP5 and PtrMTP8.2 genes showed similar expression levels in roots and stems, which were higher than those in the leaves. PtrMTP6, PtrMTP8.1 and PtrMTP8.6 have similar expression levels in stems and leaves, which were higher than those in roots (Figure 7).  We present an overview of the expression levels of all the PtrMTP genes under heavy metal toxicity relative to these under normal conditions in Table 5. In detail, we summarized the PtrMTP genes in each tissue with expression changes over four times: In root, Cd enhanced the expression of PtrMTP11.1; Cu increased the expression levels of PtrMTP8.1 and PtrMTP10.3, but decreased the expression levels of PtrMTP9; Mn repressed the expression levels of PtrMTP9 and PtrMTP10.3; Ni also repressed the expression levels of PtrMTP10.3, but Zn enhanced its expression. In stem, Cd repressed the expression levels of PtrMTP12; Co increased the expression levels of PtrMTP8.6 but decreased the expression levels of  (Figure 7 and Table 5). However, the expression levels of the PtrMTP3.1, PtrMTP3.2 and PtrMTP6 genes nearly did not change in each tissue under heavy metal toxicity ( Figure 7 and Table 5).

Evolution and Differentiation of PtrMTP Genes as well as Their Proteins Architectures
Here, a total of 22 PtrMTPs were identified and named by bioinformatics methods. Compared with that of previous report, three new MTP genesn were found in P. trichocarpa [20]. The number of MTPs in P. trichocarpa was only second to N. tabacum among the plant species in which the MTPs have been identified [5,8,[22][23][24][25][26][27][28]. However, despite the large number of PtrMTPs, the homolog of AtMTP2 in P. trichocarpa was not detected. These results indicated that gene expansion and/or gene loss may have occurred in the history of PtrMTP gene family evolution. This hypothesis was later supported by the gene duplication analysis, which revealed that ten duplication events existed in the PtrMTP gene family, among which five were segmental and another five were tandem duplications ( Figure 4, Table 2). Gene duplication in the genome may lead to the production of new genes with novel functions [29]. Thus, it would be of great interest to investigate the relationship between MTP gene family differentiation and heavy metal tolerance in different plant species in future studies.
As the most typical structural features of MTP proteins, the cation efflux domain and the modified signature sequence were detected in all of the PtrMTPs, although some other motifs/domains were not present in certain PtrMTP members. Interestingly, PtrMTP6 and PtrMTP8.6 did not possess any TMDs (Table 1), a common structure of membrane proteins [30], which suggested that these two proteins might play novel roles, i.e., other than transporters. In addition, it is notable that unlike other PtrMTP members, PtrMTP12 had the largest protein size (869 amino acids), MW (97.50 kD) and TMD number (12 TMDs) ( Table 1). This result was consistent with the characteristics of other plants MTP12 [23,25,26], indicating the distinctive biological functions and evolutionary processes of MTP12. Moreover, ZT_dimer has been recognized as the dimerization region of metal ion transporters, through which the homodimers or heterodimers of MTPs can form [31]. In this study, the ZT_dimer was detected in members of groups 6, 8 (except for PtrMTP8.6) and 9, but the question of whether this domain is correlated with their functions requires further investigation ( Figure 5).

Regulation of PtrMTP Gene Expression in P. trichocarpae
Gene expression control happens at two different levels: one is transcriptional regulation, and the other is post-transcriptional regulation. For the former, cis-acting regulatory elements (CRE) play essential roles by interacting with RNA polymerase and specific transcription factors. In this study, the CAAT-box and TATA-box, which are involved in regulating the expression frequency and initiation of transcription respectively [32], were detected in the upstream region of PtrMTP genes at a high frequency (Table 3 and Supplementary Table S5). In addition, light-responsive elements, phytohormonal responsive elements and abiotic stress and site-binding elements were also widely distributed in most of PtrMTP genes (Table 3 and Supplementary Table S5), implying that PtrMTP genes could be transcriptionally regulated by multiple stimuli.
Previous studies have shown that miRNAs play versatile roles in plant growth and development control and stresses responses [33]. Most of the 11 miRNAs identified in this study perform their functions by cleaving target mRNA, like ptc-miR2111a/b, ptc-miR172b-5p, ptc-miR172g-5p, ptc-miR6427-3p, ptc-miR6464 and ptc-miR6466-3p, while others, e.g., ptc-miR6426a, ptc-miR6426b, ptc-miR473b and ptc-miR480, carry out their functions by translation inhibition. Most of these miroRNAs have been demonstrated to play important roles against environmental stress responses. miR473 was reported only in tree species, and participates in the response to mycorrhizal symbiosis and drought in Poncirus trifoliate and populus, respectively [34][35][36]. miR2111 could be induced by Pi starvation in Arabidopsis, but it could also fulfill shoot-to-root translocation to control rhizobial infection in legume roots [37,38]. Moreover, the miR6426 are involved in the response to nutrient deficiencies and contribute to Mg-deficiency tolerance in Citrus sinensis [39]. miR172b is a key controller of the autotrophic development transition in Arabidopsis [40]. These findings imply the possible involvement of PtrMTP in abiotic and biotic stress response through post-transcriptional regulation mediated by miRNAs.

The Diverse Expression Patterns of PtrMTP Genes
Of 22 PtrMTPs, 15 exhibited organ/tissue specificity of gene expression during growth and development of poplar, and most PtrMTP genes could respond to at least two metal ions in poplar roots, stems or leaves under heavy metal stresses. The abundant cis-acting regulatory elements and different expression pattern of the PtrMTP gene imply that PtrMTP members play an important role in plant development and stress responses.
Notably, some paralogous genes from the same group, especially some tandem and segmental duplication gene pairs, showed different expression patterns when analyzed by transcriptome data and quantitative RealTime-PCR (qRT-PCR) under various metal ions stresses in different tissues and developmental stages ( Figure 6, Figure 7), such as PtrMPT1.1/PtrMPT1.2 and PtrMPT3.1 /PtrMPT3.2 gene pairs. This may due to the long-term evolution of PtrMPT genes, and resulted in the precise and specific regulation mechanism of metal homeostasis in poplar species. Notably, although the majority of tissue expression patterns of PtrMPT genes, including groups 1, 3, 4, 5, 6, 7 and 8, were consistent with the result of transcriptome data, some genes belonging to groups 9, 10, 11 and 12 showed some contradictions between the qRT-PCR results and transcriptome data (Figure 7). This inconsistency was also mentioned in our previous literature [26].
Furthermore, the expression of most PtrMTP genes could be induced by multiple metals, although some were not putative heavy metal substrates (Figure 7, Table 5). For instance, excess Cu, which is not a potential substrate for the MTP family, could induce the transcription levels of PtrMTP genes ( Figure 7, Table 5). Similar findings were also made in sweet orange and turnip [24,25]. Other than Mn, the expression of the Mn-CDFs, PtrMTP9 and PtrMTP10.3 were sharply upregulated by excess Fe and Zn, respectively ( Figure 7, Table 5). On the other hand, MTP genes may not response to their transport substrates. For example, PtrMTP1.2, PtrMTP3.1, PtrMTP3.2, PtrMTP4 and PtrMTP12, which are Zn-CDF genes, underwent no large changes under excess Zn treatment. We speculate that this may be due to the Zn treatment concentrations used in our experiment. In other words, the response of these genes may require higher or lower concentrations of Zn treatment. Certainly, gene response may also happen at the post-transcriptional level, rather than the transcription level. For example, the abundance of the Zn transporter CsMTP1 was increased four-fold in excess Zn treated cucumber roots, whereas the mRNA level of the CsMTP1 gene was not significantly changed [41]. Other studies have shown that the response of heavy metal transport happened at the post-translational level through changing protein levels, localization and the turnover of transporters [42][43][44][45]. In addition, some MTPs even displayed no response to Zn treatment at both the transcriptional and post-transcriptional levels, though it functioned in maintaining intracellular Zn homeostasis by Zn transport, such as AtMTP1 in Arabidopsis [9]. Interestingly, although the accumulation of AtMTP12 was not dependent on Zn concentration, but it could form a functional complex with another MTP, AtMTP5, to transport Zn [12]. These results indicate a complicated and multilayered regulatory mechanism underlying the response of MTPs to heavy metal substrates. Thus, it is necessary to investigate protein level changes under excess heavy metal treatment, and to identify the protein complex of PtrMTPs in future studies.

Some PtrMTPs were Co, Fe and Mn Transporters in Yeast Cells
A yeast-metal sensitivity test assay was performed to clarify the heavy metal substrates of selected PtrMTP transporters. Our results indicated that PtrMTP6, an MTP member with no predictable TMD domain, could transport three different heavy metals: Co, Fe and Mn ions (Figure 8). Similar findings were recently reported by Migocka et al., who found that CsMTP6 could affect iron and manganese homeostasis in cucumber mitochondria [46]. In addition, three PtrMTP members, PtrMTP8.1, PtrMTP9, and PtrMTP10.4, showed specific transport abilities for Mn in yeast cells (Figure 8). These results were in agreement with previous studies. CsMTP8 and OsMTP8.1 are tonoplast-localized Mn transporters which are responsible for Mn tolerance in cucumber and rice, respectively [22,47,48]. Cucumber MTP9 homologue was proved to be a plasma membrane antiporter that functions in Mn 2+ and Cd 2+ efflux from root cells [49]. Erbasol et al. found that the Golgi apparatus localized sea beet BmMTP10 was specific to Mn 2+ transport, with a role in reducing excess cellular Mn 2+ levels in yeast [50]. Nevertheless, PtrMTP8.4, unlike its paralog PtrMTP8.1, did not show any metal transport abilities in the present study, indicating a functional diversity of PtrMTP within the same group (Figure 8). CsMTP4 from cucumber is the only MTP4 homolog that has been functionally characterized to date. This protein is localized in the vacuolar membranes to sequestrate Zn and Cd into vacuole [42]. However, in our investigation, PtrMTP4 could not transport Zn and Cd as well as other heavy metals in yeast cells (Figure 8). These results suggest that homologs of MTP across plant species may have different biological functions.

Identification and Phylogenetic Analysis of the MTPs in P. Trichocarpa
Twelve MTP genes in A. thaliana were downloaded from the TAIR database (https://www. arabidopsis.org/), and were used as queries to perform in TBLASTP searches against the P. trichocarpa genome [51]. After analysis with InterProScan (http://www.ebi.ac.uk/interpro/search/sequence-search), the remaining 22 nonredundant candidates were recognized as PtrMTP proteins. Sequence similarity analysis was performed at NCBI (https://blast.ncbi.nlm.nih.gov/Blast.cgi) as described previously [26]. For phylogenetic analysis, sequences of MTP proteins from P. trichocarpa and other plant species were first aligned by the Clustal X2.1 [52]; then, phylogenetic trees were established with the Maximum Likelihood method by the MEGA6.06 software [53]. Sequences of MTPs from B. diastychon, C. sativus, N. tabacum, O. sativa, Sorghum bicolor, V. vinifera, and Z. mays were obtained from the databases, as described by Liu et al. [26]. The TBtools software (https://github.com/CJ-Chen/TBtools) was used to determine the exon-intron organization and chromosomal localization of PtrMTP genes using the genome annotation files downloaded from phytozome12 database (https://phytozome.jgi.doe.gov/pz/ portal.html) [54]. The gene duplication events of PtrMTP genes were determined by Multiple Collinearity Scan toolkit (MCScanX) [55]. DnaSP v6 software was used to compute the Ka and Ks substitution rates, as described by Rozas et al. [56].

Prediction of Cis-Acting Regulatory Elements and MicroRNA (miRNA) Target Sites of PtrMTP Genes
Promoter sequences located 1.0 kb upstream of PtrMTP genes were searched using the phytozome12 database, and related cis-acting regulatory elements were analyzed using the PlantCare database (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/). The miRNA target sites of PtrMTPs were predicted using a small RNA target analysis server (psRNATarget server: http://plantgrn.noble.org/psRNATarget/).

Tissue Expression Pattern Based on RNA-seq Data
The fragments per kilobase of exon model per million mapped reads (FPKM) of PtrMTP genes in 12 different tissues in varying developmental stages were downloaded from the Phytozome12 database. Subsequently, the normalized data (log2(FPKM+1)) were used to estimate the expression levels, and a heatmap was created by TBtools software.

Growth Condictions and Heavy Metal Treatments
Tube plantlets of P. trichocarpa were grown in a greenhouse with 16 h/8 h photoperiod (24 • C/18 • C). Thirty-day-old plants were then placed in 1/2 Hoagland solutions (pH 6.0) supplemented with different heavy metal, with concentrations of 0.1 mM CdCl 2 , 0.1 mM CoCl 2 , 0.1 mM CuSO 4 , 0.5 mM FeSO 4 , 1 mM MnSO 4 , 0.1 mM NiSO 4 and 0.5 mM ZnSO 4 , respectively, and normal 1/2 Hoagland solutions were used as the control (CK). Then, 24 h later, the roots, stems and leaves of tube plantlets were collected, and used as materials for RNA extraction.

RNA Extraction and qRT-PCR Analysis
RNA extraction, cDNA preparation and qRT-PCR were performed as suggested by Liu et al. with minor modifications [26]. Each experiment was performed with three technical replicates. Two house-keeping genes, UBQ (GenBank accession LOC7455401) and EF1α (GenBank accession LOC18100225), were used as internal reference. The primers used for qRT-PCR are listed in Supplementary Table S1. The qRT-PCR programs were as described by Gao et al. [57]. The relative expression values were calculated using the 2 -∆∆ Ct method [58].

Yeast Transformation and Growth Assay
The full coding regions of six PtrMTP genes were amplified from the cDNAs that was reverse-transcribed from the total RNA of plant leaves. The specific primers used for yeast expression vector construction are presented in Supplementary Table S2. The PCR products were then inserted into Kpn I + Xba I or Hind III + Xba I sites of the pYES2.0 vector.

Conclusions
In this study, 22 MTP members in P. trichocarpa were identified, and a systematic and comprehensive analysis of PtrMTP genes was performed. The 22 PtrMTPs were divided into three major substrate-specific clusters and seven groups. The MTP gene family in poplar underwent expansions in MTP1, MTP3, MTP8, MTP10, and MTP11 compared with those of Arabidopsis, although MTP2 might have been lost to evolutionary history. The structural characteristics of PtrMTP were similar within groups, but were diverse among different groups. The temporal and spatial expression patterns of PtrMTP genes were either similar or varied within the same group. In response to different heavy metal stresses, most PtrMTP genes were induced by at least two metal ions in roots, stems or leaves. Additionally, PtrMTP8.1, PtrMTP9 and PtrMTP10.4 were found to function as Mn transporters, and PtrMTP6 could transport three different heavy metal ions, i.e., Co, Fe and Mn in yeast cells, indicating that these proteins might play important roles in heavy metal homeostasis, detoxification and tolerance in poplar. These results will provide an important foundation for better understanding the mechanism of heavy metal transport mediated by PtrMTP proteins. In addition, our study will also provide important gene resources for the genetic modification of the heavy metal accumulation abilities of plants which can be widely used in phytoremediation.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.