Transcriptomic Analysis of Flower Development in Wintersweet (Chimonanthus praecox)

Wintersweet (Chimonanthus praecox) is familiar as a garden plant and woody ornamental flower. On account of its unique flowering time and strong fragrance, it has a high ornamental and economic value. Despite a long history of human cultivation, our understanding of wintersweet genetics and molecular biology remains scant, reflecting a lack of basic genomic and transcriptomic data. In this study, we assembled three cDNA libraries, from three successive stages in flower development, designated as the flower bud with displayed petal, open flower and senescing flower stages. Using the Illumina RNA-Seq method, we obtained 21,412,928, 26,950,404, 24,912,954 qualified Illumina reads, respectively, for the three successive stages. The pooled reads from all three libraries were then assembled into 106,995 transcripts, 51,793 of which were annotated in the NCBI non-redundant protein database. Of these annotated sequences, 32,649 and 21,893 transcripts were assigned to gene ontology categories and clusters of orthologous groups, respectively. We could map 15,587 transcripts onto 312 pathways using the Kyoto Encyclopedia of Genes and Genomes pathway database. Based on these transcriptomic data, we obtained a large number of candidate genes that were differentially expressed at the open flower and senescing flower stages. An analysis of differentially expressed genes involved in plant hormone signal transduction pathways indicated that although flower opening and senescence may be independent of the ethylene signaling pathway in wintersweet, salicylic acid may be involved in the regulation of flower senescence. We also succeeded in isolating key genes of floral scent biosynthesis and proposed a biosynthetic pathway for monoterpenes and sesquiterpenes in wintersweet flowers, based on the annotated sequences. This comprehensive transcriptomic analysis presents fundamental information on the genes and pathways which are involved in flower development in wintersweet. And our data provided a useful database for further research of wintersweet and other Calycanthaceae family plants.


Introduction
The small, evolutionarily ancient Calycanthaceae family comprises four genera, namely Calycanthus L., in North America, Idiospermum Blake, in Australia, and Sinocalycanthus Cheng & S. Y. Chang and Chimonanthus L., in China [1]. Wintersweet (Chimonanthus praecox), also known as the 'wax shrub', is a hardy, fastgrowing perennial shrub, native to China; it is dichogamous and diploid (2n = 22) [2]. It is an important deciduous aromatic plant, and it is also one of the most precious epibiotic species dating back to the Tertiary period, being classified as a Class II protected wild plant in China [3]. Wintersweet has over 1000 years' history of cultivation and, as its name indicates, it blooms particularly in winter, from late November to March in central southern and south-western China. Its unique flowering time and strong fragrance make it one of the most popular ornamental plants in China; it is appreciated as a pot plant and for cut flowers and it has a high ornamental and economic value. It has also been introduced into Korea, Japan, Europe, America, and Australia [2]. Wintersweet flowers are used in traditional Chinese medicinal preparations to treat heatstroke, vomiting, coughs and measles [4].
They have also been recognized as the source of a natural essential oil, which can be used in perfumery, cosmetics and aromatherapy [5,6]. Recently, these applications of wintersweet have received much attention in New Zealand [7].
There are several traits and properties of wintersweet flowers that are important in a commercial context. These include flower development, senescence, scent biosynthesis and emission, and resistance to biotic and abiotic stresses. Wintersweet blooms especially in winter and has a strong fragrance, and its molecular mechanisms of flower development may therefore be different from those of model species or of plants that bloom in the spring. The molecular and genetic processes that determine some of these flowering characteristics cannot be studied using model species such as Arabidopsis thaliana, or at least only to a limited extent. During the past decade, the genes and gene networks associated with important floral traits, including flower scent, color, morphology and senescence, have been identified and functionally characterized in several plants grown for their flowers, such as rose [8], carnation [9,10,11], and Antirrhinum majus [12,13]. However, the genomic resources that are available for wintersweet or for other members of the Calycanthaceae family are scant and the transcriptional changes and molecular mechanisms that control these important traits and developmental processes in wintersweet flowers are still far from being elucidated. This hampers gene discovery and seriously hinders the improvement of wintersweet as a commercially important species.
This paucity of genomic information is indicated by the fact that although a normalized cDNA library has been constructed from flowers of wintersweet and 479 unigenes have been assembled [14], nevertheless as of late October 2013 only 867 expressed sequence tags (ESTs) for Chimonanthus praecox were available on the NCBI website (http://www.ncbi.nlm.nih.gov/). Because gene identification and characterization in Chimonanthus praecox remain so limited, the large-scale discovery and characterization of functional genes via genome sequencing or exploration of the transcriptome is essential. The highthroughput capacity of the latest generation of RNA sequencing (RNA-Seq) technology provides a unique opportunity for genomic exploration and gene discovery in non-model plant species for which there are no reference genome sequence data [15,16]. The results of RNA-Seq also show high levels of reproducibility, for both technical and biological replicates [17,18]. RNA-Seq generates absolute rather than relative gene expression measurements, providing greater insight and accuracy than microarrays [19,20].
In this study, we used the Illumina RNA-Seq method to analyze the transcriptome of wintersweet flowers based on three cDNA libraries from different stages of flower development. A total of 106,995 transcripts were identified, and the sequences were annotated against the NR database using BLASTX. By using Illumina's digital gene expression platform, we investigated differential gene expression in open and in senescing flowers, and analyzed the selective expression of Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways. In particular, we identified genes of plant hormone signal transduction pathways that were differentially expressed during flower development, and also key genes of floral scent biosynthesis. To our knowledge, this is the first comprehensive transcriptomic study to identify genes and pathways that are differentially expressed during flower opening and senescence in wintersweet. It provides an important new bioinformatics resource for the further identification of genes and gene function involved in flower development in this species.

Illumina Sequencing and the Assembly of Sequence Reads
In this study, we divided wintersweet flower development into three essentially distinct stages: i) bud stage, with the flower buds enlarged and turned yellow and with a displayed petal (DP); ii) open flower stage (OF), with the flowers fully opened and emitting a strong fragrance; and iii) senescing flower stage (SF), with withering petals and loss of fragrance ( Figure 1). Three individual flowers were used to prepare one pooled RNA sample for each of the three developmental stages (DP, OF and SF). Three cDNA libraries were then constructed and subjected to Illumina deep sequencing. We obtained 21,412,928, 26,950,404, and 24,912,954 qualified Illumina reads, respectively, for the DP, OF and SF stages, giving rise to 2,113,119,501, 2,653,944,208 and 2,451,552,573 bp total residues, respectively (Table 1). After eliminating primer and adapter sequences and filtering out the low-quality reads, we pooled together all the high-quality Illumina reads from the three different developmental-stage libraries. Using Trinity software, we then combined all the reads to form a transcriptome database for wintersweet [21]. We identified 106,995 transcripts, with total residues of 127,396,344 bp . The  average length of each transcript was 1,190 bp, the shortest  sequence being 351 bp and the longest being 16,998 bp (Table 1). The sequence length distribution of transcripts is indicated in Figure 2. Most of the transcripts (40.6%) were 401,800 bp in length (Table S1).

Sequence Annotation
For annotation, the sequences were searched against the NCBI non-redundant (NR) database using BlastX, setting a cut-off Evalue of 10 25 . A total of 51,793 (48.4%) sequences showed significant similarity to known proteins in the NR database. The similarity distribution for all the search results showed that 54.7% of the matches were of high similarity, i.e. ranging from 80% to 100% similarity as reported in the BlastX results, whilst 39.2% of the matches were of similarity ranging from 60% to 80% ( Figure 3A). In a further analysis of the matching sequences, we found that 37.7% of the sequences showed closest matches with sequences from Vitis vinifera. The next-closest matches were with sequences from Theobroma cacao and Prunus persica; 12.2% of the sequences showed closest matches with sequences from Theobroma cacao, whilst 7.5% of the sequences showed closest matches with sequences from Prunus persica ( Figure 3B).

GO, COG and KEGG Classification
We added gene ontology (GO) terms using Blast2GO [22], which annotates high-score BLAST matches found to sequences in the NCBI NR proteins database. Genes of wintersweet were classified into three main GO categories and 57 sub-categories. Of the 106,995 assembled transcripts, 32,649 were successfully annotated by GO assignments, and some of them belonged to one or more of the three categories ( Figure 4). Amongst the annotated sequences, biological process categories included cellular process (18,805; 57.6%), metabolic process (18,692; 57.3%), single-organism process ( In a further analysis, the annotated sequences were subjected to a search against the Clusters of Orthologous Group (COG) database, for functional prediction and classification. Based on sequence homology, 21,893 unique sequences were assigned a COG functional classification. These sequences were classified into 25 COG categories, denoting involvement in cellular process, signal transduction, metabolism and other processes ( Figure 5). The most common category was the non-specific category of 'general function prediction only' (4,341; 19.8%), followed by replication, recombination and repair (2,360; 10.8%), transcription (2,267; 10.4%), signal transduction mechanisms (2,022; 9.2%), posttranslational modification, protein turnover and chaperones (1,326; 6.1%), and carbohydrate transport and metabolism (1,138; 5.2%). Amongst other categories represented (Table S3) were translation, ribosomal structure and biogenesis (1,115; 5.1%) and amino acid transport and metabolism (851; 3.9%). The categories of cell motility (11; 0.05%) and nuclear structure (2; 0.01%) were the least represented.
In order to identify biochemical pathways, we mapped the annotated sequences onto the KEGG database, which is an alternative approach to the categorization of gene function and which places emphasis on biochemical pathways. In total, 15,587 transcripts were assigned to 312 KEGG pathways. A summary of the findings is presented in Table S4. The largest number of sequences were those associated with metabolic pathways (4,293; 27.5%), followed by sequences that were involved in the biosynthesis of secondary metabolites (2,049; 13.2%), in microbial metabolism in diverse environments (831; 5.3%), in the biosynthesis of amino acids (495; 3.2%), in plant hormone signal transduction (467; 3.0%), and in a number of other pathways. In particular, several important pathways of secondary metabolite biosynthesis, including the phenylpropanoid, flavonoid and carotenoid biosynthetic pathways, and several signaling pathways, including the p53, mitogen-activated protein kinase (MAPK) signaling pathway, were represented. These pathway assignments provide valuable information for the investigation of specific biochemical and development processes.

Differentially Expressed Genes (DEGs) during Flower Development
To study gene expression during flower development, a genome-wide expression analysis was carried out for each of the DP, OF and SF flower developmental stages ( Figure 1). We mapped the Illumina reads from each stage onto our assembled transcriptome database. The expression of each gene was calculated using the numbers of reads mapping onto its transcripts. Comparisons of gene expression between DP and OF (DP vs OF) and between OF and SF (OF vs SF) showed that 2,706 and 1,466 genes were differentially expressed in DP vs OF and in OF vs SF, respectively (Log 2 fold changes: $2 or #22; FDR #0.05). An overall view of the expression pattern for DP vs OF is depicted in Figure 6A. This shows the expression changes for 2,706 differentially expressed genes (DEGs), ranging from a 14-fold to a -12-fold change in expression. A corresponding view for OF vs SF is depicted in Figure 6B, which shows the changes for 1,466 DEGs, ranging from a 12-fold to an -8-fold change in expression. Figure 6C illustrates that for DP vs OF there were 1,681 up-regulated genes and 1,025 down-regulated genes, whilst in the case of OF vs SF, there were 1,067 and 399, respectively.
To investigate further, the 20 most up-regulated and the 20 most down-regulated genes were each selected from both the DP vs OF and the OF vs SF comparisons. Some of these genes showed matches to sequences in the NCBI NR database, such as cytochrome P450, a-terpineol synthase, lipid-transfer protein, cysteine protease, and GDSL esterase/lipase (Tables S5; S6). For example, we found that a cytochrome P450 gene was the most up-regulated gene (13-fold) at the OF stage but was downregulated substantially at the SF stage (Table S5). Two cysteine protease genes were down-regulated with 210.9 and 28.9 fold changes respectively at the OF stage, and showed no expression at all at the SF stage (Table S5). Finally, we found that a GDSL esterase/lipase gene was induced specifically and strongly (10-fold change in expression) at the SF stage (Table S6).

Pathway Enrichment Analysis of DEGs
The KEGG Orthology Based Annotation System (KOBAS) [23] was used for the further identification of biosynthetic and other pathways and to explore in greater depth the functions of DEGs. Pathways that were found to be significantly upregulated (corrected P-value #0.05) in both DP vs OF and OF vs SF were involved in cutin, suberin and wax biosynthesis;  cyanoamino acid metabolism; pentose and glucuronate interconversions; phenylalanine metabolism; phenylpropanoid biosynthesis; plant hormone signal transduction; and starch and sucrose metabolism ( Table 2). Since wintersweet, otherwise known as 'wax shrub', blooms in winter and displays a strong resistance to biotic and abiotic stress, we focused upon the pathway of cutin, suberin and wax biosynthesis, and obtained matches to key enzymes of cutin, suberin and wax biosynthesis, such as CYP86A1, CYP86A4, CYP86B1, ECERIFERUM (CER), and alcohol-forming fatty acyl-coenzyme A reductases (FARs). The expression of these genes was found to be differentially regulated during flower development in wintersweet. In addition, the CER1 and FAR genes, which are related to wax biosynthesis, were found to be up-regulated at the OF stage and down-regulated at the SF stage ( Figure 7). We also focused upon plant hormone signal transduction pathways, apparently for the first time in wintersweet, almost all the key genes could be identified from our data ( Figure S1). We identified 92 and 72 transcripts associated with plant hormone signal transduction pathways that were differentially expressed in DP vs OF and OF vs SF, respectively ( Table 2). The DEGs were identified as phytohormone receptor or phytohormone-    responsive genes known to be induced by various phytohormone signaling pathways, including those for auxins, cytokinins (CK), abscisic acid (ABA), ethylene (ET), brassinosteroids (BR), jasmonic acid (JA) and salicylic acid (SA), as indicated in Table 3. The early auxin-responsive genes (i.e. those associated with the early phase of a response to auxins) have been grouped broadly into three major classes: the auxin/indole-3-acetic acid (Aux/IAA), the GH3, and the small auxin-up RNA (SAUR) gene families. Our data showed that the Aux/IAA, GH3 and SAUR genes were expressed differentially during wintersweet flower opening and senescence. Moreover, a number of JAresponsive genes, including JAR (JAR Jasmonate Resistant 1), JAZ (Jasmonates ZIM-domain protein) and MYC2, were all induced significantly at the OF stage, whereas the salicylic acidresponsive genes, TGA and PR1, were induced specifically and strongly at the SF stage ( Table 3).

Genes of Floral Scent Biosynthesis
Wintersweet flowers develop a strong and specific fragrance during flower opening. In this study, we obtained from our transcriptomic data the key biosynthetic genes of floral scent production, including GPPS (geranyl diphosphate synthase), MRY (myrcene synthase), GES (geraniol synthase), LIS (S-linalool synthase), FPPS (farnesyl pyrophosphate synthase), TER (aterpineol synthase), and SAMT (S-adenosyl-L-methionine: salicylic acid carboxyl methyltransferase). These included three homologous genes of SAMT and MRY and two homologous genes of TER. Based on the annotated sequences, we have therefore proposed a biosynthetic pathway for the formation of monoterpenes and sesquiterpenes in wintersweet ( Figure 8A).
The expression of genes of floral scent biosynthesis was induced at the OF stage ( Figure 8B). The expression of the LIS gene, which is responsible for a-linalool biosynthesis, was increased 7-fold at the OF stage; and it has been reported that a-linalool accounts for 36% of the total quantity of volatile compounds detected in wintersweet flowers [6]. Similarly, methyl salicylate has also been demonstrated to be one of the major volatile compounds in wintersweet flowers [6]; and the three SAMT gene homologues, responsible for the conversion of salicylic acid to methyl salicylate, were all induced significantly to varying extents at the OF stage.

Real-time Quantitative PCR Validation of RNA-Seq Results
To validate the findings from the sequencing data, an appropriate alternative methodology was chosen. Ten DEGs involved in plant hormone signal pathways were selected for validation using Real-time qPCR, with gene-specific primers designed using Primer Primer software (version 5.0), as shown in Table S7. The expression pattern of these genes at each of the three flower developmental stages is shown in Figure 9. Expression patterns determined by Real-time qPCR were consistent with those obtained by RNA-Seq, confirming the accuracy of the RNA-Seq results reported in this study.

Illumina Sequencing and Sequence Annotation
Wintersweet (Chimonanthus praecox) is mainly a garden plant and a traditional cut flower in China. However, little genomic information is available for this species. In this study, we have generated large amounts of cDNA sequence data, comprising long sequences of good quality, which will facilitate subsequent and more detailed studies. The availability of transcriptomic data for wintersweet will meet the basic requirements needed to undertake a thorough molecular-genetic and biochemical characterization of this species and its relatives.
For gene annotation, the sequences we had determined were searched against the NCBI NR database. Because of the limited genetic information available, no matches were obtained for as many as 55,202 (51.6%) of the sequences. Some of these sequences may represent novel genes in wintersweet, suggesting that flower development in wintersweet might involve some unique processes and pathways. Interestingly, the annotated sequences of wintersweet showed a higher homology to proteins of Vitis vinifera (Figure 3), which may reflect the evolutionary relationship between Chimonanthus praecox and Vitis vinifera. Chimonanthus praecox (Calycanthaceae) belongs to the order Laurales and to the clade magnoliids, which arose in evolution soon after the divergence of the eudicots. Vitis vinifera (Vitaceae) belongs to the clade rosids, which lies at the root of the eudicots, and which is closer phylogenetically to the clade magnoliids than to the order Malvales (e.g. Theobroma cacao), or to the order Rosales (e.g. Prunus persica).
In spite of the large proportion of sequences that showed no matches, a large number of transcripts were nevertheless assigned to a wide range of GO and COG classifications, which indicated that our transcriptome data represented a broad diversity of transcripts in wintersweet. KEGG predictions identified many genes associated with secondary metabolites, including genes of phenylpropanoid biosynthesis, flavonoid biosynthesis, phenylalanine metabolism and carotenoid biosynthesis; in addition, we identified genes related to three alkaloid biosynthesis pathways (Table S4). Alkaloids as a class often exhibit antibacterial and immunostimulatory properties, and this would appear to be consistent with the medicinal and health-promoting properties that are associated with wintersweet [24].

DEGs during Flower Development
Although the analysis of DEGs during flower opening and senescence has been performed using microarrays in Iris [25], Alstroemeria [26], and Rosa [27], nevertheless the availability of data relating to developmental expression profiles and comparing the transcriptomes of flowers or flower parts at different stages remains generally limited. The approach we used here identified .4000 DEGs, at a cut-off of 5% FDR, that showed differential expression patterns in the pairings, DP vs OF and OF vs SF. Further insights were obtained by examining the changes of gene expression for sets of related genes. Thus, amongst the 20 most up-regulated DEGs, we found one cytochrome P450 gene that was up-regulated significantly at the OF stage and down-regulated at the SF stage (Table S5); this is in apparent contrast to a report of increased expression of cytochrome P450 genes in senescing petals of Alstroemeria [26]. Cytochrome P450 genes are found in all kingdoms and show extraordinary diversity in their chemical reactions [28]; they are one of the largest gene families in plants, although the functions of cytochrome P450 genes in flowers are as yet unknown.
In our study, two cysteine protease genes were down-regulated significantly at the OF stage, and showed no expression at all at the SF stage (Table S5). This appears to contrast with a range of observations showing an up-regulation of cysteine proteases during flower senescence, in both ethylene-sensitive and ethyleneinsensitive flowers, such as Dianthus caryophyllus (carnation) [29], Hemerocallis spp. (daylily) [30], Alstroemeria peruviana [31], Sandersonia aurantiaca [32], Narcissus pseudonarcissus (daffodil) [33], and Gladiolus grandiflora [34]. This suggests that there may be a different function, or regulation mechanism, for cysteine protease genes in wintersweet, compared to other species.

Genes Related to Plant Hormone Signaling Pathways
Plant hormones are involved in many different processes throughout the life of a plant, including growth, development and senescence [35]. Several plant hormones such as ET, ABA,  GA have been reported to play important roles in flower opening and in senescence processes. The specific regulation of flower development is a complex process and is governed both by the endogenous levels of the hormones and by the sensitivity of the hormone response mechanisms. The expression of numerous genes has been found to be altered through the action of plant hormone signaling pathways, as well as by the diversity of interactions amongst plant hormones themselves [36].
Ethylene is one of the most important plant hormones, and it plays a major role in many flowers, such as rose, carnation and orchids, in regulating flower opening, senescence and abscission. A burst of endogenously produced ethylene in such flowers initiates senescence and coordinates the expression of genes required for the process [37]. In wintersweet, however, previous research indicated that no obvious production of ethylene could be detected from the flowers during flower development [38]. In the work we Differentially expressed genes with Log 2 fold changes $2 or #22, and FDR #0.05 were included. The genes involved in the pathway are color-coded: green boxes, genes identified in our data; and light blue boxes, genes involved in the pathway present in KEGG database but undetectable in our data. doi:10.1371/journal.pone.0086976.t003 report here, we obtained families of genes encoding ethylene receptors (ETR1, ETR2, EIN4, ERS1 and ERS2), EIN3/EILs (EIN3-like proteins), and the ethylene response factor (ERF); in no cases, however, did we find that the expression of these genes was induced significantly during flower opening and senescence. These findings therefore indicated that flower opening and senescence in wintersweet may not be dependent upon an ethylene-mediated signal transduction process triggered by the endogenous production of small amounts of ethylene. On the other hand, it has been reported that ethephon treatment of wintersweet cut flowers can accelerate the process of flower opening, indicating that wintersweet flowers are sensitive to exogenous ethylene treatment [39]. Salicylic acid (SA) or its derivatives function in a diversity of plant processes, such as seed germination, respiration, stomatal responses, senescence and responses to abiotic stress [40]. Many of the components of the SA signaling pathway, including those of signal perception, have not yet been revealed [41]. However, it is known that SA or its signaling pathway is associated with the activation of diverse groups of defense-related genes, including those encoding pathogenesis-related (PR) proteins [42]. Moreover, the non-expresser of PR genes 1 (NPR1) and transcription factors such as TGACG-motif binding factors (TGAs) have been identified as key components of the SA response [43]. In our present study, the TGA and PR1 genes were found to be upregulated specifically and significantly at the SF stage ( Table 3), suggesting that SA may be involved in the regulation of wintersweet flower senescence. The reported involvement of SA in defense responses and in reactions to abiotic stresses in plants may also be relevant, since wintersweet blooms in winter when temperatures are low.
In addition, we found that a range of auxin-related, ABArelated and JA signaling pathway-related genes isolated in our study, such as Aux/IAA, ARF, GH3, SAUR,JAR1, JAZ, MYC2, SnRK2 and ABF, were differentially regulated during flower opening and senescence in wintersweet.
Taken together, these findings suggest that flower development in wintersweet is subject to complex regulation by multiple hormones. The regulatory mechanisms by which phytohormones are involved in wintersweet flower opening and senescence are still far from elucidated. The identification in this study of genes related to plant hormone signaling pathways is therefore an essential step in understanding fully the hormonal regulation of flowering and flower development in wintersweet and other winter-flowering species.

Floral Scent Biosynthesis and its Genetic Characterization
Floral scents are unique and distinct for each plant. A distinctive floral scent plays an important role in the reproductive processes of many plants and also enhances the aesthetic properties and economic value of ornamental plants and cut flowers. Flower scents vary between plant species on account of their characteristic profiles of volatile compounds of low molecular mass, such as monoterpenes, seisquiterpenes, benzenoids, phenylpropanoids, and fatty acid derivatives [44]. Many of these volatile compounds have already been identified in a number of flower species and in addition many of the biosynthetic enzymes involved have been identified and functionally characterized [45]. Some pathways have been characterized at the molecular-genetic level, primarily in model systems such as Clarkia breweri [46][47][48][49][50], Petunia [51], rose [52,53], and Antirrhinum [50].
In a previous study, more than 30 volatile compounds had been detected in wintersweet flowers; these consisted almost exclusively of volatile terpenoids (monoterpenes and sesquiterpenes) and benzenoids [5,6]. There was little molecular-genetic characterization of the biosynthesis of these compounds in wintersweet, however; and to the best of our knowledge, only the FPPS (farnesyl pyrophosphate synthase) gene and a homologue of the SAMT(Sadenosyl-L-methionine: salicylic acid carboxyl methyltransferase) gene had been cloned [54,55].
In this study, we isolated additional key floral scent-synthesizing genes from wintersweet, based on sequence annotations and analyses of changes in gene expression during flower development. The expression of the FPPS gene was found to be induced significantly during flower opening, consistent with previous findings, and was correlated with the emission of sesquiterpenoids [54]. Terpineol synthases (TERs) have been isolated from plant species representing six genera: Magnolia grandiflora [56], Pinus taeda [57], Santalum album [58], Vitis vinifera [59], Zea mays [60] and Nicotiana species [61]. The TER enzyme of M. grandiflora produces the single product, a-terpineol, whilst the TER enzymes of the other species produce additional compounds, such as limonene, bmyrcene, a-pinene, and b-pinene [61]. In wintersweet flowers, both a-terpineol and a-pinene have been detected [4,5].
The identification of genes of floral scent biosynthesis will facilitate future molecular-genetic and physiological studies of the regulation of floral scent production during flower development in wintersweet. It may also open new opportunities for the metabolic engineering of floral volatiles in plants [62].

Conclusions
This work presents the first de novo transcriptome sequencing analysis of wintersweet flower development using the Illumina RNA-Seq method. A total of 10,699 transcripts were assembled and 51,793 sequences were annotated. Our work provides an excellent platform for future genetic and functional genomics research in wintersweet and also provides an invaluable resource for genomics studies in other members of the Calycanthaceae family. Using Illumina sequencing-based gene analysis for the identification of DEGs, we identified a large number of candidate genes that were regulated either at the OF stage or at the SF stage, or at both, and we were able to identify a number of regulated pathways. In particular, we identified DEGs that were involved in plant hormone signal transduction pathways during flower development. The results indicated that flower opening and senescence may not be dependent upon the ethylene signal transduction pathway in wintersweet; in contrast, however, SA may be involved in the regulation of flower senescence. Finally, we obtained a number of key genes of floral scent biosynthesis and on this basis we have proposed a biosynthetic pathway for monoterpenes and sesquiterpenes in wintersweet flower based on annotated sequences. In summary, this comprehensive transcriptomic analysis has provided essential information on the genes and pathways which are involved in flower development in wintersweet. It will be fundamental to the development of wintersweet cultivars of improved horticultural or ornamental quality.

RNA Extraction, cDNA Library Construction and Illumina Deep Sequencing
For each flower developmental stage, three individual flowers from individual plants were randomly chosen and total RNA samples were extracted from each flower by using an RNAprep pure Plant RNA Purification Kit (Tiangen Biotech, Beijing, China). The quality and quantity of total RNA were analyzed using an UltrasecTM 2100 pro UV/Visible Spectrophotometer (Amersham Biosciences, Uppsala, Sweden) and by gel electrophoresis. For each developmental stage, the RNA samples from the three individuals were then pooled together in equal amounts to generate one mixed sample. These three mixed RNA samples (i.e. one for each flower stage) were subsequently used in cDNA library construction and Illumina deep sequencing.
Three cDNA libraries were prepared using a TruSeq TM RNA sample preparation Kit from Illumina (San Diego, CA, USA), and using in each case 5 mg of total RNA. Messenger RNA was isolated by polyA selection with oligo (dT) beads and fragmented using fragmentation buffer; cDNA synthesis, end repair, A-base addition and ligation of the Illumina-indexed adaptors were then all performed according to the Illumina protocol. Libraries were then size-selected on 2% Low Range Ultra Agarose for cDNA target fragments of 300-500 bp; this was followed by PCR amplification using Phusion DNA polymerase (NEB) for 15 PCR cycles. Following quantification by TBS380, paired-end libraries were sequenced using Illumina HiSeq TM 2000 system (26100 bp read length) at Shanghai Majorbio Bio-pharm Biotechnology Co., Ltd. (Shanghai, China). All raw-sequence reads data were deposited in NCBI Sequence Read Archive (SRA, http://www. ncbi.nlm.nih.gov/Traces/sra) with accession number SRA106143.

Transcriptome De novo Assembly
The raw reads were cleaned by removing low-quality reads (Q value ,25), adapter sequences, reads with ambiguous bases 'N', and fragments of less than 20 bp in length, using SeqPrep (https:// github.com/jstjohn/SeqPrep) and ConDeTri_v2.0.pl software(http://code.google.com/p/condetri/downloads/detail?name = con detri_v2.0.pl). De novo assembly of the Chimonanthus praecox transcriptome in the absence of a reference genome, using the three libraries obtained from the different flower stages, was accomplished using Trinity de novo transcriptome assembly software [21].

Sequence Annotation and Classification
For annotation, the sequences were searched against the NCBI non-redundant (NR) protein database [63] using BlastX, with a cut-off E-value of 10 25 . Gene ontology (GO) terms were extracted from the annotation of high-score BLAST matches in the NCBI NR proteins database (E value #1.0610 25 ) using blast2go (http://www.blast2go.com/b2ghome), and then sorted for the GO categories using in-house perl scripts [64]. Functional annotation of the proteome was carried out by a BlastP homology search against the NCBI Clusters of Orthologous Groups (COG) database (http://www.ncbi.nlm.nih.gov/COG/). KEGG pathway annotations were performed using Blastall software against the KEGG database.

Expression Analysis
After assembling the Chimonanthus praecox transcriptome, every RNA-seq library was separately aligned to the generated transcriptome assembly, using Bowtie [65]. The counting of alignments was performed using the RSEM package [66]. The DEGs were analyzed using the R Bioconductor package, edgeR [67]. The P-value set the threshold for the differential gene expression test. The threshold of the P-value in multiple tests was determined by the value for the false discovery rate (FDR) [68]. We used ''FDR#0.05 and the absolute value of Log 2 fold change (Log 2 FC)$2'' as the threshold to judge the significance of gene expression differences. For pathway enrichment analysis, DEGs were mapped to the terms in the KEGG database by using the KOBAS program (http://kobas.cbi.pku.edu.cn/home.do). Pathways were defined as significantly enriched pathways with respect to DEGs by taking a corrected P-value #0.05 as the threshold.

Real-time Quantitative PCR Validation of RNA-Seq Data
Ten DEGs involved in plant hormone transduction pathways were chosen for validation using Real-time qPCR. The primers, designed with the software Primer Primer 5.0, are listed in Table  S7. Total RNA was extracted from flowers at each stage using an RNAprep pure Plant RNA Purification Kit (Tiangen Biotech, Beijing, China) and reverse transcribed into cDNA using a PrimeScript RT Reagent Kit with gDNA Eraser (TaKaRa, Otsu, Japan). The real-time polymerase chain reaction (PCR) was performed on 10 mL reactions using 5 mL of Ssofast EvaGreen Supermix (Bio-Rad, Hercules, CA), 0.5 mL of each primer (final concentration of 500 nM), 3.5 mL of water, and 0.5 mL of cDNA template. PCR was carried out using the Bio-Rad CFX96 RealTime PCR system (Bio-Rad, US), using a denaturation temperature of 95uC for 30 s, followed by 40 cycles of 95uC for 5 s and 59uC for 5 s. To control the specificity of the annealing of the oligonucleotides, dissociation kinetics were performed using the real-time PCR system at the end of the experiment (60 to 95uC, continuous fluorescence measurement). A comparative Ct method (2 -DDCt ) of relative quantification [69] was used to analyze the realtime quantitative PCR data, using Bio-Rad CFX Manager Software, Version 1.6. The genes for actin and tubulin were used for the calculation of relative transcript abundance [14]. The sizes of the amplified products were confirmed through gel electrophoresis. Negative controls without templates were carried out concurrently. All quantitative PCRs for each gene used three biological replicates, with three technical replicates per experiment. Figure S1 The distribution of wintersweet genes in the plant hormone signal transduction pathway based on KEGG database. The genes involved in the pathway are colorcoded: green boxes, genes identified in our data; and light blue boxes, genes involved in the pathway present in KEGG database but undetectable in our data. (TIF)