The plant materials of S. plebeia were obtained from the Medicinal Botanic Garden of China Pharmaceutical University, Nanjing, China (118.83E, 31.95N). The samples were collected at four representative stages—the basal leaf stage (stage 1, harvest parts: leaves (L1) and root (R1)); stem elongation stage (stage 2, harvest parts: leaves (L2), stem (S2), and root (R2)); flower stage (stage 3, harvest parts: leaves (L3), stem (S3), flower (F3), and root (R3)); and ripening stage (stage 4, harvest parts: leaves (L4), stem (S4), and root (R4)) ( Figure 1 ). All organs were chopped and snap frozen in liquid nitrogen, then stored at -80°C. Eight biological replicates were sampled for metabolomics assay. Three biological replicas were sampled for transcriptomic assay.

The reference standards, 6-hydroxyluteolin 7-O-β-glucoside, nepitrin, luteolin 7-O-β-glucoside, homoplantaginin, apigenin 7-O-β-glucoside, luteolin, nepetin, 6-methoxynaringenin, apigenin, hispidulin, phenylalanine, cinnamic acid, 4-coumaric acid, tyrosine, naringenin chalcone, naringenin, scutellarein, UDP-glucose, and NADPH had >98% purity and were commercially available (Push, Chengdu, China).

All samples were crushed, weighed 0.5 g into 10 mL of 70% methanol (0.02 mg/mL baicalin and 0.01 mg/mL irigenin as internal standards for UPLC-Q-TOF-MS/MS assay). The mixed solution was sonicated for 30 min. After ultrasonic extraction, the extraction solution centrifuged at 10,000 rpm for 30 min. Then, the supernatant was filtered by a 0.22µm membrane filter. The solutions were injection into the column for UPLC-Q-TOF-MS/MS and HPLC-DAD analysis.

The UPLC-Q-TOF-MS/MS analysis was using an AB SCIEX TripleToF® 5600 mass analyzer (Redwood City, CA, USA), in negative ion modes. A C₁₈ reversed phase column (50 mm × 2.1 mm, 1.5 µm, Thermo Scientific, USA) was used for UPLC analysis with 40°C column temperature. The gradient contains solvents A (0.1% formic acid in water) and solvents B (acetonitrile) with the gradient as follows: 0 min, 12% B; 3 min,16% B; 5 min, 16.5% B; 8 min, 18% B; 9 min, 27% B; 12 min, 33% B; 13 min, 50% B; and 14 min, 95% B. The flow rate was maintained at 0.4 mL/min. MS survey scan of 100-2000 Da; ion source heater, 550°C; ion spray voltage, 4500 V; and collision energy, 44 V. The Peakview Software (version 1.2.0.3, AB SCIEX, Redwood City, CA, USA) were used to analyzed MS/MS data.

The contents of 10 main flavonoids (6-hydroxyluteolin 7-O-β-glucoside, nepitrin, luteolin 7-O-β-glucoside, apigenin 7-O-β-glucoside, homoplantaginin, luteolin, nepetin, 6-methoxynaringenin, apigenin, and hispidulin) in each sample of S. plebeia were performed by external standard method, using a 1290 HPLC instrument (Agilent Technologies, Cambridge, USA) with DAD wavelength at 280 nm and 342 nm. The HPLC condition were consistent with those for the UPLC-Q-TOF-MS/MS analysis. Three replicates per sample.

Total RNA was isolated from S. plebeia using the RNA simple total RNA kit (TIANGEN, Beijing, China). A cDNA library was then constructed and sequenced on an Illumina NovaSeq 6000 platform (Illumina, San Diego, CA, USA). All raw reads were uploaded to NCBI at SRA database (PRJNA952802) and include 21 accession items (SRR24128449−SRR24128469). FPKM (the number of fragments per kilobase transcript per million fragment map) was used to calculate the gene expression level. Genes with |log2FC|≥1 and Benjamini−Hochberg-adjusted P-values<0.05 were determined as differentially expressed genes (DEGs). Using the following databases for gene functional annotation: NCBI non-redundant protein sequences (Nr), Kyoto Encyclopedia of Genes and Genomes (KEGG), NCBI non-redundant nucleotide sequences (Nt), Gene Ontology (GO), a manually annotated and reviewed protein sequence database (Swiss-Prot), a manually annotated and reviewed protein sequence database (Pfam), and Clusters of Orthologous Groups of Proteins (KOG/COG).

Based on the results of gene functional annotation, the candidate key genes PAL, C4H, 4CL, CHS, CHI, F6OMT, and UGT were selected in the homoplantaginin biosynthesis pathway. The sequences of FNS and F6H were obtained by a local BLAST algorithm-based search and relevant literature (Chen et al., 2014; Zhao et al., 2018). Phylogenetic tree construction with MEGA11 (https://megasoftware.net/). The Neighbor-Joining method was used to construct this tree with bootstrap (n=1000).

The open reading frames (ORFs) of PAL, 4CL, CHS, CHI, and F6OMT; C4H, FNS, and F6H; and UGT were individually cloned into the prokaryotic expression vectors pET28a, pET32a, and pSJ8-MBP, respectively. The primers used for construction are shown in Supplementary Material Table S1 . Then, the recombinant plasmids were transferred into E. coli BL21(DE3), respectively. The transformants were grown in 100 mL Luria-Bertani (LB) medium containing antibiotics, which was cultured at 37°C. To express of target proteins, 0.5 mM IPTG (isopropyl 1-β-D-thiogalactoside) was added to the medium at OD₆₀₀ 0.6–0.8 of the strains, and then further incubated at 16°C or 25°C for 16 h. The cells were collected by centrifugation, then re-suspended with phosphate-buffered saline, and disrupted using a sonicator on an ice-bath. The crude lysate was centrifuged (12000 rpm, 30 min) at 4°C, removed from the cell fragments. The supernatant was purified using Ni-NTA columns (Smart-Lifesciences, Changzhou, China) and further purified using protein purification system SDL-030-F2 (SePure Instruments Co., Ltd, Suzhou, China). The concentrations of protein were determined with the Bradford method. The purity of purified protein was detected by SDS-PAGE

The reaction volume was 200 µL. SpPAL reaction components included 100 mM Tris-HCl, 3 mM phenylalanine, and purified SpPAL enzyme, which were incubated at 37°C for 8 h. SpC4H reaction components included 100 mM Tris-HCl, 1 mM NADPH, 0.5 mM GSH (reduced glutathione), 50 μM cinnamic acid, and purified SpC4H enzyme, which were incubated at 28°C for 12 h. Sp4CL reaction components included 100 mM Tris-HCl, 5 mM ATP, 0.3 mM p-Coumaroyl-CoA, 5 mM MgCl₂, 50 μM cumaric acid, and purified Sp4CL enzyme, which were incubated at 37°C for 6 h. SpCHS reaction components included 100 mM Tris-HCl, 5 mM ATP, 6 mM p-Coumaroyl-CoA, 5 mM MgCl₂, 3 mM tyrosine, and purified SpCHS enzyme, which were incubated at 30°C for 12 h. SpCHI reaction components included 100 mM PBS, 50 μM naringin chalcone, and purified SpCHI enzyme, which were incubated at 25°C for 2 s. SpFNS reaction components included 100 mM Tris-HCl, 50 μM naringenin, 1 mM NADPH, 0.5 mM GSH, and purified SpFNS enzyme, which were incubated at 28°C for 12 h. SpF6H reaction components included 100 mM Tris-HCl, 50 μM apigenin, 0.5 mM GSH, 1 mM NADPH, and purified SpF6H enzyme, incubated at 28°C for 12 h. SpF6OMT reaction components included 100 mM Tris-HCl, 50 μM scutellarein, 1 mM SAM (S-adenosyl-L-methionine), 1 mM DTT (dithiothreitol), 1 mM MgCl₂, and purified SpF6OMT enzyme, which were incubated at 28°C for 12 h. SpUGT reaction components included 50 mM PBS, 100 µM hispidulin, 1 mM UDPG (uridine diphosphate glucose), and purified SpUGT enzyme, which were incubation at 30°C for 4 h. Enzyme activity was analyzed via HPLC and UPLC-Q-TOF-MS/MS.

Total RNA from leaves, roots, stems, and flowers of S. plebeia was extracted as the protocol described above. Reverse transcription was done by using HiScript II Q RT SuperMix for qPCR (Code: R223-1, Vazyme, Nanjing, China). The ChemQ SYBR qPCR Master Mix (Code: 341-02, Vazyme, Nanjing, China) was used for qRT-PCR analysis and gene-specific primer pairs ( Table S2 ). The relative expression levels of the target genes were evaluated using the 2^-△△Ct approach with Spβ-actin as the reference gene (Livak and Schmittgen, 2001). Each sample was analyzed in three biological replicates.

The genes of homoplantaginin pathway with known function (SpFNS, SpF6H1, SpF6OMT2, and SpUGT1) were inserted in different combinations into the vectors pET32a and pACYCDuet-1 downstream of T7 promoters. The 3′ end of SpFNS and SpF6H1 were fused to the 5′ end of AtCPR via the linker sequence ACTAGTGGTTCTACCTCTTCTGGTTCTGGT. The 5′ end of SpUGT1 connects the MBP sequence of pSJ8-MBP vector. All the recombinant E. coli strains used in this study are shown in Table S3 . The recombinant plasmids were transformed into E. coli after sequencing verified. The recombinant strains were cultured in LB medium with antibiotics at 37°C and 200 rpm until the OD₆₀₀ reached 0.6-0.8. Then, 0.5 mM IPTG and 100 μM naringenin were added and cultured at 20°C and 135 rpm. After 24 h, 10 mL fermentation liquor was taken and mixed with a double volume of ethyl acetate to extract the reactions. The ethyl acetate layer was dried under reduced pressure, and dissolved with 100 μL methanol for HPLC analysis.

The chemical compositions of different organs of S. plebeia during growth and development were analyzed using UPLC-Q-TOF-MS/MS method. The total ion chromatograms (TICs) of different organs and developmental stages of S. plebeia are shown in Figures 2 , S1 . Combined with information such as the parent ion, molecular formula, retention time, secondary fragmentation, and relevant literature, a total of 70 compounds including 46 flavonoids, 16 phenolic acids, seven terpenoids, and one organic acid, were accurately or tentatively identified ( Table 1 ; Figure 3 ). We identified 46 flavonoids, including 24 flavones, 17 flavonones, 4 flavanones and one flavanonol, of which 21 compounds were discovered in S. plebeia for the first time ( Table 1 ). Especially, the peaks F2 and F5 were tentatively identified as new flavanone diglycosides from S. plebeia. The characteristic product ions and neutral losses are shown in the Figure S2 . From the TICs and Venn graph ( Figures 2 , S3 ), we can see that there was small difference in the number of the identified compounds in the samples of different developmental stages, which 43 compounds were common. However, the number of the identified compounds in different organs varies greatly. For instance, 68 compounds were identified in the leaf, while only 21 compounds were identified in the root. Principal component analysis (PCA) also revealed that the stems and roots, flowers, and leaves clustered into three groups ( Figure 4A ). The distribution of compounds in different organs of S. plebeia is quite different. There was little difference in the compounds present in the roots and stems at different stages. Therefore, we selected leaves from different periods and different organs from the flowering stage for differential metabolite analysis. The groups showed clear separation in the OPLS-DA score plots of S. plebeia with satisfactory goodness of fit and statistical significance ( Figure S4 ). Twenty-four of the 70 annotated metabolites were found to be differential metabolites (DEMs) ( Table S4 ). Based on the relative contents of the DEMs in the different parts at different developmental stages, a heatmap was created ( Figure 4B ) and revealed differences. The compounds were mainly distributed in leaves and flowers and less existed in the stems and roots. Specifically, the flavonoids were the main DEMs and mainly distributed in the leaves, while terpenoids mostly existed in the flowers. Besides, the flavonoid metabolites in different developmental stages showed differences. The main DEMs were flavonoid aglycones in L1, and glycosides in L2, L3, and L4.

To further accurately explore the contents of the main flavonoids in different parts at different plant developmental stages, 10 flavonoids (6-hydroxyluteolin 7-O-β-glucoside, nepitrin, luteolin 7-O-β-glucoside, apigenin 7-O-β-glucoside, homoplantaginin, luteolin, nepetin, 6-methoxynaringenin, apigenin, and hispidulin) were quantitatively analyzed using HPLC-DAD. Results indicated that flavonoid content varied among different organs and at different growth periods of the plant ( Figures 4C, D ). Regarding the distribution of flavonoids among different organs, the flavonoids were mainly distributed in leaves and flowers and to a less extent in the stems and roots. In different growth and development periods of the plant, flavonoid content peaked at the stem elongation stage. Specifically, flavonoid content in the leaves followed a parabolic trend and peaked when the stem elongated (L2, 18.21 ± 0.33 mg/g). Similarly, the flavonoid content in the stems was also the highest in S2 (8.97 ± 0.25 mg/g). In addition, the flavonoid content in flowers was up to 16.73 ± 0.16 mg/g (F3). A fluctuating trend of increase, decrease, and increase trend was observed in the roots; however, the highest yield was only 0.94 ± 0.08 mg/g at the ripening stage (R4).

Notably, hispidulin and homoplantaginin exhibited the highest levels in all the organs at each developmental stages of the plant. Our results showed that they accumulated differently in different growth stages ( Figure S5 ). The content of hispidulin reached the highest level of 6.55 ± 1.44 mg/g in the leaves at the basal leaf stage, and dropped to about 1 mg/g from vegetative growth to reproductive stage, but remained abundant in flowers (6.07 ± 0.43 mg/g). Homoplantaginin content in leaves gradually increased from 3.63 ± 0.59 to 6.38 ± 0.68 mg/g from the basal leaf to the ripening stage. A similar trend was observed for accumulation of 6-hydroxyluteolin 7-O-β-glucoside, luteolin 7-O-β-glucoside, and nepitrin. As the glycosides are maybe more stable than aglycones in the plant (Bowles et al., 2005), aglycones might transform into glycosides during the accumulation process. The homoplantaginin content in the stem was high in the vegetative growth period, but decreased during the reproductive stage. In addition, similar to hispidulin, homoplantaginin was abundant in flowers (5.73 ± 0.59 mg/g).

To further study the mechanisms of flavonoid biosynthesis in different organs and developmental stages of S. plebeia, 21 samples were used for transcriptomic analysis. The analysis produced approximately 161.44 Gb of clean data and 58,905 unigenes (140,509 transcripts) with a mean length of 1161.286 nt and an N50 of 1754 nt, following de novo assembly. The Q30 and GC percentages ranged from 94.46% to 95.15% and from 47.3% to 49.96%, respectively ( Tables S5, S6 ). The KEGG, NR, GO, Swissprot, KOG, COG, TrEMBL, Pfam, and eggNOG databases were used to annotate the predicted protein sequences. Finally, 44,185 unigenes were functionally annotated ( Table S7 ). Results from qPCR validation were mostly consistent with those obtained using RNA-seq, indicating the accuracy of the sequencing data ( Figure S6 ). Accordingly, the RNA-seq data were used for further analysis.

We then screened the differentially expressed genes (DEGs) in different organs and developmental stages of S. plebeia. There were 7896 DEGs in F3, 6366 DEGs in R3, 4471 DEGs in S3, 5112 DEGs in L1, 3798 DEGs in L2, 4211 DEGs in L4 compared with those in L3 ( Figure 5A ; Table S8 ). The number of DEGs in different tissues was higher than that in the developmental stages, which were consistent with the trend of metabolite accumulation. The statistical significance of the differences in gene expression in S. plebeia in different tissues and developmental stages was also represented using the volcano plots ( Figure S7 ). The number of DEGs varied in different groups, and there were more special DEGs in the flowers and roots ( Figure 5B ). DEGs were further assigned to the KEGG pathway. Interesting, the metabolites involved in “flavonoid biosynthesis” and “phenylpropanoid biosynthesis” were significantly enriched in different tissues and developmental stages ( Figure S8 ).

Next, we focused on the genes involved in flavonoid biosynthesis. A total of 27 annotated flavonoid biosynthesis genes were screened from the differential transcript analysis. Based on their expression ( Table S9 ), a heatmap were drawn ( Figures 5C, D ). We found that at different stages, most flavonoid biosynthesis genes, including C4H (c128957.graph_c4), CHI (c124024.graph_c0), CHS (c82070.graph_c0), F3H (c130349.graph_c0), and DFR (c83812.graph_c0), were upregulated in L2 ( Figure 5C ). Notably, UGT (c101658.graph_c1, c115214.graph_c0, c115878.graph_c0, c118512.graph_c0, c118635.graph_c0, c127652.graph_c0, and c128814.graph_c0) and CCOMT (c103796.graph_c0 and c119249.graph_c0) were upregulated in L4 ( Figure 5D ). Specifically, UGT and CCOMT are the key enzymes involved in the synthesis of homoplantaginin and hispidulin, respectively. This was consistent with the highest concentration of homoplantaginin observed at L4. At the flowering stage, the annotated flavonoid genes were highly expressed in flowers, such as CHS (c122262.graph_c0 and c82070.graph_c0), DFR (c119231.graph_c0 and c83812.graph_c0), F3H (c115905.graph_c0, c130349.graph_c0, and c100929.graph_c1), and CHI (c124024.graph_c0). Some OMT (c107876.graph_c0 and c118570.graph_c0), CCOMT (c103796.graph_c0, c119249.graph_c0, and c113142.graph_c1), and UGT (c115214.graph_c0, c118512.graph_c0, c127652.graph_c0, c128814.graph_c0, and c111184.graph_c1) were upregulated in S3 and R3.

The biosynthetic pathways of flavonoids have been extensively studied in Arabidopsis thaliana, Scutellaria baicalensis, Carthamus tinctorius (Saito et al., 2013; Xu et al., 2020; Wang et al., 2021). Hence, the biosynthetic pathway from L-phenylalanine to homoplantaginin was predicted, comprising nine consecutive steps ( Figure 6A ). To fully screen for the possible unigenes involved in the biosynthesis of hispidulin and homoplantaginin, we mined the candidate genes by combining differential transcript analysis with transcriptome annotation. A total of 49 candidate genes, including two PAL, 17 4CL, two C4H, three CHS, two CHI, five F6OMT, and 18 UGT, were obtained. However, FNS and F6H were not annotated in the transcriptome analysis. Therefore, we searched for these two genes using the local BLAST and relevant literature. FNS has been reported to be divided into FNS I and FNS II. FNS I is mainly found in Umbelliferae, whereas FNS II is widely distributed in plants containing flavonoids. Therefore, we chose FNS I (GenBank: AY817680.1) from Petroselinum crispum (Umbelliferae) and FNS II (GenBank: KP337723.1) from Salvia miltiorrhiza (Labiatae) as templates to conduct a similarity search against the S. plebeia RNA-sequencing (RNA-Seq) datasets. However, only one gene (c127704.graph_c0) with high similarity to FNS II was found. The phylogenetic analysis demonstrates that SpFNS clustered with members of the FNS II clade and was closest in affinity to the CYP93B25 from Salvia officinalis, and FNS II from Ziziphora clinopodioides ( Figure S9 ). F6H is mainly from the CYP71D family of Leguminosae and the CYP82D family of Labiatae. Considering the phylogenetic relationships among the species, the F6H of S. plebeia is more likely to have originated from the CYP82D family. Hence, c120693.graph_c0, c125994.graph_c1, and c127863.graph_c4 were annotated as CYP82D family in our datasets, which named SpF6H1, SpF6H2, and SpF6H3. These three genes cluster in family CYP82D with CYP82D1.1 from Scutellaria baicalensis, which has been shown to have broad catalytic activity for 6-hydroxyflavonoids such as apigenin ( Figure S10 ).

Finally, 53 unigenes were selected as candidate genes ( Table S10 ). The sequences of the 53 genes were used as the queries in a BLAST algorithm-based search of NCBI databases to predict their possible functions and the open reading frame (ORF) sequences. Finally, 28 candidate genes were selected to validate their catalytic functions ( Table S1 ).

Specifically, SpPAL, SpC4H, Sp4CL1-6, SpCHS1-3, SpCHI, SpFNS, SpF6H1-2, SpF6OMT1-3, and SpUGT1-10 were expressed in E. coli, respectively. The recombinant proteins were extracted and purified for biochemical analysis in vitro. HPLC-DAD and LC-MS analyses showed that SpPAL can catalyze the formation of cinnamic acid from phenylalanine ( Figure 6B ). When AtCPR proteins and NADPH were added, SpC4H was able to catalyze the production of 4-coumaric acid from cinnamic acid ( Figure 6C ). As we do not have the standard of 4-coumaryl-CoA, we identified the reaction products of Sp4CL2 and Sp4CL5 by comparing the different peak areas of 4-coumaryl-CoA at 330 and 260 nm according to a previous report ( Figure 6D ) (Li et al, 2010). Similarly, in order to verify the enzymatic activity of SpCHS1–3, we used a binary plasmid with pCDFDuet-RtTAL-Pt4CL and pET28a-SpCHS1–3 co-expressed in E. coli. The co-expression of RtTAL and Pt4CL can directly catalyze tyrosine into 4-coumaryl CoA, providing a sufficient substrate for the catalytic reaction of SpCHS. Upon feeding with tyrosine, naringenin chalcones were found in the culture extracts ( Figure 6E ). Although naringenin chalcone can be converted to naringenin spontaneously, it can be converted almost completely instantaneously in the presence of CHI. As shown in Figure 6F , some naringenin was also generated in the control group, but naringenin chalcone was completely consumed in the sample group to generate naringenin, indicating the catalytic activity of SpCHI. FNS and F6H, which belong to the CYP450 family and require the NADPH-CPR to provide electrons, are similar to C4H. AtCPR was co-expressed with SpFNS and SpF6H1-2, respectively. SpFNS produced apigenin from naringenin ( Figure 6G ). SpF6H1 produced scutellarein from apigenin, a precursor of hispidulin ( Figure 6H ). SpF6OMT1 and SpF6OMT2 converted scutellarein into hispidulin ( Figure 6I ). Thus, SpPAL, SpC4H, Sp4CL2, Sp4CL5, SpCHS1, SpCHI, SpFNS, SpF6H1, SpF6OMT1, and SpF6OMT2 formed a complete biosynthetic pathway for hispidulin. In the final step of homoplantaginin biosynthesis, SpUGT1, SpUGT2, and SpUGT3 produced homoplantaginin from hispidulin in vitro using UDP-glucose as a sugar donor ( Figure 6J ). In sum, all the enzymes involved in the hispidulin and homoplantaginin biosynthetic pathways were identified.

To further investigate the key enzymes involved in the hispidulin and homoplantaginin pathway in the S. plebeia, we reconstructed the pathway leading from naringenin to homoplantaginin in E. coli ( Figures 7A ). Simultaneously introduced SpF6H1 and SpUGT1 into E. coli to harbor AtCPR, the recombinant E. coli HP1 was able to produce scutellarein and homoplantaginin ( Figures 7B, C ). When SpFNS and SpF6OMT2 were simultaneously introduced into E. coli, which harbored the AtCPR, the resulting recombinant E. coli HP2 was able to produce apigenin and hispidulin ( Figures 7D, E ). The pET32a vector carrying SpF6H1 and SpUGT1 and the pACYCDuet-1 vector carrying SpFNS and SpF6OMT2 were transferred into E. coli to construct HP3. HP3 was able to produce hispidulin and homoplantaginin under naringenin induction. After 48 h of cultivation, the accumulation of hispidulin and homoplantaginin reached to 5.33 and 3.86 mg/L, respectively ( Figure 7F ). Although the yield was not desirable, it could be increased through metabolic engineering and fermentation optimization to make the full use of the E. coli fermentation system in the future.