GENE-39939; No. of pages: 8; 4C: Gene xxx (2014) xxx–xxx

Contents lists available at ScienceDirect

Gene journal homepage: www.elsevier.com/locate/gene

Methods paper

Identification of differentially expressed genes in Chrysanthemum nankingense (Asteraceae) under heat stress by RNA Seq Jing Sun, Liping Ren, Yue Cheng, Jiaojiao Gao, Bin Dong, Sumei Chen, Fadi Chen ⁎, Jiafu Jiang ⁎ College of Horticulture, Nanjing Agricultural University, Nanjing 210095, China

a r t i c l e

i n f o

Article history: Received 10 July 2014 Received in revised form 14 August 2014 Accepted 4 September 2014 Available online xxxx Keywords: Chrysanthemum nankingense Heat stress RNA-Seq Stress response

a b s t r a c t The RNA-Seq platform was used to characterize the high-temperature stress response of Chrysanthemum nankingense. A set of 54,668 differentially expressed unigenes was identified. After a threshold of ratio change ≥2 and a q-value of b 0.05 were applied, the number of differentially transcribed genes was reduced to 3955, of which 765 were up-regulated and 3190 were down-regulated in response to heat stress. The differentially transcribed genes were predicted to participate in 26 biological processes, 4 cellular components, and 13 molecular functions. Among the most differentially expressed genes between the two libraries were well-recognized high-temperature responsive protein families, such as heat shock factors and heat shock proteins, various transcription factor families, and a number of RNA metabolism-related genes. Overall, the RNA-Seq analyses revealed a high degree of transcriptional complexity in early heat stress response. Some of these high-temperature responsive C. nankingense genes may prove useful in efforts to improve thermotolerance of commercial chrysanthemum. © 2014 Elsevier B.V. All rights reserved.

1. Introduction Global climate change is becoming a serious threat to the environment as well as human, plant and animal populations, with the frequency of extreme weather conditions, flooding, and typhoons on the increase in the last few decades (Parmesan and Yohe, 2003). Exposure to extremely high temperature results in heat stress, which causes disruption of many biological processes. Plant response to heat stress involves three distinct phases: triggering, maintenance and recovery (Baniwal et al., 2004; Morimoto, 1998). In the model plant Arabidopsis thaliana, the genetic basis of each of these responses has been shown to be polygenic, with the transcription level of a large number of genes being affected; at the phenotypic level, a range of physiological and biochemical responses has been identified, including altered patterns of photosynthesis (Zonia and Munnik, 2006), respiration (Davidson and Schiestl, 2001), hormone production (Larkindale and Knight, 2002; Teale et al., 2006) and antioxidation pathway (Koussevitzky et al., 2008). In rice too, prolonged exposure to high temperature has been shown to affect a wide spectrum of

Abbreviations: FDR, false discovery rate; GO, Gene Ontology; DTG, different transcribed genes; KEGG, Kyoto Encyclopedia of Genes and Genomes; PCR, polymerase chain reaction; qPCR, real-time quantitative PCR; HSE, heat shock element; SR, serine/arginine-rich; snRNP, small nuclear ribonucleoprotein; tRNA, transfer ribonucleic acid; MS medium, Murashige and Skoog medium. ⁎ Corresponding authors. E-mail addresses: [email protected] (J. Sun), [email protected] (L. Ren), [email protected] (Y. Cheng), [email protected] (J. Gao), [email protected] (B. Dong), [email protected] (S. Chen), [email protected] (F. Chen), [email protected] (J. Jiang).

morphological traits (Jagadish et al., 2010). In Guan's study, several rice class I small heat shock protein genes were found to have relatively faster response kinetics than other genes from the same family (Guan et al., 2004). The diploid species Chrysanthemum nankingense is closely related to the commercially important polyploid ornamental species Chrysanthemum morifolium (garden chrysanthemum) (Cheng et al., 2010; Yang et al., 2006). It therefore represents both a convenient genomic model for the crop, as well as harboring potentially beneficial genes in the context of chrysanthemum improvement, particularly in the area of abiotic stress tolerance (Liu et al., 2011; Yang et al., 2005). The species has not as yet, however, been extensively studied at molecular level. Here, we report an RNA-Seq based survey of gene transcription in C. nankingense, focusing on its early response to high temperature stress. RNA-Seq is a next-generation high-throughput screening technology that can be used to identify differentially expressed genes under specific biological processes (it is an upgraded version of digital gene expression [DGE]). Through a reference sequence database or characterized genome of related organisms, the identified genes can be annotated for further research (Hegedűs et al., 2009). To investigate the genetic mechanisms regulating early heat response of C. nankingense, we choose the Illumina HiSeq™ 2000 sequencing platform to perform RNA-Seq (Quantification). 3955 genes were found to be up- or downregulated under heat stress; these genes primarily comprise the response network of heat stress. Since C. morifolium has a large and complex genome of approximately 9.4 Gb (http://data.kew.org/cvalues/), its genome has not yet been sequenced. However, C. nankingense shares a close genetic relationship with C. morifolium, so this study may provide

http://dx.doi.org/10.1016/j.gene.2014.09.013 0378-1119/© 2014 Elsevier B.V. All rights reserved.

Please cite this article as: Sun, J., et al., Identification of differentially expressed genes in Chrysanthemum nankingense (Asteraceae) under heat stress by RNA Seq, Gene (2014), http://dx.doi.org/10.1016/j.gene.2014.09.013

2

J. Sun et al. / Gene xxx (2014) xxx–xxx

a useful molecular biology reference not only for C. nankingense but also for further studies on other chrysanthemum species. 2. Results 2.1. Analysis of RNA-seq data The main characteristics of two RNA-seq libraries were summarized in Table 1 and Fig. 1. The number of raw reads per library ranged from ~7.23 to ~7.25 million, and the total number of base pairs sequenced ranged from 354,456,886 to 355,484,220 (Accession No. for library SRP032828, http://www.ncbi.nlm.nih.gov/). Each raw read was sequenced from one end and the length was 50 bp. After filtering low quality reads containing adaptor and N-terminal sequence, and removing low quality sequence (i.e., reads in which b1% of the bases were uncertain), we obtained 7,233,814 and 7,254,780 clean reads, which corresponded to 354,456,886 and 355,484,220 base pairs respectively. The proportion of clean reads was N99.29% in each library (Fig. 1). Transcriptomes from both control (non-stressed) and high-temperature samples were represented by at least seven million clean reads each, which allowed for quantitative analysis of gene transcription. 2.2. Read mapping The chrysanthemum transcriptome database was used as reference. For quality evaluation of the reference gene database, the length of sequences assembled was used as a criterion for assembly success. We calculated the distribution of length of the contigs and unigenes of reference C. morifolium transcriptomes (Fig. 2). Through de novo assembly, 51,377,444 total clean reads of reference database were assembled into 196,944 contigs and 88,277 unigenes. Then, Blastx alignment (evalue b 0.00001) between unigenes/contigs and protein databases such as nr, Swiss-Prot, KEGG and COG was performed, and the results with the best alignment were used to determine sequence direction of unigenes. After directly mapping clean reads from the two C. nankingense libraries into reference unigenes, 60.01% and 64.38% of the clean reads of both library sequences were mapped to the C. morifolium transcriptome database (Table 1). Sequencing quality can be reflected by sequencing saturation data to some extent. With the increase in the amount of mapped sequence data (number of reads), the number of mapped clean reads detected also rose. The number of identified genes increased with the number of clean reads. As shown in Fig. 3a and b, when the amount of clean reads reached 6 million, the curve started to gradually level out and stabilized. The distribution of unigene coverage (i.e., the number of clean reads aligned to reference unigene sequences) is shown in Fig. 3c and d. In all, 3955 filtered unigenes (FDR ≤ 0.001 and log2Ratio ≥ 1) of varying overall length showed evidence of differential transcription (Fig. 4). Reads that were obtained at least 60,000 times were evenly distributed across the matching unigene sequence (Fig. 3e and f). The proportion of clean reads that showed a perfect match with the chrysanthemum transcriptome was 36.25% for those obtained from non-stressed (control) plants and 37.46% for those obtained from the high-temperature samples (Table 1).

2.3. GO and KEGG classification In all, the abundance of 54,668 transcripts was altered as a result of high-temperature stress, but this number was reduced to 3955 after the chosen threshold of P b 0.001 and the absolute value of log2 abundance ratio of ≥1 were applied, based on an FDR of ≤0.001 (Fig. 4b). A higher number of genes were down-regulated rather than up-regulated in response to the stress. The differentially transcribed genes (DTGs) were classified according to “biological process,” “cellular component” and “molecular function” and were distributed among 43 GO categories. Early responding DTGs reflected 26 biological processes (Fig. 5), including metabolism, transport, catabolism, response to stress and stimulus, biosynthesis, regulation and cell communication. Normal developmental processes were strongly affected by the high-temperature stress treatment. KEGG analysis identified 31 DTGs encoding photosynthesisantenna proteins; 210 involved in plant–pathogen interaction; 64 encoding the synthesis of stilbenoid, diarylheptanoid and gingerol; 55 in the regulation of flavonoid synthesis; 153 in hormone metabolism and hormone signaling; and 31 in flavone and flavonoid synthesis. 2.4. High temperature stress induced DTGs When ambient temperature rose from 22 °C to 40 °C, a complex response network of C. nankingense genes appeared. The transcription data indicated that heat stress regulated a number of transcription factors and genes associated with heat stress response. During the early stress response, 22 genes encoding heat shock proteins (Hsps) and five encoding heat shock factors (Hsfs) were induced (Table S1). The DTGs also included a number of genes encoding transcription factors AP2/EREBP, MYB, GRAS, WRKY, bHLH and bZIP (Table S2). The AP2/ EREBP genes were all down-regulated by the high temperature treatment, while two DREB members were markedly up-regulated. Among the ten differentially transcribed MYB family genes, some were upregulated whereas others were down-regulated. Similarly, the transcript abundance of eleven WRKY, five bHLH and two bZIP genes was reduced by the stress, as was that of most of the GRAS family genes. A total of 172 genes encoding various kinases (e.g., Ca-dependent, CBLinteracting, serine/threonine- and leucine-rich repeat receptor-like protein kinase) were involved in the stress response (Table S3). Eleven DTGs encoding peroxidases were all down-regulated, although to varying extents (Table S4). In addition, a number of genes involved in brassinosteroid, auxin, gibberellin and ethylene synthesis/signaling pathways were affected by exposure to high temperature (Table S5). Surprisingly, six TIR-NBS-LRR genes (which are generally associated with defense against pathogens) were repressed by two- to ten-fold (Table S6). Another set of early-responding DTGs featured genes involved in RNA metabolism, such as those encoding RNA polymerase, a pre-mRNA splicing factor, RNA helicase, RNA methyltransferase, and U1 snRNP auxiliary factor (Table S7). Finally, among the genes critically associated with photosynthesis, those encoding PSI, PSII and RuBisCO were all affected at transcriptional level by the high-temperature treatment (Table S8). 2.5. Independent qRT-PCR validation of RNA-Seq-identified DTGs

Table 1 Overview of the RNA-Seq reads acquired from non-stressed and high temperature stressed C. nankingense plants. Summary

Control

Heat stress

Total reads Total base pairs Total mapped reads Perfect match ≤2 bp mismatch Unique match Multi-position match Total unmapped reads

7,233,814 (100.00%) 354,456,886 (100.00%) 4,341,334 (60.01%) 2,622,600 (36.25%) 1,718,734 (23.76%) 2,944,682 (40.71%) 1,396,652 (19.31%) 2,892,480 (39.99%)

7,254,780 (100.00%) 355,484,220 (100.00%) 4,670,833 (64.38%) 2,717,381 (37.46%) 1,953,452 (26.93%) 3,199,140 (44.10%) 1,471,693 (20.29%) 2,583,947 (35.62%)

To test and verify the credibility of RNA-Seq results, independent experiment was performed to validate the RNA-seq results through qPCR. RNA samples at different time points (0 h, 0.5 h, 1 h, 2 h, 4 h and 6 h) were reverse-transcribed and tested. The 2−ΔΔCT method was used for analysis of relative gene expression data (Livak and Schmittgen, 2001). Eight DTGs were randomly chosen including induced and repressed genes. The expression of eight DTGs at 2 h was corroborated by the RNA-Seq-based regulation in every case (Fig. 6). During six constitutive hours of heat shock (hs) treatment, CnHsps and CnHsfs were the most up-regulated genes (Fig. 6a–c). Their expression increased dramatically at the first half hour. It's worth noting that all the members

Please cite this article as: Sun, J., et al., Identification of differentially expressed genes in Chrysanthemum nankingense (Asteraceae) under heat stress by RNA Seq, Gene (2014), http://dx.doi.org/10.1016/j.gene.2014.09.013

J. Sun et al. / Gene xxx (2014) xxx–xxx

3

Fig. 1. Summary of the quality of the RNA-Seq reads. RNA extracted from (a) non-stressed and (b) high temperature stressed C. nankingense plants. After filtering containing adaptor reads, containing N-end reads and low quality reads. The proportion of clean reads was N99.29% in each library.

of the WRKY family were down-regulated, for example CnWRKY (Unigene10297) (Fig. 6e). 3. Discussion A variety of response mechanisms have been developed by plants to adapt to abiotic stresses. Plants exposed to excessively high temperatures commonly react by inducing a number of stress response networks, a process that involves the action and interaction of many transcription factors. A feature of RNA-Seq is that it allows for a genome-wide, quantitative view of how a changing environment can affect the pattern of transcription. Thus, in C. nankingense, as in other plant species, the outcome of high temperature stress is the induction of some genes and the repression of others. 3.1. Heat shock proteins and heat shock factors were mostly up-regulated by heat stress Hsps play a central role in the response of plants to high temperature stress and their acquisition of thermotolerance. Since the original identification of Hsps in Drosophila melanogaster (Tissieres et al., 1974), this class of protein has been identified in many plants. These proteins sometimes act as stress responders, but they act as molecular chaperones in other cases (Ellis, 2000). Hsp70 family is a key player involved in the detection of accumulation of denatured proteins in stressed cells (Nover and Scharf, 1997). They are some of the most abundant hs proteins that function in an ATP-dependent manner in eukaryotic cells. This protein plays a role in the folding of nascent polypeptides released from ribosome (Hartl and Hayer-Hartl, 2002), sorting of proteins to cell organelles (Zhang and Glaser, 2002), and the formation of a bridge to ubiquitin-mediated proteasomal degradation pathway (Lee

et al., 2009). In C. nankingense, the up-regulation of CnHsp70 (Unigene20443) was accompanied by a rapid increase in the expression of CnHsp90 (CL3198.Contig3). A model of the CnHsp and CnHsf network under heat stress was proposed based on the data of present study (Fig. 7). In mammals, the interaction of Hsp90, Hsp70, Hsp27 and Hsp40 plays an important role in proteasomal degradation (Mehlen et al., 1996). Hsp90, Hsp70 and Hdj-1 also play a role in the proper folding of unfolded substrates into native proteins (Morimoto and Santoro, 1998). Hsp101 is a key component in the acquisition of thermotolerance in plants. Over-expression of Hsp101 in basmati rice could significantly improve its heat-tolerance (Katiyar-Agarwal et al., 2003). After heat shock, CnHsp101 (Unigene17914) and small CnHsp (Unigene15402) were upregulated, but to a different extent. Hsfs are pivotal components of signal transduction process and mediate the transcription of Hsp genes. In keeping with this, in this study too, we found that a specific set of Hsfs is also essential for thermotolerance (Fig. 7). The activity of HsfA1a, the master regulator of heat stress response in tomato, is repressed by Hsp70 in cytoplasm (Hahn et al., 2011). Through oligomerization with Hsp90, the HsfA1a protein is activated and binds to heat shock elements (HSE), which are present in the promoters of hs genes. This subsequently results in the expression of a series of hsrelated genes, the main ones being Hsp70, Hsp90, HsfA1a, HsfA2 and HsfB1 (Baniwal et al., 2004). However, the regulation of HsfA1a and HsfB1 was not detected in C. nankingense samples exposed to heat stress, which indicates that the up-regulation of HsfA2 and HsfA3 may play more important roles in the heat stress pathway in this species (Fig. 7). In Arabidopsis, 21 Hsfs were first systematically classified into three classes, A, B and C, based on the presence of conserved DNAbinding domains (DBDs) plus the adjacent HR-A/B region amino acid structure (Nover et al., 2001). In tomato, HsfA1a, HsfA2 and HsfB1 seem to form a regulatory network that is responsible for the expression

Fig. 2. Variation in the length of unigenes and contigs in reference gene database. Through de novo assembly, total clean reads of reference database were assembled into 88,277 unigenes (a) and 196,944 contigs (b).

Please cite this article as: Sun, J., et al., Identification of differentially expressed genes in Chrysanthemum nankingense (Asteraceae) under heat stress by RNA Seq, Gene (2014), http://dx.doi.org/10.1016/j.gene.2014.09.013

4

J. Sun et al. / Gene xxx (2014) xxx–xxx

Fig. 3. RNA-Seq analysis of each library. (a–b) Sequencing saturation of the two libraries. The number of different genes detected rose as the read number was increased; (c–d) Distribution of gene coverage analysis of each library. C: control and H: heat stress. Since reference genes have different lengths, the mapped reads have different coverage. The number of reads coverage is counted in ratio. (e–f) Randomness assessment of reads. We use the distribution of reads on the reference genes to evaluate the randomness. When the sequencing counts reached about 60,000 reads, the reads in every position were evenly distributed. C: control and H: heat stress.

of downstream hs-responsive genes (Baniwal et al., 2004; Charng et al., 2007) (Fig. 7). HsfA2 is tightly integrated within a network of interacting proteins, while HsfB1 acts mainly to co-regulate the activity of HsfA1a (Mishra et al., 2002). In HsfA2-overexpressing Arabidopsis plants, a large number of hs proteins, such as ascorbate peroxidase 2, were highly expressed compared with the wild-type plants (Nishizawa et al., 2006; Sakuma et al., 2006). In rice, expression of OsHsfA2e enhanced the tolerance of transgenic plants to environmental stresses (Yokotani et al., 2008). In this study, CnHsfA2 was shown to be a thermo-inducible gene during the early stage of heat stress (Fig. 6). Therefore, CnHsfA2 was speculated as a key heat stress response factor in C. nankingense. HsfA3 was later on proved to be regulated by DREB2A under heat stress (Schramm et al., 2008). Consistently, in this study, expression of the CnDREB2 gene family (Unigene26932, CL9872.Contig2) was found to be significantly elevated on exposure of plants to heat stress. Initial increase in the expression of Hsf genes in C. nankingense triggered the up-regulation of downstream Hsp genes; these signaling processes

and transcription control activated stress-responsive mechanisms such as the multi-chaperone network, and subsequently resulted in thermotolerance and protection of protein homeostasis (Kotak et al., 2007; Wahid et al., 2007). The performance of the Hsp and Hsf family members in C. nankingense is consistent with the transcriptional profile of Hsps and Hsfs in maturing tomatoes under heat stress (Frank et al., 2009). 3.2. RNA metabolism is involved in the response to heat stress Regulation of RNA metabolism is also necessary to maintain protein dynamics. Serine/arginine-rich (SR) proteins belong to a conserved family of splicing regulators in eukaryotes; the pre-mRNA of SR genes in A. thaliana is extensively alternatively spliced under heat stress (Palusa et al., 2007). In heat-treated seedlings of Arabidopsis, the transcription of SR30 gene was affected by heat and decreased. However, a C. nankingense homolog of CnSR gene (Unigene26030) encoding the SR

Please cite this article as: Sun, J., et al., Identification of differentially expressed genes in Chrysanthemum nankingense (Asteraceae) under heat stress by RNA Seq, Gene (2014), http://dx.doi.org/10.1016/j.gene.2014.09.013

J. Sun et al. / Gene xxx (2014) xxx–xxx

5

Fig. 4. Differentially transcribed genes following high temperature stress treatment of C. nankingense. (a) The number of DT unigenes identified in each library contrast applying a threshold of the ratio change ≥2 and a q-value of b0.05. The red/green column represents genes up/down-regulated by the stress. (b) Scatter plot of DTGs (FDR ≤ 0.001 and |log2Ratio| ≥ 1) illustrating the full set of genes sampled. Red points: up-regulated genes; Green points: down-regulated genes; Blue points: not DTGs. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

splicing factor was up-regulated by 2.87-fold under heat stress, which indicates that different pathways may exist in different species. In A. thaliana, DEAD box RNA helicase is essential for mRNA export and acts as a negative regulator of high temperature tolerance (Palusa et al., 2007). In this study, the expression of five genes encoding C. nankingense DEAD-box RNA helicases, Unigene15341, CL14189.Contig3, Unigene17117, CL5642.Contig3 and CL11096.Contig3, was decreased in the heat stressed plants (Table S8). Intriguingly, inactivation of DEADbox RNA helicase resulted in enhanced resistance to oxidative stress instead of reducing the perceived stress in prokaryotes (Briolat and Reysset, 2002). Therefore, the down-regulation of five DEAD-box RNA helicase genes was expected to be involved in the export of HSP-related mRNAs and improve thermotolerance of C. nankingense. Some genes that were affected by heat stress were related to RNA modification or processing, including RNA polymerase, RNA methyltransferase, pre-mRNA splicing factor, RNA helicase, and U1 snRNP auxiliary factor in C. nankingense (Table S8). tRNAs are responsible for the transport of amino acids during protein biosynthesis. In this study, 22 of the DTGs identified were involved in tRNA synthesis, for example, arginyl-tRNA

synthetase (Fig. 6), leucyl-tRNA synthetase and glycyl-tRNA synthetase (Table S8). The RNA-Seq results suggest that the expression levels of these tRNAs were all increased, which indicate an improved level of protein biosynthesis during the short-term heat stress. 3.3. Transcription factors were affected by heat stress The expression of several classes of transcription factors was induced by high-temperature stress in C. nankingense (Table S3). In A. thaliana, the AP2/EREBP family is represented by 147 members and divided into four subfamilies, namely AP2, DREB, ERF and RAV (Dietz et al., 2010). In C. nankingense, six CnAP2 and CnERF members were down-regulated by the heat stress, while two CnDREB2 genes were obviously up-regulated. qRT-PCR validation showed that the expression of CnAP2 (CL10743.Contig1) was decreased in the first 2 h and then increased. This result was consistent with hs-induced AtAP2 (At4G28140) expression in A. thaliana (Zeller et al., 2009). In Arabidopsis, overexpression of DREB2A resulted in the induction of Hsps in non-stressed plants (Sakuma et al., 2006). DREB2A/2B is

Fig. 5. Gene Ontology (GO) classifications of DTGs. Most consensus sequences were grouped into three major functional categories, namely biological process, cellular component, and molecular function. Red column: up-regulated genes; green column: down-regulated genes. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Please cite this article as: Sun, J., et al., Identification of differentially expressed genes in Chrysanthemum nankingense (Asteraceae) under heat stress by RNA Seq, Gene (2014), http://dx.doi.org/10.1016/j.gene.2014.09.013

6

J. Sun et al. / Gene xxx (2014) xxx–xxx

Fig. 6. The relative expression of eight DTGs at different time points. The relative expression of each gene was calculated through 2−ΔΔCT method in the form of log2 ratio.

known to bind directly to the HsfA3 promoter in the early phase of high temperature stress (Schramm et al., 2008) (Fig. 7). With respect to MYB family, two (CL3691.Contig2 and Unigene52128) of ten CnMYB genes were up-regulated while the other eight were down-regulated in response to the heat stress. In plants, different MYB transcription factors play roles as activators, repressors, or both. Moreover, MYBs have been identified as immediate targets of other regulators (http://arabidopsis.med.ohio-state.edu/). In a recent study, OsMYB55 was found to be involved in tolerance to high temperature via modulation of amino acid metabolism (Thirunavukkarasu et al., 2013). This new function of MYB transcription factor in rice led us to assume that differentially regulated CnMYBs may contribute to the thermotolerance network in C. nankingense. A total of six WRKY genes, CnWRKY1, CnWRKY6, CnWRKY21, CnWRKY23, CnWRKY30 and CnWRKY65, were implicated in

the heat stress response, as their transcription was down-regulated in response to the high temperature. A set of stress-related genes up-regulated by PtWRKY23 overexpression was found to be related to heat shock proteins (Levée et al., 2009), but this is the first time the role of CnWRKY1, 6, 21, 30, and 65 in heat stress response has been reported. In a transcriptome analysis of hsfA1a/hsfA1b double-knockout A. thaliana mutants, WRKY7 was found to be HSF-dependent and down-regulated (Busch et al., 2005). These results illustrate that CnWRKYs may participate in the HSP/HSF signaling processes associated with transcriptional reprogramming when plants encounter heat stress. With global warming on the rise, plants will have to adapt to a more stressful environment. High-temperature stress has a major impact on plant growth and reproduction, so a deeper understanding of stress response can help establish genetic strategies needed to improve the

Fig. 7. A model of CnHsp and CnHsf network under heat stress. Red frame highlighted the DTGs detected in C. nankingense. The figure shows the CnHsp and CnHsf network involved in the heat stress response and protective factors leading to thermotolerance that are described in the text. When heat stress comes, the expression of Hsp70 and Hsp90 increased dramatically. They function not only as chaperon to help the folding of protein, but also active Hsfs to bind with HSE-containing genes in the nucleus. This process promotes the expression of a series of genes including CnHsp70, CnHsp90, HsfA1a, CnHsfA2 and HsfB1 (Baniwal et al., 2004). CnHsfA2 may play a major role to promote its downstream Hsps such as CnsHsp, CnHsp101, CnHsp70 and CnHsp90 may contribute to the thermotolerance of C. nankingense in vivo.

Please cite this article as: Sun, J., et al., Identification of differentially expressed genes in Chrysanthemum nankingense (Asteraceae) under heat stress by RNA Seq, Gene (2014), http://dx.doi.org/10.1016/j.gene.2014.09.013

J. Sun et al. / Gene xxx (2014) xxx–xxx

tolerance of temperate species such as C. nankingense. The present study has succeeded in identifying a range of DTGs in C. nankingense, many of them might play a key role in the process of hs response. Overexpression of several important DTGs such as CnHsp101, CnHsf A2 and CnDREB2 may bring a high level of plant temperature tolerance. Further research is being carried on to clarify the specific mechanism of upregulated network, some of which may provide useful leads in ongoing efforts to develop plants with specific tolerance to abiotic stress. 4. Materials and methods

7

GO (http://www.geneontology.org/) database was used for functional annotation. The RPKM algorithm was used to quantify transcript abundance (Mortazavi et al., 2008). The identifications of differentially transcribed genes (DTGs) were identified through an algorithm developed by Audic and Claverie (1997) Cluster (3.0) and Java TreeView (1.1.6r2) were used to analyze transcription patterns. A false discovery rate of ≤0.001 and a log2 transcript abundance ratio of ≥1 were used as the significance thresholds in order for the differences in transcript abundance to be considered as significant. Differentially transcribed genes were subjected to GO and KEGG ontology enrichment analysis, based on the Hypergeometric distribution model.

4.1. Plant materials, stress treatments and tissue sampling C. nankingense was obtained from the Chrysanthemum Germplasm Resource Preserving Centre, Nanjing Agricultural University, China. Disinfected shoot cuttings of the same size were grown for 50 days on solid MS containing no additives (Gamborg et al., 1976), and were cultured with a 16-h light period (100 μEm− 2 s− 1, 22 °C) and an 8-h dark period (18 °C). The relative humidity was set at 68–75%. When the plants grew into 10–15 leaf stage, half of them (when they were in light period) were transferred to another 40 °C chamber for 2 h under a 2-h light period (100 μEm−2 s−1) with 68–75% relative humidity. Control plants were kept in a 22 °C chamber for the same period. After exposure to the 40 °C temperature for 2 h, the second and third true leaves were collected. At the same time, leaf samples (the 2nd and 3rd leaves) of control plants were collected. Each sample comprised two leaves each from three plants, and was stored in liquid nitrogen until RNA extraction. Total RNA was extracted using a Plant RNA Kit (TaKaRa, Japan), following the manufacturer's protocol. The extract was enriched for mRNA using oligo (dT) magnetic beads (Dynabeads; Invitrogen, Carlsbad, CA). The resulting mRNA was fragmented into ~ 200 nt pieces and converted to cDNA using random hexamerprimers. After synthesis of the complementary strand, the resulting double-stranded cDNA was purified using magnetic beads and endrepaired, following which an A nucleotide was attached to the 3′ ends. Sequencing adaptors were ligated to the fragments for PCR amplification. An Agilent 2100 Bioanalyzer (Agilent, USA) and the StepOnePlus Real-Time PCR System (ABI, USA) were then used to quantify the sample. Sequencing of the double-stranded cDNAs was carried out using a HiSeq™ 2000 device (Illumina, USA). 4.2. RNA-seq analysis The sequencing data were analyzed (base calling, quality filtering, and per base confidence scoring) using the Illumina platform maintained at the Beijing Genomics Institute (Shenzhen, China; http:// www.genomics.cn/index.php), following the manufacturer's protocols. Because the genome of chrysanthemum has not been sequenced yet, a reference transcriptome database was established and the raw sequence data were deposited in the NCBI Sequence Read Archive database (Accession number: SRP029991, http://trace.ncbi.nlm.nih.gov/ Traces/sra_sub/sub.cgi?). This transcriptome de novo assembly was carried out with the short reads assembling program Trinity (Grabherr et al., 2011). Reads with a certain length of overlap were combined to form longer fragments, which are called contigs. The different contigs were connected to obtain sequences that could not be extended on either end; such sequences were defined as unigenes. By far, this transcriptome database included most known C. morifolium unigene sequences (Li et al., 2014). Through short reads alignment program SOAPaligner/SOAP2 (Li et al., 2009), clean reads of two C. nankingense libraries (control and heat stress plants) were directly mapped onto the C. morifolium reference gene database. Reads that had more than two-nucleotide mismatches or were mapped to multiple locations were excluded. The Blast algorithm was used to assess the clean reads distribution on the reference genes and a combination of Blast, Blast2GO (v2.2.5), KEGG (http://www. genome.jp/kegg/pathway.html), NR (http://blast.ncbi.nlm.nih.gov/) and

4.3. Independent validation of RNA-Seq through quantitative real-time PCR (qRT-PCR) To see the dynamic expression process of different genes during the heat stress, an independent experiment was carried out by qRT-PCR. Uniform rooted cuttings of C. nankingense were cultured in a 2:1 mixture of garden soil and vermiculite without fertilizer supplementation. They were grown under the same condition mentioned above. When they reached the 10- to 15-leaf stage, the plants were being subjected to the heat stress (40 °C) for six consecutive hours. The second and third true leaf samples were collected at different time points: 0, 0.5, 1, 2, 4 and 6 h. Samples collected at the 2 h time point were recognized as a validation of RNA Seq results. Eight DTGs were chosen including induced and repressed genes. Total leaf RNA (1 μg) extracted from each plant was reverse-transcribed in a 10-μl reaction based on the superscript first-strand synthesis system (Invitrogen, Foster City, CA, USA). qRT-PCR was conducted according to the guidelines of the Minimum Information for Publication of Quantitative Real-Time PCR Experiments (Bustin et al., 2009). Primers (Table S9) were designed using Vector NTI v11.0 software (http://www.lifetechnologies.com/cn/) and synthesized commercially (Genery, China). qRT-PCR was performed using the SYBR Premix Ex Taq™ Kit (TaKaRa) according to the manufacturer's protocol; three technical replicates were used in an Eppendorf Real Time PCR System (Mastercycler®ep Realplex, Germany). The PCR cycles were as follows: 1 cycle of 2 min at 95 °C, followed by 40 cycles at 95 °C for 15 s, 55 °C for 15 s and 68 °C for 20 s. Following amplification, all the products were subjected to melting curve analysis. A negative control without a cDNA template was used to evaluate the overall specificity. The C. nankingense homolog of EF1α (Oda et al., 2012) was used as an internal reference gene. Supplementary data to this article can be found online at http://dx. doi.org/10.1016/j.gene.2014.09.013. Acknowledgments This work was supported by funding from the Program for New Century Excellent Talents in University of Chinese Ministry of Education (NCET-12-0890, NCET-11-0669, NCET-10-0492), the National Natural Science Foundation of China (31171987, 31272196, 31372100), the Natural Science Fund of Jiangsu Province (BK201164, BK2012773), the Fundamental Research Funds for the Central Universities (KYZ201147), and the College Graduate Research and Innovation Plan of Jiangsu Province (CXZZ12_0286). References Audic, S., Claverie, J.M., 1997. The significance of digital gene expression profiles. Genome Res. 7 (10), 986–995. Baniwal, S.K., Bharti, K., Chan, K.Y., Fauth, M., Ganguli, A., Kotak, S., Mishra, S.K., Nover, L., Port, M., Scharf, K.-D., 2004. Heat stress response in plants: a complex game with chaperones and more than twenty heat stress transcription factors. J. Biosci. 29 (4), 471–487. Briolat, V., Reysset, G., 2002. Identification of the Clostridium perfringens genes involved in the adaptive response to oxidative stress. J. Bacteriol. 184 (9), 2333–2343. Busch, W., Wunderlich, M., Schöffl, F., 2005. Identification of novel heat shock factor-dependent genes and biochemical pathways in Arabidopsis thaliana. Plant J. 41 (1), 1–14.

Please cite this article as: Sun, J., et al., Identification of differentially expressed genes in Chrysanthemum nankingense (Asteraceae) under heat stress by RNA Seq, Gene (2014), http://dx.doi.org/10.1016/j.gene.2014.09.013

8

J. Sun et al. / Gene xxx (2014) xxx–xxx

Bustin, S.A., Benes, V., Garson, J.A., Hellemans, J., Huggett, J., Kubista, M., Mueller, R., Nolan, T., Pfaffl, M.W., Shipley, G.L., 2009. The MIQE guidelines: minimum information for publication of quantitative real-time PCR experiments. Clin. Chem. 55 (4), 611–622. Charng, Y.Y., Liu, H.C., Liu, N.Y., Chi, W.T., Wang, C.N., Chang, S.H., Wang, T.T., 2007. A heatinducible transcription factor, HsfA2, is required for extension of acquired thermotolerance in Arabidopsis. Plant Physiol. 143 (1), 251–262. Cheng, X., Chen, S., Chen, F., Fang, W., Deng, Y., She, L., 2010. Interspecific hybrids between Dendranthema morifolium (Ramat.) Kitamura and D. nankingense (Nakai) Tzvel. achieved using ovary rescue and their cold tolerance characteristics. Euphytica 172 (1), 101–108. Davidson, J.F., Schiestl, R.H., 2001. Mitochondrial respiratory electron carriers are involved in oxidative stress during heat stress in Saccharomyces cerevisiae. Mol. Cell. Biol. 21 (24), 8483–8489. Dietz, K.-J., Vogel, M.O., Viehhauser, A., 2010. AP2/EREBP transcription factors are part of gene regulatory networks and integrate metabolic, hormonal and environmental signals in stress acclimation and retrograde signalling. Protoplasma 245 (1–4), 3–14. Ellis, R.J., 2000. Chaperone substrates inside the cell. Trends Biochem. Sci. 25 (5), 210–212. Frank, G., Pressman, E., Ophir, R., Althan, L., Shaked, R., Freedman, M., Shen, S., Firon, N., 2009. Transcriptional profiling of maturing tomato (Solanum lycopersicum L.) microspores reveals the involvement of heat shock proteins, ROS scavengers, hormones, and sugars in the heat stress response. J. Exp. Bot. 60 (13), 3891–3908. Gamborg, O., Murashige, T., Thorpe, T., Vasil, I., 1976. Plant tissue culture media. In Vitro Cell. Dev. Biol. Plant 12 (7), 473–478. Grabherr, M.G., Haas, B.J., Yassour, M., Levin, J.Z., Thompson, D.A., Amit, I., Adiconis, X., Fan, L., Raychowdhury, R., Zeng, Q., 2011. Full-length transcriptome assembly from RNASeq data without a reference genome. Nat. Biotechnol. 29 (7), 644–652. Guan, J.-C., Jinn, T.-L., Yeh, C.-H., Feng, S.-P., Chen, Y.-M., Lin, C.-Y., 2004. Characterization of the genomic structures and selective expression profiles of nine class I small heat shock protein genes clustered on two chromosomes in rice (Oryza sativa L.). Plant Mol. Biol. 56 (5), 795–809. Hahn, A., Bublak, D., Schleiff, E., Scharf, K.-D., 2011. Crosstalk between Hsp90 and Hsp70 chaperones and heat stress transcription factors in tomato. Plant Cell Online 23 (2), 741–755. Hartl, F.U., Hayer-Hartl, M., 2002. Molecular chaperones in the cytosol: from nascent chain to folded protein. Science 295 (5561), 1852–1858. Hegedűs, Z., Zakrzewska, A., Ágoston, V.C., Ordas, A., Rácz, P., Mink, M., Spaink, H.P., Meijer, A.H., 2009. Deep sequencing of the zebrafish transcriptome response to mycobacterium infection. Mol. Immunol. 46 (15), 2918–2930. Jagadish, S., Muthurajan, R., Oane, R., Wheeler, T.R., Heuer, S., Bennett, J., Craufurd, P.Q., 2010. Physiological and proteomic approaches to address heat tolerance during anthesis in rice (Oryza sativa L.). J. Exp. Bot. 61 (1), 143–156. Katiyar-Agarwal, S., Agarwal, M., Grover, A., 2003. Heat-tolerant basmati rice engineered by over-expression of hsp101. Plant Mol. Biol. 51 (5), 677–686. Kotak, S., Larkindale, J., Lee, U., von Koskull-Döring, P., Vierling, E., Scharf, K.-D., 2007. Complexity of the heat stress response in plants. Curr. Opin. Plant Biol. 10 (3), 310–316. Koussevitzky, S., Suzuki, N., Huntington, S., Armijo, L., Sha, W., Cortes, D., Shulaev, V., Mittler, R., 2008. Ascorbate peroxidase 1 plays a key role in the response of Arabidopsis thaliana to stress combination. J. Biol. Chem. 283 (49), 34197–34203. Larkindale, J., Knight, M.R., 2002. Protection against heat stress-induced oxidative damage in Arabidopsis involves calcium, abscisic acid, ethylene, and salicylic acid. Plant Physiol. 128 (2), 682–695. Lee, S., Lee, D.W., Lee, Y., Mayer, U., Stierhof, Y.-D., Lee, S., Jürgens, G., Hwang, I., 2009. Heat shock protein cognate 70-4 and an E3 ubiquitin ligase, CHIP, mediate plastiddestined precursor degradation through the ubiquitin–26S proteasome system in Arabidopsis. Plant Cell Online 21 (12), 3984–4001. Levée, V., Major, I., Levasseur, C., Tremblay, L., MacKay, J., Séguin, A., 2009. Expression profiling and functional analysis of Populus WRKY23 reveals a regulatory role in defense. New Phytol. 184 (1), 48–70. Li, R., Yu, C., Li, Y., Lam, T.-W., Yiu, S.-M., Kristiansen, K., Wang, J., 2009. SOAP2: an improved ultrafast tool for short read alignment. Bioinformatics 25 (15), 1966–1967. Li, H., Chen, S., Song, A., Wang, H., Fang, W., Guan, Z., Jiang, J., Chen, F., 2014. RNA-Seq derived identification of differential transcription in the chrysanthemum leaf following inoculation with Alternaria tenuissima. BMC Genomics 15 (1), 9. Liu, S., Chen, S., Chen, Y., Guan, Z., Yin, D., Chen, F., 2011. In vitro induced tetraploid of Dendranthema nankingense (Nakai) Tzvel. shows an improved level of abiotic stress tolerance. Sci. Hortic. 127 (3), 411–419.

Livak, K.J., Schmittgen, T.D., 2001. Analysis of relative gene expression data using realtime quantitative PCR and the 2−ΔΔCT method. Methods 25 (4), 402–408. Mehlen, P., Kretz-Remy, C., Preville, X., Arrigo, A.-P., 1996. Human hsp27, Drosophila hsp27 and human alphaB-crystallin expression-mediated increase in glutathione is essential for the protective activity of these proteins against TNFalpha-induced cell death. EMBO J. 15 (11), 2695. Mishra, S.K., Tripp, J., Winkelhaus, S., Tschiersch, B., Theres, K., Nover, L., Scharf, K.-D., 2002. In the complex family of heat stress transcription factors, HsfA1 has a unique role as master regulator of thermotolerance in tomato. Genes Dev. 16 (12), 1555–1567. Morimoto, R.I., 1998. Regulation of the heat shock transcriptional response: cross talk between a family of heat shock factors, molecular chaperones, and negative regulators. Genes Dev. 12 (24), 3788–3796. Morimoto, R.I., Santoro, M.G., 1998. Stress-inducible responses and heat shock proteins: new pharmacologic. Nat. Biotechnol. 16. Mortazavi, A., Williams, B.A., McCue, K., Schaeffer, L., Wold, B., 2008. Mapping and quantifying mammalian transcriptomes by RNA-Seq. Nat. Methods 5 (7), 621–628. Nishizawa, A., Yabuta, Y., Yoshida, E., Maruta, T., Yoshimura, K., Shigeoka, S., 2006. Arabidopsis heat shock transcription factor A2 as a key regulator in response to several types of environmental stress. Plant J. 48 (4), 535–547. Nover, L., Scharf, K.-D., 1997. Heat stress proteins and transcription factors. Cell. Mol. Life Sci. 53 (1), 80–103. Nover, L., Bharti, K., Döring, P., Mishra, S.K., Ganguli, A., Scharf, K.-D., 2001. Arabidopsis and the heat stress transcription factor world: how many heat stress transcription factors do we need? Cell Stress Chaperones 6 (3), 177. Oda, A., Narumi, T., Li, T., Kando, T., Higuchi, Y., Sumitomo, K., Fukai, S., Hisamatsu, T., 2012. CsFTL3, a chrysanthemum FLOWERING LOCUS T-like gene, is a key regulator of photoperiodic flowering in chrysanthemums. J. Exp. Bot. 63 (3), 1461–1477. Palusa, S.G., Ali, G.S., Reddy, A.S., 2007. Alternative splicing of pre-mRNAs of Arabidopsis serine/arginine-rich proteins: regulation by hormones and stresses. Plant J. 49 (6), 1091–1107. Parmesan, C., Yohe, G., 2003. A globally coherent fingerprint of climate change impacts across natural systems. Nature 421 (6918), 37–42. Sakuma, Y., Maruyama, K., Qin, F., Osakabe, Y., Shinozaki, K., Yamaguchi-Shinozaki, K., 2006. Dual function of an Arabidopsis transcription factor DREB2A in water-stressresponsive and heat-stress-responsive gene expression. Proc. Natl. Acad. Sci. 103 (49), 18822–18827. Schramm, F., Larkindale, J., Kiehlmann, E., Ganguli, A., Englich, G., Vierling, E., Von KoskullDöring, P., 2008. A cascade of transcription factor DREB2A and heat stress transcription factor HsfA3 regulates the heat stress response of Arabidopsis. Plant J. 53 (2), 264–274. Teale, W.D., Paponov, I.A., Palme, K., 2006. Auxin in action: signalling, transport and the control of plant growth and development. Nat. Rev. Mol. Cell Biol. 7 (11), 847–859. Thirunavukkarasu, N., Hossain, F., Mohan, S., Shiriga, K., Mittal, S., Sharma, R., Singh, R.K., Gupta, H.S., 2013. Genome-wide expression of transcriptomes and their coexpression pattern in subtropical maize (Zea mays L.) under waterlogging stress. PLoS One 8 (8), e70433. Tissieres, A., Mitchell, H.K., Tracy, U.M., 1974. Protein synthesis in salivary glands of Drosophila melanogaster: relation to chromosome puffs. J. Mol. Biol. 84 (3), 389–398. Wahid, A., Gelani, S., Ashraf, M., Foolad, M., 2007. Heat tolerance in plants: an overview. Environ. Exp. Bot. 61 (3), 199–223. Yang, N., Ren, A., Hu, W., Qian, S., Duan, J., Tian, L., 2005. Chemical constituents of Dendranthema nankingense. J. China Pharm. Univ. 36 (5), 402. Yang, W., Glover, B.J., Rao, G.Y., Yang, J., 2006. Molecular evidence for multiple polyploidization and lineage recombination in the Chrysanthemum indicum polyploid complex (Asteraceae). New Phytol. 171 (4), 875–886. Yokotani, N., Ichikawa, T., Kondou, Y., Matsui, M., Hirochika, H., Iwabuchi, M., Oda, K., 2008. Expression of rice heat stress transcription factor OsHsfA2e enhances tolerance to environmental stresses in transgenic Arabidopsis. Planta 227 (5), 957–967. Zeller, G., Henz, S.R., Widmer, C.K., Sachsenberg, T., Rätsch, G., Weigel, D., Laubinger, S., 2009. Stress-induced changes in the Arabidopsis thaliana transcriptome analyzed using whole-genome tiling arrays. Plant J. 58 (6), 1068–1082. Zhang, X.-P., Glaser, E., 2002. Interaction of plant mitochondrial and chloroplast signal peptides with the Hsp70 molecular chaperone. Trends Plant Sci. 7 (1), 14–21. Zonia, L., Munnik, T., 2006. Cracking the green paradigm: functional coding of phosphoinositide signals in plant stress responses. Biology of Inositols and Phosphoinositides. Springer, pp. 207–237.

Please cite this article as: Sun, J., et al., Identification of differentially expressed genes in Chrysanthemum nankingense (Asteraceae) under heat stress by RNA Seq, Gene (2014), http://dx.doi.org/10.1016/j.gene.2014.09.013