Abstract Background Ferula (Apiaceae), a genus of perennial herbs with significant medicinal and culinary value, derives its characteristic odor primarily from volatile organic compounds (VOCs), particularly volatile sulfur compounds (VSCs) and terpenoids. Despite their economic importance, two critical gaps hinder targeted utilization: Ferula species-specific VOCs remains undocumented, and the metabolic pathways and key genes of VSCs and terpenoids in Ferula are still unclear. Results Volatile metabolomic and transcriptomic analyses were performed on two species of Ferula—viz., Ferula sinkiangensis K. M. Shen (XJ, strong odor) and Ferula ferulaeoides (Steud.) Korov (DS, bland odor). A total of 561 VOCs were annotated, with VSCs dominating XJ (45.98%) and terpenoids prevailing in DS (44.01%). The results VSCs were the primary source of the strong odor in XJ. In total, 42,817 differentially expressed genes (DEGs) were identified between XJ and DS in this study. KEGG pathway enrichment analysis of DEGs revealed significant divergence in the biosynthetic pathways of VSCs and terpenoids between the two Ferula species. Integrative multi-omics analysis identified 87 genes linked to VSC biosynthesis and 54 regulatory genes governing terpenoid production. In XJ, elevated accumulation of VSCs was mechanistically associated with the synergistic upregulation of glutamate-cysteine ligase (GSH), γ-glutamyltransferase (GGT), and flavin monooxygenase (FMO). High levels of VSCs precursors (e.g., S-allyl-L-cysteine sulfoxide) and transcriptional activation of these genes strongly correlated with the species-specific pungent odor profile. In DS, geranylgeranyl diphosphate synthase (GGPS) and farnesyl pyrophosphate synthase (FDPS) were identified as key regulators of terpenoid biosynthesis. Conclusions The present study expands our knowledge of the volatile metabolic profiles of the two Ferula species and lays the foundation for further studies on the regulatory mechanisms of Ferula quality. Supplementary Information The online version contains supplementary material available at 10.1186/s12870-025-06969-x. Keywords: Ferula, Volatile compounds, Volatile sulfur compounds, Terpenoids, Metabolome, Transcriptome Background Ferula L., the third-largest genus in the Apiaceae family, comprises approximately 180 species distributed across southern Europe, northern Africa, Central Asia, and surrounding regions [[38]1]. The genus exhibits significant biodiversity in China, with 31 documented species, 25 of which are found in Xinjiang [[39]2]. Historically valued since antiquity, various Ferula species have been used for both medicinal and culinary purposes since the Roman era [[40]3]. Contemporary pharmacological investigations have substantiated these species’ diverse bioactivities, including antitumor, anthelmintic, antioxidant, antimicrobial, and neuroprotective effects [[41]4–[42]8]. Beyond therapeutic applications, Ferula species are commonly utilized in cooking as vegetables, flavoring agents, or food additives [[43]8]. In traditional practices, both the roots and the resin of Ferula are used as medicinal herbs and spices, while the shoots and young leaves are used as culinary ingredients [[44]9]. This study focuses on root tissues, which contain a large number of odor-related metabolites. Ferula sinkiangensis K. M. Shen (XJ), endemic to China’s Xinjiang Uygur Autonomous Region, is a critically endangered species valued for its medicinal, culinary, and aromatic applications [[45]10–[46]12]. This species emits a distinctive pungent odor resembling alliaceous plants (particularly Allium L.), with greater olfactory persistence compared to common Allium species [[47]13]. Its endangered status is attributed to growth characteristics, overharvesting, habitat degradation, and challenging environmental conditions [[48]14]. In regional medicinal practices, Ferula ferulaeoides (Steud.) Korov. (DS) is frequently used as a substitute for medicinal species despite notable organoleptic differences: DS produces a pronounced apiaceous odor (similar to Apium graveolens) but lacks the characteristic alliaceous notes of XJ [[49]15]. Current research on Ferula species predominantly focuses on resin composition and bioactivity profiles [[50]1], with limited attention given to interspecific odor variations spanning from apiaceous to alliaceous olfactory characteristics [[51]2]. Despite the conspicuous divergence observed in Ferula species, the molecular regulatory mechanisms governing these pronounced phenotypic and volatile metabolic disparities remain poorly elucidated, representing a crucial knowledge gap in Ferula-related research. Recent research indicates that sulfur-containing compounds, monoterpenes, and other volatile terpenoids dominate the volatile oils of Ferula species [[52]15–[53]19]. Notably, organosulfur derivatives act as both chemotaxonomic markers and primary determinants of their characteristic organoleptic profiles [[54]1, [55]8]. These volatile organic compounds (VOCs) in Ferula exhibit therapeutic potential in treating respiratory disorders, such as asthma, through pulmonary excretion mechanisms, and also demonstrate efficacy in alleviating gastrointestinal dysfunctions, including flatulence and gastric erosions [[56]20]. Numerous biologically active components have been isolated from the essential oils of Ferula, showing a wide range of activities, including antimicrobial, insecticidal, antioxidant, antigerminative, antitumor, and antidiabetic properties [[57]1]. Despite progress in characterizing volatile oil composition across Ferula species, substantial interspecific variability remains undocumented. This variability is attributed to the biosynthetic plasticity of these secondary metabolites, which are influenced by cultivar, growth conditions, harvest timing, pests, and diseases. Such heterogeneity has impeded the systematic elucidation of VOCs biosynthetic pathways in the genus. In this study, two Ferula varieties, XJ and DS, were selected to investigate interspecific differences in secondary metabolites. Previous research has demonstrated significant differences in essential oils composition between these two varieties. By integrating transcriptomic and metabolomic datasets from root tissues, this study advances the mechanistic understanding of species-specific secondary metabolism, thereby facilitating both taxonomic classification and sustainable exploitation of Ferula resources. Materials and methods Plant materials The abbreviations used for the roots of the two Ferula species are as follows: Ferula sinkiangensis K. M. Shen. (XJ) and Ferula ferulaeoides (Steud.) Korov. (DS). The seeds were provided by Professor Zhu Yun of the School of Pharmacy at Shihezi University. The root samples were collected in May 2023 from two-year-old cultivated Ferula plants. The collection site was Shihu Kiln, Shihezi City, Xinjiang, China (86 °11 ′ E, 44 °15 ′ N). Sample species were identified by professor Yan Ping (College of Life Science, Shihezi University, China). Voucher specimens (No. 20230321–20230530) were preserved in the herbarium (College of Life Science, Shihezi University, China). For each species, three roots were pooled to form one biological replicate, with three independent biological replicates collected per species. Samples were immediately stored at −80 °C for subsequent metabolome and transcriptome analyses. SPME-GC-MS analysis Frozen roots of XJ and DS were ground into powder in liquid nitrogen. A 500 mg sample of the powder was promptly transferred into a 20 mL headspace vial (Agilent, Palo Alto, CA, USA), containing 20 µL of internal standard ((R)-(-)-Carvone-4,4,6,6-d4, 10 µg/mL) and 2 mL of saturated NaCl solution. Vials were sealed with crimp-top caps equipped with TFE-silicone headspace septa (Agilent). For SPME analysis, vials were incubated at 60 °C for 5 min, followed by exposure of a 120 μm DVB/CWR/PDMS fiber (Agilent) to the headspace for 15 min at 60 °C for GC–MS analysis. Desorption was performed in the injection port of an Agilent 8890 GC at 250 °C for 5 min in splitless mode. Analysis utilized an Agilent 8890 GC coupled with a 7000D MS, equipped with a DB-5MS capillary column (30 m × 0.25 mm × 0.25 μm). Helium carrier gas was set to 1.2 mL/min. The oven temperature program started at 40 °C (held for 3.5 min), ramped at 10 °C/min to 100 °C, then at 7 °C/min to 180 °C, followed by a ramp of 25 °C/min to 280 °C (held for 5 min). Mass spectra were recorded in electron impact (EI) ionization mode at 70 eV. The quadrupole mass detector, ion source, and transfer line temperatures were set to 150 °C, 230 °C, and 280 °C, respectively. A custom analyte database (MetWare Biotech. Co., Ltd, Wuhan, Hubei, China) containing retention times (RTs), qualitative and quantitative ions of the target compounds was used. For each compound, one quantitative ion and two or three qualitative ions were selected for precise scanning. All ions were detected separately in time slots according to the order of peak emergence. Compounds were identified by matching RTs to reference standards and verifying the presence of diagnostic ions after background subtraction [[58]21]. Semi-quantitative analysis of volatile compounds in Ferula was performed using the internal standard method, with (R)-(-)-Carvone-4,4,6,6-d4 (10 µg/mL) as the internal standard. The t-test was used to calculate the difference in significance of the p-value using SPSS 26.0 (IBM Corporation, USA). Orthogonal partial least squares-discriminant analysis (OPLS-DA) was performed using SIMCA software (version 14.1, Umetrics, Umeå, Sweden). In two-group analysis, differentially accumulated metabolites (DAMs) were determined by|fold change| >2 and VIP (Variable Importance in Projection) > 0.8. UPLC-MS/MS analysis Root samples were processed using vacuum freeze-drying technology, followed by grinding (30 Hz, 1.5 min) into powder using a grinder (MM 400, Retsch). A 50 mg of sample powder was weighed, and 1,200 µL of pre-cooled (−20 °C) 70% methanolic aqueous solution containing the internal standard (1 mg/L 2-chlorophenylalanine) was added for extraction. The mixture was vortexed for 30 s every 30 min, repeated 6 times. After centrifugation (11,304 ×g for 3 min), the supernatant was aspirated and filtered through a 0.22 μm microporous membrane (CNW, Anpel, Shanghai, China; cat. no. SCAA-104). The final filtrate was stored in an injection vial for UPLC-MS/MS analysis. Analyses were performed using a UPLC-ESI-MS/MS system (UPLC, ExionLC™ AD; MS/MS: SCIEX). The analytical conditions were as follows: the UPLC column: Agilent SB-C18 (1.8 μm, 2.1 mm × 100 mm); the mobile phase A, pure water with 0.1% formic acid; the mobile phase B, acetonitrile with 0.1% formic acid. The gradient program started at 95% A, 5% B, linearly changed to 5% A and 95% B over 9 min, held for 1 min, then returned to 95% A and 5% B within 1.1 min and maintained for 2.9 min. The flow rate was 0.35 mL/min, the column temperature was 40 °C, and the injection volume was 2 µL. ESI source operating parameters include a source temperature of 550 °C; ion spray voltage (IS) of 5,500 V in positive mode and − 4,500 V in negative mode; Gas I (GSI), Gas II (GSII), and curtain gas (CUR) were set at 50, 60, and 25 psi, respectively. Collision-activated dissociation (CAD) was set to high. Multiple reaction monitoring (MRM) experiments were performed with collision gas (nitrogen) set to medium. The declustering potential (DP) and collision energy (CE) for individual MRM transitions were optimized. A specific set of MRM transitions was monitored for each period based on the eluted metabolites. RNA extraction, sequencing, analysis and RT-qPCR validation RNA extraction and sequencing analysis of all Ferula samples were performed by MetWare Biotechnology Co., Ltd. (Wuhan, China). The 80 mg tissue samples were ground into powder by liquid nitrogen, and the powder were transferred into 1.5 mL of preheated 65 °C CTAB-PBIOZOL reagent (Takara, Beijing, China). The samples were incubated in Thermo mixer (Thermo mixer, Eppendorf, Hamburg, Germany) for 15 minutes at 65 °C to allow for the complete dissociation of nucleoprotein complexes. After centrifugation at 12,000 ×g for 5minutes at 4 °C, 400 µL of chloroform was added to per 1.5 mL of CTAB-PBIOZOL reagent, followed by centrifugation at 12,000 ×g for 10 minutes at 4°C. The supernatant was transferred to a new 2.0 mL tube and added 700 µL of acidic phenol and 200 µL of chloroform, followed by centrifugation at 12,000 ×g for 10 minutes at 4°C. An equal volume of chloroform was added to the aqueous phase, and the mixture was centrifuged at 12,000 ×g for 10 minutes at 4°C. The supernatant was added to an equal volume of isopropyl alcohol and placed at −20°C for 2 hours for precipitation. Afterward, the mixture was centrifuged at 12,000 ×g for 20 minutes at 4°C. The supernatant was then removed. After washing with 1 mL of 75% ethanol, the RNA pellet was air-dried in a biosafety cabinet and dissolved in 50 µL of DEPC-treated water. Subsequently, total RNA was assessed for quality and quantified using NanoDrop and Agilent 2100 bioanalyzer (Thermo Fisher Scientific, MA, USA). Then, mRNA (containing polyA) was enriched from total RNA using magnetic beads with oligo (dT). Then, cDNA libraries were constructed using the NEBNext^® UltraTM RNA Library Prep Kit (NEB, USA) and sequenced on the Illumina HiSeq platform. Transcriptome assembly was performed using Trinity, with a de novo assembly approach. Corset was used to regroup relevant transcripts into ‘gene’ clusters. The assembly was based on the reads generated from the RNA-Seq data, and no reference genome was used for genome-guided assembly ([59]https://github.com/trinityrnaseq/trinityrnaseq). Gene expression levels were estimated using RNA-Seq by Expectation-Maximization (RSEM), and the Fragments Per Kilobase of transcript per Million mapped reads (FPKM) values were calculated based on the gene length. The original data were then filtered. In this context, a fuzzy base (N) represents a base that cannot be determined during sequencing and is denoted by ‘N’. Sequencing reads were removed if reads containing N content exceeded 10% of the total base count of that read. Additionally, paired reads were removed if any sequencing read contained more than 50% low-quality bases (Q ≤ 20). All subsequent analyses were based on clean reads. The functions of unigenes were annotated through homology search against the following databases: Nr, Pfam, KOG, COG, Swiss-Prot, KEGG, and GO. Differentially expressed genes (DEGs) were identified using the DESeq2 toolkit with thresholds of log2(fold change) ≥ 2 and corrected P-value ≤ 0.05 [[60]22]. To validate the quality of the RNA-Seq data, 8 DEGs were selected for RT-qPCR analysis, including three unigenes associated with the terpenoid metabolism pathway and five unigenes involved in the biosynthesis of VSCs. The primers used are listed in Supplementary Table [61]S1. The 18S rRNA served as the reference for normalizing expression levels. The 2^−ΔΔCt method was used to calculate the relative expression levels of genes [[62]23]. Results Volatile metabolites analysis of XJ and DS Overview of VOCs profiles in XJ and DS The two Ferula varieties, XJ and DS, exhibited distinct phenotypic differences (Fig. [63]1), consistent with their divergent odor profiles. XJ had a robust onion-garlic-like odor, while DS produced a lighter, more delicate scent. To explore the factors behind these odor differences, we analyzed the VOCs profiles of XJ and DS using HS-SPME-GC/MS. A total of 561 VOCs were identified and classified into 14 categories as shown in Fig. [64]2A and Table S2. These categories include 144 terpenoids, 76 esters, 70 heterocyclic compounds, 47 hydrocarbons, 44 alcohols, 39 ketones, 37 aldehydes, 25 aromatics, 23 sulfur-containing compounds, 19 acids, 15 amines, 10 phenols, 7 nitrogen-containing compounds, and 5 others. Fig. 1. [65]Fig. 1 [66]Open in a new tab Photograph of XJ (A) and DS (B) Fig. 2. [67]Fig. 2 [68]Open in a new tab Volatile profiles of XJ and DS. A Classification of metabolites. B The relative contents and the total contents of each class. C Score scatter plot of principal component analysis (PCA). D Orthogonal partial least squares-discriminant analysis (OPLS-DA). E Volcano plots of DAMs from DS vs. XJ. F Bar chart of the top 20 metabolites with differential changes in DAMs. Red represents an upregulation in DS, and green represents a downregulation in DS The total volatile content was notably higher in DS (686.43 ± 215.55 µg/g) compared to XJ (466.73 ± 104.81 µg/g) (Fig. [69]2B). The primary VOCs in XJ and DS were volatile sulfur compounds (VSCs) and terpenoids, which accounted for 45.98% of total VOCs in XJ and 44.01% in DS. The VSC content was significantly higher in XJ (214.59 ± 132.49 µg/g) than in DS (18.83 ± 9.18 µg/g), with XJ containing approximately 11.4-fold more VSCs. Conversely, DS contained over three-fold higher total terpenoids (302.13 ± 77.82 µg/g) compared to XJ (95.63 ± 3.22 µg/g). Principal component and multivariate analysis of differential VOCs In this study, two multivariate analysis methods, Principal Component Analysis (PCA) and Orthogonal Partial Least Squares Discriminant Analysis (OPLS-DA), were used to investigate the distinct distribution patterns of VOCs [[70]24]. As shown in Fig. [71]2C-D, both PCA and OPLS-DA effectively distinguished XJ and DS. PCA analysis revealed that PC1 and PC2 explained 64.98% and 14.99% of the variability, respectively. The OPLS-DA model demonstrated excellent predictive capability with R²Y = 1 and Q² = 0.988, indicating reliable model validation. Using thresholds of|fold change| >2, VIP > 0.8, and P < 0.05, we identified 490 DAMs between DS and XJ, including 329 up-regulated and 161 down-regulated (Fig. [72]2E). To further identify the key VOCs, we analyzed the top 20 most significantly altered VOCs (Fig. [73]2F). In XJ, the concentration of (E)-sec-butyl propenyl disulfide, accounting for 43.07% of the total volatile compounds, was 201.06 ± 136.12 µg/g, showing 61.91-fold higher than in DS. Among the 14 most downregulated compounds in XJ, three were terpenoids. In DS, limonene was the most abundant volatile compound, with a relative content of 40.6 ± 17.07 µg/g, which is 410.32-fold higher than in XJ. Terpinolene and α-pinene also showed marked enrichment in DS, with 268.74-fold and 160.06-fold higher than in XJ, respectively. These results underscore stark contrasts in volatile profiles between the two varieties. Non-volatile metabolites analysis of XJ and DS Analysis of non-volatile metabolites associated with VSCs metabolism To elucidate the metabolic divergences underlying odor variation, we conducted further analysis using UPLC-MS/MS to compare key sulfur-containing precursors and amino acids in XJ and DS (Table S3). A total of 1,302 non-volatile metabolites were analyzed, with 22 directly linked to the VSCs metabolic pathway. Correlation analysis revealed that 9 non-volatile metabolites were significantly correlated with (E)-sec-butyl propenyl disulfide (|r| >0.8, p-value < 0.05), with 7 showing positive correlations and 2 showing negative correlations (Table S4). Targeted metabolomics revealed that S-allyl-L-cysteine sulfoxide (alliin), a precursor to (E)-sec-butyl propenyl disulfide, exhibited 241.46-fold higher abundance in XJ compared to DS. Additionally, S-methyl-L-cysteine was also present at higher levels in XJ and showed a significant positive correlation with (E)-sec-butyl propenyl disulfide. Analysis of non-volatile terpene metabolites In this study, 71 terpenoids were detected using UPLC-MS/MS (Table S3). Among them, 30-norhederagenin, nootkatone, and roburic acid had the highest relative content in DS. Progesterone had the highest abundance in XJ. DAMs were selected according to the following thresholds: fold change ≥ 2 or ≤ 0.5 and VIP > 0.8. A total of 36 differential terpenoid metabolites were identified, including 15 up-regulated and 21 down-regulated in XJ relative to DS. Among the five terpenoids with the highest VIP values, three (nootkatone, 7-Deoxyloganic acid, and roburic acid) were upregulated in DS, while two (dehydrovomifoliol and verminoside) were downregulated in DS. Transcriptomic analysis of XJ and DS To investigate intracellular metabolic variations between XJ and DS, RNA sequencing and differential gene expression analysis were performed. The sequencing of six samples yielded high-quality clean data ranging from 6.32 Gb to 6.92 Gb per sample with minimal read loss (0.02%) and consistently high sequencing accuracy (Q30 scores: 94.06-94.34%). The GC content for each sample ranged from 42.89 to 43.01%. PCA of the transcriptome expression profiles revealed distinct clustering patterns, with the three biological replicates of each sample clustering tightly. This clear separation between XJ and DS indicates substantial inter-species transcriptional differences (Fig. [74]3A). Among the 180,744 detected genes, 42,817 genes were identified as DEGs based on the criteria of fold change ≥ 2 or ≤ 0.5 and a false discovery rate (FDR) < 0.05. The DEGs comprised 21,503 up-regulated and 21,314 down-regulated genes in XJ compared to DS (Fig. [75]3B). Fig. 3. [76]Fig. 3 [77]Open in a new tab The features of transcriptomic dataset compared between XJ and DS. A Principal component analysis (PCA) of all gene expressions between XJ and DS. B Volcano plots of DEGs from DS vs. XJ. C KEGG enrichment of DEGs in XJ and DS. The abscissa is rich factor, the ordinate is the KEGG pathway name, the number of DEGs enriched in this pathway in the size of the circle, and the p-value is marked in the color. D The numbers of DEGs in the two key VOCs formation. E GO annotations maps of DEGs in XJ and DS. The X-axis represents GO Term and the Y-axis represents the number of genes, MF: molecular functions, BP: biological processes, CC: cellular components Functional annotation through Gene Ontology (GO) analysis revealed enrichment in three categories (Fig. [78]3E). In the biological process category, the most enriched subcategories included cellular processes, metabolic processes, responses to stimuli, and biological regulation. In the cellular component category, the predominant subcategories were cellular anatomical entities and protein-containing complexes. Within the molecular function category, the most prominent subcategories included binding, catalytic activity, and transcription regulator activity. KEGG pathway analysis further identified 145 pathways, among which 19 were significantly enriched metabolic pathways (P-value < 0.05) (Fig. [79]3C, Table S5). Among the the pathways enriched to, 5 were related to the metabolism of sulfur-containing compounds, while 5 were associated with terpenoid biosynthesis (Fig. [80]3D). The most enriched pathway was cysteine and methionine metabolism (ko 00270), followed by glutathione metabolism (ko 00480) and terpenoid backbone biosynthesis (ko 00909). Metabolic changes of the key VOCs biosynthetic pathways Analysis of the VSCs biosynthetic pathways Due to the instability of most VSCs precursors, the biosynthetic pathways of VSCs in Ferula species remain largely unexplored. To explore potential differences in odor formation, we compared the intracellular metabolism of XJ and DS using UPLC-MS/MS and RNA sequencing analysis. Based on these findings, we proposed the biosynthetic pathways of VSCs in XJ and DS. Through KEGG enrichment analysis, a total of 22 non-volatile metabolites and 87 genes were identified and associated with five metabolic pathways (Fig. [81]4A, Table S6). Fig. 4. [82]Fig. 4 [83]Open in a new tab Profiles of non-volatile intermediates and the genes encoding enzymes involved in the volatile sulfur compounds biosynthesis. A Potential biosynthetic pathways of VSCs in XJ and DS. B Heatmap of the FKPM values of each identified genes. The solid and dotted arrows represented the one-step and multi-step reactions, respectively. Boxed metabolites indicate detection by UPLC-MS/MS or GC-MS. metabolites with rectangular dark yellow shadows indicate detection in UPLC-MS/MS, with red arrows indicating upregulation in XJ and green arrows indicating downregulation in XJ. The metabolite with blue shading indicates detection in SPME-GC-MS and represents downregulation in XJ, the metabolite with pink shading indicates detection in SPME-GC-MS and represents upregulation in XJ, and the dashed box represents other possible pathways Among the 22 non-volatile metabolites analyzed,11 showed significant changes. Of these, 6 were up-regulated in XJ, and 5 were down-regulated. Alliin and methiin precursors, key compounds contributing to the characteristic onion and garlic odor, were identified in both XJ and DS. The content of alliin and methiin-precursor (S-methyl cysteine) was 241.46-fold and 4.41-fold higher in XJ (Table S3). However, no isoalliin-related compounds were identified. This suggests that lower precursor metabolites in DS may lead to the reduced accumulation of onion and garlic odor compounds. We analyzed the biosynthesis pathway of VSCs by studying 87 unigenes encoding 24 enzymes involved in this pathway. As shown in Fig. [84]4B, nearly all DEGs were upregulated in XJ, except for genes involved in the methionine-homocysteine cycle, including cystathionine gamma-synthase (METB), homocysteine methylase (METE), DNA 5-cytosine methylase (DNMT), adenosylhomocysteinase (ACHY), nicotianamine synthase (NAS), and aminocyclopropanecarboxylate synthase (ACS). This differential expression pattern suggests that VSCs biosynthesis in XJ is uniquely regulated, potentially driving distinct metabolite production. Analysis of the terpenoids biosynthesis pathway To investigate the regulatory mechanisms underlying terpenoid accumulation patterns in XJ and DS, we conducted comparative expression analyses of genes involved in terpenoid biosynthesis. Terpenoid synthesis in plants occurs through two distinct biosynthetic pathways: the mevalonate (MVA) pathway and the 2-C-methyl-D-erythritol-4-phosphate (MEP) pathway [[85]25–[86]28]. Significant differences between the two pathways, with the MEP pathway primarily generating monoterpenes and diterpenes, while the MVA pathway predominantly produces sesquiterpenes and triterpenoids [[87]29, [88]30]. Our transcriptomic profiling identified 54 unigenes associated with terpenoid biosynthesis pathways (Fig. [89]5A, Table S6). Fig. 5. [90]Fig. 5 [91]Open in a new tab Profiles of the genes encoding enzymes involved in the terpenoids biosynthesis. A Biosynthetic pathways of terpenoids in XJ and DS. B Heatmap of the FKPM values of each identified genes. The arrows represented the one-step reaction. The metabolites with blue shading indicate detection in SPME-GC-MS and UPLC-MS/MS, where the red arrow indicates upregulation in XJ, the green arrow indicates downregulation in XJ, and the dashed box represents other possible pathways In the MVA pathway, 16 unigenes were predicted in Ferula, with 14 DEGs showing higher expression in DS. In the MEP pathway, 16 unigenes were predicted. Among the first four steps of the MEP pathway, the following enzymes showed higher expression in XJ: four 1-deoxy-D-xylulose-5-phosphate synthase (DXS), one 1-deoxy-D-xylulose-5-phosphate reductoisomerase (DXR), three 2-C-methyl-D-erythritol 4-phosphate cytidylyltransferase (ISPD) and two 4-diphosphocytidyl-2-C-methyl-D-erythritol kinase (ISPE). In contrast, five unigenes involved in the late stages of the MEP pathway—2-C-methyl-D-erythritol 2,4-cyclic diphosphate synthase (ISPF), two 4-hydroxy-3-methylbut-2-en-1-yl diphosphate synthases (GCPE), and two 4-hydroxy-3-methylbut-2-en-1-yl diphosphate reductases (ISPH)—exhibited increased expression in DS. The biosynthesis of terpenoid precursors involves critical enzymatic conversions. Geranylgeranyl diphosphate synthase (GGPS) and farnesyl diphosphate synthase (FDPS) mediate the formation of geranyl diphosphate (GPP) and farnesyl diphosphate (FPP) from dimethylallyl pyrophosphate (DMAPP) and isopentenyl pyrophosphate (IPP), respectively. The isopentenyl diphosphate delta isomerase (IDI) gene regulates the rate of interconversion between IPP and DMAPP. We identified two GGPS, five FDPS, and three IDI genes, with nine of these showing high expression in DS. Characterization of terpenoid synthase (TPS) revealed twelve unigenes distributed across four functional categories, including three monoterpene synthases (MTPS), three diterpene synthases (DTPS), three sesquiterpene synthases (STPS), and three triterpene synthases (TTPS). Notably, DTPS exhibited highly expressed levels in XJ, while the other synthase categories were more highly expressed in DS. Overall, our results indicate that the high expression levels for the enzymes in the terpenoid backbone pathway correlate with a high rate of terpenoid synthesis and the different expression levels of terpenoid synthase genes in XJ and DS may contribute to their different terpenoid compositions. Correlation between key volatile differential metabolites and genes In this study, we selected 87 DEGs related to key differential volatile metabolites, specifically (E)-sec-butyl propenyl disulfide, along with non-volatile metabolites and associated metabolic pathways, for Pearson correlation analysis. By calculating the correlation coefficients between genes and metabolites, we found that 54 genes exhibited significant correlations with (E)-sec-butyl propenyl disulfide (|r| >0.8, P < 0.05). Among these genes, 35 showed positive correlations and 19 showed negative correlations. Additionally, 33 genes exhibited weak correlations with (E)-sec-butyl propenyl disulfide, which did not reach statistical significance (Fig. [92]6A, Table S7). Fig. 6. [93]Fig. 6 [94]Open in a new tab Correlation networks of key volatile compounds and associated genes. A Correlation network of DEGs with (E)-sec-butyl propenyl disulfide and (B) correlation network of DEGs with key terpenoids. In the figure, |r| >0.8 and p < 0.05. Genes are represented as triangles, with red lines indicating positive correlations and blue lines indicating negative correlations Further analysis of key terpenoids (limonene, terpinolene, and α-pinene) and 54 DEGs involved in the terpenoid biosynthesis pathway revealed a gene-metabolite co-network comprising 41 gene nodes, 3 metabolite nodes, and 109 related edges (Fig. [95]6B, Table S8). RT-qPCR validation of DEGs from the RNA-Seq analysis To validate the expression patterns of key genes involved in terpenoids and VSCs biosynthesis, we performed RT-qPCR on eight unigenes selected from RNA-Seq data: DXS-4 (Cluster-18746.7), GGPS-1 (Cluster-74530.2), TTPS-3 (Cluster-39771.0), cysK-3 (Cluster-63830.0), cysE-2 (Cluster-75231.0), GGT-3 (Cluster-58909.3), GST-4 (Cluster-1353.7), and FMO-3 (Cluster-68501.0). The expression levels of these selected genes from RT-qPCR analyses were generally consistent with those deduced from their FPKM data from RNA-Seq (Fig. [96]7A-H). The correlation between the RT-qPCR and RNA-Seq measurements was evaluated and the coefficient of determination (or R-squared) was 0.829 (Fig. [97]7I). Fig. 7. [98]Fig. 7 [99]Open in a new tab A-H Expression patterns of candidate DEGs was assessed using RT-qPCR. I Gene expression correlation between RT-qPCR and FPKM values obtained from RNA-Seq data Discussion Analysis of VOCs in XJ and DS Previous studies have established that nearly all Ferula species including F. lehmannii [[100]31], F. assa-foetida [[101]32], and F. jaeschkeana [[102]33], typically produce VOCs dominated by terpenoids and sulfur-containing compounds. These VOC profiles are critical for the medicinal properties and culinary use of Ferula [[103]33]. However, previous research has primarily focused on extracting and analyzing volatile oils from Ferula [[104]34], and technological limitations have hindered comprehensive VOCs characterization and comparisons [[105]2]. To address these gaps, GC-MS detected 561 VOCs from two Ferula species, far exceeding previous reports [[106]18]. VSCs and terpenoids were identified as the predominant VOCs. VSCs typically exhibit intense sulfurous odors, whereas terpenoids produce a milder fragrance, often described as floral, citrusy, or herbal [[107]9, [108]35]. In XJ, (E)-sec-butyl propenyl disulfide dominated the total volatile profile, consistent with previous reports (10.38–64.47%) [[109]17]. This metabolite contributes to the characteristic pungency of Ferula [[110]36] and is key to the distinct odor of XJ. The contribution of VSCs and terpenoids to Ferula odor remains unclear, despite the compositional differences. Analysis of the VSCs biosynthetic pathway (Z)-sec-butyl propenyl disulfide and (E)-sec-butyl propenyl disulfide have been identified as representative compounds of Ferula [[111]31]. Similar derivatives of 1-propenyl sulfide have also been identified in Allium species [[112]36–[113]38]. Allium species synthesize specialized sulfur-containing metabolites derived from cysteine, notably S-alk(en)yl cysteine sulfoxides (CSOs). Tissue damage in Allium species activates the enzyme alliinase (EC 4.4.1.4), which catalyzes the cleavage of CSOs to generate sulfenic acids that convert into VSCs, contributing to Allium odor and bioactivity [[114]39, [115]40]. We hypothesize that Ferula forms these VSCs via a similar pathway to Allium. Studies suggest that the levels of these organic sulfur compounds depend on the regulation of the sulfur assimilation pathway genes [[116]41, [117]42]. Plants assimilate inorganic sulfur as sulfate, which is activated by ATP sulfurylase (PAPSS), converting it to adenosine 5’-phosphosulfate (APS). Subsequently, APS is converted to sulfite and ultimately to sulfide (H[2]S) through sulfite reductase (sir), and sulfide is integrated into O-acetylserine to form cysteine [[118]43]. Previous studies have characterized PAPSS gene expression in response to sulfur supply in onions [[119]44, [120]45]. Building on this, both XJ and DS showed expression of genes PAPSS, ADP-sulfurylase (HINT4), sir and APS kinase (cysC) were all expressed, indicating active sulfate reduction. However, higher expressed levels in XJ suggest enhanced sulfate uptake efficiency in its root tissues, consistent with studies linking increased sulfur supply to higher PAPSS and sir activity in onions [[121]44]. In H[2]S assimilation pathways, the cystathionine synthetase (metB) genes exhibited high expression levels in DS, while cysteine synthase (cysK) genes showed higher expression in XJ, with a 17.45-fold higher FPKM value compared to DS. Notably, the relative substrate (O-acetyl serine) was also elevated in XJ. These results indicate that H[2]S is preferentially channeled into cysteine synthesis via the serine pathway rather than the aspartate pathway in XJ. However, DS may follow the opposite trend (Fig. [122]4B). In XJ, higher expression of S-adenosylmethionine synthetase (METK) led to greater abundance of its product, S-(5’-Adenosyl)-L-methionine, and increased expression of downstream genes such as spermidine synthase (speE) and S-adenosylmethionine decarboxylase (AMD1). Conversely, lower expression levels of DNMT and AHCY suggest that more metabolic flux flowed out of the methionine-homocysteine cycle in XJ (Fig. [123]4B). In Allium, cysteine is converted to glutathione and degraded into CSOs, such as alliin or isoalliin [[124]37]. The odor of Allium arises from the enzymatic hydrolysis of the CSOs by alliinase (EC 4.4.1.4), generating VSCs along with by-products pyruvic acid and NH[3] [[125]46, [126]47]. Non-volatile metabolites analyzed in XJ and DS are similar to those in Allium. Significant differences were observed in the expression levels of CSOs biosynthesis and hydrolysis pathways between XJ and DS. Specifically, genes encoding glutamate-cysteine ligase (GSH), glutathione synthetase (GSS), glutathione S-transferase (GST), glutamate carboxypeptidase (CPG), γ-glutamyltransferase (GGT), and flavin monooxygenase (FMO) wereexpressed at a higher level in XJ than in DS, except for FMO-1. To date, only two enzymes have been identified in the glutathione-to-CSOs pathway: GGT (catalyzing the removal of glycyl group) and FMO (catalyzing S-oxygenate) [[127]48, [128]49]. These genes have been found in the Allium sativum L [[129]47], Allium hookeri [[130]50], Brassica oleracea [[131]51], and Toona sinensis [[132]35]. CSOs can convert into VSCs (e.g., 1-propenyl sulfenic acid) via alliinase (ALL) or spontaneous reactions [[133]47, [134]52]. Notably, XJ exhibited elevated levels of the key precursor alliin. Thus, the higher VSC content in XJ likely stems from both increased abundance of VSC precursors (e.g., S-methyl-L-cysteine, alliin) and enhanced expression of CSO pathway enzymes. Analysis of the terpenoids biosynthesis pathway Recent transcriptomic advances have identified regulatory genes involved in terpenoid biosynthesis [[135]53, [136]54]. GGPS catalyzes the condensation of three IPP and one DMAPP molecule into geranylgeranyl diphosphate (GGPP), upregulated in DS. GGPP is a pivotal precursor for monoterpenes and diterpenoids [[137]55]. This upregulation aligns with the observed increase in monoterpene levels in DS. Additionally, the function of key TPS genes in terpenoid biosynthesis has been elucidated [[138]56]. In Ferula, TPS genes displayed differential expression patterns, MTPS, such as α-terpineol synthase and limonene synthase, were up-regulated in DS. DTPS, including momilactone-A synthase, were down-regulated in DS. These results are consistent with the higher concentration of diterpenes found in XJ (Fig. [139]5B). Terpenoids are the predominant components in Apiaceae essential oils, contributing to herbal and woody odors [[140]57]. Our analysis revealed significant upregulation in DS of sesquiterpenes including guaiol (guaiac wood rose), alpha-bisabolene (balsamic spicy floral), and beta-bisabolene (balsamic woody). Sesquiterpene biosynthesis involves two IPP molecules and one DMAPP molecule. In this study, high expression of ACTA, HMGCS, MAK and MVAK suggests that biosynthesis of IPP and DMAPP via the MVA pathway is enhanced in DS [[141]58]. The concurrent upregulation of farnesyl diphosphate synthase (FDPS) in DS further supports sesquiterpenoid biosynthesis. FDPS catalyzes the condensation of two IPP molecules and a DMAPP molecule to form farnesyl diphosphate, which can be converted to sesquiterpenoids by cytoplasmic sesquiterpene synthase [[142]59]. The significant upregulation of ACTA, HMGCS and HMGCR in DS was consistent with the observed increase in sesquiterpenoid levels, which are known to contribute to a more pronounced aromatic odor (Fig. [143]5B). Terpenoid synthesis is controlled by key enzymes and multiple reaction [[144]60]. Previous studies revealed that the expression profiles of genes encoding key enzymes such as GGPS, DXS, MVK, and TPS, were positively correlated with terpenoid yield [[145]61, [146]62]. Consistent with these earlier findings, our data suggest that upregulated upstream genes (DXSs, ACTAs, HMGCSs, HMGCRs, MVK, MVAKs, IDIs and GGPSs) in the DS might enhance the metabolic flux toward the precursors of terpenoid biosynthesis. Conclusion In summary, this study investigated odor profiles and formation mechanisms of XJ and DS. XJ is characterized by VSCs as major volatiles and key biomarkers, while DS predominantly features terpenoids as main VOCs. The accumulation of VSCs in XJ may derive from cysteine synthesis and the availability of key precursor substances, with S-allyl-L-cysteine sulfoxide emerging as the key precursor for VSCs biosynthesis. Transcriptomic analyses indicated that high expression of genes related to the MVA and MEP pathways led to the enhanced synthesis of terpenoids in DS. Collectively, these findings provide a theoretical basis for the characterization of specific odor associated with Ferula. Supplementary Information [147]Supplementary Material 1.^ (148.8KB, docx) Acknowledgements