The P. odoratum rhizomes used for the study were planted for 3 years in a regulated experimental field in Nanling County, Wuhu, Anhui, People’s Republic of China (31°33′N, 118°38′E). Fresh P. odoratum rhizomes were collected on February 2 (winter, abbreviated as WI), May 2 (spring, abbreviated as SP), August 2 (summer, abbreviated as SU), and November 2 (autumn, abbreviated as AU) 2022 ( Figure 1A ). They were identified by Prof. Hongwei Fu (Zhejiang Sci-Tech University, China) and all samples were kept at the College of Life Science and Medicine, Zhejiang Sci-Tech University. The P. odoratum was collected in the third year after planting in the same plot with uniform fertilizer, irrigation and insect spray management. At each sampling, fresh rhizomes from a total of 20 plants were dug from a total of five points (four corners and intersections of the same plot) and processed together. After washing with purified water, the rhizomes were cut into 0.5 cm pieces, flash-frozen in liquid nitrogen and stored at -80°C for subsequent experiments.

Biological samples were dried using a freeze dryer (Scientz-100F) and then ground to a powder using an MM 400 grinder (Retsch, China). Then, 50 mg of the powder was weighed using an MS105DM electronic balance (Sartorius, Germany) and 1200 μL of 70% aqueous methanol with internal standard pre-cooled to -20°C was added. The samples were vortexed six times, centrifuged, and the filtered supernatant used as the sample for UPLC-MS/MS analysis.

Metabolomics analyses were performed by Met-ware Biotechnology Co., Ltd. (Wuhan, China), using peak area integral analysis of all metabolites of different samples by mass spectrometry and integral correction of metabolite peaks of different samples. For two-group analysis, differentially abundant metabolites (DAMs) were determined based on the P-value (P < 0.05) and absolute Log2FC (|Log2FC| ≥ 1.0). For multi-group analysis, DAMs were determined based on the P-value (P < 0.05, ANOVA). Metabolic pathways related to the DAMs were retrieved from the Kyoto Encyclopedia of Genes and Genomes (KEGG) database.

Proteins were extracted from the samples by acetone precipitation. Briefly, samples were ground to a fine powder in liquid nitrogen and homogenized, and 100 mM Tris HCl (pH 7.6), 1 mM phenylmethylsulfonyl fluoride (PMSF), and 2 mM ethylenediaminetetraacetic acid (EDTA) extraction solution were added. Then, the protein samples were boiled for 15 min, lysed by sonication on ice for 10 min, and centrifuged to obtain a clear protein solution. A 4-fold volume of frozen acetone was added to the protein solution, which was allowed to precipitate overnight at -20°C, and centrifuged at 4°C. The resulting precipitate was washed with cold acetone and dissolved in 8 M urea. Finally, the protein concentration was determined using a BCA assay kit according to the manufacturer’s instructions (Ross et al., 2004; Jing et al., 2010).

An equal amount of protein from each sample was used for trypsin digestion. To the supernatant, 8 M urea was added to 200 µl, followed by reduction with 10 mM DTT for 45 min at 37°C and alkylation with 50 mM iodoacetamide (IAM) for 15 min at room temperature in the dark. Then, a 4-fold volume of ice-cold acetone was added and the peptides were allowed to precipitate for 2 h at -20°C. After centrifugation, the protein precipitate was air-dried, re-suspended in 200 µl of 25 mM ammonium bicarbonate solution with 3 µl of 1:50 trypsin (trypsin: protein (m/m), Promega), and digested at 37°C overnight. After digestion, peptides were desalted using a C18 cartridge (IonOpticks, Australia), then dried in a vacuum concentrator, concentrated by vacuum centrifugation, and redissolved in 0.1% (v/v) aqueous formic acid (Jing et al., 2012; Florian et al., 2018).

Samples were separated using a nanoliter flow rate NanoElute UPLC system (Bruker, Germany) and then subjected to mass spectrometry data acquisition using the DDA-PASEF mode of the timsTOF Pro2 mass spectrometer (Bruker, Germany). The library search software used for sample DIA mass spectrometry data was DIA-NNN (v1.8.1), and spectral libraries were used to reanalyze the DIA data for protein quantification. The precursor ions and protein levels of the FDR were filtered at 1%, and the post-filtered data were ready for subsequent biosignature analysis (Yu et al., 2020). For proteomic analysis, three biological replicates were included for each group. Differences with P-values <0.05 were considered statistically significant.

Cluster analysis: K-means analysis was performed on differential metabolites to obtain clustering information of differential proteins based on their correlation with differential metabolites. K-Means analysis was performed using Metware Cloud, a free online platform for data analysis (https://cloud.metware.cn).

Principal component analysis (PCA): the prcomp function in R soft­ware (version R 4.2.0, https://www.R-project.org/) (R, 2022) was used for PCA analysis.

Gene Ontology (GO) enrichment: Gene Ontology enrichment of DEPs was analyzed using the topGO function in R software (version 4.2.0).

qRT-PCR was used to determine the expression of genes related to steroidal saponin and flavonoid synthesis. Three biological replicates were included in each group of samples. Primers were designed using Primer 5.0 software ( Supplementary Table S1 ), and SYBR ^®Green Pro Taq HS qPCR Kit (Accurate Biotechnology Co., Ltd, Hunan, China) was used for qRT-PCR on an ABI7500 fluorescence quantitative PCR instrument (Applied Biosystems, CA, USA). The relative expression levels were calculated using the 2^{− △△Ct} method (Wu et al., 2020).

The results of the principal component analysis showed that metabolic profiles could be clearly distinguished between the different sample groups, with biological samples clustered together ( Figure 1B ), indicating good repeatability and reliability of metabolomics analysis. In this study, a total of 1,460 known metabolites were detected in winter (February), spring (May), summer (August) and autumn (November) rhizome samples using a UPLC-MS/MS system. These included 1,237 differentially abundant metabolites (DAMs, Figure 1C ), which were categorized as amino acids and derivatives, alkaloids, flavonoids, phenolic acids, lipids, steroids, organic acids, lignans and coumarins, nucleotides and derivatives, terpenoids, quinones, and tannins ( Figure 1C ; Supplementary Table S2 ).

Pairwise comparisons revealed that spring-harvested rhizomes showed the highest number of upregulated metabolites, followed by those harvested in winter ( Figure 1D ). Based on the data of differential metabolite content and heatmap based on Z-score normalization, we found that metabolite accumulation showed different trends at different growth stages ( Figure 1E ). The differential metabolites were mainly enriched in amino acid metabolism, specifically valine, leucine and isoleucine degradation, as well as propanoate metabolism ( Figure 1F ).

In order to gain a deeper understanding of the changes in of secondary metabolites at different growth stages, six clusters with different expression trends across the growth stages were generated using the K-means algorithm ( Figure 2A ). Subcluster 1 was the largest and contained 405 metabolites, while subclusters 2-6 contained 186, 243, 179, 145 and 79 metabolites, respectively ( Figure 2B ). Subclusters 1 and 5 show the types and quantities of metabolites upregulated in spring, with the highest accumulation of amino acids and their derivatives. Compared with the other three seasons, alkaloids, phenolic acids, steroids, and organic acids also accumulated more abundantly in spring. Subclusters 2 and 6 show the types and quantities of metabolites upregulated in summer, with the highest accumulation observed for amino acids and their derivatives, followed by nucleotides and their derivatives. Subclusters 3 and 4 show the number and types of metabolites upregulated in winter, with flavonoids accumulating the most, followed by lipids. Subclusters 4 and 6 show the types and quantities of metabolites upregulated in autumn, with lipid accumulation being the highest, followed by flavonoids. Among the metabolites detected in the four seasons, steroids and flavonoids were the main bioactive compounds, and more and more pharmacological effects have been discovered, offering important guidance for the clinical utilization of herbal resources. By combining heatmap and sub-cluster analyses, we focused on steroidal compounds in subcluster 1 and flavonoid compounds in subcluster 3. Subcluster 1 had 38 steroidal compounds, including 12 spirostane saponins, 16 isospirostane saponins, 8 furostane saponins ( Supplementary Table S3 ) and 2 steroids. Subcluster 3 had 56 flavonoids, including 19 other flavonoids (14 homoisoflavanones), 12 flavones, 8 flavanones, 7 flavonols, 5 chalcones, 2 flavanonols, 2 isoflavones and 1 flavanol ( Supplementary Table S4 ). We conducted KEGG functional enrichment analysis on metabolites in clusters 1 and 3, and selected the top 20 based on the P-value ranking. Differential metabolites from cluster 1 were mainly enriched in amino acid biosynthesis and metabolism, ABC transporters, linoleic acid metabolism, nucleotide metabolism, glycerolipid metabolism, and biosynthesis of various plant secondary metabolites ( Supplementary Figure S1A ). Differential metabolites from cluster 3 were mainly enriched in the biosynthesis of flavonoids, isoflavonoid biosynthesis, amino acid biosynthesis and metabolism, plant hormone signal transduction, and biosynthesis of secondary metabolites ( Supplementary Figure S1B ).

In order to elucidate the mechanisms underlying the changes in the content of steroid and flavonoid compounds, we performed proteomic analysis of rhizomes from different growth stages of P. odoratum. By using a protein analysis system, we identified 43,849 different peptides and 7,027 proteins in all 12 samples. Based on the P < 0.05 screening criterion, 3,864 differentially expressed proteins (DEPs) were obtained. To investigate their biological functions, we categorized the identified DEPs into 24 classes using the COG database. The largest category was general function prediction only (425 DEPs), followed by posttranslational modification, protein turnover, chaperones (369 DEPs), translation, ribosomal structure and biogenesis (236 DEPs), energy production and conversion (235 DEPs), carbohydrate transport and metabolism (216 DEPs), signal transduction mechanisms (198 DEPs), intracellular trafficking, secretion and vesicular transport (172 DEPs), as well as secondary metabolite biosynthesis, transport and catabolism (134 DEPs) ( Supplementary Figure S2 ). Classification has shown that many DEPs are involved in the transport of substances and metabolic processes. Proteins with these functions may play an important role in the synthesis and transport of secondary metabolites.

We also analyzed the changes in the expression of DEPs at different growth stages, and categorized the DEPs at different growth stages into six clusters. Based on this, we focused on protein subcluster 5 and protein subcluster 1 ( Figure 3A ), which had the same trend as metabolic subcluster 1 and metabolic subcluster 3. To better analyze the functions of DEPs in subcluster 5 and subcluster 1, we performed GO analysis and selected the top 30 based on the P-value rankings of GO enrichment analysis. In the GO annotations of subcluster 5, the DEPs were assigned to 13 biological processes (mainly protein modification, response to stimuli and defense, metabolic processes, and bioregulation), 3 cellular components (plasma ectodomain, plasma membrane and endoplasmic reticulum membrane) and 14 molecular functions (mainly enzyme activity and binding activity)( Figure 3B ). In the GO annotation of subcluster 1, the DEPs were assigned to 14 biological processes (mainly protein synthesis and hydrolysis, cell cycle, and metabolic processes), 3 cellular components (THO complex, intracellular membrane-bounded organelle, and mitochondrial matrix) and 13 molecular functions (mainly enzyme activities and binding activities) ( Figure 3C ). GO enrichment analysis indicated that the expression of proteins in the rhizomes of P. odoratum at different growth stages may be related to antioxidant activity, binding, stimulatory responses, catalytic activity, metabolic processes, bioregulation, enzyme activity and regulation of molecular functions.

To assess the pathway that are likely involved in changes of secondary metabolites, DEPs from different subclusters were analyzed using the KEGG database. The DEPs in protein subcluster 5were found to be significantly enriched in 16 pathways (P < 0.05), mainly related to plant secondary metabolites, lipid metabolism and amino acid metabolism. Steroid biosynthesis (ko00100) and brassinosteroid biosynthesis (ko00905), which are associated with steroidal saponin synthesis, had significant DEPs ( Figure 4A ). The DEPs in protein subcluster 1 were found to be significantly enriched in 12 pathways (P < 0.05), mainly related to plant secondary metabolites and phytohormone metabolism. Flavonoid biosynthesis (ko00941), isoflavonoid biosynthesis (ko00943), as well as flavone and flavonol biosynthesis (ko00944) ( Figure 4B ), which are related to flavonoid synthesis, were found to have significant DEPs in subcluster 1. These results suggest that proteins from different physiological and metabolic pathways work together to regulate the accumulation of steroidal saponins and flavonoids.

In the steroid biosynthesis pathway, a total of 10 DEPs associated with steroid synthesis were identified. Among them, FDFT1, TM7SF2, CAS1, DHCR7, 3BETAHSDD and TM7SF2 were significantly upregulated in spring, and these proteins may contribute to the accumulation of steroidal saponins ( Figure 5 ; Supplementary Table S5 ). In the flavonoid biosynthesis pathway, a total of 12 DEPs associated with flavonoid synthesis were detected. Among them, HIDH, IF7MAT, PGT1, FLS, CYP75B1 and UGT73C6 were significantly upregulated in winter, and these proteins may promote the accumulation of flavonoid compounds ( Figure 6 ; Supplementary Table S6 ). The relative abundance of differential metabolites related to steroid and flavonoid metabolism differed significantly among groups. The significant differences in these DAMs at different times provide a basis for selecting the appropriate time to extract specific components.

Based on metabolomic and proteomic data, KEGG pathway enrichment was conducted for comprehensive analysis. A bubble plot was drawn using the KEGG pathways enriched in both independent datasets, showing only the top 25 pathways based on P-value ranking ( Supplementary Figure S3 ). The common pathways identified in both the proteome and metabolome included flavonoid biosynthesis, biosynthetic of various plant secondary metabolites, glycyl and dicarboxylate metabolism, alpha linolenic acid metabolism and stilbene biosynthesis. Integrated omics analysis showed that these metabolites and proteins mainly participate in primary and secondary metabolism.

In order to better understand the relationship between metabolomics and proteomics, we conducted Pearson correlation analysis on the biosynthetic pathways of steroidal saponins and flavonoids. Pearson correlation analysis of 38 steroidal saponins from metabolic subcluster 1 and the related DEPs in protein subcluster 5 showed that FDFT1, TM7SF2, DHCR7, CAS1 and 3BETAHSDD were positively correlated with the accumulation of steroidal saponins (P < 0.05), However, SMT2 was negatively correlated with steroidal saponins ( Supplementary Figure S4A ). Pearson correlation analysis of 56 flavonoid compounds in metabolic subcluster 3 and related DEPs in protein subcluster 1 showed that PGT, FLS, CYP75B1, HIDH, IF7MAT and UFT73C6 were positively correlated with flavonoid accumulation (P < 0.05). However, CHS, HCT, CHI were negatively correlated with flavonoids ( Supplementary Figure S4B ).

qRT-PCR was used to analyze the relative expression levels of genes encoding proteins involved in the biosynthesis of steroidal saponins and flavonoids ( Supplementary Figure S5 ). The relative expression levels of the genes encoding farnesyl-diphosphate farnesyltransferase (FDFT1), Δ14-sterol reductase (TM7SF2), cholestenol Δ-isomerase (EBP) and cycloartenol synthase (CAS1) related to steroidal saponin biosynthesis showed a trend of first increasing and then decreasing across the four growth stages. The relative expression levels of the genes encoding flavonol synthase (FLS) and flavonoid 3’-monooxygenase (CYP75B1), which are closely related to flavonoid biosynthesis, showed a trend of first decreasing and then increasing across the four growth stages. While 2-hydroxyisoflavanone dehydratase (HIDH) gradually decreased, chalcone synthase (CHS) and chalcone isomerase (CHI) showed a trend of first increasing and then decreasing. The expression levels of these genes were consistent with the protein levels of corresponding enzymes ( Supplementary Tables S5 , S6 ).