From the G. biloba nursery in Nanjing Forestry University (32° 04’N, 118° 48’E, Nanjing, China), 33 clones (namely 1#-33#, 18-year-old) with leaf utilization potential were selected. These clones were originated from Pizhou, Jiangsu Province, and were selected in 2005 based on seedling height and ground diameter. Afterward, these clones were transplanted into a common garden in Nanjing forestry university with uniform spacing (4m × 4m) to minimize environmental variation. In mid-April, the healthy functional leaves were harvested from every clone and immediately wrapped in tinfoil and frozen in liquid nitrogen. Leaf samples were ground into fine powder with a mortar and a pestle in liquid nitrogen. The fine-powdered samples were then stored at an ultralow temperature refrigerator for further analyses. For each clone, six branches that exhibited similar growth status and exposed to sunlight were selected. For physiological and metabolomic analyses, leaves from each branch were collected separately to form an individual sample, resulting in six biological replicates for each clone. For the transcriptomic analysis, equal amounts of leaf samples from every two branches were pooled together to form a mixed sample. Consequently, three mixed samples were obtained for each clone.

The concentration of shikimic acid was measured as described by Zelaya et al. (2011) and Joubert et al. (2023) with minor modifications. Briefly, ca. 50 mg oven-dried fine powders were mixed with 0.25 mol L^-1 hydrochloric acid for 30 min. After centrifuged (6000 g) for 10 min, the supernatant was collected and diluted three times. The reaction was initiated by adding 1% (w/v) periodate solution. Three hours later, the reaction was stopped by adding a solution containing 1 mol L^-1 NaOH and 1mol L^-1 glycinate. Shikimic acid concentration was measured spectrophotometrically at 380 nm.

Phenylalanine concentration was assayed using a phenylalanine content kit (Nanjing Jiancheng Bioengineering Research Institute Co., Nanjing, China) following the kit’s instructions.

Tryptophan and tyrosine were assayed as described previously (Park et al., 2012; Botella et al., 2023) with minor modifications. Briefly, frozen powder was homogenized in 100 µL extraction solution containing 50% ethanol and 0.1 mol L^-1 HCl. The mixture was centrifuged (13800 g, 4°C, 20 min), and the supernatant was filtered through an organic membrane (0.22 μm). The filtered supernatant was used to determine tryptophan and tyrosine with liquid chromatography-mass spectrometry (LC/MS, LTQ-XL, Thermo Scientific, Waltham, MA, USA).

The concentration of flavonoid in the leaves were analyzed as described previously (An et al., 2024; Cui et al., 2025). Briefly, oven-dried fine powder (ca. 200 mg) was packed with filter paper and washed with 100 mL petroleum ether in Soxhlet extraction to remove impurities for 8 h. Afterward, the purified packed fine powder was extracted in 10 mL methanol at 60 °C for 30 min, and this procedure was repeated twice. The extraction was transferred into a new volumetric flask, ensuring the final volume to 20 mL. The extraction was homogenized with 5% sodium nitrite, 10% aluminum nitrate, and 1 mol L^-1 NaOH for 5, 6 and 10 min, respectively. The absorbance of the mixture was determined spectrophotometrically at 510 nm.

Since shikimic acid content in 1# (referred as HS) was the highest while that in 24# (referred as LS) was the lowest ( Figure 1a ), 1# and 24# were selected to explore metabolomic analysis of shikimic acid synthesis and downstream metabolism. The clones of No.1 (HS) and No.24 (LS) originate from distinct individual trees and are not derived from mutated shoots of any single tree. Nontargeted metabolites in the leaves were analyzed based on the method described previously (Yu et al., 2022; Hong et al., 2023). Briefly, fresh fine powders (ca. 100 mg) were extracted with 500 μL of 80% methanol and were incubated on ice for 5 min. The mixture was centrifuged (15000 g, 4°C) for 20 min. The collected supernatant was diluted with mass spectrometry water until methanol content was 53%. Then, the samples were transferred to fresh tubes and centrifuged (15000 g, 4°C) for 20 min. Finally, the supernatant was collected and injected into a Vanquish UHPLC system (Thermo Fisher Scientific, Waltham, MA, USA) coupled with an Orbitrap Q ExactiveTM HF-X mass spectrometer (Thermo Fisher Scientific) for liquid chromatography tandem mass spectrometry (LC-MS/MS) analysis.

LC–MS/MS raw data were processed using the Compound Discoverer 3.1 (CD3.1, Thermo Fisher Scientific) to perform peak alignment, picking and area quantification. Subsequently, the precise qualitative and relative quantitative results were obtained by matching the peaks with mzCloud, mzVault and MassList databases. The metabolites were annotated using the KEGG database, HMDB database and LIPIDMaps database (Marco-Ramell et al., 2018). The metabolites with VIP > 1, P-value < 0.05 and log₂ (Fold change) ≥ 1 (or ≤ -1) were considered to be differentially regulated. Notably, electrospray ionization (ESI) was employed as the ion source, and was conducted in both positive and negative ion modes.

Total RNA from the leaves of HS and LS was isolated using a polysaccharide polyphenol plant total RNA kit (DP441, TianGen, Beijing, China), respectively. RNA integrity was evaluated using Agilent 2100 bioanalyzer. Subsequently, cDNA libraries were constructed and sequenced on illumina NovaSeq 6000 (Illumina, San Diego, CA, USA). The analysis of RNA sequencing data was performed as described earlier (Lu et al., 2022; Du et al., 2024). Briefly, the original sequencing data were filtered to obtain clean data. The clean reads were then aligned to the G. biloba genome (Guan et al., 2016) (Genome ID: 100613, available at http://gigadb.org/dataset/100613) using HISAT2 (v.2.0.5). StringTie software was used to assemble the new transcripts, and featureCounts (v.1.5.0-p3) was used to calculate the Fragments Per Kilobase Million (FPKM) of each gene. Significantly differentially expressed genes (DEGs) were identified on the basis of |log₂ (Fold change)| ≥ 1 and a false discovery rate (FDR) < 0.05. The clusterProfiler software (v.3.8.1) was used for gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses (Li et al., 2024b). GO terms and KEGG pathways with P-values < 0.05 were considered to be significantly enriched (Xiang et al., 2024). Three cDNA libraries were generated and sequenced for HS and LS, respectively. The sequencing data were deposited to the Genome Sequence Archive (GSA, https://ngdc.cncb.ac.cn/gsa/) under Project ID CRA025492.

DEGs were annotated and functional categorized as described by Yu et al. (2021); (Lu et al., 2024) with minor modifications. Briefly, the coding sequences of DEGs were retrieved from the G. biloba’s genomic database. The closest homologue of a G. biloba gene in Arabidopsis thaliana was identified by searching its coding sequence against the protein sequence database of Arabidopsis using translated nucleotide BLAST (BLASTX). Identifiers of Arabidopsis genes closest to these DEGs were then submitted to MapMan (http://mapman.gabipd.org/) for functional analysis.

Considering the critical roles of transcription factors (TFs) in shikimic acid biosynthesis and downstream metabolism, the Pearson correlation coefficient (R) between differentially expressed TFs and shikimic acid content was calculated. Candidate TFs were defined to have an absolute R-value greater than 0.92 and a P-value less than 0.05. To further identify key TFs regulating the accumulation of shikimic acid, R were calculated between these candidate TFs and DEGs involved in the shikimic acid biosynthesis and downstream metabolism. Significant correlations were defined as |R| > 0.95 with P-value < 0.05. Correlations between key TFs and DEGs were displayed using Cytoscape (v.3.10.1).

The differentially accumulated metabolites (DAMs) and DEGs involved in shikimic acid biosynthesis and downstream metabolism were used for the integrative analysis, and the Pearson’s correlation coefficients between them were calculated (Yu et al., 2022). A correlation was considered statistically significant if the absolute value of the |R| exceeded 0.8 with a P-value less than 0.05. Heatmaps were used to reveal the correlation between these DAMs and DEGs.

Total RNA was extracted as mentioned above. Quantitative RT-PCR (RT-qPCR) was conducted as described previously (Zhang et al., 2023). Specific primers were designed for each DEG, and Glyceraldehyde 3- phosphate dehydrogenase (GAPDH) was selected as an internal standard ( Supplementary Table S1 ).

Statgraphics (STN, St Louis, MO, USA) was employed to do statistical tests as described previously (Lu et al., 2019, 2023). The data was tested to explore the normality prior to the analysis. One-way analysis of variance (ANOVA) was employed, and the means were regarded to be significantly different if the P-value was less than 0.05 on the basis of ANOVA F-test.

The concentration of shikimic acid in the leaves varied greatly among the 33 G. biloba clones ( Figure 1a ). Notably, 1# had the highest shikimic acid level and 24# possessed the lowest content ( Figure 1a ). The shikimic acid content in HS was 119% higher than that in LS ( Figure 1a ). HS and LS were selected to further study the metabolomic and transcriptomic mechanism of shikimic acid accumulation in G. biloba leaves. The concentrations of flavonoids, phenylalanine, tryptophan and tyrosine, which are downstream metabolites of shikimic acid, in HS were 17%, 22%, 11% and 23% lower than those in LS, respectively ( Figure 1b ).

A total of 862 metabolites were identified in positive (451) and negative (411) ion modes, respectively ( Supplementary Table S2 ). Principal component analysis (PCA) demonstrated that the first two principal components explained 60.48% and 57.61% of the total variance in negative and positive ion modes, respectively, with PC1 serving as the dominant contributor in both analyses ( Figures 2a, b ). Partial least squares discriminant analysis (PLS-DA) showed clear metabolic differences between HS and LS ( Figures 2c, d ). In the PLS-DA model, PC1 accounted for 34.79% (negative mode) and 33.44% (positive mode), while PC2 explained 20.94% (negative mode) and 15.43% (positive mode) of variances, respectively ( Figures 2c, d ). For metabolites detected in negative ion mode, compared with LS, the abundances of 67 metabolites were higher, whereas those of 71 metabolites were lower in HS, respectively ( Supplementary Figure S1a ; Supplementary Table S3 ). 61 upregulated and 91 downregulated metabolites (in positive ion mode) were identified in HS vs. LS, respectively ( Supplementary Figure S1b ; Supplementary Table S3 ). According to KEGG pathway annotation, these DAMs were mainly involved in metabolic pathways, tryptophan metabolism, galactose metabolism, aminoacyl−tRNA biosynthesis, and arginine and proline metabolism ( Supplementary Figures S1c, d ).

The DAMs were divided into 10 categories, mainly including benzenoids, flavonoids, amino acids and hormones ( Figure 2e ; Supplementary Table S3 ). The abundances of most benzenoids, amino acids and hormones, were significantly lower in HS compared to those in LS ( Figure 2e ; Supplementary Table S3 ). For example, abscisic acid, indole-3-acetic acid, tyrosine, tryptophan, kynurenic acid, bergenin and 3-dehydroshikimic acid were downregulated in HS compared to LS ( Figure 2e ; Supplementary Table S3 ). The abundances of 12 flavonoid metabolites, such as luteolin, kaempferol-3-gentiobioside and wight one in HS were significantly lower when compared with those in LS ( Figure 2e ; Supplementary Table S3 ).

Based on the above physiological and metabolomic data, it was speculated that there are differences in gene expression patterns involved in shikimic acid synthesis and downstream metabolism between HS and LS. To test this hypothesis, the transcriptomes of leaves from HS and LS were sequenced. A total of six cDNA libraries from HS and LS were constructed, and 45.1-46.3 million clean reads were obtained from each library ( Supplementary Table S4 ). 91.10-91.94% of clean reads can be mapped to the G. biloba genome database ( Supplementary Table S4 ). Compared with LS, 671 upregulated and 873 downregulated genes were detected in HS ( Figure 3a ). Genes were randomly selected to validate the transcriptomic data using RT-qPCR ( Supplementary Figure S2 ).

GO enrichment analysis showed that the differentially expressed genes (DEGs) were significantly enriched in oxidoreductase activity, acting on paired donors, with incorporation or reduction of molecular oxygen (GO:0016705), transferase activity, transferring hexosyl groups (GO:0016758), xyloglucan: xyloglucosyl transferase activity (GO:0016762), hydrolase activity, acting on glycosyl bonds (GO:0016798) and polysaccharide metabolic process (GO:0005976) ( Figure 3b ). KEGG pathway enrichment analysis showed that DEGs were involved in phenylpropanoid biosynthesis, flavonoid biosynthesis, starch and sucrose metabolism, flavone and flavonol biosynthesis, and cyanoamino acid metabolism ( Figure 3c ). These results suggest that DEGs are involved in the biosynthesis and downstream metabolism of shikimic acid. MapMan was used to further classify these DEGs into functional categories, including photosynthesis, flavonoid metabolism, hormone metabolism, amino acid metabolism and transcriptional regulation ( Supplementary Table S5 ), which are closely correlated with synthesis and downstream metabolism of shikimic acid.

In the shikimic acid synthesis pathway ( Supplementary Figure S4 ), malate dehydrogenase (MDH) encoding a pivotal enzyme responsible for the conversion of l-malate into oxaloacetate (Rozova et al., 2015) and Ribose-5-phosphate isomerase (RPI) coding for a cytosolic ribose-5-phosphate isomerase that catalyzes the conversion of D-ribulose 5-phosphate to d-ribose 5-phosphate (Faria et al., 2016), were more upregulated in HS compared to LS ( Figure 4a ; Supplementary Table S5 ).

In the downstream metabolism of shikimic acid ( Supplementary Figure S4 ), a number of genes that are implicated in biosynthesis of aromatic amino acids and flavonoids were downregulated in HS compared to LS ( Figure 4b ). For instance, the transcript levels of chorismate synthase (CS) encoding an enzyme that catalyzes the conversion of 5-enolpyruvylshikimate-3-phosphate (EPSP) to chorismate was downregulated in HS compared to LS ( Figure 4b ; Supplementary Table S5 ). Similarly, the mRNA levels of 4-coumarate:CoA ligase (4CL) encoding an enzyme that catalyzes the ATP-dependent conversion of p-coumaric acid to p-coumaroyl-CoA, chalcone synthase (CHS) encoding a key enzyme that condenses p-coumaroyl-CoA with three molecules of malonyl-CoA to form naringenin chalcone, flavanone-3b-hydroxylase (F3H) encoding an enzyme that mediates the stereospecific 3β-hydroxylation of naringenin to dihydrokaempferol, and dihydroflavonol 4-reductase (DFR) encoding an enzyme that catalyzes the NADPH-dependent reduction of dihydromyricetin to leucocyanidin, were downregulated in HS compared to LS ( Figure 4b ; Supplementary Table S5 ).

In the category of hormone metabolism, the transcript levels of abscisic acid-responsive element-binding factor (ABF) and gem-relate 5 (GER5) involved in ABA signaling (Liu et al., 2019) were downregulated in HS compared to LS ( Supplementary Figure S3 ; Supplementary Table S5 ). The mRNA levels of BIG and cullin-associated and neddylation dissociated (CAND1) involved in auxin regulation (Gil et al., 2001; Chuang et al., 2004), and like aux1 3 (LAX3) and GH3.1 implicated in auxin transport (Kasahara et al., 2020), were downregulated in HS compared to LS ( Supplementary Figure S3 ).

A total of 126 differentially expressed TFs were identified in this study ( Supplementary Table S5 ). The highest abundance of TF family was MYB (18, 14.29%), followed by PTRs (12, 9.52%), C2H2 (12, 9.52%), HB (7, 5.56%) and AP2-EREBP (7, 5.56%) ( Figure 5a ). Notably, only four MYBs were upregulated, while the other 14 MYBs were significantly downregulated in HS compared to LS ( Figure 5b ). In addition, most C2H2, bHLH and ARF family genes were downregulated, while WRKYs and bZIPs were more upregulated in HS compared to LS ( Figure 5b ).

Pearson correlation analysis between the transcript levels of differentially expressed TFs and shikimic acid concentration identified 28 TFs with closely correlations (R ≥ 0.92 or ≤ -0.92) ( Supplementary Table S6 ). Among these, nine TFs, such as enhanced downy mildew 2 (EDM2) and set domain group 26 (SDG26), showed a positive correlation with shikimic acid content ( Supplementary Table S6 ). The remaining 19 TFs, such as constans-like 9 (COL9) and myb domain protein 50 (ATMYB50), displayed a negative correlation with shikimic acid content ( Supplementary Table S6 ).

Correlation analysis was performed to further calculate the correlation coefficients between these 28 TFs and DEGs participating in shikimic acid biosynthesis and downstream metabolism ( Figure 5c ). Each TF was closely correlated with 1–6 DEGs involved in shikimic acid biosynthesis and downstream metabolism ( Figure 5c ; Supplementary Table S7 ). Notably, related to AP2 11 (RAP2.11) and reproductive meristem 16 (REM16) had the most numbers of closely correlated DEGs ( Figure 5c ; Supplementary Table S7 ). RAP2.11 and REM16 were considered to be key TFs in regulating shikimic acid biosynthesis and downstream metabolism.

The correlation between DAMs and DEGs that are involved in shikimate biosynthetic and downstream metabolic pathways were analyzed ( Figure 6 ). Apigenin c-pentoside, kaempferol-3-gentiobioside and luteolin were three metabolites having the highest numbers of correlated DEGs ( Figure 6 ). Specifically, kaempferol-3-gentiobioside, apigenin c-pentoside and luteolin were positively correlated with 19, 8 and 19 DEGs, and were negatively correlated with 10, 20 and 8 DEGs, respectively ( Figure 6 ). In addition, beta glucosidase 7 (BGLU7), Arabidopsis thaliana peroxidase 34 (AtPERX34), peroxidase (PRX) and copper amine oxidase (CuAO) were found to have the largest numbers (21-22) of correlated DAMs ( Figure 6 ). Therefore, the three DAMs and four DEGs may play key roles in the shikimic acid biosynthesis and downstream metabolism in the leaves of G. biloba.

The biosynthesis and downstream metabolism of shikimic acid determine the final shikimic acid concentration in G. biloba leaves. Notably, in the upstream synthesis of shikimic acid, both MDH and RPI were significantly upregulated in HS compared to LS. Functionally, MDH catalyzes the formation of OAA (as a precursor of PEP), and RPI generates Ru5P, which can be transformed into E4P (Stover et al., 2000). E4P and PEP constitute the essential substrate pairs for shikimic acid biosynthesis (Kruger and von Schaewen, 2003). These findings suggest that the upregulation of MDH and RPI play key roles in higher shikimic acid concentration in HS vs. LS.

Previous studies have demonstrated that the attenuated downstream metabolic flux exhibit significant accumulation of shikimic acid in plants (Bochkov et al., 2011; Eloy et al., 2017). The biosynthesis of aromatic amino acids and flavonoids are two vital metabolic pathways of downstream shikimic acid metabolism. Correspondingly, the abundances of several aromatic amino acids and flavonoid-related metabolites, such as phenylalanine, tryptophan, tyrosine, dihydromyricetin, (+)-catechin and luteolin, were significantly lower in HS vs. LS. Based on the association analysis, luteolin was identified as a key metabolite. In this study, the expression levels of multiple genes, such as CS, CHS, DFR, PRX, ATPERX34, BGLU7 and CuAO, which are involved in biosynthesis of aromatic amino acids and flavonoids were markedly downregulated in HS compared to LS. Additionally, Integrated analysis of metabolomic and transcriptomic data revealed that PRX, ATPERX34, BGLU7, and CuAO may play pivotal roles in shikimate biosynthesis and downstream metabolism. In maize (Zea mays), suppression of CHS results in lower flavonoid production and higher accumulation of shikimic acid derivatives (Eloy et al., 2017). Similarly, inhibition of the CS in Petunia hybrida leads to reduced flavonoid accumulation (Zhong et al., 2020). These results suggest that the downregulation of metabolites and genes participating in downstream metabolism of shikimic acid contribute to higher concentration of shikimic acid in HS than that in LS.

Previous studies have shown that ABA and IAA could affect the metabolism of flavonoid (Jeong et al., 2004; Dong and Lin, 2021), probably regulating the shikimic acid accumulation in G. biloba leaves. For instance, exogenous application of ABA increases the levels of flavonoid-related metabolites, such as isorhamnetin-3-O-gallate and dihydromyricetin, in Rhododendron chrysanthum after UV-B radiation (Yu et al., 2024). ABA positively regulates flavonoid accumulation by upregulating F3’H under drought stress in G. biloba (Yu et al., 2022). Auxin is also reported to be positively correlated with flavonoid accumulation by upregulating CHS in Arabidopsis root galls (PÄSold et al., 2010). In this study, the abundances of ABA and IAA, and the mRNA levels of ABF, GER, BIG, CAND1, LAX3 and GH3.1, which are involved in ABA and IAA signaling and metabolism, were significantly downregulated in HS compared to LS. These results indicate that the downregulation of ABA and IAA, along with the lower expression of genes participating in ABA and IAA signaling and metabolism, play vital roles in limiting flavonoid biosynthesis, thus contributing to higher shikimic acid accumulation in HS than that in LS.

In this study, the most abundant TFs belonged to the MYB family members ( Figure 5a ). Furthermore, in the co-expression network ( Figure 5c ), four TFs—MYB6, ZP1, RAP2.11, and MTERF34—exhibited the most connections with key structural genes, suggesting their potential regulatory roles in regulating shikimic acid biosynthesis and downstream metabolism. Previous studies have revealed that MYBs, C2H2s and ERFs play key roles in regulating flavonoid biosynthesis (Li et al., 2024a), probably affecting the shikimic acid concentration in G. biloba leaves. For instance, overexpression of MYB6 brings about upregulation of DFR2, resulting in larger concentrations of anthocyanin and proanthocyanidins in poplar (Wang et al., 2019). As a C2H2 family member, ZAT could regulate flavonoid biosynthesis by interacting with the promoter regions of genes involved in flavonoid biosynthesis, such as CHS in Macadamia integrifolia under high temperature stress (Yang et al., 2023). AP2/ERF transcription factors can directly target the key genes in the flavonoid biosynthesis pathway, thereby regulating the synthesis of flavonoid in Solanum lycopersicum and Citrus (Zhao et al., 2020; Cao et al., 2024). These results suggest that TFs play key roles in regulating shikimic acid downstream metabolism, thus bringing about higher shikimic acid accumulation in HS than that in LS.