Genome-wide identification of SHMT family genes in C₃, C₃-C₄, and C₄ Salsoleae s.l. species

Genome-wide identification of SHMT gene family members

In this study, we utilized genome sequencing data for Salsola junatovii, Oreosalsola laricifolia, Soda foliosa, and Xylosalsola arbuscula, provided as unpublished genome data by Institute of Genetics and Developmental Biology, Chinese Academy of Sciences. We extracted protein and coding sequences (CDS) for these four Salsoleae species using the ‘GXF Sequence Extract’ and ‘Batch Translate CDS’ modules in TBtools (v2.154) (Chen et al., 2023). Reference protein sequences of Arabidopsis thaliana SHMTs were obtained from the TAIR database (https://www.arabidopsis.org/), and the Hidden Markov Model (HMM) of the SHMT conserved domain (PF00464) was downloaded from the Pfam database (https://www.ebi.ac.uk/interpro/entry/pfam/). Using the SHMT HMM model as a template, we employed HMMER (v3.0) to perform a whole-genome scan of these Salsoleae species to identify potential SHMTs. To validate these candidates, we constructed a local protein database for each Salsoleae species using BLAST (v2.14), with A. thaliana SHMT protein sequences as queries (E-value threshold set at 1e−5) (Ahmad et al., 2024). The candidate genes were further refined by integrating results from HMMER and BLAST, and protein sequences were extracted using the ‘Fasta Extract’ module in TBtools. All candidate genes were confirmed for structural integrity using the NCBI Conserved Domains Database (https://www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi). The SHMT numbers of these Salsoleae species were designated according to homologous genes in A. thaliana. Finally, molecular weight (Mw) and isoelectric point (pI) predictions for all SHMT protein-coding genes were conducted using the ExPASy online tool (https://web.expasy.org/compute_pi/) (Wilkins et al., 1999).

SHMT protein and gene structure analysis

Secondary structure predictions for SHMT proteins were performed using the SOPMA online tool (https://npsa.lyon.inserm.fr/cgi-bin/npsa_automat.pl?page=/NPSA/npsa_sopma.html). The three-dimensional structure of SHMT was predicted utilizing AlphaFold3 (https://alphafoldserver.com/), and the resulting models were simulated and visualized with PyMOL (v3.1.0) (Jumper et al., 2021; Mooers, 2020). The subcellular localization of SHMT proteins was inferred through the CELLO online tool (http://cello.life.nctu.edu.tw/) (Yu et al., 2006), while signal peptide predictions were performed using SignalP (v5.0) (https://services.healthtech.dtu.dk/services/SignalP-5.0/) (Nielsen et al., 2019).

For gene structure and conserved motif analysis, the locations of SHMT genes and their exon-intron structures in these species were extracted from the General Feature Format (GFF) annotation files and visualized using the ‘Gene Location Visualize from GTF/GFF’ function in TBtools. Conserved motifs within SHMT protein sequences were identified using MEME (https://meme-suite.org/meme/tools/meme), with the motif number maintained at 20 and other parameters set to default (Bailey et al., 2015). These results were visualized using the ‘Gene Structure View’ function in TBtools.

Phylogenetic analysis

SHMT protein sequences from A. thaliana, Glycine max, Solanum lycopersicum, Populus trichocarpa, Cucumis sativus, Beta vulgaris and Oryza sativa were retrieved from Phytozome (v13) and CuGenDB (http://cucurbitgenomics.org/). These sequences were aligned using ClustalW, with the Delay Divergent Cutoff (%) set to 30 (Liu et al., 2022), while all other options remained at their default settings. A phylogenetic tree comprising 71 SHMT protein sequences was constructed using the maximum likelihood method in MEGA-X (v10.1.8), applying the Jones-Taylor-Thornton (JTT) amino acid substitution model with uniform rates among sites (no discrete gamma categories or invariant sites), the Nearest-Neighbor-Interchange (NNI) heuristic search (initial tree generated automatically by NJ/BioNJ), and 1,000 bootstrap replicates; all other parameters were left at their defaults (Gao et al., 2024). The resulting phylogenetic tree was subsequently visualized using iTOL (v7.0) (https://itol.embl.de/).

SHMT gene family syntenetic analysis

The gene location information for SHMT family members in four Salsoleae species was analyzed using the ‘Gene Location Visualize from GTF/GFF’ function in TBtools (Chen et al., 2023). Gene density was calculated utilizing the ‘Gene Density Profile’ function. Synteny within each species was assessed using the ‘One Step MCScanX’ module in TBtools, with results visualized through the ‘Advanced Circos’ function. Furthermore, synteny relationships among these species, A. thaliana, and B. vulgaris were examined by downloading the genomic data of Arabidopsis and B. vulgaris from Phytozome (v13) and employing the ‘One Step MCScanX’ Synteny analysis plots were generated to illustrate the syntenic relationships of homologous SHMT genes across these species.

RNA extraction and reverse transcription-qPCR

Seeds from S. junatovii, O. laricifolia, X. arbuscula, and S. foliosa were collected in October 2023 in Xinjiang. The seeds were air-dried at room temperature for 2 weeks and subsequently stored at 4 °C in a refrigerator. Following sterilization (Guo et al., 2024), the seeds were sown on 1/2 MS solid medium (1/2 MS + 15 g L⁻¹ sucrose + 8 g L⁻¹ agar) and incubated in a growth chamber for 3 days. The growth chamber was maintained under a 14 h light/10 h dark cycle at 25 °C during the light period and 15 °C during the dark period. The light intensity was maintained at approximately 300 μmol·m⁻²·s⁻¹. Subsequently, healthy seedlings were selected and transferred to a Hoagland nutrient solution for hydroponic culture, with the solution being changed every 3 days. After 6 to 8 weeks of growth, leaf tissue samples were collected at 11:00 AM, immediately frozen in liquid nitrogen, and stored at −80 °C for subsequent total RNA extraction. Root, stem, and leaf tissue samples were collected from these Salsoleae species, with each tissue sample weighing approximately 100 mg. The samples were immediately frozen in liquid nitrogen and ground into a fine powder.

Total RNA was extracted using the TaKaRa MiniBEST Plant RNA Extraction Kit. The quality and concentration of the RNA were assessed using a NanoDrop spectrophotometer, and only RNA samples with a 260/280 ratio between 1.8 and 2.1 were selected for further analysis. To eliminate DNA contamination from the samples and synthesize complementary DNA (cDNA), we employed the PrimeScript™ RT Reagent Kit with gDNA Eraser. Reverse transcription- qPCR (RT-qPCR) primers specific to the SHMTs in these species were designed using the NCBI Primer Design Tool (https://www.ncbi.nlm.nih.gov/tools/primer-blast/index.cgi). The designed primers were submitted to Sangon Biotech for synthesis via the PAGE purification method. Quantitative RT-PCR was performed utilizing species-specific β-actin genes as internal reference controls for normalization (Zhang et al., 2019b), with the primer sequences for both SHMT and β-actin genes in each species provided (Table S1). RT-qPCR experiments were conducted with the TB Green Premix Ex Taq™ II kit. All experiments were carried out in the molecular laboratory of the research group at the Xinjiang Institute of Ecology and Geography Chinese Academy of Sciences (Zhang et al., 2019a). For each species, the average expression level of all gene members in the root tissue was used as the control, and the relative expression levels of SHMTs in the root, stem, and leaf tissues were calculated using the 2^−ΔΔCT method (Livak & Schmittgen, 2001). The expression data were log-transformed (log10), and the results were visualized using the ggplot2 package (v3.5.1).

Gene expression level analysis based on RNA-seq data

To further investigate the expression patterns of SHMT family members associated with different photosynthetic types (C₂, C₃, and C₄) across various plant groups, we selected five plant groups that include species exhibiting C₂, C₃, and C₄ photosynthesis based on prior studies. Subsequently, we downloaded 27 RNA-seq datasets corresponding to each of the 27 species for further analysis. These datasets comprise dicotyledonous plants from Alternanthera (Amaranthaceae), Flaveria (Asteraceae), Heliotropium (Boraginaceae), and Mollugo (Molluginaceae), as well as monocotyledonous plants from Neurachne (Poaceae) (Chinthapalli et al., 2000; Voznesenskaya et al., 2013; Stata et al., 2014; Wen & Zhang, 2015; Tao, Lyu & Zhu, 2016; Thulin et al., 2016; Lundgren, 2020; Lyu et al., 2020; Bernardo et al., 2023; Lauterbach et al., 2024). Additionally, newly sequenced RNA-seq data from four species of Salsoleae (Amaranthaceae) were included (Table S2). All RNA-seq data were generated using a paired-end sequencing strategy and underwent quality control using FastQC. The transcriptomes were assembled de novo using Trinity (v2.11.0) with default parameters (Haas et al., 2013). Subsequently, Cluster Database at High Identity with Tolerance (CD-HIT) (v4.8.1) was employed to cluster the transcript sequences at a 0.95 similarity threshold, effectively removing redundant transcripts (Fu et al., 2012). To assess the quality of the transcript assembly, Bowtie2 (v2.4.4) was utilized to map the reads back to their respective transcriptomes (Langmead & Salzberg, 2012), while Salmon (v1.10.3) was used for the quantitative analysis of the assembled transcripts (Patro et al., 2017). The assembled 31 transcripts were annotated using the Arabidopsis database (https://www.arabidopsis.org/), and the transcipts per million (TPM) values of the SHMT gene family members for each species were extracted to measure transcript abundance. Finally, stacked bar plots were generated using the ggplot2 package (v3.5.1) to display the expression levels of SHMT in Salsoleae plants and the five aforementioned genera. Individual bar plots were created for each SHMT member to illustrate the expression of SHMTs in plants with different photosynthetic types.

Prediction of cis-acting elements and transcription factor binding sites in promoter sequences

Prediction of cis-acting elements and transcription factor binding sites in promoter sequences located 2,000 bp upstream of the SHMT coding sequence (CDS) were extracted using the GTF/GFF3 Sequences Extract function in TBtools. These sequences were then submitted to the PlantCARE website (https://bioinformatics.psb.ugent.be/webtools/plantcare/html/) for predicting cis-acting elements (Lescot, 2002), while transcription factor binding sites were predicted utilizing the PlantTFDB database (https://planttfdb.gao-lab.org/) (Tian et al., 2020). The results were systematically organized and visualized through the ‘Gene Structure View (Advanced)’ and ‘Heatmap’ functions available in TBtools.

Genome-wide identification of SHMT gene family member

Based on the whole genome data of four Salsoleae species, we identified four to five members of each SHMT gene family. O. laricifolia and X. arbuscula each contained five SHMTs, including SHMT1, SHMT2, SHMT3, SHMT4, and SHMT7. In contrast, S. junatovii and S. foliosa each possessed four members, namely SHMT1, SHMT2, SHMT4, and SHMT7 (Table 1). All identified SHMTs contained the characteristic SHMT domain (Pfam: PF00464) (Table S3). Further analysis of the physicochemical properties of these SHMT members (Table 1) revealed that the amino acid lengths ranged from 458 (OlSHMT3) to 615 (SjSHMT7, OlSHMT7, SfSHMT7, and XaSHMT7). Notably, SHMT7 across these species exhibited the same amino acid length. Among these members, XaSHMT7 had the largest molecular weight at 68,405.64 Da, while OlSHMT3 had the smallest molecular weight at 50,244.21 Da. Additionally, the isoelectric points (pI) of these SHMTs ranged from 5.98 (OlSHMT7) to 8.86 (XaSHMT4). It was noteworthy that SHMT3 and SHMT7 were classified as acidic proteins (pI < 7).

Prediction of SHMT protein secondary structure and subcellular localization

The secondary structures of the SHMT proteins were analyzed (Table 2). The results indicated that the SHMT family across all four Salsoleae species was predominantly composed of α-helix, extended strands, β-turns, and random coils, with proportions ranging from 37.4% to 48.05%, 9.27% to 11.87%, 2.15% to 3.93%, and 38.13% to 49.27%, respectively. Subcellular localization analysis revealed that SHMT proteins were primarily distributed in the mitochondria, chloroplasts, and nuclei: SHMT1 and SHMT2 in each species were predicted to localize mainly to mitochondria, while SHMT3 and SHMT4 were predicted to localize to chloroplasts, and SHMT7 to nuclei. Signal peptide prediction indicated that none of the SHMT members possessed a typical signal peptide region. Additionally, to further analyze their structural characteristics, we predicted and visualized the three-dimensional structures of SHMT1, SHMT2, SHMT4, SHMT7, and SHMT3 (Fig. S1), with confidence scores are provided (Table S4). The results demonstrated that the overall fold and the distribution of the SHMT domain (pfam00464) were highly conserved among all four species, and each SHMT isoform assembled into a homotetrameric structure. However, distinct local structural variations were evident among different SHMT members; although the global fold was maintained, these local conformational differences may underlie the substrate-binding specificity or regulatory functions of each isoform.

Gene structure and conserved motif analysis of SHMT members

To further elucidate the functional differences among SHMTs, we analyzed their gene structures and conserved motifs. The conserved motif analysis revealed the presence of 20 conserved motifs, with lengths ranging from 8 to 50 amino acids (Table S5). Utilizing TBtools software, we categorized the SHMTs into four distinct classes based on phylogenetic analysis and conserved motifs (Fig. 1A). Specifically, SHMT1 and SHMT2 from each species were classified into Class IV, while SHMT3 from O. laricifolia and X. arbuscula was assigned to Class III. All SHMT4s were placed in Class II, and SHMT7s were categorized into Class I. Notably, SHMT proteins within the same group exhibited similar motif compositions (Fig. 1A), with the exception of SHMT3. Motifs 1, 3, 4, 5, and 8 were present in all SHMT members. Furthermore, the number of exons was consistent within each class (Fig. 1B). SHMTs in Class I, Class II, Class III, and Class IV contained 4, 4, 10, and 15 exons, respectively. These analyses suggested structural differences among SHMT members across different classes, implying potential functional differences.

Phylogenetic analysis of the SHMT gene family

To investigate the phylogenetic relationships among members of the SHMT gene family, this study selected seven representative species, including six dicot species (A. thaliana, G. max, S. lycopersicum, B. vulgaris, C. sativus, and P. trichocarpa) and one monocot species (O. sativa) (Table S6). Based on the topology of the phylogenetic tree constructed using the maximum likelihood method, the 71 SHMTs were classified into four classes (Fig. 2). Members of Class I, Class III, and Class IV were localized in the nucleus, chloroplast, and mitochondrion, respectively. Group II exhibited two predicted subcellular localizations: SHMTs from B. vulgaris and the Salsolee species were localized to the chloroplast, while other SHMTs were localized in the cytosol.

Collinearity analysis among Salsoleae species, A. thaliana and B. vulgaris

We investigated the gene synteny relationships among the four Salsoleae species, as well as their syntenic relationships with A. thaliana and B. vulgaris. No segmental or tandem duplication events were observed among the Salsoleae species. Chromosomal localization analysis (Table 1, Fig. S2) revealed that the SHMTs of S. junatovii and S. foliosa were distributed across two chromosomes, while the SHMTs of O. laricifolia and X. arbuscula were distributed on three chromosomes. Furthermore, in S. junatovii, O. laricifolia, and S. foliosa, the SHMT1 and SHMT7, as well as SHMT2 and SHMT4, were each located on the same chromosome, respectively. However, in O. laricifolia, the OlSHMT3 gene was located on a separate chromosome. X. arbuscula possessed five SHMTs, and exhibited a distinct distribution pattern. Specifically, XaSHMT3 and XaSHMT7 were located on the same chromosome, while the other genes were distributed across different chromosomes, with one gene (XaSHMT1) found in a scaffold region. This suggests that X. arbuscula may have undergone unique genomic rearrangements during its evolutionary history.

Through the analysis of the genetic relationships of SHMT genes among Salsoleae species, Arabidopsis, and B. vulgaris, we found that SHMT1, SHMT2, and SHMT4 exhibited strong collinearity across the six species. Notably, the collinearity network of SHMT1 encompassed all analyzed species (Fig. 3 and Table S7), suggesting that SHMT1 may exhibit a significant degree of functional conservation across different species. Conversely, SHMT3 and SHMT7 displayed significant species-specific evolutionary traits. Specifically, OlSHMT3 from O. laricifolia showed collinearity with SHMT3 from B. vulgaris, whereas XaSHMT3 from X. arbuscula completely lacked this conserved collinearity. The SHMT7 collinearity network was more limited, with a clear homologous relationship observed solely between S. junatovii and B. vulgaris. Furthermore, SHMT4 gene pairs were identified across all species, suggesting that this gene might have existed prior to the divergence of the ancestral species. Notably, Salsoleae species exhibited 1–2 pairs of homologous genes with Arabidopsis, while there were more homologous genes (2–4 pairs) with B. vulgaris. This may be attributed to the close phylogenetic relationship between Salsoleae and B. vulgaris, as both belong to the Amaranthaceae family, leading to greater homology in their SHMTs.

Quantitative real-time PCR analysis in different tissue

Quantitative real-time PCR (RT-qPCR) was employed to analyze the expression levels of SHMTs in the roots, stems, and leaves of various Salsoleae species (Fig. 4). The results demonstrated that all SHMTs were expressed in all three tissues, indicating a lack of strict tissue specificity. Among these SHMT members, SHMT1 exhibited the highest expression levels, with its transcripts predominantly accumulating in the leaves. Additionally, SHMT4 and SHMT7 demonstrated significantly higher expression in leaves compared to roots and stems. The expression pattern of SHMT2 varied across species, in O. laricifolia and S. foliosa, SHMT2 expression was higher in roots and stems than in leaves, whereas in S. junatovii and X. arbuscula, expression levels were elevated in stems and leaves compared to roots. Notably, SHMT3 was detected only in O. laricifolia and X. arbuscula, with markedly lower expression compared to other members. These findings suggest that distinct SHMT homologs may fulfill specialized functional roles in different organs across Salsoleae species.

Transcriptional expression analysis of SHMTs in different groups

To further investigate the expression of SHMT in Salsoleae and other plant groups, we analyzed SHMTs across four Salsoleae species and five additional genera: Alternanthera (Amaranthaceae), Flaveria (Asteraceae), Heliotropium (Boraginaceae), Mollugo (Molluginaceae), and Neurachne (Poaceae) (Figs. 5 and S3, Table. S2). The results indicated that, except for MpenSHMT1, FkocSHMT4, and McerSHMT7, which were not detected in Mollugo pentaphylla, Flaveria kochiana, and M. cerviana, respectively, all other species exhibited the expression of SHMT1, SHMT4, and SHMT7. Additionally, SHMT2, SHMT3, and SHMT6 were detected in some species. Although SHMT genes were widely distributed across various plant groups, significant differences existed in the number of gene members and their expression levels. Notably, the types and expression patterns of SHMT members varied among different species within these genera. In general, the expression levels of SHMT6 and SHMT7 were relatively low compared to other SHMT members, while SHMT2, SHMT3, and SHMT4 exhibited relatively higher expression levels in specific groups. Furthermore, SHMT1 was the most predominantly expressed member in most species and showed higher expression levels in C₂ and C₃ plants compared to C₄ plants.

Cis-acting elements and transcription factor binding site analysis

In our study, SHMT1 was identified as the predominantly expressed SHMT. Consequently, we conducted a detailed analysis of the 2,000 bp upstream promoter region of SHMT1 across four Salsoleae species. A total of 40 key cis-acting elements was detected, including 12 light-responsive elements, eight growth and development-related elements, 11 hormone-responsive elements, and nine stress-responsive elements (Fig. 6A). The types of cis-acting elements presented in the SHMT1 promoter region varied among these four species. Specifically, S. junatovii (C₃) exhibited the highest number of distinct cis-acting element types, with 29 identified, while O. laricifolia (C₂) had the fewest, with only 19. Further analysis revealed that the frequency of occurrence of different cis-acting elements also varied among these species. For instance, the Sp1 element appeared five times in S. junatovii, whereas its occurrence was lower in the other species. The Box-4 element occurred most frequently in O. laricifolia, with six occurrences. We further analyzed the positions of cis-acting elements within the SHMT1 promoter regions across various species. Our findings revealed significant variations in both the quantity of elements and their specific locations within these regions among different species. For example, OlSHMT1 exhibited a higher concentration of cis-acting elements near the transcription start site, whereas XaSHMT1 displayed fewer elements that were more dispersed (Fig. 6B).

Additionally, we predicted the transcription factor binding sites (Fig. 6C). The results indicated that X. arbuscula possessed only 11 transcription factor binding sites in its SHMT1 promoter, with the BasicHelix-Loop-Helices (BHLH) binding site being particularly prevalent, occurring a total of 62 times. In contrast, the SHMT1 promoter of S. foliosa exhibited the greatest diversity, identifying 18 different types of transcription factor binding sites. Our analysis further revealed that v-myb avian myeloblastosis viral oncogene homolog (MYB) binding sites were abundant in the SHMT1 promoters across all species. Moreover, certain transcription factor binding sites displayed distinct species specificity, for instance, APETALA2 (AP2), Cysteine-rich Polycomb-like (CPP), and Far-red Impaired Response 1 (FAR1) binding sites were detected solely in the SHMT1 promoter of S. junatovii. In summary, these findings reveal significant species-specific differences in the types, numbers, and distribution of cis-acting elements and transcription factor binding sites within the SHMT1 promoter among different Salsoleae species, suggesting that these elements may play unique roles in environmental adaptation and the regulation of photosynthesis.