A transcriptomic analysis of Chrysanthemum nankingense provides insights into the basis of low temperature tolerance

RNA-Seq libraries and reads mapping

An overview of the RNA-Seq reads derived from the six libraries using Illumina HiSeq™ 2000 platform is presented in Table 1 and Additional file 1: Figure S1. The raw sequence data have been deposited in the NCBI Sequence Read Archive (http://trace.ncbi.nlm.nih.gov/Traces/sra_sub/sub.cgi). The number of clean reads per library ranged from 7.01 to 7.47 million, and the total number of nucleotides sequenced from 343,646,114 to 365,975,267 (Accession No. for library A SRS591717; Accession No. for library B1 591719; Accession No. for library B2 591720; Accession No. for library C1 591721; Accession No. for library C2 591722; Accession No. for library CK 591679). The proportion of clean reads was over 99% in each library (Additional file 1: Figure S1). Overall, from library A to library CK, 7,270,059, 7,052,023, 7,013,186, 7,228,380, 7,299,665, and 7,468,883 clean reads were obtained respectively. Correspondingly, 356,232,891, 345,549,127, 343,646,114, 354,190,620, 357,683,585, and 365,975,267 total base-pairs (Table 1) were generated. The clean reads were mapped onto a reference gene database, which included all known C. nankingense unigene sequences. Raw sequence data were deposited in the NCBI Sequence Read Archive database under the accession number SRP041330. The proportion of unambiguously mapped reads per library ranged from 77.13% in library C2 to 81.03% (C1), and the proportion of unique matches from 57.64% (A) to 60.89% (C1). As the number of reads increased, the identification rate of new genes slowed, which indicated the saturation around seven million reads (7,013,186 in library B2, 7,468,883 in CK) (Additional file 2: Figure S2).

Quantification of transcripts and identification of differentially transcribed genes (DTGs)

The quality of the RNA-Seq dataset is assessed by gene coverage, which is the percentage of a gene covered by reads. This value is determined as the ratio of the base number in a gene covered by unique mapping reads to the total bases number of that gene. The distribution of the six libraries was presented in Additional file 3: Figure S3. In addition, transcript abundances for each gene (Additional file 4: Table S1) were calculated according to the method following Mortazavi et al. [31]. Moreover, differential transcription was identified through pair-wise comparison between various libraries, by setting a threshold FDR of 0.001 and a│log2 ratio│ of 1 based on the algorithm developed by Audic et al. [32]. From the seven comparisons, including treatment CKA (CK vs A), CKB1 (CK vs B1), CKB2 (CK vs B2), CKC1 (CK vs C1), CKC2 (CK vs C2), AC1(A vs C1), and AC2(A vs C2), the results showed that a large number of DTGs were identified (Additional files 5, 6, 7, 8, 9, 10 and 11: Table S2-8). The number of DTGs detected was as follows: treatment CKA, 3,779 (2,096 up- and 1,683 down-regulated); CKB1, 337 (250 and 87); CKB2, 718 (571 and 147); CKC1, 3,722 (2,271 and 1,451); CKC2, 4,119 (2,611 and 1,508); AC1, 194 (169 and 25); and AC2, 111 (92 and 19) (Figure 1). These results indicated that more DTGs were identified in the treatments, which underwent a prior CA (A, C1 and C2), as compared to the treatments which didn’t undergo CA (B1 and B2). In addition, a smaller number of DTGs was found in A vs C1 and A vs C2 comparisons than in CK vs B1 and CK vs B2. Moreover, on extending the freezing treatment, fewer DTG were detected in A vs C2 comparison than in A vs C1. However, contrary results were obtained in case of CK vs B1 and CK vs B2 comparisons. Based on the assumption that genes with similar expression patterns usually exhibit functional correlation, the consistency of the DTGs was checked by multiple comparisons clustering among the CKA, CKC1 and CKC2; between the CKB1 and CKB2 treatments; and between the comparisons, A vs C1 and A vs C2. A total of 2,340 DTGs were observed in the first multiple comparison clustering, out of which, only three genes behaved inconsistently (that is, showed up-regulation in one treatment and down-regulation in the other, or vice versa). Of the 2,337 consistent DTGs, 1,410 were up- and 927 were down-regulated (Additional file 12: Table S9). In the second multiple comparison clustering, 188 DTGs (142 up- and 46 down-regulated) were obtained (Additional file 13: Table S10), and all were found to be consistent. In the third comparison clustering between A vs C1 and A vs C2, 38 DTGs (37 up- and 1 down-regulated) showed consistency (Additional file 14: Table S11).

GO classification of differentially transcribed genes

In the treatment CKA, 1,535 of the 3,779 DTGs could be assigned a GO term; the equivalent number for other comparisons were as follows: treatment CKB1, 155/337; CKB2, 246/718; CKC1, 1,522/3722; CKC2, 1,691/4119; A vs C1, 60/194; and A vs C2, 30/111 (Additional file 15: Table S12). For CK vs A, 21 GO classes fell into the categories “biological process”, 12 into “cellular component” and 11 into “molecular function”. The equivalent distribution in CK vs B1 was 18, 10, and 7; in CK vs B2, 20, 11 and 9; in both, CK vs C1 and CK vs C2, 21, 12, and 11; in A vs C1, 17, 9, and 5; and in A vs C2, 11, 7, and 5. The major classes of biological process among the DTGs in CK vs A comparison were “metabolic process”, “cellular process”, “single organism process”, “response to stimulus”, “localization”, “establishment localization”, “biological regulation” and “regulation of biological process”; the predominant cellular components were “cell”, “cell part”, “organelle” and “organelle part”; and for molecular function “binding”, “catalytic activity”, “transporter activity”, “nucleic acid binding transcription factor activity” and “antioxidant activity”. Only a few genes belonged to the categories “cell killing”, “electron carrier activity”, “positive regulation of biological process”, “extracellular matrix”, “receptor activity”, “cell junction”, “protein binding transcription factor activity” and “carbon utilization”. The details of GO classification of DTGs in CK vs A, and other comparisons are presented in Figure 2. Plant hormone signal transduction pathways (mediated by either auxin or gibberellin) were well represented, particularly those associated with auxin-mediated signalling. Low temperature sensing and signalling genes influenced by Ca²⁺, as well as other protein kinases were also identified. A number of TF families, genes encoding functional proteins and post-translational regulated genes were represented.

Function annotation of DTGs using the KEGG database

Unigene KEGG annotation was aimed at DTGs from the above comparisons. In the CK vs A comparison, 1,972 DTGs were assigned to the KEGG database involved 122 pathways; for CK vs B1, 175 DTGs were assigned to 61 pathways; for CK vs B2, 302 DTGs were assigned to 75 pathways; for B1 vs B2, 100 DTGs were assigned to 46 pathways; for CK vs C1, 1,915 DTGs were assigned to 120 pathways; for CK vs C2, 2,175 DTGs were assigned to 120 pathways; for A vs C1, 90 DTGs were assigned to 46 pathways; and for A vs C2, 37 DTGs were assigned to 29 pathways. The details of the KEGG classification of the above comparisons are presented in Additional file 16: Table S13. The major pathways identified were “metabolism”, “biosynthesis of secondary metabolites”, “hormone signal transduction”, “plant-pathogen interaction”, “spliceosome”, “phenylpropanoid biosynthesis”, “endocytosis” and “protein processing in the endoplasm”, “ribosome” and “starch and sucrose metabolism”.

Genes involved in the response to low temperature

More DTGs were identified in the CK vs A, CK vs C1 and CK vs C2 comparisons than in the CK vs B1 and CK vs B2. This observation demonstrated that many genes were involved in low temperature acclimation. The DTGs revealed in the comparisons, A vs C1 and A vs C2, were found to be involved in response to freezing treatment, when the plants were exposed to prior CA. However, the DTGs obtained from the CK vs B1 and CK vs B2 comparisons participated in the response mechanism, when the plants didn’t undergo prior CA. Some of all the DTGs belonged to the group of genes involved in cold sensing and signal transduction pathways. Ca²⁺ signaling pathway, which is regarded as an important signal sensing and transduction pathway under stresses, includes many members, such as cation/calcium exchangers (CCX), calcium-binding proteins (CBP), calmodulin-like proteins (CML), CBL-interacting protein kinases (CIPK), calcium-dependent protein kinases (CDPK) and calmodulin-binding receptor-like kinases (CBRLK). The CK vs A comparison included three CCX, two CBP, five CML, three CIPK, and one CBRLK genes; the equivalent order for the CK vs C1 comparison was, respectively, three, one, five, five and one; for CK vs C2, three, zero, four, four and one; for CK vs B1, zero, one, three, zero and zero, along with one sodium/calcium exchanger family protein gene (SCE); for CK vs B2, zero, zero, five, two, zero, along with one SCE and one CDPK; for A vs C1, zero, one, two, one, zero; and for A vs C2, only one CBP was involved. In the comparisons, CK vs A, CK vs C1, and CK vs C2, three CCX (CL8580.Contig_S2_3, Unigene82641_S2_3, and Unigene-21268_S2_3), one CML (Unigene74966_S2_1), and one CDPK (Unigene5836_S2_1) were all found to be up-regulated; in addition, two CMLs (Unigene13572_S2_1 and 3462_S2_3), and one CDPK (CL8813.Contig2_S2_1) show-ed down-regulation. One CML (Unigene76667_S2_3) was found to be up-regulated in both comparisons, CK vs B1 and CK vs B2. These findings indicated that the genes related to Ca²⁺ signaling pathway played an essential role in the phenomenon of CA and freezing response without a prior CA; the complete details of DTGs involved in Ca²⁺ signalling pathway are presented in Table 2. No DTGs associated with the MAPK cascade was identified in the CK vs A comparison. In contrast, one MAPKKK and one MAPK gene were found to be differentially transcribed in CK vs C1, and one MAPKKK was represented in CK vs C2 as well as in CK vs B1. Four MAPKKK genes were observed in the CK vs B2 comparison, and one MAPKKK gene was identified in both, A vs C1 and A vs C2. The results showed that MAPK cascade, especially MAPKKK genes, participated in the plant’s response to freezing treatment, instead of CA. Therefore, all the DTGs involved in MAPK cascade were found to be up-regulated (Table 3). A higher number of genes encoding serine/threonine-protein kinases were differentially transcribed in CK vs A, CK vs C1 and CK vs C2 than in CK vs B1 or CK vs B2 (Table 4). However, no DTG encoding this specific kinase was identified in the A vs C1 or A vs C2 comparisons.

In the present study, members of various low tem-perature-responsive transcription factor (TF) families were identified; and 43, 44 and 46 such genes were found to be differentially transcribed in the comparisons, CK vs A, CK vs C1 and CK vs C2, respectively. The major TF families presented were AP2/ERF, bHLH, WRKY and TCP, along with small numbers of MYB, MYC, NAC, DOF, and the trihelix family. The CK vs B1 and CK vs B2 comparisons showed 8 and 19 TF DTGs, respectively. The number of TF DTGs identified in the CK vs B1 or CK vs B2 comparison was lesser than the comparisons, CK vs A, CK vs C1, or CK vs C2, which were involved in the process of low temperature acclimation. A larger number of TF DTGs were present in CK vs B2 than in the CK vs B1 comparison. In addition, five (3 WRKYs, 1 DREB and 1 ERF) and seven TF DTGs (1 DREB, 1 bHLH and 5 ERFs) were found in the A vs C1 and A vs C2 comparisons, respectively. The identified TF DTGs from the A vs C1 and A vs C2 comparisons, were involved in response to freezing treatment in plants with prior exposure to CA. The differently transcribed TFs of the CKA treatment, the A vs C1, and A vs C2 comparisons are presented in Table 5.

Other relevant classes of protein, which featured as DTGs products, were dehydrin, LEA (late embryogenesis abundant) proteins, heat shock proteins (HSPs). The proteins involved in post-transcriptional regulation, such as ribosomal proteins and a DEAD-box ATP-dependent RNA helicase (RH), were particularly found in the CK vs A, CK vs C1 and CK vs C2 comparisons. Yet, with respect to the above classes of proteins, only one LEA was found in both, A vs C1 and A vs C2, while two HSPs were identified in the A vs C1 comparison. The details of classes of protein of the CKA treatment are presented in Table 6.

Verification of differential transcription using quantitative real time PCR (qPCR)

To further verify the expression profiles of genes in our Illumina RNA-Seq results, we have performed a selection of 15 DTGs for their key roles in response to low temperature by qRT-PCR, these incorporated genes encoding serine/threonine-protein kinase, LEA protein, dehydrin, a gibberellin-regulated protein, a jasmonate ZIM-domain protein, and a DEAD-box ATP-dependent RH, along with a selection of TFs (WRKY, DREB, AP2, bHLH and DOF). The qPCR outcomes in each case correlated closely with the transcript abundances estimated from the RNA-Seq output (Figure 3).

Global patterns of transcription in response to low temperature

The information available on the molecular basis of the response of Chrysanthemum to low temperature is still meagre. However, the development of next-generation sequencing technology provides a straightforward method for the identification of genes involved in this process, and we can try to elaborate the molecular mechanism underlying the response to low temperature. Over 77% of the reads in each of the six RNA-Seq libraries corresponded with known transcripts (Table 1), a proportion which was as high as that achieved in a similar study of Anthurium [30]. The less than 23% of reads were unmapped probably as a result of unidentified transcripts [33]. Around 4,000 DTGs were identified in each of the CKA, CKC1 and CKC2 treatments (Figure 1). As was also the case for Camellia sinensis [28], the majority of DTGs involved up- rather than down-regulation by the stress. In both species, many of the DTGs comprised genes associated with low temperature sensing or signal transduction, low temperature-responsive TFs, stabilization of the plasma membrane and osmosensing-responsiveness.

Low temperature sensing and signaling genes

Low temperature stress-induced signals are directed to various pathways [8]. Ca²⁺ is well recognized as a messenger in stress signalling [34], and is sensed by proteins of three main classes: CDPKs, CaMs and CBLs [35, 36]. A further important pathway is the mitogen-activated protein kinases (MAPKs) cascade [8], while a number of receptor-like protein kinases (RLKs) are known to be responsible for perceiving changes in the external environment and transducing the appropriate signal [37, 38]. The Arabidopsis and rice genomes harbor, respectively, 50 and 32 CML-encoding genes [39]. In rice, the transcription of the calmodulin-like OsMSR2 gene is significantly up-regulated by a series of stresses including low temperature in different tissues at different developmental stages, and its heterologous expression in A. thaliana has suggested that the gene affects salinity and drought tolerance in an ABA-dependent manner [40]. In Camellia sinensis exposed to low temperature, five calmodulin genes, two CDPK genes and one CBL gene were identified to be involved in signal transduction [28]. Here, five Unigenes resembling CML (Unigene74966_S2_1, 13572_S2_1, 3462_S2_3, 76667_S2_3 and 42277_S2_1) were identified as significant DTGs in the CKA treatment, and three of these were also differentially transcribed in both the CKC1 and CKC2 treatments.

CBL-interacting protein kinases (CIPKs), which specifically interact with CBLs, are thought to act as sensors since they lack any enzymatic activity [41]. Three DTGs with homology to CIPK (Unigene5836_S2_1, CL4456.Contig1_S2_3 and CL8813.Contig2_S2_1) were identified in the CKA treatment. Both the CKC1 and CKC2 treatments also featured the same two DTGs (Unigene5836_S2_1 and CL8813.Contig2_S2_1), as well as the other three CIPK homologues (CL2470.Contig1_S2_1, Unigene82472_S2_1, Unigene71605_S2_1) in the former treatment, and two CIPK homologues (CL794.Contig2_S2_1, CL4456.Contig1_S2_3) in the latter. A number of CDPKs have been proven to participate in rapid abiotic stress and immune signaling responses [42]. Transgenic Arabidopsis heterologously expressing the Populus euphratica gene PeCPK10 show an enhanced level of freezing tolerance, perhaps through the transgene's enhancement of the transcript abundance of the abiotic stress-responsive genes RD29B and COR15A [43]. Only one CDPK homologue (CL4878.Contig1_S2_1) was identified as a DTG in the CKB2 treatment, but no other homologues featured as DTGs in any of the other treatments. However, in A vs C1 comparison, with respect to Ca²⁺ signalling pathway, only two CML genes (Unigene42277_S2_1 and 86214_S2_3) and one CIPK gene (CL2470.Contig1_S2_1) were identified as DTGs. Four MAPKKK genes (Unigene34028_S2_3, Unigene36211_S2_1, Unigene8656_S2_1, and Unigene1-412_S2_1) were differently transcribed in the CKB2 treatment. Unigene36211_S2_1 also featured in the CKB1 and CKC2 treatments, while a fifth MAPK gene (Unigene-36598_S2_1) was found to transcribe differently in the CKC1 treatment. No one gene of MAPK cascade was identified in the CKA treatment. In addition, many DTGs encoding serine/threonine-protein kinases were found in the CKA, CKC1 and CKC2 treatments; however, no gene was detected in the comparisons, A vs C1 and A vs C2. These findings provided evidence for the crucial role of Ca²⁺ in the low temperature acclimation process in C. nankingense, and further proved that the MAPK pathway and serine/threonine-protein kinases are more strongly involved in the response to freezing.

Major classes of TF involved in the response to low temperature

Transcriptional regulation of stress-responsive genes is a vital component of the response to both abiotic and biotic stress [44]. Five major TF classes (AP2/ERF, bHLH, WRKY, TCP and MYB) were identified as DTGs in the treatments involved in a process of low temperature acclimation. The AP2/ERF family, which a large group of plant-specific transcription factors, has been sub-divided into AP2, RAV, ERF and DREB TFs [45]. DREBs control the ABA-independent transcription of low temperature responsive genes in A. thaliana [10]; the AtDREB1 sub-family harbors six members [46], of which DREB1A/CBF3, DREB1B/CBF1 and DREB1C/CBF2 are the ones which respond most rapidly to low temperature. A. thaliana plants constitutively expressing any one of these TFs display a heightened tolerance to freezing, drought and salinity [45]. The DREB2B TF present in the desert-adapted plant Eremosparton songoricum has been shown to enhance the tolerance of both yeast and tobacco against a variety of abiotic stresses. The constitutive expression in tobacco of EsDREB2B promotes the accumulation of proline in response to abiotic stress (including low temperature) [47]. Here, the DREB2 homologue Unigene73473_S2_1 was up-regulated upon exposure to low temperature acclimation. It was speculated that the DREB2 may be involved in the accumulation of proline in response to low temperature. It is well known Proline accumulates in many plant species in response to environmental stress [48]. The constitutive expression of an AP2 TF has been shown to improve the tolerance of A. thaliana to low temperature, as well as to drought and high temperature [49]. Here, three AP2-like genes (Unigene27271_S2_3, Unigene27661_S2_3 and CL1514.Contig4_S2_1) were also all up-regulated by CA. A number of other TFs, belonging to the WRKY, bHLH, TCP, MYB, MYC, Trihelix and b-ZIP families were also among the DTGs identified in treatments involving low temperature acclimation. The class of b-ZIP transcription factors, ABRE binding proteins (AREBs or ABFs), can bind to ABRE and activate ABA-dependent gene expression when plants are exposed to low temperature [10]. Many researches have also indicated that three families of transcription factors: WRKY, bHLH and MYB closely related to plant cold stress [50–52]. Five of six WRKYs, which were detected as DTGs, were down-regulated in CK vs A comparison. However, three WRKYs were all up-regulated in A vs C1, and no WRKY as DTG was found in A vs C2. These findings suggested that WRKY family played essential and different roles in CA and early freezing response. Most ERFs also showed down-regulated in CK vs A comparison, while, one ERF was detected up-regulated as DTG in A vs C1; on extending the freezing treatment, five ERFs as DTGs were all up-regulated in A vs C2. The ERFs performing under different temperature treatment showed its roles in CA and cold tolerance. Few reports to date have indicated that TCP transcription factors involve in the response to cold stress.

Low temperature responsive genes related to post-transcriptional regulation

Post-transcriptional regulation (pre-mRNA processing, mRNA stabilization and mRNA export from the nucleus) has been implicated in the process of low temperature acclimation [7]. DEAD-box RHs are intimately associated with RNA-mediated processes and are related to various RNA metabolism events, including RNA synthesis to RNA degradation by means of catalyzing the ATP-dependent unwinding of local RNA secondary structures [53]. The transcription of the genes encoding these proteins is known to be regulated by stress in both bacteria and plants [54–57]. In A. thaliana, RH25 has been associated with enhanced freezing tolerance, probably through its function as an RNA chaperone [58]. The product of RCF1, a low temperature-inducible RH gene, is important for low temperature-responsive gene regulation and low temperature tolerance in plants through maintenance of normal pre-mRNA splicing instead of regulating mRNA export like a previously reported DEAD-box RH (LOS4) regulates mRNA export [59]. It has been demonstrated that the functional roles and RNA chaperone activity related to intron splicing in mitochondrial and chloroplast [59–62]. Here, respectively, eleven, nine and twelve genes encoding DEAD-box RHs were up-regulated in the CKA, CKC1 and CKC2 treatments, but none were differentially transcribed in either the CKB1 or the CKB2 treatment; this was taken to imply that DEAD-box RHs are activated during the process of low temperature acclimation. In future, a major task is to clear how RNA chaperones recognize substrate RNAs and how they work with other proteins to regulate post-transcriptional RNA metabolism in response to developmental and environmental condition [53].

Genes encoding functional proteins

Besides protein factors involved in further regulation of signal transduction and gene expression such as transcription factor and protein kinase that probably function in stress response, various functional proteins featured as products of the DTGs, in particular, LEA protein and dehydrin, LEA proteins are accumulated during the late stage of seed maturation and under moisture deficient conditions, and act to protect higher plants from damage caused by abiotic stress. When the maize gene ZmLEA3 is expressed in tobacco and yeast, it improves the plant's level of tolerance against both osmotic and oxidative stress [63]. The wheat LEA protein gene WCI16 expressed heterologously in A. thaliana enhances freezing tolerance [64]. Here, seven LEA protein genes were among the DTGs identified in the CKA treatment, five in both the CKC1 and CKC2 treatment, none were present in either the CKB1 or the CKB2 treatment, and only one was identified in both comparisons of A vs C1 and A vs C2. The implication was that LEA proteins probably enhance low temperature tolerance through their participation in the low temperature acclimation process. Dehydrins constitute a group of plant proteins involved in tolerance to low temperature and drought [65]. The thranscrips of genes encoding three dehydrins in E. globulus accumulate strongly in the stem and leaf tissue of acclimated plants, compared to non-acclimated [66]. Here, four genes encoding dehydrins were identified as DTGs in the CKA, CKC1 and CKC2 treatments. But none in either the CKB1 or CKB2 treatments, suggesting that, as with the LEA protein genes, their contribution to low temperature tolerance is expressed during the low temperature acclimation process. HSPs were initially identified from their involvement in the response to high temperatures, but it is now recognized that many of them also respond to low temperature. It is well known HSPs can protect plants against stress by means of reestablishing normal protein conformation and thus cellular homeostasis. Five major HSP families have been conservely defined, based on their molecular weight: HSP100s, HSP90s, HSP70s, HSP60s and small HSPs. Five HSP90-encoding DTGs (all up-regulated) featured in the CKA, CKC1 and CKC2 treatments. HSP90s function as molecular chaperones during signal transduction, cell cycling, the stress response, and in protein folding, degradation and transport [67, 68]. In A. thaliana, expression of Hsp90 is developmentally regulated, but it also responds to high and low temperature, as well as salinity [67]. It is speculated that Hsp90 as molecular chaperones play an important role in signal transduction and stress management in C. nankinginse. A gene HSP70- encoding DTGs (all down-regulated) featured in the CKA, CKB2, CKC1 and CKC2 treatments, which not found in the CKB1 treatment; and a number of small HSPs, belonging to the 15.7, 22.7 and 23.6 families were also among the DTGs identified in treatments involving low temperature response.

In conclusion, it was clear that a large number of genes were induced when exposed to low temperature. The complex gene networks involved a set of interactions between cold sensing and signaling genes, genes related to post-transcription regulation responding to cold stress, transcription factors and certain functional proteins. In the present study, the DTGs identified as candidate genes in response to low temperature, require further investigation for a complete understanding of the molecular basis of cold response in C. nankingense. This information will prove beneficial in molecular breeding programs for excellent chrysanthemum varieties.

Plant material

The accession of C.nankingense utilized is conserved by the Chrysanthemum Germplasm Resource Preserving Centre, Nanjing Agricultural University, China. Plants raised from tissue-cultured plantlet were grown on MS medium (16 h photoperiod, 22°C/18°C day/night temperature, 70% relative humidity). Four week old plants were subjected to one of the following temperature treatments: (A) 4°C for one week, (B1) -5°C for 1 h, (B2) -5°C for 2 h, (C1) 4°C for one week, followed by -5°C for 1 h, (C2) 4°C for one week, followed by -5°C for 2 h. CK plants were harvested without any additional treatment. For each treatment, the leaves and stems of three seedlings were sampled, and two samples were harvested. 30 μg of total RNA were pooled in equal amounts from the two biological replicates for subsequent RNA- Seq.

RNA isolation and cDNA library construction

Six separate libraries were prepared. The samples from the six treatments (A, B1, B2, C1, C2 and CK) were snap-frozen in liquid nitrogen and ground to a fine powder. Total RNA was extracted using a Total RNA Isolation System (Takara, Japan) according to the manufacturer’s instructions. When the quality of the resulting RNA was verified using a 2100 Bioanalyzer RNA Nano chip device (Agilent, Santa Clara, CA, USA), all six extractions delivered an RNA integrity number value of > 8.0, and a 28S:18S ratio >1.5. After checking for the absence of contamination by protein (A260/A280 nm ratios) and reagent contamination (A260/A230 nm ratios) by a Nanodrop ND-430 1000 spectrophotometer, the extractions were selected based on 28S/18S rRNA band intensity (1.5:1 ~ 2:1) and spectroscopic A260/A280 nm readings between 1.8 and 2.2, A260/A230 nm readings greater than 2.0. 10 μg RNA was pooled from each of the three sampled plants. The total RNA preparation was treated with RNase-free DNase I (Takara, Japan) to degrade any possible DNA, and mixed with oligo (dT) coated magnetic beads to concentrate the polyA mRNA.

The mRNA was fragmented into short fragments ~200 nt pieces by incubation in a fragmentation buffer under elevated temperature (The Beijing Genomics Institute). The first strand of cDNA was then synthesized by priming with random hexamer, and the second strand was generated with buffer, dNTPs, RNase H and DNA polymerase I. The double strand cDNA was purified with a QiaQuick PCR extraction kit and resolved with EB buffer for end repair and addition of single nucleotide A. Finally, sequencing adaptors were ligated to the fragments. Following agarose gel electrophore, suitable fragments were selected as templates for PCR.

RNA-Seq

The library was sequenced using an Illumina HiSeq^™2000 located at the Beijing Genomics Institute (Shenzhen, China; http://www.genomics.cn/index). The data were deposited in the US National Center for Biotechnology Information (NCBI) Sequence Read Archive (SRA, http://www.ncbi.nlm.nih.gov/Traces/sra; [69] under accession number (SRP041138). Raw data were saved as .fastq files. Data filtering is performed to obtain “clean reads” for further analysis. Clean reads were obtained by removing adaptor sequences, reads in which the percentage of unknown bases (N) was greater than 10% and low quality reads. The clean reads were mapped onto the reference sequences using SOAP (2.21) software [70]. A maximum of two mismatches was allowed in the alignment. The NCBI non-redundant protein (Nr) database (http://www.ncbi.nlm.nih.gov) and the Swiss-Prot protein database (http://www.expasy.ch/sprot) were used for blast search and annotation using an E-value cut off of 10^-5. Functional annotation by gene ontology terms (GO, http://www.geneontology.org) was analyzed using the Blast2GO program [71]. The Kyoto Encyclopedia of Genes and Genomes Pathway (KEGG; http://www.genome.jp/kegg), the major public pathway-related database [72] was also used to predict and classify possible functions. The RPKM reads (clean reads per kilo base per million) method [31] was used to estimate transcript abundance on the base of eliminating the influence of different gene length and sequencing discrepancy. Therefore, the RPKM values can be directly used for comparing the difference of gene expression among samples.

Identification of differentially expressed genes

To compare the differences in gene expression, the method of an algorithm developed by Audic et al. [32] was used to identify DTGs. The criteria applied were an FDR (false discovery rate) less than 0.01 and an absolute value of log2 ratio of at least 1. Then, the DTGs were subjected to GO and KEGG Ontology (KO) enrichment analysis on the base of a hypergeometric test.

qPCR validation of differential transcription

Total RNA was isolated from leave and stem of plants subjected to the various treatments described above. Contaminating DNA was removed by treating with RNase-free DNase I and the first cDNA strand was synthesized from 1 μg total RNA using PrimeScript® Reverse Transcriptase (Takara, Dalian, China) and an oligo (dT) primer, according to the manufacturer’s instructions. qPCRs were performed in an Eppendorf Real Time PCR System (Mastercycler®ep realplex, Germany) using a SYBR Premix Ex Taq™ Kit (Takara), according to the manufacturer’s protocol. Gene-specific primers were designed using Primer5 software (sequences given in Table 7). Each 20 μL qPCR contained 5 μL diluted cDNA, 100 nM of each primer, and 10 μl SYBR Green PCR master mix, and was exposed to an initial denaturation (95°C/2 min), followed by 40 cycles of 95°C/15 s, 60°C/15 s, 72°C/15 s. After amplication, all results were screened to verify a single peak melting curve for the specificity of the amplifications. Three biological replicates were performed for each sample. Relative transcript abundance was obtained by including the C. nankingense EF1α gene as the reference, and was based on the 2^-ΔΔCT method [73].