The plant materials, cultivated under natural conditions without exposure to extreme drought, plant diseases, and insect pests, were collected from Guanglong village (N: 27°31′, E: 110°40′; 1320 m above sea level), Shaoyang city, Hunan Province, China. For this experiment, the original flower buds of the wild-type (WT) and xianglei-type (XL) L. macranthoides Hand-Mazz (LM) at seven growth stages were selected between June and July 2022 for sampling( Figure 1A ). Eight biological replicates of flowers from WT and XL (WF and XF) were randomly selected from the upper parts of the branches on a clear day. Three independent biological replicates were used for RNA library construction, another three biological replicates were used for UHPLC−MS/MS analysis, and the remaining samples were stored for further analysis; each biological replicate comprising flowers tissues, randomly collected from each branch, was frozen using liquid nitrogen and stored at -80 °C.

Metabolite extraction and analysis were performed by Shanghai BIOTREE Biological Technology Co., Ltd. (Shanghai China). A 10 mg crushed aliquot of each freeze-dried samples were accurately weighed then transferred to Eppendorf tubes, which containing 500 μL extraction solution (75% methanol in water, precooled at -40°C, containing the internal standard 2-amino-3-(2-chloro-phenyl)-propionic acid). The resulting supernatant was filtered through a 0.22 μm microporous membrane, subsequently diluted 10-fold with the extraction solution containing internal standard, UHPLC separation was conducted via an EXIONLC system (Sciex). Chromatographic separation was achieved by injecting 2 μL sample into an Waters Acquity UPLC HSS T3 column (1.8 μm 2.1x100 mm) with the flow rate set at 400 μL/min, column temperature at 40°C and autosampler temperature at 4°C, respectively. The mobile phases consisted of A and B (A:0.1% formic acid in water, v/v; B:acetonitrile) Separation was conducted with the following gradient: 2% B for 0.5 min,; 2% to 50% B, in 9.5 min,; 50% to 95% B in 1min,; 95% B for 2 min,; 98% to 2% B in 0.1 min,; and 2% B over 15 min,;. The mass spectrometry analysis was performed in MRM mode using Sciex QTrap 6500+(Sciex Technologies), with temperature set at 400°C, and the curtain gas pressure set at 35 psi. The IonSpray voltage was set as 5500 for positive polarity, and 4500 V for negative polarity; The ion source gas pressure was 1:60 psi and 2:60 psi with a DP at 100 V for both positive and negative polarity.

SCIEX Analyst Work Station Software (Version 1.6.3) was employed for all MS data collection and quantitative analysis of target compounds. The peak detection and annotation were performed via an in-house R program and MS/MS database. The relationships among samples was implemented by principal component analysis (PCA) following R package models. Afterwards, the metabolites were annotated and classified to facilitate the analysis of active ingredient biosynthesis via the KEGG compound database. Metabolites with VIP > 1, p values < 0.05 and twofold-change were regarded as significantly differentially accumulated metabolites (DAMs) between groups. The metabolite enrichment analysis and pathways analysis to which the significantly regulated metabolites mapped were integrated using MetaboAnalyst (https://www.metaboanalyst.ca/).

RNA sequencing was performed on the Illumina NovaseqTM 6000 platform, taking the total RNA extracted from frozen petals of wild-type and xianglei-type LMs as samples. After assembly and annotation, the transcripts per million (TPM) values were calculated via SaLmon (version 0.8.2) to quantify gene expression, and DESeq was adopted to determine the differentially expressed genes (DEGs) based on the following parameters: p value<0.05 and absolute |log2(fold-change)| ≥2. To identify the DEGs that were significantly enriched in GO terms and metabolic pathways (p value<0.05), functional analyses was accomplished by performing GO (Gene Ontology) and KEGG enrichment on the annotated DEGs.

The Weighted gene coexpression network construction based on the RNA-seq data via the WGCNA Rpackage (Zhang and Horvath, 2005; Langfelder and Horvath, 2008) was constructed to investigate the genes related to the biosynthesis of flavonoid, phenolic acid and terpenoids during flower development. The resulting adjacency matrix was converted by unsigned TOMtype using the default parameters. Eigengene values were used to calculate the correlation between metabolites and gene modules which then visualized by Cytoscape v3.9.1.

The correlation coefficients among key genes and metabolites were analyzed via the Pearson’s correlation algorithm (Bishara and Hittner, 2012). The integrated genes-metabolite regulatory network was constructed based on KEGG pathway.

Sixteen unigenes including structural gene, hormone-related gene and TF were randomly selected for RT−qPCR evaluation with 18S as endogenous control. The gene-specific RT−qPCR primers were designed via PrimerQuest™ Tool and the sequences were listed in Supplementary Table S12 . RT−qPCR analysis performed on the Bio-Rad CFX96 real-time PCR detection system with three biological replicates. The PCR reaction system and procedure reference to the previous papers (Long et al., 2021) and the 2^-ΔΔCT method (Zhao et al., 2021) was used to calculate the relative expression level of selected gene.

The wild-type and xianglei-type LMs were categorized into seven developmental stages (wild-type, WF1–WF7; xianglei-type, XF1–XF7) ( Figure 1A ) considering previous classification standards for flowering periods, which were based on harvest time, degree of development and use color, size, and appearance as descriptive indicators (Long et al., 2023). The corollas of the two varieties gradually grew at a different rate from the early developmental to the aging stage, the color of the flower buds significantly changed at different stages of development as well. The corolla of wild-type LM gradually expanded from young buds and unfolded in the fifth stage of flowering (WF5), with a color change from green to white, followed by a transformation to yellow. The corolla of XL LM never unfolds from young buds to withering, and the color gradually changes from green to yellow-green, finally turning yellow.

Next, VIP>1, p-value < 0.05 and two-fold change were used to determine differentially accumulated metabolites (DAMs) between the groups ( Figure 2A ). A total of 647 significant DAMs were subsequently grouped into eight clusters via Mfuzz clustering based on their content ( Supplementary Figure S6 ). When the number of DAMs was compared with that of both WF and XF for each flower development stage, venn diagrams were generated, which revealed that 12 flavonoids, 11 phenolic acids and 7 terpenoid compounds accumulated differently between the wild and xianglei variety and different flowering periods ( Figure 2B ). The number of DAMs in the F3, F4, and F5 comparison groups was only half of that in F1, F2, F6, and F7, indicating that the differences in metabolic accumulation between the two varieties of LM were mainly in the early and late stages of flower development, which is consistent with the PCA results. The samples were further divided into the prefloral (F1, F2, and F3), mid-floral (F4 and F5), and late-floral (F6 and F7) development stage for differential analysis and KEGG analyses. KEGG pathway enrichment analyses. revealed that 566 DAMs were annotated. The flavonoid biosynthesis, arginine and proline metabolism, ABC transporter, plant hormone signal transduction, nucleotide metabolism, purine metabolism, and D-amino acid metabolism pathways were significantly enriched ( Figures 2C–H ).

After removing low-quality reads, ambiguous reads, and adapter sequences, we initially generated 259.93 Gb of data from 42 libraries, and each library yielded between 33.05 and 57.26 million clean reads ( Supplementary Table S3 ). The Q20%, Q30%, and GC% results indicated that the throughput and sequencing quality were sufficiently high to warrant further analysis ( Supplementary Table S4 ). After de novo assembly, a total of 107,737 unigenes with an N50 of 1602 bp were detected in the flower samples ( Supplementary Table S4 ). In addition, the assembled 107,737 sequences were all annotated in the database, and 51757, 33397, 45457, 45742, 56402, and 47389 unigenes were annotated respectively ( Supplementary Table S5 ).

Through differential expression analysis between the two L. macranthoidesis varieties at different flowering periods, a total of 34976 DEGs were identified between the sample groups, with the number of DEGs ranged from 137 to 17985 ( Figure 3A ). There were 14,507 DEGs in the seven WF comparison groups, 16,014 DEGs in XF, and 32,075 DEGs in XF vs. WF. Compared to WF, there were more DEGs in XF, and most of them showed an upward trend, which was reflected in each comparison group. Furthermore, among the WF, XF, and XF vs WF groups, 5,776 common DEGs were indicated using Venn diagrams ( Figure 3B ). These common 5,776 DEGs exhibited varied expression patterns, showing a closer relationship( Figures 3C, D ). Among the WFs, the most abundant DEGs were found in the WF6 vs. WF7 comparison group, whereas those in the XF5 vs XF6 comparison group were in the XF. However, the number of DEGs in the early bud (F1 and F2) comparison between the two varieties was significantly lower, suggesting that the two varieties differed significantly in terms of middle and late flower development.

KEGG-based enrichment analysis of all the DEGs revealed the top 20 KEGG pathways which were significantly enriched ( Figure 3E ). These pathways included metabolic pathways (ko01100), biosynthesis of secondary metabolites (ko01110), flavonoid biosynthesis (ko00941), terpenoid backbone biosynthesis (ko00900), sesquiterpenoid and triterpenoid biosynthesis (ko00909), phenylpropanoid biosynthesis (ko00940), etc. GO enrichment analysis of these DEGs revealed that they were widely distributed in three functional groups, i.e., biological processes, molecular functions, and cellular components ( Figure 3F ). According to the GO analysis, the response to acid chemical was the largest group of biological processes. The genes related to “cell periphery”, “cell wall”, and “external encapsulating structure” were predominant in the category of “cellular components”. In the molecular function category, the genes were mostly involved in oxidoreductase activity, catalytic activity, transmembrane transporter activity, and transporter activity.

To reveal the molecular mechanisms underlying the different accumulations of active ingredients in LM at different stages, we conducted transcriptome sequencing and correlated the results with the metabolomic data. We identified 128 enriched structural genes that regulate 34 enzymes (CHS, CHI, FNS, F3’,5’H, FLS, F3H, DFR, ANS, ANR, PAL, C4H, 4CL, CCR, CAD, HCT, C3’H, F5H, COMT, CCoAOMT, ACAT, HMGCS, HMGCR, MVK, MVD, IDI, DXS, DXR, MCT, CMK, MDS, HDS, HDR, TPS, FPS, SS, and SE), of which 16 genes regulate 9 enzymes involved in flavonoid biosynthesis, anthocyanin biosynthesis, flavonoid and flavonol biosynthesis, and isoflavonoid biosynthesis; 46 genes regulate 8 enzymes involved in phenolic acid biosynthesis; and 66 genes regulate 17 enzymes involved in the biosynthesis of terpenoid backbone biosynthesis and sesquiterpenoid and triterpenoid biosynthesis ( Supplementary Table S6 ). To identify the genes closely related to flavonoid, phenolic acid and terpenoid biosynthesis, all 128 enriched structural genes (DEGs) and differentially accumulated metabolites (DAMs), including 165 flavonoids, 34 phenolic acids, and 134 terpenoids detected at different flower developmental stages, were selected for integrative analysis. Using Pearson’s correlation analysis, we evaluated the correlation coefficients (r ≥ 0.8) between the DEGs and DAMs during flower development ( Figure 4 ). In flavonoid biosynthesis, 2 CHIs (TRINITY_DN46726_c4_g1, TRINITY_DN43583_c0_g1), ANSs (TRINITY_DN28844_c0_g1), 4 CHSs (TRINITY_DN46824_c0_g4, TRINITY_DN46824_c0_g2, TRINITY_DN46824_c0_g5, TRINITY_DN46824_c0_g6), DFRs (TRINITY_DN48211_c1_g5), and F3Hs (TRINITY_DN41729_c3_g1) were positively correlated with hesperetin, chalconaringenin, naringenin, amentoflavone, kaempferol, safflor yellow A, and butin. Nine candidate DEGs, including PAL (TRINITY_DN41142_c2_g1), 4CL (TRINITY_DN44156_c3_g1), C3’H (TRINITY_DN50056_c3_g3), 3 CCR (TRINITY_DN48406_c0_g4, TRINITY_DN53148_c0_g2, TRINITY_DN53072_c0_g4), CAD (TRINITY_DN42601_c1_g4), COMT (TRINITY_DN54977_c4_g2), and HCT (TRINITY_DN44906_c4_g2), were identified, which are positively connected by 5 phenolic acid compounds (isoferulic acid, methyl caffeate acid, dihydroconiferyl alcohol, and syringarin). Furthermore, 3 ACAT, DXR, HDS, 3 HMGCS, MDS, MVD, 2 MVK, SS, and 5 TPS genes were identified as key genes correlated with 14 terpenoids involved in terpenoid biosynthesis. Among them, 30 pairs were positively correlated, and 2 pairs were negatively correlated. To further explore the mechanism of flavonoid, phenolic acid, and terpenoid accumulation in LMs, markedly upregulated or downregulated DAMs and key structural genes were mapped onto the flavonoid, phenolic acid, and terpenoid metabolic pathways, which are depicted in Figures 4A–C , respectively. The expression levels of the genes encoding enzymes that catalyze each step of the biosynthetic pathway are shown in Supplementary Figure S4 .

Hormone-related genes, including 105 IAA-related genes (auxin), 18 CTK-related genes (cytokinin), 10 GA-related genes (gibberellin), 16 ABA-related genes (abscisic acid), 19 ETH-related genes (ethylene), 15 BR-related genes (brassinosteroid), 7 JA-related genes (jasmonic acid), and 10 SA-related genes (salicylic acid), were further analyzed to explore the role of plant hormones in the biosynthesis of active ingredients at different developmental stages of L. macranthoides,. The results revealed 77 genes exhibiting significantly differentially expression between the two breeds, with IAA-related genes accounting for 59.7%. The heatmap ( Supplementary Figure S5 ), cluster profiles ( Figure 5A ) and bar chart ( Figure 5B ) based on relative expression show the distributions of these hormone-related genes in eight different expression trend gene sets. Subsequent trend analysis revealed that a total of the 200 genes were distributed across 8 clusters. Clusters 1 (19 genes) and 7 (33 genes) were highly expressed only in WF4 and XF1, respectively; Clusters 2 (20 genes) and 8 (20 genes) were highly expressed mainly in WF4, 5, 6, 7, and the middle and late flower stages in WF. Conversely, 17 genes in Cluster 3 were highly expressed mainly in XF4, 5, 6, 7, and the middle and late flower stages in XF; the expression of 24 genes in Cluster 4 was significantly lower in XF than in WF; notably, the expression of 31 genes in Cluster 6 showed the same trend in both WF and XF, which was highly expressed in early flowering stages F1, 2, and 3 but decreased in the middle and late stages; 27 genes in Cluster 5 presented similar trends in WF. The enrichment results indicate that IAA-related genes (ARF/AUX/GH3/IAA/SAUR/LAX) were enriched mostly in Clusters 1, 2 and 7, which play a role mainly in the late stage of flower development, with some playing a specific role in XF1; CTK-related genes (AHP/AHK/ARR) were enriched mostly in Clusters 5 and 6, indicating potential functions in the early stage of flower development; 11 ETH-related genes (CTR1/EIN2/ERS1/ETR/ERF/EIL/MKK) were enriched in Clusters 3 and 8, showing high expression in the middle and late stages (F4, 5, 6, 7) of WF and XF, closely related to their function in flower maturation and aging processes; and GA-related genes (GID/GAI), ABA-related genes (ABF/PP2C/PYL), BR-related genes (BAK/CYCD/BRL), JA-related genes (MYC2/JAZ/COI1), and SA-related genes (NPR/TGA/TGAL) presented expression trends in Clusters 6, 4, 7, 4, and 5, respectively. Cluster analysis revealed a close connection between these genes and LM flower development. However, these findings reflect the direct or indirect regulation of hormone-related genes and the accumulation of active substances during flower development in LM.

For coexpression network analysis, the WGCNA package was used. After filtering, the top 50% of genes with minimal MADs, outlier genes and samples were removed. We obtained a total of 18,082 genes to construct the module. Modules with a distance of less than 0.5 were merged, and a total of 15 modules were generated ( Figure 6 ). The gray module represent a gene set that cannot be attributed to any module. Five coexpression modules, including M1, M2, M4, M5, and M10, were specifically correlated with individual flower stages. Among them, M4 and M5 were strongly correlated with WF6 (r = 0.99) and XF6 (r = 1), respectively, indicating that the gene sets in these modules are specifically expressed in their corresponding stages. In addition, analysis revealed that more than 80.5% of the key structural genes were duplicated in the WGCNA modules. Subsequently, Pearson correlation analysis revealed that modules M9, M10, and M14 were co-expressed with 20 abundant metabolites related to flavonoid biosynthesis; modules M3, M9, and M12 were coexpressed with 20 abundant metabolites related to phenolic acid biosynthesis. Modules M3, M6, and M9 were co-expressed with 20 abundant metabolites related to terpenoid biosynthesis. These findings indicated that genes related to the regulation of key metabolite accumulation in these modules had the strongest correlation. Furthermore, we mapped a regulatory network of TFs and key structural genes, revealing that the expression of key structural genes (such as CHS, CHI, F3H, DFR, ANS, PAL, 4CL, COMT, ACAT, HMGCR, MVK, DXR, and TPS) was closely associated with TF families (such as MYBs, ARFs, bHLH, GRASs, HD-ZIPs, and B3) ( Supplementary Figure S6 ).

To confirm the RNA-seq results, RT−qPCR was performed to evaluate the expression of 16 genes involved in the biosynthetic pathway. Our observations revealed that MYB (TRINITY_DN43595_c7_g1), bZIP (TRINITY_DN56880_c3_g1), and bHLH (TRINITY_DN44673_c3_g1) TFs were significantly upregulated in WFs, and exhibited the lowest expression levels in XF ( Figure 7 ). The RT-qPCR results of the structural genes significantly correlated with the RNA-seq results. The expression levels of CHS (TRINITY_DN46824_c0_g4), HMGCR (TRINITY_DN46117_c1_g1), DFR (TRINITY_DN48211_c1_g5) and F3H (TRINITY_DN41729_c3_g1) in XF were 2 to 3 fold higher than those in WF, indicating that these genes are potentially involved in the regulation of biosynthetic pathways. Additionally, we detected the expression levels of the JA-related gene MYC (TRINITY_DN49103_c1_g1) and the IAA-related gene GH3 (TRINITY_DN51174_c1_g2) in WF and XF, which was accurate. These findings indicated that the relative expression levels were generally in line with the RNA-Seq data, confirming the accuracy and reliability of the transcriptomic profiling results.

Flavonoids are widely present in most green plants and play essential roles in plant growth, development, and stress response processes. For a long time, flavonoids have long been considered indispensable sources of nutrients and medicinal components for the human body (Roy et al., 2022; Yao et al., 2004). In this study, we analyzed 165 flavonoids among the seven stages of flower development in wild-type and Xianglei-type LMs, demonstrated that more than half of the compounds showed greater accumulation in XF than in WF. Most anthocyanins, including delphinidin-3-O-glucoside, pelargonidin-3,5-O-diglucoside chloride, malvidin 3-O-glucoside, malvidin 3-glucoside, delphinidin, cyanidin, procyanidin B1, and procyanidin B2, were highly abundant in WF4 and XF4. The variation and high accumulation of anthocyanins in the F4 stage may crucially regulate the development of the yellow color in LMs, resulting in phenotype changes in color during the development of floral organs. Plant flavonoids biosynthesis begins with the phenylpropane biosynthesis pathway, which is regulated by multiple key enzymes (Liu et al., 2021). Naringenin is a key precursor substance for the biosynthesis of many flavonoid components. Under the catalysis of CHS and CHI, p-coumaroyl-CoA produces flavonoids, including naringenin, quercetin, and chalconingenin, whose levels decrease consistently. In the pathway, the ionic abundances of the terminal compounds in the pathway, including rutin, luteolin, cynaroside, and tricetin initially increased and then peaked at F4 and F5 before decreasing, mirroring the trend of their upstream counterparts. This may be induced by the reduced accumulation of secondary metabolites upstream of the pathway, which redirects metabolic ion flow and consequently increases the accumulation of downstream metabolites. Correlation analysis revealed that key rate-limiting enzyme-encoding genes, such as PAL, CHS, CHI, FNS, ANS and ANR, are involved in the metabolism of flavones and CQA in LM. We detected significant differences in the expression of structural genes, which are similar to those involved in the metabolic accumulation of flavonoids. Among them, F3’,5’H (TRINITY_DN50067_c0_g1) was significantly elevated in the WF4, WF5, and WF6 stages but was stable at relatively low levels in XF, which reflected consistency with the accumulation patterns of rutin, luteolin, cynaroside, and tricetin. ANS (TRINITY_DN28844_c0_g1) and ANR (TRINITY_DN37089_c0_g1) was more highly expressed in the prophase stages of F1, F2, and F3, which was accompanied by an increased abundance of anthocyanins at F4 and F5. These findings indicate that these genes play important roles in regulating the accumulation of secondary metabolites during different developmental periods. However, their mechanisms of action and regulatory mechanisms require further investigation. Additionally, our network analysis revealed that the TFs that presumably regulated flavonoids in LMs, among which 25 bHLHs and 24 MYBs were enriched. Numerous studies have shown that the MYB and bHLH transcription factors can act as transcriptional repressors or activators that regulate plant growth, development, and stress response processes by mediating of flavonoid biosynthesis. For instance, the overexpression of CmMYB8 in Chrysanthemum morifolium significantly inhibits the expression of genes such as CHS, CHI, F3H, F3’H, and DFR, and significantly downregulates components such as rutin, isorhamnetin, quercetin, and kaempferol (Wang et al., 2021b; Zhu et al., 2020). After the transformation of Erigeron breviscapus EbbHLH80 into tobacco (N. tabacum), the total flavonoids content in tobacco was significantly increased, resulting in the differential enrichment of 98 flavonoid components, and significant upregulation of gene expression including CHS, CHI, F3H, FLS, ANS, and DFR (Gao et al., 2023). MYBs typically form complexes with bHLH and WD40 proteins to regulate gene transcription. The expression of key structural genes in rice (O. sativa) is effectively activated by MBW members, impacting the anthocyanin content (Sun et al., 2022). These results indicate variable TF associated complex regulatory mechanisms underlying flavonoid accumulation, and the specific mechanism. Further research can validate the synergistic regulations of different active components in LM.

Phenolic acids are the main active ingredient of LM with multiple therapeutic and pharmacological activities such as cardioprotective, antioxidant, antiviral, antibacterial, anti-inflammatory, and antitumor activities, which are closely associated with the traditional effects of clearing away heat and detoxification (Chen et al., 2024). Our previous findings suggest that the content of chlorogenic acid accounts for over 62% of the total phenolic acid content in LM flowers and that the contents of chlorogenic acid, isochlorogenic acid A, and isochlorogenic acid C comprise over 90% of the total phenolic acid content (Hu et al., 2016). In the Chinese Pharmacopoeia 2020, it is considered as an important quality control indicator for the medicinal herb Lonicerae Flos (Chinese Pharmacopoeia Commission, C.P, 2020). The phenolic acid components in LM are mostly isomers of chlorogenic acid; the primary metabolites of chlorogenic acid include caffeic acid, quininic acid and ferulic acids; the products of further metabolism are mostly small molecules, such as catechin and hydroxybenzoic acid derivatives. We identified 43 phenolic acid compounds in the two varieties of LM. During flower development, the levels of most phenolic acid compounds, including 5-O-caffeoylshikimic acid, caffeic acid, quinic acid, and ferulic acid, as well as methylcinnamate, peaked in F4. The content of isochlorogenic acid B was the highest in F1 and then decreased. Chlorogenic acid gradually increased from F1, reaching its highest level at F4, and then gradually decreased. The isochlorogenic acid A content changed irregularly and then gradually decreased. The variation trend of phenolic acid content during the development of LM flowers varied depending on the types of phenolic acid components, which may be related to their mutual transformation. In plants, phenolic acids are primarily obtained via the phenylpropanoid biosynthetic pathway. PAL, C4H, and 4CL enzymes contribute to the synthesis of various phenolic acids (Białecka-Florjańczyk et al., 2018). Different phenolic acid derivatives are synthesized by the CAD, CCR, and COMT enzymes, and HCT acts as a key rate-limiting enzyme in CGA metabolism, suggesting that phenolic acid biosynthesis genes are closely associated with phenolic acid accumulation. In this study, all 46 key genes (PAL, C4H, 4CL, HCT, CAD, CCR, COMT, F5H, and C3’H) and TFs (r > 0.80) were found to be differentially expressed at seven developmental stages in WT and XL. These TFs, including 9 NACs, 6 MYBs, 6 bHLHs, 3 ERFs, 6 WRKYs, 10 C2H2s and other family members, have been reported to regulate the biosynthesis of phenolic acids such as chlorogenic acid and salvianolic acid in response to hormone signaling (Liu et al., 2022; Cao et al., 2021; Shi et al., 2022; Guan et al., 2023). According to the correlation network analysis, the key rate-limiting enzymes LmPAL1 (TRINITY_DN41142_c2_g1) and LmPAL2 (TRINITY_DN52595_c1_g2) were positively regulated by LmMYB37 (TRINITY_DN55478_c3_g1) and LmWRKY12 (TRINITY_DN47707_c4_g4); LmbHLH57 (TRINITY_DN56449_c0_g1), LmNAC12 (TRINITY_DN48409_c3_g8), LmNAC16 (TRINITY_DN49974_c1_g3) and LmWRKY13 (TRINITY_DN48226_c2_g2) negatively regulated the Lm4CL1 (TRINITY_DN44156_c3_g1), LmCCR6 (TRINITY_DN53148_c0_g2) and LmCOMT2 (TRINITY_DN54977_c4_g2) in LM. These TFs are hypothesized to regulate the expression of structural genes through activation or inhibition; hence, they potentially participate in the biosynthesis and accumulation of phenolic acids in LM.

Previous research indicates that triterpenes and their glycosides are characteristic compounds of LM. More than 50 triterpenoid saponins and their aglycones have been identified in plants belonging to the Lonicera genus (Fang et al., 2020). Structural analysis identified dipsacoside B, macranthoidin B, macranthoidin A, Hederagenin, cauloside C, Akebia saponin D, and α-Hederin as oleanane-type saponins (Liu et al., 2013). The triterpenoid saponins in LM exhibit significant anticancer, anti-solid and nonsolid tumor, antiviral, free radical scavenging, antioxidant, and liver protection activities both in vitro and in vivo and thus have substantial developmental potential (Fan et al., 2018; Chen et al., 2023, Chen et al., 2020). In this study, 134 terpenoids were tentatively identified by metabolite profiling using UHPLC-MS/MS. Our results reflected the rapid accumulation of most triterpenoids in LM at stages F1, F2, and F3, with the peak in stage F4, followed by a reduced level. Terpenoids, comprising more than 556 odor compounds, constitute the most diverse class of plant floral volatiles (Dötterl and Gershenzon, 2023). The main types of plant floral compounds are monoterpenes and sesquiterpenes, which play important roles in the presentation of different plant aromas, such as in herbaceous plants including Lilium sp (Yang et al., 2022), Hedychium (Ke et al., 2021), and woody plants including Rosa damescena (Shi and Zhang, 2022), Osmanthus fragrans L (Sheng et al., 2021), Jasminum sambac (L.) Ait (Wang et al., 2022d). Monoterpenoid and sesquiterpenoid components with differential accumulation at different flowering stages were identified in the two varieties, mainly in the mid-to-late stages of flower growth characterized by flowers movements and opening. We speculated that this accumulation pattern is the possibly associated with variable aroma during development. Although triterpenoid saponins in different plants share a common MVA and MEP skeleton synthesis pathways, the post-modification processes from 2,3-oxidized squalene to the triterpenoid saponin formation involve different patterns, resulting in multiple complex structures. Multiple key enzymes derived from the MEP and MVA pathways, which included synthases, intermediates, and transferases enzymes which involved in the synthesis and metabolism pathways of terpenoids were identified using RNA-seq data. These genes exhibit distinct tissue-specific expression patterns during the flower development cycle. TPS genes play crucial roles in plant terpene metabolism by modifying the structural types of various terpenes. The catalytic formation of terpenoids by TPS also requires further modifications by glycosyltransferases, cytochrome P450 monooxygenases, methyltransferases, acetyltransferases and short-chain dehydrogenases to produce a diverse variety of terpenoids (Cheng et al., 2020; Sawai and Saito, 2011). Most terpenoid-related CYPs are considered members of the CYP71 subfamily and many are enzymes that catalyze terpenoid oxidation. UGT 71, 73, 74, 85, and 91 clans have been reported participate in saponin biosynthesis (Rahimi et al., 2019), among which the UGT73 clan is considered as candidate gene crucial for the biosynthesis of oleanane-type triterpenoid saponin owing to the catalytic activity for the production of α - hederin from cauloside A (Yin et al., 2023). This function forms the evolutionary basis for the diversity of terpene skeletons. In this study, 15 TPS genes, 9 CYP450s and 2 UGTs genes involved in skeleton modification were identified in the WGCNA gene module. These genes constitute the core gene set and play an irreplaceable role in the LM terpenoid biosynthesis and metabolic pathways. Several studies have confirmed that members of multiple TF families, including the MYB, bHLH, WRKY, ERF/AP2, bZIP, and ARF TF family proteins, participate as regulatory factors in terpenoid biosynthesis (Huang et al., 2019; Xu et al., 2022; Man et al., 2023). In some species, MYC2 TFs belonging to the bHLH family positively regulate the gene expressions by directly binding to G-box-like motifs in the TPS and CYP promoters. In Panax notoginseng, ginsenoside biosynthesis is regulated by TFs, including WRKY, AP2, DELLA, ERF3, and MYC2, through the JA pathway via PgLOX6 (Rahimi et al., 2016). In this study, several TFs, including 12 bHLHs, 8 MYBs, 4 AP2/ERF family members and others, which are co-expressed with terpenoid pathway genes and play important regulatory roles, were overexpressed in the WGCNA module, suggesting that these candidate TFs may participate in the biosynthesis of terpenoids in LM.

Previously, several active ingredients including flavoids, terpenoids, alkaloids, saponins, and phenolic acids, were isolated and identified from medicinal plants owing to their crucial roles in the development of drugs and the promotion of human health (Wang et al., 2023). Numerous reports reflected that the transcriptional regulation of the biosynthesis of active ingredients in medicinal plants involves complex and multilayered mechanisms, enhancing its importance in molecular pharmacology research (Zheng et al., 2023; Shi et al., 2024). A complex network of upstream transcriptional regulators comprises different plant hormones that are considered key signals for plant growth and developmental processes, including root elongation, tillering, flowering, and fruit ripening (Mazzoni-Putman et al., 2021). As macroscopic upstream regulatory factors, plant hormones usually indirectly regulate the transcription of transcription factors such as MYBs, bHLHs, and ERFs through specific signaling pathways (Tong et al., 2024; Wu et al., 2021). Transcription factors regulate the activation or inhibition of downstream target genes, thereby regulating the expression of multiple downstream target genes with multiple functions. These functional TFs belong to multiple families and are induced by various hormones, which have been shown to regulate the biosynthesis of active ingredients. For instance, both auxin and ethylene can increase the transcription of AtCHS, AtCHI, and AtFLS during seed development in Arabidopsis, thereby increasing the anthocyanins and flavonols contents. However, AtF3’H is significantly increased by auxin only (Lewis et al., 2011). Extensive research has shown that JA plays an important role in regulation of saponin accumulation. The TF PgbHLH28 promotes saponin content by inducing the expression of PgHMGR2 and PgDXS2, and mediates JA signaling in Platycodon grandiflorus (Zhang et al., 2024). The expression of SmMYB111 is significantly downregulated by MeJA, GA, or SA treatment (Li et al., 2018). Many hormone-related genes involved in the biosynthesis of LM active ingredients were identified using transcriptional data analysis. KEGG analysis also revealed that plant hormone signal transduction pathways were significantly enriched in multiple comparison groups, suggesting that plant signal transduction potentially plays an important role in the flower development of L. macranthoides, simultaneously affecting the accumulation of active ingredients. Through a combined analysis of transcription and metabolism, we further identified hormone-related genes and key transcription factors involved in the regulatory network associated with active ingredients, thus forming a bridge between hormone signaling TF structural genes. A network diagram of several hormones involved in the transcription factor regulation of downstream active component structural genes is shown in Figure 8 .