Transcriptomic Analysis of Cadmium Stressed Tamarix hispida Revealed Novel Transcripts and the Importance of Abscisic Acid Network.

Cadmium (Cd) pollution is widely detected in soil and has been recognized as a major environmental problem. Tamarix hispida is a woody halophyte, which can form natural forest on the desert and soil with 0.5 to 1% salt content, making it an ideal plant for the research on response to abiotic stresses. However, no systematic study has investigated the molecular mechanism of Cd tolerance in T. hispida. In the study, RNA-seq technique was applied to analyze the transcriptomic changes in T. hispida treated with 150 μmol L-1 CdCl2 for 24, 48, and 72 h compared with control. In total, 72,764 unigenes exhibited similar sequences in the Non-redundant nucleic acid database (NR database), while 36.3% of all these unigenes may be new transcripts. In addition, 6,778, 8,282, and 8,601 DEGs were detected at 24, 48, and 72 h, respectively. Functional annotation analysis indicated that many genes may be involved in Cd stress response, including ion bonding, signal transduction, stress sensing, hormone responses and ROS metabolism. A ThUGT gene from the abscisic acid (ABA) signaling pathway can enhance Cd resistance ability of T. hispida by regulating the production of ROS under Cd stress and inhibit absorption of Cd. The new transcriptome resources and data that we present in this study for T. hispida may facilitate investigation of molecular mechanisms governing Cd resistance.

Introduction

Heavy metal pollution in soil has become a worldwide problem. It not only inhibits crop growth and reduces yield and quality but also poses a considerable threat to human health (Wu et al., 2007). The heavy metal cadmium is a biologically non-essential element. Cd is highly toxic and migratory with a biological half-life of 10 to 30 years (Suwazono et al., 2009). Therefore, Cd can easily enter the human body through the food chain and accumulate in the body, causing injury to the kidneys, lungs, liver, testicles, brain, bones, and blood system (Hamada et al., 1991).

High concentration of Cd is toxic to plants, but some plants grow under high Cd stress without exhibiting toxic effects. Approximately 10–33% of Arabidopsis halleri subsp. gemmifera can accumulate more than 100 mg kg–1 of Cd in contaminated soil (Bert et al., 2002). Yang et al. (2004) showed no reduction in shoot and root dry matter yields when S. alfredii were grown at Cd levels of 200 μmol L–1 in nutrient solution. Under natural conditions, the aboveground part of Noccaea caerulescens can accumulate up to 164 mg kg–1 Cd (Baker et al., 1994). T. praecox is a hyperaccumulator plant of zinc, cadmium and lead, and the aboveground part can accumulate up to 5,030 mg kg–1 of Cd (Vogel-Mikus et al., 2005). Under Cd pollution level of 25 mg kg–1, the Cd content in stems and leaves of S. nigrum exceeded 100 mg kg–1, and it was greater in shoots than in roots (Wei et al., 2005). These studies suggested that these super enriched plants with good cadmium tolerance can provide the theoretical basis for the study of plant remediation of Cd-contaminated soil.

At present, it is generally believed that accumulation of Cd in plants is primarily reflected in two aspects. On the one hand, at the cellular level, Cd primarily accumulates in the vacuoles and apoplasts of plants. On the other hand, at the organ level, this process is manifested in the epidermal cells, subepithelial cells and epidermal hairs of plants. According to Küpper et al. (2000), mustard mesophyll cells are important sites for Cd accumulation. In addition, Salt and Wagner (1993) reported deposit of large amount of Cd in the leaf epidermis and epidermis hairs in mustard. Recently, researchers found that the vacuolar membrane of rapeseed and A. thaliana play an important role in regulating the ion channel protein activity of NO3– and Cd, and the vacuolar compartmentalization and cell wall fixation of Cd may be the main physiological reasons for the difference in Cd toxicity resistance between Cd-resistant cultivar Z11 and Cd-sensitive cultivar W10 of rapeseed, which provides a means of synergistically improving the NUE and Cd toxicity of rapeseed (Zhang et al., 2019b). In A. thaliana, the defensive protein AtPDF2.5 may chelate cytoplasmic Cd and mediate its efflux, promote Cd accumulation in apoplasts, and regulate plant detoxification and accumulation of Cd (Luo et al., 2019). Therefore, accumulation of Cd in cell walls, vacuoles, epidermal cells or epidermal hair is likely to be one of the ways in which plants achieve detoxification.

In recent years, with of availability of transcriptional data, a growing body of knowledge regarding the genetic basis underlying Cd stress physiological processes has greatly increased our understanding of the molecular mechanism of Cd transcription and toxicity in some Cd hyperaccumulating plants, such as A. halleri (Herbette et al., 2006), Brassica juncea (Farinati et al., 2010), S. alfredii (Gao et al., 2013) and Noccaea caerulescens (Halimaa et al., 2014; Milner et al., 2014). At the same time, the molecular mechanisms of Cd stress on some cultivated plants such as cabbage (Brassica oleracea subsp. capitata f. alba) (Ba̧czek-Kwinta et al., 2019), pea (Pisum sativum L.) (Rodriguez-Serrano et al., 2009), barley (Hordeum vulgare L.) (Cao et al., 2014), rice (Oryza sativa L.) (Oono et al., 2014), tobacco (Nicotiana tabacum L.) (Martin et al., 2012), ramie (Boehmeria nivea L.) (Liu et al., 2015) and pakchoi (Brassica chinensis L.) (Zhou et al., 2016) were studied.

Tamarix hispida is a woody halophyte that grows in arid and semiarid regions. In a previous study, the transcriptome of T. hispida treated with NaHCO3 was constructed and analyzed to detect the response of T. hispida to alkaline treatment (Wang et al., 2013). Some transcription factors, including ThNAC7, ThCRF1, ThZFP1, and ThbHLH1, involved in the process of reducing ROS to confer salt or osmotic tolerance in transgenic plant have been cloned (Zang et al., 2015; Ji et al., 2016; Qin et al., 2017; He et al., 2019). There were studies in which multiple T. hispida genes enhanced tolerance to Cd. For example, the transfer of the metallothionein gene ThMT3, increased resistance to Cd in transgenic tobacco and yeast (Yang et al., 2011; Zhou et al., 2014). Overexpression of vacuolar membrane H+-ATPase c subunit gene ThVHAc1 improved Cd tolerance of Saccharomyces cerevisiae, Arabidopsis, and T. hispida. ThWRKY7, a possible upstream gene of ThVHAc1, exhibited similar expression patterns as ThVHAc1 under CdCl2 treatment and improved Cd tolerance in T. hispida (Gao et al., 2011; Yang et al., 2016).

Investigating transcriptomic response of Cd-stressed leaves would be particularly useful for furthering the genetic improvement of T. hispida to Cd stress. To elucidate the initial perception mechanism in response to Cd stimuli in T. hispida leaves, we examined gene expression changes at different time points and identified Cd-specific regulatory networks. This study helps to elucidate the mechanism of Cd tolerance in T. hispida and provides a useful reference for further exploration in woody plants and used for remediation of heavy metals (Cd) from contamination soils.

Materials and Methods

Plant Materials and Cadmium Treatments

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 CdCl2 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.

Determination of Cadmium Concentration and H2O2-Related Physiological Indices

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 HNO3. 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 H2O2 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 (H2O2) 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).

Ribonucleic Acid Extraction, Sequencing and de novo Assembly

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).

Normalization of Gene Expression Levels and Identification of Differentially Expressed Genes

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 (| log2Ratio| ≥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[1].

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.

Validation of Differentially Expressed Genes With Quantitative Real-Time Polymerase Chain Reaction

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 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.

Sequence Annotation, Functional Classification, and Biological Pathway Analysis

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.

Cloning of ThUGT Gene in Abscisic Acid Signaling Pathway and It’s Cadmium Resistance Function Analysis

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 CdCl2 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 H2O2 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.

Results

Changes of Cadmium and H2O2 Concentrations in Tamarix hispida Subjected to Cadmium Stress

In this study, the Cd and H2O2 contents were detected after T. hispida was treated with 150 μmol L–1 CdCl2. 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 H2O2 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).

FIGURE 1

Cd stress analysis of T. hispida. The leaves Cd (A) and H2O2 (B) concentrations in T. hispida at control or under 150 μmol L– 1 CdCl2 treated for 24, 48, or 72 h. (C) DAB staining of T. hispida leaves under control or after cadmium stress. * (P < 0.05) indicate signification difference compared with control.

Pairwise Comparisons of Transcriptome Between Control and Cadmium Stressed Leaves

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).

FIGURE 2

DEGs of T. hispida under cadmium treatments. (A) Gene number analysis of DEGs between the cadmium treated transcriptomes compared with the control. (B) Venn diagrams of these DEGs.

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.

FIGURE 3

Verification of eight selected DEGs by qRT-PCR. Comparison of RNA-seq data (Red point) with qRT-PCR data (Blue point). The normalized expression level (FPKM; expected number of Fragments Per Kilobase of transcript sequence per Millions base pairs sequenced) of RNA-seq is indicated on the y-axis to the left. The relative qRT-PCR expression level of selected DEGs is shown on the y-axis to the right. β-actin was used as the internal control. Three biological replicates were used.

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.

FIGURE 4

GO terms and KEGG pathways involved in ROS and hormones production analysis. GO analysis of biological process terms (A), molecular function terms (B,C) KEGG pathways involved in hormones between the cadmium treated transcriptomes compared with the control. The x-axis in (A,B) indicates the percentage of DEG snumbers vs. background gene numbers in each GO term. The x-axis in (C) indicates the percentage of DEGs numbers vs. background gene numbers in each KEGG pathway.

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.

TABLE 1

Top 30 KEGG pathways based on the percentage of DEGs in 24 h vs. control.

Term	Total	24 h vs. control		48 h vs. control		72 h vs. control
		DEGs	%	DEGs	%	DEGs	%
Glucosinolate biosynthesis	7	2	28.6	1	14.3	2	28.6
Photosynthesis – antenna proteins	48	10	20.8	1	2.1	10	20.8
Cutin, suberin, and wax biosynthesis	77	7	9.1	2	2.6	9	11.7
Brassinosteroid biosynthesis (BRs)	31	2	6.5	0	0.0	1	3.2
Flavonoid biosynthesis	73	4	5.5	1	1.4	11	15.1
Phenylpropanoid biosynthesis	372	19	5.1	12	3.2	49	13.2
Stilbenoid, diarylheptanoid, and gingerol biosynthesis	61	3	4.9	0	0.0	10	16.4
Cyanoamino acid metabolism	161	7	4.4	6	3.7	13	8.1
Phenylalanine metabolism (SA)	123	5	4.1	1	0.8	17	13.8
Nitrogen metabolism	104	4	3.9	4	3.9	11	10.6
Arachidonic acid metabolism	57	2	3.5	1	1.8	5	8.8
Phenylalanine, tyrosine, and tryptophan biosynthesis	206	7	3.4	8	3.9	9	4.4
alpha-Linolenic acid metabolism (JA)	184	6	3.3	4	2.2	10	5.4
Tryptophan metabolism (IAA)	134	4	3.0	1	0.8	7	5.2
Alanine, aspartate and glutamate metabolism	241	7	2.9	4	1.7	22	9.1
Ascorbate and aldarate metabolism	214	6	2.8	7	3.3	19	8.9
Pentose phosphate pathway	215	6	2.8	7	3.3	15	7.0
Amino sugar and nucleotide sugar metabolism	408	11	2.7	6	1.5	24	5.9
Cysteine and methionine metabolism (ETH)	342	9	2.6	9	2.6	23	6.7
Terpenoid backbone biosynthesis	190	5	2.6	3	1.6	4	2.1
Taurine and hypotaurine metabolism	76	2	2.6	0	0.0	7	9.2
Steroid biosynthesis	123	3	2.4	1	0.8	8	6.5
Plant hormone signal transduction	633	15	2.4	11	1.7	30	4.7
Butanoate metabolism	128	3	2.3	1	0.8	2	1.6
Photosynthesis	131	3	2.3	3	2.3	12	9.2
Carbon fixation in photosynthetic organisms	306	7	2.3	8	2.6	28	9.2
Plant-pathogen interaction	482	11	2.3	13	2.7	20	4.2
Glutathione metabolism	264	6	2.3	8	3.0	18	6.8
Pentose and glucuronate interconversions	221	5	2.3	4	1.8	13	5.9
Sphingolipid metabolism	133	3	2.3	4	3.0	5	3.8

Hormone-related KEGG pathways have a gray background.

Analysis of Temporal Expression Pattern of Genes

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.

FIGURE 5

Analysis of temporal expression pattern of genes. (A) Heat map analysis of DEGs. (B) Expressive trend pattern analysis of DEGs. The value of ratios take the base 2 logarithm were used to analysis each gene‘s expression trend in every group. For each gene, ratios = FPKM of the gene in sample/FPKM value of the gene in control.

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).

TABLE 2

The top 20 GO molecular function based on number of common DEGs.

Gene Ontology Molecular Function	DEG_item
protein binding	88
ATP binding	52
DNA binding	45
zinc ion binding	33
nucleic acid binding	31
oxidoreductase activity	28
protein kinase activity	24
metal ion binding	21
structural constituent of ribosome	20
RNA binding	18
transcription factor activity, sequence-specific DNA binding	17
calcium ion binding	14
transmembrane transporter activity	14
catalytic activity	13
hydrolase activity, hydrolyzing O-glycosyl compounds	12
GTP binding	10
electron carrier activity	9
GTPase activity	9
heme binding	9
nucleotide binding	8

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.

Differentially Expressed Genes Involved in Hormone Biosynthetic Pathways

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.

ThUGT-Overexpressing Tamarix hispida Increased Abscisic Acid Content and Cadmium Tolerance

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 CdCl2 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).

FIGURE 6

The expression levels of ThUGT gene in transiently transformed T. hispida under different CdCl2 stress conditions.

FIGURE 7

Cadmium tolerance analysis of ThUGT gene. (A) qRT-PCR analysis (B) staining analysis and (C) physiological indicators determination of ThUGT transgenic T. hispida under 100 μmol L– 1 CdCl2 stress. * (P < 0.05) indicate signification difference compared with control.

The staining analysis showed that the overexpressing plants were stained lighter after CdCl2 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 CdCl2 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 H2O2 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.

Discussion

When plants are under heavy metal stress, it can produce a series of responses to relieve the toxic effects of heavy metals. Some reports suggested that the production of glutathione, phytochelatin, and metal chelates or chaperones, which can bind to heavy metal ions and then transport them out of cells. In addition, antioxidant enzymes, such as superoxide dismutase (SOD), peroxidase (POD), glutathione reductase (GR), ascorbate peroxidase (APX), and glutathione peroxidase (GPX), were produced to scavenge oxygen free radicals generated by heavy metals (Emamverdian et al., 2015).

In this work, we used RNA-seq to explore the time course of the response mechanism in T. hispida under Cd stress. Our study showed that 114,292 unigenes were obtained after the transcriptome data were spliced and 36.3% of the unigenes were novel transcripts under Cd stress. Further analysis found that 6,778 DEGs were detected at 24 h and then increased to 8,282 and 8,601 at 48 and 72 h, respectively (Figure 2A), which indicated that an increasing number of DEGs were induced or activated under Cd stress in T. hispida.

H2O2 is a product of aerobic metabolism of cells, and its production is increased under various stresses. It not only has the effect of damaging biological macromolecules and thereby harming cells, it is also an important signal molecule that induce the expression of a series of defense genes in the cell, improve the activity of protective enzymes to remove active oxygen, prevent its excessive accumulation under adversity conditions, and protect plants from damage (Moller, 2001; Mittler, 2002). Interestingly, Cd content increased with increasing treatment time (Figure 1A), while the H2O2 content peaked at 24 h (Figures 1B,C). At the same time, we also found that many studies reported H2O2 functions a secondary messenger during plant development and defense to abiotic stress (Avshalumov et al., 2003; Jiang et al., 2012; Khalili et al., 2014; Liu et al., 2016; Saxena et al., 2016). Therefore, we suggest that H2O2 may be serving as a secondary messenger to induce the expression of stress-related proteins in the early stage of cadmium stress and initiate the development of systemic acquired resistance in T. hispida. This plant synthesizes many ROS-clearing genes, thereby effectively eliminating excess H2O2 to help T. hispida tolerate Cd stress. Consistent with this property, the results of the DEGs at 24 h showed that 1/6 DEGs participated in the redox reaction, the transport of ions, and the synthesis of signal substances (Supplementary Table 3), which have important functions in plant stress resistance.

In previous studies, multiple transcription factors were reported to be involved in the scavenging of ROS or increasing SOD and POD activities in T. hispida to improve the salt tolerance or osmotic stress ability of plants (Zang et al., 2015; Ji et al., 2016; Qin et al., 2017; He et al., 2019). Many studies on transcriptome analyses of Cd-treated plants have found that it is mainly related to the pathways of “ROS-scavenging enzymes” (Guo et al., 2017; Gupta et al., 2017), “cell wall alternation and strengthening” (Fan et al., 2011; Wan and Zhang, 2012; Xu et al., 2015), “lipid oxidation” (Guo et al., 2017), “auxin biosynthesis and metabolism” (Yue et al., 2016), and “nitric oxide-mediated homeostasis” (Zuccarelli et al., 2017).

In T. hispida, the GO analysis showed in the biological pathway that “oxidation–reduction process” was the top enriched term, and “oxidoreductase activity” in the molecular function also ranked in the top 10 (Figures 4A,B). These results indicated that the expression of antioxidant and redox homeostasis-related genes play important role in response to Cd stress in T. hispida.

Abscisic acid, a widely known phytohormone involved in the plant response to abiotic stress, plays a vital role in mitigating Cd2+ toxicity in herbaceous species. Studies have found that when plants were exposed to Cd, endogenous ABA levels were increased in plant cells (Sharma and Kumar, 2002). Several other studies demonstrated that the application of ABA can reduce Cd accumulation in crops (Hsu and Kao, 2003, 2005; Uraguchi et al., 2009; Fan et al., 2014). Fan et al. (2014) reported that ABA treatment correlates with the downregulation of ABA-inhibited IRON-REGULATED TRANSPORTER 1 (IRT1) to decrease Cd accumulation. Through interaction with MYB49, the ABI5 represses MYB49 binding to the downstream genes bHLH38, bHLH101, HIPP22, and HIPP44, which result in the inactivation of IRT1 and reduced Cd uptake (Zhang et al., 2019a).

The results of KEGG pathway analysis showed that seven main phytohormones participate in Cd stress in T. hispida. Among the phytohormones, the input numbers vs. background gene numbers for ABA pathways increased over time (Figure 4C), with ratios of 0.9, 2.7, and 3.5% being observed at 24, 48, and 72 h, respectively. By GO function analysis, we found that seven DEGs were involved in the ABA signaling pathway, three of which were upregulated (Supplementary Figure 1). In particular, the expression of ThUGT gene was clearly induced. The ABA content in transient overexpression ThUGT T. hispida plants was significantly increased after cadmium stress treatment, but it was less compared with the control (Figure 7C). Moreover, increase in the H2O2 content and Cd content was less in OE plants than the control after Cd stress (Figure 7C), which was consistent with the results reported in previous studies, indicating that the ThUGT gene has the ability to remove ROS under Cd stress and enhance the tolerance of T. hispida to Cd stress.

Uridine diphosphate-glucosyltransferases (UGTs) are a family of proteins involved in physiological responses to the inactivation of many glycosylation hormones (Kleczkowski et al., 1995). The overexpression of UGT74E2 gene in A. thaliana can improve the resistance to drought and salt stress by regulating the ABA dynamic balance (Tognetti et al., 2010). ABA β-glucosyltransferase belongs to the UGT family and is a key enzyme in the ABA catabolism binding pathway. This enzyme plays an important role in maintaining the normal physiological level of ABA. The ABA β-glucosyltransferase gene AtUGT71B6 in A. thaliana regulates intracellular ABA balance (Priest et al., 2006; Dong et al., 2014). The Phaseolus vulgaris PvABAGT gene regulates the ABA balance and stress response in adversity during bean development (Xu et al., 2002).

In Beta vulgaris, the ROS produced by a plasma membrane NADPH oxidase may act as a signal to induce BvGT (UGT family gene) expression after wounding and bacterial infiltration (Sepulveda-Jimenez et al., 2004). Oxidative stress and conditions that promote cell death could induce the expression of glucosyltransferase genes and produce transportable glucosides that function as ROS scavengers (Mazel and Levine, 2002). Therefore, we concluded that ThUGT may play a crucial role in the ABA conjugation pathway and in adaptation to Cd stress by inducing ROS scavengers to deduce H2O2.

Conclusion

The present study identified novel transcripts, gene structures, and DEGs in T. hispida under Cd stress. In total, 114,292 unigenes were identified. The transcripts identified in this study will serve as a valuable genomic resource for future studies. Among these genes, a large number of genes were related to ROS clearance and hormone signals, which may facilitate the analysis of ABA signaling pathways. The overexpression of related a gene ThUGT reduced the accumulation of Cd in T. hispida under high CdCl2 stress. These results will help to establish a foundation for future research to improve cadmium tolerance in T. hispida and other plants.

RNA-Seq Data Submitted to Public Database

The data presented in the study are deposited in the NCBI repository, accession number PRJNA795701 (https://www.ncbi.nlm.nih.gov/bioproject/PRJNA795701).

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author Contributions

P-LW and C-QG conceived and designed the experiments and wrote the manuscript. P-LW and Y-YW performed the experiments. P-LW and X-JL analyzed the data. All authors provided editorial advice.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.