Whole-Transcriptome Analysis of Differentially Expressed Genes in the Vegetative Buds, Floral Buds and Buds of Chrysanthemum morifolium

Background Chrysanthemum morifolium is an important floral crop that is cultivated worldwide. However, due to a lack of genomic resources, very little information is available concerning the molecular mechanisms of flower development in chrysanthemum. Results The transcriptomes of chrysanthemum vegetative buds, floral buds and buds were sequenced using Illumina paired-end sequencing technology. A total of 15.4 Gb of reads were assembled into 91,367 unigenes with an average length of 739 bp. A total of 43,137 unigenes showed similarity to known proteins in the Swissprot or NCBI non-redundant protein databases. Additionally, 25,424, 24,321 and 13,704 unigenes were assigned to 56 gene ontology (GO) categories, 25 EuKaryotic Orthologous Groups (KOG) categories, and 285 Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways, respectively. A total of 1,876 differentially expressed genes (DEGs) (1,516 up-regulated, 360 down-regulated) were identified between vegetative buds and floral buds, and 3,300 DEGs (1,277 up-regulated, 1,706 down-regulated) were identified between floral buds and buds. Many genes encoding important transcription factors (e.g., AP2, MYB, MYC, WRKY, NAC and CRT) as well as proteins involved in carbohydrate metabolism, protein kinase activity, plant hormone signal transduction, and the defense responses, among others, were considerably up-regulated in floral buds. Genes involved in the photoperiod pathway and flower organ determination were also identified. These genes represent important candidate genes for molecular cloning and functional analysis to study flowering regulation in chrysanthemum. Conclusion This comparative transcriptome analysis revealed significant differences in gene expression and signaling pathway components between the vegetative buds, floral buds and buds of Chrysanthemum morifolium. A wide range of genes was implicated in regulating the phase transition from vegetative to reproductive growth. These results should aid researchers in the study of flower-time regulation, breeding and molecular biology in chrysanthemum.

Chrysanthemum morifolium. A wide range of genes was implicated in regulating the phase transition from vegetative to reproductive growth. These results should aid researchers in the study of flower-time regulation, breeding and molecular biology in chrysanthemum.

Introduction
Asteraceae belongs to the euasterids II clade of the core eudicots, one of the largest families of flowering plants. The head-like inflorescence (capitulum) of Asteraceae resembles a single large flower, and this trait is considered to be the key innovation behind the evolutionary success of the Asteraceae [1]. Thus, the trait has been the subject of many phylogenetic and biological studies. In recent years, molecular genetic research on flower development has been restricted to Gebera hybrida and Helianthus annuus in Asteraceae. Chrysanthemum morifolium, another species of Asterceae, is an important crop that is widely cultivated, and indeed, complete year-round production of chrysanthemums is necessary to meet the demands of a growing market [2]. However, the flowering mechanisms and molecular pathways of chrysanthemums involved in flower development have not been well characterized to date due to the lack of genomic information [2].
In temperate regions, the duration of the daily light period, also called the photoperiod, is one of the most important factors controlling flowering time [3]. In Arabidopsis, many genes involved in daylength response have been identified through molecular-genetic approaches, including the genes encoding regulatory proteins for the regulation of flowering, light signal transduction pathways and circadian clock function, such as CRYPTOCHROME2/FHA (CRY2), GIGANTEA (GI), TIMING OF CAB 1 (TOC1) and CONSTANS (CO) [3]. Some important downstream regulatory genes have been identified in the photoperiod pathway. For example, the crucial regulator LEAFY (LFY) functions in floral specification. LFY exhibits weak expression in vegetative tissues, and its expression is significantly up-regulated after receiving photoperiodic signals regulated by the FT pathway and gibberellins [4,5]. Many MADS-Box transcription factors members, along with APETALA2 (AP2), also play important roles in floral initiation, such as SUPPRESSOR OF CONSTANS1 (SOC1), SHORT VEGETA-TIVE PHASE (SVP), LEAFY (LFY), APETALA1 (AP1), GIGANTEA (GI), AGAMOUS (AG) and PISTILLATA (PI). After the initial expression of LFY, AP1 is expressed throughout the floral meristem, and together with AP2, specifies sepal identity [4,6]. Various feedback mechanisms are used to regulate the actions of these genes in floral meristem specification and in the transition to flowering that occurs by promoting floral initiation and repressing vegetative growth [4]. The classic ABC model includes A-class genes (AP1, AP2), which determine sepal identity, A-class genes together with B-class genes (AP3, PI), which determine petals, B-class genes together with a C-class gene (AG), which specify stamens, and the C-class gene, which specifies carpel identity [7].
Chrysanthemum is a typical short-day plant. Under short-day treatment, chrysanthemum undergoes a major physiological change and morphological programs as it transits from a vegetative growth period to a reproductive period and finally blooms. This process is regulated by a combination of endogenous and environmental signals. By morphological observation, three developmental periods can be identified during the developmental process from vegetative growth to reproductive growth, including vegetative bud growth, flower bud differentiation, and bud development. In chrysanthemum, only a handful of important transcription factors involved in the flower developmental pathway have been isolated and analyzed, including the homologous genes APETALA1, SEPALLATA3, and FRUITFULL [8]. Therefore, the identification of novel genes and molecular regulatory pathways with important functions will be necessary to improve year-round production techniques in chrysanthemum.
Transcription expression profiling, data assembly, and analysis can provide broad and deep insights into gene regulatory networks and biological pathways because they can reveal the genes downstream of key transcription factors in related pathways [9]. The use of RNA-Seq in the studies of chrysanthemum and related species has been demonstrated [10][11][12]. However, to date, no transcriptome information concerning the different flower developmental periods of chrysanthemum has been reported. In this study, transcriptional sequencing and analysis of vegetative buds (the apical buds during vegetative growth), floral buds (the apical buds during flower differentiation) and buds of chrysanthemum were performed using Illumina assembly technology and RNA-Seq quantification analysis. Based on transcriptional sequencing and analysis, we first identified differentially expressed genes (DEGs) between the vegetative buds, floral buds and buds to uncover novel genes involved in regulating the phase transition from vegetative to reproductive growth. Next, we identified important regulatory genes involved in the photoperiod pathway and the control of flower organ identification to create a list of candidate genes for studying flowering-time regulation in chrysanthemum. These results will be helpful for elucidating the molecular mechanisms of flower development and will contribute to the development of techniques for studying flowering-time regulation, breeding and molecular biology in chrysanthemum.

Illumina sequencing and assembly
The cDNA libraries of vegetative buds, floral buds and buds were sequenced using the Illumina HiSeq 2000 system. After stringent quality checks and data cleaning, 58,530,798, 41,424,050 and 54,035,520 clean reads were left for vegetative buds, floral buds and buds, respectively. The average proportion of clean reads for these libraries was approximately 93%.
A total of 123,354 contigs were assembled based on the high-quality reads, with a total size of 67,027,473 bp, an N50 of 641 bp and an average contig length of 543 bp. Then, scaffolds were constructed between the contigs via the paired-end relationships between the reads. A total of 95,280 scaffolds were obtained, with an N50 of 818 bp and an average length of 721 bp. We filled the intra-scaffold gaps and constructed a non-redundant unigene set from all three of the assembled datasets using CAP3 software. Finally, a total of 91,367 high-quality unigene sequences with lengths greater than 400 bp were obtained, with an average unigene length of 739 bp (Table 1).

Gene annotation and functional classification
Of the 91,367 unigenes, 43,137 (47.21% of the total) were aligned to the Nr protein (date 2014.03) and Swissprot protein databases (date 2014.03) using an E-value threshold of < 1e-5, which meant that 48,230 unigenes (52.79% of the total) had no Swiss-Prot annotations because of missing chrysanthemum genome and expressed sequence tag (EST) information ( Table 1).
The international standardized gene functional classification system Gene Ontology (GO) provided three ontologies-molecular functions, cellular components and biological processesthat were useful for gene annotation and analysis. Based on the Nr annotation and Swissprot protein databases, 25,424 unigenes were classified into 56 functional GO categories using Blast2GO software [13]. A total of 19,161 unigenes were classified into 25 biochemical process categories, 10,063 unigenes were classified into 17 cellular component categories, and 22,539 unigenes were classified into 14 molecular function categories (Fig 1). For each of the three main GO classifications, the 'metabolic process' (in 'biological process'), 'cell part' (in 'cell component'), and 'binding' (in 'molecular function') terms were dominant in level 2. And 17,815 (92.98%), 7244 (71.98%), and 15,516 (68.84%) unigenes were assigned to 'metabolic process' (GO:0008152), 'cell part' (GO:0044464) and 'binding' (GO:0005488), respectively. It was indicated that complex metabolic processes regulated by a wide range of genes that also interact with one another in cell parts were activated during the development of flower in chrysanthemum. In addition, we also noted high percentages of genes from the 'cellular process' (GO:0009987), 'cell' (GO:0005623), and 'catalytic activity' (GO:0003824) categories, which included 13,907 (72.58%), 7244 (71.98%) and 14,876 (65.99%) unigenes, respectively. By contrast, there were relatively few genes from the 'cell killing' (GO:0001906), 'extracellular matrix part' (GO:0044420), and 'protein tag' (GO:0031386) categories in each of the three main GO classifications (Fig 1).
Within the 'biological process' category, we identified 13 terms related to flower development in level 6, including In total, 113 unigenes were assigned to these 13 terms, of which 86 were annotated as uncharacterized or predicted proteins and the remainder showing homology to F-box family genes and GIGAN-TEA etc. (S1 Table).
As a eukaryote-specific version of the Clusters of Orthologous Groups (COG) tool, euKaryotic Orthologous Groups (KOG) is used to identify orthologous and paralogous proteins, providing a way to identify Joint Genome Institute (JGI)-predicted genes by KOG classification or ID. The annotated sequences were further searched for genes involved in KOG classifications for evaluating the completeness of our transcriptome library and the effectiveness of our annotation process. Of 43,137 Nr hits, 24,321 sequences were assigned to KOG classifications. Among the 25 KOG categories, the cluster for 'Signal transduction mechanisms' (3487, 14.34%) represented the largest group, followed by 'General function prediction only' (2798, 11.50%) and 'Posttranslational modification, protein turnover, chaperones' (2453, 10.09%). The 'Nuclear structure' (108, 0.44%), 'Extracellular structures' (104, 0.43%), and 'Cell motility' (8, 0.032%) categories represented the smallest groups (Fig 2).
To further analyze the chrysanthemum transcriptome, all unigenes were compared with the Kyoto Encyclopedia of Genes and Genomes (KEGG) database using BLASTx with an E-value threshold of < 1e-5. Of the 43,137 unigenes, 13,704 had significant matches in the database and were assigned to 285 KEGG pathways ( Table 2). The most represented pathways were 'Metabolic pathways' (3426 members), 'Biosynthesis of secondary metabolites' (1881   Transcriptome Analysis of Chrysanthemum morifolium

Comparison of transcriptomes between vegetative buds, floral buds and buds
The common set of unigenes between vegetative buds, floral buds and buds. The number of unigenes with an RPKM value of > 0.3 that were shared by vegetative buds and floral buds, floral buds and buds, and vegetative buds and buds were 80,769, 81,794, and 80,984, respectively. Vegetative buds, floral buds and buds shared 78,561 unigenes in common (Fig 3). By contrast, 496, 1294 and 3128 unigenes showed specific expression in vegetative buds, floral buds and buds, respectively. During development from vegetative growth to flower differentiation and bud formation, an increasing number of genes are expressed in chrysanthemum.
DEGs between vegetative buds, floral buds and buds. The resulting reads were mapped to our reference transcriptome; statistical significance values were taken to be reliable with an RPKM value 2 in at least one of the three samples. Transcriptomes were compared, and differentially expressed genes (DEGs) were analyzed in three samples: vegetative buds, floral buds and buds. To determine which of the unigenes were differentially expressed, we filtered with an FDR 0.05 and a |log2 (ratio)| 2; 1,876 DEGs were significantly changed in expression between vegetative buds and floral buds, of which 1,516 DEGs were up-regulated and 360 DEGs were down-regulated in the floral buds ( Fig 4A).
A total of 2,356 DEGs (983 down-regulated and 1,373 up-regulated in buds) were identified between buds and vegetative buds (Fig 4A). In addition, we also identified 174 DEGs specifically expressed in buds relative to floral buds, which contained many genes involved in pistil and stamen development (S6 Table).  To identify unigenes with significant changes in expression, the correlations in gene expression between the three samples were studied using an algorithm developed from the correlation scatter plot. As shown in Fig 4B, the difference in gene expression was smallest between vegetative buds and floral buds, yielding a strong correlation in gene expression (R 2 = 0.907). By contrast, the difference in gene expression was the largest between floral buds and buds, leading to a poor correlation in gene expression (R 2 = 0.831, Fig 4C).
Among the samples from three different developmental periods in chrysanthemum, the number of up-regulated genes in floral buds relative to vegetative buds (1516) was greater than the number in buds relative to floral buds (1277). Fewer genes (360) were down-regulated in floral buds relative to vegetative buds than the number down-regulated (1706) in buds relative to floral buds. Indeed, there were only 90 common DEGs among vegetative buds, floral buds and buds (Fig 5). In total, 1301 DEGs in vegetative buds and floral buds were also differentially expressed in floral buds and buds. By contrast, 397 DEGs in vegetative buds and floral buds were also differentially expressed in buds and vegetative buds (Fig 5). We also found that 1367 DEGs in floral buds and buds were differentially expressed in buds and vegetative buds (Fig 5).
Up-regulated unigenes during flower differentiation. By comparing the transcriptomes of vegetative buds, floral buds and buds, we found that many unigenes were dramatically upregulated in floral buds and down-regulated in both vegetative buds and buds, suggesting that these genes play an important role during flower differentiation.
The expression levels of many important transcription factor genes, including the homologs of AP2, MYB, MYC, WRKY, NAC and CRT, were dramatically up-regulated in floral buds (Table 6 and Fig 6). Interestingly, we identified one CRT-binding gene (CRTCL) that was highly expressed in floral buds, weakly expressed in vegetative buds, and not expressed at all in buds (Table 6). One homolog of AP2 (AP2CL3) was assigned to "organ morphogenesis" (GO:0009887) under the 'biological process' category. Furthermore, many other homologs of these transcription factors were assigned to the 'nucleus' (GO:0005634), 'DNA binding' (GO:0003677) and 'regulation of transcription, DNA-dependent' (GO:0006355) categories under 'cellular component,' 'molecular function' and 'biological process,' respectively. The genes that were up-regulated during flower differentiation also included many genes that are responsive to phytohormones or are involved in phytohormone synthesis, as well as sugar synthesis/transport genes, kinase-like genes, ubiquitin-like genes, and resistance genes; numerous uncharacterized protein coding genes were also identified. These genes were enriched for the 'carbohydrate metabolic process' (GO:0005975), 'signal transduction' (GO:0007165), 'protein kinase activity' (GO:0004672) and 'defense response' (GO:0006952) categories and are listed in S7-S10 Tables.

Identification of genes involved in the photoperiod pathway in chrysanthemum
Chrysanthemum is a typical short-day plant; it can flower in response to a single short day. We searched for homologous genes of the important regulators involved in the photoperiod pathway in chrysanthemum. In Arabidopsis, many genes required for the daylength response have been identified using molecular-genetic approaches, among which some members encode regulatory proteins specifically playing roles in the regulation of flowering; others encode components in light signal transduction pathways or pathways related to circadian signaling, including CRYPTOCHROME (CRY), PHYTOCHROME (PHY), LATE  GIGANTEA (GI) and others [3]. In this study, many unigenes from the chrysanthemum transcriptome were identified as homologous genes of the photoreceptor and circadian clock components in the photoperiod pathway (Fig 7).  (Fig 7). In the photoperiod pathway, many MADS box genes also play key roles in promoting floral meristem identity, including SUPPRESSOR OF CONSTANS1 (SOC1), SHORT VEGETATIVE PHASE (SVP), LEAFY (LFY), APETALA1 (AP1), GIGANTEA (GI), AGAMOUS (AG) and PISTILLATA (PI). APETALA2 (AP2) also play key roles in promoting floral meristem identity. And homologs of these MADS box genes and AP2 were also identified (Fig 7). LEAFY (LFY) is a vital regulator of the specification of floral meristem identity. This gene is initially expressed very early throughout the presumptive floral meristem, and its expression leads to a cascade of transcriptional activities that control floral meristem development and meristem identity [4]. We found one homolog of LFY (LFYCL) that was significantly up-regulated in floral buds. SOC1 is an upstream regulatory gene of LFY. Two homologs of SOC1, SOC1CL1 and SOC1CL2, were identified in chrysanthemum, and their expression levels were significantly up-regulated in the floral bud stage. AP1 is 'A-class' like gene; LFY directly activates AP1 transcription [4,14]. In Arabidopsis, AP1 and LFY are expressed throughout the floral meristem and are thought to regulate floral initiation via cell-to-cell signaling [15]. AP1CL1 and AP1CL2, which are  homologous to AP1, were initially expressed in floral buds and exhibited higher expression levels in buds in chrysanthemum. The third homolog of AP1, AP1CL3, exhibited higher expression in vegetative buds, suggesting that AP1CL3 is probably not involved in floral meristem identity. AP2 is the second 'A-class' gene that is not a MADS-box transcription factor. Many homologs of AP2 were identified in chrysanthemum (Fig 7), and the majority were significantly up-regulated in floral buds. The B-class genes of most core eudicot species include three different lineages: PI, euAP3 and TM6. TM6-like genes appear to have been lost in Arabidopsis and Antirrhinum [16]. In chrysanthemum, we identified homologs for all three genes PI, euAP3 and TM6 (PICL, AP3CL, TM6CL). In addition, we identified one homolog of the Cclass gene AG (AGCL) in chrysanthemum. The homologs of AP3, PI, TM6 and AG, were initially expressed in the floral buds, although they showed the highest expression in buds. Annotation information for unigenes involved in the photoperiod pathway in chrysanthemum is listed in S11 Table. Verification of gene expression profiles by qRT-PCR To further confirm the expression data of the genes in the Illumina sequencing analyses, we randomly selected 20 unigenes for qRT-PCR evaluation using the vegetative bud, floral bud and bud samples that were originally used for RNA-Seq (Fig 8). Among the 20 selected unigenes, 10 were annotated as known proteins (MADS-boxCL1, MADS-boxCL2, MADS-boxCL3, GDEF1CL, AP2/ ERFCL3, CRTCL, MYBCL1, DoubleWRKYCL, GA2O2CL, ETHR3CL) and 10 were annotated as uncharacterized or predicted proteins (UnknownC1, UnknownC2, UnknownC3, UnknownC4, UnknownC5, UnknownC6, UnknownC7, UnknownC8, UnknownC9, UnknownC10). The expression of 4 genes (GDEF1CL, MADS-boxCL1, MADS-boxCL2 and MADS-boxCL3) was only weakly detected in vegetative buds and floral buds by qRT-PCR (Fig 8A). Consistent with this, these genes were also assigned minimal expression levels by RNA-Seq analysis (S12 Table). The expression patterns of the 20 selected genes were nearly identical, and the gene expression levels detected by Illumina sequencing analysis were consistent with the qRT-PCR results.

Discussion
For the chrysanthemum transcriptome, approximately 16.5 Gb of data were generated and assembled into 91,367 unigenes with an average length of 739 bp. However, only 43,137 unigenes (47.21% of the total) were annotated in the Swiss-Prot database, and the remaining 48,230 unigenes (52.79% of the total) had no Swiss-Prot annotations due to the lack of chrysanthemum genomic and EST information. The assembly unigenes in this study were assigned to an extensive range of GO categories, KEGG pathways and KOG classifications, which indicated that various transcripts are involved in chrysanthemum flower development. Comparative analysis of the transcriptomes revealed a wealth of valuable information. As shown in Fig 3, vegetative buds, floral buds and buds shared 78,561 common unigenes, whereas only 496, 1294 and 3128 unigenes were specifically expressed in vegetative buds, floral buds and buds, respectively. Therefore, the majority of unigenes were expressed during all three developmental stages. However, during flower bud differentiation and bud development, many unigenes were up-regulated and a small number of unigenes were turned off. Many unigenes were significantly up-regulated during flower differentiation, including many important transcription factors shown in Table 6. Previous studies have shown that the AP2/ERF transcription factors are expressed in response to developmental and environmental stimuli; some also played roles downstream of the ethylene, biotic, and abiotic stress signaling pathways in model plants [17]. WRKY proteins play regulatory roles in plant defense responses [18]. MYC/MYB transcription factors are involved in the ABA-dependent pathway to up-regulate abiotic stressresponsive genes [19,20]. NAC proteins are plant-specific transcription factors that function in diverse developmental processes and stress responses, though the biological functions of most NAC family genes remain unknown [21]. In the chrysanthemum transcriptome, AP2, WRKY, MYB, MYC, NAC genes are significantly up-regulated in floral buds relative to vegetative buds and buds, indicating that these genes play important roles in chrysanthemum flower development and should be addressed in future studies (Table 6). Interestingly, two CRT genes were identified as highly expressed in floral buds. Both CRT genes (CRT/DRECL, CRTCL) showed minimal expression in vegetative buds, and one CRT gene (CRTCL) was not detected in buds. It is known that CRT genes are important regulators of plant reproductive events [22]. However, our findings suggest that CRT genes likely play a vital role in regulating chrysanthemum flower bud differentiation.
Among the important transcription factors mentioned above, some are involved in the same biological processes as genes involved in phytohormone signaling, the photoperiod pathway and organ determination. Genes involved in signal transduction as well as sugar, kinase, ubiquitin and defense pathways were also up-regulated in floral buds relative to vegetative buds and buds. These up-regulated genes were enriched for the GO categories 'carbohydrate metabolic process,' 'protein kinase activity,' 'plant hormone signal transduction,' and 'defense response,' among others. These findings indicate that genes involved in carbohydrate metabolic processes, protein kinase activities, plant hormone signal transduction, defense responses and other related processes play a role in floral differentiation in chrysanthemum. This hypothesis is also consistent with previous studies showing that 1) various phytohormones play an important role in the regulation of flowering in both long-and short-day plants; 2) the soluble sugar and starch contents of leaves and buds are temporarily increased in Sinapis during the induction of flowering; and 3) a calcium-dependent protein kinase is involved in photoperiodic flower induction in Pharbitis nil. [23,24,25].
In addition, there were many unigenes coding uncharacterized proteins that are also worthy of attention. Some of these were assigned to flower development-related GO terms, and many others were significantly up-regulated in floral buds, suggesting that they are involved in regulating chrysanthemum flower development. Further studies should be conducted to reveal the exact roles of these uncharacterized genes in the regulation of flower development.
The chrysanthemum floral meristem initiates in response to a single short day; the photoperiod is considered to be the most important factor controlling chrysanthemum flowering time. We identified many genes (CRY, PHY, LHY, FKF1, TOC1, ZTL, GI, CO, FT) encoding components of light signal transduction pathways or related to circadian rhythms. LFY expression leads to a series of transcriptional activities that control floral meristem formation and meristem identity [4]. The homologous gene (LFYCL) in chrysanthemum was identified. In addition, the homologs of AP2 and many MADS box genes involved in floral meristem identity were identified in chrysanthemum. According to the ABC model, A function genes (AP1, AP2) specify the identity of the sepal, A function genes together with B function genes (AP3, PI, and TM6) specify the petals, B function genes together with C function genes (AG) specify the stamen, and C function genes specify the carpel. However, the capitulum of chrysanthemum is composed of hundreds of individual florets which are specialized in structure and function, indicating that the regulatory mechanism of A, B and C function genes should be further studied in this species. The homologous genes of A, B and C function genes were identified in this study. These genes are important candidate genes for molecular cloning and functional analysis of flowering regulation genes in chrysanthemum.

Conclusion
In this study, comparative transcriptome analysis revealed significant differences in gene expression and signaling pathways among the vegetative buds, floral buds and buds in Chrysanthemum morifolium. We found that genes involved in specific biological processes were activated and considerably up-regulated during the process of flower bud differentiation, including genes coding for important transcription factors, as well as genes involved in carbohydrate metabolic processes, protein kinase activities, plant hormone signal transduction and defense responses. These findings indicate that the phase transition from vegetative to reproductive growth is a complex biological process regulated by a wide range of factors, during which time many physiological processes, such as sugar metabolism and phytohormone regulation, become activated. The identification of homologous genes involved in the photoperiod pathway will likely contribute to studies of flowering-time regulation in chrysanthemum.
The transition from vegetative growth to the reproductive period in chrysanthemum is regulated by a large number of genes involved in complex molecular pathways. These results represent a first step towards illuminating the molecular mechanisms of flower development in chrysanthemum, and they will provide abundant genomic resources and new candidate genes for studies of flower-time regulation, breeding and molecular biology in this plant species.

Plant materials and RNA extraction
The tissues (vegetative buds, floral buds, buds) used in the study were obtained from a groundcover chrysanthemum variety (Chrysanthemum morifolium 'Fenditan', a hybrid of chrysanthemum varieties) cultivated in a greenhouse under a 16-hr light/8-hr dark cycle for 180 days and then an 8-hr light/16-hr dark cycle at 23°C in Beijing Forestry University (116.3°E, 40.0°N). Under long-day treatment, 70 vegetative buds were collected between 9:00-12:00 a.m. When the short-day treatment began, 70 apical buds were collected between 9:00-12:00 a.m. every week until visible floral buds were formed. After floral buds were formed, 30-50 buds were collected between 9:00-12:00 a.m. every week until the ray florets started to show color. All plant tissues collected were placed immediately in liquid N 2 and stored at −70°C until RNA was extracted. The individual plants used for collecting materials were all derived by tissue culture from the C. morifolium genotype 'Fenditan.' Total RNA was extracted using an RNAisomate RNA Easyspin Isolation System, and a NanoDrop ND2000 instrument was used to assess the RNA quantity and quality.

Illumina sequencing, de novo assembly and functional annotation
To reveal differences in the chrysanthemum transcriptome before and after flower differentiation, we collected apical buds during the adult vegetative phase before short-day treatment. However, flower differentiation and bud development are continuous developmental process. Therefore, to obtain more complete transcriptome information, we pooled equivalent quantities of total RNA isolated from different developmental stages of flower differentiation and bud development, and the pooled total RNA samples were used for sequencing. In total, 10 μg total RNA was collected for each sample for sequencing. We conducted Illumina sequencing using an Illumina HiSeq 2000 system according to the manufacturer's protocols with paired-end 2×100-nt multiplexing according to the manufacturer's instructions.Briefly, poly (A) mRNA was purified using magnetic beads with oligo (dT) and was then fragmented into small pieces. After the first-strand cDNAs were synthesized using random hexamer primers, the second strand was synthesized using appropriate buffers, dNTPs, RNase H, and DNA polymerase I. After purification using the QiaQuick PCR extraction kit, the double-stranded cDNAs were resolved using EB buffer for end reparation and poly (A) addition. We then connected the short fragments with sequencing adapters. Finally, a cDNA library was constructed, which was sequenced using an IlluminaHiSeq 2000. Low-quality reads were removed using an in-house Perl script. The transcriptome sequences were then assembled with clean reads. First, the distinct contigs were assembled using the short reads with the software program CLC Genomic Workbench 5.5. Next, scaffolds between the contigs were constructed by employing the paired-end relationships between the reads. Finally, we filled the intra-scaffold gaps and constructed a non-redundant unigene set from the assembled datasets using CAP3 software [26].

Statistical analyses
The expression levels of the unigenes were calculated by RNA-Seq quantification analysis as the number of reads per kilobase of an exon region in a given gene per million mapped reads (RPKM) [28]. The transcriptomes of the three samples were used as a reference to screen and analyze the differentially expressed unigenes. A rigorous algorithm was used to identify DEGs using the method of Audic et al. (1997) [29]. The false discovery rate (FDR) was used to determine the threshold of the P-value in multiple tests and analyses [30]. A FDR of < 0.05 and an absolute value of log2 (ratio) 2 were the thresholds used to define significant differences in gene expression [28]. The DEGs with minimum two-fold changes in expression were subjected to differential gene expression analysis.

Real-time quantitative PCR verification by qRT-PCR
Total RNA was extracted from the vegetative buds, floral buds and buds as described above. After treatment with DNase, the total RNA was subjected to reverse transcription to produce cDNA using a reverse transcription system. Real-time RT-PCR was conducted using the PikoReal real-time PCR system. Each reaction was performed in a total reaction mixture volume of 20 μL, which contained 2 μL of first-strand cDNA as template. The amplification program was: 30 s at 95°C, 40 cycles of 5 s at 95°C and 30 s at 60°C. Gene-specific primers (S12 Table) were designed to provide a relative quantification of each gene. All of the real-time RT-PCR experiments were conducted using three biological replicates, and all replicates were measured in triplicate. We analyzed the relative expression levels using the 2-ΔΔCt method, and the protein phosphatase 2A (PP2Acs) gene of C. morifolium was used as a reference gene [31].
Supporting Information S1

Author Contributions
Conceived and designed the experiments: HL MS QZ. Performed the experiments: HL. Analyzed the data: HL. Contributed reagents/materials/analysis tools: HL DD HP TC JW QZ. Wrote the paper: HL.