Utilizing HiFi sequencing technology with rigorous quality control, we generated 33,799,059,475 bp of data ( Supplementary Figure S2 ; Supplementary Table S1 ). Comparative analysis of the NT library identified the top 10 species as green plants, confirming no contamination in the sequencing data ( Supplementary Figure S3 , Supplementary Table S2 ). The assembled genome size was 482.30 Mb, with a Contig N50 of 30.97 Mb ( Supplementary Table S3 ). Chromosome-level assembly was achieved via Hi-C technology, yielding a 100% chromosome loading rate and pseudomolecular maps for 16 chromosomes. The Hi-C contact map’s lack of significant noise outside the diagonal region indicates high assembly quality ( Figure 1B ; Supplementary Figure S5 , Supplementary Table S4 ). BUSCO analysis showed 99.29% genome completeness with a high proportion of single-copy sequences ( Supplementary Figure S6 , Supplementary Table S6 ). CEGMA further assessed core gene accuracy and completeness ( Supplementary Figure S7 ). GC-depth analysis indicated a GC content of approximately 36.00% and a sequencing depth of 68.39X, affirming the absence of contamination ( Supplementary Figure S8 , Supplementary Table S7 ). These results collectively demonstrate the genome assembly’s high reliability and completeness, especially for conserved genes ( Table 1 ).

Utilizing biotool software (https://github.com/zxgsy520/biotool), we identified the presence of telomere sequences across all 16 chromosomes. Notably, telomeres were located at both the 5’ and 3’ ends of 14 of these chromosomes. The genome analyzed is characterized as a T2T-sized genome, devoid of chromosomal gaps ( Figure 1B ; Supplementary Table S8 ). The structural annotation of the genome revealed a total of 46,759 genes, with an average gene length of 2,734.57 bp, an average coding sequence (CDS) length of 1,172.84 bp, an average exon length of 253.88 bp, and an average intron length of 431.45 bp. Functional annotation of the protein-coding genes was conducted using various databases, including NR, KEGG, KOG, Swissprot, and GO ( Supplementary Figure S9 , Supplementary Table S9 ). The evaluation results from BUSCO ( Supplementary Table S10 ) indicated that 97.88% of gene elements were complete, which serves as an indirect measure of the completeness and reliability of the predicted results. Furthermore, our annotated findings included statistical data on non-coding RNAs ( Supplementary Table S11 ). The genome contains 60,840 simple sequence repeats (SSR), which represent 0.24% of the total genomic length. Scattered repeats (transposable elements, TE) constituted 38.17% of the genome, while tandem repeats (TR) and other repeat types accounted for 2.61%. Specifically, long terminal repeats (LTR) comprised 13.77%, DNA transposon factors made up 0.2%, long interspersed nuclear elements (LINE) represented 0.06%, and short interspersed nuclear elements (SINE) accounted for 0.26%. It is noteworthy that the majority of LTRs identified were of the Copia and Gypsy types, which accounted for 8.06% and 5.14%, respectively ( Supplementary Table S12 ; Table 1).

In this study, we selected 20 dicotyledonous plant species, including Hd, and constructed a phylogenetic tree utilizing 225 single-copy orthologous genes, with Dioscorea alata serving as the outgroup for comparative analysis. Our analysis yielded a total of 30,563 gene families, and by integrating the clustering information of each species’ gene family ( Supplementary Figure S10 , Supplementary Table S13 ), we determined that Hd possesses the highest number of unique gene families, totaling 14,641. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses conducted using clusterProfile v3.14.0, revealed that these genes are predominantly enriched in signaling pathways associated with vital regulatory processes, including cell apoptosis, protein digestion and absorption, glycolysis/gluconeogenesis, and RNA degradation. Notably, the enrichment of UDP glycosyltransferase activity (GO: 0008194) and terpene synthesis activity (GO: 0010333) underscores their significant role in the biosynthesis of iridoids and their glycosides in Hd ( Supplementary Figures S12A, B ). Furthermore, Hd exhibited a substantial expansion of gene families, comprising 3,348 expansion genes and 1,132 contraction genes ( Supplementary Figure S11 ). Functional annotation indicated that the expanded gene families were particularly rich in Gene Ontology (GO) terms. Specifically, the entries for chitinase activity (GO: 0004568), chitin catalytic process (GO: 0006032), and chitin binding (GO: 0008061) are all associated with the function of chitinase. Chitinase is capable of hydrolyzing chitin, a key component of the cell walls of plant pathogenic fungi, thereby exhibiting antibacterial properties and directly eliminating pathogenic fungi while also triggering plant defense mechanisms (Vaghela et al., 2022). Additionally, UDP glycosyltransferase activity (GO: 0008194) and acetyltransferase activity (GO: 0016747) facilitate the post-modification processes of iridoids ( Supplementary Figure S12C ).

The most recent whole genome duplication (WGD) event in Hd has been identified through phylogenetic tree and collinearity analyses, in conjunction with previously documented WGD events in various species ( Figure 2A ) (Yu et al., 2017; Xie et al., 2019; Wang et al., 2021; Fan et al., 2022; Wu et al., 2022; Yao et al., 2022). The core eudicotyledonous common gamma whole genome triplication (γWGT) event is estimated to have occurred approximately 150 million years ago (MYA), with a peak synonymous substitution rate (Ks) of around 2.0. The Ks curve indicates that the peak Ks value (0.4) for Hd was observed subsequent to its divergence from Coffea canephora, and later than the WGD event in Sesamum indicum, which occurred approximately 48 MYA. Furthermore, the Ks peak for Olea europaea (approximately 0.3) corresponds to a WGD event dated to about 25 MYA, which is marginally more recent than that of Hd ( Figure 2C ; Supplementary Figure S13 ). The WGD time for Hd has been calculated using the formula: T=Ks/2r, yielding an estimate of approximately 32 MYA. The presence of multiple homologous gene pairs within the genomic structure of Hd substantiates the occurrence of recent WGD events in this species ( Figure 2B ). Notably, there has been no subsequent WGD event in Coffea canephora and Coffea eugenioides following the gamma WGT event. Consequently, a classical collinearity model further elucidates the 1:2 relationship between the genomic regions of Coffea and Hd ( Figure 2D ).

In this study, we conducted extensive targeted metabolomics analyses utilizing the roots, stems, leaves, flowers, and fruits of Hd. The results of principal component analysis (PCA) indicated a high degree of consistency among the five biological replicates within each group, while also demonstrating clear differentiation among the five tissues based on the conditions of PC1 (30.6%) and PC2 (21.4%) ( Figure 3A ). We identified 30 iridoid metabolites classified as Class II monoterpenoids ( Supplementary Table S13 ) and elucidated their relative concentrations and tissue distribution through hierarchical clustering heat maps ( Figure 3B ). The findings revealed that these 30 iridoids were broadly distributed across various tissues, with 10 iridoids specifically accumulating in leaves, while only Mussaenosidic acid and Loganic acid exhibited significant accumulation in flowers. The three most abundant iridoids–Deacetyl asperuloside, Deacetyl asperulosidic acid, and Asperulosidic acid –were found to be highly concentrated in the roots, leaves, and fruits, respectively. Figure 3 illustrates the intrinsic correlations among these 30 iridoids, which can be categorized into five distinct clusters based on their correlation (a, b, c, d, e). In fact, the correlation cluster aligns with the tissue distribution of the iridoid metabolites; cluster a (leaf accumulation) is positively correlated with cluster b (flower accumulation) and cluster c (stem accumulation), with the exception of Deacetyl asperulosidic acid methyl ester ( Figures 3B, C ). Conversely, these clusters exhibit negative correlations with clusters e and f.

In this investigation, transcriptome sequencing was performed on root, leaf, flower, and fruit tissues of Hd, utilizing six biological replicates for each tissue type. The Q30 quality score for each sample exceeded 92%, with GC content varying between 41.76% and 46.23%, yielding an average of approximately 44.53% ( Table 2 ). To comprehensively elucidate the functions of differentially expressed genes, we conducted Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses by establishing a series of tissue comparison combinations: a) Root and Stem; b) Root and Leaf; c) Fruit and Root; d) Flower and Root; e) Stem and Leaf; f) Leaf and Fruit; g) Flowers and Fruit; h) Flower and Stem; i) Stem and Fruit; j) Flower and Leaf. The analysis identified the number of up-regulated and down-regulated genes within each comparison group ( Figure 4A ). The functional enrichment outcomes from the GO and KEGG analyses highlighted associations with terpenoid biosynthesis. In the GO analysis of up-regulated genes, the differentially expressed genes across each comparison group were found to be involved in monooxygenase activity, hydrolase activity (specifically hydrolyzing O-glycosyl compounds), and UDP-glycosyltransferase activity. Notably, only 20 differentially expressed genes in comparison group e were enriched for terpene synthase activity, which may correlate with the accumulation of iridoids in the leaves ( Figure 4B ). The KEGG analysis of up-regulated genes revealed that all comparison combinations were enriched in the biosynthesis of various plant secondary metabolites, with comparison group c primarily enriched in terpenoid backbone biosynthesis, while groups d, h, and j were predominantly focused on monoterpenoid biosynthesis. Additionally, several comparison groups also emphasized other terpenoid-related pathways, including sesquiterpenoid and triterpenoid biosynthesis, ubiquinone and other terpenoid-quinone biosynthesis, as well as carotenoid biosynthesis ( Figure 4C ). These findings indicate that multiple genes in Hd are involved in catalyzing the biosynthesis of monoterpenes and other terpenes through diverse pathways, with the biosynthesis of monoterpenes and their precursors primarily localized in the roots, stems, and leaves.

Previous research has demonstrated that the MVA and MEP pathways regulate the production of iridoid precursors, with terpenoids utilizing these pathways to synthesize IPP and subsequently activate the monoterpene-iridoid biosynthesis pathway via GPP (Geranyl pyrophosphate) synthesis (Bohlmann et al., 1998; Ye et al., 2019). In this study, 30 genes associated with the MVA/MEP pathways were identified, including AACT(1), HMGS(2), HMGR(3), MVK(1), PMK(3), MVD(2), and MEP pathway genes DXS(4), DXR(2), ISPE(2), ISPF(2), ISPG(2), ISPH(3), along with two IDI genes ( Figure 5A ). Heat maps revealed that these genes were predominantly expressed in roots, leaves, and flowers. In flowers, genes associated with the MVA pathway (1 HMGR, 1 HMGS, 2 PMK) and the MEP pathway (1 DXS, 1 ISPE, 1 ISPF) were highly expressed. The stem showed a predominance of the MEP pathway, including 1 DXS and 1 ISPH. Notably, the root exhibited high expression of both MVA pathway genes (1 AACT, 2 HMGR, 1 HMGS, 2 MVD, 1 MVK, 1 PMK) and MEP pathway genes (2 DXS, 2 DXR, 2 ISPG, 1 ISPH, 1 ISPE), along with 2 IDI genes. These findings suggest that the root is a crucial site for the formation of terpenoid precursors via the MVA/MEP pathways in Hd, facilitating the accumulation of synthetic raw materials. Despite the low gene enrichment in leaves and fruits, we observed high expression of 12 MEP pathway genes in leaves and 4 MVA pathway genes in fruits. These findings suggest that the MEP pathway is predominantly active in stems and leaves, whereas the MVA pathway is active in fruits ( Figure 5B ).

The iridoid biosynthetic pathway involves the cyclization process of the iridoid backbone, catalyzed by IRIS, which converts GPP into iridodial and 8-epi-iridodial. This process is regulated by the expression of 28 genes encoding key enzymes, including GES, G8/10H, 8/10HGO, and IRIS. Additionally, a total of 65 coding genes for enzymes in two distinct iridoid pathways were identified: Pathway 1, involving IO, 7DLGT, 7DLH, and LAMT, and Pathway 2, involving UGT, ALDH, F3H, and 2HFD ( Figure 5C ). A heat map of gene expression revealed that genes involved in iridoids biosynthesis are widely distributed across five tissues, with cyclization genes predominantly expressed in roots and stems. Notably, multiple G8/10H and 8/10HGO genes were enriched in all tissues, underscoring their role in regulating iridoids backbone synthesis. Additionally, Pathway 1 genes were more highly expressed in stems, whereas Pathway 2 genes were more prominent in flowers. This widespread gene enrichment contributes to the tissue-specific distribution of iridoid metabolites in Hd. However, our transcriptome data did not detect DCH in Pathway 2, suggesting that geniposidic acid in metabolites is not synthesized via Pathway 2 but may interact with Pathway 1 ( Figure 5D ).

To validate the reliability and authenticity of transcriptome data, we conducted qRT-PCR analysis to confirm the expression of key genes involved in the iridoid biosynthesis pathways across roots, stems, leaves, flowers, and fruits. We selected 13 coding genes for qRT-PCR validation, primarily from the MEP pathway, the cyclization process, and pathway1. The results indicated that most gene expression patterns aligned with the transcript data, except for Hd_27967 (ISPE), Hd_21103 (ISPH), and Hd_30653 (G8/10H) ( Figure 6A ). Specifically, MEP pathway genes Hd_44192 (DXS), Hd_28360 (DXR), and Hd_37217 (ISPG) exhibited high expression in roots. Cyclization process genes Hd_05671 (GES) and Hd_29763 (8/10HGO) were highly expressed in roots and fruits, respectively. Additionally, the two LAMT-encoding genes, Hd_09207 (LAMT) and Hd_07967 (LAMT), showed distinct expression patterns, with high expression in fruit and root, respectively ( Figure 6A ).

We selected eight iridoids, based on three with high relative content, that share similar chemical structures and biosynthetic pathways. We then examined the correlation between the relative content of these iridoid metabolites and the expression levels of enzyme-coding genes using qRT-PCR ( Figure 6B ). Significant positive correlations were found in eight pairs, while 15 pairs showed significant negative correlations. Deacetyl asperloside, the most abundant iridoids, was significantly positively correlated with the genes DXR, ISPE, ISPG, ISPH, and 2HFD (P<0.05). These genes were negatively correlated with mussaenosidic acid, loganic acid, geniposidic acid, geniposide, and deacetyl asperulosidic acid. Additionally, asperulosidic acid was significantly negatively correlated with 7DLH (P<0.05). These genes collectively regulate iridoids biosynthesis through multiple coordinated interactions ( Figure 6B ).

We investigated the post-modification processes of eight iridoid metabolites in Hd, focusing on pathways involving methyltransferase (MT), O-acetyltransferase (OAT), and Cytochrome P450s ( Figure 7A ). Loganic acid O-methyltransferase (LAMT), a key methyltransferase, typically methylates the -OH group of the iridoid C-4 carboxyl group, as corroborated by several studies (Petronikolou et al., 2018; Kang et al., 2021; Hao et al., 2023). MT requires two kinds of coding genes with varying expression patterns for the biosynthesis of Geniposide and Deacetyl asperulosidic acid methyl ester. We identified 17 genes annotated as LAMT in the transcriptome that exhibit expression patterns similar to those of Geniposide and Deacetyl asperulosidic acid methyl ester ( Figure 7B ; Supplementary Table S14 ). Additionally, we identified 20 OAT coding genes consistent with the enrichment pattern of Asperulosidic acid, the catalytic product of OAT ( Figure 7C ; Supplementary Table S14 ). The CYP450 gene family plays a crucial role in terpenoid biosynthesis, with the CYP71 family specifically implicated in monoterpenoid biosynthesis through hydroxylation of specific metabolites (Drew et al., 2013; Yuting et al., 2020). Our comprehensive analysis of the CYP71 gene cluster in Hd revealed that this family includes two key enzymes, CYP71D55 and CYP71AT96, with notable gene coding activity. Heat map analysis indicated that CYP71D55 is predominantly expressed in the root, while CYP71AT96 is expressed across various tissues ( Figure 7D ). Additionally, CYP71BE52 and CYP71BE79 exhibit high expression in the root, leaf, and fruit, respectively. Notably, despite multiple coding genes, there is no evidence supporting that CYP71AT96 is involved in hydroxylation regulation. Consequently, excluding CYP71AT96, we identified 18 genes with expression patterns similar to Deacetyl asperulosidic acid methyl ester and Deacetyl asperulosidic acid, corresponding to two enzyme types: CYP71D55 and CYP71BE52 ( Supplementary Table S14 ).

To validate the gene’s function, we integrated the screened gene with known functional enzymes (LAMT, OAT, CYP71D55, and CYP71BE52) from other plants into a bootstrap consensus gene tree using the Neighbor-Joining method (Felsenstein, 1985; Saitou and Nei, 1987). Sequence similarity analysis revealed that among the LAMT candidate coding genes, only Hd_18862 showed homology with CrLAMT (B2KPR3.1/AGX93063.1), HpLAMT (UYZ56985.1), and MsLAMT (WIE96233.1), with homologies of 75.75%, 78.21%, and 78.81%, respectively ( Figure 7E ). For OAT candidate genes, Hd_41896 clustered with CrOAT (Q9ZTK5.1), while Hd_22877, Hd_33633, and Hd_33631 clustered with CrSAT (A0A2P1GIW7.1), showing homologies of 32.76%, 43.09%, 46.31%, and 48.18%, respectively ( Figure 7F ). Within the CYP71 family, only Hd_18118 and Hd_18119 clustered with CYP71D55 (A6YIh8.7), with homologies of 61.89% and 62.11%, respectively ( Figure 7G ). These results suggest that, except for OAT, the genes exhibit high homology with known amino acid sequences, indicating potential functional consistency. We identified three genes, including one LAMT (Hd_18862) and two CYP71D55 (Hd_18118 and Hd_18119), which display distinct expression patterns ( Figure 7H ).

The Hd samples analyzed in this study were sourced from a traditional Chinese medicine cultivation site in Xinzhou City, Hubei Province. These samples were confirmed to be genuine Hd specimens. High-quality DNA was extracted from fresh plant leaves using the QIAGEN Genomic testing kit ( Figure 2A ).

We employed the HiFi Reads technology from the PacBio sequencing platform, known for its long read length and high accuracy, to enhance mutation detection, expedite assembly, and improve the continuity, accuracy, and completeness of genome assembly (Wenger et al., 2019). Initially, DNA quality was assessed using 0.75% agarose electrophoresis, Nanodrop, and Qubit methods. Subsequently, genomic DNA underwent cleavage with g-TUBE (Covaris, USA), followed by enrichment and purification of target DNA fragments using magnetic beads. DNA polymerase was utilized for repairing fragmented DNA and linking stem-loop sequencing adapters at both ends. Target fragments were then screened and purified using the gel cutter BluePippin (Sage Science, USA) to generate the library. The size of library fragments was assessed using the Agilent 2100 Bioanalyzer (Agilent Technologies, USA). Subsequently, samples were loaded onto the nanopores of the PacBio Sequel II series sequencer for real-time single-molecule sequencing to acquire raw sequencing data. Finally, the HQRF (High-Quality Region Finder) in the Smrtlink software was employed to identify the longest region of single enzyme activity and eliminate low-quality regions based on SNR. Read data exceeding 1000bp in length were selected for assembly.

Sequencing the Hi-C library directly enables the acquisition of comprehensive chromatin interaction data essential for constructing chromosome-level genomic super scaffolds (Belton et al., 2012). Fresh leaves of Hd were dissected into 2 cm segments and fixed with 2% formaldehyde for genomic DNA extraction to build a Hi-C fragment library. Quality control of the Hi-C raw data involved filtering out low-quality sequences, bridging sequences, and sequences shorter than 30 bp with a quality score <20. Utilizing bowtie2 (v2.3.2) facilitated obtaining uniquely mapped paired-end reads, followed by the identification and retention of valid interacting read pairs from these uniquely mapped reads using HiC-Pro (v2.8.1), which also eliminated invalid read pairs. Subsequently, LACHESIS was employed for further aggregation, sorting, and orientation of scaffolds onto chromosomes, employing parameters including CLUSTER, MIN-RESITES=100, CLUSTER-MAX_LINK_DENSITY=2.5, CLUSTER NONFORMATIVE RATIO=1.4, and ORDER.

The genome assembly of third-generation sequencing data was conducted using Hifiasm software to achieve high-precision genomes. Decontamination and redundancy analyses were performed based on existing genomic background data to assess genome quality. Specifically, the Smrtlink software, developed independently from the PacBio platform, was utilized to identify the longest region of single enzyme activity through HQRF. Subsequently, low-quality regions of the genome were filtered using SNR. The CCS software was employed to convert subreads into HiFi reads for quality control, followed by filtering reads longer than 1000 bp for direct assembly. The quality of the assembled high-precision genome by Hifiasm was evaluated using BUSCO, CEGMA, and genomic quality assessment methods such as sequence consistency evaluation and GC deep analysis.

Concatenated repeat sequences were identified using GMAA30 and TRF (Tandem Repeat Sequence Finder) with default parameters (Benson, 1999; Wang and Wang, 2016). Duplicate sequences were searched using MITE hunter software (Han and Wessler, 2010), LTR finder (Xu and Wang, 2007) and ltr harperst (Ellinghaus et al., 2008). RepeatMask was employed to detect repetitive sequences across the genome. Gene structures were defined by integrating de novo prediction, homologous annotation, and transcriptome prediction methods with EVM. Genes containing transposable elements were analyzed using TransposonPSI. Prediction of non-coding RNA sequences, including rRNA, snRNA, and miRNA, was conducted by comparing submissions to the Rfam database using cmscan in Inferna software. Functional annotation of protein-coding genes was performed using BLASTP and InterproScan software. Protein sequences were compared to the NR, KEGG, SwissProt, GO, and KOG databases to merge protein annotation information.

To investigate the evolutionary patterns of Hd, we identified protein sequences of dicotyledonous plants, including Hd, using Orthofinder v2.5.4 software and constructed a phylogenetic tree based on single-copy genes. Subsequently, we aligned the sequences of each single-copy gene family using MAFFT v7.20535 (Katoh et al., 2009), followed by the conversion of aligned protein sequences into codon alignments using the PAL2NAL v14 program (Suyama et al., 2006). For model selection, IQ-TREE was employed, identifying the JTT+F+R4 model as optimal. A maximum likelihood (ML) phylogenetic tree was constructed with 1000 bootstrap replicates, achieving a bootstrap support rate of 100%. Dioscorea alata was designated as the outgroup, and divergence times were estimated using the MCMCTREE package in PAML v4.9i software.

To identify whole-genome duplication (WGD) events in Hd, we initially employed diamond v2.0.337 (Buchfink et al., 2015) to align gene sequences across species and detect orthologous gene pairs. Subsequently, collinearity maps of genes were generated using WGDI software. Utilizing the blast method for protein homology alignment and MCScanX for homologous gene pair identification, we computed the synonymous substitution rate (Ks) and constructed Ks curves to analyze synonymous mutation rates within and between species. Integration of these analyses allows for inference of WGD events in the focal species.

Samples of Hd from various plant parts (roots, stems, leaves, flowers, fruits) were selected for metabolic analysis (n=6). Metabolites were detected using a comprehensive metabolomics approach based on UPLC-MS/MS on an Agilent Ultra High Performance Liquid Chromatography System (ExionLC™ AD) with a self-constructed database. The chromatographic column employed was an Agilent SB-C18 (1.8 µm * 2.1 mm * 100 mm), with a flow rate of 0.35 mL/min. The mobile phase consisted of ultrapure water and acetonitrile, both containing 0.1% formic acid, with an initial gradient elution ratio of 5% ultrapure water. Plant samples were freeze-dried, ground into powder, and 50 mg of the sample was weighed. Subsequently, 1200 μl of pre-cooled 70% methanol water internal standard extraction solution was added, followed by vortexing and centrifugation (12000 rpm, 3 min). The supernatant was aspirated, and a 0.22 μl sample solution was filtered through a microporous filter membrane for analysis.

The mass spectrometry parameters were set as follows: Electrospray Ionization (ESI) temperature was maintained at 500°C; Ion spray voltage (IS) was set at 5500 V in positive ion mode and -4500 V in negative ion mode; The ion source gases I (GSI), II (GSII), and curtain gas (CUR) were adjusted to 50, 60, and 25 psi, respectively. Collision-induced ionization settings were optimized to high levels. The QQQ scan operated in Multiple Reaction Monitoring (MRM) mode with medium collision gas (nitrogen). Declustering potential (DP) and collision energy (CE) were fine-tuned for each MRM ion pair. A specific set of MRM ion pairs corresponding to eluting metabolites was monitored in each period. Data analysis was performed using Analyst 1.6.3 software, and the relative content of each component was determined using peak area normalization. The samples’ metabolites were analyzed qualitatively and quantitatively using mass spectrometry based on a local metabolic database. Detection employed multi-response monitoring mode (MRM) to identify characteristic ions for each substance through a quadrupole mass analyzer. Signal strength (CPS) of these ions was recorded. Subsequently, the mass spectrometry data was processed using MultiQuant software for peak integration and correction. The peak area of each chromatographic peak reflected the relative content of the corresponding metabolite. Integrated data of all peak areas were saved. To compare metabolite content differences among samples, chromatographic peaks were adjusted based on metabolite retention time and peak type to ensure qualitative and quantitative accuracy.

Total RNA was extracted from roots, stems, leaves, flowers, and fruits of Hd using the HiPure Plant RNA Kit (Magen, Guangzhou, China) for RNA sequencing. Raw sequencing data were processed using Fastp software (https://github.com/OpenGene/fastp) to filter low-quality reads. The quality of the cleaned data was assessed using FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc). Clean reads were then aligned to the reference genome using HISAT2 software. Gene expression levels were quantified as FPKM (fragments per kilobase of transcript per million mapped reads). Differential expression analysis was conducted using DESeq2 software, with thresholds set at |log2(FoldChange)| ≥ 2 and adjusted p-value (Padj) ≤ 0.05.

Pairwise comparisons were performed across different tissue types, and differentially expressed genes (DEGs) were subjected to Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses, with significance determined by corrected p-values (Padj < 0.05). Data visualization was performed using the ChiPlot bioinformatics platform (https://www.chiplot.online/). In GO and KEGG analyses, “count” refers to the number of genes annotated in the respective enriched terms.

To identify the candidate genes involved in the post-biosynthetic modification processes (e.g., methylation, acetylation, hydroxylation) of iridoids, we constructed clustering heatmaps for the candidate genes of three key enzymes (LAMT, OAT, CYP71) based on their expression patterns across five tissues. A phylogenetic tree of the candidate genes was generated using MEGA11 software. The protein sequences of the candidate genes were aligned with those from other plant species (Catharanthus roseus, Hyoscyamus muticus, Hamelia patens, Mitragyna speciosa) using Clustal with default parameters. Subsequently, a phylogenetic tree was constructed using the neighbor-joining method with 1,000 bootstrap replicates. Clustering with known functional genes was performed to infer evolutionary relationships and functional similarities. Additionally, BLASTP analysis (https://blast.ncbi.nlm.nih.gov/Blast.cgi) was conducted to evaluate the sequence similarity between the candidate genes and reference enzymes. The term “Per.ident” refers to the percentage of sequence identity in alignments, which serves as a metric for assessing alignment consistency. We propose that a “Per.ident” value exceeding 60% provides evidence, to some extent, that the candidate gene may possess a function similar to that of the reference protein.

RNA quality assessment was conducted prior to extracting 1μg of total RNA from each sample. The extracted total RNA was then reverse-transcribed into cDNA using the PrimeScript™ 1st strand cDNA Synthesis Kit.Add the following components to a tube placed in an ice bath: 1 μg of total RNA, 1 μl of Oligo(dT) primer (50 μM), and 1 μl of dNTP Mix (10 mmol/L) were combined with RNase-free dH2O to a total volume of 10 μl. The mixture was then incubated at 65 °C for 5 minutes followed by rapid cooling on ice. Subsequently, add the following components to the same tube in the ice bath: 10 μl of Template RNA Primer Mixture (obtained from the previous step), 4 μl of 5×Reaction Buffer, 0.5 μl of RNase Inhibitor (40 U/μl), 1 μl of MMLV RT (200 U/μl), and 20 μl of RNase-free dH2O. The reaction mixture was then incubated at 42 °C for 30-60 minutes and terminated by heating at 95 °C for 5 minutes. For the PCR reaction setup, prepare a mixture consisting of 10 μl of 2× SYBR real-time PCR premixture, 0.4 μl of 10 μM PCR specific primer F, 0.4 μl of 10 μM PCR specific primer R, 1 μl of cDNA, and 20 μl of RNase-free dH2O. Perform an initial denaturation at 95 °C for 5 minutes, followed by 30 cycles of denaturation at 95 °C for 15 seconds, annealing at 60 °C for 30 seconds. The relative content was determined using the 2^-ΔΔCt analysis method with the formula: Ratio = 2− [CT (test)− CT (calibrator)] (test/calibrator).