A method for mining key genes for the synthesis of terpenoid compounds in Zanthoxylum bungeanum leaves based on third-generation full-length and second-generation transcriptome sequencing

Abstract

The present invention discloses a method for mining key genes for the synthesis of terpenoid compounds in Zanthoxylum armatum leaves based on third-generation full-length and second-generation transcriptome sequencing, including: determining the components and contents of terpenoid compounds in the leaves of Zanthoxylum armatum DC. var. ferrugineum Rehd. and its bud mutant variety Rongchang thornless Zanthoxylum armatum at three different developmental stages; extracting total RNA from the leaves to construct a sequencing library, performing data analysis and functional annotation. CDS, SSR, lncRNAs and As events were predicted, differentially expressed genes were analyzed, and differentially expressed genes of structural genes and transcription factors related to terpenoid biosynthesis were screened; a co-expression gene network was constructed, and the FPKM values of hub genes were analyzed for correlation with the contents of terpenoid compounds, and genes positively correlated with the contents of terpenoid compounds were selected as genes related to their synthesis. The present invention successfully constructed a high-quality transcriptome of Zanthoxylum armatum, revealing 12 key genes that may be involved in regulating the synthesis of terpenoid compounds, laying a foundation for analyzing the molecular mechanism of terpenoid biosynthesis in Zanthoxylum armatum leaves.

Description

Technical Field

[0001] The present invention relates to the technical fields of plant molecular biology and plant genetic engineering, and particularly relates to a method for mining key genes for the synthesis of terpenoid compounds in Zanthoxylum bungeanum leaves based on third-generation full-length and second-generation transcriptome sequencing. Background Art

[0002] Zanthoxylum bungeanum is a deciduous small tree of the genus Zanthoxylum in the Rutaceae family and is an important economic crop tree species, widely planted in Asian countries such as China, South Korea, and Japan. There are approximately 250 species of Zanthoxylum in the world, and Zanthoxylum armatum is one of the main cultivated varieties in China and has a long cultivation history in China. Zanthoxylum armatum is also known as Sichuan green pepper. Its fruits, peels, leaves, bark and other tissues are rich in essential oils and have a strong spicy aroma. Terpenoid compounds are the main flavor components of Zanthoxylum bungeanum and are widely used in medicine and food, having important economic value.

[0003] There are two biosynthetic pathways of terpenoid compounds in plants, namely the mevalonic acid (MVA) and 2-methyl-D-erythritol 4-phosphate (MEP) pathways. Pharmacological studies have shown that terpenoid compounds have various biological activities such as anti-inflammatory, antibacterial, antioxidant, and cytotoxicity, thus attracting extensive attention from researchers. In biological control, terpenoid compounds can act on the insect nerve cell membrane receptor to change the ion channel and inhibit the nerve conduction of insects, thereby poisoning and repelling insects. For humans, terpenoid compounds can be used to improve diabetes, promote the sensitivity of cervical cancer cells to chemotherapy drugs, and scavenge free radicals. It has been found that there are significant differences in the essential oil components and contents of Zanthoxylum bungeanum leaves at different developmental stages. Therefore, according to the contents of different terpenoid compounds in the leaves of Zanthoxylum bungeanum at different developmental stages, the leaf developmental stage for essential oil extraction is determined to achieve the maximum utilization rate of Zanthoxylum armatum leaves and improve the yield of Zanthoxylum bungeanum.

[0004] The related research on Zanthoxylum bungeanum leaves mainly focuses on the identification and determination of active substances such as antioxidants, flavonoids, lignans, and amides. However, there has been no report on using SMRT-seq technology to screen candidate genes to determine the molecular mechanism of terpenoid biosynthesis in Zanthoxylum bungeanum leaves. Zanthoxylum armatum has a very complex genetic background and a large genome. At present, the lack of high-quality genome data has hindered the further research on this plant. Summary of the Invention

[0005] The object of the present invention is to provide a method for mining key genes for the synthesis of terpenoid compounds in Zanthoxylum bungeanum leaves based on third-generation full-length and second-generation transcriptome sequencing. The present invention uses Zanthoxylum armatum (ZY) and its bud mutant variety Rongchang thornless Zanthoxylum bungeanum (WC) as research materials, combines transcriptome sequencing (third-generation and second-generation) with the determination of terpenoid compounds for analysis, to identify the key genes for the synthesis of terpenoid compounds in Zanthoxylum armatum, and to clarify the biosynthetic pathway and regulatory mechanism of terpenoid compounds, improve the Zanthoxylum bungeanum transcriptome data resources, provide gene resources for the later research on the cloning, expression analysis and functional research of key enzyme genes in the terpenoid compound synthesis and metabolism pathway of Zanthoxylum bungeanum leaves, at the same time lay a foundation for analyzing the molecular mechanism of the biosynthesis of terpenoid compounds in Zanthoxylum bungeanum leaves, further strengthen the development and utilization of Zanthoxylum bungeanum leaves, and also provide genetic data for using genetic engineering technology to increase the content of terpenoid compounds in Zanthoxylum bungeanum leaves.

[0006] To achieve the above object, the technical solution adopted by the present invention is:

[0007] A method for mining key genes for the synthesis of terpenoid compounds in Zanthoxylum bungeanum leaves based on third-generation full-length and second-generation transcriptome sequencing, characterized in that the method comprises:

[0008] Step 1, determining the components and contents of terpenoid compounds in three different development stages of Zanthoxylum armatum and Rongchang thornless Zanthoxylum bungeanum, namely young leaves, mature leaves and old leaves, and comparing the differences;

[0009] Step 2, respectively extracting the total RNA of young leaves, mature leaves and old leaves of Zanthoxylum armatum and Rongchang thornless Zanthoxylum bungeanum, and extracting the total RNA of the mixed samples of young leaves, mature leaves and old leaves of Zanthoxylum armatum and Rongchang thornless Zanthoxylum bungeanum, sequencing respectively and constructing a full-length transcriptome and a transcriptome dataset to obtain raw data;

[0010] Step 3, correcting the raw data obtained in Step 2 to obtain high-quality sequences, filtering the high-quality sequences to obtain clean data, and performing assembly processing on the clean data to obtain the full-length transcriptome database and transcriptome database of Zanthoxylum armatum and Rongchang thornless Zanthoxylum bungeanum;

[0011] Step 4, predicting alternative splicing candidate events and simple sequence repeats in the full-length transcriptome database and transcriptome database obtained in Step 3, and verifying the authenticity of the alternative splicing candidate events;

[0012] Step 5, predicting the coding region sequences and lncRNAs of the transcript sequences in the full-length transcriptome database and transcriptome database obtained in Step 3;

[0013] Step 6: Perform functional annotation on the genes in the transcriptome database in Step 3 and predict transcription factors, identify differentially expressed genes among the samples, and screen for enzyme genes and transcription factors related to terpenoid synthesis based on the functional annotation and differential gene information;

[0014] Step 7: Predict the target genes of the lncRNAs in Step 5, and perform regulatory network analysis on the lncRNAs of the genes screened in Step 6 according to the lncRNA target gene prediction information;

[0015] Step 8: Perform weighted gene co-expression network analysis on the genes in Step 6, construct a gene regulatory network diagram through visual network construction, and screen for hub genes with high connectivity to terpenoid metabolism;

[0016] Step 9: Perform correlation analysis on the significantly different terpenoids in Step 1 and the hub genes in Step 8, and select the genes that are correlated with the terpenoid content as the key genes for terpenoid synthesis.

[0017] Furthermore, in Step 1, a gas chromatography-mass spectrometry (GC-MS) instrument is used to determine the components and relative contents of terpenoids in the young leaves, mature leaves, and old leaves of Zanthoxylum armatum and Rongchang thornless Zanthoxylum.

[0018] Furthermore, in Step 2, the PacBio RS II real-time sequencing platform and the Illumina Hiseq4000 platform are respectively used for full-length transcriptome and transcriptome sequencing and to construct a sequencing data set.

[0019] Furthermore, in Step 3, the original sequences are converted into ROI (Reads of Insert) sequences to obtain full-length sequences and non-full-length sequences. The ICE (Iterative isoform-clustering) algorithm is used to perform clustering analysis on the RoI (Reads of Insert) sequences from the same transcript to obtain consensus sequences, and the non-full-length sequences are used to correct the obtained consensus sequences to obtain high-quality sequences.

[0020] Furthermore, in Step 3, the method for filtering the high-quality sequences is as follows: remove the data with rRNA repeats and adapter-containing data, and remove the low-quality data with a quality value Q ≤ 20 and more than 50% bases to obtain clean data.

[0021] Further, in step 4, the BLAST software is used to perform pairwise alignment on all sequences in the full-length transcriptome database and the transcriptome database, and the alignment results meet the following conditions: 1) The lengths of both sequences are greater than 1000 bp, and there are two HSPs in the alignment; 2) The alternative splicing Gap is greater than 100 bp and is at least 100 bp away from the 3' / 5' end; 3) Sequences allowing 5-bp overlap of all alternative transcripts are considered candidate alternative splicing events.

[0022] Further, in step 4, total RNAs of Zanthoxylum armatum and Z. simulans var. inermis are used to prepare cDNAs, specific primers are designed based on the open reading frames of the sequencing database, and the authenticity of alternative splicing candidate events is verified by reverse transcription polymerase chain reaction.

[0023] Further, in step 4, transcripts longer than 500 bp are selected, and the MISA software (http: / / pgrc.ipk-gatersleben.de / misa / ) is used to predict simple sequence repeats (SSRs).

[0024] Further, in step 5, the TransDecoder software is used to predict the coding region sequences of transcript sequences;

[0025] Further, in step 5, four methods, namely CPC analysis, CNCI analysis, pfam protein domain analysis, and CPAT analysis, are used to predict lncRNAs. Transcripts with a length exceeding 200 nt and having more than two exons are selected as lncRNA candidates.

[0026] Further, in step 6, the genes in the transcriptome database of Zanthoxylum are subjected to sequence alignment and functional annotation in the Nr, Swiss-Prot, Pfam, KEGG, GO, and COG databases using the BLASTN software, where E-value < 10 -5 ; then the iTAK software is used to predict and classify transcription factors in the database; then the DESeq2 analysis is used to analyze the differentially expressed genes in the leaves at three different developmental stages of Z. armatum and Z. simulans var. inermis. The Benjamini-Hochberg method is used to correct the significant p-values obtained from the original hypothesis test, and a false discovery rate less than 0.01 and a fold change ≥ 2 are used as the screening criteria for differentially expressed genes; structural genes and transcription factors related to terpene compound synthesis are screened based on the annotation information and differentially expressed gene information.

[0027] Further, in step 7, the OmicShare (https: / / www.omicshare.com / tools / Home / Soft / cytoscape) was used to analyze and display the regulatory relationship between lncRNAs and candidate genes.

[0028] Further, in step 8, the Biomarker Cloud Platform (https: / / international.biocloud.net / zh / software / tools / detail / small / 8a8300b253cf73e70153d16368250f32) was used to perform weighted gene co-expression network analysis (WGCNA), and the OmicShare (https: / / www.omicshare.com / tools / Home / Soft / cytoscape) was used to construct a visual network to construct a gene regulatory network (GRN) to identify hub genes with high connectivity to terpenoid metabolism.

[0029] Further, the method further includes: verifying the full-length transcriptome and transcriptome sequencing results, including: randomly selecting the lncRNAs predicted in step 5 and the differentially expressed genes screened in step 6, using qRT-PCR technology to measure the expression levels of lncRNAs and differentially expressed genes, and calculating their correlation with the sequencing FPKM value data. If the expression levels of lncRNAs, the expression levels of differentially expressed genes, and the sequencing FPKM value data are correlated, the sequencing results are accurate and reliable.

[0030] The present invention also provides an application of the above method in screening genes related to the biosynthesis of terpenoid compounds in Zanthoxylum bungeanum leaves.

[0031] The present invention also provides genes related to the biosynthesis of terpenoid compounds in Zanthoxylum bungeanum leaves mined by the above method, and the genes include: AACT1, AACT2, HMGR4, AP2 / ERF5, AP2 / ERF6, bZIP1, bZIP3, bZIP5, HDR2, HDR3, HDR4, and bHLH10.

[0032] Compared with the prior art, the beneficial effects of the present invention are as follows: Based on the combined analysis method of full-length transcriptome and transcriptome sequencing, the present invention has successfully constructed the full-length transcriptome and transcriptome databases of Zanthoxylum armatum and its bud mutation variety Rongchang thornless Zanthoxylum armatum. Based on the sequencing data and the determination and analysis of terpenoids, 31 terpenoids have been identified and 12 key genes that may be involved in the synthesis of terpenoids in Zanthoxylum armatum leaves have been mined, providing new ideas for analyzing the biosynthesis of terpenoids in Zanthoxylum armatum leaves, laying a theoretical foundation for improving the content of terpenoids in Zanthoxylum armatum leaves through genetic engineering means, further strengthening the development and utilization of Zanthoxylum armatum leaves, and ultimately increasing the benefits of Zanthoxylum armatum. BRIEF DESCRIPTION OF THE DRAWINGS

[0033] Figure 1 It is the statistical result of ORF and SSR information in Example 1 of the present invention. Figure 1 -a is the quantity and length distribution of the protein sequences encoded by complete ORFs. Figure 1 -b is the SSR statistics in the transcript database.

[0034] Figure 2 It is the analysis of lncRNA in Example 1 of the present invention, where Figure 2 -a is the Venn diagram of the number of lncRNAs predicted by four analysis methods. Figure 2 -b is the analysis of the expression pattern of lncRNA.

[0035] Figure 3 It is the statistical classification of NR annotation species in Example 1 of the present invention.

[0036] Figure 4 It is the statistical result of the number of types of transcription factor families in Example 1 of the present invention.

[0037] Figure 5 It is the statistical result of DEGs and transcript annotation information among samples in Example 1 of the present invention.

[0038] Figure 6 It is the functional annotation result of DEGs in Example 1 of the present invention, where Figure 6 -a is the statistical classification of Go annotation. Figure 6 -b is the scatter plot of KEGG enrichment.

[0039] Figure 7 It is the qRT-PCR verification of DEGs in Example 1 of the present invention.

[0040] Figure 8 It is the PCR verification of AS events in Example 1 of the present invention.

[0041] Figure 9This is the WGCNA analysis of the leaves of Zanthoxylum armatum and Rongchang thornless Zanthoxylum bungeanum at different developmental stages in Example 1 of the present invention, where Figure 9 -a is the hierarchical clustering tree of co-expression modules based on WGCNA analysis, Figure 9 -b is the correlation analysis between transcription modules and tissue samples;

[0042] Figure 10 is the regulatory network diagram of candidate lncRNAs and corresponding target genes involved in the regulation of terpenoid synthesis in Example 1 of the present invention, where Figure 10a is the regulatory network diagram of candidate lncRNAs involved in the regulation of terpenoid synthesis and target genes MDS1, HDR4, HDR3, AP2 / ERF32, AP2 / ERF34, AP2 / ERF31, AP2 / ERF28, DXR6, DXR1, DXR2, bHLH11, HDR2, Figure 10b is the regulatory network diagram of candidate lncRNAs involved in the regulation of terpenoid synthesis and target genes HDS6, bHLH28, DXR11, DXR13, AP2 / ERF63, DXR12, DXR15, WRKY37, WRKY35, bHLH29, HMGR4, bZIP18, bZIP19, AP2 / ERF35, bZIP1, bZIP3, bZIP5, AACT1, AACT2, AP2 / ERF6, AP2 / ERF5, Figure 10c is the regulatory network diagram of candidate lncRNAs involved in the regulation of terpenoid synthesis and target genes AACT8, HMGS4, HMGS6, bHLH10, AP2 / ERF36, AP2 / ERF38;

[0043] Figure 11 This is the expression heat map and transcriptional regulatory network related to the tissue-specific module in Example 1 of the present invention, where Figure 9 -a, 9-b, 9-c, and 9-d in the heat map are significantly overexpressed in ZY-I, WC-Y, ZY-O, and WC-O for the genes in this module, respectively;

[0044] Figure 12 This is the correlation analysis between the expression levels of 40 genes and the relative contents of 18 terpenoids in Example 1 of the present invention. Detailed implementation manners

[0045] Next, the technical solutions of the present invention will be clearly and completely described in conjunction with the embodiments in the present invention. Obviously, the described embodiments are only a part of the embodiments of the present invention, rather than all of the embodiments. All other embodiments obtained by those of ordinary skill in the art based on the embodiments in the present invention without creative efforts belong to the scope of protection of the present invention.

[0046] Example 1

[0047] This example provides a method for mining key genes for the synthesis of terpenoid compounds in Zanthoxylum bungeanum leaves based on third-generation full-length and second-generation transcriptome sequencing, which is as follows:

[0048] 1. Determination of the components and contents of terpenoid compounds in Zanthoxylum bungeanum leaves by gas chromatography-mass spectrometry

[0049] Weigh 0.4 g of fresh leaves of Zanthoxylum armatum (ZY) and Rongchang thornless Zanthoxylum bungeanum (WC) at three different developmental stages (young leaves I, mature leaves M, and old leaves O) respectively, grind them into powder under liquid nitrogen, add 1 ml of methyl tert-butyl ether (containing tetradecane as an internal standard), shake for 2 hours, centrifuge at 15,000 rpm for 10 minutes, and finally collect the supernatant as the sample for injection. Use gas chromatography-mass spectrometry (GC-MS) to detect the components and contents of terpenoid compounds in different samples. The chromatographic column is Agilent DB-5MS (30 m × 0.25 mm, 0.25 μm), and the program is as follows: 40°C for 2 minutes, heat up to 220°C at a rate of 5°C / min, hold at 220°C for 5 minutes, heat up to 280°C at a rate of 5°C / min, and hold at 280°C for 5 minutes. Helium is used as the carrier gas with a flow rate of 1.0 ml / min, the injection volume is 2 μl, and the injector temperature is set at 250°C. The transfer line temperature and ion source temperature are set at 250°C and 230°C respectively. Three biological replicates and six technical replicates are set for analysis. Based on the measurement results of the internal standard and the data in the NIST mass spectrometry library, use R software (http: / / cran.r-project.org / ) and TagFinder to process the chromatographic data, and determine the statistical significance through Duncan's multiple range test. Determine the relative contents and components of terpenoid compounds in Zanthoxylum bungeanum leaves, and the measurement results are shown in Table 1.

[0050] Table 1 Components and relative contents of terpenoid compounds in leaves of Zanthoxylum armatum (ZY) and Rongchang thornless Zanthoxylum bungeanum (WC) at three different developmental stages (young leaves I, mature leaves M, old leaves O)

[0051]

[0052]

[0053] Note: The data are the average values of biological replicates (n = 6). Different letters at the end of the same row of data indicate significant differences in the relative contents of terpenoid compounds in different samples (P < 0.05).

[0054] According to the measurement results in Table 1, a total of 31 major compounds were identified, including monoterpenoids, sesquiterpenoids, diterpenoids, triterpenoids, and other compounds. The results showed that sesquiterpenoids were the main terpene compound type in ZY, with the highest content of elemene, and the relative content of sesquiterpenoids in ZY was higher than that in WC; monoterpenoids were the main terpene compound type in WC, with the highest content of linalool, and the relative content of monoterpenoids in WC was higher than that in ZY. Obviously, the substance contents were different between the two varieties and among different organs, and the total terpene content in ZY was significantly higher than that in WC.

[0055] 2. RNA Extraction, Library Preparation, and Transcriptome Sequencing

[0056] (1) Extraction of Total RNA from Zanthoxylum bungeanum

[0057] Taking 36 Zanthoxylum armatum DC. (ZY) and 36 Zanthoxylum simulans Hance var. inermis (WC) as research objects, with every six plants as one replicate, for a total of six replicates. Pest-free leaf samples at three different developmental stages, young leaves (10 days after germination, WC-I, ZY-I), mature leaves (60 days after germination, WC-M, ZY-M), and old leaves (120 days after germination, WC-O, ZY-O), were collected, frozen with liquid nitrogen, and stored in a -80 °C refrigerator for later use.

[0058] The three different sample types at the three stages of the six replicates were mixed in equal amounts to obtain the mixed analysis samples of WC and ZY for full-length transcript sequencing. Transcriptome sequencing libraries were constructed using the total RNA of leaves at three different developmental stages of WC and ZY. With every six plants as one replicate, every two replicates were combined into a composite sample as a biological replicate, for a total of three biological replicates. The total RNA of the mixed samples of WC and ZY and the total RNA of leaves at three different developmental stages of WC and ZY were extracted using the TaKaRa MiniBEST Plant RNA Extraction Kit (TaKaRa, China). The purity of the RNA samples was measured using a NanoDrop 2000 microspectrophotometer, and the integrity and concentration of the RNA samples were detected using an Agilent 2100 RNA Nano 6000 bioanalyzer. OD 260 / OD 280 values between 1.6 - 1.8 were better for subsequent experiments.

[0059] (2) Full-Length Transcriptome Library Preparation and Transcriptome Sequencing

[0060] After the total RNA passed the test in step (1), it was used to construct a full-length transcriptome library and a transcriptome library. The library construction and sequencing were completed by Biomarker Technologies Corporation, and the raw sequences (raw reads) were obtained by sequencing. The raw sequences of the third-generation full-length sequencing were converted into RoI (Reads of Insert) sequences to obtain full-length sequences and non-full-length sequences; the RoI sequences from the same transcript were clustered using the ICE (Iterative isoform-clustering) algorithm, and the RoI with similar sequences were clustered into a cluster, and a consensus sequence was obtained for each cluster; the consensus sequence obtained was corrected (polishing) using the non-full-length sequences to obtain high-quality sequences; the sequences with differences only in the 5'-terminal exons of the high-quality sequences and the same other exons were merged, and the longest sequence among them was taken as the unigene transcript, that is, the final transcript sequence. The sequences with high similarity were merged using the CD-HIT software, and the redundant sequences in the high-quality transcripts were removed using CD-HIT to obtain 148,445 non-redundant transcript sequences. The rRNA repeats, reads containing adapters, poly-N, and low-quality data with a quality value Q≤20 and more than 50% bases were removed from the original data through an internal Perl script to obtain clean reads, and a total of 117.65G of clean data (Clean Data) was obtained.

[0061] 3. Prediction of coding sequence (CDS), simple sequence repeat (SSR) and lncRNA

[0062] (1) CDS prediction

[0063] The TransDecoder software was used to predict the coding sequence and its corresponding amino acid sequence of the transcript sequence. Based on information such as the length of the open reading frame (Open Reading Frame, ORF), the log-likelihood score, and the ratio of the amino acid sequence to the protein domain sequence in the Pfam database, reliable potential coding sequences (Coding Sequence, CDS) were identified from the transcript sequence. A total of 143,122 CDSs were predicted in WC and ZY, including 114,205 complete CDSs. The length distribution of the protein sequences encoded by the predicted complete CDSs is as Figure 1 shown in a. The proteins encoded by the open reading frame contain 0-2,700 amino acids, and the proteins containing 100-200 amino acids (22,591, 19.78%) are the most, followed by the proteins containing 200-300 amino acids (20,749, 18.17%).

[0064] (2) Prediction of SSRs

[0065] Transcripts longer than 500 bp were selected and the frequency and distribution of SSRs were analyzed using the MISA tool (http: / / pgrc.ipk-gatersleben.de / misa / ). MISA identified seven types of SSRs, namely mononucleotide, dinucleotide, trinucleotide, tetranucleotide, pentanucleotide, hexanucleotide, and compound SSRs. SSR analysis was performed on 145,505 sequences, and the results are as Figure 1 b. A total of 105,465 SSRs were detected. In addition, 66,239 sequences contained SSRs, and 25,256 of these sequences contained more than one SSR. Among all the SSRs, mononucleotides were found to be the most abundant (48,877), followed by trinucleotides (13,688), dinucleotides (10,700), hexanucleotides (980), tetranucleotides (946), and pentanucleotides (249).

[0066] (3) Prediction of lncRNAs

[0067] Four analytical methods, namely CPC analysis, CNCI analysis, pfam protein domain analysis, and CPAT analysis, were combined to screen the coding potential of transcripts for predicting lncRNAs. Transcripts encoding more than 100 amino acids were filtered using the minimum length and exon number threshold method. Transcripts longer than 200 nucleotides and with more than two exons were selected as lncRNA candidates. The results are as Figure 2 shown in a. A total of 4,719 lncRNAs were predicted among the sample tissues.

[0068] (4) qRT-PCR verification of lncRNAs

[0069] Six lncRNAs were randomly selected from the lncRNAs predicted in step (3). Specific primers were designed using Primer 5.0 software based on the transcript sequences. Dual internal reference genes were used, and the primer sequences of the genes are shown in Table 2. Sample RNA was reverse transcribed into cDNA using the Hiscript II QRT SuperMix for qPCR (+gDNA wiper) kit. The cDNA was diluted 10-fold and used as the qPCR template. qPCR was performed using the BioEasy Master Mix (SYBR Green, High Rox) kit. Real-time fluorescence quantitative PCR was carried out on a Bio-Rad Mini OpticonTM real-time PCR mini cycler. The operation steps are as follows:

[0070] ①Genomic DNA removal: 4×gDNA wiper Mix, 4 μL; Total RNA, 1 pg - 1 μg; RNase Free dH2O, To 16 μL. Pipette gently to mix well. Reaction program: 42°C for 2 min.

[0071] ②Reverse transcription reaction: 5×qRT SuperMix II, 4 μL; Reaction solution from step ①, 16 μL. Pipette gently to mix well. Reaction program: 50°C for 15 min, 85°C for 2 min, store the template cDNA at -20°C.

[0072] ③qPCR reaction: Reaction system 20 μL: 2×ChamQ Universal SYBR qPCR Master Mix, 10.0 μL; Forward Primer (10 uM), 0.4 μL; Reverse Primer (10 uM), 0.4 μL; Templete cDNA, 2 μL; ddH2O, To 20 μL. Reaction program: 95°C, 1 min; 95°C, 15 s; 60°C, 1 min, for a total of 35 cycles, add melting curve program, 95°C, 1 min; 65°C, 1 min; 95°C, 20 s; 30°C, 1 min.

[0073] ④Each sample was performed with three biological replicates and three technical replicates, using ultrapure water as a negative control, and two internal references were used. The internal reference genes were reference gene 1 (Actin) and reference gene 2 (Ubiquitin - conjugating enzyme). The 2 -ΔΔCT method was used to convert the CT value of each amplification into the relative expression level for analyzing the correlation between qRT - PCR and RNA - seq results.

[0074] The expression of lncRNA is as Figure 2 shown in b, and the correlation between the qRT - PCR results and FPKM values of each lncRNA was analyzed. The results showed that the qRT - PCR results of lncRNA were significantly correlated with the FPKM values (R 2 > 0.5391), indicating the reliability of the sequencing data.

[0075] Table 2 Sequence information of designed primers

[0076]

[0077]

[0078]

[0079] 4. Functional annotation and identification of transcription factors (TFs)

[0080] (1) Functional annotation

[0081] Annotation information was obtained by using BLASTX (E-value ≤ 1e-5). The 148,445 transcript sequences obtained in step 2 were aligned to eight public databases, namely Nr (Non-redundant database), Swiss-Prot, KEGG (Kyoto Encyclopedia of Genes and Genomes), GO (Gene Ontology), and COG (Clusters of Orthologous Group). A total of 142,829 transcripts had annotation information. 142,829 transcripts were successfully annotated, and the annotation results are shown in Table 3. By aligning to the non-redundant protein (Nr) database, the species distribution of the best-matching results was obtained ( Figure 3 ). The transcripts showed the closest match to Citrus sinensis (69,539, 48.74%), followed by Citrus clementina (54,193, 37.99%) and Theobroma cacao (2,070, 1.45%), all belonging to the Rutaceae family.

[0082] Table 3 Statistics of the number of annotated transcripts

[0083]

[0084] (2) Identification of transcription factors

[0085] Transcription factors (TFs) were predicted by iTAK (http: / / bioinfo.bti.Cornell.edu / tool / ITak), and different types of TFs were classified and counted. In all tissue samples, a total of 20,145 TFs were predicted and divided into different transcription factor families. The results are shown in Figure 4 and show the top 20 families with the most TF members. The family with the most members is RLK-Pelle_DLSV, with a total of 723 gene members; followed by C3H transcription factors (693), bHLH (565), and GRAS (539).

[0086] 5. Differential expression gene analysis (DEGs)

[0087] (1) Identification of DEGs

[0088] Differential expression analysis between sample groups was performed using DESeq2. During the differential expression analysis, the recognized and effective Benjamini-Hochberg method was used to correct the significant p-values obtained from the original hypothesis test to reduce the false positives brought about by independent statistical hypothesis tests on the expression values of a large number of genes. During the screening process, a false discovery rate (FDR) of less than 0.01 and a fold change (FC) of greater than or equal to 2 were used as the screening criteria.

[0089] The results are as Figure 5 shown. A total of 96,931 DEGs were identified in the leaves at three different developmental stages of Zanthoxylum armatum and Z. armatum var. inermis. The number of DEGs was the largest among different varieties at the same developmental stage. Between WC-M_vs_ZY-M, 70,800 DEGs were identified ( Figure 5 e); between WC-I_vs_ZY-I, 65,781 DEGs were identified ( Figure 5 c); between WC-O_vs_ZY-O, 63,225 DEGs were identified ( Figure 5 f).

[0090] (2) Annotation of DEGs

[0091] The 96,931 DEGs obtained in step (1) were aligned with the COG, GO, KEGG, KOG, Pfam, Swiss-Prot, eggNOG, and Nr databases to obtain DEG annotation information. The statistical results of the functional annotation information of DEGs between different tissues are as Figure 5 shown in j. Among them, 65,206 DEGs were successfully annotated to 49 types in GO. In the biological process category, the DEGs annotated to metabolic process, cellular process, and single-organism process were the most; in the cellular component category, the DEGs annotated to cell part, cell, and organelle were the most; in the molecular function category, the DEGs annotated to catalytic activity and binding were the most (Figure 6a). The results of the KEGG pathway enrichment analysis of DEGs are as Figure 6 shown in b. A total of 26,606 DEGs were annotated to the KEGG pathways, and the pathways of photosynthesis (ko00195), DNA replication (ko03030), and homologous recombination (ko03440) were significantly enriched (q < 0.1).

[0092] (3) Screening of candidate genes involved in the synthesis of terpenoids in Zanthoxylum leaves

[0093] Combined with the gene function annotation information in Step 4 and the enzyme gene information in the biosynthesis and metabolism of terpenoids in plants, using the enzyme EC number and enzyme gene name as the search basis, gene sequences were screened, and gene sequences with incomplete ORFs were further excluded. In addition, through BLAST alignment on the NCBI website, it was verified whether the screened genes were the corresponding enzyme genes searched for. The statistics of the initially screened gene numbers are shown in Table 4, and 36, 84, and 21 genes involved in the MVA, MEP pathways, and downstream branch points were screened respectively.

[0094] Table 4 Statistics of the number of genes involved in the biosynthesis of terpenoids

[0095]

[0096]

[0097] In addition, four TFs that may be related to the synthesis of terpenoids were screened, namely bHLH (103), bZIP (62), AP2 / ERF (159), and WRKY (95). Combining the genes that may be involved in terpenoid biosynthesis, a total of 560 genes (FPKM≥10) were screened. Combining the information of 65,206 DEGs in Step (1), 526 of the 560 genes were DEGs and were used as the objects for subsequent analysis.

[0098] (4) Analysis of the expression patterns of DEGs

[0099] Among the 65,206 DEGs in Step (1), 30 DEGs were randomly screened, and according to the gene sequences, specific primers were designed using Primer5.0 software, and two internal reference genes were used. The primer sequences of the genes are shown in Table 2. The specific implementation method is the same as that in Step 3(4), and the results are as Figure 7 , the RNA-Seq data was significantly correlated with the qRT-PCR results (R 2 >0.7, p<0.05), indicating the accuracy of the sequencing data.

[0100] 6. Alternative splicing (AS) analysis and reverse transcription polymerase chain reaction (RT-PCR)

[0101] (1) Prediction of AS and identification of AS candidate events in the terpenoid biosynthesis pathway

[0102] Predict AS candidate events based on the transcripts after redundancy removal of the three-generation non-parametric transcriptome. Use the BLAST software to perform pairwise alignment of all sequences. Sequences meeting the following criteria are considered candidate AS events: 1) The lengths of both sequences are greater than 1000 bp, and there are two HSPs (High-scoring Segment Pairs) in the alignment; 2) The variable splicing gap is greater than 100 bp, and it is at least 100 bp away from the 3' / 5' end; 3) An overlap of 5 bp of all variable transcripts is allowed. In WC and ZY, 6647 and 3795 AS events were found respectively. Among the 526 DEGs in step 5(3), 16 genes had AS events, including structural genes for terpenoid biosynthesis [1-deoxy-D-xylulose 5-phosphate synthase (DXS), 1-deoxy-D-xylulose 5-phosphate reductoisomerase (DXR), 4-hydroxy-3-methylbut-2-enyl diphosphate synthase (HDS), farnesyl diphosphate synthase (FPPS), and germacrene synthase (GDS)] and TFs (bHLH, WRKY, and AP2 / ERF).

[0103] (2) AS candidate events in the terpenoid biosynthesis pathway and verification by RT-PCR

[0104] Among the 16 AS event sequences predicted in step (1), 6 AS event sequences were randomly selected, and specific primers were designed using Primer 5.0 software according to the open reading frame (ORF) of the sequences (Table 2). Using the cDNA reverse transcribed from the total RNA of the mixed samples of WC and ZY as a template, RT-PCR reactions were carried out for verification. The results are as Figure 8 shown, and the AS events truly exist.

[0105] 7. Weighted gene co-expression network analysis (WGCNA)

[0106] According to the expression level of each gene, using the Biomarker Cloud platform (https: / / international.biocloud.net / zh / software / tools / detail / small / 8a8300b253cf73e70153d16368250f32), perform WGCNA analysis on the 560 genes in step 5(3) to study the relationship between terpenoids and tissue samples. The results are as Figure 9 shown, Figure 9 In the dendrogram of a, 6 different modules were identified in the main branches, which are represented by different colors. Pearson correlation coefficient analysis showed that the 6 modules were associated with different tissues, and 4 co-expression modules were highly correlated with a single specific sample (r≥0.8) ( Figure 9b). The brown module was specifically related to WC-O (r = 0.98), and 79 genes in this module were highly aggregated in WC-O. The yellow module was related to WC-I (r = 0.94), and 48 genes in this module were highly aggregated in WC-I. The blue module was related to ZY-O (r = 0.88), and 96 genes in this module were aggregated in ZY-O. Finally, the green module was related to ZY-I (r = 0.87), and 52 genes in this module were highly accumulated in ZY-I.

[0107] 8. lncRNAs Targeting Candidate Genes for Terpenoid Synthesis

[0108] The target genes of lncRNAs were predicted by analyzing the correlation of lncRNA and mRNA expression levels among samples, and 4719 lncRNA sequences predicted in step 3(3) were used for target gene prediction. Combining the genes highly aggregated in the transcriptional modules in step 7 and the results of target gene prediction analysis, the lncRNA-regulated mRNA network diagram was further constructed using the Biomarker Cloud Platform (https: / / international.biocloud.net / zh / software / tools / detail / small / 8a8300e260593ad30160596333f10000). The results are shown in Figure 10. A total of 219 lncRNAs (FPKM≥10) targeted 39 candidate genes, including 18 structural genes and 21 transcription factors, indicating that these lncRNAs play important roles in terpenoid synthesis.

[0109] 9. Identification of Transcriptional Regulatory Modules and Key Genes Involved in Terpenoid Synthesis

[0110] (1) Visual Network Analysis of Transcriptional Regulatory Modules Involved in Terpenoid Synthesis

[0111] The genes highly aggregated in the transcriptional modules in step 7 were subjected to visual network analysis on the OE Biotech Cloud Platform (https: / / www.OmicShare.com / Tools / Home / Soft / cell cape) to find highly connected hub genes. The results are as Figure 11 shown. The heatmap shows that the genes in the green module have high expression levels in ZY-I samples ( Figure 11 a), and the network diagram shows that the highly connected hub genes include structural genes (AACT1, AACT2, and HMGR4) and TFs (bZIP1, bZIP3, bZIP5, AP2 / ERF5, and AP2 / ERF6). Figure 11a). Genes in the yellow module had relatively high expression levels in the WC-I sample. The hub genes with high connectivity were AACT8, HMGS4, HMGS6, AP2 / ERF35, AP2 / ERF36, and AP2 / ERF38( Figure 11 b). Genes in the blue module had relatively high expression levels in the ZY-O sample. There were a total of 14 hub genes with high connectivity, including 8 structural genes (HDR1, HDR2, MDS1, HDR3, HDR4, DXR1, DXR2, and DXR6) and 6 TFs (AP2 / ERF28, AP2 / ERF31, AP2 / ERF32, AP2 / ERF34, bHLH10, and bHLH11)( Figure 11 c). Genes in the brown module had relatively high expression levels in the ZY-O sample. The hub genes with high connectivity included HDS6, DXR11, DXR12, DXR13, and DXR15, as well as TFs: bZIP18, bZIP19, bHLH28, bHLH29, WRKY35, WRKY37, and AP2 / ERF63( Figure 11 d). A total of 40 highly connected hubs were screened from the above 4 transcription modules for subsequent analysis.

[0112] (2) Correlation analysis of the content of significantly different terpenoids between hub genes and samples

[0113] Perform a correlation analysis on the 40 hub genes in step (1) and 18 significantly different terpenoids. According to the gene expression levels and the relative contents of terpenoids, perform a correlation analysis on the OmicShare cloud platform (https: / / www.omicshare.com / tools / Home / Soft / ica2), and use TBtools software to draw the figure as shown. The results are as follows Figure 12As shown, there is a significant correlation between the expression levels of AACT1, AACT2, HMGR4, AP2 / ERF5, AP2 / ERF6, bZIP1, bZIP3 and bZIP5 genes and the relative contents of various terpenoids, including elemene (r > 0.93, P < 0.01), isoeugenol (r > 0.91, P < 0.01), naphthalene, decahydro-4a-methyl-1-methylene-7-(1-methylethenyl)-[4aR-(4a,αα,2,6-octadien-1-ol, 3,7-dimethyl-, acetate (r > 0.78, P < 0.07), α-pinene (r > 0.76, P < 0.08), 6-octen-1-ol, 3,7-dimethyl-, (R)-(r > 0.75, P < 0.09), linalool (r > 0.72, P < 0.1); the expression level of HDR2 is significantly correlated with the relative contents of (R)-3,7-dimethyl-6-octenal (r > 0.87, P < 0.03) and α,α-dimethyl-4-methylenecyclohexanol (r > 0.81, P < 0.05). The expression levels of HDR3 (r > 0.90, P < 0.02), HDR4 (r > 0.81, P < 0.05) and bHLH10 (r > 0.82, P < 0.05) are significantly correlated with the relative content of (R)-3,7-dimethyl-6-octenal. Finally, 12 key candidate genes (AACT1, AACT2, HMGS4, AP2 / ERF5, AP2 / ERF6, bZIP1, bZIP3, bZIP5, HDR2, HDR3, HDR4, bHLH10) were screened out to participate in the regulation of the biological accumulation of terpenoids in Zanthoxylum bungeanum leaves.

[0114] As described above, it is only a preferred specific embodiment of the present invention, but the protection scope of the present invention is not limited thereto. Any changes or substitutions that can be easily conceived by those skilled in the art within the technical scope disclosed by the present invention should be covered by the protection scope of the present invention.

Claims

1. A method for mining key genes for the synthesis of terpenoid compounds in Zanthoxylum bungeanum leaves based on third-generation full-length and second-generation transcriptome sequencing, characterized in that, The method includes: Step 1: Determine the components and contents of terpenoids in three different development stages of Zanthoxylum armatum and Rongchang thornless Zanthoxylum bungeanum, namely young leaves, mature leaves and old leaves, and compare the differences. Step 2: Extract the total RNA from the young leaves, mature leaves and old leaves of Zanthoxylum armatum and Rongchang thornless Zanthoxylum bungeanum respectively, and extract the total RNA from the mixed samples of young leaves, mature leaves and old leaves of Zanthoxylum armatum and Rongchang thornless Zanthoxylum bungeanum. Sequence them respectively and construct full-length transcriptomes and transcriptome datasets to obtain raw data. Step 3: After correcting the raw data obtained in Step 2, obtain high-quality sequences. After filtering the high-quality sequences, obtain clean data. Assemble the clean data to obtain the full-length transcriptome databases and transcriptome databases of Zanthoxylum armatum and Rongchang thornless Zanthoxylum bungeanum. Step 4: Predict alternative splicing candidate events and simple sequence repeats in the full-length transcriptome databases and transcriptome databases obtained in Step 3, and verify the authenticity of the alternative splicing candidate events. Step 5: Predict the coding region sequences and lncRNAs of the transcript sequences in the full-length transcriptome databases and transcriptome databases obtained in Step 3. Step 6: Perform functional annotation and prediction of transcription factors on the genes in the transcriptome databases in Step 3, identify differentially expressed genes between samples, and screen for enzyme genes and transcription factors related to terpenoid synthesis based on functional annotation and differential gene information. Step 7: Predict the target genes of the lncRNAs in Step 5, and perform regulatory network analysis on the lncRNAs of the genes screened in Step 6 according to the lncRNA target gene prediction information. Step 8: Perform weighted gene co-expression network analysis on the genes in Step 6, and construct a gene regulatory network diagram through visualization to screen for hub genes with high connectivity to terpenoid metabolism. Step 9: Perform correlation analysis on the significantly different terpenoids in Step 1 and the hub genes in Step 8, and select the genes correlated with the terpenoid content as the key genes for terpenoid synthesis.

2. The method according to claim 1, wherein In Step 3, convert the original sequences into RoI sequences to obtain full-length sequences and non-full-length sequences. Use the ICE algorithm to perform clustering analysis on the RoI sequences from the same transcript to obtain consensus sequences, and use the non-full-length sequences to correct the obtained consensus sequences to obtain high-quality sequences.

3. The method according to claim 1, wherein In Step 3, the method for filtering high-quality sequences is: remove the rRNA repeats and data containing adapters, and remove the low-quality data with a quality value Q ≤ 20 and more than 50% bases to obtain clean data.

4. The method according to claim 1, wherein In Step 4, use the BLAST software to perform pairwise alignment on all sequences in the full-length transcriptome databases and transcriptome databases. The alignment results meet the following conditions: 1) The lengths of both sequences are greater than 1000 bp, and there are two HSPs in the alignment; 2) The alternative splicing Gap is greater than 100 bp and at least 100 bp away from the 3' / 5' ends; 3) Sequences allowing a 5-bp overlap of all alternative transcripts are considered candidate alternative splicing events.

5. The method according to claim 1, characterized in that In step 5, four methods, namely CPC analysis, CNCI analysis, pfam protein domain analysis and CPAT analysis, are used to predict lncRNAs, and transcripts with a length exceeding 200 nt and more than two exons are selected as lncRNA candidates.

6. The method according to claim 1, wherein In step 6, the genes in the transcriptome database of Zanthoxylum bungeanum were subjected to sequence alignment and functional annotation using the BLASTN software in the Nr, Swiss-Prot, Pfam, KEGG, GO, and COG databases, where E-value < 10 -5 ; then the iTAK software was used to predict and classify the transcription factors in the database; then DESeq2 was used to analyze the differentially expressed genes in the leaves of Zanthoxylum armatum and Zanthoxylum bungeanum var. inermis at three different developmental stages, and the Benjamini-Hochberg method was used to correct the significant p-values obtained from the original hypothesis test, with a false discovery rate of less than 0.01 and a fold change ≥ 2 as the screening criteria for differentially expressed genes; according to the annotation information and differentially expressed gene information, the structural genes and transcription factors related to terpenoid synthesis were screened.

7. The method according to claim 1, characterized in that, In step 7, the target genes of lncRNAs are predicted by the method of analyzing the correlation of the expression levels of lncRNAs and mRNAs among samples.

8. The method according to claim 1, wherein The method further includes: verifying the full-length transcriptome and the results of transcriptome sequencing, including: randomly selecting the lncRNAs predicted in step 5 and the differentially expressed genes screened in step 6, using qRT-PCR technology to measure the expression levels of the lncRNAs and the differentially expressed genes, and calculating their correlation with the FPKM value data of the sequencing. If the expression levels of the lncRNAs, the expression levels of the differentially expressed genes and the FPKM value data of the transcriptome sequencing are correlated, the sequencing results are accurate and reliable.

9. Use of the method according to any one of claims 1-8 in screening genes related to the biosynthesis of terpenoid compounds in Zanthoxylum bungeanum leaves.

10. A gene related to the biosynthesis of terpenoid compounds in Zanthoxylum bungeanum leaves mined by the method according to any one of claims 1-8, characterized in that, The genes include: AACT1, AACT2, HMGR4, AP2 / ERF5, AP2 / ERF6, bZIP1, bZIP3, bZIP5, HDR2, HDR3, HDR4, and bHLH10.

Images