Seeds of T. hispida sourced from The Turpan Desert Botanical Garden (Xinjiang, 293 China) were germinated in plastic pots containing a mixture of turf peat and sand (1:1 v/v) under constant photoperiod conditions (14/10 h light/dark) with a light intensity of 1,500∼2,000 lx at temperature (24 ± 1°C). After culturing for 3 months in a greenhouse, at least 800 healthy seedlings of similar size (9 cm in height) were selected for Cd treatment. Based on preliminary test results (Gao et al., 2011), 150 μmol L^–1 CdCl₂ was used to irrigate the seedlings. At the same time, samples irrigated with fresh water were treated as control. After 24, 48 or 72 h treatment (each treatment contained three separate repeats with at least 200 seedlings), the leaves were washed with clean water and frozen in liquid nitrogen immediately, after which they were stored at −80°C for subsequent experiments. Each sample contains three replicates.

To detect the Cd concentration of the samples, leaves from samples (each sample contains at least 20 seedlings) containing three replicates for each control, 24, 48 and 72 h were dried at 72°C to a constant weight and then digested with HNO₃. Subsequently, the Cd ion content was determined using ICP-OES 5110 VDV (Agilent Instruments Inc., CA, United States) (Gdbrijel et al., 2009; Han et al., 2016). The H₂O₂ content was detected by a hydrogen peroxide assay kit (Nanjing Jiancheng Bioengineering Institute), and detailed operating procedures were carried out according to the manufacturer‘s instructions. At the same time, hydrogen peroxide (H₂O₂) content in the leaves after Cd stress was detected by 3,3-diaminobenzidine (DAB) staining. Briefly, the above mentioned samples were placed in PBS (pH 7.0) solution containing 1 mg mL^–1 DAB and treated in the dark for 12 h at 37°C. After exposure for 1 h, samples were decolorized with ethanol, and finally, the seedlings were photographed (Daudi and O’Brien, 2012).

Total RNA was extracted from the leaf tissues of each sample by CTAB method (Jiang and Zhang, 2003). The degree of degradation of RNA samples was verified by RNase-free agarose gel electrophoresis. The RNA concentration was detected by qubit, and the RNA integrity was accurately detected with an Agilent 2100 Bioanalyzer. Equal quantities of high-quality RNA from the samples was used for the subsequent RNA sequencing.

Complementary DNA libraries were constructed for each of the samples and sequenced on the Illumina HiSeq 2000 platform (Illumina Inc., CA, United States). The original sequenced reads or raw reads containing reads with adapters or low quality were filtered to obtain clean reads. Specifically, reads with adapters were removed, unknown bases (N bases) over 10% and/or low-quality reads (the number of bases with a mass value of Qphred ≤ 20 accounts for more than 50% of the total reads) were removed from each data. Then, the clean reads of the twelve samples with high quality were spliced to construct unique sequences as the reference sequences using the Trinity package (Grabherr et al., 2011). The quality of transcripts were estimated by the value of FPKM (expected number of fragments per kilobase of transcript sequence per million base pairs sequenced).

The clean reads of each sample were remapped to reference sequences using RSEM software (Li et al., 2009). RSEM counts the results of the bowtie comparison and further obtains the number of read counts for each sample that was aligned to each gene and performed FPKM conversion to analyze gene expression levels. For genes with more than one alternative transcript, the longest transcript was selected to calculate the FPKM.

To infer transcriptional changes over time under Cd stress conditions, differentially expressed genes (DEGs) after 24, 48, and 72 h of Cd treatment were identified by comparing the expression levels with control. The false discovery rate (FDR) was calculated to adjust the threshold of p-value to correct for multiple testing (Rajkumar et al., 2015). Transcripts with a minimal four-fold difference in expression (| log₂Ratio| ≥2) and an FDR ≤ 0.001 were considered differentially expressed between two time points (Audic and Claverie, 1997). For convenience, DEGs with higher expression levels at 24, 48, and 72 h than control were denoted as ‘‘upregulated,’’ whereas those with the opposite were denoted as ‘‘downregulated.’’ At the same time, Venn diagrams of these differentially expressed genes were made to distinguish the differences between them using online software¹.

Short Timeseries Expression Miner (STEM) version 1.3.8 was used to analyze expression pattern (Ernst and Bar-Joseph, 2006). To further explore the temporal expression patterns, K-means clustering was applied to the identified DEGs. The DEGs belonging to the same cluster have expression patterns similar to each other. For each genotype, the clustering profiles of DEGs with p-values < 0.05 were considered to be significantly different from the reference group.

Eight genes were randomly selected for quantitative real-time RT-PCR (qRT-PCR) to determine the expression patterns revealed by RNA sequencing. RNA was extracted from T. hispida leaves of Cd treated samples (control, 24, 48, and 72 h). The PrimeScript™ RT reagent Kit (Takara) was used for first-strand cDNA synthesis. Primers for qRT-PCR analysis are listed in Supplementary Table 1. β-actin (FJ618517) was used as an internal control (Vandesompele et al., 2002). qRT-PCR was performed using a real-time PCR instrument (qTOWER 2.0) (analytik jena, Jena, Germany). The reaction mixture (20 μl) consisted of 10 μl of TransStar^R Top Green qPCR SuperMix (TRANS), 2 μl of cDNA template (equivalent to 500 ng of total RNA), 0.5 μmol L^–1 of forward and reverse primers. The reaction procedure was as follows: one cycle at 95°C for 3 min, followed by 45 cycles of 95°C for 30 s, 58°C for 15 s, 72°C for 30 s. Three independent replicates were performed to ensure the reproducibility of results. Expression levels of these genes were determined according to the 2^–ΔΔ(Ct) (Livak and Schmittgen, 2001). Then, qRT-PCR results and sequencing results were analyzed together to verify the accuracy of the sequencing results.

The unigenes were analyzed for functional annotation and functional classification. After splicing, the unigene sequences were compared with the protein database by Blast, and the annotation included the NCBI Non-redundant nucleic acid database (NR) and the Swiss-port protein sequence database (Swiss-Prot) with a threshold of e-value < 0.00001. The GO annotation information was obtained by Blast2GO analysis based on the NR annotation information (Conesa et al., 2005) and classification of all unigenes was performed by WEGO (Ye et al., 2006). At the same time, the Kyoto Encyclopedia of Genes and Genomes (KEGG) was performed to further characterize the metabolic pathways and biological functions of DEGs in the transcriptome.

The abscisic acid (ABA) signaling pathway gene ThUGT was successfully cloned using F: ATGGCTTCAGAATCCCATGAT and R: TTAGTTAATCCGGCCAC CTTT as primer, then an overexpression vector (ThUGT-pROKII) was constructed using F: GCTCTAGAATGGCTTCAGAATCCCATGAT and R: CGGGGTACCTTAG TTAATCCGGCCACCTTT following Wang et al. (2020). Afterward, the overexpression vector strain was transiently transformed into T. hispida according to the method of Ji et al. (2014) and treated with 100 μmol L^–1 CdCl₂ for 24 h, and the pROKII empty vector transiently transformed seedlings were used as control. Each treatment contained three separate repeats with at least 45 seedlings. Then the expression levels of ThUGT in transient overexpression T. hispida and control seedlings were analyzed by qRT-PCR. At the same time, the Cd content, ABA content (SenBeiJia Biological Technology Co., Ltd., Nanjing, China) and H₂O₂ content (Nanjing Jiancheng Bioengineering Institute) of the samples were determined. DAB, NBT and Evans Blue staining were also performed on each sample (Zhang et al., 2011). Each experiment was repeated three times.

In this study, the Cd and H₂O₂ contents were detected after T. hispida was treated with 150 μmol L^–1 CdCl₂. The results showed that Cd concentration showed an upward trend after Cd stress compared with control. Especially at 72 h, the concentration reached a peak with an absolute value of 41.6 mg kg^–1, which was 227 times that of control (Figure 1A). The H₂O₂ concentration was also significantly increased after Cd stress and reached the highest level at 24 h (20.2 times that of control) (Figure 1B). At the same time, the dense DAB staining of the leaves of T. hispida was observed at 24 h (Figure 1C).

Using a cutoff of four-fold difference in gene expression, 3,505, 3,983, and 4,443 upregulated genes and 3,273, 4,299, and 4,158 downregulated genes were identified at 24, 48, and 72 h, respectively, compared with those in the control (Figure 2A). Among these genes, 2,673, 3,058, and 3,077 DEGs at 24, 48, and 72 h had no similar sequences in the NR database, respectively (Supplementary Table 2). Interestingly, the expression of many DEGs changed significantly only at a certain time point. For example, there were 5,012 DEGs at 72 h and there were 1,069 DEGs common to all three treatment time points compared with the control (Figure 2B).

To validate the expression data obtained from RNA sequencing, eight genes were randomly selected from the identified DEGs to perform qRT-PCR analysis. The results showed a strong correlation between the RNA sequencing and qRT-PCR data (Figure 3), which supports the reliability of the expression results generated by RNA sequencing.

Through GO and KEGG pathway enrichment analysis, the function of the DEGs were characterized. GO annotation suggested that biological processes and molecular functions related to ROS functions and biosynthetic and metabolic processes were enriched among the DEGs at different time points (Figures 4A,B). In biological processes, biosynthetic, and metabolic processes, DEGs were enriched at all three time points (Figure 4A). In addition, the molecular function term “phenylalanine ammonia-lyase activity” was enriched at 24 and 72 h. The term “methylenetetrahydrofolate reductase NAD(P)H activity” occurred at all three time points (Figure 4B). These results indicated that these genes or proteins participate in hormone and ROS metabolism play crucial roles in the T. hispida response to Cd stress.

The Kyoto Encyclopedia of Genes and Genomes pathway enrichment analysis results showed that 90, 91, and 105 pathways were categorized from the pairwise comparisons between 24 h vs. control, 48 h vs. control, and 72 h vs. control, respectively. Based on these results, the DEGs involved in the biosynthetic pathways of seven hormones and one pathway for “plant hormone signal transduction” were enriched (Figure 4C). At the same time, we observed that the number of genes involved in six hormone synthetic pathways were the highest at 72 h, whereas brassinosteroid biosynthesis pathways peaked at 24 h. Interestingly, six hormone synthetic pathways were included in the top 30 pathways from 24 h vs. control (Table 1). In the top 20 pathways in 24 h vs. control, five hormone synthesis pathways were found. These results demonstrated that the expression of genes involved in hormone synthesis may play an important role in T. hispida response to Cd stress.

Hierarchical clustering produced six groups with similar expression trends to those of K-means clustering. Specifically, 1,069 qualifying genes were categorized into six groups (referred to as G1, G2, G3, G4, G5, and G6), comprising 240, 110, 384, 142, 164, and 29 genes, respectively (Figure 5), of which the G2, G3, G4, G5, and G6 groups mainly showed an upregulated trend under Cd treatment.

The expression of G1 group showed a downregulation trend, and reached its lowest value at 48 h, then slightly increased at 72 h. The G2 group mainly reached a peak at 24 h, then showed a downward trend and slightly increased at 72 h. The G3 group showed an upregulation trend at the beginning (24 h) followed by stable expression in the subsequent stage. The G4 group showed a continuous increasing trend at 24 h and 48 h. Interestingly, the G5 group showed an increasing expression trend during the whole stage and peaked at 72 h. The G6 group showed an upregulation trend at 24 and 48 h and a slight decrease at 72 h. From these results, we can assume that the genes in G3 and G5 group can be rapidly induced by Cd stress, and the genes in G2, G4, and G6 group showed a time-dependent trend in the process. In contrast, the genes in G1 group showed inhibition trend in response to cadmium stress.

In the GO functional annotation analysis, the number of genes in each group based on their biological pathway, molecular function and cellular component were counted. The results showed that all selected DEGs were involved in 282 biological pathways. Among them, the number of DEGs involved in protein binding, ATP binding, and DNA binding were 88, 52, and 45, respectively (Table 2).

In total, the K-means/hierarchical clustering and GO function annotation results indicated that the genes related to hormones were significantly affected during Cd stress in T. hispida leaves. Hence, the DEGs involved in the metabolism of hormones and their signaling pathways were further explored systematically.

To further explore the genes involved in hormone biochemical pathways following Cd stress treatment, the DEGs involved in the seven hormone biosynthesis or metabolism KEGG pathways (ABA, ETH, IAA, SA, GA, BRs, and JA) were analyzed (Table 1). The number of genes in the ABA, SA, ETH, IAA, GA, and JA biosynthetic or metabolism pathways primarily increased, while BR biosynthesis showed decreasing trends (Figure 4C).

Under Cd stress, 7 DEGs were identified in the ABA-related pathway, including ABA β-glucosyltransferase, zeta-carotene desaturase (PLN02487), isoprenoid biosynthesis enzymes (IBe), carotene beta-ring hydroxylase (PLN02738), antheraxanthin epoxidase/zeaxanthin epoxidase (PLN02927) and cytochrome P450. Three of them were induced, especially ThUGT (ABA β-glucosyltransferase), which appeared at 24 and 72 h with 3.7 and 4.2-fold increased, respectively. In contrast, the expression of PLN02487, PLN02738, PLN02927, and cytochrome P450 showed reduced expression trend.

Thirty DEGs in ethylene biosynthetic pathway were detected. Among the metabolic processes, there were nine beta-eliminating lyase genes, four hypothetical proteins, two 1-aminocyclopropane-1-carboxylate oxidases and two 5-methyltetrahydropterin glutamate homocysteine methyltransferases. The expression of these 17 genes were upregulated. In contrast, the expression of genes involved in homocysteine S-methyltransferase, 1-aminocyclopropane-1-carboxylate synthase, ARD/ARD’ family and S-adenosyl-methionine synthase were downregulated.

Hormones content, antioxidant activities, and downstream signals are induced in response to environmental stresses to control cell stability and mitigate the negative effects of ROS (Abdullah et al., 2021; Heidari et al., 2021). These results indicated that the key regulatory components of the biosynthetic pathways in ABA and ethylene changed significantly during the T. hispida response to Cd stress. In addition, there are also many genes involved in BRs, JA, GA, IAA, and SA, also responded to Cd stress.

To further explore the genes involved in hormone pathways following Cd stress treatment, the ABA-related KEGG pathway after each stress treatment point was analyzed. The results showed that ABA pathway-related genes mainly participated in phytoene, lycopene, zeaxanthin, abscisate and lutein biosynthetic processes (Supplementary Figure 1). ThUGT, a predicted ABA β-glucosyltransferase gene, was one of the DEGs in the ABA signaling pathway, and the expression level was induced (3.7 and 4.2 times of the control at 24 and 72 h, respectively). Therefore, ThUGT was transiently transformed into T. hispida. Pre-experimental results showed that in case of treatment with 100 μmol L^–1 CdCl₂ after 24 h, the relative expression level of ThUGT gene had the largest difference with the normal conditions (Figure 6). The results of qRT-PCR showed that the expression level of ThUGT was significantly higher in the overexpressing plants than in the control plants (Figure 7A).

The staining analysis showed that the overexpressing plants were stained lighter after CdCl₂ stress treatment compared with the control plants (Figure 7B). The results obtained with physiological indicators showed that the ABA content in ThUGT-overexpressing T. hispida and the control plants was effectively increased under 100 μmol L^–1 CdCl₂ treatment, while in the ThUGT transgenic plants, it increased less than in the control (Figure 7C). At the same time, the Cd ion content in both the overexpression plants and the control plants increased significantly after Cd treatment, but the Cd ion content in the overexpressing plants was significantly less than that of the control (Figure 7C). The H₂O₂ content in overexpression plants was lower than that of the control after Cd treatment (Figure 7C). These results suggested that the ThUGT transgenic plants eliminated more ROS and inhibited the absorption of Cd to a certain extent by T. hispida under cadmium stress, thereby enhancing the cadmium resistance of T. hispida.