Biological samples were freeze-dried using a vacuum freeze-dryer (Scientz-100F). The freeze-dried samples were then crushed using a mixer mill (MM 400, Retsch) with a zirconia bead for 1.5 minutes at 30 Hz. A 100 mg portion of the lyophilized powder was dissolved in 1.2 ml of 70% methanol solution. The mixture was vortexed for 30 seconds every 30 minutes, repeated six times, and stored at 4°C overnight. After centrifugation at 12,000 rpm for 10 minutes, the extracts were filtered through a 0.22 μm pore size filter (SCAA-104; ANPEL, Shanghai, China) before undergoing UPLC-MS/MS analysis (Luo et al., 2013).

For the rutin standard solution, 0.05 g of rutin was dissolved in anhydrous ethanol in a 50 mL volumetric flask. Aluminum nitrate solution was prepared by dissolving 17.60 g of aluminum nitrate nonahydrate in water, and potassium acetate solution was prepared by dissolving 9.814 g of potassium acetate in water, with both solutions adjusted to 100 mL. Calibration was performed by placing 1.0 to 5.0 mL of the rutin solution in separate 50 mL volumetric flasks and filling them to 15 mL with ethanol. Subsequently, 1.00 mL each of the aluminum nitrate and potassium acetate solutions was added, and the final volume was adjusted with water. After standing for 1 hour, absorbance was measured at 420 nm. C. paliurus leaves were freeze-dried for 24 hours, and then 1.50 g of the dried sample was ground. The sample was extracted with 30 mL of ethanol at 40°C using ultrasound for 1 hour, with shaking every 20 minutes. The solution was then filtered and brought to a final volume of 50 mL with ethanol. For analysis, 1.0 mL of this extract was diluted to 15.0 mL with ethanol, and the flavonoid content was measured using the standard curve (Liu et al., 2016).

The sample extracts were analyzed using a UPLC-ESI-MS/MS system (UPLC, SHIMADZU Nexera X2; MS, Applied Biosystems 4500 Q TRAP). The UPLC utilized an Agilent SB-C18 column (1.8 µm, 2.1 mm x 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). The gradient program started with 95% A and 5% B, transitioning linearly to 5% A and 95% B within 9 minutes, maintaining this ratio for 1 minute. The program then reverted to 95% A and 5% B over 1.1 minutes, held for 2.9 minutes. The flow rate was 0.35 ml/min, the column oven was set to 40°C, and the injection volume was 4 μl. The effluent was connected to an ESI-triple quadrupole-linear ion trap (QTRAP)-MS. LIT and triple quadrupole (QQQ) scans were acquired using an AB4500 Q TRAP UPLC/MS/MS System with an ESI Turbo Ion-Spray interface, operating in both positive and negative ion modes, controlled by Analyst 1.6.3 software (AB Sciex). The ESI source parameters were: source temperature 550°C, ion spray voltage 5500 V (positive)/-4500 V (negative), ion source gas I and II at 50 and 60 psi, respectively, curtain gas at 25 psi, and high collision-activated dissociation (CAD). Instrument tuning and mass calibration were performed with 10 and 100 μmol/L polypropylene glycol solutions in QQQ and LIT modes. QQQ scans were acquired as MRM experiments with nitrogen collision gas set to medium. Specific MRM transitions were monitored for each period based on the eluted metabolites (Zheng et al., 2022).

In the transcriptome data collected across four different harvesting periods from August to November, we focused on the expression analysis of differentially expressed genes (DEGs) in three specific comparisons: 8MML_vs_9MML, 9MML_vs_10MML, and 10MML_vs_11MML. DEGs were identified by comparing expression levels between adjacent months using DEGseq software, with the screening criteria set to a p-value of less than 0.01 and a log2 fold change greater than 1.5. The screening results are presented in the form of volcano plots, as shown in Figures 3A, C, E . The log2 fold change values on the horizontal axis indicate upregulated genes (greater than 0) and downregulated genes (less than 0), while the vertical axis represents the significance of differential gene expression between the two sample groups. Specifically, in the 8MML_vs_9MML group, we identified 1,226 upregulated and 1,184 downregulated genes ( Figure 3A ). In the 9MML_vs_10MML group, 5 genes were upregulated and 38 genes were downregulated ( Figure 3C ). In the 10MML_vs_11MML comparison, 413 genes were upregulated and 59 genes were downregulated ( Figure 3E ).

Further analysis was performed using Gene Ontology (GO) enrichment for the DEGs identified in the 8MML_vs_9MML, 9MML_vs_10MML, and 10MML_vs_11MML comparisons, with the results illustrated in Figures 3B, D, F . The GO enrichment analysis highlighted significant biological processes, with most differentially expressed genes (DEGs) annotated to processes primarily involved in metabolic and cellular functions. This underscores the essential role of these processes in regulating flavonoid biosynthesis and related metabolic pathways across the different seasonal stages. In the biological process category, DEGs were predominantly enriched in metabolic and cellular processes, suggesting that these genes are crucial in regulating various metabolic pathways during the seasonal fluctuations in flavonoid production. The cellular component category showed a high concentration of DEGs associated with the cell membrane, particularly in processes related to the formation of cell membrane components. This finding indicates the importance of cellular structure and membrane dynamics in maintaining metabolic functions during the seasonal transitions. Moreover, in the molecular function category, DEGs were most enriched in binding and catalytic activities, further emphasizing the role of these genes in driving key biochemical reactions involved in flavonoid biosynthesis. These results not only highlight the key molecular functions associated with flavonoid metabolism but also provide a deeper understanding of the regulatory mechanisms underlying the seasonal changes observed in C. paliurus. The GO enrichment analysis thus reveals important insights into the molecular framework governing the plant’s adaptive responses, reinforcing the seasonal modulation of flavonoid biosynthesis.

Subsequently, KEGG pathway enrichment analysis was conducted for the DEGs identified in the 8MML_vs_9MML, 9MML_vs_10MML, and 10MML_vs_11MML comparisons, with the results depicted in Supplementary Figures S1A–C . In the 8MML_vs_9MML group, a total of 569 DEGs were enriched, with a significant concentration in metabolism-related pathways, including plant-pathogen interaction, starch and sucrose metabolism, and notably, 9 DEGs related to the flavonoid biosynthesis pathway ( Supplementary Tables S6 ), including CYP90A1, cinnamoyl-CoA reductase, VSR6, Flavonoid 3’-monooxygenase, CYP714C2, CYP71AN24, Chalcone and stilbene synthases, CYP87A3 and CYP86A22. In contrast, in the 9MML_vs_10MML group, 38 DEGs were enriched, with the majority being related to plant-pathogen interaction, and 2 DEGs were annotated to the flavonoid biosynthesis pathway ( Supplementary Tables S7 ), including CYP71D9 and shikimate O-hydroxycinnamoyltransferase. Similarly, in the 10MML_vs_11MML group, 420 DEGs were enriched, with a strong emphasis on metabolic and organismal systems pathways, including 4 DEGs related to the flavonoid biosynthesis pathway ( Supplementary Tables S8 ), including stemmadenine O-acetyltransferase, flavonoid 3’-monooxygenase, Isoflavone reductase and trans-cinnamate 4-monooxygenase.

The comprehensive analysis of DEGs in C. paliurus leaves across the critical periods from August to November identified numerous genes implicated in both the upstream biosynthesis and downstream metabolism of flavonoids. Key enzymes in the flavonoid biosynthesis pathway, such as Chalcone Synthase (CHS), Chalcone Isomerase (CHI), and Flavanone 3-Hydroxylase (F3H), were prominently upregulated during the transition from August to September, suggesting an early surge in flavonoid production. Additionally, Dihydroflavonol 4-Reductase (DFR) and Anthocyanidin Synthase (ANS) were significantly upregulated from September to October, correlating with increased production of anthocyanins and other flavonoid derivatives. Notably, Flavonol Synthase (FLS), which facilitates the conversion of dihydroflavonols to flavonols, was also upregulated in this period, indicating a diversification of flavonoid types. In the later months, from October to November, genes such as Leucoanthocyanidin Reductase (LAR) and Flavonoid 3’-Hydroxylase (F3’H) were downregulated, indicating a decrease in the synthesis of proanthocyanidins and other flavonoid metabolites as the leaves transitioned into winter. Additionally, genes involved in the regulation and modification of flavonoids, including UDP-glucose-3-O-glucosyltransferase (UFGT) and Cytochrome P450 enzymes like CYP75B1, which are critical for the hydroxylation and glycosylation of flavonoids, exhibited varying expression patterns, underscoring the complex regulation of flavonoid metabolism. This dynamic expression of flavonoid-related genes highlights the intricate balance between biosynthesis and metabolism that is modulated in response to seasonal changes, offering valuable insights into the optimal harvesting period and potential genetic targets for enhancing flavonoid content in C. paliurus leaves ( Supplementary Tables S6 - S8 ).

To analyze gene-metabolite correlations related to flavonoid biosynthesis, we conducted Weighted Gene Co-expression Network Analysis (WGCNA) using flavonoid-related metabolites and transcriptomic data from C. paliurus leaves collected between August and November. We selected an optimal soft threshold of β=19, ensuring a scale-free network with a fitting index above 0.8 (Langfelder and Horvath, 2008). This threshold allowed us to construct a robust co-expression network ( Figures 4A, B ).

We identified 20 distinct gene modules, each visualized by a unique color in Figure 4C . Highly correlated genes were grouped within the same module, while the gray module contained genes without significant correlation patterns. The lightsteelblue1 module had the most genes (2,479), while the thistle1 module had the fewest (41). Figure 4C illustrates the correlations between modules, with blue shades indicating weaker correlations and red indicating stronger ones. As shown in Figures 4D, E , we examined correlations between gene modules and specific flavonoids detected in the metabolomic analysis, including compounds like naringin, geranolin, kaempferol-3-O-arabinoside, kaempferol-3-O-rhamnoside, and naringin-7-O-glucoside. The tan module demonstrated a high correlation (correlation coefficient > 0.6) with key flavonoids, such as kaempferol-3-O-glucoside, kaempferol-4’-O-glucoside, and kaempferol-3,7-di-O-glucoside. Other modules, including thistle1 and plum2, were associated with naringenin and specific kaempferol derivatives, such as kaempferol-3-O-(2”-p-coumaryl) galactoside and kaempferol-3-O-(6”-p-coumaryl) galactoside (Xie et al., 2015).

To further investigate regulatory genes, we constructed an interaction network of genes within the tan module using Cytoscape, focusing on the top 200 genes with the highest co-expression weights. As visualized in Figure 4F , we identified ten key genes, CpBBX (newGene_23931), CpIP5P3 (newGene_15401), CpLHY (newGene_27003), CpSGPP (newGene_29931), CpCIPK3 (newGene_19518), CpAPF1 (newGene_24020), CpHSFB3 (newGene_10195), CpLOC5396 (newGene_2447), CpMANES (newGene_24731), and CpLOC7593 (newGene_28372), as central regulators within this network (Shannon et al., 2003). Their high connectivity suggests pivotal roles in regulating flavonoid biosynthesis. For example, Cyclocarya_paliurus_newGene_27003 and newGene_29931 were notably active in October, with corresponding upregulation of metabolites like quercetin-3-O-(2”-acetyl) rhamnoside and luteinin-7-O-rutinoside in November. This pattern suggests that these genes may contribute to the synthesis of flavonoids as seasonal conditions change.

Comprehensive differential expression analysis revealed several key genes involved in the flavonoid biosynthesis pathway in C. paliurus, with notable seasonal variations across the studied months. Chalcone Synthase (CHS), a critical enzyme responsible for initiating flavonoid biosynthesis, exhibited peak expression in August, aligning with the highest flavonoid content observed during this period. This suggests that CHS serves as a major regulatory enzyme driving flavonoid production in response to favorable environmental conditions. Similarly, Chalcone Isomerase (CHI), which converts naringenin chalcone to naringenin, displayed elevated expression in August and September, ensuring an adequate supply of naringenin as a precursor for downstream flavonoid biosynthesis (Wan et al., 2024). Flavanone 3-Hydroxylase (F3’H), critical for diversifying flavonoid structures, peaked in September, while Dihydroflavonol 4-Reductase (DFR) showed heightened expression in both August and September, supporting anthocyanin biosynthesis during this period. Flavonol Synthase (FLS), which catalyzes the conversion of dihydroflavonols into flavonols such as quercetin and kaempferol, also exhibited peak expression in August, coinciding with the accumulation of these antioxidant compounds ( Figures 5A, B ) (Jaakola, 2013).

Figure 5A illustrates the flavonoid biosynthetic pathway in C. paliurus, mapping the transformation of precursors like cinnamoyl-CoA and p-coumaroyl-CoA into diverse flavonoids through the activity of key enzymes, including CHS, CHI, F3’H, FLS, Anthocyanidin Synthase (ANS), and Anthocyanidin Reductase (ANR) (Winkel-Shirley, 2001; Hichri et al., 2010). The heat maps integrated into this pathway visually represent the differential expression of these genes, with color gradients highlighting the upregulation and downregulation patterns in response to developmental and environmental stimuli (Falcone Ferreyra et al., 2012). These expression patterns underline the dynamic regulation of flavonoid biosynthesis, driven by the plant’s adaptive responses to oxidative stress, UV radiation, and other environmental challenges, while simultaneously contributing to human health benefits through antioxidant and anti-inflammatory properties (Wang et al., 2024). Genes associated with the biosynthesis of procyanidins, such as Leucoanthocyanidin Reductase (LAR), DFR, ANS, and ANR, were also dynamically expressed. These enzymes play essential roles in synthesizing catechins and epicatechins, as shown in Figure 5B , which emphasizes the significance of flavonol glycosides and procyanidins as key flavonoid classes with diverse biological functions (Xie et al., 2015; Xu et al., 2024). Notably, the heat maps reveal seasonal changes in gene expression, with green and red denoting down- and upregulation, respectively, reflecting the plant capacity to modulate flavonoid synthesis in response to environmental or developmental cues (Hichri et al., 2010).

From August to November, dynamic shifts in the expression of flavonoid-related genes were observed. Genes such as CHS, CHI, and F3’H exhibited prominent upregulation during the transition from August to September, corresponding to an early surge in flavonoid production. In September and October, elevated expression levels of DFR and ANS coincided with increased production of anthocyanins, while FLS activity diversified flavonoid profiles by facilitating the synthesis of flavonols. By late autumn (October to November), genes such as LAR and F3’H were downregulated, marking a decline in proanthocyanidin synthesis as the plant transitioned into dormancy. Additionally, genes involved in flavonoid regulation and modification, including UDP-glucose-3-O-glucosyltransferase (UFGT) and Cytochrome P450 enzymes like CYP75B1, exhibited fluctuating expression patterns, underscoring the complexity of flavonoid metabolism.

The interaction of 57 genes linked to 16 enzymes, as noted in the KEGG pathway (ko00941), highlights the extensive regulatory control underlying flavonoid biosynthesis. These regulatory networks enable C. paliurus to adapt to seasonal changes, bolstering its resilience to environmental stressors while optimizing flavonoid production for both plant defense and human health benefits (Falcone Ferreyra et al., 2012). The findings presented in Figures 5A, B ; Supplementary Tables S6 - S8 not only emphasize the intricate balance between flavonoid biosynthesis and metabolism but also provide valuable insights for optimizing the harvesting period and identifying potential genetic targets for enhancing flavonoid content in C. paliurus.

The total number of flavonoid substances in the leaves of C. paliurus was 144, among which 39 major flavonoid substances with high content accounted for more than 80% of the total flavonoid content from August to November ( Table 1 ). Therefore, we selected these 39 flavonoid substances for the analysis of gene-flavonoid relationships and seasonal variations, Pearson correlation studies were carried out between them and the function genes of the flavonoid synthesis pathway ( Supplementary Table S9 ) and the significantly different transcription factors in DEGs ( Supplementary Table S10 ). The regulation of flavonoid biosynthesis in C. paliurus was thoroughly examined through a comprehensive Pearson correlation analysis, which linked the expression of functional biosynthetic genes and regulatory transcription factors to seasonal variations in flavonoid content. The data from Figures 6 , 7 provide critical insights into how these genes exhibit coordinated fluctuations across different months, significantly influencing the production of key flavonoids such as quercetin and kaempferol glycosides. This seasonal regulation sheds light on the plants intricate molecular responses to environmental stressors and climatic changes ( Supplementary Tables S11 - S12 ).

Beginning with the analysis of functional genes as illustrated in Figure 6 , it is evident that several core biosynthetic genes play pivotal roles in driving the flavonoid biosynthesis pathway. Among these, chalcone synthase (CHS) stands out as the key enzyme catalyzing the initial step in flavonoid biosynthesis, converting malonyl-CoA and p-coumaroyl-CoA into chalcones. The strong positive correlation between CHS expression and flavonoid accumulation, particularly in late summer, suggests that this gene is crucial for the overall regulation of flavonoid biosynthesis in C. paliurus. This observation is consistent with previous studies conducted in white apruce, where it was revealed that accumulation of CHS mRNA in needle tissue following mechanical wounding, or application of signal molecules, such as jasmonic acid or methyl jasmonate (Richard et al., 2000).

Similarly, the role of chalcone isomerase (CHI), which converts chalcones into flavanones, is highlighted by its high expression during the peak flavonoid accumulation period. This gene plays a vital role in facilitating downstream biosynthesis, particularly in the synthesis of flavonols. The expression patterns observed in C. paliurus align with findings in Glycine max, where CHI activity increased under environmental stress, enabling the plant to boost flavonoid production as a means of stress adaptation (Seehaus and Tenhaken, 1998). Furthermore, flavanone 3’-hydroxylase (F3’H), which catalyzes the conversion of flavanones into dihydroflavonols, demonstrated a strong correlation with both quercetin and kaempferol derivatives. This gene is integral to the overall balance of flavonol types produced, and its seasonal expression pattern mirrors those observed in other species such as Malus domestica, where F3’H plays a key role in regulating flavonoid biosynthesis during fruit development (Torres et al., 2020).

Moving further down the biosynthetic pathway, dihydroflavonol 4-reductase (DFR) exhibited significant correlations with anthocyanin content, particularly with cyanidin-3-O-glucoside. This suggests that DFR plays a vital role in anthocyanin biosynthesis during the late summer months, likely as part of the plant’s adaptive response to high light conditions. Such findings are corroborated by studies in Vitis vinifera (grape), where ubiquitin ligase VvWRKY24 could directly interact with the promoters of dihydroflavonol-4-reductase (DFR) to inhibit proanthocyanidins biosynthesis (Zhao et al., 2024). In parallel, flavonol synthase (FLS), which catalyzes the final conversion of dihydroflavonols into flavonols such as quercetin and kaempferol, also demonstrated strong seasonal expression patterns, correlating with the high levels of flavonols observed in late summer. This observation is consistent with research on Arabidopsis, where FLS was shown to be upregulated in response to environmental stimuli, enhancing the synthesis of flavonols for drought inducing (Ma et al., 2024).

Turning to the analysis of transcription factors in Figure 7 , it is apparent that the regulatory machinery governing flavonoid biosynthesis is equally complex and dynamic. Members of the MYB family, including MYB11, MYB12, and MYB111, exhibited strong positive correlations with quercetin and kaempferol glycosides, indicating their central roles in regulating flavonol synthesis. These transcription factors have been widely documented for their involvement in controlling flavonoid biosynthesis, particularly in response to environmental stresses such as UV radiation. In Arabidopsis thaliana, for instance, similar MYB transcription factors have been shown to directly regulate the expression of flavonol biosynthetic genes, enhancing flavonol accumulation under high-light conditions (Stracke et al., 2007).

MYB75 (PAP1) and MYB90 (PAP2) were closely associated with anthocyanin production, particularly cyanidin-3-O-glucoside. These transcription factors are well-known for their roles in regulating anthocyanin biosynthesis in a wide range of plant species. For example, in Vitis vinifera, in ‘Malbec’ grape skins, VviMyb4a and VviMyb4b likely regulate the early genes in the anthocyanin biosynthesis pathway, while VviMybC2-L2 and VviMyb4-like are co-expressed at later stages of the process. Additionally, VviMybC2-L1 and VviMybC2-L3 are activated during two distinct phases of berry development. Their initial induction occurs towards the end of the green stage, suggesting a potential cooperative role in repressing proanthocyanidin synthesis. These genes are then reactivated during maturation, potentially inhibiting anthocyanin production at this later stage (Muñoz et al., 2019). Furthermore, the role of bHLH and WD40 transcription factors, which often act as co-regulators in the MYB-bHLH-WD40 complex, is also highlighted in this study. TT8, a member of the bHLH family, was strongly correlated with proanthocyanidin accumulation, further supporting its function in regulating later stages of the flavonoid biosynthetic pathway, particularly in response to seasonal changes. Studies in Arabidopsis have shown that TT8 plays a critical role in regulating seed coat pigmentation and proanthocyanidin biosynthesis, further validating its importance in flavonoid biosynthesis (Baudry et al., 2004). Another significant finding is the negative correlation between MYB4 and flavonoid content, particularly in late summer. This suggests that MYB4 acts as a repressor of flavonoid biosynthesis, likely through the downregulation of CHS expression to prevent overaccumulation of flavonoids under high-light conditions. This regulatory function has been well-documented in other species, such as Arabidopsis, where MYB4 negatively regulates flavonoid biosynthesis by repressing the expression of early biosynthetic genes (Jin et al., 2000).

Additionally, transcription factors such as WRKY33 and NAC demonstrated significant correlations with flavonoid levels, particularly kaempferol glycosides. WRKY33 is known to play a critical role in regulating plant responses to biotic and abiotic stresses, and its positive correlation with flavonoid levels suggests that it may upregulate flavonoid biosynthesis as part of the plant’s defense mechanisms. Studies in Oryza sativa have similarly highlighted the role of WRKY33 in enhancing flavonoid production under stress conditions, further supporting its function in C. paliurus (Jan et al., 2022). The correlation between EIN3, ABI5, and anthocyanin content suggests that hormonal signaling pathways, particularly ethylene and abscisic acid (ABA), are also involved in the regulation of flavonoid biosynthesis in C. paliurus. The positive correlations between these transcription factors and flavonoid content in late summer and early autumn indicate that they may play a role in modulating the plant transition to dormancy by regulating flavonoid production. Similar findings have been reported in Glycine max, where ethylene and ABA signaling pathways were shown to regulate flavonoid accumulation in response to environmental stresses (Jan et al., 2022).