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

Background A major constraint affecting the quality and productivity of chrysanthemum is the unusual period of low temperature occurring during early spring, late autumn, and winter. Yet, there has been no systematic investigation on the genes underlying the response to low temperature in chrysanthemum. Herein, we used RNA-Seq platform to characterize the transcriptomic response to low temperature by comparing different transcriptome of Chrysanthemum nankingense plants and subjecting them to a period of sub-zero temperature, with or without a prior low temperature acclimation. Results Six separate RNA-Seq libraries were generated from the RNA samples of leaves and stems from six different temperature treatments, including one cold acclimation (CA), two freezing treatments without prior CA, two freezing treatments with prior CA and the control. At least seven million clean reads were obtained from each library. Over 77% of the reads could be mapped to sets of C. nankingense unigenes established previously. The differentially transcribed genes (DTGs) were identified as low temperature sensing and signalling genes, transcription factors, functional proteins associated with the abiotic response, and low temperature-responsive genes involved in post-transcriptional regulation. The differential transcription of 15 DTGs was validated using quantitative RT-PCR. Conclusions The large number of DTGs identified in this study, confirmed the complexity of the regulatory machinery involved in the processes of low temperature acclimation and low temperature/freezing tolerance. Electronic supplementary material The online version of this article (doi:10.1186/1471-2164-15-844) contains supplementary material, which is available to authorized users.


Background
Chrysanthemum (Chrysanthemum morifolium) is a popular ornamental plant worldwide [1,2]. Chrysanthemum plants are susceptible to damage when exposed to prolonged periods of low temperature; therefore, improving their tolerance to cold stress is perceived as an important breeding goal. The high chromosome number and polyploidy of major ornamental species complicate the genetics and capacity for gene discovery [3]. Hence, C. nankingense has been considered as a convenient genomic model due to its simple diploid nature. In addition, it displays better tolerance to low temperature as compared to the ornamental polyploid species [4].
Temperature is a major determinant of the geographical distribution and length of the growing season in most plant species [5,6]. However, episodes of low temperature during the growing season cause a substantial loss in the yield of many temperate crops. During chrysanthemum production, specifically in China, extreme low temperatures in early spring and winter, unusual freezing temperatures during late cold spring, and sudden frosts in fall often lead to growth arrest and block flower buds or inflorescence, which in turn result in significant economic losses every year [7]. Temperate plant species can acquire the ability to withstand a prolonged period of sub-zero temperature if they are previously exposed to a period of low temperature above 0°C; this phenomenon is known as low temperature acclimation [5,8,9]. However, the molecular basis of cold acclimation (CA) and low temperature/freezing tolerance in chrysanthemum has not yet been explored. In the present study, we aim to fish out candidate genes underlying the process of CA and response to low /freezing temperature, which will help to elucidate the molecular basis of the cold response in C. nankingense, and then improve the chrysanthemum varieties cold tolerance.
In Arabidopsis thaliana and some of the winter cereals, a variety of physiological, biochemical, and molecular changes are known to occur during the low temperature acclimation process [9]. In these species, the physical state of the plasma membrane has been shown to be an important determinant of the plant's ability to sense changes in the air temperature [10][11][12][13][14]. Membrane rigidification leads to an increase in the cytosolic concentration of the Ca 2+ ion [15], which is regarded as a major regulator of low temperature responsive factor. Certain Ca 2+ -dependent protein kinases have also been recognized as positive regulators [16]. The mitogen-activated protein kinase (MAPK) cascade participates in the low temperature signalling and low temperature tolerance [17]. Three candidate chilling response genes, encoding MAPKKK (MAPK kinase kinase), CLC-D (chloride channel D) and RLK (receptor-like protein kinase) homologues are all up-regulated following chilling stress in Maize [18]. The CBF low temperatureresponse pathway is well conserved across a diversity of species and is a significant component of tolerance to sub-zero temperatures [7,19]. The CBF pathway is positively regulated by the circadian clock components CCA1 and LHY [20], while the INDUCER OF CBF EX-PRESSION 1 (ICE1) protein acts upstream of CBF in the low temperature response pathway [21]. A high expression of osmotically responsive gene (HOS) 1 acts as a negative regulator of the low temperature response [21,22]. Salicylic acid (SA) is also involved in the response to low temperature stress [23,24]. RNA processing and nucleocytoplasmic transport play crucial roles in plant stress [25].
The RNA-Seq platform has become a key technology for quantifying the transcriptional response in nonmodel organisms or those with genome characteristics extremely difficult to whole-genome sequencing [26,27]. In tea (Camellia sinensis), this approach has enabled the recognition of 1,770 differentially transcribed genes (DTGs) induced during low temperature acclimation [28]. A similar transcriptomic analysis of Jatropha curcas identified over 3,000 genes as being up-or down-regulated by low temperature [29]. While in Anthurium sp., the method of digital gene expression enabled 39 low temperature-inducible transcription factors (TFs) to be identified [30].
Here, the RNA-Seq platform based on Illumina NGS technology was used to characterize the transcriptomic response to low temperature by comparing the different transcriptome of C. nankingense plants subjected to periods of sub-zero temperature with or without a prior low temperature acclimation, with a view to gaining a deeper insight into the molecular basis of this physiological adaptation. In addition, the identified candidate genes will be useful for improving adaptation to low temperature and enhancing productivity and geographical distribution.

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;  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  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 Figure 1 The numbers of DTGs identified in comparisons between pairs of libraries.

GO classification of differentially transcribed genes
In  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 2+ , as well as other protein kinases were also identified. A number of TF families, genes encoding functional proteins and posttranslational regulated genes were represented.  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 2+ 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), CBLinteracting protein kinases (CIPK), calcium-dependent protein kinases (CDPK) and calmodulin-binding receptorlike kinases (CBRLK  (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 temperature-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   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).

Discussion
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 uprather than down-regulation by the stress. In both species,

Low temperature sensing and signaling genes
Low temperature stress-induced signals are directed to various pathways [8]. Ca 2+ 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.Con-tig1_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.Con-tig1_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 2+ 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, Uni-gene36211_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 2+ 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 subfamily 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 desertadapted 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 AP2like 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][51][52]   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][55][56][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][60][61][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.

Conclusions
An overview of the many changes to the C. nankingense transcriptome induced by exposure to low temperature has been provided. Many of the DTGs identified were involved in predictable classes of gene, encoding, for example, low temperature sensors and signaling molecules, TFs and certain functional proteins associated with the stress response. The large number of DTGs identified confirms the complexity of the regulatory machinery involved in low temperature acclimation and low temperature/freezing tolerance. Establishing which of these many DTGs reflect the primary response to low temperature is necessary before the molecular basis of the response can be fully elaborated; such genes would represent candidates for intervention in the breeding of chrysanthemum varieties endowed with a heightened tolerance to low temperature stress.

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 nonredundant 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 RNasefree 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 Table 7 Primers of quantitative reverse transcription-polymerase chain reaction for validation of RNA-Seq data  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].

Additional files
Additional file 1: Figure S1. Composition of raw reads in the six RNA libraries. "Clean" reads refers to those remaining after the removal of adaptor sequences, reads in which the proportion of missing bases was >10% and reads in which low quality (≤5) bases represented >50% of the reads. The numbers in parentheses indicate the percentage of each type of read present.
Additional file 2: Figure S2. Sequencing saturation analysis in the six libraries (A, B1, B2, C1, C2 and CK). The numbers of new genes detected rose as the read number was increased, but not beyond a threshold around 7,000,000.
Additional file 3: Figure S3. Distribution of gene coverage in the six libraries.
Additional file 4: Table S1. The transcription level of each unigene derived from the number of relevant reads recovered in the four libraries.
The "GeneLength" column gives the length of exon sequence.
Additional file 5: Table S2. Genes differentially transcribed in the comparison between libraries CK and A. The criteria applied for assigning significance were: P-value < 0.05, FDR ≤ 0.001, and estimated absolute |log2 Ratio Additional file 9: Table S6. Genes differentially transcribed in the comparison between libraries CK and C2. The criteria applied for assigning significance were: P-value < 0.05, FDR ≤ 0.001, and estimated absolute |log2 Ratio(C2/CK) | ≥1. Genes listed in descending order of absolute | log2 Ratio(C2/CK) |. GeneIDs retrieved from the Chrysanthemum nankingense Reference Sequence Database. Annotation of unigene sequences performed using BlastX (E <10). The "GeneLength" column gives the length of exon sequence. CK-and C2-expression: frequency of unigene transcripts in libraries CK and C2, respectively. CK-and C2-RPKM: reads per kb per million reads for each unigene in libraries CK and C2, respectively. Log2 Ratio (C2/CK) : the ratio between the RPKM in CK and the RPKM in C2. Up-Down-Regulation (C2/CK), P-value and FDR of each gene are also shown. KEGG: annotation according to the KEGG database by BLAST. Blast nr: identification of homologues in GenBank. GO Component, GO Function and Go Process: ontology information of Cellular Components, Molecular Function and Biological Processes of Gene-corresponding GO terms. "-": no hit.
Additional file 10: Table S7. Genes differentially transcribed in the comparison between libraries A and C1. The criteria applied for assigning significance were: P-value < 0.05, FDR ≤ 0.001, and estimated absolute |log2 Ratio(C1/A) | ≥1. Genes listed in descending order of absolute |log2 Ratio(C1/A) |. GeneIDs retrieved from the Chrysanthemum nankingense Reference Sequence Database. Annotation of unigene sequences performed using BlastX (E <10). The "GeneLength" column gives the length of exon sequence. A-and C1-expression: frequency of unigene transcripts in libraries A and C1, respectively. A-and C1-RPKM: reads per kb per million reads for each unigene in libraries A and C1, respectively. Log2 Ratio(C1/A) : the ratio between the RPKM in A and the RPKM in C1. Up-Down-Regulation (C1/A), P-value and FDR of each gene are also shown. KEGG: annotation according to the KEGG database by BLAST. Blast nr: identification of homologues in GenBank. GO Component, GO Function and Go Process: ontology information of Cellular Components, Molecular Function and Biological Processes of Gene-corresponding GO terms. "-": no hit.