Method for revealing molecular mechanism of zanthoxylum bungeanum leaf terpenoid biosynthesis based on metabonomics and whole transcriptome analysis

Abstract

The invention belongs to the technical field of plant molecular biology and plant genetic engineering, and particularly relates to a method for revealing a biosynthesis molecular mechanism of zanthoxylum bungeanum leaf terpenoids based on metabonomics and whole transcriptome analysis. Through multi-omics analysis on the Chinese prickly ash leaves, the terpenoids are determined to be main contributing substances of the fragrance of the Chinese prickly ash leaves. Key candidate regulatory factors are screened and identified in combination with correlation analysis of gene expression level and terpenoid content, and the factors may participate in synthesis regulation of zanthoxylum bungeanum maxim terpenoid. On the basis of the key candidate regulatory factors, a regulatory network involving 9 DEmRNAs, 8 DEmiRNAs, 8 DElncRNAs and 14 DEcircRNAs is constructed, lncRNA and circRNA in the network are competitively combined with miRNA, expression of the key candidate regulatory factors is regulated in a targeted mode, and then the key candidate regulatory factors cooperatively participate in synthesis regulation of the zanthoxylum bungeanum maxim terpenoid.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of plant molecular biology and plant genetic engineering, and particularly relates to a method for revealing the molecular mechanism of terpenoid biosynthesis in Zanthoxylum bungeanum leaves based on metabolomics and whole transcriptome analysis. BACKGROUND

[0002] Zanthoxylum bungeanum (Z. bungeanum) as an important economic tree species of Zanthoxylum, is a typical representative of Chinese characteristic spice crops, and has important medicinal value and edible value. Studies have shown that volatile terpenoids rich in Z. bungeanum leaves and pericarp are the main contributing substances of its unique aroma characteristics. As an important secondary metabolite in plants, terpenoids not only play a key role in plant growth and development and environmental adaptation, but also have significant pharmacological activities such as antibacterial, anti-inflammatory and antioxidant activities. Zanthoxylum armatum The biosynthesis process of terpenoids is a highly complex regulatory network involving multiple levels of regulatory mechanisms. From the perspective of biosynthetic pathways, terpenoids in Z. bungeanum are mainly synthesized through the cytoplasmic mevalonate (MVA) and plastidic methylerythritol phosphate (MEP) pathways. These two pathways provide basic precursor substances for terpenoid synthesis, and then the backbone structure is formed through the catalysis of terpene synthase (TPS), and finally the diverse terpenoids are formed through the modification of cytochrome P450 monooxygenase (CYP450s) and glycosyltransferase (UGTs) enzymes.

[0003] However, current research on Z. bungeanum mainly focuses on the extraction and identification of volatile components and the analysis of pharmacological activities in two directions, and the key molecular regulation mechanism of aroma formation, especially the terpenoid synthesis regulation network based on the synergistic action of metabolomics and whole transcriptome (mRNA, lncRNA, circRNA and miRNA), still lacks systematic and in-depth research.

[0004] Therefore, we propose a method for revealing the molecular mechanism of terpenoid biosynthesis in Z. bungeanum leaves based on metabolomics and whole transcriptome analysis, hoping to solve the deficiencies in the prior art.

[0005] SUMMARY The purpose of the present application is to provide a method for revealing the molecular mechanism of terpenoid biosynthesis in Z. bungeanum leaves based on metabolomics and whole transcriptome analysis, which solves the problems in the prior art.

[0006] The purpose of the present application is to provide a method for revealing the molecular mechanism of terpenoid biosynthesis in Z. bungeanum leaves based on metabolomics and whole transcriptome analysis, which solves the problems in the prior art.

[0007] The present application is realized by the following technical scheme: A method for revealing the molecular mechanism of terpenoid biosynthesis in Z. bungeanum leaves based on metabolomics and whole transcriptome analysis, comprising the following steps: S1, volatile metabolome analysis: The volatile metabolites were detected by GC-MS, and the detected metabolites were subjected to hierarchical cluster analysis (HCA), principal component analysis (PCA) and K-means cluster analysis to screen differential metabolites (DTM), and the key flavor substances with rOAV≥1 were identified based on the relative odor activity value; S2, RNA extraction and sequencing: The total RNA of the three groups of leaves of Zanthoxylum planispinum, i.e., young leaves, mature leaves and old leaves, was extracted by using RNA Pure Plant Kit, and the RNA quality was evaluated by NanoDrop 2000 and Agilent Bioanalyzer 2100 system; chain-specific RNA-seq library and small RNA sequencing library were constructed respectively, high-throughput sequencing was performed, and the original sequencing data were obtained and processed through BMKCloud platform; S3, analysis of mRNA, lncRNA, miRNA and circRNA: The sequencing data of step S2 were subjected to analysis of mRNA, lncRNA, miRNA and circRNA; S4, construction of terpenoid synthesis regulation network: Based on the key flavor substances of terpenoids identified in step S1 and the differential expressed RNAs (DEG) screened in step S3, the key candidate genes involved in terpenoid synthesis were screened by correlation analysis; based on the ceRNA hypothesis, the interaction relationship between the key candidate genes and DEmiRNAs, DElncRNAs and DEcircRNAs was analyzed, and a complete ceRNA regulation network was constructed; S5, verification: The expression levels of DEmRNAs, DElncRNAs and DEmiRNAs screened in step S3 were verified by qRT-PCR, and the reliability of the sequencing data was confirmed.

[0008] Further, the young leaves are Zanthoxylum planispinum leaves YL 10 days after germination, the mature leaves are Zanthoxylum planispinum leaves ML 60 days after germination, and the old leaves are Zanthoxylum planispinum leaves SL 120 days after germination.

[0009] Further, in step S1, the GC-MS is Agilent Model 8890-7000D, and each sample is subjected to three biological repeats; in the differential metabolite screening, the DTM of ML_vs_YL, SL_vs_ML and SL_vs_YL groups are 721, 136 and 714 respectively, and are enriched to secondary metabolite synthesis, monoterpene biosynthesis, sesquiterpene and triterpene biosynthesis pathways through metabolic pathway annotation.

[0010] Further, in step S2, the ribosomal RNA is removed using Ribo-Zero rRNA Removal Kit when constructing the strand-specific RNA-seq library, and NEBNext® Ultra™ Directional RNA Library Prep Kit is used to construct the strand-specific library; and when constructing the miRNA library, TruSeq Small RNA Sample Prep Kit is used.

[0011] Further, in step S3, the mRNA and lncRNA analysis is as follows: the sequencing sequences are aligned to the reference genome using Hisat2 (https: / / doi.org / 10.6084 / m9.figshare.14400884.vl), the transcripts with a length > 200 nt and containing ≥ 2 exons are screened as lncRNA candidates, and the lncRNAs are identified by four calculation methods of CPC2, CNCI, Pfam and CPAT; the FPKM values of the mRNA and lncRNA are calculated using StringTie, and the DESeq R software package is used to screen the differentially expressed mRNA (DEGs) and lncRNA (DElncRNAs) with the standard of FC ≥ 2 and Pvalue < 0.01; the cis-target genes of the lncRNA are predicted based on the positional relationship (the adjacent genes within the range of 100 kb upstream and downstream), the trans-target genes are predicted based on the expression correlation, and the above target genes are annotated for functional information by Nr, Pfam, KOG / COG, Swiss-Prot, KO and GO six public databases.

[0012] Further, in step S3, the miRNA analysis is as follows: the known and new miRNAs are identified using miRDeep2 software, the target genes are predicted using TargetFinder, the expression amount is normalized using TPM algorithm, and the differentially expressed miRNAs (DEmiRNAs) are screened according to |log2(FC)| ≥ 1.00 and FDR ≤ 0.01.

[0013] Further, in step S3, the circRNA analysis is as follows: the circRNAs are identified using find_circ software, the expression level is represented by junction reads count and normalized by SRPBM, the differentially expressed circRNAs (DEcircRNAs) are screened according to FC ≥ 2 and P < 0.05, and the miRNA-circRNA binding pairs are identified using Target Finder software.

[0014] Further, the regulatory network in step S4 comprises 9 DE mRNAs, 8 DE miRNAs, 8 DE lncRNAs and 14 DE circRNAs, wherein the lncRNA and the circRNA relieve the inhibition of the miRNA on the key gene of terpenoid synthesis by competitively binding the miRNA.

[0015] Further, the qRT-PCR verification in step S5 adopts AceQ Universal SYBR qPCR MasterMix and HiScript III RT SuperMix to detect the lncRNA and the mRNA, adopts miRNA Unimodal SYBR qPCR Master Mix to detect the miRNA, takes U6 as a reference gene, and calculates the relative expression level by 2 -ΔΔCt The relative expression level is calculated by the method, and the verification result is consistent with the trend of the sequencing data.

[0016] Compared with the prior art, the present application has the following advantages: 1. The present application jointly analyzes volatile metabolome and whole transcriptome (including mRNA, lncRNA, circRNA and miRNA) data, and analyzes the mechanism of formation of main aroma substances of Zanthoxylum bungeanum Maxim. through three research dimensions: (1) identifying key aroma active substances based on the dynamic characteristics of volatile metabolome; (2) mining key structural genes and regulatory transcription factors of terpenoid synthesis; (3) constructing a ceRNA regulatory network to reveal the regulatory effect of non-coding RNA on terpenoid synthesis. The research results will provide new insights into the molecular mechanism of terpenoid biosynthesis of Zanthoxylum bungeanum Maxim., and also provide a theoretical basis for quality improvement and efficient utilization of functional components of Zanthoxylum bungeanum Maxim..

[0017] 2. The present application carries out whole transcriptome sequencing on Zanthoxylum bungeanum Maxim. leaves in three development periods of young leaves, mature leaves and old leaves, aiming to mine information related to the synthesis regulation of terpenoids. The differentially expressed mRNA, lncRNA, circRNA and miRNA are identified, and candidate key regulatory factors such as CYP450, bHLH, bZIP, UTG, GDS, LMS, HMGS and DXS are screened, which are related to the change trend of the content of terpenoids. Further, a ceRNA-miRNA-target regulatory network related to terpenoid synthesis of Zanthoxylum bungeanum Maxim. leaves is constructed, and it is clear that the lncRNA and the circRNA can act as ceRNA, regulate the expression of UTG11, DXS3 and other candidate genes of terpenoid synthesis by competitively binding the miRNA, and participate in the terpenoid biosynthesis pathway. The present application reveals the molecular regulation mechanism of terpenoid synthesis of Zanthoxylum bungeanum Maxim. leaves, identifies potential mRNAs and ncRNAs involved in terpenoid synthesis, and provides theoretical support and research reference for in-depth exploration of terpenoid metabolism of Zanthoxylum bungeanum Maxim. and regulation of secondary metabolism of spice plants.

[0018] 3. This invention, through multi-omics analysis of Sichuan pepper leaves, clarified that terpenoids are the main contributors to the aroma of Sichuan pepper leaves. Among the 1468 volatile metabolites detected, 320 were terpenoids, forming the core components of aroma substances. Combining gene expression level and terpenoid content correlation analysis, key candidate regulatory factors, including transcription factors such as WRKY and bHLH, CYP450 and UGTs family genes, and GDS, LMS, HMGS, and DXS, were screened and identified. These factors may participate in the synthesis and regulation of Sichuan pepper terpenoids. Based on these key candidate regulatory factors, a regulatory network involving 9 DEmRNAs, 8 DEmiRNAs, 8 DElncRNAs, and 14 DEcircRNAs was constructed. In this network, lncRNAs and circRNAs competitively bind to miRNAs, targeting and regulating the expression of key candidate regulatory factors, thereby synergistically participating in the synthesis and regulation of Sichuan pepper terpenoids. In summary, this study clarified the core role of terpenoids in the aroma formation of Sichuan pepper leaves and elucidated the molecular network regulating their synthesis, providing a theoretical basis and genetic resources for a deeper understanding of the terpenoid metabolism mechanism in Sichuan pepper and the quality improvement of spice plants. Attached Figure Description

[0019] Figure 1 Volatile metabolomics analysis of leaves at three different developmental stages of Sichuan pepper; Figure 2 Flavoromics analysis for DTM; Figure 3 mRNA analysis of Sichuan pepper leaves; Figure 4 Pearson correlation analysis network diagram for gene expression levels and relative metabolite expression; Figure 5 Analysis of lncRNA in Sichuan pepper leaves; Figure 6 Analysis of miRNAs in Sichuan pepper leaves; Figure 7 To identify and analyze circular RNAs in leaves and to construct a ceRNA network that regulates terpene biosynthesis; Figure 8 RT-PCR validation of miRNAs, lncRNAs, and mRNAs in Sichuan pepper leaves. Detailed Implementation

[0020] To further explain the present invention, the following specific embodiments are described.

[0021] plant materials Based on the color of the peel, Sichuan peppercorns are divided into red and green varieties. The material used in this invention is green Sichuan peppercorns (…). Zanthoxylum armatum var. novemfoliusThe samples were planted in a standardized Sichuan pepper plantation in Rongchang District, Chongqing (105.93°E, 29.36°N). To ensure biological reproducibility and data reliability, 36 genotype-consistent Sichuan pepper trees were randomly selected and divided into 6 groups, with 6 trees in each group serving as independent biological repeat sequences. Pest-free leaves were collected at different developmental stages: Stage 1 was 10 days after budding (denoted as YL), Stage 2 was 60 days after budding (denoted as ML), and Stage 3 was 120 days after budding (denoted as SL). After collection, the samples were immediately frozen in liquid nitrogen and stored in an ultra-low temperature freezer at -80°C in preparation for subsequent experimental analysis.

[0022] Example 1 1. Volatile metabolomics analysis A suitable amount of fresh sample was ground in liquid nitrogen, and approximately 500 mg was weighed into a 20 mL headspace vial containing a saturated NaCl solution (Agilent, Palo Alto, CA, USA). The vial was shaken for 5 min at a constant temperature of 60 °C. A 120 µm DVB / CWR / PDMS extraction head was inserted into the sample headspace vial, and headspace extraction was performed for 15 min. The sample was then desorbed at 250 °C for 5 min, and then separated and determined by GC-MS (Model 8890, 7000D; Agilent). Each sample was subjected to three biological replicates.

[0023] GC-MS conditions: Column, DB-5MS (30 m × 0.25 mm × 0.25 μm, Agilent J&W Scientific, Folsom, CA, USA); Carrier gas, high-purity helium; Flow rate, 1.2 mL / min; Injector temperature, 250°C, splitless, solvent delay 3.5 min. Temperature program: 40°C, hold for 3.5 min, ramp to 100°C at 10°C / min, ramp to 180°C at 7°C / min, and finally ramp to 280°C at 25°C / min, hold for 5 min. Ionization mode, electron collision (EI); Ionization energy, 70 eV; Detector, quadrupole mass analyzer; Temperature settings: quadrupole, 150°C; Ion source, 230°C; transfer line, 280°C; Selected ion monitoring (SIM) mode for identification and quantification of target analytes. Differential metabolite analysis: Differential metabolites were screened based on variable importance (VIP) in the projection (VIP > 1) and absolute log2FC (|log2FC| ≥ 1.0). The identified metabolites were then annotated using the KEGG compound database (http: / / www.kegg.jp / kegg / compound / ).

[0024] Previous research by the inventors of this invention has shown that the content of terpenoid compounds in Zanthoxylum armatum leaves varies significantly at different developmental stages (Liu X, Tang N, Xu F, Chen Z, Zhang X, Ye J, Liao Y, Zhang W, Kim SU, Wu P, Cao Z. SMRT and Illumina RNA sequencing reveal the complexity of terpenoid biosynthesis in Zanthoxylum armatum. Tree Physiol. 2022 Mar 9;42(3):664-683.). Given the unique numbing and aromatic flavor of Zanthoxylum armatum leaves, this study systematically analyzed the composition and changes of volatile metabolites in leaves at three developmental stages—young leaves, mature leaves, and old leaves—using volatile metabolomics sequencing. A total of 1468 volatile metabolites were identified, covering 15 categories of substances, with terpenoids being the most numerous (320 types), followed by esters (234 types). Hierarchical cluster analysis (HCA) and principal component analysis (PCA) showed differences in metabolites among the three groups (…). Figure 1 A and B). Further screening yielded 854 differentially expressed metabolites (DTMs), of which 49 metabolites showed significant differences across all three comparison groups. Figure 1 C). Specifically, 721 (461 down-regulations; 260 up-regulations) DTMs were identified between ML_vs_YL, SL_vs_ML, and SL_vs_YL, respectively; 136 (115 down-regulations; 21 up-regulations) DTMs were identified between SL_vs_YL and SL_vs_YL, respectively; SL_vs_YL identified 721 (461 down-regulations; 260 up-regulations) DTMs were identified between SL_vs_YL and SL_vs_YL, respectively; SL_vs_YL identified 136 (115 down-regulations; 21 up-regulations) DTMs were identified between SL_vs_YL and SL_vs_YL, respectively; SL_vs_YL identified 714 (538 down-regulations; 17 Figure 1 D). Among them, the number of differentially expressed metabolites in the old leaf and mature leaf comparison group was significantly lower than that in the other two groups, indicating that there are large differences in aroma substances between young leaves and mature and old leaves of Sichuan pepper, while the aroma substance composition between mature leaves and old leaves is relatively stable. K-means cluster analysis divided the 854 differentially expressed metabolites into 7 categories, and the results showed that most differentially expressed metabolites were present in higher amounts in young leaves. Figure 1E), which may be related to the active synthesis of volatiles during the rapid growth and development of young leaves. Metabolic pathway annotation revealed that 66 out of 854 DTMs were annotated to different metabolic pathways, mainly enriched in secondary metabolite biosynthesis (ko01110), metabolic pathways (ko01100), monoterpene biosynthesis (ko00902), sesquiterpenoid and triterpenoid biosynthesis (ko00909), and tropane, pipeidine and pyridine alkaloid biosynthesis (ko00960), further confirming the important role of terpenoids in the aroma formation process of Sichuan pepper leaves.

[0025] Figure 1 (A) Hierarchical cluster analysis of all volatile metabolites, with each metabolite represented by a column; (B) Principal component analysis (PCA) of volatile metabolites; (C) Venn diagram of shared DTMs among the ML_vs_YL, SL_vs_ML, and SL_vs_YL groups; (D) Number of upregulated and downregulated DTMs in the ML_vs_YL, SL_vs_ML, and SL_vs_YL comparison groups; (E) Cluster analysis of K-Means DTM expression patterns.

[0026] Relative odor activity (rOAV) is an important method for identifying key flavor components in food and elucidating their contribution to overall aroma, based on the sensory threshold of compounds. Generally, when rOAV ≥ 1, it indicates that the compound directly affects the flavor of the sample. Of the 1468 volatile metabolites, 675 were flavor contributors (rOAV > 1), belonging to 14 categories. Terpenes were the most numerous (105, 15.5%), followed by esters (94, 13.9%), heterocyclic compounds (76, 11.3%), alcohols (71, 10.5%), and ketones (62, 9.2%). Notably, terpenes (105) were the most numerous flavor contributors in all three different developmental stages of Sichuan pepper leaves (Table 1), confirming them as the core components of Sichuan pepper leaf aroma. To further explore the formation mechanism of aroma compounds in Sichuan pepper, sensory flavor annotation was performed on differential metabolites. A total of 442 differential metabolites were annotated, including 166 key flavor compounds with rOAV > 1, mainly terpenoids (50 species). Comparative analysis revealed that in the ML_vs_YL and SL_vs_YL groups, sweet, green, and fruity were the top three flavor characteristics. Figure 2 A and Figure 2 C), corresponding to the most abundant differential metabolites; in the SL_vs_ML group, the top three flavor characteristics are fruity, green and citrus ( Figure 2 B). These flavor characteristics are formed by the synergistic effect of multiple compounds, creating a complex flavor network. Figure 2 Terpenoids are the core contributors to the unique aroma of Sichuan pepper.

[0027] Figure 2 (A) Radar chart of sensory flavor characteristics analysis using DTM; (B) Correlation network diagram between sensory flavor characteristics and DTM. Pink circles represent sensory flavor characteristics, and green circles represent different metabolites. The larger the pink circle, the more differentially expressed metabolites are annotated for sensory flavor characteristics; the larger the green circle, the more sensory flavor characteristics are annotated by DTM.

[0028] Table 1. Number of compounds with rOAV > 1 in Sichuan pepper leaves at different developmental stages

[0029] 2. RNA extraction, library construction, and RNA sequencing Total RNA was extracted from Sichuan pepper leaves at different developmental stages using the RNA prep Pure Plant Kit (Tiangen, Beijing, China). RNA quality was assessed using a NanoDrop 2000 (Thermo Fisher Scientific, Wilmington, DE) and an Agilent Bioanalyzer 2100 system (Agilent Technologies, CA, USA), and three biological replicates were prepared for each stage. For mRNAs, lncRNAs, and circRNAs, ribosomal RNA was removed using the Ribo-Zero rRNARemoval Kit (Epicentre, Madison, WI, USA) and analyzed using NEBNext® Ultra. TM Strand-specific RNA-seq libraries were constructed using the Directional RNA Library Prep Kit; for miRNAs, small RNA sequencing libraries were constructed using the TruSeq Small RNA Sample Prep Kit (Illumina). All raw reads were further processed using the bioinformatics analysis platform BMKCloud (www.biocloud.net).

[0030] 3. mRNA analysis, lncRNA identification, and target gene prediction For lncRNAs and mRNAs, Hisat2 was used to align the pruned sequences to a reference genome (https: / / doi.org / 10.6084 / m9.figshare.14400884.v1) to identify known transcripts and predict new ones. For unknown transcripts, those longer than 200 nt with more than two exons were selected as lncRNA candidates, and potential lncRNAs were further screened using four computational methods: CPC2 / CNCI / Pfam / CPAT. Then, antisense and sense lncRNAs were selected using the SenseRNA software. StringTie (version 1.3.1) was used to calculate the FPKMs of lncRNAs and mRNAs in the samples, and the DESeq R package (1.10.1) was used to analyze the differentially expressed lncRNAs (DElncRNAs) and mRNAs (DEmRNAs) between pairs of groups (FC ≥ 2 and P value < 0.01).

[0031] In the regulation of neighboring gene expression by lncRNAs, the prediction of cis-target mRNAs is mainly based on the positional relationship between the lncRNA and the gene. Neighboring genes within a 100kb range upstream and downstream of the lncRNA are typically considered its cis-target genes. Trans-target genes are identified by analyzing the expression correlation between lncRNAs and genes. Subsequently, functional information annotation of the aforementioned target genes was performed based on six public databases, including Nr (NCBI non-redundant protein sequences), Pfam (Proteinfamily), KOG / COG (Clusters of Orthologous Groups of proteins), Swiss-Prot (Amanually annotated and reviewed protein sequence database), KO (KEGG Orthologdatabase), and GO (Gene Ontology).

[0032] 3.1 Analysis of the influence of mRNA on the synthesis of terpenoids To elucidate the molecular mechanism of Sichuan pepper aroma formation, this invention performed eukaryotic RNA-Seq analysis on Sichuan pepper leaves at three different developmental stages. Nine libraries were sequenced using the Illumina HiSeq platform, yielding a total of 61.83 Gb of clean data. The GC content of each transcriptome ranged from 44.14% to 44.54%, and the Q30 value reached 95.43% to 96.08% (Table 2). The sequencing reads were mapped to the Sichuan pepper genome map using TopHat software for alignment. The mapping rate for each library was ≥83.46%, indicating that the data quality met the requirements for subsequent analysis. Gene expression levels were analyzed using the FPKM method. In the ML_vs_YL, SL_vs_ML, and SL_vs_YL groups, 10286 (5577 upregulated, 4709 downregulated), 757 (421 upregulated, 336 downregulated), and 11796 (6550 upregulated, 5246 downregulated) differentially expressed mRNAs (DEGs) were identified, respectively. Figure 3 A). The results showed that gene expression was most active in the young leaf stage, and the gene expression pattern changed dynamically with the leaf development process. The number of differentially expressed genes was the largest between old and young leaves, suggesting that they may be involved in the specific regulation of growth and senescence-related pathways.

[0033] Table 2 Evaluation of Sample Sequencing Data

[0034] Further analysis based on the annotation information revealed that 117 DEGs (Derivative Genetic Organisms) were screened in the terpene biosynthesis pathway (Sesquiterpenoid and triterpenoid biosynthesis: ko00909, Terpenoid backbone biosynthesis: ko00900, Diterpenoid biosynthesis: ko00904, Monoterpenoid biosynthesis: ko00902) (Table 3). Figure 3 B), and most of these differentially expressed genes are located downstream of metabolic pathways ( Figure 3 (C) indicates that these structural genes may be involved in the synthesis of terpenoids in Sichuan pepper leaves. Transcription factors can regulate the production of secondary metabolites by controlling the transcription initiation of genes in metabolic pathways. Current research has found that the main families of transcription factors involved in terpenoid synthesis include AP2 / ERF, bHLH, WRKY, and bZIP. In addition, different types of terpenoids are also modified by CYP450s and UGTs. In this invention, a total of 382 DEGs that may be involved in terpenoid synthesis were screened, including WRKY (72), bHLH (75), AP2 / ERF (37), bZIP (10), CYP450s (114), and UGTs (74). Correlation analysis was performed using the FPKM values ​​of 254 candidate DEGs with FPKM values ​​greater than or equal to 10 in all three samples and the relative contents of 50 key terpenoid flavor compounds (rOAV > 1). The analysis revealed that 41 terpenoids were regulated by 91 mRNAs (|r| > 0.9 and p < 0.01). Among these, CYP450-6, bHLH14, bZIP9, bZIP10, CYP450-27, CYP450-46, UTG52, and GDS7 were regulated by alpha-Pinene, Naphthalene, 1,2,3,5,6,8a-hexahydro-4,7-dimethyl-1-(1-methylethyl)-, (1S-cis)-, and 1,3,6-Octatriene, 3,7-dimethyl-, (Z)-, Humulene, (1R)-2,6,6-Trimethylbicyclo[3.1.1]hept-2-ene, (1S)-2,6,6-Trimethylbicyclo[3.1.1]hept-2-ene, o-Cymene, p-Cymene, beta.-Ocimene, and trans-.beta.-Ocimene were all significantly positively correlated (r > 0.999 and p < 0.01). Figure 4This provides data support for identifying key regulatory genes in the synthesis of terpenoids in pepper leaves.

[0035] Figure 3 (A) Number of all DEGs in the ML_vs_YL, SL_vs_ML and SL_vs_YL comparison groups; (B) Venn diagram of differentially expressed structural genes in the terpene biosynthetic pathway in the three control groups; (C) Heatmap analysis of DEG expression levels in the terpene biosynthetic pathway.

[0036] Figure 4 The data shows associations with correlation coefficients |r| > 0.9 and p < 0.001, where circular nodes represent metabolites and triangular nodes represent genes. Red lines indicate positive correlations (r > 0.9) and green lines indicate negative correlations (r < -0.9).

[0037] Table 3. Number of differentially expressed genes involved in terpene biosynthesis

[0038] 3.2 Analysis of the influence of lncRNAs on the synthesis of terpenoids In addition to mRNAs, a total of 14,226 lncRNAs were identified from RNA sequencing data using a combination of four analytical methods: CPC analysis, CNCI analysis, pfam protein domain analysis, and CPAT analysis. Figure 5 A). The average transcript length of mRNAs is greater than that of lncRNAs, and approximately 80% of lncRNA transcripts are less than 1000 bp in length. Figure 5 B); and approximately 96% of lncRNAs contain only 1-2 exons, a significantly higher proportion than mRNAs (B). Figure 5 C). Furthermore, the FPKM value distribution characteristics of lncRNAs and mRNAs are similar ( Figure 5 D). In the three comparison groups of ML_vs_YL, SL_vs_ML, and SL_vs_YL, 614 (368 upregulated, 246 downregulated), 154 (73 upregulated, 81 downregulated), and 672 (388 upregulated, 284 downregulated) differentially expressed lncRNAs were identified, respectively. Figure 5 E).

[0039] Previous studies have shown that lncRNAs can regulate mRNA expression levels by binding to adjacent target genes, thereby indirectly participating in the regulation of plant life activities. This invention predicts lncRNA target genes based on the complementary base pairing principle between lncRNAs and mRNAs, identifying 145 lncRNAs with potential target genes. Combined with the analysis results of transcripts related to terpene synthesis, further analysis of lncRNAs targeting 91 candidate genes regulating terpene biosynthesis revealed that 145 lncRNAs target 40 candidate key genes, including 11 CYP450s, 6 UTGs, 7 bHLHs, 1 bZIP, and 16 structural genes in metabolic pathways. Figure 5 (F), among which the GDS5 (Zardc21201) gene was targeted and regulated by the most lncRNAs, totaling 52. These results indicate that lncRNAs may play an important indirect regulatory role in the biosynthesis of terpenoids in Zanthoxylum bungeanum by targeting key genes involved in terpenoid biosynthesis.

[0040] Figure 5 (A) Different types of lncRNAs in Sichuan pepper leaves; (BD) Comparison of lncRNAs and mRNAs in terms of transcription length (B), number of exons (C), and FPKM value (D); (E) Number of DE lncRNAs in the ML_vs_YL, SL_vs_ML, and SL_vs_YL control groups; (F) Network diagram of candidate mRNAs targeting terpene biosynthesis.

[0041] 4. Identification of miRNAs and prediction of their target genes For miRNAs, the filtered sequences were compared with known miRNAs using miRDeep2 software (version 2.0.5) to detect new miRNAs. TargetFinder (version 1.6) was used to predict the target genes of all miRNAs. The expression levels of miRNAs in each sample were statistically analyzed, and the expression levels were normalized using the TPM (Transcripts Per Million) algorithm. |log2(FC)| ≥ 1.00 and FDR (False Discovery Rate) ≤ 0.01 were used as the criteria for screening differentially expressed miRNAs (DEmiRNAs).

[0042] In this study, a total of 519 miRNAs were identified through deep sequencing, including 192 known miRNAs and 327 novel predicted miRNAs. Most clean miRNA reads were between 21 and 24 nt in length. Figure 6A), consistent with the typical length characteristics of plant miRNAs. Differential expression analysis showed that 199, 194, and 257 differentially expressed miRNAs were identified in the ML_vs_YL, SL_vs_ML, and SL_vs_YL groups, respectively, for a total of 332 DEmiRNAs. Figure 6 (B) indicates that miRNAs may play an important regulatory role in the physiological processes of Sichuan pepper leaves during the three developmental stages. Of all the miRNAs, 493 miRNAs target 6891 mRNAs (B). Figure 6 C). To elucidate the regulatory mechanism of miRNAs on terpene biosynthesis, and in conjunction with the analysis results of transcripts related to terpene biosynthesis, further analysis of miRNAs was conducted on 91 candidate genes regulating terpene biosynthesis. The results showed that 9 candidate mRNAs (SQLE2, SQLE3, DXS3, LMS3, CYP450-42, CYP450-55, CYP450-81, UTG11, WRKY20) were targeted by 8 miRNAs (novel_miR_13, novel_miR_138, novel_miR_62, csi-miR396b-5p, novel_miR_41, csi-miR171e-3p, csi-miR171g-3p, novel_miR_302), and all were differentially expressed miRNAs. Figure 6 D). Since miRNAs are mainly regulated through mRNA degradation, they show a negative correlation with target gene expression. Specifically, novel_miR_13 is negatively correlated with the expression levels of target genes SQLE2 and SQLE3; novel_miR_138 is negatively correlated with the expression levels of target genes DXS3 and LMS3; csi-miR171e-3p and csi-miR171g-3p are negatively correlated with the expression level of UTG11. Figure 6 D). In summary, these results indicate that miRNAs play a regulatory role in the biosynthesis of terpenoids by targeting key enzyme genes and transcription factors in the terpenoid biosynthesis pathway.

[0043] Figure 6 In the diagram: (A) Length distribution of all miRNA sequences; (B) Number of DE miRNAs in the ML_vs_YL, SL_vs_ML, and SL_vs_YL groups; (C) Number of target genes (mRNAs) of miRNAs; (D) Network diagram of candidate mRNAs for miRNA targeting terpene biosynthesis. Black arrows indicate positively correlated regulation; green "T" symbols indicate negatively correlated regulation.

[0044] 5. Identification and analysis of Circular RNA (circRNA) For circRNAs, the find_circ software was used to predict and identify them. circRNA expression levels were characterized by junction read counts and normalized using the SRPBM (Spliced ​​Reads Per Billion Mapping) normalization method. In the screening of differentially expressed circRNAs (DEcircRNAs), a Fold Change (FC) ≥ 2 and P < 0.05 were used as the criteria for significant difference. Given that circRNA sequences typically contain multiple miRNA binding sites, a miRNA target gene prediction strategy was employed, and potential miRNA-circRNA binding pairs were identified using the Target Finder software.

[0045] circRNAs, as a newly emerging class of non-coding RNA molecules, play a central role in miRNA-mediated gene regulatory networks. By specifically binding to intracellular miRNAs, circRNAs can relieve the inhibitory effect of miRNAs on target genes, thereby playing a crucial regulatory role in important biological processes such as plant growth and development and responses to environmental stress. This invention systematically analyzed the expression characteristics of circRNAs in leaves at different developmental stages, identifying a total of 1434 circRNA molecules. Structurally, exonal and intergenic types were predominant, with a relatively low proportion of intronic types. In terms of length distribution, most circRNA molecules were concentrated in the range of ≥ 3000 bp or < 800 bp. Figure 7 A). Chromosomal localization analysis showed that the largest number of circRNAs (308) were not anchored to chromosomes (unChr), followed by chromosome 2 (chr2, 70) and chromosome 5 (chr5, 59). Figure 7 B). A total of 43 differentially expressed circRNAs (DEcircRNAs) were identified in all comparison groups. Figure 7 C). In the ML_vs_YL comparison group, 28 DEcircRNAs were identified (13 upregulated and 15 downregulated); the SL_vs_ML comparison group had only 2 upregulated DEcircRNAs; while the SL_vs_YL comparison group contained 27 DEcircRNAs (17 upregulated and 10 downregulated). Figure 7 D). Notably, 13 common DEcircRNAs were found in the SL_vs_YL and ML_vs_YL comparison groups. Figure 7E). These results suggest that circRNAs may be closely related to the developmental stage of leaves, and that the expression changes of some circRNAs are consistent between old and young leaves, and between mature and young leaves, suggesting that these circRNAs may play an important role in regulating the synthesis of terpenoids during leaf development.

[0046] 6. Construction of a ceRNA network regulating terpene synthesis To systematically elucidate the regulatory mechanisms of mRNA and non-coding RNA (ncRNA) in the synthesis of terpenoids in Sichuan pepper leaves, based on the ceRNA (endogenous competitive RNA) hypothesis, this invention focuses on the target miRNAs of key candidate mRNAs in terpenoid synthesis, exploring their interactions with circRNAs and lncRNAs, and constructing a complete ceRNA regulatory network. The study found that 437 miRNAs can bind to circRNAs. Further reverse analysis showed that among 43 DEcircRNAs, 32 could be targeted and bound by 197 miRNAs. Among them, seven DEmiRNAs (csi-miR171e-3p, csi-miR171g-3p, novel_miR_138, novel_miR_62, novel_miR_13, novel_miR_41, and novel_miR_302) target and regulate UTG11, DXS3, LMS3, SQLE2, SQLE3, CYP450-42, CYP450-81, and WRKY20 terpene synthesis candidate genes. Figure 7 F). These 7 DEmiRNAs can bind to 92 lncRNAs, involving a total of 8 DElncRNAs.

[0047] Analysis of the constructed ceRNA regulatory network revealed that the network involves a total of 9 mRNAs, 8 miRNAs, 8 lncRNAs, and 14 circRNAs. Figure 7F). Among them, two lncRNAs, MSTRG.29917.1 and MSTRG.17202.11, and circRNA, Superscaffold15:3273121|3315109, can competitively bind to novel_miR_138, indirectly regulating the expression of DXS3 and LMS3 genes; five circRNAs (such as Superscaffold15:3273121|3315109, etc.) indirectly affect the regulation of CYP450-42 gene by binding to novel_miR_62. In addition, 6 circRNAs (Superscaffold29:46805745|46806278, etc.) can bind to noveL_miR_302; 4 lncRNAs (MSTRG.64631.1, etc.) and 3 circRNAs (unanchor7782:2349|22526, etc.) can competitively bind to novel_miR_13; 2 lncRNAs (MSTRG.42768.2 and MSTRG.42198.1) can competitively bind to csi-miR396b-5p; lncRNA (Superscaffold29:46805745|46806278, etc.) can bind to noveL_miR_302; 4 lncRNAs (MSTRG.64631.1, etc.) can bind to noveL_miR_302; 3 lncRNAs (Superscaffold29:46805745|46806278, etc.) can bind to noveL_miR_302; 4 lncRNAs (MSTRG.64631.1 ...4 lncRNAs (MSTRG.64631.1, etc. The circRNA (ld10:72517029|72616745) can competitively bind to two miRNAs (csi-miR171g-3p and csi-miR171e-3p); three lncRNAs (Superscaffold2:78932959|78949713, etc.) can competitively bind to novel_miR_41. These ncRNAs indirectly regulate the expression of key candidate genes for terpene synthesis, such as WRKY20, SQLE2, SQLE2 and SQLE3, CYP450-55, UTG11, and CYP450-81, by competitively binding to miRNAs. These results indicate that circRNAs and lncRNAs can competitively bind to miRNAs, thereby relieving the inhibitory effect of miRNAs on key genes for terpene synthesis and playing an important role in the regulation of terpene compound synthesis in *Zanthoxylum bungeanum* leaves, providing a new perspective for a deeper understanding of the regulatory mechanisms of terpene metabolism.

[0048] Figure 7(A) Sequence length of the identified circular RNAs and their distribution on chromosomes (B); (C) Heatmap analysis of expression levels of all DEcircRNAs; (D) Bar graph showing the number of upregulated and downregulated CircleRNAs in the ML_vs_YL, SL_vs_ML, and SL_vs_YL groups; (E) Venn diagram showing the number of overlapping DEcircRNAs between different control groups; (F) ceRNA network, which may regulate the biosynthesis of terpenoids in Sichuan pepper leaves. Regulatory targeting relationships are represented by nodes of different shapes and connections. All RNAs in the network are differentially expressed.

[0049] 7. qRT-PCR validation of DEmRNAs, DElncRNAs, and DEmiRNAs These three RNA samples were isolated from leaves of *Zanthoxylum bungeanum* at stages 1 (YL), 2 (ML), and 3 (SL). AceQ Universal SYBR qPCR Master Mix and HiScript III RT SuperMix (Vazyme Biotech, Nanjing) were used for qPCR of lncRNAs and mRNAs (+gDNA wiper). Forward-specific primers were designed based on the mature miRNA sequence, while reverse primers used universal primers (via tail A) provided in the miRNA first-strand cDNA synthesis kit. The upstream and downstream primers for the exogenous gene U6 included in the miRNA Unimodal SYBR qPCR Master Mix (Vazyme Biotech, Nanjing) were used as reference primers. qRT-PCR reactions of miRNAs were performed using the miRNA Unimodal SYBR qPCR Master Mix (Vazyme Biotech, Nanjing). All miRNA, lncRNA, and mRNA-specific primers and reference primers are shown in Table 4. qRT-PCR was performed on three independent biological replicas, each containing three technical replicas of each sample. The relative expression level of genes is determined by using 2 -ΔΔCt The method calculates its cycle threshold (Ct) value by normalizing it to the Ct value of the reference gene.

[0050] Table 4. MiRNA, lncRNA, and mRNA-specific primers and reference gene primers.

[0051] To verify the accuracy of the expression profiles in the sequencing data, the expression levels of 9 DE mRNAs, 7 DE miRNAs, and 8 DE lncRNAs were detected by qRT-PCR. The expression changes of the selected mRNAs, lncRNAs, and miRNAs were highly consistent with the trends obtained by qRT-PCR. Figure 8 These results confirm the reliability of the RNA sequencing data.

[0052] Figure 8 In the diagram: Error bars represent the standard errors of three biological replication methods. The bar graph and line graph represent the qRT-PCR and TPM / FPKM values ​​of the gene, respectively. 2 The values ​​represent the correlation between qRT-PCR and TPM / FPKM values. Significance of differences is indicated by different lowercase letters (p < 0.05).

[0053] This invention, through multi-omics analysis of Sichuan pepper leaves, clarified that terpenoids are the main contributors to the aroma of Sichuan pepper leaves. Among the 1468 volatile metabolites detected, 320 were terpenoids, forming the core components of aroma substances. Combining gene expression level and terpenoid content correlation analysis, key candidate regulatory factors, including transcription factors such as WRKY and bHLH, CYP450 and UGTs family genes, and GDS, LMS, HMGS, and DXS, were screened and identified. These factors may participate in the synthesis and regulation of Sichuan pepper terpenoids. Based on these key candidate regulatory factors, a regulatory network involving 9 DEmRNAs, 8 DEmiRNAs, 8 DElncRNAs, and 14 DEcircRNAs was constructed. In this network, lncRNAs and circRNAs competitively bind to miRNAs, targeting and regulating the expression of key candidate regulatory factors, thereby synergistically participating in the synthesis and regulation of Sichuan pepper terpenoids. In summary, this invention clarifies the core role of terpenoids in the aroma formation of Sichuan pepper leaves, elucidates the molecular network regulating their synthesis, and provides a theoretical basis and genetic resources for a deeper understanding of the terpenoid metabolism mechanism in Sichuan pepper and the quality improvement of spice plants.

[0054] The specific embodiments described above further illustrate the purpose, technical solution, and beneficial effects of the present invention. It should be understood that the above description is only a specific embodiment of the present invention and is not intended to limit the scope of protection of the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the scope of protection of the present invention.

Claims

1. A method for revealing the molecular mechanism of terpene biosynthesis in Zanthoxylum bungeanum leaves based on metabolomics and whole transcriptomics analysis, characterized in that, Includes the following steps: S1. Volatile metabolomics analysis: Volatile metabolites were detected by GC-MS. Hierarchical cluster analysis, principal component analysis and K-means cluster analysis were performed on the detected metabolites to screen for differential metabolites. Key flavor compounds with rOAV ≥ 1 were identified based on relative odor activity values. S2. RNA extraction and sequencing: Total RNA was extracted from three groups of leaves (young leaves, mature leaves, and old leaves) of *Zanthoxylum bungeanum* using the RNA Pure Plant Kit. RNA quality was evaluated using NanoDrop 2000 and Agilent Bioanalyzer 2100 systems. Strand-specific RNA-seq libraries and small RNA sequencing libraries were constructed, and high-throughput sequencing was performed to obtain raw sequencing data, which were then processed using the BMKCloud platform. Analysis of S3, mRNA, lncRNA, miRNA, and circRNA: The sequencing data from step S2 were analyzed for mRNA, lncRNA, miRNA, and circRNA. S4. Construction of the regulatory network for terpene synthesis: Based on the key terpenoid flavor compounds identified in step S1 and the differentially expressed RNAs screened in step S3, key candidate genes involved in terpenoid synthesis were screened through correlation analysis; based on the ceRNA hypothesis, the interaction between key candidate genes and DEmiRNAs, DElncRNAs, and DEcircRNAs was analyzed to construct a complete ceRNA regulatory network. S5. Verification: The expression levels of DEmRNAs, DElncRNAs, and DEmiRNAs screened in step S3 were validated using qRT-PCR to confirm the reliability of the sequencing data.

2. The method for revealing the molecular mechanism of terpene biosynthesis in Zanthoxylum bungeanum leaves based on metabolomics and whole transcriptomics analysis according to claim 1, characterized in that, The young leaves are the leaves of the Sichuan pepper plant (YL) 10 days after germination, the mature leaves are the leaves of the Sichuan pepper plant (ML) 60 days after germination, and the old leaves are the leaves of the Sichuan pepper plant (SL) 120 days after germination.

3. The method for revealing the molecular mechanism of terpene biosynthesis in Zanthoxylum bungeanum leaves based on metabolomics and whole transcriptomics analysis according to claim 1, characterized in that, In step S1, the GC-MS used was an Agilent Model 8890-7000D, and each sample was subjected to three biological replicates. In the differential metabolite screening, the DTMs of the ML_vs_YL, SL_vs_ML, and SL_vs_YL groups were 721, 136, and 714, respectively, and were enriched to the secondary metabolite synthesis, monoterpene biosynthesis, sesquiterpene, and triterpene biosynthesis pathways through metabolic pathway annotation.

4. The method for revealing the molecular mechanism of terpene biosynthesis in Zanthoxylum bungeanum leaves based on metabolomics and whole transcriptomics analysis according to claim 1, characterized in that, In step S2, when constructing the strand-specific RNA-seq library, the Ribo-Zero rRNA Removal Kit was used to remove ribosomal RNA, and the NEBNext® Ultra™ Directional RNALibrary Prep Kit was used to construct the strand-specific library; when constructing the miRNA library, the TruSeq Small RNA SamplePrep Kit was used.

5. The method for revealing the molecular mechanism of terpene biosynthesis in Zanthoxylum bungeanum leaves based on metabolomics and whole transcriptomics analysis according to claim 1, characterized in that, The mRNA and lncRNA analysis described in step S3 is as follows: Hisat2 is used to align sequencing sequences to a reference genome; transcripts > 200 nt in length and containing ≥ 2 exons are screened as lncRNA candidates; lncRNAs are identified using four computational methods: CPC2, CNCI, Pfam, and CPAT; StringTie is used to calculate the FPKM values ​​of mRNA and lncRNA; differentially expressed mRNAs and lncRNAs are screened using the DESeq R software package with FC ≥ 2 and Pvalue < 0.01 as criteria; cis-target genes of lncRNAs are predicted based on positional relationships, and trans-target genes are predicted based on expression correlations; functional information of the above target genes is annotated using six public databases: Nr, Pfam, KOG / COG, Swiss-Prot, KO, and GO.

6. The method for revealing the molecular mechanism of terpene biosynthesis in Zanthoxylum bungeanum leaves based on metabolomics and whole transcriptomics analysis according to claim 1, characterized in that, The miRNA analysis described in step S3 involves using miRDeep2 software to identify known and new miRNAs, using TargetFinder to predict target genes, normalizing expression levels using the TPM algorithm, and screening differentially expressed miRNAs based on |log2(FC)| ≥ 1.00 and FDR ≤ 0.

01.

7. The method for revealing the molecular mechanism of terpene biosynthesis in Zanthoxylum bungeanum leaves based on metabolomics and whole transcriptomics analysis according to claim 1, characterized in that, The circRNA analysis described in step S3 involves identifying circRNAs using the find_circ software, characterizing expression levels by junction read counts and normalizing them using SRPBM, screening differentially expressed circRNAs based on FC ≥ 2 and P < 0.05, and identifying miRNA-circRNA binding pairs using the Target Finder software.

8. The method for revealing the molecular mechanism of terpene biosynthesis in Zanthoxylum bungeanum leaves based on metabolomics and whole transcriptomics analysis according to claim 1, characterized in that, The regulatory network described in step S4 includes 9 DEmRNAs, 8 DEmiRNAs, 8 DElncRNAs and 14 DEcircRNAs, wherein lncRNAs and circRNAs competitively bind to miRNAs, thereby relieving the inhibitory effect of miRNAs on key genes for terpene synthesis.

9. The method for revealing the molecular mechanism of terpene biosynthesis in Zanthoxylum bungeanum leaves based on metabolomics and whole transcriptomics analysis according to claim 1, characterized in that, In step S5, the qRT-PCR validation used AceQ UniversalSYBR qPCR Master Mix and HiScript III RT SuperMix to detect lncRNA and mRNA, and miRNA Universal SYBR qPCR Master Mix to detect miRNA, using U6 as a reference gene. The assay was performed using 2... -ΔΔCt The method calculates relative expression levels, and the validation results are consistent with the trends in sequencing data.

Images