Two-year-old Gastrodia elata Bl. f. elata (red stem, GR), Gastrodia elata Bl. f. glauca S. Chow (black stem, GB), and Gastrodia elata Bl. f. glauca S. Chow and Gastrodia elata Bl. f. elata (hybridization of GR and GB, GR×GB) were collected across four Chinese provinces (Anhui, Hubei, Guizhou, and Yunnan), a total of 8 groups, each included three replicates, and each replicate contained one plant. Detailed sampling information is presented in Supplementary Table S1 . All samples were collected by the same individual in November and were identified by Prof. Huasheng Peng (Anhui University of Chinese Medicine). After collection, they were cleaned of surface soil, blotted dry with filter paper, sliced, immediately flash-frozen in liquid nitrogen, and subsequently stored at –80 °C for further analysis.

The soluble sugar content and composition were determined using MetWare (Wuhan, China, http://www.metware.cn/). A total of 20 mg powder was diluted to 500 μL methanol: isopropanol: water (3:3:2, v/v/v), vortexed for 3 min, and sonicated for 30 min. The extract was centrifuged at 12,000 rpm under 4 °C for 3 min. The supernatant (50 μL) was mixed with 20 μL internal standard (1000 μg/mL) and evaporated under a nitrogen gas stream. The evaporated samples were transferred to a lyophilizer for freeze-drying. The residue was used for further derivatization. The derivatization method was as follows: the sample was mixed with 100 μL solution of methoxyamine hydrochloride in pyridine (15 mg/mL), incubated at 37 °C for 2 h, and 100 μL of N,O-bis (trimethylsilyl) trifluoroacetamide (BSTFA) was added. The mixture was incubated at 37 °C for 30 min after vortex-mixing. The mixture was analyzed using gas chromatography-mass spectrometry (GC-MS) following dilution to an appropriate concentration (Gómez-González et al., 2010; Medeiros and Simoneit, 2007; Chen et al., 2016).

The Agilent 8890 gas chromatograph coupled to a 5977B mass spectrometer with a DB-5MS column (30 m length × 0.25 mm i.d. × 0.25 μm film thickness; J&W Scientific, USA) was employed for GC-MS analysis of soluble sugars. Helium was used as a carrier gas at a flow rate of 1 mL/min. Injections were made, with an injection volume of 1 μL, in split mode with a split ratio of 5:1. The oven temperature was held at 160 °C for 1 min, and then raised to 200 °C at 6 °C/min, 270 °C at 10 °C/min, 300 °C at 5 °C/min, and finally 320 °C at 20 °C/min and held for 5.5 min. All samples were analyzed in the selective ion monitoring mode. The ion source and transfer line temperatures were 230 °C and 280 °C, respectively (Medeiros and Simoneit, 2007; Chen et al., 2016). Principal component analysis (PCA) and hierarchical cluster analysis (HCA) were used to analyze the accumulation patterns of the metabolomic data. The HCA results of the samples and metabolites are presented as heat maps accompanied by dendrograms. The fold change (FC) ≥ 2 or fold change (FC) ≤ 0.5 was set as the screening condition for differentially accumulated metabolites (DAMs).

As previously described, with minor modification, the phenolic content was extracted from dried G. elata samples and measured (Chen et al., 2016). Briefly, the dried powder (0.5 g) of each G. elata sample was extracted in 25 mL of 50% methanol via ultrasonication for 2 h. The phenols in the extract were chromatographed and detected using an Acquity-UPLC™ (Waters, MA, USA) system equipped with a quaternary pump, degasser, autosampler, and photodiode array detector (PDA). Phenols were separated on an ACQUITY UPLC^® T3 column (1.8 μm, 2.1 mm × 100 mm) (Waters, MA, USA). The mobile phase consisted of water containing 0.05% phosphoric acid (A) and acetonitrile (B) with a flow rate of 0.2 mL/min at 25 °C. The linear gradient elution program was as follows: 0~2 min, 5%~7% B; 2~3 min, 7%~10% B; 5~7 min, 20%~24% B; 7~9 min, 15%~18% B; 9~12 min, 18%~20% B; 12~14 min, 80%~76% B; 14~19 min, 24%~27% B; 19~21 min, 27%~15% B; 21~24 min, 15%~5% B. The injection volume was 2 μL, and the chromatograms were acquired at 220 nm.

AH-GR, HB-GR, GZ-GR, AH-GR×GB, HB-GR×GB, YN-GR×GB, GZ-GB and YN-GB were used for RNA-seq analysis. Total RNA was extracted by ethanol precipitation and CTAB-PBIOZOL, and genomic DNA was removed from the total RNA using DNase I. Subsequently, total RNA was identified and quantified using a Qubit fluorescence quantifier and Qsep400 high-throughput biofragment analyzer. The cDNA libraries were constructed using MetWare Biotechnology Co., Ltd. (Wuhan, China). The final library was amplified with phi29 before being uploaded to make DNA nanoballs (DNB) with more than 300 copies of a molecule, and the DNBs were loaded onto the sequencing chip and sequenced on the BGI sequencing platform.

To obtain high-quality clean data, fastp software was used for strict quality control of the data (Chen et al., 2018). The adapters of reads are removed simultaneously, along with low-quality reads (including those with an N ratio greater than 10% and those with Q ≤ 20 bases accounting for more than 50% of the entire read). The HISAT2 software was used to compare the clean reads with the reference genome to obtain localization data of the reads on the reference genome (Kim et al., 2015). StringTie was used to assemble the reads and reconstruct the transcriptomes for subsequent analyses (Pertea et al., 2015). Finally, annotation data were obtained by sequence alignment of the genes with the Kyoto Encyclopedia of Genes and Genomes (KEGG), Gene Ontology (GO), Non-Redundant Protein Database (NR), Swiss Protein Database (Swiss-Prot), Translated European Molecular Biology Laboratory Nucleotide Sequence Database (TrEMBL), and EuKaryotic Orthologous Groups (KOG) databases using Diamond, with an E-value of 1e^-5.

The Stringtie maximum flow algorithm and fragments per kilobase of transcript per million fragments mapped (FPKM) were used to calculate gene expression levels. DESeq2 (Love et al., 2014) was used for differentially expressed genes (DEGs) screening based on the count values of genes in each sample. The false discovery rate (FDR) was obtained by correcting the probability of hypothesis testing (P-value) for multiple hypothesis testing using the Benjamini-Hochberg method. DEGs were screened for Fold Change ≥ 2 and FDR < 0.01. The total numbers of DEGs, upregulated genes, and downregulated genes in each group were counted and visualized using OmicShare tools (https://www.omicshare.com/tools/) and EVenn (http://www.ehbio.com/test/venn/#/).

Phylogenetic trees were constructed with MEGA 6.0 software using the neighbor-joining method with 1,000 bootstrap replicates. The amino acid sequence of the UGT from Arabidopsis thaliana was obtained from the CAZy database (http://www.cazy.org/) (Ross et al., 2001). Several UGT proteins that have been confirmed to encode the function of catalyzing the conversion of p-hydroxybenzyl alcohol to gastrodin in Escherichia coli and Saccharomyces cerevisiae were used to construct phylogenetic trees, namely UGT73B6 (GenBank: AY547304) (Bai et al., 2014, 2016; Zhao et al., 2021), UGT73B6FS (Bai et al., 2016), AsUGT (GenBank: Q9AR73) (Yin et al., 2022), AsUGTsyn (Yin et al., 2022), AtUGT (GenBank: KAG7614719.1) (Xia et al., 2024), CtUGT (GenBank: AJT58578.1) (Xia et al., 2024), PgUGT (GenBank: AIE12483.1) (Xia et al., 2024), and RgUGT (GenBank: KAG5565971.1) (Xia et al., 2024). The constructed phylogenetic tree was visualized and identified using the iTOL tool (https://itol.embl.de/).

Multiple sequence alignments were performed using the ESPript online program (https://espript.ibcp.fr/ESPript/cgi-bin/ESPript.cgi). A three-dimensional model of GeUGTs was generated using the SWISS-MODEL server and visualized using PyMOL 3.0. Molecular docking was performed using AutoDock 1.5.7 to predict the binding of the glycosyl donor UDPG to the GeUGT protein, thereby forming the GeUGT-UDPG complex. Subsequently, p-hydroxybenzyl alcohol was docked onto the GeUGT-UDPG complex. The protein structure modeling and molecular docking are outlined in Supplementary Tables S2 and S3 .

To validate the reliability of the RNA-seq results for determining gene expression, qRT-PCR was performed. Based on our preliminary studies and literature survey, the β-actin gene was selected as the reference gene (Yuan et al., 2018; Shan et al., 2021). Twelve DEGs involved in phenolic compound and carbohydrate metabolism were selected for qRT-PCR validation. Gene-specific primers were designed using the Primer 5.0 software and synthesized by General Biology Co., Ltd. (Anhui, China; Supplementary Table S4 ). qRT-PCR analysis was conducted using the BioRad CFX96 Real-Time PCR Detection System (Bio-Rad Laboratories) following the manufacturer’s protocol (LABLEAD, China), including optimized cycling conditions and reaction mixture proportions. Relative gene expression levels were quantified using the 2^−ΔΔCt method. All experiments were performed in triplicate to ensure data reproducibility.

Microsoft Excel and GraphPad Prism 9.30 software were used to statistically analyze the data. The statistical analysis method used was one-way ANOVA. The TB tools software 2.331, Metware Cloud (https://cloud.metware.cn), and OmicShare tools (https://www.omicshare.com/tools/) were used for visualization and analysis.

The morphology and appearance of G. elata tubers vary considerably across different cultivars. The GR exhibits an elongated cylindrical shape, whereas the GB is characterized by a distinct oval form. The GR×GB effectively combines the favorable morphological traits of both parental lines, demonstrating a more robust and plump tuber morphology than its parents ( Figure 1A ). Soluble sugars in G. elata tubers from different producing areas and cultivars were analyzed using GC-MS. The total ion overlap diagram chromatograms demonstrated excellent instrument stability ( Supplementary Figure S1A ). PCA revealed a clear separation among the samples ( Supplementary Figure S1B ). A total of 19 soluble sugar components were identified ( Supplementary Table S5 ), among which sucrose, glucose, D-fructose, and D-xylose were the major soluble sugar components in G. elata, collectively accounting for over 98% of the total soluble sugar content. Sucrose (50.91–80.81 mg/g) was the most abundant across all samples ( Supplementary Table S5 ). Comparatively, the content of glucose, D-fructose, and D-xylose in GR and GB was higher than that in GR×GB ( Figure 1B ). DAMs primarily consisted of monosaccharides and disaccharides. Specifically, 3, 5, 3, 1, 3, and 3 DAMs were detected in the comparison groups HB-GR vs. AH-GR, GZ-GR vs. AH-GR, GZ-GR vs. HB-GR, HB-GR×GB vs. AH-GR×GB, AH-GR×GB vs. YN-GR×GB, and YN-GB vs. GZ-GB, respectively ( Supplementary Table S6 ), and identified a DAM (xylitol) in the GR from three regions ( Supplementary Figure S2 ). Additionally, in the comparative groups from the same regions of Anhui, Hubei, Guizhou, and Yunnan, there were 7, 3, 3, and 8 DAMs ( Supplementary Table S6 ), respectively.

Phenolic compounds are the major class of bioactive constituents in G. elata (Shan et al., 2021). In all examined samples, parishin A, parishin B, and parishin E were present at relatively high concentrations, whereas p-hydroxybenzyl alcohol and parishin C were detected at comparatively lower levels ( Supplementary Table S7 ). In the GR group, the gastrodin content in AH-GR and GZ-GR was significantly higher (p < 0.05) than that in HB-GR. The p-hydroxybenzyl alcohol content in AH-GR was significantly higher than that in GZ-GR, while the parishin E content in HB-GR was significantly higher (p < 0.05) than that in AH-GR. In the GR×GB group ( Figure 2 ), significant differences were observed in the contents of parishin A and parishin E between the HB-GR×GB and YN-GR×GB. In the GB group, the content of p-hydroxybenzyl alcohol in YN-GB accessions was higher than that in GZ-GB ( Figure 2 ). Differences were also observed in the content of specific compounds among different varieties of G. elata from the same origin. The gastrodin and p-hydroxybenzyl alcohol content in AH-GR were significantly higher (p < 0.05) than those in AH-GR×GB, the gastrodin content in GZ-GR was significantly higher than that in GZ-GB, and the p-hydroxybenzyl alcohol content in YN-GB was significantly higher than that in YN-GR×GB ( Figure 2 ).

A total of 24 cDNA libraries were subjected to transcriptome sequencing and yielded between 43,248,992 and 56,780,560 clean reads; the Q30 and Q20 percentages for all libraries exceeded 92.54% and the GC content percentages were > 45.43% ( Supplementary Table S8 ). The correlation coefficient between samples exceeded 0.85 ( Supplementary Figure S3 ). The PCA plot shows that PC1 and PC2 account for 26.4% and 18.1% of the variance, respectively ( Supplementary Figure S4 ). Overall, the results of transcriptome analysis were of high quality and suitable for further research. Of the 23,162 identified genes, 18,415 (79.51%) were successfully annotated using at least one of the following databases: Nr, SwissPort, KOG, Trembl, Pfam, GO, or KEGG ( Supplementary Table S9 ; Supplementary Figure S5 ). Among the 18,347 genes annotated in the NR database, 87.88% showed homologous sequence matches with eight species of the Orchidaceae family, with the highest homology observed with Dendrobium catenatum ( Supplementary Figure S6 ). Among the 15,504 genes annotated with GO terms, the categories were predominantly “cellular process” and “metabolic process” in biological processes, “cellular anatomical entity” in cellular components, and “binding” and “catalytic activity” in the molecular function categories ( Supplementary Figure S7 ). Additionally, 14,161 genes were mapped to 147 pathways using KEGG annotation ( Supplementary Figure S8 ).

The analysis also revealed that the number of DEGs in the three pairwise comparisons within the GR group from different producing areas was substantially higher than that in the GR×GB and GB groups, reaching 3,255, 4,866 and 4,037, respectively. The GR×GB groups contained 1,791, 608 and 2,983 DEGs, while the GB group showed the lowest number, with a total of only 617 ( Figure 3A ). Furthermore, the number of DEGs in the four comparison groups of different varieties within the same origin was generally low. Particularly, the comparison between HB-GR×GB and HB-GR yielded only common 50 DEGs ( Supplementary Figure S9 ). Based on the identified differences from the comparisons, we conducted further analysis on G. elata of the same variety from different producing areas. Venn diagrams illustrated 682 common DEGs in the three GR comparison groups ( Figure 3B ), and 93 in the three GR×GB comparison groups ( Figure 3C ). Further analysis of the DEGs across different groups was conducted to investigate their functional categories using the KEGG database. The highest number of DEGs was identified in the comparison AH-GR_vs._GZ-GR within the GR group ( Figure 3D ), with an enrichment of 139 KEGG pathways. Conversely, the fewest DEGs were identified in the AH-GR×GB_vs._YN-GR×GB within the GR×GB group ( Figure 3E ), with 159 DEGs enriched in 93 KEGG pathways. In the GB comparison group, GZ-GB_vs._YN-GB ( Supplementary Figure S10 ) revealed enrichment of 177 DEGs across 102 KEGG pathways. Seven comparison groups demonstrated significant enrichment in the “metabolic pathway,” “secondary metabolite biosynthesis,” and “plant-pathogen interaction” pathways. Notably, pathways related to phenolic metabolism, such as phenylpropanoid biosynthesis, phenylalanine metabolism, and the biosynthesis of phenylalanine, tyrosine, and tryptophan, as well as glucose metabolic pathways, including starch and sucrose metabolism, galactose metabolism, and glycolysis/gluconeogenesis, were enriched.

In our previous study, we postulated the biosynthetic pathway of phenolic compounds in G. elata (Shan et al., 2021) ( Figure 4B ). Phenolic metabolite analysis in G. elata revealed significant variations in phenolic content across different producing areas and varieties. Based on these findings, we next investigated the molecular mechanisms underlying the differential accumulation of phenolic compounds in G. elata. Gastrodin, p-hydroxybenzaldehyde, and p-hydroxybenzyl alcohol are synthesized from phenylalanine via the phenylpropanoid pathway, whereas the citric acid biosynthesis pathway is implicated in the biosynthesis of parishin, we identified a total of 11 enzymes: Phenylalanine ammonialyase (PAL), Cinnamate4-hydroxylase (C4H), 4-comarate coenzymeA ligase (4CL), Shikimate O-hydroxycinnamoyltransferase (HCT), 5-O-(4-coumaroyl)-D-quinate3’-monooxygenase (C3H), Caffeoyl-CoA O-methyltransferase (CCoAOMT), Alcoholdehydrogenase (ADH), Glycosyltransferases (GT), Pyruvatedehydrogenase/pyruvate dehydrogenase complex (aceE), Dihydrolipoyllysine-residue acetyltransferase (DLAT) and Citrate symthase (CS) corresponding to 96 genes ( Supplementary Table S10 ). The phenylpropanoid biosynthesis pathway serves as the fundamental route for the synthesis of various phenolic compounds in plants, with its upstream steps primarily responsible for generating universal precursors. For our target metabolites, that is, gastrodin, p-hydroxybenzyl alcohol, and parishin-type compounds, the downstream ADH and GT genes represent more critical targets for investigation. The expression of ADH1 and ADH4 was lower in AH-GR than AH-GR×GB, whereas the expression of ADH5, ADH6, and ADH11 was higher ( Figure 4B ), which may be attributed to the higher content of p-hydroxybenzyl alcohol and gastrodin in AH-GR. Correlation clustering analysis between six phenolic compounds and eleven ADH genes indicated a significant positive correlation (p < 0.05) of ADH5 and ADH11 with the accumulation of p-hydroxybenzyl alcohol and gastrodin ( Figure 4C ), suggesting that ADH5 and ADH11 are likely key genes involved in the biosynthesis of p-hydroxybenzyl alcohol. The content of gastrodin and p-hydroxybenzyl alcohol in AH-GR was significantly higher than that in AH-GR×GB ( Figure 2 ). We identified a total of 46 GT genes; there were 12 differentially expressed GTs (GT1, GT4, GT5, GT11, GT16, GT17, GT20, GT22, GT30, GT31, GT35, and GT37) between AH-GR×GB and AH-GR ( Figure 4D ).

Correlation clustering analysis between six phenolic compounds and 46 GT genes revealed ( Figure 5A ) that Group A (GT3, GT38, GT11, GT16, GT24, GT7, GT39, GT9, GT25, GT41, GT45, GT4, GT8, GT35, GT37, GT46, and GT18) and Cluster I (parishin A, parishin B, parishin C, p-hydroxybenzyl alcohol, and gastrodin), which showed positive correlations (p < 0.05), negatively correlated with Cluster II (parishin E). In Group C, GT19 and GT1 significantly positively correlated with parishin E (p < 0.05). These findings provide key candidate genes and a solid basis for further functional studies on glycosyltransferases in G. elata. To gain a comprehensive understanding of the evolutionary relationships and classify the 46 UGT proteins in G. elata, a phylogenetic tree including 25 UGT proteins from Arabidopsis thaliana (Groups A to N) and 5 from Zea mays (Groups O to Q) was constructed using MEGA 7.0. The phylogenetic tree revealed that the UGT proteins in G. elata could be classified into 13 groups, with no sequences identified in groups B, H, N, or P ( Figure 5B ). Among them, group D contained the largest number of UGT proteins of G. elata, with a total of 15 members, followed by group E, which contained 9 members. Groups G, L, and A consisted of five, four, and three members, respectively. Groups C, F, K, M, and O each contained only one member. Notably, group D contained two UGT members (CtUGT and UGT73B6FS), whereas group E contained 4 UGT members (AtUGT, AsUGTsyn, PgUGT, and RgUGT), all of which have been shown to catalyze the conversion of p-hydroxybenzyl alcohol to gastrodin in vitro. The 24 UGT genes in groups D and E of G. elata may also have similar catalytic functions.

Based on the phylogenetic analysis, four GeUGTs (GT16, GT17, GT18, and GT44) from group E and two GeUGTs (GT24 and GT42) from group D, which showed close evolutionary relationships with functionally characterized UGT genes, were selected for sequence and structural analyses. Multiple sequence alignments revealed that these six GeUGTs shared 48.81% sequence identity with the known functional UGTs. All identified GeUGTs contained a highly conserved 44-amino acid sequence at their C-termini, designated the plant secondary product glycosyltransferase (PSPG) motif, which is characteristic of plant secondary metabolism glycosyltransferases. Notably, the catalytic residues His (H) and Asp (D) are highly conserved among the UGTs ( Supplementary Figure S11 ). Subsequent tertiary structure modeling demonstrated that all six GeUGT structures were monomeric proteins. In the structural models, the PSPG box is highlighted in red, while the catalytic residues are represented as spheres, with histidine (His) in blue and aspartate (Asp) in yellow ( Figure 6 ).

Molecular docking was used to investigate the interaction between the protein and ligand and predicted the stability and binding affinity of the docked complex ( Supplementary Table S11 ). The molecular docking analysis revealed that p-hydroxybenzyl alcohol binds to the GT16-UDPG complex with a binding energy of -5.043 kcal/mol. Additionally, His48 formed hydrogen bonds with UDPG at a distance of 3.3 Å ( Figure 7A ). Similarly, p-hydroxybenzyl alcohol exhibited stronger binding affinity (-5.326 kcal/mol) to the GT17-UDPG complex, where the catalytic residue His25 formed shorter hydrogen bonds (2.2 Å) with UDP-glucose ( Figure 7B ). The key conserved residues in the PSPG box of GT16, GT17, GT18, and GT44 formed hydrogen bonds with UDP-glucose (UDPG), facilitating substrate recognition and binding ( Figures 7A-C, F ). However, in GT24 and GT42, neither the catalytic residue His nor the key residues within the PSPG box formed interactions with UDPG and p-hydroxybenzyl alcohol ( Figures 7D, E ).

To confirm the accuracy of the transcriptome sequencing results, six DEGs (PAL4, CCoAMT1, C4H, 4CL2, ADH4, UGE1, and INV8) involved in phenolic and soluble sugar metabolic pathways, and six GT genes (GT16, GT17, GT18, GT24, GT42, and GT44) were randomly selected for validation using qRT-PCR analysis. As shown in Supplementary Figure S12 , the expression profiles of the 12 genes were consistent in the qRT-PCR and FPKM values, confirming the accuracy and repeatability of the transcriptome analysis.