RNA was extracted from the two lily variants using ethanol precipitation combined with the CTAB-PVIOZOL method. The extracted RNA was dissolved in 50 µL of DEPC-treated sterile water. Quality and quantity assessments were performed using a Qubit fluorometer and a Qsep400 high-throughput bio-fragment analyzer. Qualified total RNA samples from L. cernuum (labeled as R-t) and L. cernuum var. album (labeled as W-t) underwent transcriptome sequencing, each with three biological replicates, resulting in six sample groups. Sequencing was outsourced to Wuhan Maiwei Metabolic Biotechnology Co., Ltd., utilizing the Illumina platform.

To ensure high-quality data, raw reads were filtered using fastp (v0.23.2) to remove adapter-containing sequences. Clean reads were assembled into transcripts using Trinity (v2.12.2), followed by clustering and de-redundancy with Corset (v1.09) (https://github.com/trinityrnaseq/trinityrnaseq). Coding sequences (CDS) were predicted from the Trinity-assembled transcripts using TransDecoder (TransDecoder5.3.0) (https://github.com/TransDecoder/), yielding corresponding amino acid sequences. De-redundant transcript sequences were annotated by aligning them with databases such as Nr (http://www.ncbi.nlm.nih.gov/), GO (http://www.geneontology.org), SwissProt (http://www.expasy.ch/sprot/), KEGG (http://www.genome.jp/kegg/), and EggNog (http://eggnogdb.embl.de/) using DIAMOND (v5.3.0). Transcription factors were predicted using iTAK software (Zheng et al., 2016).

Differential expression analysis between sample groups was conducted using DESeq2 to identify genes exhibiting significant differences in expression between the two biological conditions. To account for multiple hypothesis testing, the Benjamini-Hochberg method was applied to adjust the p-values, resulting in the calculation of the false discovery rate (FDR). Genes were then filtered based on the criteria of |log2 fold change| ≥ 1 and FDR < 0.05 to identify significantly differentially expressed genes.

qRT-PCR was employed to validate the transcriptional levels of differentially expressed genes (DEGs) in the petals of both lily variants. Reactions were performed using SYBR Premix Ex Taq (Tli RnaseH Plus, TaKaRa), with cDNA (diluted 20-fold) serving as the template. Protocols and reaction conditions adhered to the manufacturer’s instructions. Each qRT-PCR assay was conducted with three biological replicates, and relative gene expression levels were calculated using the 2^−ΔΔCT method. Subsequently, a one-way analysis of variance (ANOVA) was performed to compare gene expression differences across groups. To mitigate the risk of false positives from multiple comparisons, Tukey’s Honest Significant Difference (HSD) test was applied for pairwise comparisons, while the Benjamini-Hochberg method was utilized to control the false discovery rate (FDR). Gene-specific primers for qRT-PCR were designed using Primer Express (v3.0), with common lily actin (LhActin, NCBI accession number: JX826390) serving as the internal control. Details of the specific and internal control primers are listed in Table 1 .

The lily petals stored in a -80°C freezer were taken out and ground into powder in liquid nitrogen. This process was performed for three biological replicates of each lily type, resulting in a total of six samples. The L. cernuum samples were labeled as R-GC, and the L. cernuum var. album samples were labeled as W-GC. 500 mg of the powder was immediately transferred to a 20 mL headspace vial (Agilent, Palo Alto, CA, USA) containing a NaCl saturated solution. The vials were sealed using crimp-top caps with TFE-silicone headspace septa (Agilent). For the SPME analysis, each vial was placed at 60°C for 5 minutes, then a 120 µm DVB/CWR/PDMS fiber (Agilent) was exposed to the headspace of the sample for 15 minutes at 60°C.

After sampling, desorption of the VOCs from the fiber coating was carried out in the injection port of the GC apparatus (Model 8890; Agilent) at 250°C for 5 minutes in splitless mode. The identification and quantification of VOCs were performed using an Agilent Model 8890 GC coupled with a 7000D mass spectrometer (Agilent), equipped with a 30 m × 0.25 mm × 0.25 μm DB-5MS (5% phenyl-polymethylsiloxane) capillary column. Helium was used as the carrier gas at a linear velocity of 1.2 mL/min. The injector temperature was maintained at 250°C. The oven temperature was programmed as follows: starting at 40°C (held for 3.5 minutes), increasing at 10°C/min to 100°C, at 7°C/min to 180°C, and at 25°C/min to 280°C, where it was held for 5 minutes. Mass spectra were recorded in electron impact (EI) ionization mode at 70 eV. The quadrupole mass detector, ion source, and transfer line temperatures were set to 150°C, 230°C, and 280°C, respectively. The MS was operated in selected ion monitoring (SIM) mode for the identification and quantification of analytes.

Differential metabolites between L. cernuum and L. cernuum var. album were initially identified based on the Variable Importance in Projection (VIP) scores derived from the OPLS-DA model. These metabolites were then further refined using the false discovery rate (FDR) obtained from univariate analysis. The screening criteria for significant differences were defined as VIP > 1 and either a fold change ≥ 2 or ≤ 0.5.

The lily petals stored in a -80°C freezer were taken out and ground into powder in liquid nitrogen. This procedure was performed for three biological replicates of each lily type, resulting in a total of six samples. The L. cernuum samples were labeled as R-MS, and the L. cernuum var. album samples were labeled as W-MS. Using vacuum freeze-drying technology, the biological samples were placed in a lyophilizer (Scientz-100F) and then ground (30 Hz, 1.5 minutes) into powder using a grinder (MM 400, Retsch). Next, 50 mg of sample powder was weighed using an electronic balance (MS105DM) and 1200 μL of -20°C pre-cooled 70% methanolic aqueous internal standard extract was added (for samples weighing less than 50 mg, the extract was added at a rate of 1200 μL per 50 mg sample). The mixture was vortexed for 30 seconds every 30 minutes, for a total of six times. After centrifugation (12,000 rpm, 3 minutes), the supernatant was aspirated, filtered through a microporous membrane (0.22 μm pore size), and stored in an injection vial for UPLC-MS/MS analysis.

The sample extracts were analyzed using a UPLC-ESI-MS/MS system (UPLC, ExionLC™ AD, https://sciex.com.cn/) coupled with a tandem mass spectrometry system (https://sciex.com.cn/). The analytical conditions were as follows: UPLC column, Agilent SB-C18 (1.8 µm, 2.1 mm × 100 mm); the mobile phase consisted of solvent A (pure water with 0.1% formic acid) and solvent B (acetonitrile with 0.1% formic acid). Sample measurements were performed using a gradient program starting at 95% A and 5% B. Over 9 minutes, a linear gradient to 5% A and 95% B was programmed, and this composition was maintained for 1 minute. Subsequently, the composition was adjusted to 95% A and 5% B within 1.1 minutes and held for 2.9 minutes. The flow rate was set to 0.35 mL/min; the column oven was set to 40°C; the injection volume was 2 μL. The effluent was alternately connected to an ESI-triple quadrupole-linear ion trap (QTRAP)-MS.

The ESI source operation parameters were as follows: source temperature 550°C; ion spray voltage (IS) 5500 V (positive ion mode)/-4500 V (negative ion mode); ion source gas I (GSI), gas II (GSII), and curtain gas (CUR) were set at 50, 60, and 25 psi, respectively; and the collision-activated dissociation (CAD) was set to high. QQQ scans were acquired as MRM experiments with collision gas (nitrogen) set to medium. The DP (declustering potential) and CE (collision energy) for individual MRM transitions were optimized. A specific set of MRM transitions was monitored for each period according to the metabolites eluted during that period.

The ESI source operation parameters were as follows: source temperature 550°C; ion spray voltage (IS) 5500 V (positive ion mode)/-4500 V (negative ion mode); ion source gas I (GSI), gas II(GSII), curtain gas (CUR) were set at 50, 60, and 25 psi, respectively; the collision-activated dissociation(CAD) was high. QQQ scans were acquired as MRM experiments with collision gas (nitrogen) set to medium. DP(declustering potential) and CE(collision energy) for individual MRM transitions was done with further DP and CE optimization. A specific set of MRM transitions were monitored for each period according to the metabolites eluted within this period.

The procedure for differential metabolite screening was identical to that outlined in Section 2.3.3.

Unsupervised PCA (principal component analysis) was performed using the ‘prcomp’ statistics function within R (www.r-project.org). The data was unit variance scaled before performing the unsupervised PCA.

The HCA (hierarchical cluster analysis) results of samples and metabolites were presented as heatmaps with dendrograms, while Pearson correlation coefficients (PCC) between samples were calculated using the ‘cor’ function in R and presented as heatmaps only. Both HCA and PCC were performed using the R package ‘ComplexHeatmap’. For HCA, normalized signal intensities of metabolites (unit variance scaling) were visualized as a color spectrum.

For two-group analysis, differential metabolites were determined by VIP (VIP > 1) and absolute Log2FC (|Log2FC| ≥ 1.0). VIP values were extracted from OPLS-DA results, which also contained score plots and permutation plots, generated using the R package ‘MetaboAnalystR’. The data was log-transformed (log2) and mean-centered before OPLS-DA. To avoid overfitting, a permutation test (200 permutations) was performed.

Identified metabolites were annotated using the KEGG Compound database (http://www.kegg.jp/kegg/compound/), and annotated metabolites were then mapped to the KEGG Pathway database (http://www.kegg.jp/kegg/pathway.html). Pathways with significantly regulated metabolites were then fed into MSEA (metabolite sets enrichment analysis), and their significance was determined by hypergeometric test p-values.

Through transcriptomic analysis of the petals of two lily color varieties, R-t generated 46,676,986 clean reads, and W-t generated 42,606,872 clean reads. The error rate was kept below 0.02%. The Q30 values ranged between 94.70% and 95.31%, and the GC content ranged between 49.47% and 50.60%. The PCA plot showed that PC1 and PC2 contributed 63.07%, with the two lily varieties clearly separated ( Figure 2a ). The annotation rate reached 73.81% after comparing the reads with six databases: KEGG, NR, SwissProt, KOG, GO, and Pfam, indicating that the transcriptomic analysis is reliable and can be used for subsequent analysis. Differential gene screening of the obtained clean reads revealed 12,174 differentially expressed genes (DEGs) between R-t and W-t, with 6,055 genes downregulated and 6,119 genes upregulated ( Figure 2b ). According to K-means clustering analysis, the DEGs were divided into 10 subclasses. Compared to R-t, W-t showed higher gene expression levels in subclasses 1, 6, 7, and 8, while in subclasses 2, 3, 4, 5, 9, and 10, W-t exhibited lower gene expression levels compared to R-t ( Figure 2c ). These genes may play a crucial role in flower color change.

The GO annotation statistics of differential genes revealed several key biological processes and molecular functions related to flower color variation during pigment biosynthesis. In the Biological Process (BP) category, differential genes were mainly concentrated in Cellular and Metabolic Processes. The Cellular Component (CC) category data indicated significant changes in molecular components within organelles related to pigment synthesis and storage, possibly implying regulation of plastid formation or function in petals. In the Molecular Function (MF) category, genes involved in catalytic activation and protein binding were highly enriched, suggesting their potential roles in regulating anthocyanin biosynthesis and transport ( Figure 3a ). The KEGG annotation of differential genes showed distinct differential expressions in multiple metabolic pathways. These differential genes were mainly concentrated in Flavonoid biosynthesis (Ko00941), Biosynthesis of secondary metabolites (Ko01110), and Metabolic pathways (Ko01100) in the Metabolism process. Among these, the RichFactor value was the highest in the Flavonoid biosynthesis pathway. In this pathway, 58 differential genes were upregulated and 13 were downregulated in R-t compared to W-t ( Figure 3b ). In flavonoid biosynthesis, these differentially expressed genes mainly belong to the transferase family, reductase family, synthetase family, isomerase family, and oxidoreductase family. Key genes include flavanone 3-hydroxylase and chalcone synthase. Among them, genes related to the chalcone and stilbene synthases, N-terminal domain account for 26.8%, playing a crucial role in the biosynthesis of flavonoids and anthocyanin pigments.

This study performed an overlay analysis of the total ion chromatograms (TIC) obtained from mass spectrometry detection of various quality control samples. The results showed a high degree of overlap in the TIC curves, with consistent retention times and peak intensities. This consistency indicates that the mass spectrometry signals remained stable across different detection times for the same sample, confirming their reliability for subsequent analysis. A total of 540 VOCs were detected in the petals, which were classified into 16 categories. Among these, terpenoids were the most abundant, with 97 types accounting for 18% of the total. This was followed by esters, with 84 types accounting for 15.6%, and heterocyclic compounds, with 80 types accounting for 14.8%. Together, these three classes of compounds accounted for 48.4% of the total ( Figure 4a ). Comparing the fold changes in VOCs quantification between W-GC and R-GC, 120 of the 540 VOCs showed changes. Compared to R-GC, 110 VOCs were significantly downregulated in W-GC, while only 10 VOCs were significantly upregulated ( Figure 4b ).

KEGG annotation was performed on the 120 VOCs, and 21 of these VOCs were annotated to 11 metabolic pathways. The classification of these results is shown in Figure 5 . Among the annotated VOCs, 12 were concentrated in the Metabolic pathways process, accounting for 57.14%, followed by the Biosynthesis of secondary metabolites process, which included 5 significantly different VOCs, accounting for 23.81% ( Figure 5a ). Based on the results of the differential metabolites, KEGG pathway enrichment analysis was performed. The metabolic pathways were ranked and displayed according to their P-value, from smallest to largest across the 11 metabolic processes ( Figure 5b ). The results indicated that in the Phenylpropanoid biosynthesis process, the enrichment index of volatile differential metabolites was relatively high, and the difference was extremely significant.

By comparing the differential metabolites identified based on the screening criteria in W-GC and R-GC and their annotated sensory flavor characteristics, a radar chart was drawn for the top 10 sensory flavors with the highest number of annotations ( Figure 6 ). The results showed that the differential VOCs between W-GC and R-GC were mainly concentrated in sweetness and floral aroma. In sweetness, 22 VOCs were annotated, including Acetophenone, 4’-hydroxy-, BenzAldehyde, BenzAldehyde, 2-methyl-, Hydrocinnamic Acid and Methyleugenol, with five VOCs being annotated to the Metabolic pathways, Phenylalanine metabolism, and Phenylpropanoid biosynthesis pathways, all of which were reduced in content. In floral aroma, 21 VOCs were annotated, including BenzeneacetAldehyde, Benzeneacetic acid, Eugenol, and trans-Isoeugenol, with four VOCs also annotated to the Metabolic pathways, Phenylalanine metabolism, and Phenylpropanoid biosynthesis pathways, all of which were also reduced in content. This suggests that during the transition of L. cernuum from red to white, some of the sweetness and floral aroma characteristics change.

The outermost circle names represent sensory flavor characteristics. The numbers corresponding to the green points indicate the number of occurrences of each sensory flavor characteristic, i.e., the number of differential metabolites annotated to that sensory flavor characteristic.

The overlay analysis of the TIC revealed a high degree of overlap in the curves for metabolite detection, with consistent retention times and peak intensities. This demonstrates excellent signal stability, confirming that the data is reliable and suitable for subsequent analysis. A total of 1,243 metabolites were detected in the petals, which were classified into 13 categories. Among these, Flavonoids were the most abundant, with 319 metabolites accounting for 25.7%, followed by Phenolic acids with 174 metabolites, accounting for 14.0% ( Figure 7a ). Comparing the fold change in metabolite quantification between R-MS and W-MS, there were 419 significantly different metabolites (DMs) among the 1,243 metabolites in the petals. Compared to R-MS, W-MS had 235 metabolites that were significantly down-regulated, of which 119 were Flavonoids, accounting for 50.64%. On the other hand, 184 metabolites were significantly up-regulated ( Figure 7b ), among which Lipids were the most abundant, reaching 50 metabolites, accounting for 42.02%.

KEGG pathway enrichment analysis was performed on the 419 DMs ( Figure 8a ), and the top 20 pathways ranked by P-value were displayed from smallest to largest. The results indicated that DMs related to the synthesis of anthocyanins and flavonoids were mainly concentrated in the Arginine biosynthesis (ko00942) and Flavone and flavonol biosynthesis (ko00944) pathways. Cluster analysis of the DMs in these two metabolic pathways revealed that in Arginine biosynthesis ( Figure 8b ), Delphinidin-3,5,3’-Tri-O-glucoside and Delphinidin-3-O-glucoside (Mirtillin) showed significantly increased content in W-MS, while they were significantly decreased in R-MS. Conversely, Cyanidin 3-O-sophoroside, Pelargonidin-3-O-glucoside, and Peonidin-3-O-glucoside showed significantly decreased content in W-MS and significantly increased in R-MS. In the Flavone and flavonol biosynthesis pathway ( Figure 8c ), Vitexin-2’’-O-rhamnoside, Kaempferol-3-O-rhamnoside (Afzelin) (Kaempferin), Quercetin-3-O-rhamnoside (Quercitrin), and Kaempferol-3-O-rutinoside (Nicotiflorin) showed significantly increased content in R-MS and significantly decreased in W-MS.

In this study, the correlation between DEGs (Differentially Expressed Genes) and VOCs (Volatile Organic Compounds) was analyzed. By calculating the Pearson correlation coefficient between DEGs and VOCs, we selected those with an absolute correlation coefficient greater than 0.8 and a p-value less than 0.05. Then, based on the KEGG annotation results for genes and volatile metabolites, the relevant data were screened ( Figure 9a ). The genes and VOCs in the Phenylpropanoid biosynthesis (ko00940) and Phenylalanine metabolism (ko00360) pathways were selected for analysis. In these two pathways, when comparing R-GC and W-GC, the content of Eugenol and Methyleugenol in W-GC significantly decreased. There were 66 genes associated with these two metabolites, 59 of which were common to both Eugenol and Methyleugenol. Comparing these 59 genes with the PFam database, they mainly belong to three protein families: Peroxidase (13 genes), Aromatic amino acid lyase (10 genes), and AMP-binding enzyme (7 genes). Comparing with the Nr database, the 59 associated genes were mainly annotated to five gene families: Peroxidase Family (8 genes), Cinnamoyl-CoA Reductase Family (8 genes), 4-Coumarate–CoA Ligase Family (4 genes), Phenylalanine Ammonia-Lyase Family (3 genes), and Caffeoyl-CoA O-Methyltransferase Family (3 genes). All 59 associated genes in W-GC were down-regulated ( Figure 9b ).

The correlation between DEGs and DMs in L.cernuum was analyzed using their quantitative values. By calculating the Pearson correlation coefficient between DEGs and DMs, we selected those with an absolute correlation coefficient greater than 0.8 and a p-value less than 0.05. Then, based on the KEGG annotation results for genes and metabolites, relevant data were screened. As shown in Figure 10 , the number of DEGs and DMs in the Anthocyanin biosynthesis (ko00942) pathway reached a highly significant level. Additionally, DEGs and DMs in the Phenylpropanoid biosynthesis (ko00940) and Flavonoid biosynthesis (ko00941) pathways were also selected for analysis.

When comparing R-MS and W-MS, in the Phenylpropanoid biosynthesis pathway, 37 genes and the metabolites Chlorogenic acid (3-O-Caffeoylquinic acid), Sinapinaldehyde, and p-Coumaric acid showed changes, with both gene expression and metabolite content being down-regulated in W-MS. These 37 genes belong to the Phenylpropanoid Pathway Enzymes, Carbohydrate Metabolism Enzymes, and Aldehyde Dehydrogenase Family categories, with the Phenylpropanoid Pathway Enzymes mainly including the Phenylalanine Ammonia Lyase (PAL), 4-Coumarate–CoA Ligase (4CL), and Cinnamoyl-CoA Reductase (CCR) families. It is worth noting that these 37 genes are also associated with the VOCs Eugenol and Methyleugenol.

In the Flavonoid biosynthesis pathway, a larger number of genes showed changes. A total of 71 genes and 7 metabolites, including Aromadendrin (Dihydrokaempferol), Catechin, Chlorogenic acid (3-O-Caffeoylquinic acid), Epicatechin, Hesperetin, Hesperetin-7-O-neohesperidoside (Neohesperidin), and Naringenin-7-O-glucoside (Prunin), showed changes. Compared to R-MS, 58 of the 71 genes in W-MS were down-regulated, and 13 were up-regulated. Among the 7 metabolites, the content of Aromadendrin and Hesperetin increased, while the other 5 metabolites decreased. These 71 genes mainly belong to the Flavonoid Biosynthesis Enzymes and Phenylpropanoid Pathway Enzymes categories, with the Flavonoid Biosynthesis Enzymes primarily including Chalcone Synthase (CHS), Chalcone Isomerase (CHI), Flavanone 3-Hydroxylase (F3H), Flavonol Synthase (FLS), Flavonoid 3’-Hydroxylase (F3’H), Dihydroflavonol 4-Reductase (DFR), Anthocyanidin Synthase (ANS), and Anthocyanidin Reductase (ANR).

In the Anthocyanin biosynthesis pathway, 7 genes and 9 metabolites showed changes. Compared to R-MS, 6 of the 7 genes in W-MS were down-regulated, and 2 were up-regulated. Among these 7 genes, 5 belonged to UDP-glucose: anthocyanidin 3-O-glucosyltransferase (UFGT), while 2 were unclassified. The 9 metabolites can be categorized into 4 types, with the content of Cyanidin, Pelargonidin, and Peonidin decreasing in W-MS, but the content of Delphinidin increasing.

The horizontal axis represents the KEGG pathway names, and the vertical axis represents the significance test p-value for enrichment in the pathway. Red and green indicate metabolomics and transcriptomics, respectively.