In this study, two Lilium genotypes, Lilium ‘Siberia’ and Lilium ‘Novano’ were used as plant materials. The bulbs were cultured in plastic pots (20 cm diameter, 20 cm height) containing medium comprised of peat and vermiculite at a ratio of 2 to 1 in the greenhouse at the Science Park of Beijing University of Agriculture under a 16/8-h light/dark 25/20°C cycle (Hu et al., 2016). The plants were irrigated every 3 days and were supplied with full-strength Hoagland's nutrient solution every 2 weeks (Hu et al., 2016). The flowering period was divided into 4 stages, i.e., the bud stage (BS), half-bloom stage (HS), full-bloom stage (FS), and late-bloom stage (LS). The floral scents were collected and analyzed during these 4 stages. Petals were also collected during these 4 stages and were immediately frozen in liquid nitrogen prior to further analysis. The FS petals were used to detect the transcriptome. Three biological replicates were collected per stage.

Using dynamic headspace sampling the floral scent was collect (Hu et al., 2016). An individual flower was placed in a Reynolds oven bag (16 × 17.5 in), which releases and absorbs few volatiles. A stainless steel tube (0.25 × 3.5 in) containing Tenax-GR (60–80 mesh, Chrompack) was used as the volatile trap, and care was taken to ensure the flower was not touched (Hu et al., 2016). A portable air sampler (QC-1; Beijing Municipal Institute of Labor Protection, China) served as the pump, and air filtered by charcoal was pumped into the bag. The volatiles were collected for 20 min at a flow rate of 300 mL·min⁻¹ (Hu et al., 2015, 2016). Then, the stainless steel tubes were sealed and placed in a refrigerator.

The automated thermal desorption—gas chromatography/mass spectrometry (ATD-GC/MS) technique was used to analyze the floral scent. The floral scent collected in the stainless steel tube was desorbed by heating in an ATD (Auto Thermal Desorber, TurboMatrix 650, PerkinElmer) at 260°C for 10 min and was then cryofocused in a cold trap where the temperature was maintained at −25°C for 3 min. Subsequently, the cold trap was rapidly heated to 300°C, which was maintained for 5 min, to transport the volatiles to a GC (Clarus 600, Perkin Elmer). The GC was equipped with a capillary DB-5MS column (30 m × 0.25 mm i.d. with a 0.25-μm film thickness). Helium was used as the carrier gas. The GC was programmed at 40°C for 2 min, 4°C·min⁻¹ up to 160°C, followed by 20°C·min⁻¹ up to 270°C, and held at 270°C for 3 min. The MS (Clarus 600T, Perkin Elmer) was operated in EI ionization mode at 70 eV, and a mass scan range of 29–600 amu was monitored. The interface and ion source temperatures were 250°C and 220°C, respectively.

According to the retention indices, the compounds were identified by searching the NIST08 and WILEY library in the TurboMass Ver5.4.2 software. α-pinene (Fluka, USA) was used as an external standard to determined the release amounts of volatile components by dissolution in ethyl acetate with different solution concentrations (Hu et al., 2009, 2015, 2016), and μg·h⁻¹ was used as the unit.

Total RNA was extracted from the lily petals at FS using RNAiso Plus (TransGen Biotech, Beijing, China). The quality and quantity of purified RNA were examined using a NanoDrop ND-1000 UV/Visible spectrophotometer (Wilmington, DE, USA). The RIN (RNA integrity number) values (>8.0) of these samples were assessed using an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA) for gel electrophoresis. High-quality RNA was used in cDNA library construction and Illumina deep sequencing.

The construction of the libraries and the RNA-Seq were performed by the Biomarker Biotechnology Corporation (Beijing, China). mRNA was enriched and purified with oligo(dT)-rich magnetic beads and then broken into short fragments. Using these cleaved mRNA fragments as templates, first- and second-strand cDNA were synthesized. The resulting cDNAs were then subjected to end-repair and phosphorylation using T4 DNA polymerase and Klenow DNA polymerase. Subsequently, an ‘A’ base was inserted as an overhang at the 3′ ends of the repaired cDNA fragments, and Illumina paired-end solexa adaptors were subsequently ligated to these cDNA fragments to distinguish the different sequencing samples. The products of the ligation reaction were purified and selected on 2% agarose gel for downstream enrichment. Then the purified cDNA template was enriched by PCR amplification. Finally, 6 libraries were sequenced using an Illumina HiSeq™ 2500 (SRA accession number: SRP112400). Data analysis and base calling were performed using the Illumina instrument software.

The raw reads were first filtered to identify clean reads by removing the reads with only adaptor sequences and unknown nucleotides >5%, as well as low quality reads. The clean reads were then assembled de novo using the Trinity platform based on the parameters “K-mer = 25, group pairs distance = 300” (Li et al., 2016). The short reads were first assembled into longer contigs on the basis of their overlapping regions. Then, different contigs were further recognized by mapping clean reads back to the corresponding contigs based on their paired-end information, and then the gene sequences were obtained. Finally, the potential gene sequences were clustered to gain unigenes using the TGI tool (Pertea et al., 2003). In order to annotate, the unigene sequences were searched against a series of protein databases using BLASTx, with a cut-off E-value of 10⁻⁵, including the NCBI Cluster of Orthologous Groups of proteins (COG), the Gene Ontology (GO), the Kyoto Encyclopedia of Genes and Genomes (KEGG), the NCBI eukaryotic Orthologous Groups (KOG), the Swiss-Prot, NCBI non-redundant (Nr), and the Pfam. The deduced amino acid sequences of unitranscripts were demanded to be longer than 70% of the corresponding sequences. If a unigene met the criteria, it was assumed to contain a near full-length contig. Alternatively, targeted assembly was performed to obtain even greater coverage of the respective genes. All reads in the databases examined were mapped to the reference sequences, and the mapped reads were then assembled using clustering and CAP3.

To evaluate the depth of coverage, all usable reads were realigned to each unigene using SOAPaligner and then normalized into RPKM values (reads per kb per million reads). The unigenes showed differential expression levels among different samples were calculated based on the ratio of the RPKM values. In order to compute the significance of the difference in transcript abundance, the threshold of the P-value in multiple tests was idenfied using the false discovery rate (FDR) control method. Only the unitranscripts with an absolute value of log2 ratio ≥2 and an FDR significance score <0.001 were used for subsequent analysis. If FDR was lower than 0.05 and the highest RPKM (reads per kilobase per million reads) of the unigene was twice that of the lowest one, the unigene was considered as differentially expressed gene (Zheng et al., 2015).

The unigenes were subjected to quantitative real-time PCR (qRT-PCR) with specific primers (Supplementary Table 1). SYBR Green PCR Master Mix was used to detect the PCR products on a 7900 HT Sequence Detection System (Applied Biosystems). QRT-PCR was performed using the SYBR Premix Ex Taq Kit (TaKaRa) according to the manufacturer's protocol. The results were normalized to the expression level of the constitutive actin gene, and a comparative Ct method (2^−ΔΔct) of relative quantification was used to evaluate the quantitative variation (Zheng et al., 2015). There were three independent biological replicates and three technical replicates of each biological replicate for each sample. The samples in qRT-PCR validation were the same as those for Illumina sequencing.

Figure 1 showed the total amounts of floral scent emitted from Lilium ‘Siberia’ (Hu et al., 2016) and Lilium ‘Novano’ at different flowering stages. These values initially increased and then decreased with flower development, showing similar patterns. The maximum amounts released from Lilium ‘Siberia’ and Lilium ‘Novano’ occurred at FS, which were 7.85- and 6.90-fold higher, respectively, than those at BS. However, the amount released from Lilium ‘Siberia’ was significantly higher than that from Lilium ‘Novano’ (P < 0.05). At HS and FS in particular, the amounts released from Lilium ‘Siberia’ were ~11.19- and 8.32-fold higher, respectively, compared with those from Lilium ‘Novano’.

We have detected 7 main categories of volatile components in the floral scents of Lilium ‘Siberia’ (Hu et al., 2016), including terpenoids, benzenoids and derivatives, alkanes, alcohols, aldehydes, ketones, and esters, which were also found in the floral scent of Lilium ‘Novano’ (Figure 2). The release amounts of these components also increased initially and then decreased during the flowering stage. The maximum release amount occurred at FS, showing a pattern similar to the total release amount. Among the volatile components of Lilium ‘Siberia’, the terpenoid compound showed the highest release amount (Hu et al., 2016). The release amounts of terpenoid compounds at different flowering stages accounted for over 72% of the total release amount from Lilium ‘Siberia’ (Table 1), but in Lilium ‘Novano’, the relative release amounts of terpenoid compounds did not exceed 16% (Table 1).

Except for alkane and ketone, the amounts of the 5 other categories of components released from Lilium ‘Siberia’ (Hu et al., 2016) were higher than those released from Lilium ‘Novano’ at different flowering stages. At BS, HS, FS, and LS, the amounts of terpenoid compounds released from Lilium ‘Siberia’ were ~7.44-, 11.19-, 8.33-, and 2.94-fold higher, respectively, compared with Lilium ‘Novano’. Compared to Lilium ‘Novano’ (Figure 2), Lilium ‘Siberia’ emitted higher amounts of aldehydes (Hu et al., 2016) at BS and HS. The amounts of esters released from Lilium ‘Siberia’ (Hu et al., 2016) were also higher than those from Lilium ‘Novano’ at BS, HS, and FS, i.e., 2.60-, 2.19-, and 1.75-fold, respectively (Figure 2). In the floral scent of Lilium ‘Novano’, benzenoid and derivative and ester accounted for the high relative release amounts. So the analysis on release amounts of different components indicated that the significant differences in the amounts of the terpenoid compounds between Lilium ‘Siberia’ and Lilium ‘Novano’ mainly contributed to the difference in floral scent.

The terpenoid components were also analyzed, and monoterpenes were found to be the dominant components in Lilium ‘Siberia’ (Hu et al., 2016). A total of 7 monoterpenes were detected in the floral scent of Lilium ‘Siberia’, including α-pinene, myrcene, ocimene, linalool, limonene, E,E-2,6-dimethyl-1,3,5,7-octatetraene, and 2,6-dimethyl-3,7-octadiene-2,6-diol (Hu et al., 2016). We found that the release amounts of the 7 monoterpenes all peaked at FS (Hu et al., 2016). In the floral scent of Lilium ‘Novano’ 4 monoterpenes, including myrcene, ocimene, linalool, and E,E-2,6-dimethyl-1,3,5,7-octatetraene, were found (Figure 3). Among these 4 compounds, E,E-2,6-dimethyl-1,3,5,7-octatetraene and ocimene were released in the highest amounts, but the maximum values did not exceed 6.00 μg·h⁻¹ at FS. The release amount of myrcene was the lowest among the 4 monoterpenes, which only peaked at 0.68 ± 0.22 μg·h⁻¹.

To understand the molecular basis of the difference in the floral scent between Lilium ‘Siberia’ and Lilium ‘Novano’, the flowers at different developmental stages were used to build 6 libraries for high-throughput sequencing. We obtained a total of 29.24 Gb of raw data for the 6 samples (Table 2). We discarded low-quality reads, which contained adapters and unknown or low-quality bases, and after stringent quality checks and data cleaning, a total of 116.07 Mb of clean reads were obtained (Table 2). The GC (guanine + cytosine) contents of these samples were 47.26–50.23%, with an average of 49.27%. The average Q20 and Q3 percentage reached 95.46% and 90.59%, respectively (Table 2).

Based on the high quality reads, 2,864,438 contigs were assembled with an N50 length of 121 bp and an average length of 81 bp, including 94.67% contigs smaller than 200 bp. Using paired-end joining and gap-filling, the contigs were further assembled into 229,128 scaffolds with an N50 length of 1,472 bp and an average length of 876 bp, including 20,354 scaffolds larger than 2,000 bp. The de novo assembly yielded 124,233 unigenes with an N50 length of 986 bp and an average length of 615 bp (Table 3). Of these unigenes, 84.26% (104,682) were shorter than 1,000 bp, 10.48% (13,201) ranged from 1,000 to 2,000 bp, and the remaining 5.26% (6,530) were longer than 2,000 bp (Table 3). Among the unigenes longer than 2,000 bp, there were 2,212 unigenes whose lengths exceeded 3,000 bp, accounting for 1.78% of all unigenes.

All unigenes were aligned to 7 protein databases including COG, GO, KEGG, KOG, Swiss-Prot, and Nr using BLASTx with an E-value threshold of 10⁻⁵ and Pfam using HMMER with an E-value threshold of 10⁻¹⁰. As shown in Table 4, of 35,749 unigenes annotated, 34,717 (97.11%) unigenes presented significant BLASTx matches in the Nr database. There were 14,642 unigenes whose lengths were longer than 1000 bp, accounting for ~42.18% of 34,717 unigenes. Based on comparison against the Swiss-Prot database, 23,338 (65.28%) unigenes had significant matches, and ~49.65% (11,588) of unigenes were longer than 1,000 bp. In the Pfam and GO databases, 18,726 (52.38%) and 17,644 (49.36%) unigenes were also found to have significant matches respectively, and 6,210 (17.37%) unigenes were similar to proteins in the KEGG database.

To further evaluate the completeness of our transcriptome library and the effectiveness of our annotation process, we searched the annotated sequences for genes with COG (cluster of orthologous groups) classifications, and 9,183 unigenes were assigned to the COG classification (Figure 4). Among the 25 COG categories, the cluster for “General function prediction only” (2,295, 24.99%) represented the largest group, followed by “Replication, recombination and repair” (1,359, 14.80%), “Transcription” (1,155, 12.58%), and “Signal transduction mechanisms” (910, 7.19%). The “Secondary metabolites biosynthesis, transport and catabolism” category about which we were concerned accounted for 4.38% (402) of the sequences.

Gene ontology (GO) was also used to classify the functions of the predicted unigenes. Based on the sequence homology, 17,644 sequences were categorized into 52 functional groups (Figure 5). The assigned functions of the unigenes covered a broad range of GO categories. The unigenes were assigned to three main categories including cellular component, molecular function, and biological process categories (Figure 5). Of these, the biological process category constituted the majority followed by the cellular component and molecular function categories. In the biological process category, the metabolic process (10,841, 61.44%) and cellular process (9,160, 51.92%) indicated that some important metabolic activities occurred in Lilium petals, including monoterpene biosynthesis. A total of 2,465 (13.97%) and 1456 (8.25%) unigenes were assigned to biological regulation and the developmental process, respectively. In the cellular component category, cell part (8,417, 47.70%), cell (8,343, 47.29%), and organelle (6,900, 39.11%) were prominently represented. In the molecular function category, binding (2,904, 16.46%) and catalytic activity (2,340, 13.26%) represented the majority. In addition, 1,057 (5.99%) unigenes were involved in transporter activity.

Based on the comparison against the KEGG database that records networks of molecular interactions and reactions, 6,210 unigenes were assigned to 118 KEGG pathways (Supplementary Table 2). These annotations provide a valuable resource for investigating the processes, functions, and pathways involved in floral scent biosynthesis.

In the comparison of floral scents between Lilium ‘Siberia’ and Lilium ‘Novano’, we found a significant difference in the amounts of monoterpenes released (Figure 3), which might be the key reason resulting in the fragrance difference. Thus, we were interested in the gene expression of the monoterpene biosynthesis. Genes involved in the pathways of terpenoid backbone biosynthesis and monoterpenoid biosynthesis are related to the production of monoterpenoid compounds in plants.

Through mapping to the KEGG reference pathways, a total of 39 annotated unigenes were assigned to the pathway of terpenoid backbone biosynthesis (Supplementary Table 3). Though no unigenes were mapped to the monoterpene biosynthesis pathway, based on GO function classification combined with the annotation of the Swissprot and nr databases, we found that 7 unigenes were annotated to monoterpene synthases mediating the biosynthesis of monoterpenes (Supplementary Table 4).

Differential gene expression of the two Lilium plants was analyzed at FS. We performed a statistical analysis on the genes with an RPKM value ≥2 to reduce false positives and false negatives. We filtered the data using an FDR ≤ 0.001 and |log2 (ratio)| ≥ 2. There was a significant difference in gene expression between the two Lilium plants in the heat-map (Supplementary Figure 1). The expression of 6,496 DEGs was found to be significantly changed between Lilium ‘Siberia’ and Lilium ‘Novano’, including 2,702 up-regulated unigenes and 3,794 down-regulated unigenes (NFS vs. SFS), among which 4,739 unigenes were annotated (Figure 6; Supplementary Table 5).

The enrichment factor of DEGs in the terpenoid backbone biosynthesis pathway was lower than in other pathways except sulfur metabolism, which indicated significant enrichment of DEGs in this pathway (Supplementary Figure 2). Among the DEGs annotated by the KEGG database, 16 unigenes were assigned to the terpenoid backbone biosynthesis (Table 5), and the sequence data were shown in Supplementary Table 6. Except DXS-2, these DEGs, including DXS-1, DXR, HDS, HDR, IDI1, GPS, and GGPS, in the MEP pathway upstream of monoterpene biosynthesis showed upregulated expression in Lilium ‘Siberia’ compared to Lilium ‘Novano’. The expression level of DXS-1 in Lilium ‘Siberia’ was far higher (>6,000-fold) than that in Lilium ‘Novano’ (Figure 7). The expression level of GPS in Lilium ‘Siberia’ was 819-fold higher than that in Lilium ‘Novano’. In addition, the genes, DHDDS-1, DHDDS-2, SDS, and GGDR, which did not mediate monoterpenoid biosynthesis downstream of the MEP pathway presented downregulated expression in Lilium ‘Siberia’ (Table 5). In Lilium ‘Siberia’, the expression of DHDDS-1 and DHDDS-2 could not be detected (Figure 7). The gene expression of HMGS and HMGR-2 in upstream of IPP biosynthesis was also upregulated in the MVA pathway (Table 5). The relative mRNA concentration of HMGR-2 in Lilium ‘Siberia’ was nearly 31-fold higher than that in Lilium ‘Novano’ (Figure 7).

Among the 7 DEGs annotated to monoterpene synthases, two unigenes annotated to the ocimene synthase gene (OCS, c65063.graph_c0) and myrcene synthase gene (MYS, c66796.graph_c0) respectively were upregulated in Lilium ‘Siberia’ (Figure 8). The relative mRNA concentration of the OCS in Lilium ‘Siberia’ almost reached 2,500, which was far higher than that in Lilium ‘Novano’. The relative mRNA concentration of the OCS in Lilium ‘Novano’ did not exceed 1. The expression level of MYS showed a similar pattern, which was nearly 12.91-fold in Lilium ‘Siberia’ higher than in Lilium ‘Novano’.

Seven DEGs involved in monoterpene biosynthesis were selected for qRT-PCR analysis to validate the expression profiles obtained by RNA-Seq. The expression patterns of DXS-1, DXR, GPS, HMGR-2, OCS, and MYS of Lilium ‘Siberia’ and Lilium ‘Novano’ at the different flowering stages were shown in Figure 9. The expression of these genes exhibited a similar pattern in both genotypes, i.e., an initial increase and then a decrease with flower development (Figure 9). High expression levels were measured at HS and FS, which was consistent with the amounts of monoterpenes released (Figure 3).

Additionally, it was found that the relative expression levels of these 6 genes in Lilium ‘Siberia’ were higher than those in Lilium ‘Novano’ at the 4 stages, especially DXS-1 and OCS, whose relative expression levels in Lilium ‘Siberia’ at FS were 9526- and 2015-fold higher, respectively, than those in Lilium ‘Novano’. This result was basically consistent with the RNA-Seq results.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.