Based on the full-length transcriptome assembly results, the CDS sequence, protein sequence, and annotation file of the genome were obtained (bioproject_accession: PRJNA1028662). The Pfam numbers PF16916 and PF01545 of MTP were obtained from the InterPro (https://www.ebi.ac.uk/interpro/ database to download its Hidden Markov Model (HMM) file. The HMM file of MTP was used as input, and hmmsearch (v3.3.2) software (Johnson et al., 2010) was used to search the tulip protein sequence file with E-value < 0.05.

Secondly, the sequence files of MTP proteins in rice and A. thaliana were downloaded from Phytozome (https://phytozome-next.jgi.doe.gov/) database as submission files. Blastp (v2.13.0) (Kelley et al., 2004) was used to search for potential MTP proteins in tulip protein files, and the results of hmmsearch and blastp were combined to obtain candidate gene ID. The protein sequences of candidate genes were extracted and submitted to Pfam, CDD (https://www.ncbi.nlm.nih.gov/Structure/cdd/cdd.shtml) and SMART (http://smart.embl-heidelberg.de) database for domain confirmation to obtain the determined TgMTP protein.

By comparing with A. thaliana and rice MTP family members and naming TgMTP family members, the physicochemical properties of TgMTP family members were analyzed by Expasy (https://web.expasy.org/protparam/) software. Cell-PLoc-2 (http://www.csbio.sjtu.edu.cn/bioinf/Cell-PLoc-2/) software was used to predict the subcellular localization of TgMTP.

All MTP protein sequences of tulips were used as input files, and the conserved motifs of TgMTP family were analyzed by MEME software (version: v5.0.5, http://meme-suite.org/) (Bailey et al., 2006). The parameter motif search number is 10, the minimum motif amino acid length is 6, and the maximum motif amino acid length is 100. The R package gggenomes (https://thackl.github.io/gggenomes/) was used to visually analyze the conserved motif and domain of TgMTP, and TBtools was used to visualize the CDS and UTR regions. The signal peptide of TgMTP protein was predicted by SignalP (v5.0) software (https://services.healthtech.dtu.dk/service.php?SignalP-5.0).

Ka/Ks analysis, a prevalent and crucial tool in bioinformatics, plays a significant role in examining the evolutionary dynamics of nucleic acid molecules. Within the realm of genetics, the Ka/Ks ratio signifies the balance between the nonsynonymous substitution rate (Ka) and the synonymous substitution rate (Ks) in two protein-coding genes. This metric serves as an indicator of selective pressure acting on these genes, using KaKs _ Calculator (v2.0) (Wang et al., 2010) software for analysis.

The MTP family protein sequences of A. thaliana, rice, tomato and soybean were downloaded from the Phytozome database, and the phylogenetic tree was drawn using MEGA (v10) software (Kumar et al., 2018). Clustal W (Thompson et al., 1994) was selected as the comparison method, and the phylogenetic tree was constructed by Neighbor-joining method. The model was p-distance, the missing data method was Partial deletion, the cutoff was 50%, and Bootstrap was set to 1000. We use Evolview3 (https://www.evolgenius.info/evolview/) to beautify the tree.

The protein secondary structure of TgMTP was predicted by sopma (https://npsa-prabi.ibcp.fr/cgi-bin/npsa_automat.pl?page=/NPSA/npsa_sopma.html) website, and the protein tertiary structure of TgMTP was predicted by SwissModel (https://swissmodel.expasy.org/). The clusterProfiler (Yu et al., 2012) in R was used to analyze the GO (Gene Ontology) annotation of TgMTP.

The psRNATarget online tool (https://www.zhaolab.org/psRNATarget/) was use to analysis the miRNA target gene network of TgMTP gene family. Additionally, since miRNA research in tulip remains scarce, we selected the model plant A. thaliana, which has the most extensive miRNA databases, to analyze the miRNAs regulating TgMTP genes. Cytoscape (Smoot et al., 2011) was used to visualize the interaction networks.

Tulip cultivar “Dow Jones” was exposed at 5°C for 4 weeks and then plant them in the soil (growing media 4.0, Shangdao Co., Ltd, Jinan, China). After treatment, untreated tulip bulbs were used as CK group (0 d), and tulip bulbs were washed and placed in the medium containing different concentrations of metal ions (for FeSO₄: 1 mM; for MnSO₄: 5 mM; for ZnSO₄: 1 mM and for CdCl₂: 5 mM) in Hoagland nutrient solution. They were then cultivated in a conventional water-fertilizer incubator in soil with a 25°C/21°C day-night temperature and 16 h day^-1 photoperiod. Bulbs subjected to these conditions were sampled at 1, 3, 5, and 7 days and total RNA was extracted using the BioRun Plant Total RNA Extraction Kit (Biorun Co., Ltd, Wuhan, China). cDNA was synthesized, followed by qRT-PCR to measure relative expression. Relative gene expression was calculated using the 2^-ΔΔCt method. Graphpad prism 9 software (https://www.graphpad.com/features) was used for plotting and statistical analysis. Primers were designed at NCBI and are listed in Supplementary Table 1 .

The heterogeneic expression of TgMTP7.1 in yeast BY4741 and ycf1Δ were used to investigate whether TgMTP7.1 can improve the Cd tolerance of yeast. The ORF of TgMTP7.1 was amplified from tulip cDNA and inserted into the yeast expression vector pYES2 under the control of a galactose-inducible promoter. The resulting pYES2-TgMTP7.1 and empty vector pYES2 were transformed into yeast cells. Transgenic yeast strains were cultured in synthetic complete (SC) medium lacking uracil with 2% (w/v) galactose to induce gene expression. For Cd tolerance assays, the exponentially growing yeast cultures were normalized to an OD₆₀₀ of 0.2, and serial dilutions (10⁰, 10¹, 10², 10³) were prepared (Liu et al., 2019) 2 μL of each dilution was spotted onto SC-Ura agar medium containing 5 mM CdCl₂. Plates were incubated at 30°C for 3-5 days before comparing growth differences. To ascertain metal tolerance, the growth of yeast (as measured by the OD₆₀₀ value) was quantified at various intervals, employing methodologies delineated by Peng (Peng et al., 2018).

Based on homologous sequence alignment and conserved domain analysis, a total of 11 TgMTP gene family members were identified in this study, named as TgMTP1, TgMTP2, TgMTP3, TgMTP4, TgMTP5, TgMTP6.1, TgMTP6.2, TgMTP6.3, TgMTP7.1, TgMTP7.2 and TgMTP7.3 ( Table 1 ). The amino acid lengths (a. a. length) of the 11 TgMTP proteins ranged from 256 aa (TgMTP6.1) to 484 aa (TgMTP7.3). The theoretical isoelectric points (pI) and molecular weights (MW) of TgMTP proteins were 5.23 (TgMTP3 and TgMTP4) - 8.33 (TgMTP5) and 27576.85 Da (TgMTP6.1) - 51803.91 Da (TgMTP7.3), respectively. The grand average of hydropathicity (GRAVY) values distributed from -0.001 (TgMTP7.3) to 0.291 (TgMTP5). Except for TgMTP7.3 with a GRAVY of -0.001, all other members had GRAVY values greater than 0, indicating that the TgMTP family proteins are mostly hydrophobic, with only one member TgMTP7.3 being hydrophilic.

The phylogenetic relationships of TgMTPs with MTP members from soybean, tomato, rice and A. thaliana were displayed in the phylogenetic tree ( Figure 1 ). The 63 MTP proteins were clustered into three major clades, designated as Mn-CDF, Fe/Zn-CDF and Zn-CDF subfamilies. The Mn-CDF clade contained TgMTP3 and TgMTP4. The Fe/Zn-CDF subfamily included 6 TgMTP members (TgMTP6.1-6.3, TgMTP7.1-7.3). TgMTP1, TgMTP2 and TgMTP5 belonged to the Zn-CDF group. Each clade also incorporated the homologous MTPs from the other plant species. The phylogenetic analysis revealed that the 11 TgMTPs were distributed among the 3 subfamilies together with respective orthologs. This suggests that the divergence of TgMTPs into different subfamilies preceded the speciation of tulip, soybean, tomato, rice and A. thaliana. The results provide insights into the evolutionary relationships and putative functional specificities of TgMTP transporters based on their classification into Mn-CDF, Fe/Zn-CDF and Zn-CDF subfamilies.

The predicted secondary structure composition of the 11 TgMTP proteins is summarized in Figure 2A . Four types of secondary structure elements were analyzed, including alpha helix, beta turn, random coil and extended strand. Among the TgMTPs, TgMTP3 and TgMTP4 contained the highest percentages of alpha helices at 63.27%, while TgMTP5 was enriched in extended strands at 12.47%. TgMTP7.1 exhibited the lowest content of random coils at 21.04%, whereas TgMTP7.3 showed the highest random coil percentage of 31.20%. The remaining TgMTPs had moderate compositions of the four secondary structures. In general, alpha helix appeared to be the predominant element in most TgMTPs, accounting for 38.96%-63.27% of their sequences. However, the proportions of alpha helices, beta turns, random coils and extended strands varied among different TgMTP members. These compositional divergences in secondary structures may contribute to functional specificities of the TgMTP transporters.

The predicted tertiary structures of the 11 TgMTP proteins are presented in Figure 2B . Homology modeling was performed using SWISS-MODEL with the best-match templates. The tertiary structures of TgMTPs displayed the typical architecture of MTP/CDF transporters, consisting of 6 transmembrane helices (TM1-6) and connecting loops.

The subcellular localizations of the 11 TgMTP proteins were predicted using Cell-PLoc-2 ( Figure 2C ). The results showed TgMTP1, TgMTP2, TgMTP3, TgMTP4, TgMTP5, TgMTP7.1, TgMTP7.2 and TgMTP7.3 were localized to the vacuolar membrane. TgMTP6.2 and TgMTP6.3 were predicted to target the plasma membrane. Whereas TgMTP6.1 was identified as a cytoplasmic protein.

The phylogenetic relationship of TgMTP genes were shown in Figure 3A , the phylogenetic tree illustrates a high correspondence between evolutionary relationships and conserved motif as well as domain architectures. To better understand the diversity and similarity of TgMTP genes, conserved motifs in TgMTP proteins were analyzed using the MEME software ( Figure 3B ). The results revealed 10 conserved motifs. TgMTP5 contained the fewest motifs with only motif 4 among all TgMTP family members. Motif 1 was only present in TgMTP members belonging to the Mn-CDF subfamily and TgMTP6.1, 6.2, 6.3 of the Fe/Zn-CDF subfamily. Motif 2 existed in all TgMTPs of the Fe/Zn-CDF and Zn-CDF subfamilies but not in the Mn-CDF subfamily. Motif 3 was specific to the Fe/Zn-CDF subfamily members. Motif 7 was only found in TgMTP6.1, 6.2 and 6.3 of the Fe/Zn-CDF subfamily. Motif 4 and 6 were absent only in TgMTP6.1, 6.2 and 6.3. Motif 5 was unique to TgMTP7.1, 7.2 and 7.3 of the Fe/Zn-CDF subfamily. Motifs 9 and 10 were present in all TgMTP members except the Mn-CDF subfamily. The sequences of all motifs are shown in Supplementary Figure 1 .

Figure 3C displays two major protein domains (PF16916 and PF01545) in the TgMTP family. PF16916 only existed in the Mn-CDF clade and TgMTP6.1-6.3, while PF01545 distributed across all members.

Gene structure analysis ( Figure 3D ) showed that TgMTP members of the same clade had higher consistency in CDS length, such as TgMTP6.1, TgMTP6.2 and TgMTP6.3. Moreover, as the identification was based on full-length transcriptome, all members displayed only one CDS region and UTR region, without intron region.

The Ka/Ks analysis revealed varied selection pressures among TgMTP gene pairs ( Table 2 ). The Ka/Ks ratio of TgMTP1 and TgMTP2 was 0.166849, significantly lower than 1, indicating purifying selection between these two members. TgMTP6.1 and TgMTP6.2 had a Ka/Ks ratio of 0.565191, while TgMTP6.1 and TgMTP6.3 showed a Ka/Ks ratio of 0.680997, slightly higher than 0.5, suggesting relatively weak positive selection. Notably, an extremely low Ka/Ks value of 0.0700296 was observed between TgMTP7.1 and TgMTP7.2, evidencing robust purifying selection. TgMTP7.2 and TgMTP7.3 also exhibited a Ka/Ks ratio of 0.194371, lower than 1, implying purifying selection as well. Collectively, most TgMTP gene pairs were under various extents of purifying selection, reflecting functional differentiation and constraints among family members. However, a relaxation of selective pressure was observed between TgMTP6.1/6.2 and TgMTP6.1/6.3. Overall, this analysis provided insights into the molecular evolution patterns and selection forces shaping diversity and conservation within the TgMTP gene family.

Gene ontology (GO) analysis was performed to predict the functional roles of the 11 TgMTP members. As shown in Figure 4 , the significantly enriched GO terms could be categorized into biological process (BP), cellular component (CC) and molecular function (MF). For BP, the major terms included metal ion transport, metal ion homeostasis, transition metal ion transport and sequestering of metal ion. Other notable terms were intracellular chemical homeostasis, chemical homeostasis and seedling development. In the CC category, key enriched terms were plasma membrane bounded cell projection and root hair. The predominant MF terms were metal ion transmembrane transporter activity, monoatomic cation transmembrane transporter activity and monoatomic ion homeostasis. These results indicate that the TgMTP family is associated with metal transport, homeostasis and sequestration, consistent with their predicted roles as metal cation transporters. The GO terms related to root and plasma membrane localization further suggest their functional importance in roots for metal acquisition and compartmentalization. This GO analysis provides an overview of putative biological roles of the TgMTP transporters and will inform further investigations into their metal transport capabilities and physiological functions in tulip.

We predicted the potential miRNA binding sites of TgMTPs using psRNATarget (Dai and Zhao, 2011) In total, 30 miRNAs were identified that may regulate TgMTPs ( Figure 5 ; Supplementary Table 1 ). Most miRNAs had multiple TgMTP targets, such as miR159a, which could target 3 TgMTP genes ( Figure 5 and Supplementary Table 1 ). In contrast, some TgMTPs could be targeted by several miRNAs. For instance, TgMTP7.1 could be targeted by 6 miRNAs, including ath-miR159a, ath-miR159b-3p and ath-miR159c ( Figure 5 ; Supplementary Table 1 ). The relationships between TgMTPs and miRNAs warrant further investigation.

The relative expression profiles of 11 TgMTP genes in response to Fe²⁺, Mn²⁺, Cd²⁺ and Zn²⁺ treatments are shown in Figure 6 . The transcript levels were examined by qRT-PCR at 0, 1, 3, 5 and 7 days after metal exposure. Significant variations in expression were observed among TgMTP genes by metal treatments over time. Overall, most TgMTPs were markedly induced by four metal (except TgMTP4 under Cd treatment), especially at later time points for most TgMTP genes. For instance, TgMTP7.1 displayed over 30-fold upregulation after Cd treatment for 1, 3, 5 and 7 days compared to 0 day. Substantial stimulation of transcript levels was also evident for TgMTP3, TgMTP6.1 and TgMTP6.3 following Cd exposure. For Fe²⁺ treatment, 72.7% (8/11) of the TgMTP genes exhibited over 10-fold upregulation in expression. Similarly, greater than 10-fold increased expression was observed for 54.5% (6/11) and 81.8% (9/11) of the genes under Zn²⁺ and Mn²⁺ exposure, respectively. These results demonstrate that most TgMTPs are responsive to heavy metal stresses. However, the percentage of TgMTP genes that could be markedly induced (with relative expression elevated by over 10-fold) by Cd²⁺ treatment was 45.4% (5/11). The data reveals the potential roles of these TgMTP genes in metal tolerance and accumulation in tulip. Further functional analyses of the responsive transporters such as TgMTP7.1 may provide insights into heavy metal detoxification mechanisms in this species.

The TgMTP7.1 gene was selected for functional validation as it exhibited the most pronounced upregulation upon Cd treatment when compared with other TgMTP genes. The TgMTP7.1 was heterologous expressed in the parental strain BY4741 and a yeast mutants that are highly sensitive to Cd (ycf1Δ), respectively. As shown in Figure 7 , yeast cells transformed with empty vector pYES2 displayed severely inhibited growth on the medium containing 5 mM CdCl₂, whether in plates ( Figure 7A ) or liquid medium ( Figure 7B ). In contrast, ycf1Δ expressing pYES2-TgMTP7.1 exhibited significantly enhanced tolerance to Cd toxicity. At the 10² dilution factor, the survival rate of the pYES2 transformants (ycf1Δ+pYES2) was 0, and at the 10³ dilution factor, the survival rate of the pYES2 transformants (BY4741+pYES2) was 0. Whereas yeast cells with TgMTP7.1 (ycf1Δ+pYES2-TgMTP7.1) grew readily at dilutions of 10¹ and 10², even the dilutions of 10³, the ycf1Δ+pYES2-TgMTP7.1 transformants showed a survival. These results demonstrate that TgMTP7.1 expression could effectively rescue the Cd sensitivity of yeast cells. The increased Cd tolerance conferred by TgMTP7.1 heterologous expression indicates it possesses Cd transport activity and likely functions in vacuolar sequestration of Cd in plants.