Integrated Metabolomic and Transcriptomic Analyses Reveal the Coumarin Biosynthesis Pathway and Key Regulatory Genes in the Pericarp of Zanthoxylum

TL;DR

Researchers combined metabolomic and transcriptomic data to identify genes involved in coumarin production in two types of Sichuan pepper.

Contribution

The study identifies key regulatory genes and metabolic pathways responsible for coumarin variation in Zanthoxylum pericarp.

Findings

24 coumarins were identified, with most being more abundant in green Sichuan pepper.

56 candidate genes were found to be strongly correlated with scopoletin and scopolin accumulation.

Overexpression of CCoAOMT, COMT, and F6H increased scopoletin and scopolin levels in transient assays.

Abstract

Coumarins in the pericarp of Zanthoxylum contribute to the characteristic numbing–aromatic flavor and are associated with diverse bioactivities. To characterize coumarin divergence between two Zanthoxylum materials, mature pericarps of Dahongpao Z. bungeanum (red Sichuan pepper) and Z. planispinum var. dingtanensis (green Sichuan pepper) were analyzed by widely targeted UPLC–ESI–MS/MS metabolomics integrated with transcriptome sequencing. This approach enabled joint profiling of metabolites and transcripts to identify genes associated with material-specific coumarin accumulation. Across the two materials, 583 metabolites were detected, with flavonoids, phenolic acids, and alkaloids as the predominant classes. Among these, 24 coumarins were identified, and most showed significantly higher abundance in green Sichuan pepper than in red Sichuan pepper. Pathway enrichment analysis indicated that differentially accumulated coumarins were mainly associated with the phenylpropanoid biosynthesis pathway, consistent with coordinated metabolic and transcriptional regulation. The integration of metabolite abundance with gene expression patterns identified 56 candidate genes strongly correlated with scopoletin and scopolin accumulation. To evaluate functional relevance, CCoAOMT, COMT, and F6H were cloned and transiently overexpressed in Nicotiana benthamiana. Transient expression assays showed that overexpression of each gene increased scopoletin and scopolin, supporting their involvement in coumarin biosynthesis. Collectively, these results clarify molecular determinants of coumarin variation between the two materials and highlight candidate genes for quality improvement and metabolic engineering.

Figures

Full text

1. Introduction

Sichuan pepper broadly refers to several economically important Zanthoxylum species and varieties used in China as spices and medicinal resources, including Zanthoxylum bungeanum, Zanthoxylum schinifolium, and Zanthoxylum planispinum var. dingtanensis, which belongs to the genus Zanthoxylum in the family Rutaceae. These species are widely used as both flavoring agents and traditional medicinal materials. Previous phytochemical investigations have shown that the pericarp of Sichuan pepper contains a wide array of secondary metabolites, including essential oils, amides, alkaloids, flavonoids, lignans and coumarins [1,2,3,4]. Multiple simple coumarins and furanocoumarins, including 7-oxo-7H-furo [3,2-g]chromen-9-dimethyl carbamate, have been isolated from the pericarp of red-fruited Sichuan pepper [5]. More than 70 coumarins, predominantly simple coumarins with a smaller proportion of furanocoumarins, have been identified in green Sichuan pepper (Z. schinifolium), including 7-methoxycoumarin, isoscopoletin and fraxinol [3]. Coumarins concentrated in the pericarp of Sichuan pepper are thought to contribute to its aromatic profile and may contribute to its pharmacological effects [2]. In particular, simple coumarins such as scopoletin and its glucoside scopolin are frequently detected in Zanthoxylum pericarps and are considered representative metabolites that connect the phenylpropanoid metabolism to pericarp aroma and bioactivity. In addition to simple coumarins, furanocoumarins (e.g., bergapten- and xanthotoxin-type compounds) have also been reported in some Zanthoxylum materials, highlighting the chemical diversity of coumarins in Sichuan pepper pericarp [3,5]. Simple coumarins refer to compounds characterized by the 1,2-benzopyran-2-one core scaffold, which are further substituted on the benzene ring with groups such as hydroxyl, methoxy, methylenedioxy, or prenyl. Importantly, these compounds do not contain any additional fused heterocyclic systems, meaning they lack fused furan or pyran rings. In contrast, coumarins that have a fused furan ring on the coumarin scaffold are designated as ‘furanocoumarins’, while those with a fused pyran ring are referred to as ‘pyranocoumarins’. This structure-based classification is extensively utilized in the study and review of plant coumarins, providing a clear framework for distinguishing between coumarin subclasses in Sichuan pepper and elucidating their potential biosynthetic origins [3,5].

Coumarins extracted from Sichuan pepper exhibit a broad spectrum of pharmacological activities, including anti-inflammatory, antioxidant, analgesic, anticancer, antithrombotic, cardioprotective and neuroregulatory effects. They constitute a major chemical basis for the pharmacological actions of Sichuan pepper [3]. In addition to cytotoxic and apoptosis-inducing effects in tumor cells, coumarin and its derivatives show antimicrobial, antispasmodic, and anti-platelet (antiaggregation) activities. They have also been reported to exert anti-rheumatic effects in arthritis models [6,7]. Beyond the volatile oils and amide compounds that dictate the characteristic numbing and aromatic flavor of Sichuan pepper, coumarins also contribute to its distinctive taste profile. Furthermore, increasing attention has been paid to the potential of coumarin and related phenylpropanoid derivatives from Sichuan pepper in the management of immunological, digestive, and cardiovascular problems [1,6].

Plant coumarins are generally derived from the phenylpropanoid pathway. In the early steps of this pathway, key enzymes such as phenylalanine ammonia-lyase (PAL), cinnamate 4-hydroxylase (C4H), and 4-coumarate: CoA ligase (4CL) convert phenylalanine into hydroxycinnamoyl-CoA derivatives. Subsequently, ortho-hydroxylation, cis-trans isomerization, and lactonization generate the characteristic 2H-1-benzopyran-2-one skeleton [8,9]. Among these, the primary rate-limiting enzyme in simple coumarin biosynthesis is feruloyl-CoA 6′-hydroxylase (F6′H), a dioxygenase that is dependent on ferredoxin. Its structure and function have been confirmed in species such as Clematis, Manihot esculenta, and Arabidopsis thaliana. Coumarin accumulation is often positively associated with F6′H expression levels [2,9,10,11]. Nevertheless, compared with these model plants, research on the coumarin biosynthesis and its key structural and regulatory genes in Sichuan pepper pericarp remains limited. Existing work has largely focused on isolation/identification of individual coumarins and bioactivity evaluation, whereas systematic analyses of coumarin metabolic networks across different Sichuan pepper pericarps are still scarce [4,12].

Zanthoxylum planispinum var. dingtanensis is a green Sichuan pepper variety widely distributed in the karst plateau canyon regions of Guizhou Province and is characterized by rapid growth and strong stress tolerance. Dahongpao, an important red Sichuan pepper cultivar of Z. bungeanum, is widely cultivated in China. These two cultivars show significant differences in pericarp aroma and numbing intensity. To elucidate coumarin-related biosynthetic pathways in Sichuan pepper pericarp and to prioritize key candidate genes, we performed integrated transcriptomic and metabolomic analyses of mature pericarp (fresh weight) from Z. planispinum var. dingtanensis and Dahongpao Z. bungeanum. Heterologous overexpression was employed to validate the functions of selected candidate genes. This study aims to provide a theoretical basis for targeted genetic improvement of Sichuan pepper and to elucidate the molecular mechanisms underlying its aromatic and therapeutic properties.

2. Materials and Methods

2.1. Plant Materials

Dahongpao Zanthoxylum bungeanum (red Sichuan pepper) and Zanthoxylum planispinum var. dingtanensis (green Sichuan pepper) were used as experimental materials (Figure 1A). One-year-old seedlings of Z. bungeanum were sourced from Hanyuan County, Sichuan Province, China, whereas one-year-old seedlings of Z. planispinum var. dingtanensis were obtained from Zhenfeng County, Qianxinan Buyi and Miao nationality Autonomous Prefecture. Both Zanthoxylum species were planted in 2019 at the same cultivation base in Zhenfeng County, Guizhou Province, China (105°38′36″ E, 25°24′10″ N). In 2023, five-year-old plants with vigorous growth and comparable field conditions (management practices, light exposure, and soil moisture) were selected from this site for subsequent analyses. Nine plants per species were sampled; three plants were pooled as one biological replicate, yielding three biological replicates in total. During sampling, fruits with uniform ripeness, similar color and shape, and no mechanical damage were collected from four directions of each plant (40 fruits per plant; 120 fruits per replicate). After rinsing with sterile water and air-drying, the fruits were promptly dehusked. The pericarps were immediately transferred into prechilled 50 mL centrifuge tubes, snap-frozen in liquid nitrogen, transported to the laboratory on dry ice, and stored at −80 °C until further analysis. Tissue-cultured tobacco seedlings were provided by the Mianyang Academy of Agricultural Sciences (Mianyang, China).

2.2. Reagents

Scopolin (CAS No. 531-44-2, purity ≥ 98%) and scopoletin (CAS No. 92-61-5, purity > 98%) standards were purchased from Shanghai Yuanye Biotechnology Co., Ltd. (Shanghai, China). These standards were used for metabolite identification and quantification, including the preparation of calibration curves, as described in the corresponding analytical procedures. In addition, the pB1121 vector, Escherichia coli strain DH5α, and Agrobacterium tumefaciens strain GV3101 were obtained from BioMADE (St. Paul, MN, USA)

2.3. Metabolite Extraction and UPLC–MS/MS Analysis

Pericarp samples were removed from −80 °C storage and kept on ice prior to extraction. The biological samples were freeze-dried using a vacuum freeze-dryer (SCIENTZ-100F (Ningbo Scientz Biotechnology Co., Ltd., Ningbo, China)) and then ground to a fine powder with a mixer mill (MM400, Retsch (Haan, Germany) at 30 Hz for 1.5 min under liquid nitrogen. Approximately 100 mg of powder was transferred to a 2 mL microcentrifuge tube and extracted with 1.2 mL of pre-chilled 70% (v/v) methanol. The mixture was vortexed thoroughly and extracted by shaking at 4 °C. After incubation, the samples were centrifuged at 12,000 rpm for 10 min, and the supernatant was filtered through a 0.22 μm organic membrane filter prior to UPLC–MS/MS analysis. Scopolin standard (CAS No. 531-44-2, purity ≥ 98%; Shanghai Yuanye Bio (Shanghai, China)) was used for instrument calibration and metabolite quantification/verification, as appropriate. Metabolite profiling of the two Sichuan pepper pericarp samples was performed using UPLC–MS/MS as described previously [13,14]. Metabolite detection and date acquisition were conducted by Maiwei Metabolic Biotechnology Co., Ltd. (Wuhan, China).

2.4. Metabolome Data Analysis

Following identification, peak alignment, and normalization of the raw mass spectrometry data, principal component analysis (PCA) and orthogonal partial least squares discriminant analysis (OPLS-DA) were used to assess metabolic differences and samples clustering. Differential metabolites (DMs) were identified using both of the following criteria: fold change (FC) ≥ 2.0 or |log2FC| ≤ 1 and a false discovery rate (FDR)-adjusted p value (q value) < 0.05 [15]. The KEGG database was utilized to annotate DMs. Metabolite set enrichment analysis was performed using MetaboAnalyst (MSEA module), and pathways with FDR-adjusted p values < 0.05 were considered significantly enriched [15].

2.5. RNA Extraction and Transcriptome Sequencing

Total RNA was extracted from Sichuan pepper pericarp tissues using standard plant RNA isolation procedures. The integrity of the RNA and potential genomic DNA contamination were assessed by agarose gel electrophoresis. RNA purity was evaluated using a NanoPhotometer spectrophotometer (OD260/280 and OD260/230) (Implen, Munich, Germany), and RNA concentration was quantified using a Qubit 2.0 fluorometer (Thermo Fisher Scientific, Waltham, MA, USA). RNA integrity was further evaluated using an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). mRNA was enriched based on the poly(A) tail using Oligo(dT) beads, then fragmented with fragmentation buffer, and reverse-transcribed into first-strand cDNA using random hexamer primers, followed by second-strand cDNA synthesis [15]. Purified double-stranded cDNA underwent end repair, A-tailing, and adapters ligation, followed by purification and size selection using AMPure XP (Beckman Coulter, Brea, CA, USA) beads and PCR enrichment to generate the final cDNA libraries. Prior to sequencing, library quality was assessed by Qubit 2.0 quantification, insert-size evaluation using an Agilent 2100 Bioanalyzer, and qPCR based quantification (effective library concentration > 2 nM) [15]. Qualified libraries were pooled and sequenced on an Illumina HiSeq platform (Illumina, Inc., San Diego, CA, USA) using sequencing-by-synthesis (SBS). Raw image data were base-called and converted to FASTQ format using CASAVA (Illumina, Inc., San Diego, CA, USA) to generate raw reads. These raw reads were stringently filtered to obtain clean reads by (i) removing adapter-containing reads; (ii) discarding read pairs in which either read contained >10% ‘N’ bases; and (iii) discarding read pairs in which low-quality bases (Q ≤ 20) accounted for >50% of the read length [15]. Sequencing was performed by Maiwei Metabolic Biotechnology Co., Ltd.

2.6. Transcriptome Data Analysis

Clean reads were de novo assembled using Trinity (v2.6.6) [16], and the resulting transcripts were hierarchically clustered with Corset (v1.07) [16]. The longest sequence in each cluster was retained to generate the unigene set for downstream analyses. Clean reads from each sample were subsequently mapped to the unigene reference, and transcript abundance was estimated using RSEM with Bowtie2 [16,17,18]. Expression levels were normalized to FPKM (Fragments Per Kilobase of transcript per Million fragments mapped). For functional annotation, unigene sequences were searched against the KEGG, NR, Swiss-Prot, GO, COG/KOG, and TrEMBL databases using BLAST (2.16.0), and Pfam annotation was performed using HMMER; annotation statistics were summarized accordingly. Differential expression analysis was conducted using DESeq2 (v1.22.2) based on raw read counts generated by featureCounts. Genes with |log2FC| ≥ 1 and FDR < 0.05 were considered differentially expressed [16,17,18].

2.7. Integrated Metabolome and Transcriptome Analysis

Differentially accumulated metabolites (DAMs) and differentially expressed genes (DEGs) were independently mapped to the KEGG database. Pathways that were significantly enriched in both datasets were extracted to delineate coumarin-related functional modules. Candidate genes were subsequently screened based on their pathway annotations and differential expression patterns between ‘green Sichuan pepper’ and ‘red Sichuan pepper’. To quantitatively assess gene–metabolite relationships, the relative abundances of representative coumarins (scopoletin and scopolin) were paired with unigene expression levels (FPKM) across biological replicates, followed by Pearson correlation analysis. Statistically significant associations (p < 0.05) were retained and assembled into a correlation network to identify genes showing significant associations with scopoletin and/or scopolin.

2.8. Validation of Gene Expression by qRT-PCR

Quantitative reverse transcription polymerase chain reaction (qRT-PCR) was performed using the SYBR^®^ Premix Ex Taq™ II Reagent Kit (Takara Bio Inc., Kusatsu, Shiga, Japan). The internal reference gene was the Zanthoxylum polyubiquitinase gene (UBQ), as reported by Ruijie et al. [19]. Primers were designed using Primer 5.0 (Table S1) with the following conditions: annealing temperature 58–65 °C, primer length 18–24 bp, and amplicom length of 80–150 bp. Relative gene expression was quantified using the 2^−^ΔΔCt method, with reverse-transcribed cDNA as the template [20]. The total reaction volume was 20 μL (Table S3). Cycling conditions were as follows: 95 °C for 2 min; 95 °C for 15 s, 60 °C for 30 s, and 60 °C for 30 s. Melting curve analysis was performed at 95 °C for 15 s, 60 °C for 1 min, 95 °C for 15 s, and 60 °C for 15 s. The cycling program was repeated for 36 cycles, and the reactions were held at 4 °C.

2.9. Candidate Gene Cloning and Overexpression in Tobacco

The three candidate genes, CCoAOMT (Cluster-21146.94405), COMT (Cluster-21146.25660), and F6H (Cluster-21146.64498), were significantly associated with scopoletin and scopolin levels. The complete primer sequences for these three genes are provided in Table S2. The CaMV 35S promoter was used to drive expression of the cloned CDS in a plant overexpression vector. After sequence verification, the recombinant vector was introduced into Agrobacterium tumefaciens strain GV3101. The Agrobacterial suspension carrying the recombinant vector was infiltrated into Nicotiana benthamiana leaves (See Supplementary Materials for details). After infiltration, the plants were maintained under standard growth conditions for 48–72 h, after which the infiltrated leaves were harvested for RNA extraction and qRT-PCR to confirm transgene expression. Metabolite were then extraction as described previously [20,21]. Changes in coumarin content in tobacco leaves, including scopoletin and its glycoside scopolin, were determined using UPLC-MS/MS.

3. Results

3.1. Metabolomics Data Analysis

3.1.1. Sample Stability Assessment and Quantitative Analysis

In both positive and negative ion modes, the total ion chromatograms (TIC) of quality control (QC) samples exhibited a high degree of overlap (Figure S1), indicating stable instrument performance. In negative ion mode, the response peaks of each QC sample showed consistent retention times and peak shapes, with sharp, non-tailing peaks. This observation indicates that neither ionization efficiency nor chromatographic separation performance drifted during the run. In positive ion mode, early-eluting peaks exhibited considerable overlap; moreover, the locations and intensities of the most abundant peaks remained consistent across QC samples. Later-eluting signals were stable, showing no significant signal attenuation or retention-time shift. These results demonstrate good repeatability and stability of the ion source and chromatographic system, supporting their suitability for subsequent metabolite measurements. Using UPLC-ESI-MS/MS, six pericarp samples of green Sichuan pepper and red Sichuan pepper were analyzed, and 583 metabolites were identified and grouped into seven classes (Figure 1B). The classes comprised 246 flavonoids, 150 phenolic acids, 96 alkaloids, 46 lignans and coumarins, 19 tannins, 18 terpenoids, and 8 other compounds. A total of 24 coumarin derivatives were detected, including scopoletin, scopoletin-7-O-glucoside, 6-methylcoumarin, and suberosin (Figure 1C).

3.1.2. PCA and PLS-DA Analysis

The OPLS-DA model showed clear separation structure between the two groups in the predictive component(p1), with R^2^X = 0.855, R^2^Y = 1.000, and Q^2^ = 0.997 (Figure 2A). The model captured differences in metabolite accumulation in Sichuan pepper pericarp and exhibited strong predictive performance. Principal Component Analysis (PCA) revealed intrinsic compositional differences between the two pericarp types across major metabolic classes (Figure 2B). PC1 and PC2 together exceeded >84% of the total variance, indicating that these components captured the major sources of metabolic variation. Green Sichuan pepper and red Sichuan pepper samples were distinctly separated along PC1 axis, while samples within each group clustering tightly, suggesting high within-group consistency of metabolite profiles. The hierarchical clustering heatmap showed marked differences in metabolite distributions between the two pericarp types, indicating systematic shifts in accumulation across multiple metabolite categories. Green Sichuan pepper and red Sichuan pepper formed two distinct clusters in the standardized expression of most metabolites, with particularly pronounced differences observed in flavonoids, phenolic acids, and lignan-coumarin compounds (Figure 2C).

3.1.3. Identification and Analysis of Differential Metabolites

Significant differences were observed in secondary metabolites of the pericarps of two Zanthoxylum species. Differential metabolites were mainly flavonoids (103 upregulated and 58 downregulated metabolites), phenolic acids (39 upregulated and 51 downregulated metabolites), alkaloids (35 upregulated and 21 downregulated metabolites), and lignans and coumarins (14 upregulated and 18 downregulated metabolites) (Figure 3A).

To further investigate differences in coumarin constituents between the pericarps of the two Zanthoxylum samples, we conducted a volcano plot analysis of coumarin-derived secondary metabolites (Figure 3B). The results showed that among the 24 coumarin-related metabolites, 16 differed significant between DHPHJ and DTHJ, comprising eight upregulated and eight downregulated compounds, whereas the remaining eight metabolites showed no significant changes. Comparison with green Sichuan pepper, red Sichuan pepper showed a bidirectional shift in its coumarin profile: several methoxylated and furanocoumarin-related metabolites accumulated significantly, whereas certain coumarins and their conjugated derivatives decreased markedly. These findings suggest a clear redistribution of metabolic flux within the coumarin network between the two sample groups. Hierarchical clustering was performed for the 24 coumarin compounds across both samples (Figure 1C). The clustering results indicated that the concentrations of nine compounds, including 4-Hydroxycoumarin, 7-Methoxycoumarin, Ayapin, and 5,7-Dimethoxycoumarin, were significantly higher in DHPHJ. In contrast, the remaining 15 compounds, such as Aurapten, Nodakenin, Esculin, and Scopoletin-7-O-glucoside (Scopolin), formed a major cluster and showed higher concentrations in green Sichuan pepper. This pattern indicates that most coumarin compounds were present at relatively higher concentrations in green Sichuan pepper.

3.1.4. KEGG Enrichment Analysis

In the comparison between green Sichuan pepper and red Sichuan pepper, 363 differential metabolites were identified and mapped to the KEGG database, and 32 metabolic pathways were enriched. The top 20 enriched KEGG pathways are shown in Figure 3C, indicating that the most significantly overrepresented pathways are closely related to secondary metabolism. Notably, biosynthesis of secondary metabolites and flavonoid biosynthesis exhibited the strongest enrichment signals. In contrast, enrichment in phenylpropanoid biosynthesis and related upstream amino acid pathways, such as phenylalanine metabolism and biosynthesis, was comparatively weaker. Importantly, the 363 differential metabolites were predominantly enriched in the phenylpropanoid biosynthesis pathway (ko00940), the biosynthesis of secondary metabolites (ko01110), and flavonoid biosynthesis (ko00941). Specifically, the coumarin compounds scopoletin and scopolin showed significant enrichment within phenylpropanoid biosynthesis.

3.2. Transcriptomic Analysis

3.2.1. Quality Assessment of Sequencing

DataHigh-throughput sequencing was performed on six pericarp samples of Z. bungeanum using the Illumina HiSeq platform, with three biological replicates for each of green Sichuan pepper and red Sichuan pepper. Separate cDNA libraries were constructed for each sample. The number of clean reads ranged from 4.26 × 10^7^ to 5.28 × 10^7^, per sample, and sample yielding over 6 Gb of clean data. The average GC content exceeded 43%, with Q20 values surpassing 97% and Q30 values exceeding 93%, while the sequencing error rate was maintained at 0.03% (Table 1). These quality metrics indicated that the data were suitable for downstream transcriptome analysis. Using Trinity, a total of 191,959 unigenes were assembled and annotated against seven major databases: KEGG, NR, SwissProt, Trembl, GO, KOG, and Pfam. The annotation counts were as follows: 96,043 (50.03%); 136,059 (70.88%); 93,201 (48.55%); 134,046 (69.83%); 82,909 (43.19%); 108,677 (56.61%); and 86,321 (44.97%).

Following PCA (Figure 4A), the green Sichuan pepper and red Sichuan pepper groups separated clearly along PC1, with PC1 and PC2 explaining 85.54% and 4.35% of the total variation, respectively. Both groups showed tight within-group clustering, indicating high reproducibility among biological replicates. The sample correlation heatmap (Figure 4B) showed within-group correlation coefficients of approximately 0.99 for red Sichuan pepper and 0.98 for green Sichuan pepper, whereas between -group correlations were significantly lower (0.63–0.67), suggesting substantial transcriptional differences between the two materials. Furthermore, box plots of FPKM expression distributions (Figure 4C) reinforce the consistency of expression across samples, indicating an adequate sequencing depth that is suitable for differential expression analysis.

3.2.2. Expression Characteristics of DEGs

A total of 68,619 differentially expressed genes (DEGs) were identified between green Sichuan pepper and red Sichuan pepper, with 37,721 upregulated and 30,898 downregulated (Figure 5A). By mapping these DEGs to the KEGG database, we characterized metabolic pathways represented in Sichuan pepper pericarp and identified DEGs associated with coumarin-related biosynthesis. The DEGs were significantly enriched in pathways related to metabolism, biosynthesis of secondary metabolites, and phenylpropanoid biosynthesis (Figure 5B).

Thirteen genes were selected for qRT-PCR validation. The RNA-seq heatmap results indicated that key genes involved in the phenylpropanoid and coumarin-related pathways exhibited distinct group-specific expression patterns (Figure 5C). Multiple transcripts showed higher expression in green Sichuan pepper, whereas another subset was preferentially expressed in red Sichuan pepper, highlighting substantial transcriptional divergence between the two pericarp materials. These findings further suggest that both the upstream steps of the phenylpropanoid pathway and the downstream modification processes of coumarin biosynthesis may be differentially regulated. To further validate the transcriptomic data, qRT-PCR analysis was conducted using SUBA as the reference gene (Figure 5D). The qRT-PCR results corroborated the expression trends observed in the RNA-seq data, and the relative expression changes between red Sichuan pepper and green Sichuan pepper showed high concordance with the transcriptomic results for each gene. This concordance between qRT-PCR and RNA-seq supports the reliability of the expression patterns of these candidate genes.

3.3. Reconstruction of the Coumarin Biosynthetic Pathway in Zanthoxylum Pericarp and Screening of Candidate Genes

By integrating KEGG enrichment results from metabolomics and transcriptomics, we constructed a biosynthetic pathway framework for the simple coumarins scopoletin and scopolin in the pericarp of Z. bungeanum (Figure 6A). The upstream pathway initiates with the PAL-catalyzed deamination of phenylalanine, resulting in the formation of cinnamic acid, which is subsequently converted to p-coumaroyl-CoA through the action of Cinnamate-4-hydroxylase (C4H) and 4-coumarate-CoA ligase (4CL). Enzymatic processes involving hydroxycinnamate O-methyltransferase (HCT), Cinnamate-3-hydroxylase (C3′H), Caffeoyl shikimate esterase (CSE), Caffeoyl-CoA O-methyltransferase (CCoAOMT), and caffeic acid O-methyltransferase (COMT) yield various hydroxycinnamic acid derivatives. The key ortho-hydroxylation step is mediated by cytochrome P450 enzymes (CYP73A) and 2-oxoglutarate-dependent dioxygenases (F6H), resulting in intermediates that readily undergo lactonization. Subsequently, enzymes such as TOGT1 and BGLU facilitate the glycosylation–deglycosylation cycle, culminating in the formation of scopoletin and scopolin.

Based on pathway localization and differential expression analysis, a total of 102 coding genes potentially involved in the coumarin biosynthesis of Zanthoxylum were identified. This includes 9 PAL, 7 4CL, 4 C3′H, 9 CCoAOMT, 17 COMT, 8 CYP73A, 6 F6H, 19 HCT, 11 TOGT1, 13 BGLU, and 3 CSE genes (Figure 6B). Most of these genes showed significantly higher expression in DTHJ than in DHPHJ, consistent with the higher overall coumarin content in DTHJ. Further correlation analysis between the contents of scopoletin and scopolin and the expression levels of these genes identified 56 candidate genes that were highly correlated with the contents of these two coumarins, including PAL, 4CL, CCoAOMT, and COMT. The correlation network indicated that some midstream enzymes (CCoAOMT, COMT, and F6H) exhibited high correlation coefficients with both coumarins, thus displaying ‘hub’ node characteristics within the network (Figure 6C).

3.4. Functional Validation of CCoAOMT, COMT, and F6H

Based on the correlation network and expression patterns, three candidate genes—CCoAOMT (Cluster-21146.94405), COMT (Cluster-21146.25660), and F6H (Cluster-21146.64498)—were selected for gene cloning and heterologous expression analyses. The constructed overexpression vectors were introduced into tobacco leaves via Agrobacterium-mediated transient infiltration (Figure S2). PCR verification confirmed that constructs carrying the three target genes were detectable in tobacco leaves after transient infiltration (Figure S3), and transient overexpression was achieved (Figure S4). UPLC-based quantitative analysis demonstrated that, compared with the control, the contents of scopoletin and scopolin significantly increased in tobacco leaves overexpressing each of the three genes, with more pronounced increases observed in the COMT- and F6H-overexpression lines (Figure 7). Coupled with their expression profiles in Zanthoxylum pericarp and the metabolite-gene correlation analysis, these results indicate that the three genes play crucial catalytic roles in the biosynthesis of coumarins in Zanthoxylum, particularly in the formation of scopoletin and its glycoside.

4. Discussion

This study analyzed coumarin metabolites and their synthesis-related genes in the mature pericarp of two Sichuan pepper cultivars. Significant variations were identified in the overall secondary metabolite composition between green Sichuan pepper and red Sichuan pepper, with green Sichuan pepper pericarp demonstrating higher overall abundance of simple coumarins. This conclusion is consistent with previous reports of differences in volatile oils, flavonoids, and amides among Sichuan pepper cultivars, and further supports the observation that green Sichuan pepper contains higher coumarin levels than red Sichuan pepper [22,23,24,25].

The coumarin biosynthesis framework proposed in this study is fundamentally consistent with the phenylpropanoid branch pathways previously discovered in plants of the Apiaceae, Rutaceae, and Brassicaceae families. This process originates with PAL-C4H-4CL, catalysed by enzymes including HCT, C3′H, and CCoAOMT/COMT to create various hydroxycinnamic acid derivatives. Subsequently, crucial ortho-hydroxylation reactions are catalysed by 2-oxoglutarate-dependent dioxygenases such as F6′H, leading to intermediates that undergo lactonisation and glycosylation [22,23,24,25]. Multiple gene family members including PAL, 4CL, HCT, C3′H, CCoAOMT, COMT, F6H, TOGT1 and BGLU were similarly discovered in Sichuan pepper, demonstrating that coumarin production in this plant follows the normal phenylpropanoid metabolic pathway. The diversity and abundance of coumarins in Sichuan pepper pericarp differ markedly from those reported in species such as Arabidopsis thaliana and Angelica sinensis, in which roots are the principal medicinal tissues [25,26,27]. For example, in A. sinensis, simple coumarins are primarily concentrated in roots, whereas Sichuan pepper pericarp exhibits a distinct coumarin profile consistent with its use as a seasoning [28]. Variations in the division of labour across various plants or tissues in defensive functions, signal transduction, and interactions with microbes could be the cause of this specificity discrepancy. This pattern may also be influenced by synergistic defense networks involving other secondary metabolites, including volatile oils and amides [29].

Although green Sichuan pepper DTHJ and red Sichuan pepper differ in coumarin content, the expression patterns of major enzyme-encoding genes within the CCoAOMT, COMT, and F6H families also show substantial divergence. Research reveals that these enzymes hold critical locations within the phenylpropanoid metabolic pathway, where even modest alterations in their expression can drastically alter the distribution of carbon flux among diverse products such as flavonoids, coumarins, and lignin [30]. In this study, most structural genes associated with coumarin synthesis exhibited higher expression in green Sichuan pepper, while several genes involved in lignin biosynthesis showed relatively higher expression in red Sichuan pepper. This pattern suggests that carbon flux allocation may differ between the two materials, potentially favouring coumarin and flavonoid biosynthesis in green Sichuan pepper; however, this interpretation remains hypothetical and requires further experimental validation, such as lignin quantification, total phenolic measurements, or metabolic flux analysis. CCoAOMT, COMT, and F6H exhibited strong positive correlations with scopoletin and scopoletin glycoside within the network and demonstrated the capacity to enhance coumarin content in heterologous systems. This observation is consistent with reports in other species identifying F6′H and related methyltransferases as key enzymes in coumarin biosynthesis [24]. However, the presence of numerous members within the TOGT1 and BGLU gene families in Sichuan pepper, coupled with the significant correlation of multiple members with coumarin content, suggests that glycosylation-deglycosylation cycles may play a more pivotal role in regulating the dynamic equilibrium of coumarin in pericarp. Although such cycles have been reported less frequently in species in which roots are the primary site of coumarin accumulation, further functional evidence is required to substantiate this regulatory role.

5. Conclusions

Metabolomic analysis of mature pericarps from green Sichuan pepper and red Sichuan pepper identified a total of 583 metabolites, including 24 coumarins that differed significant between the two Zanthoxylum materials. The pericarp of DTHJ showed a higher overall abundance of simple coumarins. Integrating metabolomic and transcriptomic datasets identified 102 structural genes potentially involved in coumarin metabolism; among these, 56 candidates showed significant correlations with scopoletin glycoside accumulation. Transient overexpression assays conducted in Nicotiana benthamiana revealed that CCoAOMT (Cluster-21146.94405), COMT (Cluster-21146.25660), and F6H (Cluster-21146.64498) significantly enhanced the accumulation of scopoletin and scopolin, thereby supporting their roles in coumarin biosynthesis. These functional data were obtained from heterologous transient expression, and genetic validation in Zanthoxylum remains to be established. Future studies employing gene knockdown or knockout techniques (e.g., RNAi or CRISPR/Cas9), in vitro enzymatic assays, and absolute metabolite quantification using internal standards will be essential to further validate gene functions and refine the regulatory framework of coumarin biosynthesis.