Integrating Transcriptomic and GC-MS Metabolomic Analysis to Characterize Color and Aroma Formation during Tepal Development in Lycoris longituba

Lycoris longituba, belonging to the Amaryllidaceae family, is a perennial bulb bearing flowers with diverse colors and fragrance. Selection of cultivars with excellent colored and scented flowers has always been the breeding aim for ornamental plants. However, the molecular mechanisms underlying color fading and aroma production during flower expansion in L. longituba remain unclear. Therefore, to systematically investigate these important biological phenomena, the tepals of L. longituba from different developmental stages were used to screen and analyze the metabolic components and relevant genes. Utilizing the Illumina platform, a total of 144,922 unigenes were obtained from the RNA-Seq libraries. Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis indicated that the phenylpropanoid biosynthesis and flavonoid biosynthesis pathways might play important roles during color and aroma changes. Metabolomic analysis identified 29 volatile organic components (VOCs) from different developmental stages of L. longituba tepals, and orthogonal partial least-squares discriminate analysis (OPLS-DA) revealed that trans-β-ocimene—a terpene—was the most important aroma compound. Meanwhile, we found the content of anthocyanin was significantly reduced during the tepal color fading process. Then, we identified two dihydroflavonol-4-reductase (DFR) and three terpene synthase (TPS) genes, for which expression changes coincided with the production patterns of anthocyanins and trans-β-ocimene, respectively. Furthermore, a number of MYB and bHLH transcription factors (TFs) which might be involved in color- and aroma-formation were also identified in L. longituba tepal transcriptomes. Taken together, this is the first comprehensive report of the color and fragrance in tepals of L. longituba and these results could be helpful in understanding these characteristics and their regulation networks.


Introduction
Lycoris longituba, commonly known as Chinese tulip, is a bulbiferous species of the Amaryllidaceae family and distributed in central eastern China [1]. It can tolerate extremes of drought, waterlogging and shade, as well as poor soil conditions. Its plentiful flower colors, large tepals, elegant fragrance, and some medicinal potential make it a popular ornamental plant [2].

Anthocyanin Level in the Different Tepal Development Stages of Lycoris longituba
In the small bud stage, the red color of L. longituba 'Pink' tepal was very deep, and then the color intensity was significantly decreased with the rapid elongation of tepals, as shown in Figure 1a. In contrast, the tepal color of L. longituba 'White' was always white, as shown in Figure 1a. As shown in Figure 1b, the anthocyanin content in L. longituba 'Pink' was dramatically reduced from S1-P to S3-P, and nearly no anthocyanin was detected in S3-W, as shown in Figure 1b. These results suggested that content changes of anthocyanin could be the main reason that led to the tepals' red color fading of L. longituba 'Pink'.

Anthocyanin Level in the Different Tepal Development Stages of Lycoris longituba
In the small bud stage, the red color of L. longituba 'Pink' tepal was very deep, and then the color intensity was significantly decreased with the rapid elongation of tepals, as shown in Figure 1a. In contrast, the tepal color of L. longituba 'White' was always white, as shown in Figure 1a. As shown in Figure 1b, the anthocyanin content in L. longituba 'Pink' was dramatically reduced from S1-P to S3-P, and nearly no anthocyanin was detected in S3-W, as shown in Figure 1b. These results suggested that content changes of anthocyanin could be the main reason that led to the tepals' red color fading of L. longituba 'Pink'.

Transcriptome Sequencing and de Novo Assembly
Twelve total RNA samples were isolated from different L. longituba tepal developmental stages S1-P, S2-P, S3-P, and S3-W. These RNA samples were at concentrations of about 200-500 ng/μL with OD260/280 ≥ 1.9 and the RNA Integrity Numbers (RINs) of 8.6-10.0 were used for cDNA library construction. The Illumina HiSeqTM 4000 platform was used to obtain the dataset of 12 cDNA libraries. About 663.25 million raw sequencing reads with a length of 150 bp were generated, and after discarding the low-quality reads, we obtained about 85.23% (565.28 million) clean reads. For all 12 samples, the quality score above 20 (Q20) was ~98.20% and the GC percentages were 45.55-46.88%. Using Trinity software, the de novo assembly totally generated 144,922 unigenes, of average length 941 bp, from the twelve tepal transcriptomes, as shown in Table 1. In this research, the N50 was determined to be 1527 bp, which indicated that the quality of sequence assembly was good. All raw high throughput sequence data have been deposited in the NCBI Sequence Reads Archive (SRA) with the accession number PRJNA490415.

Transcriptome Sequencing and de Novo Assembly
Twelve total RNA samples were isolated from different L. longituba tepal developmental stages S1-P, S2-P, S3-P, and S3-W. These RNA samples were at concentrations of about 200-500 ng/µL with OD260/280 ≥ 1.9 and the RNA Integrity Numbers (RINs) of 8.6-10.0 were used for cDNA library construction. The Illumina HiSeqTM 4000 platform was used to obtain the dataset of 12 cDNA libraries. About 663.25 million raw sequencing reads with a length of 150 bp were generated, and after discarding the low-quality reads, we obtained about 85.23% (565.28 million) clean reads. For all 12 samples, the quality score above 20 (Q20) was~98.20% and the GC percentages were 45.55-46.88%. Using Trinity software, the de novo assembly totally generated 144,922 unigenes, of average length 941 bp, from the twelve tepal transcriptomes, as shown in Table 1. In this research, the N50 was determined to be 1527 bp, which indicated that the quality of sequence assembly was good. All raw high throughput sequence data have been deposited in the NCBI Sequence Reads Archive (SRA) with the accession number PRJNA490415.

Functional Classification of Genes during Tepal Development Stages
The assembled unigenes were annotated using blastx against NCBI nonredundant protein (Nr), Swiss-Prot, Kyoto Encyclopedia of Genes and Genomes (KEGG), and Cluster of Orthologous Groups (COG) protein sequence databases with an E-value ≤10 −5 . A total of 85,563 (59.04%) unigenes could be annotated while other unigenes had no significant BLAST hit, as shown in Table 2, indicating that numerous new genes specific to L. longituba are still functionally unknown and need to be further studied in the future. Remarkably, the plant species with the top three numbers of blastx hits were Elaeis guineensis (23,200 transcripts; 29.62%), Phoenix dactylifera (19,045 transcripts; 24.31%), and Musa acuminata (6512 transcripts; 8.31%), as shown in Table S1. These results implied that the assembled L. longituba transcripts shared similarity with transcripts from several monocotyledons and were reliable. Based on sequence homology, gene ontology (GO) assignment analysis was performed. Of the 85,563 annotated unigenes, 44,813 (52.37%) sequences were assigned into three main categories (biological process, cellular components, and molecular function), which could be further distributed under 58 GO terms, as shown in Figure S1. Metabolic process, cell, and catalytic activity were the most highly enriched GO terms in biological process, cellular components, and molecular function categories, respectively.

Identify Differentially Expressed Unigenes between Tepal Transcriptomes
To identify differentially expressed unigenes (DEGs) during tepal color fading, the unigenes that were differentially expressed between developmental stages and opening tepals with different colors were compared and shown in Figure 2. Among the four comparisons, the smallest number of DEGs was between the S1-P and S2-P libraries (4674), of which 2330 were up-regulated and 2344 were down-regulated, and the largest number of DEGs was between the S3-P and S3-W libraries (9958), with 5786 up-regulated and 5463 down-regulated unigenes. As the tepals faded from S2-P to S3-P, 8110 unigenes were differentially expressed, with 3526 up-regulated and 4584 down-regulated unigenes. In the S1-P vs. S3-P comparison, 11,024 DEGs were detected, including 5335 up-regulated and 5689 down-regulated unigenes. The comparison between S1-P and the other tepal developmental stages (S2-P and S3-P) showed that the number of up and down-regulated unigenes were both significantly increased as the tepals developed.

KEGG Pathway Enrichment
According to KEGG pathway enrichment analysis (p-value < 0.05), the phenylpropanoid biosynthesis (96 DEGs, ko00940), flavonoid biosynthesis (58 DEGs, ko00941), flavone and flavonol biosynthesis (23 DEGs, ko00944), and anthocyanin biosynthesis (4 DEG, ko00942) pathways which related to color formation were significantly different as compared to S3-P and S3-W. Interestingly, except for the anthocyanin biosynthesis and flavonol biosynthesis pathways, the phenylpropanoid biosynthesis and flavonoid biosynthesis pathways were also identified in the S1-P vs. S2-P, S2-P vs. S3-P, and S1-P vs. S3-P comparisons. These results suggested that the phenylpropanoid biosynthesis and flavonoid biosynthesis pathways could play critical roles during red color fading and aroma formation. The statistically enriched pathways between each two transcriptomes are shown in Figure  S2.

Validation of the Gene Expression Profiles by qRT-PCR
To validate the transcription profile revealed by RNA-Seq data, the expression levels of 21 genes from the flavonoid biosynthesis pathway were also assessed using qRT-PCR, as shown in Figure 3a. Linear regression analysis was used to obtain the overall correlation coefficient between RNA-Seq and qRT-PCR data, which showed a good correlation (R = 0.89) between these two data, as shown in Figure 3b, indicating that the 12 transcriptomics data were reliable.

KEGG Pathway Enrichment
According to KEGG pathway enrichment analysis (p-value < 0.05), the phenylpropanoid biosynthesis (96 DEGs, ko00940), flavonoid biosynthesis (58 DEGs, ko00941), flavone and flavonol biosynthesis (23 DEGs, ko00944), and anthocyanin biosynthesis (4 DEG, ko00942) pathways which related to color formation were significantly different as compared to S3-P and S3-W. Interestingly, except for the anthocyanin biosynthesis and flavonol biosynthesis pathways, the phenylpropanoid biosynthesis and flavonoid biosynthesis pathways were also identified in the S1-P vs. S2-P, S2-P vs. S3-P, and S1-P vs. S3-P comparisons. These results suggested that the phenylpropanoid biosynthesis and flavonoid biosynthesis pathways could play critical roles during red color fading and aroma formation. The statistically enriched pathways between each two transcriptomes are shown in Figure S2.

Validation of the Gene Expression Profiles by qRT-PCR
To validate the transcription profile revealed by RNA-Seq data, the expression levels of 21 genes from the flavonoid biosynthesis pathway were also assessed using qRT-PCR, as shown in Figure 3a. Linear regression analysis was used to obtain the overall correlation coefficient between RNA-Seq and qRT-PCR data, which showed a good correlation (R = 0.89) between these two data, as shown in Figure 3b, indicating that the 12 transcriptomics data were reliable.

Metabolome Analysis of Lycoris longituba Tepal Development Stages by GC-MS
To investigate the volatile metabolic components changes of L. longituba tepals in the opening processes, we obtained the GC-MS total ion current (TIC) chromatograms for nine L. longituba tepal samples from three typical developmental stages (S1-P, S2-P, and S3-P), as shown in Figure 1a. The obvious differences of chromatographic peaks were observed between sample groups, and the retention times were fairly consistent and reproducible, as shown in Figure 4a. In this study, a total of 29 metabolites were identified in our sample libraries across all samples, as shown in Table 3.
To assess the volatile components profile changes during L. longituba tepal development, the orthogonal partial least-squares discriminate analysis (OPLS-DA) plot was generated from the GC-MS metabolite data of S1-P, S2-P, and S3-P tepals and showed clear metabolic differences between two stages. Remarkably, S1-P, S2-P, and S3-P tepals could be completely separated sufficiently by use of two principal components. The first principle component (PC1, accounting for 49.89%) and the second component (PC2, accounting for 32.96%) of the variation in the data could separate all three types of tepals with no outliers, as shown in Figure 4b. The contribution of each variable to PC1 and PC2 was also calculated by giving each variable a weight value. The top two core differential metabolites of PC1 and PC2 discrimination were caryophyllene and trans-β-ocimene, as shown in Figure 4c. Interestingly, the trans-β-ocimene, which had a high content in S3-P, was not detected in S1-P and S2-P, as shown in Table 3.

Metabolome Analysis of Lycoris longituba Tepal Development Stages by GC-MS
To investigate the volatile metabolic components changes of L. longituba tepals in the opening processes, we obtained the GC-MS total ion current (TIC) chromatograms for nine L. longituba tepal samples from three typical developmental stages (S1-P, S2-P, and S3-P), as shown in Figure 1a. The obvious differences of chromatographic peaks were observed between sample groups, and the retention times were fairly consistent and reproducible, as shown in Figure 4a. In this study, a total of 29 metabolites were identified in our sample libraries across all samples, as shown in Table 3.

Analysis of Candidate Genes Related to Color and Fragrance Metabolics
To explore the genetic regulation of L. longituba tepal color fading and aroma emission, the genes which have been reported to be involved in these two metabolic pathways were selected. With the development of tepals, several anthocyanin biosynthesis structural genes had the lowest expression levels in S3-P, such as the DFRs, CHIs, CHS2, F3'H1, and FLSs, as shown in Figure 3a. Especially, the DFR-annotated unigenes (DFR2-1 and DFR2-2) which had more than a 7-fold down-regulated expression level in S1-P vs. S3-P, as shown in Figure 3a. Five TPS genes were also identified from the DEGs, and interestingly three (Unigene81776, CL6106.Contig2, and Unigene3859) of them were predominantly expressed in S3-P, as shown in Figure 5c.
The spatial and temporal expression of pigment structural and aroma genes were usually controlled by transcription factors from MYB and bHLH [14,25]. In this study, 35 MYBs and 29 bHLHs with fragments per kilobase per million fragments (FPKM) ≥ 5 were identified from DEGs, as shown in Figure 5a,b. Among them, six MYBs and four bHLHs (in black frames) had the similar down-regulated expression trends with DFR2-1 and DFR2-2, and four MYBs and one bHLH (in red frames) had the similar up-regulated expression patterns with the above three TPSs, as shown in Figure 5. To assess the volatile components profile changes during L. longituba tepal development, the orthogonal partial least-squares discriminate analysis (OPLS-DA) plot was generated from the GC-MS metabolite data of S1-P, S2-P, and S3-P tepals and showed clear metabolic differences between two stages. Remarkably, S1-P, S2-P, and S3-P tepals could be completely separated sufficiently by use of two principal components. The first principle component (PC1, accounting for 49.89%) and the second component (PC2, accounting for 32.96%) of the variation in the data could separate all three types of tepals with no outliers, as shown in Figure 4b. The contribution of each variable to PC1 and PC2 was also calculated by giving each variable a weight value. The top two core differential metabolites of PC1 and PC2 discrimination were caryophyllene and trans-β-ocimene, as shown in Figure 4c. Interestingly, the trans-β-ocimene, which had a high content in S3-P, was not detected in S1-P and S2-P, as shown in Table 3.

No.
Name S1-P S2-P S3-P  down-regulated expression trends with DFR2-1 and DFR2-2, and four MYBs and one bHLH (in red frames) had the similar up-regulated expression patterns with the above three TPSs, as shown in Figure 5.

Discussion
Flower color and fragrance are considered to be two critical characters that influence plants' ornamental value and insect pollination. The pigmentation and aroma formation processes involve multiple gene expression networks and complex biochemical pathways. However, the underlying molecular regulation mechanism of flower color fading and fragrance synthesis in L. longituba remains to be uncovered. To date, transcriptome sequencing technology has been widely used in exploring the critical metabolic pathways and functional genes that are involved in color biosynthesis

Discussion
Flower color and fragrance are considered to be two critical characters that influence plants' ornamental value and insect pollination. The pigmentation and aroma formation processes involve multiple gene expression networks and complex biochemical pathways. However, the underlying molecular regulation mechanism of flower color fading and fragrance synthesis in L. longituba remains to be uncovered. To date, transcriptome sequencing technology has been widely used in exploring the critical metabolic pathways and functional genes that are involved in color biosynthesis and aroma formation in various species, even those that lack reference genomes. In this study, the tepal transcriptomes in small bud, medium bud, and opening stages of L. longituba were compared. The phenylpropanoid biosynthesis, flavonoid biosynthesis, and flavone and flavonol biosynthesis pathways which related to the color and aroma formation were highlighted between each two transcriptomes by KEGG pathway enrichment analysis, as shown in Figure S2. This result implied that these pathways could play critical roles in the two ornamental characters formation in L. longituba tepals.
It has been demonstrated that the different content of anthocyanidin was a critical reason for red color alteration in plant organs [9,35]. In this research, tepal color was altered from deep red to pink accompanying L. longituba flower bud development, as shown in Figure 1a. Remarkably, the anthocyanidin level of the tepal was significantly reduced during the color fading processes, as shown in Figure 1b. This result suggested that the decreasing of anthocyanidin accumulation could be the directly responsible for the color fading of L. longituba.
Dihydroflavonol-4-reductase (DFR), which is a downstream enzyme in anthocyanin biosynthesis, could regulate the critical rate-limiting step of anthocyanin biosynthesis processes by catalyzing dihydroflavonols to leucoanthocyanidins [4,36]. Recently, a series of DFR gene homologs have been identified in various plant species such as in Rosa rugosa [4], Camellia sinensis [37], Populus trichocarpa [38], and Calibrachoa hybrida [39]. The previous works have confirmed that overexpression of DFR genes within tobacco and petunia could accelerate the anthocyanin accumulation and promote the red coloration of flowers [4,39]. Here, two LlDFRs (DFR2-1 and DFR2-2) were isolated from DEGs. Both of them had the highest expression levels in the S1-P stage (deep red bud), then significantly reduced in the S2-P stage (pink red bud), which were positively correlated with the content variation trend of anthocyanins in L. longituba tepals, as shown in Figures 1b and 3a. Previous studies have demonstrated that the expressions of anthocyanin biosynthesis structural genes are controlled by the MYB, bHLH, and WD40 transcription factors [40], such as PhAN11 in Petunia hybrida [41], PeMYB11 in Phalaenopsis spp. [10], MdMYB10 and MdbHLH3 in Malus domestica [11,42], CmMYB6 in Chrysanthemum [43], AtMYBL2 in Arabidopsis thaliana [44], FaMYB1 in Fragaria ananassa [45], as well as PpMYB16 and PpMYB111 in Prunus persica [46]. In this study, several MYB and bHLH genes, which had similar expression patterns with anthocyanidins and DFR genes, were identified in these transcriptomes, as shown in Figure 5. This evidence indicated that the down-regulated expression of DFR genes during the tepal developmental processes might be controlled by the above transcription factors, and these TFs could directly control the flower color fading.
Recently, GC-MS technology has been successfully used to analyze qualitative and quantitative differences in aroma-related metabolomics, and some important volatile terpenes which play dominant roles in the floral scent formation have been well demonstrated, such as in Lilium [47], Chimonanthus praecox [18], and Osmanthus fragans [19]. Here, a total of 29 floral volatile organic compounds (VOCs) and 6 volatile terpenes were identified from three developmental stages of L. longituba tepals, as shown in Figure 4. In previous studies, linalool has been verified as one of the critical aroma metabolites in various plants [18,19]. Conversely, linalool was undetected in the tepals of L. longituba, as shown in Table 3, which indicated that the composition of L. longituba floral aroma could be very particular. Interestingly, some VOCs were specifically emitted in the different stages; especially the trans-β-ocimene and β-myrcene which belong to the monoterpenes, which were only detected in the strongest aroma period of opening flowers, as shown in Table 3. Meanwhile, OPLS-DA loading values showed that the top two critical components were caryophyllene and trans-β-ocimene, as shown in Figure 4c, and the content of trans-β-ocimene was much higher than caryophyllene in the full-blooming floral tepals (S3-P), as shown in Table 3. Taken together, we concluded that the trans-β-ocimene could be the most important aroma compound during the L. longituba floral fragrance emission.
It has been demonstrated that the TPSs could control the β-ocimene biosynthesis at the transcriptional level [48][49][50]. In this study, the expression patterns of three TPSs were in excellent agreement with the emission trend of trans-β-ocimene in L. longituba tepals, as shown in Figure 5c and Table 3, implying that the three TPS members might directly influence the protein abundance of ocimene synthase. Some members of the MYB and bHLH families have also been confirmed that could control the floral aroma formation, such as the CpMYC2 [18], AtMYC2 [24], HcMYB1, and HcMYB2 [25]. In this research, several MYB and bHLH genes which had the consistent expression patterns with trans-β-ocimene and TPS genes were identified, as shown in Figure 5a,b, suggesting that these genes could be involved in the regulation of aroma formation.
The color fading and aroma formation during the L. longituba tepal development could be two very complex and dynamic processes. While many candidate genes have been identified in this work, the functions of these members still need to be investigated in our future studies.

Plant Materials
Plants of L. longituba 'Pink' and 'White' were grown in the Lycoris Experimental Plantation of Nanjing Forestry University in Nanjing, China. The L. longituba tepals of three typical developmental stages were selected based on tepal length, which were S1-P (15 ± 5 mm tepals at small bud stage), S2-P (55 ± 5 mm tepals at medium bud stage), S3-P (90 ± 5 mm tepals at opening stage) of 'Pink', and S3-W (90 ± 5 mm tepals at opening stage) of 'White', as shown in Figure 1a. At each stage, the fresh samples of tepals were harvested and immediately frozen in liquid nitrogen and stored at -80 ºC.

Anthocyanin Level Measurement
Freeze-dried tepals were finely ground and 0.2 g was extracted with 2 mL acidic methanol (0.1% hydrochloric acid) at 4 ºC in darkness for 12 h, mixing the extract up and down every 6 h. Then, the extract was centrifuged at 12,000 rpm for 10 min, the supernatant was diluted 4 times with acidic methanol, and the absorbances were tested spectrophotometrically at 530 and 657 nm. Finally, total anthocyanin was defined using the equation: Q = (A 530 − 0.25 × A 657 ) × FW − 1, where Q = total anthocyanins; A 530 = absorption at 530 nm; A 657 = absorption at 657 nm; FW = fresh weight of tepals (g). Three biological replicates were performed for each group.

RNA-Seq and de Novo Assembly
Total RNA from three biological replicates was independently isolated from the tepals of L. longituba using RNAiso Reagent (Takara, Otsu, Japan) according to the previously described method [51]. The RIN value of total RNA sample was examined with Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, USA), and the concentration was assessed using NanoDrop (Thermo Scientific, Waltham, MA, USA).
Twelve cDNA libraries which consisted of separate RNA samples from L. longituba 'Pink' tepals of three different developmental stages (S1-P, S2-P, and S3-P) and L. longituba 'White' tepals (S3-W) were prepared using the TruSeq RNA Sample Preparation Kit (Illumina, San Diego, CA, USA) according to the manufacturer's instructions. Firstly, the enriched mRNA were fragmented into short fragments and reverse transcripted into first-strand cDNA. Then after second-strand cDNA synthesis, end repair, adapter ligation, and PCR amplification, the cDNA library products were sequenced on an Illumina HiSeq TM 4000 (Illumina, San Diego, CA, USA) instrument using paired-end sequencing technology by staff at Gene Denovo Biotechnology Corporation (Guangzhou, China). After the raw reads removing the adapter sequences, reads with more than 20% low-quality bases (quality value <20) and ambiguous nucleotides (denoted with an "N" in the sequence trace) of the raw reads, the high-quality clean reads were accomplished assembled using Trinity software [52]. In order to avoid the interference of alternative splicing transcripts, only the longest transcript was taken as the unigene in this research.

Sequence Annotation and Gene Expression Difference Analysis
The assembled unigenes functional annotations were performed through a blastx search against four public protein databases including NCBI nonredundant protein (Nr), Swiss-Prot, Kyoto Encyclopedia of Genes and Genomes (KEGG), and Cluster of Orthologous Groups (COG) of proteins, with an E-value of less than 1e-5. The Blast2GO software (http://blast2go.bioinfo.cipf.es/) was used to acquire GO terms of the unigenes.
The expression level of unigenes was calculated using fragments per kilobase per million fragments mapped (FPKM) method [53]. In this study, the false discovery rate (FDR) ≤ 0.001 and the absolute value of log 2 Ratio ≥ 2 were taken as the threshold for significantly differential expression of unigenes. All differentially expressed unigenes (DEGs) were mapped to the KEGG pathway database and the numbers of unigenes for every KEGG Orthology (KO) term were calculated. Significantly enriched KO terms from the set of DEGs were identified by comparing the observed DEG count to the expected count of the genes involved in a given pathway with a random distribution of the L. longituba tepal transcriptome using the formula of the hypergeometric test [54].

qRT-PCR Analysis
Using TransScript One-Step gDNA Removal and cDNA Synthesis SuperMix kit (Transgene, Beijing, China), 2 µg of total RNA was reverse transcripted to the first-strand cDNA on the basis of the manufacturer's instructions. The SYBR Premix Ex TaqTM II kit (Takara) was used to perform qRT-PCR in the ABI 7500 Fast Real-Time PCR System (Applied Biosystems, Applied Biosystems, Cheshire, UK according to the manual's description [55]. The relative expression level of genes was calculated by the 2 −∆∆CT method. To ensure reliability, the expression data of each gene was obtained from three independent biological replications (each biological replication included three technical replications). The data are shown as mean values ± SE (standard error). All primer pairs were designed by Primer 5 software and listed in Table S2, and the specificity of them was assessed by the sequencing of amplified qRT-PCR products. The ELF gene was used as an internal reference control [56].

GC-MS Analysis
Fresh tepals at three different stages (S1-P, S2-P, and S3-P), defined by the size of the flower as shown in Figure 1a, were picked from the plants at the same time as samples that were collected for the above transcriptome studies. Sampling was replicated four times, and the samples were quickly put into polyethylene bags impermeable to gases, kept in the ice-box, and analyzed immediately. Headspace solid phase microextraction (SPME) combined with GC-MS was used to determine the identity and quantity of the fragrance volatiles. Tepals (0.3 g) were placed in a 4 mL solid-phase microextraction vial (Supelco Inc, Bellefonte, PA, USA), 1 µL of 1000x diluted ethyl caprate (Macklin Inc, Shanghai, China) was added, and vials were capped with a 65 µm DB-5 ms extraction head (Supelco Inc). The oven temperature was programmed at 60 • C for 2 min, increasing at 5 • C/min to 150 • C, then increasing at 10 • C/min to reach 250 • C, followed by maintaining the temperature of the transfer line at 250 • C and helium was used as the carrier gas at a linear velocity of 1.0 mL/min. Mass detector conditions on MS were carried out according to our previous method [57]. The quantities of the volatile aroma compounds were calculated by normalizing the peak-areas and volatile compounds were first identified using the NIST98 database (Agilent). SIMCA-P 11.5 software (Umetrics AB, Umea, Sweden) was selected to test the differences in the metabolite levels of L. longituba tepals from different developmental stages by orthogonal partial least-squares discriminate analysis (OPLS-DA).

Supplementary Materials:
The following are available online http://www.mdpi.com/2223-7747/8/3/53/s1. Figure S1: GO classification of DEGs among different samples in L. longituba, Figure S2: The KEGG enrichment of DEGs among the different samples, Table S1: Species distribution of the BLAST hits for all homologous sequences, Table S2: Primers used for qRT-PCR.

Conflicts of Interest:
The authors declare no conflict of interest.