This study was conducted at the Jiaojiang Campus of Taizhou University, Zhejiang Province, China. Seeds of S. glauca and S. salsa were procured from the Qiaohua Tamarix Planting Co., Ltd., Shandong, China. In April 2024, the seeds of the two halophytes were surface-sterilized with a 0.5% NaClO solution for 15 minutes, washed three times with deionized water, and sown in prepared soil to germinate. They were kept in the greenhouse at a temperature of 25/20 °C with a relative humidity of 60%, under 16 h day/8 h night conditions.

In April 2024, saline-alkaline soil was collected from saline-alkali land in the Taizhou Economic Development Zone in Zhejiang, China (121°34′55″ E, 28°31′29″ N). This area is located on the south-eastern coast and has a typical subtropical monsoon climate. This has resulted in patterns of low-lying terrain and strong tidal action, as well as saline-alkali land, which have been influenced by evaporation and precipitation. The average annual temperature is 18 °C, and the annual rainfall is 1400 to 1800 mm. The soil was air-dried, sieved (2 mm mesh), and homogenized to remove coarse particles to facilitate plant growth.

The peat was supplied by Taizhou Agricultural Trade Co., Ltd. in Zhejiang, China. It was made from moss residues and had a pH value between 4 and 6. It contained over 95% organic matter and had a conductivity of 0.1 mS cm⁻¹ (Fu et al., 2021).

To evaluate the effects of peat and halophytes in the phytoremediation of saline-alkaline soil, two halophytes (S. glauca and S. salsa) were cultivated under three different peat concentrations (0, 6, and 18 g/kg soil) following previous studies (Duan et al., 2023; Fu et al., 2021). The low concentration (6 g/kg) was derived by converting the optimal field application rate of wood peat for saline-sodic soil (19.57–23.95 t·ha^-1, Duan et al., 2023) using a standard conversion coefficient (1 t·ha^-1 ≈ 0.385 g/kg soil), yielding a converted pot concentration range of 7.53–9.21 g/kg. The high concentration (18 g/kg) was consistent with the pot-validated effective concentration (18.5 g/kg, Fu et al., 2021), which has been proven to significantly improve degraded soil properties and promote plant growth. Soil properties were identical to those reported in our previous study (Wang et al., 2025b). Uniformly grown seedlings were transplanted into plastic pots (10.5 cm×14 cm), each containing 1 kg of pretreated field soil, with one plant per pot. Plants were grown in a greenhouse under controlled conditions (28 °C light for 16 h, 18 °C dark for 8 h, and approximately 45% relative humidity). All pots were watered daily with the same amount of tap water.

Leaf chlorophyll content was determined using the 80% acetone method. The extracted solution was analyzed using a multifunctional enzyme marker (Readmax 1900, Flash, Shanghai, China) at wavelengths of 645 and 663 nm. The contents of chlorophyll a (C_a) and chlorophyll b (C_b) were calculated as follows: C_a = 12.72×A663-2.59×A645, C_b = 22.88×A645-4.67×A663 (Lichtenthaler and Wellburn, 1982). Three leaves per plant from four randomly selected plants were measured repeatedly.

Four fresh leaves per plant were sampled, and every three plants were pooled to form one biological replicate. Each treatment included three biological replicates. All samples were immediately frozen in liquid nitrogen and then stored at −80 °C until use.

Total RNA was extracted from 80 mg of leaf tissue using the ethanol precipitation protocol and the CTAB PBIOZOL reagent. RNA concentration was quantified using a NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, CA, USA), and RNA quality was assessed using the Bioanalyzer 2100 system (Agilent Technologies, CA, USA). Sequencing libraries were prepared using NEB Next^® Ultra^™ RNA Library Prep Kit for Illumina^® (NEB, USA) following the manufacturer’s instructions. Subsequently, libraries were sequenced on an Illumina Nova Seq 6000 (Illumina Inc., San Diego, CA, USA) to produce 150 bp paired-end reads (Novogene, China).

Initially, clean reads were obtained by removing reads containing adapters, poly-N, and low-quality reads from raw reads. At the same time, Q20, Q30, and GC content of the clean reads were calculated, ensuring high-quality clean reads for downstream analyses. Then, clean reads were used for transcriptome de novo assembly using Trinity (v2.6.6) with min_kmer_cov set to three and min_glue set to four, along with all the other parameters being set to default, producing reference sequences for the downstream analysis (Grabherr et al., 2011). The longest transcripts within each cluster were regarded as a unigene to eliminate redundant sequences. The function of these unigenes were annotated in the following databases: NCBI non-redundant protein sequences (Nr, e-value = 1e-5), NCBI non-redundant nucleotide sequences (Nt, e-value = 1e-5), Protein family (Pfam, e-value = 0.01), Clusters of Orthologous Groups of proteins (KOG/COG, e-value = 1e-5), A manually annotated and reviewed protein sequence database (Swiss-Prot, e-value = 1e-5), KEGG Ortholog database (KO, e-value = 1e-5), and Gene Ontology (GO, e-value = 1e-6) databases.

Differentially expressed genes (DEGs) between groups were identified using the DESeq 2 (1.26.0) (Love et al., 2014) with the thresholds of p<0.05 and |log₂(FoldChange)| ≥ 1. GO and KEGG pathway enrichment analyses were performed using the cluster Profile R package v3.10.1 to infer the putative functions of the DEGs. GO terms and KEGG pathways with p < 0.05 were considered significantly enriched (Mao et al., 2005). TBtools was used to create gene expression heatmaps (Chen et al., 2020). Raw transcriptome sequencing data have been deposited in the NCBI Sequence Read Archive (SRA) under BioProject accession number PRJNA1189186.

After two months, the entire plants were harvested and washed with sterilized distilled water. Then, the plants were carefully separated into leaves, stems, and shoots. Each component was then dried at 70 °C until a constant weight was reached, after which the biomass was accurately measured using a precision balance.

The RNA samples used for RT-qPCR were consistent with those for RNA-seq. Total RNA was converted to complementary DNA (cDNA) using FastKing RT Kit (With gDNAase) (Tiangen, China) according to the manufacturer’s instructions. The qPCR was performed on a CFX Connect Real-time system (Bio-Rad, USA) via a FastFire qPCR premix (SYBR Green) kit (Tiangen, China) followed previously established protocols (Zhou et al., 2023). Three biological replicates and three technical replicates were analyzed with the EF1 gene as an internal control. Relative expression was calculated using the 2−ΔΔCT method (Livak and Schmittgen, 2001). The primers for RT-qPCR are listed in Supplementary Table 1.

Statistical data analysis was performed using the Statistical Product and Service Solution (SPSS) software (v16.0). The effects of peat addition on chlorophylls and growth parameters were analyzed using one-way analysis of variance (ANOVA). Subsequently, Duncan’s multiple range test was performed at a 0.05 confidence level to examine differences among treatments.

Peat addition significantly increased chlorophyll a and b contents of S. glauca, while it had no significant effect on the chlorophyll a content in S. salsa (Figures 1a, b; Table 1). A significant interaction was observed between peat addition and plant species for chlorophyll b content (Table 1), indicating that the effect of peat on chlorophyll b accumulation varies with plant species. At the low peat concentration (6 g/kg), no significant effects were detected on the growth of either S. glauca or S. salsa. In contrast, the high peat concentration (18 g/kg) significantly increased aboveground, belowground, and total biomass in both species (Figures 1c–e). S. glauca exhibited significantly higher aboveground, belowground, and total biomass than S. salsa. High peat concentration also significantly increased root length and root surface area in S. glauca (Figures 1f, g), although no significant effect was observed on plant height (Figure 1h). Furthermore, a significant interaction between peat addition and plant species was found for both aboveground and total biomass (Figures 1c, e; Table 1).

A total of 116.77 Gb clean reads were generated, with an average of 6.5 Gb per library, a mean Q30 percentage of 94.39%, and a mean GC percentage of 43.09% (Supplementary Table 2). As genomes of S. glauca and S. salsa were not available, the clean reads were assembled using Trinity, resulting in 81,510 unigenes for subsequent analysis. The assembled unigenes ranged in length from 301 to 15,968 bp, with an average length of 1,154 bp, and an N50 value of 1,777 bp (Supplementary Figure 1a). All unigenes were aligned against seven databases based on sequence similarity using BLASTX. As a result, 49.48%, 41.3%, 40.68%, 40.68%, 39.58%, 20.91%, and 17.23% of the unigenes were annotated in the NR, Swiss-Prot, GO, PFAM, NT, KO, and KOG databases, respectively. In total, 55,197 (67.71%) unigenes were annotated in at least one of the databases, and 4,219 (5.17%) unigenes were annotated across all seven databases (Supplementary Figure 1b).

Principal component analysis (PCA) and principal coordinates analysis (PCoA) revealed that peat addition caused significant changes in gene expression profiles in both S. glauca (Figure 2a; Supplementary Figures 2a, c) and S. salsa (Figure 2b; Supplementary Figures 2b, d), with a significant difference between low and high peat treatments (p < 0.05). DEGs were identified based on the criteria of |log₂(FoldChange)| ≥ 1 and p < 0.05. In S. glauca, 2,631 DEGs (1,579 upregulated and 1,052 downregulated) and 13,379 DEGs (10,034 upregulated and 3,345 downregulated) were identified under low and high peat concentrations, respectively. In S. salsa, 2,284 DEGs (1,490 upregulated and 794 downregulated) and 5,501 DEGs (3,267 upregulated and 2,234 downregulated) were identified under low and high peat concentrations, respectively (Figure 2c). In both species, the number of DEGs increased with the increasing levels of peat concentrations, with more genes being upregulated and fewer downregulated. High peat levels caused more DEGs than low levels. Notably, only 73 DEGs (0.4%) were between the two species across both low and high levels of peat addition (Supplementary Figure 2c). In summary, peat addition has exerted significant and concentration-dependent regulatory effects on the gene expression of S. glauca and S. salsa, and the transcriptomic responses of the two species to this treatment are highly species-specific. This difference may reflect their distinct molecular strategies for adapting to peat-amended environments.

Additionally, KEGG and GO enrichment analyses of the DEGs showed distinct responses to peat addition between S. glauca and S. salsa (Figure 2d; Supplementary Figure 3). In S. glauca, upregulated DEGs were mainly enriched in growth-related pathways. At low peat concentrations, genes related to ribosome, tricarboxylic acid (TCA) cycle, and nitrogen metabolism pathways were activated. In contrast, high peat concentrations induced genes related to fatty acid biosynthesis, ribosome biogenesis in eukaryotes, and porphyrin and chlorophyll metabolism pathways. In S. salsa, the upregulated DEGs were mainly enriched in secondary metabolite-related pathways. Both low and high peat treatments activated genes involved in linoleic acid metabolism, flavonoid biosynthesis, and cutin, suberine, and wax biosynthesis. High peat concentrations further induced the biosynthesis of phenylpropanoids, anthocyanins, and other secondary metabolites (Figure 2d).

In S. glauca, the downregulated DEGs were mainly enriched in amino acid degradation and metabolism, as well as secondary metabolites-related pathways. Specifically, both low and high peat concentrations caused the expression of genes related to the degradation of valine, leucine, and isoleucine; the metabolism of cysteine, methionine, arginine, proline, and cyanoamino acid; as well as the biosynthesis of phenylpropanoids and diterpenoids. In S. salsa, the downregulated DEGs were mainly enriched in growth-related pathways. Both peat treatments caused the expression of genes involved in photosynthesis, carbon fixation in photosynthetic organisms, etc. (Figure 2e).

A total of 30 cation transporter genes associated with Na⁺ homeostasis and K⁺ absorption were identified among DEGs, covering 7 functional families: 12 CHX (cation/H⁺ exchanger), 9 VHA (vacuolar H⁺-ATPase), 3 AKT1 (inward rectifying K⁺ channel), 2 NHX (vacuolar Na⁺/H⁺ antiporter), 2 HKT (high-affinity K⁺ transporter), 1 SOS1 (plasma membrane Na⁺/H⁺ antiporter), and 1 SKOR (stelar K⁺ outward rectifier). The expression patterns of these transporters showed pronounced species-specific responses to peat addition. In S. glauca, 5 VHA genes were significantly upregulated under low peat concentration, whereas no other transporter families exhibited differential expression. Under high peat concentration, 19 transporters were differentially expressed, including 10 downregulated genes (8 CHX, 1 SOS1) and 9 upregulated genes (7 VHA, 1 SKOR, 1 AKT1). In S. salsa, multiple transporter families were upregulated under the low peat treatment, including 5 CHX, 2 AKT1, 2 NHX, 2 VHA, and 1 HKT. Under high peat concentration, all 9 detected transporters were upregulated, covering CHX (3), AKT1 (2), NHX (2), VHA (1), and HKT (1) (Supplementary Figure 2d).

Heatmap analysis further revealed that VHA family genes were consistently upregulated in both species under peat addition, indicating a conserved role in salt tolerance (Figure 2f). In contrast, CHX family genes displayed opposite expression trends between the two species: they were downregulated in S. glauca under the high peat but upregulated in S. salsa under both low and high peat concentrations. Notably, SOS1, a key gene mediating Na⁺ extrusion or xylem Na⁺ retrieval in salt-stressed plants, was specifically downregulated in S. glauca under the high peat condition, which may reflect a reduced demand for active Na⁺ efflux under ameliorated saline conditions. On the contrary, NHX and HKT, the key genes responsible for vacuolar Na⁺ sequestration and long-distance Na⁺ transport, were specifically upregulated in S. salsa (Figure 3b). The species-specific expression patterns of K⁺ channel genes (AKT1 and SKOR) further highlighted divergent ion homeostasis strategies between the two species: under the high peat treatment, S. glauca exhibited upregulated expression of both AKT1 and SKOR, which facilitated K⁺ homeostasis maintenance; in contrast, S. salsa showed upregulation of AKT1 alone, with no differential expression of SKOR under peat treatments, suggesting a distinct mode of K⁺ regulation (Figure 2f).

In S. glauca, the upregulated DEGs were mainly enriched in growth-related pathways, such as nitrogen metabolism and the TCA cycle. Based on KEGG enrichement analysis and Swiss-prot annotations, a total of 20 DEGs involved in nitrogen metabolism pathway were identified, including seven carbonic anhydrase (CA) genes, three formamidase genes, four glutamate synthase (GLT1) genes, two glutamine synthetase (GLUL) genes, two nitrate reductase (NR) genes, one glutamate dehydrogenase (GLUD1-2), and one ferredoxin-nitrite reductase (NirA) gene (Figure 3a). Expression analysis revealed that nine of these genes were upregulated under the low peat treatment, and 15 were upregulated under the high peat treatment. Only one gene, Cluster-16558.28329_formamidase, was downregulated under high peat addition but upregulated under the low peat condition. Moreover, four genes (Cluster-16558.12787_GLT1, Cluster-16558.13182_GLUD1-2, Cluster-16558.13629_CA, and Cluster-16558.13839_GLUL) were consistently upregulated under both low and high peat concentrations (Figure 3a; Supplementary Table 3).

A total of 30 DEGs involved in TCA cycle were identified, including five isocitrate dehydrogenase (IDH) genes, four succinyl-CoA synthetase alpha subunit (LSC) genes, three pyruvate dehydrogenase E1 component (aceE) genes, three dihydrolipoyllysine-residue acetyltransferase (DLAT) genes, three malate dehydrogenase (MDH) genes, three succinate dehydrogenase (ubiquinone) flavoprotein subunit (SDH) genes, two aconitate hydratase (ACO) genes, and one gene each for dihydrolipoyl dehydrogenase (DLD), citrate synthase (CS), ATP citrate (pro-S)-lyase (ACLY), 2-oxoglutarate dehydrogenase E2 component (dihydrolipoamide succinyltransferase) (DLST), 2-oxoglutarate dehydrogenase E1 component (OGDH), Fumarse (FumA) and pyruvate carboxylase (pyc) (Figure 3b). Of these, 16 genes were upregulated under the low peat treatment, and all 30 showed upregulation under the high peat condition. Notably, 16 genes, including three MDH, two ACO, two aceE, two LSC, and one each of OGDH, SDH, DLD, DLST, DLAT, ACLY, and CS, were upregulated under both low and high peat concentrations (Figure 3b; Supplementary Table 4).

In S. salsa, the upregulated DEGs were mainly enriched in secondary metabolites-related pathways, including flavonoid, phenylpropanoid, and anthocyanin biosynthesis. According to KEGG enrichment analysis and Swiss-prot annotations, a total of 43, 20, and four DEGs were identified as being involved in phenylpropanoid, flavonoid, and anthocyanin biosynthesis, respectively (Figure 4a). Among the phenylpropanoid-related genes, six were upregulated under low peat addition, and 42 were upregulated under high peat addition. Five of which were upregulated across both low and high peat additions, including three PRDX6 and two C3’H (Figure 4a; Supplementary Table 5). For flavonoid biosynthesis, eight and 19 DEGs were upregulated under low and high peat levels, respectively, with six genes upregulated under both conditions. These included two CHS and two C3’H, and one gene each for F3H, FLS, and CHI (Figure 4a; Supplementary Table 6). Regarding anthocyanin biosynthesis, four DEGs, all encoding BZ1, were upregulated specifically under high peat addition (Figure 4a; Supplementary Table 7).

In contrast, the downregulated DEGs in S. salsa were mainly enriched in growth-related pathways. Both low and high levels of peat addition caused the expression of genes related to photosynthesis, carbon fixation in photosynthetic organisms, etc. A total of 37 DEGs involved in the photosynthesis pathway were identified, including 16 DEGs related to photosystem I (PsaD, PsaE, PsaF, PsaG, PsaH, PsaK, PsaL, PsaN, and PsaO), 13 DEGs associated with photosystem II (PsbB, PsbM, PsbQ, PsbP, PsbQ, PsbR, PsbS, PsbW, and PsbY), four DEGs related to F-type ATPase (b, delta, and gamma), three DEGs involved in photosynthetic electron transport (PetE and PetF), and one DEG from the cytochrome b6/f complex (PetA). Additionally, 12 DEGs involved in the photosynthesis-antenna proteins pathway were identified, such as Lhca1, Lhca2, Lhca4, Lhca5, Lhcb1, Lhcb2, Lhcb4, Lhcb5, and Lhcb7 (Figure 4b).

Peat addition of peat suppressed the degradation and metabolism of amino acids in both S. glauca and S. salsa. Based on KEGG enrichment analysis and Swiss-prot annotations, 18 and 15 DEGs in S. glauca were involved in valine, leucine, and isoleucine degradation under low and high peat levels, respectively. Among them, 12 DEGs were downregulated at both peat levels, including three 3-methylcrotonyl-CoA carboxylase alpha subunit (MCCC), two aldehyde dehydrogenase (NAD+) (ALDH), two 2-oxoisovalerate dehydrogenase E2 component (dihydrolipoyl transacylase) (DBT), and one each of branched-chain amino acid aminotransferase (ilvE), isovaleryl-CoA dehydrogenase (IVD), 2-oxoisovalerate dehydrogenase E1 component subunit alpha (BCKDHA), and alanine-glyoxylate aminotransferase 2 (AGXT2) (Supplementary Figure 4a; Supplementary Table 8). In S. salsa, 10 and 18 DEGs were identified under low and high peat treatments, respectively. Seven of these were downregulated under both levels of peat addition, including three MCCC, and one each of ALDH, DBT, IVD, and HPD1 (3-hydroxyisobutyrate/3-hydroxypropionate dehydrogenase) (Supplementary Figure 4b; Supplementary Table 9).

TFs such as WRKY, bHLH, bZIP, NAC, MYB, and AP2/ERF play crucial roles in regulating plant responses to abiotic stress. To better understand TFs involved in nitrogen metabolism, the TCA cycle, and the biosynthesis of flavonoids, phenylpropanoids, and anthocyanins, as well as amino acid degradation and photosynthesis, we conducted a correlation analysis between these TFs and relevant DEGs. This analysis identified 166 TFs, including 42 bHLH, 38 MYB, 30 bZIP, 29 AP2/ERF, 17 NAC, and 10 WRKY (Supplementary Table 10), which are closely associated with the key metabolic pathways mediating saline-alkaline tolerance and growth. In S. glauca, 25 TFs were detected under low peat addition and 82 under high peat addition, with 12 TFs shared between the two conditions (Supplementary Figure 5). Notably, 16 bHLH, 16 bZIP, 6 MYB, 3 WRKY, and 2 AP2/ERF TFs were significantly correlated with 18 nitrogen metabolism-related genes (e.g., GLT1, GLUL, CA), which are critical for nitrogen uptake and assimilation (Figure 5a; Supplementary Table 11). Additionally, 14 bHLH, 10 bZIP, five MYB, and one AP2/ERF TFs were tightly linked to 27 TCA cycle genes (e.g., IDH, MDH, LSC, Figure 6b; Supplementary Table 12), which provide energy precursors for growth. These TFs likely act as positive regulators of the growth-promoting pathway in S. glauca, directly driving the upregulation of nitrogen metabolism and TCA cycle genes under high peat concentration, thereby supporting biomass accumulation and root trait improvement.

In S. salsa, 19 TFs were identified under low peat addition and 86 under high peat addition, also with 12 shared TFs (Supplementary Figure 5). A total of 20 MYB, nine NAC, seven AP2/ERF, three bHLH, two WRKY, and two bZIP TFs were significantly correlated with 30 genes involved in the biosynthesis of flavonoids, phenylpropanoids, and anthocyanins (e.g., PRDX6, C3’H, BZ1)—the core stress-mitigating pathways in this species (Figure 5c; Supplementary Table 13). For the downregulated photosynthesis-related genes (e.g., Lhca, Lhcb), nine AP2/ERF, four bHLH, two NAC, two bZIP, and one MYB TFs were found to be correlated (Figure 5d; Supplementary Table 14), suggesting these TFs may mediate the transcriptional repression of photosynthetic pathways, allowing metabolic resources to be redirected toward secondary metabolite synthesis.

Notably, only one TF (Cluster-16558.33486_MYB) was present across all four libraries (Supplementary Figure 5), indicating that species-specific TF repertoires are responsible for the distinct transcriptional strategies of S. glauca and S. salsa. Furthermore, TFs involved in amino acid degradation suppression were identified in both species: in S. glauca, five AP2/ERF, two MYB, two bHLH, and one WRKY TFs were correlated with 20 genes involved in valine, leucine, and isoleucine degradation (e.g., MCCC, ALDH, IVD, Supplementary Figure 4c; Supplementary Table 15); in S. salsa, 13 AP2/ERF, nine MYB, eight bHLH, six NAC, six bZIP, and two WRKY TFs were linked to 20 amino acid degradation-related genes (e.g., MCCC, ALDH, DBT, Supplementary Figure 4d; Supplementary Table 16). This suggested that peat addition-induced suppression of amino acid degradation is also transcriptionally regulated by conserved and species-specific TF families, which helps conserve nitrogen resources and reduce energy consumption under salt stress. Collectively, these results revealed that high peat concentration modulates the expression of species-specific TFs, which in turn target the core metabolic pathways (growth-promoting pathways in S. glauca and stress-mitigating pathways in S. salsa) and the conserved amino acid degradation pathway.

To validate the reliability of the RNA-seq data and clarify the expression patterns of key genes involved in core metabolic pathways, we performed RT-qPCR on 7 representative genes, including those associated with nitrogen metabolism (Cluster-16558.13629, CA), TCA cycle (Cluster-16558.8660, SDH; Cluster-16558.14519, CS), amino acid degradation (Cluster-16558.35715, IVD; Cluster-16558.35550, DBT), and secondary metabolite biosynthesis (Cluster-16558.37546, CHS; Cluster-16558.45693, PRDX6, Figure 6). The results showed that the relative expression levels of all detected genes by qPCR were completely consistent with their FPKM values from RNA-seq. Specifically, in S. glauca, growth-promoting pathway genes (CA, SDH, and CS) exhibited concentration-dependent upregulation (Figures 6a–c), their expression in high-peat (JN3) was significantly higher than in low-peat (JN1) and control (JN0), while amino acid degradation genes (IVD, and DBT) were downregulated (Figures 6d, e). In S. salsa, the expressions of secondary metabolite biosynthesis genes (CHS, and PRDX6) were notably upregulated, peaking in high peat (YN3). In contrast, the expressions of amino acid degradation genes (IVD, and DBT) were also downregulated (Figures 6g, h). Collectively, these results confirmed RNA-seq accuracy and validated the concentration-dependent regulation of key pathways and species-specific transcriptional strategies of the two Suaeda species under peat addition.