Physiological and transcriptome response to cadmium in cosmos (Cosmos bipinnatus Cav.) seedlings

To date, several species of Asteraceae have been considered as Cd-accumulators. However, little information on the Cd tolerance and associated mechanisms of Asteraceae species Cosmos bipinnatus, is known. Presently, several physiological indexes and transcriptome profiling under Cd stress were investigated. C. bipinnatus exhibited strong Cd tolerance and recommended as a Cd-accumulator, although the biomasses were reduced by Cd. Meanwhile, Cd stresses reduced Zn and Ca uptake, but increased Fe uptake. Subcellular distribution indicated that the vacuole sequestration in root mainly detoxified Cd under lower Cd stress. Whilst, cell wall binding and vacuole sequestration in root co-detoxified Cd under high Cd exposure. Meanwhile, 66,407 unigenes were assembled and 41,674 (62.75%) unigenes were annotated in at least one database. 2,658 DEGs including 1,292 up-regulated unigenes and 1,366 down-regulated unigenes were identified under 40 μmol/L Cd stress. Among of these DEGs, ZIPs, HMAs, NRAMPs and ABC transporters might participate in Cd uptake, translocation and accumulation. Many DEGs participating in several processes such as cell wall biosynthesis, GSH metabolism, TCA cycle and antioxidant system probably play critical roles in cell wall binding, vacuole sequestration and detoxification. These results provided a novel insight into the physiological and transcriptome response to Cd in C. bipinnatus seedlings.

Cd accumulation and distribution. The Cd concentration was not observed in all samples under 0 μmol/L Cd stress ( Table 1). The Cd concentrations of all samples increased significantly with increasing Cd concentration. Cd accumulated highly in the roots, followed by the stems and the leaves (Table 1). Translocation factor (TF) values of the stems ranged from 0.56-0.64 was higher than those of the leaves ranged from 0. 19-0.29 (Table 1), suggesting that most of Cd in aboveground was sequestrated in the stems.
In order to understand whether different Cd stresses exhibited different Cd detoxifications or toxicity, we analyzed the Cd subcellular distribution mainly in the roots under these three Cd stresses. Under 40 μmol/L Cd stress, more than 80% Cd was accumulated in the soluble fraction, approximate 15% Cd was accumulated in the cell wall fraction, and only 5% Cd was transported into the organelle fraction (Fig. 3). Although Cd in soluble fraction was dramatically reduced with the increasing Cd concentrations, more than 55% Cd was still sequestrated in this fraction when treated with 120 μmol/L Cd (Fig. 3). Meanwhile, Cd in cell wall fractions were significantly increased with the increasing Cd concentrations, up to 40% Cd was binding in cell wall fraction when treated with 120 μmol/L Cd. These results indicated that the sequestration of Cd into soluble fraction is the main Cd detoxification mechanism under low Cd stress, while the sequestration of Cd into soluble fraction and the binding of Cd in the cell wall fraction represent a coaction for Cd detoxification with the increasing Cd concentrations.
Effects of Cd on Zn, Ca, and Fe concentrations in C. bipinnatus. After 9 days of treatments, the Cd stresses significantly decreased the uptake of Zn in the roots when compared with CK (Fig. 4A). Interestingly, Zn concentration in the stems and leaves were mainly increased, although some decreased at 120 μmol/L Cd in stems ( Fig. 4B and C). These results indicated that Cd inhibited the uptake of Zn in the roots, but promoted the translocation of Zn from root to shoot. Meanwhile, Cd stresses significantly decreased the Ca concentration in the stems and the roots ( Fig. 4D and E), but did not affect the Ca concentration in the leaves (Fig. 4F). Cd increased Fe concentrations in roots and stems leaves, although the differences were not significant ( Fig. 4G and H). But it significantly increased the Fe concentrations in leaves (Fig. 4I). These results indicated that Cd treatment may differentially affect the uptake of metal nutrients in C. bipinnatus.  (Fig. 5B). Meanwhile, the POD activity in the leaves and roots were dramatically increased ( Fig. 5C and D). Interestingly, Cd stresses did not affect the CAT activity in the leaves (Fig. 5E) and roots (except of 40 μmol/L, Fig. 5F). Except of 120 μmol/L Cd in the leaves, the SOD activity in the leaves and roots were increased by all three Cd treatments ( Fig. 5G and H). Although 40 μmol/L Cd did not affect the GR activity, 80 and 120 μmol/L Cd significantly increased the GR activity in leaves and roots ( Fig. 5I and J). These The growth of C. bipinnatus exposed to Cd. A: the fresh weight of the plants; B: the dry weight of the plants, and C: the root length. Values were means ± standard deviation (n = 3); values followed by different lowercase letters show significant differences at P < 0.05 .  results indicated that lower Cd treatment (40 μmol/L) did not cause severe oxidative stresses, while higher Cd treatment (80-120 μmol/L) produced more oxidative damages.
Transcriptome sequence and de novo assembly. Cd stresses changed the accumulation of some metal nutrients, the subcellular distribution of Cd, and the concentration and activity of some biochemical indexes. Meanwhile, some significantly changes were induced by 40 μmol/L Cd and no obvious toxic symptom was observed. Root samples treated with 40 μmol/L Cd was interestingly used to transcriptome analysis to reveal molecular response. 14.9 Gb nucleotides were generated (Table 2), which was deposited in the Sequence Read Archive (SRA) database with the accession numbers SRR3546768 and SRR3546769. 66,407 unigenes with the means length of 817 bp and N50 value of 1,344 bp were assembled. Among these assembled unigenes, the length of 18,491 unigenes (27.8% of all the unigenes) was more than 1,000 bp (Table 2). These results suggested that RNA sequencing and assembled unigenes had well quality and could be used for further transcriptome analysis.
Functional annotation and classification. 41,674 (62.76%) unigenes were functionally annotated in at least one of the five databases: GO, KEGG, COG, Swissprot and NR (Table 3). 24,733 (37.24%) unigenes were not annotated in any public database. Among these annotated unigenes, 15,481 unigenes were classified into 25 COG categories. In detail, the major group of COG was 'general functions prediction only' , followed by 'translation, ribosomal structure and biogenesis' , 'transcription' and 'replication, recombination and repair' (SFig. 1). 24,639 unigenes were classified into three major categories of GO classification. 'Cell' , 'organelle part' and 'cell part' represented the largest proportion in the cellular component category, while 'catalytic activity' and 'cell part' represented the most abundant categories in molecular function category. Moreover, the most abundant categories were 'metabolic process' , 'cellular process' and 'single-organism process' in biological process category (SFig. 2). A total of 18,496 unigenes were annotated in the KEGG database and were classified into 128 KEGG pathways (STable 2). Briefly, 'ribosome' pathway (ko03010) contains the most abundant unigenes, followed by 'carbon metabolism' (ko01200), 'biosynthesis of amino acids' (ko01230), and 'protein processing in endoplasmic reticulum' (ko04141). All sequences and functional information were deposited in the NCBI Transcriptome Shotgun Assembly database with accession number GEZQ00000000.  Noteworthy DEGs and metabolic pathways related to Cd uptake, transportation and detoxification. As shown in Fig. 6, a network associated to Cd uptake, transport, translocation, and detoxification were aggregated. Cd in extracellular could be obstructed by cell wall structures for toxicity reduction, while part of Cd would be transported by some metal transporters. When Cd entered into the C. bipinnatus cells, various metabolic processes could be induced for Cd detoxification. GSH, MTs and organic acid would be bound with Cd and then sequestrated into vacuoles. Cd also induced unigenes of antioxidant enzymes for oxidative defense. Finally, some Cd-chelates would be uploaded into xylem for transportation by some metal transporters. Briefly, this network was mainly composited by the several following processes. Firstly, total 49 metal transporters (34 up-regulated and 15 down-regulated) including zinc transporters (ZIPs), ATP-binding   Table 3. Result of unigne annotation. Figure 6. The presumed transcriptional network related to Cd uptake, translocation, and detoxification in C. bipinnatus root. The red arrows represent up-regulated genes, while the green arrows represent down-regulated genes. The dotted red boxes represent noteworthy mechanism pathways. cassette (ABC) family members, heavy-metal ATPases (HMAs) like HMA3, HMA2, the family of natural resistance-associated macrophage protein (NRAMP) members such as NRAMP2, NRAMP3, yellow stripe-like (YSL) family members, and plant cadmium resistance protein PCR2 were observed ( Table 4), suggesting that these metal transporters might participate in Cd uptake, transport and translocation. 29 DEGs (14 up-regulated and 15 down-regulated) involved in sulfate and GSH metabolism such as gene encoding adenylyl-sulfate kinase, sulfate adenylyl-transferase, serine acetyltransferase, adenylyl-sulfate reductase, sulfite reductase, serine acetyl-transferase, cysteine synthase, ornithine decarboxylase, spermidine synthase, glutathione S-transferase, and γ-glutamyl-transpeptidase 3 were regulated by Cd, suggesting the potential biosynthesis of GSH and PCs and chelation of Cd (Table 4, SFig. 5). 7 DEGs (1 up-regulated and 6 down-regulated) were metal chelates such as metallothioneins (MTs), suggesting MTs also play role in Cd chelation (Table 4). 30 DEGs (13 up-regulated and 17 down-regulated) were involved in cell wall metabolism, such as UDP-glucose 6-dehydrogenase1 (UGDH1), UDP-glucose pyrophosphorylase2 (UGP2), glucose-6-phosphate isomerase (GPI), sucrose synthase, pectinesterase, UDP-glucuronate-4-epimerase, fructokinase, beta-fructofuranosidase, xyloglucan endotransglucosylase/ hydrolase protein, expansin, glucan 1,3-alpha-glucosidase, cellulose synthase, callose synthase, and laccase, suggesting the compositions of cell wall were changed, mainly focused on pectin, callose and cellulose (Table 4). In addition, 12 DEGs (1 up-regulated and 11 down-regulated) participated in phenylpropanoid metabolism were regulated by Cd (Table 4). 17 DEGs (9 up-regulated and 8 down-regulated) were participated in tricarboxylic acid cycle (TCA) pathway, which potentially enhanced glucose metabolism and produced more organic acids such as citrate and malate for Cd binding (Table 4, SFig. 6). 11 DEGs (4 up-regulated and 7 down-regulated) were related to nitrogen metabolism, including genes encoding high-affinity nitrate transporter, nitrate reductase, ferrddoxin-nitrite reductase, glutamine synthetase, and glutamate synthase [NADPH/NADH] (Table 4). Moreover, 14 DEGs (4 up-regulated and 10 down-regulated) involved in antioxidant system, like superoxide dismutase, catalase isozyme, phospholipid hydroperoxide glutathione peroxidase, monodehydroascorbate reductase, and ascorbate peroxidase were differentially regulated by Cd (Table 4).

RT-qPCR validation.
To confirm the differential expression profiles of DEGs identified from RNA-Seq analysis, a total of 14 candidate DEGs were randomly selected from RNA-Seq and their expression levels in CK and Cd were examined by quantitative RT-PCR. As expected, the expression pattern of those unigenes obtained from qRT-PCR was similar with the differential expressions from RNA-Seq (Fig. 7).

Discussion
Under Cd stress, if the plants accumulate more than 100 μg/g Cd in dry aerial tissue, the value of TF is more than 1, and with normal growth, these plants are recommended as Cd hyper-accumulators [37][38][39] . However, most plants exhibit Cd toxicity when the leaves accumulate more than 5-10 μg/g Cd 40  Thus, C. bipinnatus should be a Cd accumulator, which would be potentially used for phytoremediation under mild Cd stress. However, the mechanism of strong tolerance and high Cd accumulation of C. bipinnatus was unknown. The physiological parameters and transcriptome analysis would help us revealing the mechanism. Normally, Cd was intake, translocated and accumulated using other metal transporters, such as Zn, Fe and Cu 41 . Although ZIP genes mainly transport Zn, some ZIPs participate in Cd transport in Arabidopsis and Thlaspi caerulescens 42 . We found that Cd up-regulated ZIP1, ZIP3, ZIP5, ZIP7, and ZTP29 ( Table 4), suggesting that these genes were involved in Cd and Zn transport, thus the Zn concentrations in roots were reduced (Fig. 4A). HMA2, a Cd/Zn transporter, loads Cd/Zn into xylem 43,44 . The increased expression level of HMA2 in Arabidopsis 45 , rice 44 , and barley 43 induced Cd or Zn xylem uploading for translocation to the shoot. HMA3, another P-type ATPase gene, segregates Cd or Zn into the root vacuolar to limit the Cd xylem loading 46,47 . Down-regulation of HMA3 resulted in a decreased concentration of Cd in the root 48 . Certainly, when HMA3 is localized at the leaf vacuolar, it also transports Cd into the leaf vacuolar, finally producing Cd hyper-accumulator of Thlaspi caerulescens 16 . Therefore, in the present study, up-regulation of HMA2 and down-regulation of HMA3 not only implied most of Cd was uploaded to aerial tissues, which resulted in high root-shoot translocation (Table 1), but also suggested a certain amount of Zn was uploaded to shoots, so that higher Zn concentrations accumulated in aerial tissues under Cd treatment ( Fig. 4B and C). NRAMP (nature resistance associated with microphage) family members display poor selectivity towards divalent mental cations, which are responsible for heavy metal ions uptake and transport 49 . Previous studies have found that NRAMP1, NRAMP3, NRAMP4 and NRAMP6 transport Fe and Cd [49][50][51][52] . NRAMP2 and NRAMP3 involved in metal efflux from the vacuole 49,53,54 . In Arabidopsis, increased expression levels of NRAMP3 result in an increased metal output from the vacuole 55 . Two NRAMP family genes (NRAMP2 and NRAMP3) were both up-regulated by Cd in our research, suggesting that they involved in Fe or Cd efflux from vacuole, thereby leading to high accumulation of Fe in aerial tissues ( Fig. 4H and I). Higher Fe concentration in shoots alleviated the Cd toxicity in Arabidopsis 56 . Therefore, the increased concentration of Fe in leaves of C. bipinnatus may also be associated with the detoxification of the plant. Additionally, NRAMP2 and NRAMP3 may also have a role in Cd transport processes. YSL genes participate in Fe-nicotinamide (Fe-NA) complex root-to-shoot transport 57 . Cd mediated the expression of YSL3 and YSL7, suggesting Cd affected Fe transport in the plant. Previous study also found that YSL3 in Solanum nigrum and YSL7 in Brassica juncea were induced by Cd as well 58,59 . However, it remains unknown whether YSL family genes are involved in transport of Cd-NA complexes transport, and further study must be conducted to analyze the function of YSL genes under Cd stress.
The subcellular distribution of Cd in the root is associated with the accumulation, translocation, and detoxification of Cd 60 . Cd bound in the cell wall fraction is an important mechanism for Cd tolerance 61 Continued of Cd was bunded in root cell wall of C. bipinnatus, especially at high Cd concentration treatment (Fig. 3). Cell wall is comprised of polysaccharide (including cellulose, semi-cellulose and pectate) and protein 64,65 . Containing abundant hydroxyl (OH-) for metal binding 66 , cellulose and semi-cellulose are both essential components of primary and secondary cell walls of higher plants 31,67 . In addition, cellulose synthase plays important role in cellulose formation, while xyloglucan endotransglucosylase (XTH) are involved in cell wall extension by cutting loosened xyloglucan strands and integrating new xyloglucans into the cell walls 68 . Present study found that the genes encoding cellulose synthase and XTH are up-regulated by Cd, suggesting the synthase and remobilization of cellulose and semi-cellulose play critical role in Cd binding. Moreover, the existence of pectin enhances the binding capacity of cell wall 69,70 . The synthase of pectin is associated with glucose metabolism. Several DEGs involved in UGP2, GPI, beta-fructofuranosidase and pectinesterase were up-regulated by Cd treatment (Table 4, SFig. 7), suggesting that Cd induce formation of pectin, thereby enhancing the capacity of Cd accumulation in cell walls, resulting high Cd tolerance of C. bipinnatus. Moreover, callose functions as a mechanical barrier to prevent ions penetration 52,[71][72][73] , and the multi-copper-containing glycoprotein laccases involved in cell wall lignification 74 . The unigenes encoding callose synthase and laccase were up-regulated by Cd treatments in our study, implied that C. bipinnatus accumulated callose deposition and enhance cell wall lignification in root to prevent Cd entering the protoplasts under Cd stress. These results indicated that cell wall obstruction is one of important detoxification mechanism in C. bipinnatus. After Cd entered into cytoplasm, it would be bound with metal chelates and then sequestrated into vacuoles to reduce their toxicity. The result of Cd subcellular distribution demonstrated that a large proportion of Cd was found in soluble fractions, suggesting the vacuole was the predominant detoxification sink for Cd in C. bipinnatus root. Vacuole possesses abundant sulphur-rich peptides such as GSH and PCs [75][76][77] , and organic acids 78 , which binds heavy metals and decreases their migration to reduce toxicity. Interestingly, unigenes involved in glutathione (GSH) metabolism, including serine acetyltransferase, cysteine synthase, ornithine decarboxylase, spermidine synthase, glutathione S transferase, γ-glutamyltranspeptidase, and unigenes from ABCC family were up-regulated under Cd treatment (Table 4). Meanwhile, ABCC1, ABCC2 and ABCC3 are major vacuolar PC-Cd transporters in other plants 75,79 , which were also up-regulated. Thus, Cd in cytoplasm turned into GSH (PC)-toxic compounds, finally transported by ABCC transporters without displaying cytotoxic to plant cell [80][81][82] . Vacuole sequestration of Cd also plays main role in Phytolacca Americana 83 and Arachis hypogaea 84 . Moreover, GSH is also one of important antioxidants in plants. AsA -GSH cycle system can be able to eliminate ROS in many plants. Genes encoding enzymes involved in AsA-GSH cycle like glutathione peroxidase, monodehydroascorbate reductase and L-ascorbate peroxidase were up-regulated by Cd. Similarly, the activity of GR under Cd stress significantly increased compared with CK. These results indicated that enhancement of AsA-GSH cycle improved Cd tolerance of C. bipinnatus.
MDA is the product of lipid peroxidation, and its concentration reflects the degree of oxidative damage. In our study, 40 μmol/L Cd did not increase the MDA concentrations in leaves, while MDA concentration increased in roots of C. bipinnatus with higher accumulation of Cd ( Fig. 5A and B), indicating that C. bipinnatus had strong tolerance under lower Cd but suffered cellular oxidative stress at higher Cd. Activities of antioxidant enzymes are induced by oxidative stress, and increased antioxidant level prevent oxidative damages 85 . Previous studies illustrated their functions in scavenging ROS in plants 86,87 . The enzyme SOD alters O 2 − to H 2 O 2 and oxygen 88 . CATs convert H 2 O 2 to water and molecular oxygen, while PODs have a more elevated affinity to H 2 O 2 than CATs 89 . The coordination between different enzymes can alleviate oxidative stress in the plant. In our study, compared with CK, the activity of POD and SOD in leaves were significantly increased under 40 μmol/L Cd treatment (Fig. 5C and G), while other two enzymes did not show significant changes ( Fig. 5E and I), suggested that the activities of SOD and POD possessed sufficient capacity to scavenging ROS under lower Cd treatment. Therefore, the MDA concentration did not increase under 40 μmol/L Cd. However, generation of ROS were increased with the increased of Cd concentrations, exceeding the limits of POD and SOD scavenging ability (Fig. 5I). Meanwhile, the GR activity increased, complementing the ability of ROS scavenging. In addition, two unignenes encoding peroxidase were up-regulated under Cd treatment. These results demonstrated that antioxidative enzymes indeed play an important role in Cd detoxification and enhance tolerance of C. bipinnatus under adverse environment.
Summary, C. bipinnatus was recommended as a "Cd-accumulator" that would be potentially used for phytoremediation under mild Cd stress. Subcellular distribution of Cd displayed different detoxification mechanisms under different levels of Cd stress. C. bipinnatus initiated diverse defense and detoxify response to keep strong tolerance when treated with Cd stress. RNA-Seq analysis revealed that ZIPs, NRAMPs, HMAs, and ABC transporters were involved in Cd uptake, translocation and accumulation. Meanwhile, several processes such as cell wall biosynthesis, glutathione (GSH) metabolism, TCA cycle and the antioxidant system probably played critical roles in cell wall binding, vacuole sequestration and detoxification.

Materials and Methods
Plant culture and Cd treatment. Cosmos seeds (Cosmos bipinnatus Cav.) were sterilized with 2% NaClO for 20 min then rinsed with deionized water. The sterilized seeds were germinated on clean sand at 25 °C. After 2 weeks, the uniform seedlings were transplanted into plastic pots with 2.5 L half-strength Hoagland nutrient (60 plants per pot, pH 6.5) for 7 days and this was then replaced with full Hoagland nutrient solution. The plastic pots were randomly divided into 4 groups, each in triplicate. The four groups were treated with control, 40, 80, and 120 μmol/L CdCl 2 , respectively. All plants grew in a growth chamber with a daily temperature of 25 °C, a relative humidity of 70% and a photon flux density of 500 μmol/m 2 ·S. Leaf, stem and root samples of all the treatments were collected on the 9th day. The root samples from 10 plants (10 plants per replicate, three biologic repeats) were collected and rapidly frozen in liquid nitrogen and then stored at −80 °C for RNA extraction.
Phenotype characterization. On the 9th day after treatment, the plant height and root length were measured (12 plants per biological replicate, three biologic repeats). The fresh weight of root, stem and leaf were also weighted. After weighing, all tissues were then dried at 80 °C for two days for dry weight calculation and metal concentration measurement.

Measurement of metal concentrations and calculation of translocation factor (TF).
The concentration of Cd, Zn, Fe, and Ca was measured as described by Wang et al. 90 with some modifications 90  The superoxide dismutase (SOD) activity was determined according the method described earlier 94 . The reaction mixture consisted of 50 mmol/L Tris-HCl buffer (pH 7.8), 0.1 mmol/L EDTA, 0.1 mmol/L nitroblue tetrazolium (NBT), 13.37 mmol/L methionine, and 0.1 mmol/L riboflavin and enzyme extract. The reaction was initiated by adding the riboflavin. The mixture was first placed under light then transferred into darkness immediately and the absorbance recorded at 560 nm. One unit of SOD activity was defined as the amount of enzyme that inhibited 50% of NBT photoreduction.
The catalase (CAT) activity was assayed in a reaction mixture containing 2.9 mL 50 mmol/L Tris-HCl buffer (pH 7.0), 50 μL 750 mmol/L H 2 O 2 , and 50 μL enzyme extract as per the method of Aebi (1984) 95 . Activity was measured by following the decomposition of H 2 O 2 at 240 nm.
The peroxidase (POD) activity was determined according to the guaiacol method 96 with some modifications. The reaction mixture was 50 mmol/L Tris-HCl buffer (pH 7.0) containing 0.1 mmol/L EDTA, 10 mmol/L guaiacol, 5 mmol/L H 2 O 2 and 100 μL enzyme extract. The reaction was initiated by adding the extract. Guaiacol oxidation was determined based on an increase in the absorbance at 470 nm. One unit of POD activity was expressed as units (μmol guaiacol decomposed per minute) per mg of fresh weight (FW).
The glutathione reductase (GR) activity was assayed as described by Foyer and Halliwell (1976) with some modifications 97 . The reaction mixture consisted of 450 μL of the enzyme extract, 2.34 mL 50 mmol/L Tris-HCl buffer (containing 0.1 mmol/L EDTA, 5 mmol/L MgCl 2 pH 7.5), 60 μL 10 mmol/L NADPH and 150 μL 10 mmol/L oxidized glutathione (GSSG). The reaction was initiated by adding the extract, NADPH, and GSSG. The NADPH oxidation rate was determined by recording the decrease in absorbance at 340 nm. The GR activity was expressed as the amount of enzyme needed to oxidize 1 μmol of NADPH /min· mg FW. Subcellular distribution of Cd in the C. bipinnatus root. Cd subcellular distribution was determined according to Su et al. 84 with some modifications 84 . The frozen root samples (1 g) were ground into powder with a pre-cold extraction buffer [50 mmol/L Tris-HCl buffer solution (pH 7.5), 250 mmol/L sucrose, 1.0 mmol/L DTE (C 4 H 10 O 2 S 2 ) and 5.0 mmol/L ascorbic acid]. The homogenate was centrifuged at 4000 r/min for 15 min and the precipitate was designated as a cell wall fraction consisting mainly of cell walls and cell wall debris. The supernatant solution was further centrifuged at 16000 r/min for 45 min. The resultant deposit and supernatant solution were designated as the organelle-containing fraction and the soluble fraction, respectively. All fractions were dried and then digested in 5 mL HNO 3 . The Cd concentrations in the different fractions were analyzed by FAAS. RNA extraction. The total RNA of each root sample (CK, 40 μmol/L Cd treatment) was extracted by using the Quick RNA isolation Kit (Huayueyang Biotech Co., Ltd., Bejing, China) according to the instruction manual. RNase-free DNasel (TaKaRa Biotech Co., Ltd., Dalian, China) was used for removing residual DNA in the extracted RNA. The quality of the total RNA sample was measured by 1% agarose gels, and the concentrations of the total RNA samples were assayed with an Agilent 2011 Bioanalyzer (Agilent Technologies, Inc., Santa Clara, CA, USA).

Library construction and illumina sequencing.
High-quality RNA samples from C. bipinnatus root were prepared for cDNA library construction and sequencing. The mRNA was purified from total RNA using oligo (dT) magnetic beads and poly (A) tails. The RNA sequencing libraries were generated using the TruSeq RNA sample Prep Kit (Illumina, San Diego, CA) with multiplexing primers, according to the manual. The cDNA library was constructed with average inserts of 250 bp, with non-stranded library preparation. The QIAquick PCR extraction kit (Qiagen, Inc., Hilden, Germany) was used for cDNA purifying. The short cDNA fragments were subjected to end repair, adapter ligation, and agarose gel electrophoresis filtration. Subsequently, the appropriate fragments were selected as templates for PCR amplification. Sequencing was performed via a paired-end 125 cycle rapid run on 2 lanes of the Illumina HiSeq. 2500 system, generating pairs of reads of great quality as intended.
Transcriptome assembly. Adapter-related and low-quality reads including ambiguous reads ('N'), duplicated sequences were removed from the raw reads to obtain the clean reads. Trinity software (http://trinityrnaseq.sourceforge.net/) was used for the de novo assembled transcriptomes. In brief, the contigs were formed by combining the certain overlap length into long fragments without N (contigs) and then they were clustered using the TGICL software to produce unigines (without N) and finally the redundancies were removed to obtain non-redundant unigenes 98 .
Unigene functional annotation. A series of databases and software were used for putative unigenes annotations. BLAST software 99 was used to align the unigene with the NR 100 , Swiss-Prot 101 , GO 102 , COG 103 , and KEGG databases 104 (E-value ≤ 1E −5 ) to retrieve protein functional annotations based on sequence similarity. The ESTScan software was used to decide the sequence direction of the unigenes that could not be aligned to any of the above databases 105 . Functional categories of putative unigenes were grouped using the GO and KEGG databases.
Differential expression analysis. FPKM values were used to compare gene expression differences between the two samples. The DESeq package was used to obtain the base mean value for identifying DEGs. FDR ≤0.01 and the absolute values of log2 ratio ≥1 were set as the thresholds for the significance of the gene expression difference between the two samples.
Real-time quantitative (qRT-PCR) validation of partial DEGs. qRT-PCR was performed in a 96-well plate with the CFX-96 real-time system (Bio-Rad, CA, USA). Each reaction of 15 μL contained 6.3 μL (30 ng/μL) cDNA, 0.6 μL (4 pmol/μL) for each forward and reverse primer, and 7.5 μL iQ SYBR Green Supermix (Bio-Rad, CA, USA). Each cDNA sample was amplified in triplicates. The PCR reaction conditions were 95 °C for 5 min, 39 cycles of 95 °C for 15 s, 56 °C for 30 s, and 72 °C for 10 s, followed by the generation of a dissociation curve by increasing the temperature starting from 65 °C to 95 °C to check for the specificity of amplification. Actin was used to standardize the transcript levels in each sample. The relative expression level was calculated with the 2 −△△CT formula 106 . The primers that were designed and used in the RT-qPCR analyses are shown in STable 1.