The CM samples gathered from the medicinal botanical garden at Hubei University of Chinese Medicine were identified and authenticated by Professor Gong Ling of Hubei University of Chinese Medicine. Fresh CM samples were immediately imaged after collection and separated into petal and non-petal tissues for subsequent analysis, as shown in Figure 1 . Non-petal samples were collected from white, yellow, and gold CM, and were designated as groups WCNP, YCNP, and GCNP, respectively. The petal samples from these colors were designated as groups WCP, YCP, and GCP. Fresh samples were stored at -80°C prior to metabolomics analysis and RNA extraction. High-performance liquid chromatography-grade acetonitrile and methanol were acquired from Merck Chemicals, Darmstadt, Germany.

The samples were first freeze-dried to a constant weight and then ground into a fine, uniform powder. To begin the extraction process, 100 mg of the powdered sample was weighed and extracted with 3 mL of an 80% methanol solution using ultrasonic extraction for 10 minutes. This extraction step was repeated twice: first with 3 mL of a 50% methanol solution and then with 3 mL of a 95% methanol solution. The three resulting extracts were combined, and the mixture was centrifuged at 12,000 rpm for 10 minutes. After centrifugation, the supernatant was carefully filtered and prepared for further analysis. Curcumin was employed as an internal standard in the samples at a final concentration of 200 ng mL^-1. After filtration, the quality control (QC) sample was prepared by combining equal aliquots from all individual samples.

For the analysis of the extracts, an ACQUITY UPLC H-Class system from Waters (Milford, MA, USA) was utilized, featuring a Waters ACQUITY UPLC BEH C18 column (2.1 × 100 mm, 1.7 μm). The injection volume was 2.0 μL, the flow rate was maintained at 0.3 mL min^-1, and the column temperature was set to 40°C. The mobile phase consisted of two solvents: mobile phase A, a solution of water and formic acid in a 1000:1 (v/v) ratio, and mobile phase B, acetonitrile. The chromatographic gradient was programmed as follows: 0 min, 5% (B); 12 min, 35% (B); 18 min, 80% (B); 22 min, 95% (B); 25 min, 95% (B); 26 min, 5% (B); and 30 min, 95% (A) and 5% (B).

Mass spectrometry analysis was carried out using a Waters Xevo G2-XS QTOF system equipped with an electrospray ionization source. The parameters were set according to our previously published reports (Zeng et al., 2023, 2024). The optimal parameters were as follows: desolvation temperature of 500°C, cone voltage of 20 V, capillary voltage of 3 KV, source temperature of 100°C, desolvation gas flow rate of 1000 L h^-1, collision energy range of 30 to 40 eV, and cone gas flow rate of 50 L h^-1. The mass range for full scans was set from m/z 50 to 1500 Da, with a scan duration of 1.0 s. Data acquisition was conducted in MS^E mode, employing both positive and negative ion electrospray modes. Metabolomics analysis predominantly utilized positive ion mode, while structural determination of metabolites was achieved with both ion modes.

RNA extraction and Illumina sequencing of the samples were primarily conducted by MetWare Biotechnology Co., Ltd. (Wuhan, Hubei, China). Total RNA was extracted from plant samples using ethanol precipitation and the CTAB-PBIOZOL method, then dissolved in DEPC-treated water. RNA integrity and concentration were assessed using a Qubit fluorescence quantifier and Qsep400 biofragment analyzer (Bioptic Inc., Taiwan, China). PolyA-tailed mRNAs were enriched with Oligo (dT) magnetic beads, fragmented, and reverse-transcribed into first-strand cDNA using random hexamer primers. Strand-specific second-strand cDNA synthesis was performed with dUTPs to ensure strand specificity, followed by end repair, A-tailing, and sequencing adapter ligation. The cDNA was size-selected (250-350 bp), PCR-amplified, purified, and quantified using Qubit 4.0 (Thermo Fisher Scientific, Massachusetts, USA) and Q-PCR (Bio-rad, California, USA). Libraries were pooled based on effective concentrations and sequenced on an Illumina platform, generating 150 bp paired-end reads via sequencing-by-synthesis. Raw sequencing data were filtered with fastp, and clean reads were assembled using Trinity. Corset was then used to remove redundant isoforms from the assembled transcripts. CDS prediction was conducted with TransDecoder, and gene function annotation was performed using DIAMOND and HMMER across databases. The fragments per kilobase of transcript per million mapped reads (FPKM) for each gene were subsequently calculated based on the gene length and the number of reads mapped to it. Transcript expression levels were quantified using RSEM, with differential expression analysis conducted using DESeq2 and edgeR, followed by Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses.

The raw data from UPLC-QTOF-MS/MS were acquired using MassLynx v4.1 (Waters, Milford, MA, USA). Following data acquisition, the open-access software MS-DIAL (Version 4.9) was utilized for comprehensive data processing, encompassing peak detection, alignment, spectral deconvolution, identification, and normalization (Tsugawa et al., 2015). The retention time for the collected peaks was set between 1 and 30 minutes. Peaks with a minimum height of 3000 and m/z values ranging from 100 to 1500 Da were selectively retained, guided by the expected component range, while the m/z value for fragment ions was set between 50 and 1500 Da. Mass tolerance parameters were established at 0.015 Da for MS and 0.02 Da for MS/MS. To determine the elemental composition of the peaks, common positive ion adducts such as [M+H]⁺, [M+Na]⁺, [M+K]⁺, [M+H-H₂O]⁺, and [2M+H]⁺ were employed. In negative ion mode, adducts like [M-H]^-, [M+HCOOH-H]^-, and [M-H₂O-H]^- were used to assist in adduct correction alongside the positive ion adducts. Feature alignment across all samples was performed with a retention time tolerance of 0.1 and an MS tolerance of 0.015. Furthermore, the deconvolution value and MS/MS abundance cutoff were set at 0.6 and 100, respectively, to facilitate peak deconvolution. The response values of the peaks were normalized using internal standards and the LOWESS method. Finally, to ensure that the screened metabolites accurately represented each group, peaks were required to be present in at least one group, with every sample within that group exhibiting the corresponding peak.

In the search for characteristic differentially accumulated metabolites (DAMs) between two groups, we utilized variable importance in projection (VIP) values exceeding 1.25 and fold change (FC) thresholds greater than 1.5 or less than 0.67 to identify potential distinctive peaks. In the transcriptome analysis, to ensure the genes were both representative and meaningful, we applied the following criteria to eliminate false positives: (1) the average FPKM value of the gene in at least one group exceeds 1; (2) the gene’s detection rate across all samples is greater than 0.167; (3) the gene’s detection rate in any individual group is above 0.667. For the differentially expressed genes (DEGs) identified in the comparison between the two groups, the initial requirement was that genes had an FPKM > 1 in at least one group and were detected in at least two samples within any group. Subsequently, based on FPKM quantitative data, the fold change (FC) for the DEGs had to exceed 2 or fall below 0.5, with a q-value between the two groups required to be < 0.05. In the co-enrichment analysis, to identify pathways that represent both metabolomics and transcriptomics, we selected pathways containing at least two DAMs and conducted a comprehensive analysis of the DEGs upstream and downstream of the DAMs within these pathways.

The SIMCA 14.1 software (Umetrics AB, Umeå, Sweden) was employed for principal component analysis (PCA) and orthogonal partial least squares-discriminant analysis (OPLS-DA). Data visualization was carried out using Origin 2021 (OriginLab Corporation, Northampton, MA, USA) and R version 4.3.2 (R Foundation for Statistical Computing, Vienna, Austria). Statistical analyses were performed with SPSS version 26.0 (SPSS, Inc., Chicago, IL, USA). The significance of DAMs was assessed with a P-value < 0.05 using the Mann-Whitney U test, while DEGs were evaluated with a q-value < 0.05 using the Benjamini-Hochberg test. Qualitative analysis of DAMs was conducted by comparing fragment ions using literature references (Chen et al., 2017), MS-FINDER version 3.60 (Tsugawa et al., 2016), and public databases such as the Human Metabolome Database (https://www.hmdb.ca) and MassBank (https://massbank.eu/MassBank/). Quantitative determination of metabolites was achieved using Quanlynx software version 4.1 (Waters, Milford, MA, USA). For the co-enrichment analysis of metabolomics and transcriptomics, metabolic pathways of DAMs and DEGs were primarily constructed with reference to KEGG.

In this study, a total of 6704 features were acquired for all samples using MS-DIAL based on UPLC-QTOF-MS/MS. As shown in Figure 2A , the PCA model was applied for unsupervised analysis of all samples, including quality controls. The first two principal components accounted for 63.30% and 9.74% of the variance among the samples, respectively. The results revealed clear intra-group clustering, demonstrating strong reproducibility within each sample group. Additionally, the WCNP and YCNP groups displayed a higher degree of global similarity to each other. The clustering of the QC sample further confirmed the robust stability of instrument throughout the analysis.

Furthermore, supervised OPLS-DA was employed to identify DAMs between the two groups in both petal or non-petal areas (see Supplementary Figures S1 and S2 ). The scatter plots derived from OPLS-DA revealed clear segregation between each pair of groups. All pairwise comparisons yielded R2Y and Q2 scores consistently exceeding 0.9. Additionally, the results of the 999 permutation tests showed that the blue regression line of Q2 intersected the vertical axis below zero, confirming the absence of overfitting in the original OPLS-DA model. Ultimately, 90 DAMs (P < 0.05) were characterized, with a higher number identified in the petals than in the non-petal areas. The qualitative information of these metabolites, including retention time, actual m/z mass of the quasi-molecular ion peak, tentative identification, molecular formula, secondary fragments, and classification details, was provided in Supplementary Table S1 . The FC and VIP values of DAMs across different comparison groups were listed in Supplementary Table S2 . When comparing various petal groups, WCP vs. YCP had the fewest number of DAMs, while WCP vs. GCP had the most ( Figure 2B ). Meanwhile, WCNP vs. GCNP had the highest number of DAMs among the various non-petal groups. As illustrated in Figure 2C , the identified DAMs were mainly categorized into 6 categories: amino acids, flavonoid glycosides, flavonoids, lipids, other glycosides, and others. Among these, flavonoid glycosides, lipids, and flavonoids were predominantly classified in the petal parts, whereas flavonoid glycosides and lipids were primarily classified in the non-petal parts.

To elucidate variations in the expression levels of these DAMs across different morphologies of CM, we conducted a relative quantitative analysis of the corresponding groups with the ratio of peak area between DAMs and internal standards, as shown in Figure 3 . In the petals of CM, the following metabolite changes were observed: 8 up-regulated and 21 down-regulated in the WCP vs. YCP, 14 up-regulated and 20 down-regulated in YCP vs. GCP, and 12 up-regulated and 25 down-regulated in WCP vs. GCP. Similarly, in the non-petals of CM, 16 up-regulated and 6 down-regulated metabolites were identified in the WCNP vs. YCNP, 5 up-regulated and 17 down-regulated in YCNP vs. GCNP, and 11 up-regulated and 16 down-regulated in WCNP vs. GCNP.

As highlighted in Figure 3 , most flavonoid glycosides showed higher concentrations in white CM, both in non-petal and petal parts. However, within the flavonoids, the petal part of white CM exhibited lower contents, while the non-petal part showed higher contents. Additionally, most lipids within DAMs were found at lower concentrations in both WCP and WCNP compared to the other two groups. Conversely, in the petals, most lipids were more abundant in YCP, whereas in the non-petal parts, they were predominantly higher in GCNP.

To gain a more comprehensive understanding of the non-petal and petal parts across different phenotypes, transcriptomic sequencing and variance analysis of DEGs were employed to investigate the underlying molecular mechanisms. After data acquisition and transcript assembly, a total of 231670 raw sequencing reads were retained. To reduce false positives while preserving valid data, 181404 reads were selected for further transcriptomic analysis. The PCA results revealed that both yellow and gold CM showed similar trends in non-petal and petal parts compared to white CM ( Figure 4A ). Moreover, each sample displayed a distinct clustering pattern within its respective group. The first two principal components accounted for 23.16% and 13.49% of the total variance among the samples, respectively. DEGs between the two groups were identified using criteria of a FC greater than 2 or less than 0.5, and a q-value less than 0.05. As illustrated in Figure 4B , the comparisons between WCP and GCP, and WCNP and GCNP, revealed the highest number of DEGs in the petal and non-petal groups, respectively. Additionally, fewer transcriptional differences were observed between yellow and gold CM.

To further illustrate the differences between the groups, a heat map was generated from the quantitative data of DEGs, as detailed in Supplementary Figure S3 . Statistical analysis of DEG expression trends between the two groups, presented in Figure 4C , revealed that most DEGs exhibited a decreasing trend, particularly in comparisons between WCNP vs. YCNP and WCNP vs. GCNP.

To systematically explore the biological functions of DEGs in petal and non-petal parts across different flower colors, we performed KEGG pathway enrichment analysis on the DEGs from each group separately. As shown in Figure 5 , the top fifteen pathways enriched by DEGs in each comparison group were selected and visualized in a bubble chart.

These pathways fall into several categories, including cellular processes, environmental information processing, genetic information processing, metabolism, and organismal systems, with metabolism representing the largest portion (78.95%). Within the metabolism category, pathways related to carbohydrate and lipid metabolism were the most prevalent. Moreover, the three metabolic pathways with the highest number of DEGs were phenylpropanoid biosynthesis, pentose and glucuronate interconversions, and glycerophospholipid metabolism.

To investigate the genetic regulation of CM flower color, this study conducted a comprehensive analysis of pathways involving DEGs and DAMs, building on prior metabolomic and transcriptomic data. In order to identify pathways representing both the transcriptome and metabolome, only metabolic pathways containing at least two DAMs were retained for further analysis, with DEGs included only if located upstream or downstream of the DAMs. The KEGG codes and pathways, along with the numbers and lengths of the DEGs involved in the core pathways, were displayed in Supplementary Table S3 .

As shown in Figure 6 , glycerophospholipid metabolism and glycerolipid metabolism, both primarily involving lipids, were the primary pathways associated with CM flower color. In petal tissues, the majority of DEGs in WCP exhibited the lowest expression compared to other groups, particularly EPT1, DAD1, MGD, and psd ( Figure 6A ). However, at the metabolite level, only MGDG, MGMG, and LysoPC (sn-1) showed reduced levels in WCP. Additionally, compared to YCP, a greater number of DAMs in GCP displayed a downward trend, including LysoPC (sn-2), PC, PE, and LysoPE (sn-2).

Likely, in non-petal tissues, a significant number of DEGs, including MGLL, EPT1, LPCAT1_2, and DAD1, showed reduced expression in WCNP ( Figure 6B ). Unlike the DAMs in petal tissues, the majority of DAMs in WCNPs exhibited significantly lower levels, particularly LysoPC (sn-1), LysoPC (sn-2), LysoPA (sn-1), and LysoPE (sn-2). Moreover, a substantial number of DAMs in GCP showed an upward trend compared to YCP, which contrasts with the pattern observed in the petal tissues.