Whole-transcriptome analysis of differentially expressed genes in the ray florets and disc florets of Chrysanthemum morifolium

Chrysanthemum morifolium is one of the most important global cut flower and pot plants, and has been cultivated worldwide. However, limited genomic resources are available and the molecular mechanisms involved in the two morphologically distinct floret developmental cycles in chrysanthemum remain unclear. The transcriptomes of chrysanthemum ray florets, disc florets and leaves were sequenced using Illumina paired-end sequencing technology. In total, 16.9 G reads were assembled into 93,138 unigenes with an average length of 738 bp, of which 44,364 unigenes showed similarity to known proteins in the Swissprot or NCBI non-redundant protein databases. Additionally, 26,320, 22,304 and 13,949 unigenes were assigned to 54 gene ontology (GO) categories, 25 EuKaryotic Orthologous Groups (KOG) categories, and 280 Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways, respectively. A total of 1863 differentially expressed genes (DEGs) (1210 up-regulated and 653 down-regulated) were identified between ray florets and disc florets, including genes encoding transcription factors and protein kinases. GO and KEGG pathway enrichment analyses were performed on the DEGs to identify differences in the biological processes and pathways between ray florets and disc florets. The important regulatory genes controlling flower development and flower organ determination, as well as important functional genes in the anthocyanin biosynthetic pathway, were identified, of which two leucoanthocyanidin dioxygenase-encoding genes showed specific expression in ray florets. Lastly, reverse transcription quantitative PCR was conducted to validate the DEGs identified in our study. Comparative transcriptome analysis revealed significant differences in patterns of gene expression and signaling pathways between ray florets and disc florets in Chrysanthemum morifolium. This study provided the first step to understanding the molecular mechanism of the differential development of ray florets and disc florets in chrysanthemum, and also provided valuable genomic resources for candidate genes applicable for the breeding of novel varieties in chrysanthemum.


Background
Chrysanthemum morifolium is one of the most important global cut flower and pot plants, and has been cultivated worldwide [1]. Breeders have attempted to breed more chrysanthemum varieties with novel flower types and colors through cross breeding and transgenic biotechnology [1]. Capitulum is the main ornamental part of C. morifolium, and the typical inflorescence structure of chrysanthemum is composed of two morphologically distinct florets; namely, ray florets and disc florets. Ray florets are ligulate and zygomorphic, with a showy corolla (petals) and aborted stamens, and play an important role in attracting pollinators. The central disc florets with fertile pollens are radially symmetrical and hermaphroditic, and are used for reproduction in chrysanthemum (Additional file 1). However, the molecular mechanism involved in the development of two morphologically distinct florets has not been well-characterized to date due to the lack of genomic information.
In Arabidopsis, significant progress has been made toward understanding the molecular mechanisms involved in flower development [2]. ABCE models have revealed A-class genes that specify sepal identity, A-class genes together with B-class genes that specify petals, B-class genes together with C-class genes that determine stamens, C-class genes that determines carpel identity, and E-class genes that were thought to pre-establish the floral context so that other organ identity genes could function [3,4]. Many MADS-Box genes including AP1, AP3, PI, and AG, as well as AP2 genes, were shown to be important regulators involved in flower organ specification [4,5]. In chrysanthemum, only a few transcription factors in the flower developmental pathway were isolated and analyzed, such as homologs of APETALA1, SEPALLATA3, and FRUITFULL [6].
Flavonoids are ubiquitous secondary metabolites that have important functions in plant physiology [7]. Anthocyanin is a class of flavonoids responsible for the orange-toblue colors identified in many flowers, leaves, fruits, seeds, and other tissues, and it can attract pollinators and protect against damage from UV irradiation [8]. The production of anthocyanins provide attractive colors for chrysanthemums and improve the ornamental value of chrysanthemum, especially for cut flowers that require bright colors to decorate flower bouquets and baskets. As another class of abundant flavonoids, flavonols also have diversified functions in the plant kingdom, e.g. as effective UV filters, as regulators of auxin transport, as phytoalexins, or as pigments in flowers and fruit [9,10]. Due to the significant differences in corollaceous colors between ray florets and disc florets of chrysanthemum, it is necessary to perform the comparison of pigments between them. And the molecular mechanism for different flower colors of ray florets and disc florets should be further explored.
In addition, transcriptome profiling of chrysanthemum ray florets was performed using RNA-Seq to investigate light-induced anthocyanin biosynthesis [11]. Interestingly, in chrysanthemum, the ray florets of many varieties exhibit bright colors, such as red, pink, and purple. However, the disc florets show only yellow pollen. The molecular mechanism controlling the flower color difference between ray florets and disc florets remains unknown.
Transcription expression profiling, data assembly, and analysis increase our understanding of the gene regulation networks and biological pathways, can reveal genes downstream of key transcription factors involved in related pathways or networks, and may ultimately explain specific processes [12]. The use of RNA-Seq in studies of chrysanthemum and related species has been demonstrated [11,[13][14][15][16][17]. However, at this time no comparative transcriptome information on the ray florets and disc florets of chrysanthemum has been reported. Here, transcriptional sequencing and comparative analysis of chrysanthemum ray florets, disc florets, and leaves were performed using Illumina assembly technology and RNA-Seq quantification analysis. Based on transcriptional sequencing and analysis, we identified DEGs between ray florets and disc florets to reveal important regulators controlling differential development of ray and disc florets and novel genes involved in regulating stamen development in disc florets. Next, we identified important regulatory genes involved in controlling flower development and flower organ identification, as well as important functional genes in the anthocyanin synthesis pathway to create a list of candidate genes for functional analysis of flowering regulation in chrysanthemum. Our study provides valuable genomic resources and candidate genes for the breeding of novel varieties, and increases our understanding of the developmental molecular mechanisms of both ray florets and disc florets in chrysanthemum.

Illumina sequencing and assembly
The cDNA libraries of ray florets, disc florets and leaves were sequenced using the IlluminaHiSeq™ 2000 system. After stringent quality checks and data cleaning, 46,952,994, 51,333,542 and 68,328,428 clean reads were generated for ray florets, disc florets and leaves, respectively. The average proportion of clean reads for the library was 96.2 %.
A total of 124,264 contigs were assembled based on the high-quality reads, with a total size of 68,266,284 bp, an N50 of 648 bp and an average contig length of 549 bp. We then constructed scaffolds between the contigs via the paired-end relationships between the reads. A total of 96,942 scaffolds were obtained, with an N50 of 819 bp and an average length of 721 bp. We filled the intra-scaffold gaps and constructed a non-redundant unigene set from three of the assembled datasets using CAP3 software. Finally, a total of 93,138 high-quality unigene sequences were obtained with an average length of 738 bp (Table 1).

Gene annotation and functional classification
Of the 93,138 unigenes, 44,364 (47.63 % of the total) were aligned to the Nr protain database (date 2014.03) and Swissprot protein database (date 2014.03) using an E-value threshold of < 1e-5, which meant that 48,774 unigenes (52.37 % of the total) had no Swiss-Prot annotation because of the lack of chrysanthemum genome and EST information ( Table 1). The average length of the sequences with significant matches was 927 bp, while that of those without matches was 567 bp. The match efficiency was 89.63 % for unigenes longer than 2000 bp. However, the match efficiency decreased to 31.97 % for those between 200 and 600 bp. Therefore, the sequence lengths affected the annotation.
Gene Ontology (GO) consisted of three ontologies; molecular functions, cellular components, and biological processes, which facilitated gene annotation and analysis. A total of 26,320 unigenes were classified into 54 functional categories and 19,004 GO terms using the Blast 2GO software, among which 7201 unigenes were assigned to 25 GO categories and 10,282 terms in biochemical processes, 13,285 unigenes were assigned to 13 categories and 5890 terms in molecular functions and 5834 genes were assigned to 16 categories and 2832 terms in cellular component. As shown in Fig. 1, in 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, respectively, indicating that numerous metabolic activities were activated during the development of ray and disc florets, which were regulated by a wide range of genes that interacted within cell parts. In addition, we observed a high percentage of genes from the 'cellular process' , 'membrane' , and 'catalytic activity' categories, but few from 'cell killing' , 'cell junction' , and 'metallochaperone activity' in each of the three main GO classifications (Fig. 1).
Within the 'biological process' category, we identified 25 terms related to flower development in level 6, including 'flower development' (GO: 0009908), 'floral organ development' (GO: 0048437), 'regulation of flower development' (GO: 0009909), 'floral whorl development' (GO: 0048438), and 'photoperiodism, flowering' (GO: 0048573) ( Table 2). In total, 154 unigenes were assigned to these 25 terms, of which 116 were annotated as uncharacterized or predicted proteins and the remainder showed homology to F-box family genes, SET domain protein -coding genes, and GIGANTEA (Additional file 2). In addition, as shown in Table 3, we identified 19 terms that were related to the pollen-pistil interaction, recognition of pollen, pollen development, and stamen development in the 'biological process' category. In addition, 295 unigenes were assigned to these 19 terms, of which 166 were annotated as uncharacterized or predicted proteins and the remainder showed homology to serine/threonine-protein kinase and ARK3-like protein coding genes (Additional file 3). These genes assigned to the terms related to flower and stamen development were important candidate genes that may regulate flower development, specification of floral organ identity, and stamen development in chrysanthemum.
The Kyoto Encyclopedia of Genes and Genomes (KEGG) Pathway is a collection of manually drawn pathway maps representing the networks of molecular interactions in the cells and variants of these pathways specific to particular organisms. To further analyze the transcriptome of chrysanthemum, all unigenes were compared with the KEGG database using BLASTx with an E-value threshold of < 1e-5.
Of the 44,364 unigenes, 13,949 had significant matches with at least one KEGG pathway in the database and were assigned to 280 KEGG pathways ( Table 4). The most represented pathways were 'metabolic pathways' (3475 members) and 'biosynthesis of secondary metabolites' (1887 members), followed by 'microbial metabolism in diverse environments' (741 members), 'biosynthesis of amino acids' (482 members), and 'carbon metabolism' (471 members). In addition, 381 unigenes were assigned to "plant hormone signal transduction".

Comparison of the transcriptomes of ray florets and disc florets
The set of unigenes common to disc florets and leaves The number of unigenes with an RPKM value of > 0.3 that were shared by ray florets and disc florets were 82,017. (Fig. 3). By contrast, 2414 and 5293 unigenes showed specific expression in ray florets and disc florets, respectively. Therefore, more genes were expressed in disc florets than in ray florets, since stamens that were aborted in ray florets developed normally in disc florets.

DEGs between ray florets and disc florets
Transcriptomes of ray florets and disc florets were compared, and we mapped the resulting reads to the reference transcriptome. Statistical significance was considered reliable with an RPKM value ≥ 2 in at least one of the three samples. DEGs were filtered with an FDR ≤ 0.05 and a |log2 (ratio)| ≥ 2; only DEGs with a minimum of a two-fold change in expression were used in the analysis of differential gene expression. A total of 1863 DEGs (1210 up-regulated in disc florets and 653 down-regulated in disc florets) were identified between ray and disc florets (Fig. 4a). The correlation of gene expression between ray florets and disc florets was examined using an algorithm developed from the correlation scatter plot (Fig. 4b).
Among the 653 up-regulated unigenes in ray florets relative to disc florets, 247 were annotated to characterized proteins, including many important transcription factors factor-coding genes, functional genes in the anthocyanin biosynthetic pathway, and genes related to pigments or aromatic constituents synthesis (Additional file 4). Of the 1210 genes up-regulated in disc florets relative to ray florets, 387 were annotated to characterized proteins, including some important transcription factors, and several protein kinases (Additional file 5). A total of 132 DEGs were specifically expressed in disc florets relative to ray florets, including WUS and pectinesterase-coding genes, which may be involved in the regulation of stamen development. The annotation information on the DEGs specifically expressed in disc florets was provided in Additional file 6.
GO and KEGG pathway enrichment analyses were performed on the DEGs to identify differences in biological processes and pathways between ray florets and disc florets. In total, 990 DEGs were enriched in 697 GO categories. In the 'biological process' category, the dominant terms were the following: 'oxidation-reduction process' (GO: 0055114), 'metabolic process' (GO: 0008152), and 'carbohydrate metabolic process' (GO: 0005975). In the 'molecular function' category, the most representative terms were the following: 'hydrolase activity' (GO: 0016787), 'oxidoreductase activity' (GO: 0016491), and 'nucleotide binding' (GO: 0000166). Finally, in the 'cellular component' category, the most representative terms were the following: 'membrane' (GO: 0016020), 'integral to membrane' (GO: 0016021), and 'cell wall' (GO: 0005618). Thus, the physiological and biochemical activities involved in oxidation-reduction processes, metabolic processes, and carbohydrate metabolic processes differed between ray florets and disc florets. A total of 399 DEGs were enriched in 204 KEGG pathways. The dominant pathways were the following: 'metabolic pathways' (ko01100), 'biosynthesis of secondary metabolites' (ko01110), and 'starch and sucrose metabolism' (ko00500), which indicated that there were significant differences in physiological and biochemical processes involved in metabolic pathways, biosynthesis of secondary metabolites, as well as starch and sucrose metabolism between ray florets and disc florets. Interestingly, for many of the enriched KEGG pathways, all DEGs were up-regulated in disc florets relative to ray florets, as was the case for 'phagosome' (ko04145), 'arginine and proline metabolism' (ko00330), 'progesterone-mediated oocyte maturation' (ko04914), and many signaling pathways such as 'oxytocin signaling pathway' (ko04921), 'PI3K-Akt signaling pathway' (ko04151), 'Rap1 signaling pathway' (ko04015), 'neurotrophin signaling pathway' (ko04722), and 'Hippo signaling pathway' (ko04390) ( Table 5). Therefore, these signaling pathways may play important roles in disc floret development. The enriched GO terms and KEGG pathways are listed in Additional files 7 and 8.

Important transcription factors differentially expressed between ray florets and disc florets
A total of 19 important transcription factors family genes were dramatically differentially expressed between ray florets and disc florets. As shown in Table 6, 13 transcription factor family genes were significantly up-regulated in ray florets relative to disc florets, including nine TCP family members (four TCP-like genes and five CYC-like genes) and four other transcription factor family genes (CRT, AP2/ ERF, BHLH and DOF). And three CYC-like genes (CYC2CM2, CYC2CM4, CYC3CM1) showed more than 14-fold higher transcription levels in ray florets compared with disc florets, which indicated that these three CYC-like genes may play important roles during the development of ray florets. By contrast, six transcription factor family genes were remarkably up-regulated in the disc florets relative to ray florets, including MADS-box, WUS, NAC and MYB family genes (Table 6). Interestingly, one WUS family gene (WUSCL) was specially expressed in the disc florets. These transcription factors may have important functions in the differential development of ray and disc flowers.

Identification of important genes controlling flower development and organ determination in Chrysanthemum
Flowering is a complex process controlled by environmental conditions and developmental regulation. Several signaling pathways involved in flowering control have been identified in Arabidopsis, including the photoperiod pathway, the GA pathway, the vernalization pathway, and the autonomous pathway [18]. These different signaling pathways are known to converge on the regulation of the same floral development genes. As shown in Fig. 5, the gene regulatory networks involved in the different signaling pathways in Arabidopsis were outlined, and the homologs of these important regulators were identified in chrysanthemum (Additional file 9).  In most core eudicot species, B class genes include three different lineages: PI, euAP3 and TM6. However, TM6 like genes seem to have been lost in Arabidopsis and Antirrhinum [20]. In this study, the homologs of PI (PICL1, PICL2, GLOCL), euAP3 (AP3CL1, AP3CL2) and TM6 (TM6CL) were identified in C. morifolium. We identified one homolog of C-class gene AG (AGCL) and two homologs of E class gene SEP (SEPCL1, SEPCL2) in chrysanthemum. WUSCHEL (WUS) is known to maintain the stem cell activity. In the floral meristems, WUS binds to adjacent sites of LFY in the AG regulatory region, and both WUS and LFY promote the up-regulation of AG [2]. KNUCKLES (KNU) encoded a C2H2 zinc finger protein, which is required to repress the expression of WUS in the floral meristem, and WUS expression disappears by stage 6 in the process of normal floral development [2]. The expression of KNUCKLES (KNU) is induced by AG. One homolog of WUS (WUSCL) and one homolog of KNU (KNUCL) were identified. We found that WUSCL was specifically expressed in disc florets, which may be important for the differential development of ray and disc flowers. In addition, several other genes known to participate in controlling floral development in chrysanthemum were   identified. Annotation information on the unigenes controlling floral development and organ determination in chrysanthemum is provided in Additional file 9.

Identification of important functional genes in the anthocyanin biosynthetic pathway and pigments in chrysanthemum
To explore the molecular basis of the difference in flower colors between ray florets and disc florets, we identified important functional genes involved in the chrysanthemum anthocyanin biosynthetic pathway and analyzed the expression differences of these genes between ray florets and disc florets (Fig. 6). As shown in The majority of these genes were significantly upregulated in ray florets relative to disc florets. Only one CHI gene (CHICM1), two F3H genes (F3HCM1 and F3HCM2), and two AT genes (ACTCM1 and ACTCM2) were weakly up-regulated in disc florets relative to ray florets. Interestingly two LDOX genes were specifically expressed in the ray florets, which was verified using qRT-PCR (Additional file 10, Fig. 6). Annotation information for these genes in the anthocyanin biosynthetic pathway is listed in Additional file 10.
The variety of flower colours was determined by the pigment types and contents in flowers. Qualitative analysis of pigments in the flowers of pink, red and purple chrysanthemums was performed by HPLC (Additional file 11). Both anthocyanins and flavonols were detected in the ray florets of pink, red and purple chrysanthemums (Fig. 7). However, no anthocyanins was detected in the disc florets of these flowers above mentioned. (Fig. 7a). As shown in Fig. 7a

Discussion
For the chrysanthemum transcriptome, 93,138 assembly unigenes were obtained, and these unigenes were assigned to a wide range of GO categories, KEGG pathways and KOG classifications, which indicated that a great variety of transcripts are involved in chrysanthemum flower development. We identified unigenes that were annotated to the GO terms related to flower development and stamen development in the 'biological process' category, of which many members encoded uncharacterized proteins. These provided important candidate genes for the discovery of novel regulators controlling flower and stamen development in chrysanthemum. Thus, our work provides valuable information and sequence resources for identifying the genes that regulate the development of ray and disc florets in C. morifolium.

Significant differences in gene expression and signaling pathways between ray florets and disc florets form the molecular basis of their morphological and functional distinction
Comparative analysis of the transcriptomes between ray florets and disc florets revealed valuable information of candidate genes and their expression patterns in chrysanthemum. As shown in Fig. 3, ray florets and disc florets shared 82,017 common transcripts, whereas 2414 and 5923 transcripts were specifically expressed in ray florets and disc florets, respectively. More DEGs (1210) were significantly up-regulated in disc florets than (653) in ray florets. Thus, although the majority of unigenes were expressed in ray florets and disc florets, the differences in gene Specific transcription factor genes were predicted to be important regulators during the development of disc florets in chrysanthemum In this study, we identified several transcription factors genes showing significantly up-regulated expression, including MADS-box, NAC, MYB and WUS transcription factor genes. As master regulators of floral organ determination, MADS-box transcription factors have played important roles during the evolution of flowering plants [5]. Thus, the MADS-box genes up-regulated in disc florets (MADSCM1 and SVPCM1) may have important functions in the development of disc florets. As a major group of plant-specific transcription factors, the NAC family genes play important roles in plant development, including seed development, embryo development, shoot apical meristem formation, fiber development, leaf senescence, and cell division etc. [21,22]. We identified two NAC family genes (NACCM1, NACCM2) that showed more than fivefold increased transcript levels in disc florets compared with ray florets, which indicated NACCM1 and NACCM2 were important regulators mediating the development of disc florets. MYB transcription factors involved in the ABA-dependent pathway were known to up-regulate abiotic stress-responsive genes [23]. However, one MYB gene (MYBRCM1) showed up-regulated expression in disc florets relative to ray florets, indicating that MYB transcription factors probably also play a regulatory role during the development of disc florets. Therefore, besides A, B and C function genes, some other transcription factor (NAC, MYB, WUS) genes may also regulate the development of disc florets in chrysanthemum. Furthermore, these transcription factors genes may have important regulatory roles in the differential development of ray and disc florets, which are required to be explored in the future studies.

Many DEGs specifically expressed in disc florets may control pollen development
We identified 132 DEGs that were specifically expressed in disc florets relative to ray florets, among which 30 DEGs were annotated to known proteins, including pectinesterase, copper transporter, metallothionein-like protein, calcium-binding protein Calnexin, profiling, and Amb a 1-like protein, and the others coded proteins of unknown function (Additional file 6). In Arabidopsis, pectinesterase genes are required for pollen separation during normal floral development, which is expressed in anther tissues shortly after meiosis is complete [24]. Thus, it was shown that pectinesterase-coding genes specifically expressed in disc florets may play a vital role in pollen development in chrysanthemum. In Arabidopsis, Copper Transporter COPT1 is expressed in pollen and COPT1 antisense plants showed pollen morphological abnormalities, which suggests that COPT1 is a copper transporter involved in pollen development, and pollen development is one of the most sensitive processes to copper depletion throughout the entire plant [25]. Thus, the classically described symptoms of male sterility caused by pollen defects in plants grown in copper-deficient soils have been accepted [25]. In our study, one copper transporter gene was expressed specifically in disc florets, which indicated that copper transporter genes might be involved in chrysanthemum pollen development. Mmetallothioneins (MTs) consist of a type of low-Mr, Cys-rich proteins, which bind heavy metals and are widely distributed in eukaryotic and prokaryotic organisms [26]. In this study, MTs were predicted to be involved in related signal transduction pathways during the development of stamens in disc florets of chrysanthemums, along with many other DEGs specifically expressed in the disc florets including calciumbinding protein, calnexin, and profiling-coding genes. In addition, three types of pollen allergen, calcium-binding protein, profiling, and Amb a 1-like protein show specific expression in the disc florets, which provides insight into the allergenic potential of chrysanthemum flowers [27].

TCP-like genes may mediate the development of ligulate and zygomorphic corolla of ray florets in chrysanthemum
Transcription factors play important roles in the regulation of plant growth and development. As plant-specific transcription factors, TCP transcription factors control cell cycle in angiosperms, resulting in differential growth at meristems and in individual organs; and TCP transcription factors have been proved to regulate morphological traits during plant evolution, including the complex architecture of Asteraceae inflorescence, composed of various types of flowers [28,29]. Studies in some core eudicotssuch as Antirrhinum majus (Plantaginaceae), Helianthus annuus, Gerbera hybrida and Senecio vulgaris (Asteraceae) -have shown that CYC2-clade TCPs play prominent roles in establishing flower symmetry by specifying identity to the dorsal (adaxial) region of the flower; and these genes are also involved in the control of flower-type differences in H. annuus, G. hybrida and S. vulgaris of Asteraceae [20,[30][31][32][33]. Thus far, little information is available regarding the TCP-like genes in chrysanthemum.
In this study, we identified nine TCP transcription factor family genes that showed higher expression in ray florets relative to disc florets; in particular, the CYC-like genes were significantly up-regulated in ray florets. This result was consistent with studies on CYCLOIDEA/ TEOSINTE BRANCHED1 (CYC/TB1) gene family in G. hybrida and H. annuus, which indicated that all sunflower and gerbera CYC2 clade genes showed upregulated expression in marginal ray (and trans) flowers.
In Asteraceae, previous studies concerning flower symmetry were restricted to CYC2-clade TCPs. However, we also identified other TCP-like genes with up-regulated expression in ray florets of chrysanthemums. Therefore, 1) the up-regulated expression of CYC-like genes contributes to the development of ligulate and zygomorphic corolla in ray florets of chrysanthemum; 2) in addition to CYC2-clade genes, other TCP-like genes may also play a role in the control of flower-type differences in chrysanthemum. Further studies are required to explore the exact roles of these TCP-like genes in the regulation of ray floret development in chrysanthemum.

A new version of the ABC model for chrysanthemum is required
In the present study, we outlined the gene regulatory networks controlling flowering involved in the different signaling pathways in Arabidopsis and identified the homologs of these important regulators involved in floral meristem identity and the regulation of flower development, including AP2 and many MADS box genes (Fig. 5). The classic ABC model was developed based on floral mutant phenotypes in Arabidopsis and Antirrhinum, which postulates that three different gene classes (or functions) specify the identity of the flower organs in the four floral whorls and later E-function genes required for all other functions were added to the model [5]. The homologs of A-, B-, C-, and E-class genes were identified based on the protein annotation of unigenes in chrysanthemum, including one homolog of the TM6 lineage (Bclass genes), which is not found in Arabidopsis and Antirrhinum. However, different from the traditional model plants, the capitulum of chrysanthemum is composed of two morphologically distinct florets (ray florets and disc florets). Thus, further studies are required to explore the regulatory mechanism of A, B, and C function genes in chrysanthemum. Studies in gerbera have shown that 1) flower-type specific MADS-box protein complexes play a central role in differential development of ray and disc flowers; 2) B class MADS-box genes; GGLO1 (PI/GLOlike gene) and GDEF2 (euAP3-like gene), mediate the classical B-function since they determine petal and stamen identities; 3) GDEF1 (B-class gene of TM6 lineages) is not involved in determining petal identity, and it may independently regulate stamen development of late stages [20,34]. In chrysanthemum, CDM111 (A function gene) is the functional equivalent to AP1, and CDM44 (E function gene) is the functional equivalent to SEP3 [6]. However, further studies are required to explore the functions of A-, B-, C-, and E-class genes in flower organ determination and differential development of ray and disc florets in chrysanthemum. Furthermore, we identified a further two regulators of flower development together with MADS-box genes; WUS (WUSCL) and KNU (KNUCL), whose function in regulation of flower development in Asteraceae remains unclear. WUS, which binds to the adjacent sites of LFY in the AG regulatory region, promotes the up-regulation of AG [2]. KNU functions as an inhibitor of WUS in the floral meristem. Interestingly, we found that WUSCL was specifically expressed in the disc florets. Recently, it has been illuminated that the prototypic WOX-family member, the WUS gene, is a bifunctional transcription factor acting as a repressor in stem-cell regulation and as an activator in floral patterning [35]. The WUS gene is sufficient to return differentiating cells to stem cells to maintain stem cells, and its expression is repressed at stage 6 in an AGdependent manner [36]. Thus, in our study it was predicted that 1) the specific expression of WUS in the disc floret region contributes to formation of the indefinite inflorescence in chrysanthemum; 2) WUS is involved in flower organ determination together with AG. Further studies are required to explore the functions of WUS in flower development and organ determination in chrysanthemum.
Anthocyanidins are not synthesized in disc florets as two LDOX genes are not expressed We analyzed important genes in the chrysanthemum anthocyanin biosynthetic pathway and identified two LDOX genes that are highly expressed in ray florets but are not expressed in disc florets. Leucoanthocyanidindioxygenase (LDOX) is also known as anthocyanidinsynthase (ANS), belongs to the OGD family, and catalyzes the synthesis of corresponding colored anthocyanidins. High expression of LDOX genes in the ray florets promoted anthocyanidin synthesis and accumulation. The qualitative analysis of pigments in the flower of chrysanthemum confirmed that anthocyanins are present only in the ray florets, while the disc florets contained no anthocyanidins, which was consistent with LODX-like gene specific expression in ray florets. Anthocyanins in the ray florets provide vibrant colors as visible signals to attract insects for pollination.

Conclusions
In our study, comparative transcriptome analysis revealed significant differences in gene expression and signaling pathways between ray florets and disc florets in C. morifolium. We identified unigenes that were annotated to the GO terms related to flower development and stamen development in the 'biological process' category, which provided important candidate genes for the discovery of novel regulators controlling flower and stamen development in chrysanthemum. We identified homologs of the important regulators of flower development and organ determination, including the homologs of A-, B-, C-, and E-class genes as well as two regulators controlling flower development together with MADSbox genes: WUS (WUSCL) and KNU (KNUCL). In addition, WUSCL, which is specifically expressed in the disc florets, may play important roles in flower development and organ determination. The important functional genes in the anthocyanin synthesis pathway were analyzed, and we identified two LDOX genes that showed specific expression in ray florets, which explains why anthocyanidins were not synthesized in disc florets. These findings indicate that the differential development of two morphologically distinct florets (ray florets and disc florets) is a complex biological process regulated by a wide range of factors, and further studies are required to explore the gene regulatory networks involved in this process. This study represents the first step in exploring the molecular mechanism of the differential development of ray florets and disc florets in chrysanthemum, and it provided valuable genomic information and candidate genes for the breeding of novel chrysanthemum varieties.

Plant materials and RNA extraction
The tissues (ray florets, disc florets and leaves) used in the study were obtained from a ground-cover chrysanthemum variety (C. morifolium 'Fenditan' , a hybrid of chrysanthemum varieties) cultivated in a greenhouse under an 8-h light/16-h dark cycle at 23°C in Beijing Forestry University (116.3°E, 40.0°N). After marginal ray florets (central disc florets) of the flower buds were fully formed, about 60-120 marginal ray florets (central disc florets) and 3-4 fully expanded leaves were collected between 9:00-12:00 a.m. once per week until the flowers were in full bloom. The collected plant materials were placed immediately in liquid N2 and stored at −80°C until RNA extraction. Total RNA was extracted using the RNeasy Plant Mini Kit (Qiagen, Cat. #74904). RNA quantity and quality were assessed using a NanoDrop ND2000 instrument (Thermo Scientific).

Illumina sequencing, and de novo assembly functional annotation
Ray florets and disc floret development are continuous developmental processes. Therefore, to obtain complete transcriptome information, we pooled equivalent quantities of total RNA isolated from different developmental stages of ray florets and disc florets, and the pooled total RNA samples were used for sequencing. Thus, leaves were collected at the same time as the samples of ray florets and disc florets. In total, 10 μg of total RNA were collected for each sample for sequencing. We conducted llumina sequencing using an Illumina HiSeq™ 2000 system (Illumina, San Diego, CA, USA) in Shanghai (ShoBiotechnology Corporation (SBC), Shanghai, China). After poly (A) mRNA was purified and fragmented into small pieces, the first strand cDNAs were synthesized using random hexamers primers, after which the second strand was synthesized. Double-strand cDNAs were purified with Qia-Quick PCR extraction kit (Qiagen) and resolved using EB buffer for end reparation and addition of a poly (A) tail. The short fragments were then connected with sequencing adapters. Briefly, a cDNA library was constructed with average insert sizes of 300-500 bp. We conducted cDNA sequencing using the Illumina HiSeq™ 2000 system according to the manufacturer's protocols, with paired end 2 × 100 nt multiplex. After removing the low-quality reads, the transcriptome sequences were assembled into distinct contigs using the short reads with the software CLC Genomic Workbench 5.5 software (CLC Bio, Denmark). Scaffolds were then constructed between contigs employing the paired-end relationships between the reads. The intra-scaffold gaps were filled and a nonredundant unigene set was constructed from all three assembled datasets using CAP3 software [37].

Analysis of chrysanthemum transcriptome sequencing results
The expression level of each unigene was calculated using RNA-Seq quantification analysis as the number of reads per kilobase of exon model per Million mapped reads (RPKM) [39]. The chrysanthemum transcriptome from the three samples was employed as reference for the screening and analysis of differentially expressed unigenes due to the unavailability of existing data. A rigorous algorithm was used to identify differentially expressed genes based on the method of Audic et al. [40]. We used the false discovery rate (FDR) to confirm the threshold of the P-value in multiple tests and analyses [41]. An FDR of < 0.05 and the absolute value of log2 (ratio) ≥ 2 were used as thresholds to define significant differences in gene expression [39]. Only the differentially expressed genes (DEGs) with a minimum of a twofold change in expression were used in the analysis of gene differential expression.

Gene expression analysis based on qRT-PCR
Total RNA was extracted from the ray florets, disc florets, and leaves as described above. Total RNA was treated with DNase (Promega, USA), and then subjected