Yingting Zhang
Jiebing Cui
Hailiang Hu
Jin Xu- Key Laboratory of Forest Genetics and Biotechnology of Ministry of Education, Co-innovation Center for Sustainable Forestry in Southern China, College of Forestry, Nanjing Forestry University, Nanjing, China
Chinese cedar (Cryptomeria fortunei) is a tree species with important ornamental, medicinal, and economic value. Terpenoids extracted from the essential oil of C. fortunei needles have been considered valuable ingredients in the pharmaceutical and cosmetic industries. However, the possible gene regulation mechanisms that limit terpenoid biosynthesis in this genus are poorly understood. Here, we adopted integrated metabolome analysis, transcriptome, small-RNA (sRNA), and degradome sequencing to analyze the differences in terpenoid regulatory mechanisms in two different overwintering C. fortunei phenotypes (wild-type and an evergreen mutant). A total of 1447/6219 differentially synthesized metabolites (DSMs)/unigenes (DEGs) were detected through metabolome/transcriptome analyses, and these DSMs/DEGs were significantly enriched in flavonoid and diterpenoid biosynthesis pathways. In C. fortunei needles, 587 microRNAs (miRNAs), including 67 differentially expressed miRNAs (DERs), were detected. Among them, 8346 targets of 571 miRNAs were predicted using degradome data, and a 72-miRNA-target regulatory network involved in the metabolism of terpenoids and polyketides was constructed. Forty-one targets were further confirmed to be involved in terpenoid backbone and diterpenoid biosynthesis, and target analyses revealed that two miRNAs (i.e., aly-miR168a-5p and aof-miR396a) may be related to the different phenotypes and to differential regulation of diterpenoid biosynthesis. Overall, these results reveal that C. fortunei plants with the evergreen mutation maintain high terpenoid levels in winter through miRNA-target regulation, which provides a valuable resource for essential oil-related bioengineering research.
Introduction
Cryptomeria fortunei (Chinese cedar), an endemic to China, is an evergreen tree species of the family Taxodiaceae. It has become a main garden and timber species in southern China due to its beautiful tree shape, rapid growth and great adaptability, and it is important in ecological environment construction and wood production. Furthermore, it has high value as a traditional medicinal material. Research on the chemical components of C. fortunei essential oil is an active field of medicinal research (Lin et al., 2005; Xie et al., 2011). Studies have shown that compounds extracted from C. fortunei essential oil have excellent detoxification, antiulcer, antioxidant, antifungal, and insecticidal properties, among others (Matsunaga et al., 2000; Wang et al., 2006; Cheng et al., 2012; Kim et al., 2013); therefore, they have long been used as raw materials in the plant protection, pharmaceutical, and cosmetic industries. Several classes of chemical components, including (oxygenated) mono- (oxygenated), sesqui-, and (oxygenated) diterpenoids, have been shown to be present in the essential oil of C. fortunei (Xie et al., 2011). In general, the functional properties of C. fortunei essential oils benefit from the presence of a very large number of active ingredients, especially terpenes, terpene alcohols, and flavonoids (Cheng et al., 2009). According to previous reports, within Cryptomeria, there are certain differences in the chemical composition of its essential oils among different species, trees of different ages and different places of production, and different tree organs or tissues (Lin et al., 2005; Cheng et al., 2009; Xie et al., 2011). Therefore, it is necessary to understand the mechanism underlying terpenoid synthesis in C. fortunei and to enhance its production overall.
Among the plant specialized metabolites, terpenoids are the most diverse and abundant; to date, more than 75,000 terpenoids have been discovered (Fujihashi et al., 2018). Although the structures of plant terpenoids vary greatly, the synthetic precursors of all plant terpenoids are the five-carbon (C5) compound isopentenyl diphosphate (IPP) and its isomer dimethylallyl diphosphate (DMAPP), which undergo rearrangement and cyclization reactions resulting in the synthesis of various terpenoid skeletons. Terpenoids can be divided into hemi- (C5), mono- (C10), sesqui- (C15), di- (C20), tri- (C30), tetra- (C40), and polyterpenoids (C > 40) (Ashour et al., 2010) based on the number of C5 units contained in the molecule. Studies have shown that terpenoid synthesis involves two biosynthetic pathways: the mevalonic acid (MVA) pathway, through which sesqui- and triterpenoids are produced, and the 1-deoxy-D-xylosyl-5-phosphate (DOXP)/2-C-methyl-D-erythritol-4-phosphate (MEP) pathway, through which monoterpenoids, diterpenoids, tetraterpenoids, zeatin, and other terpenoid quinones are synthesized (Vranová et al., 2013). In the life activities of plants, terpenoids are important in signal transduction, resistance to (a)biotic stress, and biological interactions (Izumi and Hirata, 1997; Grassmann et al., 2002).
MicroRNAs (miRNAs) are endogenous non-coding RNAs with a length of 21–25 nt and play a vital regulatory role in plant growth and development and adaptation to environmental changes. Specifically, they are important in flower transformation and development, senescence, and (a)biotic stress responses (Lu and Huang, 2008; Song et al., 2013; Xu et al., 2013). In recent years, as high-throughput technology has developed, small-RNA (sRNA) sequencing has been widely used in the identification of plant miRNAs, including those of the camphor tree (Cinnamomum camphora) (Chen et al., 2019), Masson pine (Pinus massoniana) (Ye et al., 2020), and ginkgo (Ginkgo biloba) (Jia et al., 2020). In addition, miRNAs play a very important role in regulating gene expression through their ability to shear target messenger RNAs (mRNAs) and to inhibit the translation of targeted mRNAs (Bartel, 2004). Computer algorithms such as TargetScan (Grimson et al., 2007), MiRanda (Doron et al., 2008), and PicTar (Krek et al., 2005) are the methods most commonly used to identify miRNA targets. However, the miRNA targets that are currently known are very limited and do not provide a sufficient basis for prediction, and the identification of predicted candidate genes is relatively cumbersome; therefore, it is difficult to achieve high-throughput and large-scale identification. With the emergence of degradome sequencing (degradome-seq), high-throughput identification of plant miRNA target genes has become possible; therefore, this method has been widely used to predict plant miRNA target genes. For example, C. camphora (Chen et al., 2019), P. massoniana (Ye et al., 2020), and kidney bean (Phaseolus vulgaris) (Parreira et al., 2021) have all been shown to harbor a variety of miRNA targets. In recent years, the role of miRNAs and their targets in regulating the biosynthesis and accumulation of specialized metabolites such as flavonoids, terpenoids, and alkaloids in various plants has been reported. For example, in Arabidopsis thaliana, miR156 targets squamosa promoter binding protein-like (SPL9) and negatively regulates (E)-β-caryophyllene biosynthesis (Yu et al., 2015), and in bitter gourd (Picrorhiza kurroa), miR4995 targets 3-deoxy-7-phosphoheptulonate synthase (DHS) and affects the production of picroside-I by regulating terpenoid biosynthesis (Vashisht et al., 2015). However, there are few reports concerning the regulation of terpene synthesis by miRNAs and their transcription factors (TFs)/target genes in C. fortunei; therefore, further research is needed to better understand the regulatory role of miRNAs in the production of terpenoids.
In this study, using two phenotypes of C. fortunei needles [normal (wild-type with yellowish-brown needles in winter, YWt) and an evergreen mutant (GM)] in winter as materials, transcriptome, and metabolome analyses were used to identify differentially synthesized metabolites (DSMs) and unigenes (DEGs), respectively. These DSMs/DEGs were significantly enriched in terpenoid-related pathways. Then, through analysis of the sRNAs and the degradome, miRNAs present in C. fortunei needles and their targets were further identified, and the role of the miRNA-target module in terpenoid biosynthesis was analyzed. Our results reveal the presence of differential molecular mechanisms of terpenoid biosynthesis in different C. fortunei phenotypes in winter; this finding helps direct the selection of C. fortunei for essential oil extraction and provides a theoretical basis for its further genetic improvement to increase its terpenoid content in winter.
Materials and Methods
Plant Material and Extraction of Total RNA
Two representative C. fortunei clones with different overwintering phenotypes [#X1 (GM) and #3 (YWt)] cultivated in the Garden Experimental Teaching Center of Nanjing Forestry University (32°04’41″N, 118°48’43″E, Nanjing, Jiangsu, China) were selected as experimental materials. For each clone, needle samples were collected from among the upper and middle needles of the secondary lateral branches of 5-year-old trees. According to our previous research and related references, February is the period in which Cryptomeria shows the greatest phenotypic differences (Han et al., 2003). Therefore, the needle samples were collected on February 15, 2020, and all samples were immediately stored at −80°C. Three biological replicates were used as materials for transcriptomic analysis and sRNA sequencing, and six biological replicates were used for metabolome analysis (quantitative analysis of all metabolites).
Total RNA of C. fortunei was extracted using the RNeasy Plant Mini Kit (Qiagen, Hilden, North Rhine-Westphalia, Germany) following the manufacturer’s instructions. RNA integrity was evaluated with an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, United States), and only RNA with an RNA integrity number (RIN) value >7.0 was used for high-quality sequencing.
De novo Transcriptome Sequencing
Transcriptome sequencing was conducted by OE Biotech Co., Ltd. (Shanghai, China). To obtain clean reads (CRs), adaptor sequences, low-quality reads, and poly-N sequences were removed from the raw data using Trimmomatic version 0.36 (Bolger et al., 2014) with the following parameters: LEADING: 3, TRAILING: 3, SLIDINGWINDOW: 4:15, and MINLEN: 50. The CRs were then assembled de novo into transcripts using Trinity version 2.4.0 (Grabherr et al., 2013) in paired-end mode, with RF and SS_lib_type as the parameters. The longest transcript was selected as a unigene according to the sequence length and similarity and used in the subsequent functional annotation and expression calculation.
The unigenes were annotated through alignment with the National Center for Biotechnology Information (NCBI)’s non-redundant (Nr), eukaryotic orthologous groups of proteins (KOG), and SwissProt databases using blastx (Altschul et al., 1990) with a threshold E-value of 10–5. These resulting SwissProt IDs were mapped to the Gene Ontology (GO) database, and GO annotations of the unigenes were obtained. Finally, the unigenes were compared to the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (Kanehisa et al., 2008) to obtain pathway information.
The fragments per kilobase of exon model per million mapped fragments (FPKM) values (Trapnell et al., 2010) of each unigene were calculated using bowtie2 version 2.3.3.1 (Langmead and Salzberg, 2012). The functions estimateSizeFactors and nbinomTest in DESeq version 1.26.0 (Anders and Huber, 2016) were used to identify DEGs, and p-values < 0.05 and fold change >2 or <0.5 were set as their thresholds. KEGG and GO analyses of these DEGs were performed using R according to the hypergeometric distribution to predict the biological functions and pathways mainly affected by the DEGs. Hierarchical cluster analysis of these DEGs was performed to explore the expression patterns of the related transcripts.1
sRNA Library Construction and Sequencing
High-quality CRs were obtained after removing adaptor sequences, sequences <15 or >41 bp in length, reads containing N bases, and low-quality sequences from the raw reads (RRs). The length distribution of CRs was determined to preliminarily evaluate the sRNA distribution of the sample. The CRs were compared to the transcriptome to determine the comparison rate. Non-coding RNAs were annotated as ribosomal RNAs (rRNAs), transfer RNAs (tRNAs), small nuclear RNAs (snRNAs), and small nucleolar RNAs (snoRNAs). The RNAs were aligned and then subjected to a blastn search (Altschul et al., 1990) against the Rfam version 10.02 and GenBank databases.3
Known miRNAs were identified through alignment with the miRBase database version 21.0.4 Unannotated sRNAs were analyzed using miRDeep2 version 2.0.0.8 (Friedlnder et al., 2012) to predict novel miRNAs. The miRNA expression levels were quantified as transcripts per kilobase million (TPM) values. DESeq version 1.26.0 (Anders and Huber, 2016) was used to standardize the expression of miRNA in each sample and to calculate the fold change, and the negative binomial distribution test method was used to test the significance of the difference in miRNA expression data. Differentially expressed miRNAs (DERs) were identified with the threshold of p-value < 0.05 and | log2FoldChange| > 1.
Degradome-seq and Target Gene Identification
To construct degradation libraries, two 30-μg aliquots of total RNA extracted from the GM and YWt phenotypes, respectively, of C. fortunei were purified using the Total RNA Purification Kit (LC Science, Houston, TX, United States). The mRNAs were connected to 3′ and 5′ adaptors, reverse-transcribed using a mixture of biotinylated random primers and mRNAs, and amplified by polymerase chain reaction (PCR) to complete library construction. These constructed libraries were sequenced using an Illumina HiSeq 2500 sequencer (Illumina, San Diego, CA, United States). The single-end read length was 1 × 50 bp, and the sequencing data contained more than 5 million reads.
The raw sequencing data were filtered using Illumina Pipeline version 1.5. The miRNA targets and their cleavage sites were predicted using the CleaveLand program version 4.0 (Addo-Quaye et al., 2009) and mapped to the C. fortunei transcriptome data. Based on the abundance and characteristics of superaccumulated single genes in the C. fortunei transcriptome data, a degradome density file was established. According to the abundance of cleavage site reads (CSRs), miRNA targets can be divided into five categories as the reliability decreases: category 0 (1 maximum CSR), category 1 (>1 maximum CSR), category 2 (CSR abundance > median transcript abundance), category 3 (CSR abundance ≤ recorded median transcript abundance), and category 4 (1 RR containing the cleavage site). These target mRNA sequences that matched an sRNA sequence from the sequenced species were predicted by TargetFinder software (Lavorgna et al., 1999) with default parameters. To further study miRNAs and their targets, KEGG and GO analyses were performed on all targets.5 Subsequently, these interaction relationships were applied to construct the entire miRNA–gene interaction network, and the network was then visualized in Cytoscape version 3.5.1.6
Metabolomics Study
An accurately weighed 80-mg needle sample was transferred to a 1.5-mL microcentrifuge tube and extracted using 20 μL of 1.5 mM 2-chloro-L-phenylalanine and 1 mL of 70% methanol. The precooled sample was ground at 60 Hz for 2 min, sonicated at ambient temperature for 0.5 h, allowed to stand at −20°C for 20 min, and centrifuged at 13,000 rpm at 4°C for 10 min. Then, 300 μL of the supernatant was evaporated to dryness and reconstituted with 400 μL of 20% methanol. The sample was vortexed well, precooled at 4°C for 2 min and centrifuged for 10 min. Finally, 150 μL of the supernatant was aspirated using a syringe, filtered, and transferred to a liquid chromatography (LC) injection bottle, and stored at −80°C for liquid chromatography–mass spectrometry (LC–MS) determination.
An LC–MS system consisting of an AB ExionLC ultra-high-performance LC instrument in tandem with an AB TripleTOF 6600 plus high-resolution mass spectrometer was used to analyze the metabolic spectrum in positive and negative ESI modes. An ACQUITY UPLC HSS T3 (100 mm × 2.1 mm, 1.8 μm) column was used, and water containing 0.1% formic acid and acetonitrile (v/v) was used as a binary mobile phase for gradient separation. Sample MS signal acquisition was performed in two ESI modes, and quality control (QC) samples were introduced throughout the entire analysis process to evaluate repeatability. The raw data collected by UNIFI 1.8.1 were processed by the metabolomics processing program Progenesis QI version 2.3 (Non-linear Dynamics, Newcastle, United Kingdom), and identification was then performed based on http://www.hmdb.ca/, http://www.lipidmaps.org/ and a self-built database. For the extracted data, the resulting matrix was further reduced through removing peaks for which values were missing in >50% of the samples. Finally, both ion mode data were merged into a table that was imported into the R ropls package.
After mean centring and Pareto variance scaling, principal component analysis (PCA) and (orthogonal) partial least squares-discriminant analysis [(O)PLS-DA] were performed to visualize the metabolic differences between the experimental groups. In the model score graph, the T2 area of Hotelling shown in the oval shape defines the 95% confidence interval of the model variables. The variable importance in the projection (VIP) ranks the overall contribution of each variable to the OPLS-DA model; VIP >1 indicates a potential biomarker. The default sevenfold cross-validation and 200-response permutation testing methods were used to examine the quality of the model to prevent overfitting. In addition, DSMs were selected based on a combination of the statistical significance threshold of the VIP value obtained from the OPLS-DA model and the p-value of the normalized peak area of the two-tailed Student’s t-test. Metabolites with VIP values >1.0 and p-values < 0.05 were considered DSMs. In addition, DSMs were mapped to the KEGG database to obtain their common pathway information. The KEGG IDs of the DSMs were used to perform enrichment analysis to obtain the enrichment results of the metabolic pathways. A hypergeometric test was then applied to find pathway entries that were significantly enriched in DSMs compared with the entire background metabolites. The bubble chart of KEGG pathway enrichment analysis was visualized based on https://cloud.oebiotech.cn/task/.
Validation of DERs
Nine DERs were selected for determination of expression levels using quantitative real-time PCR (qRT-PCR). Based on these obtained sequences, gene-specific primers were designed (Supplementary Table 1).
MicroRNAs were extracted from 0.1 g needle samples using the miRNA Isolation Kit (Tiangen Biotech Co., Nanjing, Jiangsu, China) following the manufacturer’s instructions. Qualified samples that had been evaluated for integrity, purity, and concentration were used to synthesize first-strand cDNAs using the miRcute Plus miRNA First-Strand cDNA Synthesis Kit (Tiangen Biotech Co.). qRT-PCR was performed on the same instrument as that used by Zhang et al. (2021) with the following amplification procedure: 95°C predenaturation for 15 min followed by 40 cycles of template denaturation and annealing; a dissociation curve was generated at 60–95°C. The 2–Δ Δ Ct method was used to calculate the relative expression levels with U6 and cln-miR6725 as references (Supplementary Table 1). All experiments were performed in triplicate, and each sample was examined in three technical replicates.
Results
Metabolomic Profiling Analysis
To investigate the metabolic differences between the C. fortunei phenotypes, we analyzed the non-targeted metabonome of C. fortunei needles using LC–MS. The total ion current of 12 samples showed large peak capacity and excellent retention time reproducibility (Supplementary Figure 1), indicating that the difference between biological replicates and the error caused by the instrument were both small. In these needles, a total of 9251 compounds with definite formulas were identified (MTBLS3183), and the parameters for metabolite annotation were described according to the guideline by Fernie et al. (2011), including 5815 lipid-like molecules, accounting for 62.86% of all metabolites. Among them, 173 metabolites were identified as terpenoids, including 87 sesqui-, 35 di-, 32 mono-, 18 tri-, and 1 polyterpenoid(s) (Supplementary Table 2). PCA and (O)PLS-DA accurately divided all the samples into two different clusters, reflecting the obvious difference between green and yellow needles (Figures 1A–C). The accuracy of the OPLS-DA model was also tested using sevenfold cross-validation (Figure 1D).
Figure 1. Multivariate statistical analysis of the metabolome and differentially synthesized metabolites (DSMs). Plots of principal component analysis (PCA) (A), partial least squares discriminant analysis (PLS-DA) (B), and orthogonal PLS-DA (OPLS-DA) (C). (A,B) The x- and y-axes represent the first (PC1) and the second (PC2) principal components, respectively; (C) the x- and y-axes represent predictive principal components and orthogonal principal components, respectively. (D) The 200-response sorting tests of the OPLS-DA model. (E) Hierarchical clustering heatmap of DSMs. Each column and row represent a sample and a differentially expressed unigene (DEG), respectively, and the colors indicate the expression levels of DSMs. GM, evergreen mutant; YWt, yellowish-brown needles in winter; GM_1, GM_2, GM_3, GM_4, GM_5, and GM_6 represent six repetitions.
A total of 1447 DSMs were detected by combining single- and multidimensional analysis; these DSMs included 31 terpenoids (10 di-, 6 mono-, and 15 sesquiterpenoids) (Supplementary Table 3), and the synthesis of most metabolites was enhanced in GM (Figure 1E and Supplementary Table 4). In addition, 111 DSMs were enriched in 56 metabolic pathways, of which the top three were “arachidonic acid metabolism” (ath00590), “flavonoid biosynthesis” (ko00941), and “diterpenoid biosynthesis” (ko00904) (Figure 2). Therefore, we preliminarily speculated that C. fortunei that display different phenotypes may differentially accumulate these metabolites.
Figure 2. Kyoto Encyclopedia of Genes and Genomes analysis of differentially synthesized metabolites (DSMs). The x- and y-axes represent the enrichment factor and the pathway term, respectively. The colors and sizes of the dots represent the significance (p) and the number of metabolites, respectively.
Assembly, Annotation, and Enrichment Analysis of the C. fortunei Transcriptome
To identify genes expressed in C. fortunei needles, six transcriptome libraries (GM_1, GM_2, GM_3, YWt_1, YWt_2, and YWt_3) were generated using total RNA extracted from three biological replicates of each phenotype. A total of 304.79 M RRs were obtained, and the sequencing results were submitted to the NCBI database under registration number PRJNA697258 (SAMN17672908–17672913). The low-quality sequences were filtered out, and in total, 302.75 M CRs were obtained; the Q30 (%) was >95.57%, and the gas chromatography (GC) content was 44.19–44.33% (Supplementary Table 5). Then, 55,839 single unigenes (submission number: SUB10083705, Supplementary Tables 6, 7) were assembled using Trinity software. The N50 value and the average length were 1219.86 and 1990 bp, respectively. The unigenes expressed in the GM and YWt samples clearly aggregated into different clusters, indicating that the biological reproducibility was good and that the samples could be used for subsequent analysis (Supplementary Figure 2). Finally, through BLAST comparison with seven public databases, a total of 27,769 single unigenes were annotated, of which 27,144 unigenes that significantly matched the Nr database accounted for the highest percentage of the total (48.61%) (Supplementary Tables 6, 8).
In total, 6219 DEGs were detected when the GM and YWt libraries were compared; these included 3096 upregulated and 3123 downregulated DEGs (Supplementary Figure 3). Classification of these DEGs by GO function yielded 1950 DEGs, including 913 upregulated and 1037 downregulated DEGs. The enrichment histogram showed that the most enriched category was “biological process” (BP), followed by “cellular component” (CC) and “molecular function” (MF) (Figure 3A). In BP, the most enriched subcategory was “biological regulation,” followed by “metabolic process” and “response to stimulus”; “cell,” “cell part,” and “organelle”; and “binding,” “catalytic activity,” and “transporter activity” were the three most enriched subcategories in CC and MF, respectively. Directed acyclic graphs were used to display the results of the GO structure, and the GO terms were enriched in the following order: CC, MF, and BP (Figure 3B). A total of 669 DEGs were enriched in KEGG pathways; the most significantly enriched pathway was “carbohydrate metabolism,” followed by “biosynthesis of other specialized metabolites” (Figure 3C); 805 DEGs were enriched for 112 KEGG terms, and the most highly enriched ones were “phenylpropanoid and flavonoid biosynthesis,” “glutathione metabolism,” “plant-pathogen interaction,” “cyanoamino acid metabolism,” “starch and sucrose metabolism,” “brassinosteroid biosynthesis,” “fructose and mannose metabolism,” and “diterpenoid biosynthesis” (Figure 3D).
Figure 3. Functional classification of differentially expressed unigenes (DEGs). (A) Comparison of the distribution of DEGs at the GO level. Red/green indicates GO entries enriched for up-/downregulated DEGs; the x-axis is the name of the entry, and the y-axis represents the number of unigenes in the corresponding entry and its percentage. (B) Directed acyclic graphs of three main categories. The 10 most significant nodes are represented by rectangles; the colors of the rectangles represent the significance of enrichment, and the significance increases as the color changes from yellow to red. The basic information on each node is displayed in the corresponding graph. (C) KEGG distribution map of DEGs. The x-axis shows the up- and downregulated DEGs annotated to each pathway and all these DEGs annotated to the KEGG pathway to the total number (%); the y-axis indicates the name of the pathway, and the number on the right side of the column indicates the number of DEGs annotated to each pathway. (D) KEGG enrichment analysis of DEGs. The x- and y-axes represent the enrichment score and the pathway term, respectively. The colors and sizes of the dots represent the significance (p) and the number of unigenes, respectively.
High-Throughput sRNA Sequencing
To identify the DERs related to the phenotypic differences between the two C. fortunei strains, six independent sRNA libraries were generated in February, and the sequencing results were submitted to the NCBI database; the sequence read file registration number is PRJNA720228 (SAMN18644350–18644352, 18644359–18644361). A total of 254.59 M RRs were obtained, and 21.00–26.99 M CRs were obtained in each library after removal of low-quality reads, poly-N sequences, adaptor sequences, and reads <15 or >41 bp in length (Supplementary Table 9). A total of 120,677,469 reads were mapped to the aforementioned transcriptome; approximately 0.146% (208,203) of the sRNA reads were mapped to non-coding RNAs (i.e., rRNAs, tRNAs, snRNAs, and cis-regulatory elements) in the Pfam database, and an average of 465,479.33 localization reads in each library were recognized as known miRNAs (Supplementary Table 9). The most abundant sequence among the C. fortunei sRNAs was the 21-nt sRNA tag (Supplementary Figure 4).
A total of 587 miRNAs were detected; 439 were expressed in both phenotypes, and 67 and 81 were expressed in only GM or YWt, respectively (Supplementary Figure 5A). A total of 391 miRNAs were identified as known miRNAs, 359 of which were distributed in 40 miRNA families (Supplementary Figure 5B). Sixty-seven DERs showed significant differences between the two phenotypes, including 29 downregulated and 38 upregulated miRNAs in GM (Supplementary Figure 6). Nine DERs were selected for qRT-PCR, and the expression profile of each miRNA was consistent with the sRNA-seq data (Figure 4), verifying the reliability of the sRNA results.
Figure 4. Quantitative real-time PCR (qRT-PCR) validation of selected DERs. The x- and y-axes represent the sample (GM, evergreen mutant; YWt, yellowish-brown needles in winter) and the transcripts per kilobase million (TPM) or relative expression level, and the values shown represent the mean ± SD; n = 3.
miRNA Targets Identified Through Degradome-seq
In total, 84,600,973 RRs were obtained from both GM and YWt degradomes, representing 9,473,983 and 31,301,265 unique RRs, respectively (Supplementary Table 10). After removal of the 3′ adaptor and shorter cut label reads (<15 nt), 7,223,211 (76.24% of the unique mappable reads) and 5,795,795 (18.52%) unique reads were successfully mapped to 42,337 (75.82% of all 55,839 reference transcriptomes) and 41,195 (73.77%) single genes, respectively (Supplementary Table 10).
A total of 8346 target sites were identified for 571 miRNAs (Supplementary Figure 7); these included 3910 targets of 375 known miRNAs and 4436 targets of 196 novel miRNAs. In GM, 2389 (28.62%) targets belonged to category 0, 1, or 2, and 3644 targets were classified into category 3 or 4, while in YWt, 2019 (24.19%), the targets belonged to category 0, 1, or 2, and 3068 targets were classified into the other two categories (Supplementary Table 11). It is worth noting that a certain miRNA was involved in regulating multiple targets. For example, gma-miR396h, gma-miR6300, novel116_star, novel132_mature, novel20_star, novel59_mature, novel72_star, and osa-miR414 had more than 80 targets. In contrast, single genes are usually regulated by multiple miRNAs. For example, TRINITY_DN8796_c0_g1_i1_2 may be regulated by aly-miR166b-5p, aly-miR396b-5p, gma-miR396h, and hbr-miR396a. In addition, some targets were annotated as TFs in C. fortunei, such as APETALA2 (AP2) and v-myb avian myeloblastosis viral oncogene homolog (MYB), indicating that C. fortunei miRNAs may have the ability to target TF gene families.
Enrichment Analyses of miRNA Targets
The functions of all targets were divided into three categories that included 54 GO terms by GO analysis. In BP, “cellular process,” “single-organism process,” and “metabolic process” had the highest number of target genes; in CC, the targets were mainly involved in “cell” and “cell part,” and in MF, the targets were mostly divided into “binding” and “catalytic activity” components (Figure 5).
Figure 5. Gene Ontology classification of targets in C. fortunei. The x-axis is the entry name, and the y-axis represents the number of genes in the corresponding entry.
Based on the KEGG analysis, these targets were matched to five categories: “genetic information processing,” “metabolism,” “organismal systems,” “cellular processes,” and “environmental information processing”; these categories represented 19 metabolic pathways. Among these pathways, “global and overview maps” (676), “carbohydrate metabolism” (254), and “translation” (239) were the most abundant (Supplementary Figure 8A). In addition, 72 and 48 targets were annotated to “metabolism of terpenoids and polyketides (TP)” and “biosynthesis of other specialized metabolites,” respectively (Supplementary Figure 8A). Notably, metabolism was the category with the highest number of annotated pathways. In metabolism, “genetic information processing and organismal systems,” “sesquiterpenoid and triterpenoid biosynthesis” (ko00909, 55.6%), “proteasome” (ko03050, 43.1%), and “circadian rhythm-plant” (ko04712, 48.9%) had the highest proportions of targets in each pathway (Supplementary Figure 8B).
Comprehensive Analysis of TP
The metabolome, transcriptome, and miRNA target analyses described above showed that terpenoid and polyketide biosynthesis was significantly enriched. Therefore, we next focused on analyzing these metabolic pathways.
The miRNA targets related to the metabolism of TP mainly included nine main pathways, among which “terpenoid backbone biosynthesis,” “carotenoid biosynthesis,” and “diterpenoid biosynthesis” (Table 1) were more enriched; targets involved in “terpenoid backbone biosynthesis” (23, 31.94%), “carotenoid biosynthesis” (16, 22.22%), and “diterpenoid biosynthesis” (8, 11.11%) accounted for more than 65% of the total targets (Figure 6 and Table 1). We measured the levels of metabolites related to TP and found that except for the metabolites related to carotenoid biosynthesis, which were reduced in GM, most of the metabolites were significantly enhanced in GM (Table 2). In particular, in the diterpenoid biosynthesis pathway, the levels of gibberellin A14 and gibberellin A53 aldehydes exceeded the sum of the levels of the other 20 metabolites, and the levels of these metabolites in GM were 14.14 and 8.46 times those in YWt, respectively (Table 2).
Table 1. MicroRNAs target functions involved in the metabolism of terpenoids and polyketides.
Table 2. Metabolites in terpenoids and polyketides.
Figure 6. MicroRNA regulatory network involved in the metabolism of terpenoids and polyketides. The red squares indicate metabolic pathways, the yellow triangles indicate genes, and the blue circles indicate miRNAs. The dotted lines with arrows indicate how the miRNAs regulate specific genes. The abbreviation of each gene is as follows: AACT, acetyl-CoA C-acetyltransferase; ALDH, aldehyde dehydrogenase; BAS1, PHYB activation tagged suppressor 1; CHYB, beta-ring hydroxylase; CISZOG, cis-zeatin O-glucosyltransferase; CKX, cytokinin dehydrogenase; CPD, cytochrome P450 family 90 subfamily A polypeptide 1; CPSent, ent-copalyl diphosphate synthase; crtB, 15-cis-phytoene synthase; crtZ, beta-carotene 3-hydroxylase; CYP735A, cytokinin trans-hydroxylase; CYP90D1, 3-epi-6-deoxocathasterone 23-monooxygenase; CYP97C1, carotene epsilon-monooxygenase; DET2, steroid 5-alpha-reductase; DHDDS, ditrans, polycis-polyprenyl diphosphate synthase; DMAPP-tRNA transferase, tRNA dimethylallyltransferase; DWF4, steroid 22-alpha-hydroxylase; DXS, 1-deoxy-D-xylulose-5-phosphate synthase; ECA4, Ca2+-transporting ATPase; FDFT1, farnesyl-diphosphate farnesyltransferase; FDPS, farnesyl diphosphate synthase; FNTA, protein farnesyltransferase/geranylgeranyltransferase type-1 subunit alpha; FOLK, farnesol kinase; GDS, (–)-germacrene D synthase; gcpE, (E)-4-hydroxy-3-methylbut-2-enyl-diphosphate synthase; gmer, GDP-4-keto-6-deoxy-D-mannose-3,5-epimerase-4-reductase; GGPS, geranyl diphosphate synthase; HMGCR, hydroxymethylglutaryl-CoA reductase; HMGCS, hydroxymethylglutaryl-CoA synthase; ispD, 2-C-methyl-D-erythritol 4-phosphate cytidylyltransferase; ispH, 4-hydroxy-3-methylbut-2-en-1-yl diphosphate reductase; ispS, isoprene synthase; lcyB, lycopene beta-cyclase; lcyE, lycopene epsilon-cyclase; MD, (+)-neomenthol dehydrogenase; miaA, tRNA dimethylallyltransferase; MVD, diphosphomevalonate decarboxylase; NCED, 9-cis-epoxycarotenoid dioxygenase; PDS, 15-cis-phytoene desaturase; GPS, geranyl-diphosphate specific/all-trans-non-aprenyl-diphosphate synthase; SQLE, squalene monooxygenase; VDE, violaxanthin de-epoxidase; ZDS, zeta-carotene desaturase; ZEP, zeaxanthin epoxidase.
Based on these findings, we inferred that diterpenoid biosynthesis (including biosynthesis of the terpenoid backbone) is important in the formation of evergreen mutations in GM. Forty-one miRNA targets were related to terpenoid backbone and diterpenoid biosynthesis (Figure 7A). Thirty-one target pairs, including eight known miRNA targets [atr-miR396d-geranylgeranyl diphosphate synthase (GGPS), gma-miR6300-hydroxymethylglutaryl-CoA reductase (HMGCR), 2 gra-miR482-farnesyltransferase/geranylgeranyltransferase type-1 subunit alpha (FNTA), mtr-miR5261-geranyl-diphosphate specific/all-trans-non-aprenyl-diphosphate synthase (GPS), pab-miR394a-isoprene synthase (ispS), vvi-miR159a-ditrans, polycis-polyprenyl diphosphate synthase (DHDDS), and vvi-miR172a-diphosphomevalonate decarboxylase (MVD)], were related to terpenoid backbone biosynthesis. Among them, the DEG ispS (TRINITY_DN6836_c0_g1_i1_1) cleaved by novel106_mature and pab-miR394a was identified as participating in “isoprene synthesis” (Figure 7B). Twenty-three novel miRNA targets were related to terpenoid backbone biosynthesis, which included only contains three DERs (novel59_mature, novel106_mature, and novel5_star), and they were all highly expressed in GM; however, their targets were not significantly expressed (Figure 7B). Ten target pairs were identified as being involved in diterpenoid biosynthesis, namely, 6 ent-copalyl diphosphate synthase (ent) target genes cut by aly-miR168a-5p, aof-miR396a, gma-miR395d, novel28_mature, novel29_mature, and novel40_star and 4 (13E)-labda-7,13-dien-15-ol synthase target genes cut by novel59_mature, pta-miR1314, and cme-miR828 (Figure 7C). Among them, aly-miR168a-5p, aof-miR396a, and novel40_star were significantly upregulated (by 361.635, 29.540, and 20.516%, respectively) in YWt compared with GM, while others (gma-miR395d, novel28_mature, novel29_mature, novel59_mature, pta-miR1314, and cme-miR828) were significantly upregulated in GM (28.144–52.122%). It is worth noting that the targets of aly-miR168a-5p, aof-miR396a, novel59_mature, and cme-miR828 were significantly upregulated (by 40.144–163.608%) in GM compared with YWt.
Figure 7. Analysis of miRNA targets involved in terpenoid and polyketide biosynthesis (TPB). (A) miRNA targets involved in TPB. miRNA cleavage is indicated by solid lines with red “T” arrows. miRNA targets involved in TPB (B) and diterpenoid biosynthesis (C). The pseudo-color bar at the bottom of the middle represents the fragments per kilobase of exon per million (FPKM) value, GM, evergreen mutant; YWt, yellowish-brown needles in winter. The abbreviation of each gene/metabolite is as follows: atoB, acetyl-CoA acetyltransferase; CMK, 4-(cytidine 5-diphospho)-2-C-methyl-D-erythritol kinase; CPSent, ent-copalyl diphosphate synthase; crtB, 15-cis-phytoene synthase; DHDDS, ditrans, polycis-polyprenyl diphosphate synthase; DMAPP, dimethylallyl diphosphate; DXR, 1-deoxy-D-xylulose 5-phosphate reductoisomerase; DXS, 1-deoxy-D-xylulose-5-phosphate synthase; FDPS, farnesyl diphosphate synthase; FNTA, protein farnesyltransferase/geranylgeranyltransferase type-1 subunit alpha; FOLK, farnesol kinase; FPP, farnesyl diphosphate; FPPS, farnesyl pyrophosphate synthase; gcpE, (E)-4-hydroxy-3-methylbut-2-enyl-diphosphate synthase; GGPP, geranylgeranyl diphosphate; GGPS, geranylgeranyl diphosphate synthase; GPP, geranyl pyrophosphate; HDR, 4-hydroxy-3-methylbut-2-enyl diphosphate reductase; HMBPP, (E)-4-hydroxy-3-methyl-but-2-enyl pyrophosphate; HMGCR, hydroxymethylglutaryl-CoA reductase; HMGCS, hydroxymethylglutaryl-CoA synthase; IDI, isopentenyl diphosphate isomerase; IPP, isopentenyl diphosphate; iSPD, 2-C-methyl-D-erythritol 4-phosphate cytidylyltransferase; iSPH, 4-hydroxy-3-methylbut-2-en-1-yl diphosphate reductase; iSPS, isoprene synthase; MCS, 2-C-methyl-D-erythritol 2,4-cyclodiphosphate synthase; MEP/DOXP, 2-C-methyl-D-erythritol 4-phosphate/1-deoxy-D-xylulose 5-phosphate; MVD, diphosphomevalonate decarboxylase; MVK, mevalonate kinase; Non-aprenyl-PP, non-aprenyl pyrophosphate; Phytyl-PP, phytyl pyrophosphate; PMK, phosphomevalonate kinase; GPS, geranyl-diphosphate specific/all-trans-non-aprenyl-diphosphate synthase.
To date, mainly six TF families [AP2/ethylene responsive factor (ERF), basic helix-loop-helix (bHLH), MYB, NAM, ATAF, and CUC (NAC), WRKY and basic leucine zipper (bZIP)] have been reported in plants in relation to the metabolism of terpenoids (Xu et al., 2019). In the study, only three TF families (12 AP2-like factors, 61 MYBs, and 2 WRKYs) were identified (Supplementary Table 12); among them, only five miRNAs (aly-miR172e-3p, cme-miR319c, novel137_mature, osa-miR159d, and zma-miR159e-3p) were differentially expressed. These findings suggest that these miRNA targets may play a vital role in differential diterpenoid biosynthesis and biological adaptation.
Discussion
Cryptomeria fortunei, a very important industrial essential oil tree species in China, contains an especially large amount of natural terpenes (Yao et al., 2008; Xie et al., 2011). miRNAs play an important role in the life activities of plants by regulating gene expression (Bulgakov and Avramenko, 2015). In addition, a role for miRNAs in regulating plant specialized metabolites has been reported in an increasing number of plant species (Xu et al., 2013; Vashisht et al., 2015; Gupta et al., 2017; Liu et al., 2017). However, the possible regulatory role of C. fortunei miRNA in terpenoid biosynthesis is still unclear. Therefore, we used metabonome, transcriptome, sRNA, and degradome-seq to explore the potential miRNA regulatory network of terpenoids in C. fortunei needles, with the goal of improving the level of terpenoid synthesis in winter.
Terpenoid Metabolites Were Differentially Produced in the Two Studied Phenotypes of C. fortunei
In C. fortunei needles, a total of 173 terpenoid metabolites were identified; these were primarily sesquiterpenoids but also included a variety of di-, mono-, and triterpenoids (Supplementary Table 2). Previous phytochemical studies of the branches, needles, heartwood, cones, and bark of C. fortunei have also isolated a variety of terpenoids, including mono-, sesqui-, and diterpenoids (Yao et al., 2008; Wu and Zhao, 2010), supporting the idea that terpenoids are present in large amounts in C. fortunei. However, there are obvious differences in the types and levels of compounds reported to be present. For example, Xie et al. (2011) identified only 68 terpenoids from the essential oil of C. fortunei using GC–MS. These differences in terpenoid metabolites may be related to the plant genetic background, environment, tissues/organs, sampling time, and many other factors.
We identified 31 terpenoids that were differentially biosynthesized in the two C. fortunei phenotypes analyzed in this study (Supplementary Table 3). Both metabolome and transcriptome functional analyses showed that DEGs/DSMs were enriched for diterpenoid biosynthesis (ko00904) (Figures 2, 3D). Therefore, we speculated that each of these two phenotypes of C. fortunei may have its own unique diterpenoid metabolic mechanisms. However, their specific biological regulatory mechanisms need to be further analyzed.
miRNA TF Targets May Regulate Terpenoid Synthesis
A number of studies have shown that the production of terpenoid substances is strictly regulated by microribonucleic acids through TFs. We identified three TF families (AP2, MYB, and WRKY), most of which were MYBs and 35 MYBs were cut by miR828/858 (Supplementary Table 12). Other studies have also found that most of the target genes of miR858/828 belong to MYBs, especially R2R3-MYB TFs. For example, up to 20 R2R3-MYBs are involved in cotton fiber development (Guan et al., 2014), and most of them are related to the metabolic pathways of phenylpropanes. However, few studies have examined the synthesis of terpenoids. In loblolly pine (Pinus taeda), PtMYB14 positively regulates the production of multiple sesquiterpenes (Bedon et al., 2010); in A. thaliana, AtMYB TFs increase terpenoid volatiles, not only by activating the transcription of the trehalose-6-phosphate synthase (TPS) gene but also by enhancing the metabolic flux to the isoprenoid pathway, while MsMYB negatively regulates the accumulation of monoterpenes and interferes with the synthesis of sesquiterpenes and diterpene derivatives in spearmint (Reddy et al., 2017). Therefore, whether a specific MYB regulates terpenoid synthesis and how to regulate it (e.g., negative/positive regulation, which structural gene to regulate) remains to be further studied. In addition, in Artemisia annua, the AP2 TFs AaERF1/2 participate in actively regulating expression of the sesquiterpene synthase gene (Yu et al., 2012). In Asian cotton (Gossypium arboreum), GaWRKY1 activates the transcription of the sesquiterpene synthase gene encoding cadinene synthase (CAD1), which leads to the biosynthesis of sesquiterpene gossypol (Xu et al., 2004), which is similar to our research results. Other TFs have also been found to be involved in the biosynthesis of terpenes. For example, in Arabidopsis, miR156 regulates the content of active ingredients such as sesquiterpenes by downregulating its target gene SPL9 (Yu et al., 2015). Therefore, it is still necessary to enrich the sample tissues and the cedar leaves of different growth stages in follow-up studies to further explore the mechanism by which miRNA-TF regulates terpenoid biosynthesis.
miRNA Target Genes Play an Important Role in Regulating Diterpenoid Synthesis
MicroRNAs usually partially or completely match their target genes, leading to endonuclease digestion or translation inhibition and thereby regulate the expression of the target genes and play important roles. Therefore, searching for target genes is of great value for elucidating the biological functions of miRNAs. Degradome-seq, which offers high throughput and coverage, has been suggested to significantly increase the probability of capturing target miRNA fragments and to provide a statistically significant number of all target sites for each miRNA (Pantaleo et al., 2010). It has therefore been widely used in research on plant miRNA targets (Xu et al., 2013; Chen et al., 2019; Ye et al., 2020). miRNA targets were identified by degradome-seq; through GO and KEGG analysis, these miRNA targets were found to be involved in a variety of molecular biological regulatory pathways, including terpenoid and polyketide biosynthesis (Figure 5 and Supplementary Figure 8). Further miRNA target and metabolite analysis showed that terpenoid backbone biosynthesis, diterpenoid biosynthesis and carotenoid biosynthesis play a vital role in terpenoid biosynthesis in C. fortunei (Figure 6 and Tables 1, 2). In addition, most of the orange-color-related DSMs in carotenoid biosynthesis were significantly upregulated in YWt, and most of the remaining terpenoid and polyketide biosynthesis-related metabolites were upregulated in GM (Table 2). As a common precursor, geranylgeranyl diphosphate (GGPP) has three main routes of metabolism: (1) via a reaction catalyzed by phytoene synthase (PSY), GGPP is further condensed to phytoene (C40) and participates in carotenoid biosynthesis; (2) via a reaction catalyzed by geranylgeranyl reductase (GGR), GGPP eventually forms the phytol side chain of chlorophyll; and (3) via a reaction catalyzed by diterpene synthase, various diterpene compounds, including gibberellin A and some phytoalexins, are formed. Therefore, we inferred that in C. fortunei terpenoid backbone biosynthesis, the precursor GGPP is allocated to different metabolic branches, such as carotenoid and diterpenoid biosynthesis; that is, YWt contains more carotenoids than GM, and GM contains more diterpenoids than YWt, resulting in differences in needle color between GM and YWt due to differential synthesis of diterpenoids (Figure 8).
Figure 8. Simplified model of aly-miR168a-5p- and aof-miR396a-CPSent regulating diterpenoid biosynthesis. Red/blue text indicates up-/downregulated genes; the sizes of the arrows indicate the amount of metabolic flow; and pseudo-color bar at the bottom of the middle represents the fragments per kilobase of exon per million (FPKM) value. CPSent, ent-copalyl diphosphate synthase; GGPP, geranylgeranyl diphosphate; GM, evergreen mutant; YWt, yellowish-brown needles in winter.
One miRNA can regulate hundreds of targets, and one target can also be regulated by multiple miRNAs (Licursi et al., 2019), thus constructing a complex interlocking miRNA-target regulatory network (Figure 6). Among them, we identified 41 miRNA targets synthesized via C. fortunei terpenoid backbone biosynthesis and diterpenoids. For many modules, the expression pattern of miRNA was opposite to that of its target genes, indicating that these miRNAs could negatively regulate their targets (Figures 7B,C). In terpenoid backbone biosynthesis, we identified multiple known miRNAs (atr-miR396d, gma-miR6300, gra-miR482, mtr-miR5261, pab-miR394a, vvi-miR159a, and miR172a) that regulated structural gene (such as GGPS, HMGCR, and FNTA) (Figure 7A). In Chinese toon (Toona sinensis) sprouts, tsi-miR408-3p, tsi-miR171a, tsi-miR390b-5p, and tsi-miR-1260 target these unigenes, such as HMGCR, MEP, and 1-deoxy-D-xylulose-5-phosphate synthase (DXS) (Zhao et al., 2021); only three miRNA-enzyme gene pairs related to terpenoid biosynthesis have been identified in tea trees (Camellia sinensis) (Zhao et al., 2018). There are some differences in miRNA targets, which may be related to species specificity, but these also confirm that miRNAs are indeed involved in the biosynthesis of terpenoids. We also identified multiple novel miRNA targets in this pathway (Figure 7A), which is similar to the phenomenon in T. sinensis sprouts (Zhao et al., 2021). Novel miRNAs are a type of speculative miRNA, which indicates that these genes may be the basis of specific biological functions in C. fortunei, and further experiments are needed to determine the functions of these novel miRNAs. Interestingly, the differentially expressed target gene ispS (TRINITY_DN6836_c0_g1_i1_1) cut by novel106_mature and pab-miR394a was identified as participating in “isoprene synthesis” (Figure 7). The interpretation of miRNA regulatory functions is inseparable from its target genes. In land plants, the target gene ispS is located in the chloroplast (Sasaki et al., 2005), where it catalyzes the conversion of the intermediate substrate DMAPP to isoprene (Silver and Fall, 1991). In our previous research, we found that the chloroplasts in GM were basically intact, whereas the chloroplast structure in YWt was destroyed in winter. Therefore, it was also confirmed that the synthesis of ispS was indeed inhibited in YWt. Previous studies have shown that isoprene, a volatile C5 terpenoid compound, is synthesized by plants and that the biological release of isoprene by plants has a certain protective effect (Kksal et al., 2010; Oku et al., 2015). Therefore, we inferred that GM could produce more isoprene to enhance its adaptability to low temperatures.
In the diterpenoid biosynthesis modules, that is, in aly-miR168a-5p- and aof-miR396a-CPSent and in novel59_mature- and cme-miR828-(13E)-labda-7,13-dien-15-ol synthase (CPSKSL1), the target genes were upregulated in GM (Figure 7C). CPSent is a key enzyme in the gibberellin biosynthesis pathway in plants (Otomo et al., 2004) and in the phytoalexin synthesis pathway. It catalyzes the cyclization of GGDP to ent-copalyl diphosphate (CDP) (Xu et al., 2010). Gibberellin is a very important plant hormone in higher plants. It regulates stem elongation, leaf expansion, flower differentiation and seed and fruit development, breaks dormancy, promotes germination, reduces organ shedding, and participates in many other physiological processes (Lang, 1957). Phytoalexins are small-molecular-weight compounds with antimicrobial activity that are synthesized and rapidly accumulate in the body when plants are attacked by biological or non-biological threats. The phytoalexins that have been isolated and identified to date can be roughly divided into terpenes, flavonoids, flavanones, and alkaloids. The phytoalexin components of these diterpenoids enhance the disease resistance of plants and their resistance to related exogenous damage to some extent. Therefore, we inferred that compared with YWt, the low expression of aly-miR168a-5p and aof-miR396a in GM promoted high expression of its target CPSent genes (TRINITY_DN30372_c1_g1_i15_2 and TRINITY_DN29820_c0_g1_i21_2), which in turn promoted the accumulation of gibberellins such as gibberellin A12, A14, and A3, gibberellin A53 aldehyde and gibberellin A8 catabolite (Figure 8). Therefore, GM exhibited reduced dormancy and remained green, a condition that ensured that it could produce more diterpenoids (Figure 8). The bifunctional CPSKSL1 gene product catalyzes the formation of the inner double bond isomer of cobalt cyclopentadiphosphate (CPP) and subsequently replaces diphosphate with a hydroxyl group to form labda-7,13E-dien-15-ol (Mafu et al., 2011). Although the latter compound is a known plant metabolite, its biosynthesis has not yet been studied, and further experimental verification is needed.
miRNA Targets Play a Role in Regulating Flavonoid Biosynthesis to Some Extent
Differentially synthesized metabolites, DEGs, and miRNA targets were all significantly enriched in the flavonoid metabolism pathway (Figures 2, 3D and Supplementary Figure 9). (+)-Catechin, epicatechin, luteoforol, and pentahydroxyflavanone were significantly upregulated in YWt; in that phenotype, the level of epicatechin was 11.81-fold higher than that in GM and much higher than the levels of other metabolites; afzelechin, homoeriodictyol, and pinocembrin were significantly upregulated in GM (Supplementary Table 13). These results showed that the two phenotypes of C. fortunei exhibited significant differences in flavonoid metabolism.
A total of 22 miRNA targets were related to flavonoids (Supplementary Figure 9). We further selected target pairs with different target genes, that is, ata-miR395b-3p-flavonol synthase (FLS, TRINITY_DN24453_c0_g1_i3_1, and TRINITY_DN29464_c0_g1_i4_2), novel116_star-chalcone synthase (CHS, TRINITY_DN26891_c0_g2_i2_3), novel20_star-flavanone 3-hydroxylase gene (F3H, TRINIT Y_DN19610_c1_g1_i2_1), novel55_mature-leucoanthocyanidin dioxygenase (LDOX, TRINITY_DN27808_c0_g1_i2_3), and lja-miR398-3p-leucoanthocyanidin reductase (LAR, TRINITY_DN23847_c0_g1_i4_1), and we found that only LAR was upregulated in GM, whereas the remaining target genes were upregulated in YWt. Studies have shown that miR395 acts on target genes and negatively regulates the key flavonoid gene FLS, which ultimately affects the accumulation of flavonoids (Li et al., 2018). Therefore, we believe that FLS is highly expressed in YWt and that this expression promotes flavonoid accumulation. Novel116_star, novel20_star, and novel55_mature are novel predicted miRNAs, with 204, 126, and 16 predicted targets, respectively, most of which are unknown. These genes showed no homology with genes from other species, indicating that they may support specific functions in C. fortunei. In addition, few studies have examined the role of miR398 in regulating flavonoid biosynthesis, and further experiments are needed to determine the functions of these miRNAs.
Conclusion
In summary, we analyzed the transcriptome, metabolome, sRNA, and degradome of two different overwintering C. fortunei phenotypes. The results of the metabolome and transcriptome analyses showed that DSMs and DEGs were enriched in the diterpenoid pathway. The miRNAs and their targets in C. fortunei needles were then identified, and the miRNA targets involved in terpenoid biosynthesis were analyzed. The module (aly-miR168a-5p- and aof-miR396a-CPSent) may be a potential key factor in the mechanism of differential diterpenoid synthesis in winter. In addition, we also identified three miRNA-TF families, that may be related to terpenoid biosynthesis, but further experimental verification is required. These findings provide reference information regarding diterpenoid regulatory mechanisms in C. fortunei and valuable resources for bioengineering research related to plant production of essential oils in winter.
Data Availability Statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi.nlm.nih.gov/genbank/, PRJNA697258; https://www.ncbi.nlm.nih.gov/genbank/, PRJNA720228.
Author Contributions
JX conceived and designed the experiments. JX and YZ performed the experiments and wrote the manuscript. JC, HH, JYX, and JY analyzed the data. All authors read and approved the final manuscript.
Funding
This work was supported by China’s State Forestry Administration, Forestry Public Welfare Industry Research (201304104), and the Priority Academic Program Development of Jiangsu Higher Education Institutions (PAPD).
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s Note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2021.740755/full#supplementary-material
Abbreviations
AP2, APETALA2; BP, biological process; CC, cellular component; CPSent, ent-copalyl diphosphate synthase; CRs, clean reads; CSRs, cleavage site reads; DEGs, differentially expressed unigenes; DSMs, differentially synthesized metabolites; DERs, differentially expressed miRNAs; GGPP, geranylgeranyl diphosphate; GO, Gene Ontology; ispS, isoprene synthase; KEGG, Kyoto Encyclopedia of Genes and Genomes; LC–MS, liquid chromatography-mass spectrometry; MF, molecular function; miRNAs, microRNAs; MYB, v-myb avian myeloblastosis viral oncogene homolog; (O)PLS-DA, (orthogonal) partial least squares-discriminant analysis; qRT-PCR, quantitative real-time polymerase chain reaction; RRs, raw reads; sRNA, small RNA; TFs, transcription factors; TP, terpenoids and polyketides; VIP, variable importance in the projection.
Footnotes
- ^ https://cloud.oebiotech.cn/task/
- ^ http://www.sanger.ac.uk/software/Rfam
- ^ http://www.ncbi.nlm.nih.gov/genbank/
- ^ http://www.mirbase.org/
- ^ https://www.omicshare.com/
- ^ https://cytoscape.org/
References
Addo-Quaye, C., Miller, W., and Axtell, M. J. (2009). CleaveLand: a pipeline for using degradome data to find cleaved small RNA targets. Bioinformatics 25, 130–131. doi: 10.1093/bioinformatics/btn604
Altschul, S. F., Gish, W., Miller, W., Myers, E. W., and Lipman, D. J. (1990). Basic local alignment search tool. J. Mol. Biol. 215, 403–410. doi: 10.1016/S0022-2836(05)80360-2
Anders, S., and Huber, W. (2016). Differential Expression of RNA-Seq Data at the Gene Level–the DESeq Package. Available Online at: http://mirrors.ustc.edu.cn/bioc/3.4/bioc/vignettes/DESeq/inst/doc/DESeq.pdf [accessed Jan 12, 2016]
Ashour, M., Wink, M., and Gershenzon, J. (2010). “Biochemistry of terpenoids: monoterpenes, sesquiterpenes and diterpenes”, in Annual Plant Review ed. M. Wink (Blackwell, Singapore: Wiley), 258–303.
Bartel, D. P. (2004). MicroRNAs: genomics, biogenesis, mechanism, and function. Cell 116, 281–297. doi: 10.1016/S0092-8674(04)00045-5
Bedon, F., Bomal, C., Caron, S., Levasseur, C., Boyle, B., Mansfield, S. D., et al. (2010). Subgroup 4 R2R3-MYBs in conifer trees: gene family expansion and contribution to the isoprenoid- and flavonoid-oriented responses. J. Exp. Bot. 61, 3847–3864. doi: 10.1093/jxb/erq196
Bolger, A. M., Marc, L., and Bjoern, U. (2014). Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics 30, 2114–2120. doi: 10.1093/bioinformatics/btu170
Bulgakov, V. P., and Avramenko, T. V. (2015). New opportunities for the regulation of secondary metabolism in plants: focus on microRNAs. Biotechnol. Lett. 37, 1719–1727. doi: 10.1007/s10529-015-1863-8
Chen, C., Zhong, Y., Yu, F., and Xu, M. (2019). Deep sequencing identifies miRNAs and their target genes involved in the biosynthesis of terpenoids in Cinnamomum camphora. Ind. Crop. Prod. 145:111853. doi: 10.1016/j.indcrop.2019.111853
Cheng, S. S., Chua, M. T., Chang, E. H., Huang, C. G., Chen, W. J., and Chang, S. T. (2009). Variations in insecticidal activity and chemical compositions of leaf essential oils from Cryptomeria japonica at different ages. Bioresource Technol. 100, 465–470. doi: 10.1016/j.biortech.2007.11.060
Cheng, S. S., Lin, C. Y., Chung, M. J., and Chang, S. T. (2012). Chemical composition and antitermitic activity against coptotermes formosanus shiraki of Cryptomeria japonica leaf essential oil. Chem. Biodivers. 9, 352–358. doi: 10.1002/cbdv.201100243
Doron, B., Manda, W., Aaron, G., Marks, D. S., and Chris, S. (2008). The microRNA.org resource: targets and expression. Nucleic Acids Res. 36, D149–D153. doi: 10.1093/nar/gkm995
Fernie, A. R., Aharoni, A., Willmitzer, L., Stitt, M., Tohge, T., Kopka, J., et al. (2011). Recommendations for Reporting Metabolite Data. Plant Cell 23, 2477–2482. doi: 10.1105/tpc.111.086272
Friedlnder, M. R., Mackowiak, S. D., Li, N., Chen, W., and Nikolaus, R. (2012). MiRDeep2 accurately identifies known and hundreds of novel microRNA genes in seven animal clades. Nucleic Acids Res. 40, 37–52. doi: 10.1093/nar/gkr688
Fujihashi, M., Sato, T., Tanaka, Y., Yamamoto, D., Nishi, T., Ueda, D., et al. (2018). Crystal structure and functional analysis of large-terpene synthases belonging to a newly found subclass. Chem. Sci. 9, 3754–3758. doi: 10.1039/c8sc00289d
Grabherr, M. G., Haas, B. J., Yassour, M., Levin, J. Z., and Amit, I. (2013). Trinity: reconstructing a full-length transcriptome without a genome from RNA-Seq data. Nat. Biotechnol. 29, 644–652. doi: 10.1038/nbt.1883
Grassmann, J., Hippeli, S., and Elstner, E. F. (2002). Plant’s defence and its benefits for animals and medicine: role of phenolics and terpenoids in avoiding oxygen stress. Plant Physiol. Bioch. 40, 471–478. doi: 10.1016/S0981-9428(02)01395-5
Grimson, A., Farh, K. H., Johnston, W. K., Garrett-Engele, P., Lim, L. P., and Bartel, D. P. (2007). MicroRNA targeting specificity in mammals: determinants beyond seed pairing. Mol. Cell 27, 91–105. doi: 10.1016/j.molcel.2007.06.017
Guan, X., Pang, M., Nah, G., Shi, X., Ye, W., Stelly, D. M., et al. (2014). miR828 and miR858 regulate homoeologous MYB2 gene functions in Arabidopsis trichome and cotton fibre development. Nat. Commun. 5:3050. doi: 10.1038/ncomms4050
Gupta, O. P., Karkute, S. G., Banerjee, S., Meena, N. L., and Dahuja, A. (2017). Contemporary understanding of miRNA-Based regulation of secondary metabolites biosynthesis in plants. Front. Plant Sci. 8:374. doi: 10.3389/fpls.2017.00374
Han, Q., Shinohara, K., Kakubari, Y., and Mukai, Y. (2003). Photoprotective role of rhodoxanthin during cold acclimation in Cryptomeria japonica. Plant Cell Environ. 26, 715–723. doi: 10.1046/j.1365-3040.2003.01008.x
Izumi, S., and Hirata, T. (1997). Signal transduction against the terpenoid-stimulation in the liverwort cells. Chem. Lett. 26, 1117–1118. doi: 10.1246/cl.1997.1117
Jia, Z., Zhao, B., Liu, S., Lu, Z., and Wang, L. (2020). Embryo transcriptome and miRNA analyses reveal the regulatory network of seed dormancy in Ginkgo biloba. Tree Physiol. 41, 571–588. doi: 10.1093/treephys/tpaa023
Kanehisa, M., Araki, M., Goto, S., Hattori, M., Hirakawa, M., Itoh, M., et al. (2008). KEGG for linking genomes to life and the environment. Nucleic Acids Res. 36, 480–484. doi: 10.1093/nar/gkm882
Kim, S. H., Lee, S. Y., Hong, C. Y., Gwak, K. S., and Choi, I. G. (2013). Whitening and antioxidant activities of bornyl acetate and nezukol fractionated from Cryptomeria japonica essential oil. Int. J. Cosmetic Sci. 35, 484–490. doi: 10.1111/ics.12069
Kksal, M., Zimmer, I., Schnitzler, J. P., and Christianson, D. W. (2010). Structure of isoprene synthase illuminates the chemical mechanism of teragram atmospheric carbon emission. J. Mol. Biol. 402, 363–373. doi: 10.1016/j.jmb.2010.07.009
Krek, A., Grün, D., Poy, M. N., Wolf, R., Rosenberg, L., Epstein, E. J., et al. (2005). Combinatorial microRNA target predictions. Nat. Genet. 37, 495–500. doi: 10.1038/ng1536
Lang, A. (1957). The effect of gibberellin upon flower formation. Proc. Natl. Acad. Sci. U. S. A. 43, 709–717. doi: 10.1073/pnas.43.8.709
Langmead, B., and Salzberg, S. L. (2012). Fast gapped-read alignment with Bowtie 2. Nat. Methods 9, 357–359. doi: 10.1038/nmeth.1923
Lavorgna, G., Guffanti, A., Borsani, G., Ballabio, A., and Boncinelli, E. (1999). TargetFinder: searching annotated sequence databases for target genes of transcription factors. Bioinformatics 15, 172–173. doi: 10.1093/bioinformatics/15.2.172
Li, H., Lin, Y., Chen, X., Bai, Y., Wang, C., Xu, X., et al. (2018). Effects of blue light on flavonoid accumulation linked to the expression of miR393, miR394 and miR395 in longan embryogenic calli. PLoS One 13:e191444. doi: 10.1371/journal.pone.0191444
Licursi, V., Conte, F., Fiscon, G., and Paci, P. (2019). MIENTURNET: an interactive web tool for microRNA-target enrichment and network-based analysis. BMC Bioinformatics 20:545. doi: 10.1186/s12859-019-3105-x
Lin, H. Y., Cheng, S. S., and Chang, S. T. (2005). Chemical composition and antifungal activity of essential oils from different tissues of Japanese Cedar (Cryptomeria japonica). J. Agric. Food Chem. 53, 614–619. doi: 10.1021/jf0484529
Liu, J., Yuan, Y., Wang, Y., Jiang, C., Chen, T., Zhu, F., et al. (2017). Regulation of fatty acid and flavonoid biosynthesis by miRNAs in Lonicera japonica. RSC Adv. 7, 35426–35437. doi: 10.1039/C7RA05800D
Lu, X. Y., and Huang, X. L. (2008). Plant miRNAs and abiotic stress responses. Biochem. Biophys. Res. Commun. 368, 458–462. doi: 10.1016/j.bbrc.2008.02.007
Mafu, S., Hillwig, M. L., and Peters, R. J. (2011). A novel labda-7,13e-dien-15-ol-producing bifunctional diterpene synthase from Selaginella moellendorffii. Chembiochem 12, 1984–1987. doi: 10.1002/cbic.201100336
Matsunaga, T., Hasegawa, C., Kawasuji, T., Suzuki, H., Saito, H., Sagioka, T., et al. (2000). Isolation of the antiulcer compound in essential oil from the leaves of Cryptomeria japonica. Biol. Pharm. Bull. 23, 595–598. doi: 10.1248/bpb.23.595
Oku, H., Inafuku, M., Ishikawa, T., Takamine, T., Ishmael, M., and Fukuta, M. (2015). Molecular cloning and biochemical characterization of isoprene synthases from the tropical trees Ficus virgata, Ficus septica, and Casuarina equisetifolia. J. Plant Res. 128, 849–861. doi: 10.1007/s10265-015-0740-9
Otomo, K., Kenmoku, H., Oikawa, H., König, W. A., Toshima, H., Mitsuhashi, W., et al. (2004). Biological functions of ent- and syn-copalyl diphosphate synthases in rice: key enzymes for the branch point of gibberellin and phytoalexin biosynthesis. Plant J. 39, 886–893. doi: 10.1111/j.1365-313X.2004.02175.x
Pantaleo, V., Szittya, G., Moxon, S., Miozzi, L., Moulton, V., Dalmay, T., et al. (2010). Identification of grapevine microRNAs and their targets using high−throughput sequencing and degradome analysis. Plant J. 62, 960–976. doi: 10.1111/j.1365-313X.2010.04208.x
Parreira, J. R., Cappuccio, M., Balestrazzi, A., Fevereiro, P., and Araújo, S. (2021). MicroRNAs expression dynamics reveal post-transcriptional mechanisms regulating seed development in Phaseolus vulgaris L. Hortic. Res. 8:18. doi: 10.1038/s41438-020-00448-0
Reddy, V. A., Wang, Q., Dhar, N., Kumar, N., Venkatesh, P. N., Rajan, C., et al. (2017). Spearmint R2R3-MYB transcription factor MsMYB negatively regulates monoterpene production and suppresses the expression of geranyl diphosphate synthase large subunit (MsGPPS.LSU). Plant Biotechnol. J. 15, 1105–1119. doi: 10.1111/pbi.12701
Sasaki, K., Ohara, K., and Yazaki, K. (2005). Gene expression and characterization of isoprene synthase from Populus alba. FEBS Lett. 579, 2514–2518. doi: 10.1016/j.febslet.2005.03.066
Silver, G. M., and Fall, R. (1991). Enzymatic synthesis of isoprene from dimethylallyl diphosphate in aspen leaf extracts. Plant Physiol. 97, 1588–1591. doi: 10.1104/pp.97.4.1588
Song, Y., Ma, K., Ci, D., Zhang, Z., and Zhang, D. (2013). Sexual dimorphism floral microRNA profiling and target gene expression in andromonoecious poplar (Populus tomentosa). PLoS One 8:e62681. doi: 10.1371/journal.pone.0062681
Trapnell, C., Williams, B. A., Pertea, G., Mortazavi, A., Kwan, G., Baren, M., et al. (2010). Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nat. Biotechnol. 28, 511–515. doi: 10.1038/nbt.1621
Vashisht, I., Mishra, P., Pal, T., Chanumolu, S., Singh, T. R., and Chauhan, R. S. (2015). Mining NGS transcriptomes for miRNAs and dissecting their role in regulating growth, development, and secondary metabolites production in different organs of a medicinal herb, Picrorhiza kurroa. Planta 241, 1255–1268. doi: 10.1007/s00425-015-2255-y
Vranová, E., Coman, D., and Gruissem, W. (2013). Network analysis of the MVA and MEP pathways for isoprenoid synthesis. Annu. Rev. Plant Biol. 64, 665–700. doi: 10.1146/annurev-arplant-050312-120116
Wang, S. Y., Lai, W. C., Chu, F. H., Lin, C. T., Shen, S. Y., and Chang, S. T. (2006). Essential oil from the leaves of Cryptomeria japonica acts as a silverfish (Lepisma saccharina) repellent and insecticide. J. Wood Sci. 52, 522–526. doi: 10.1007/s10086-006-0806-3
Wu, J., and Zhao, W. M. (2010). New sesquiterpene and triterpene from the fruits of Cryptomeria fortunei. J. Asian Nat. Prod. Res. 12, 382–387. doi: 10.1080/10286021003785565
Xie, Y., Wang, J., Yang, F., and Lei, C. (2011). Comparative analysis of essential oil components of two Cryptomeria species from China. Ind. Crop. Prod. 34, 1226–1230. doi: 10.1016/j.indcrop.2011.04.016
Xu, M., Hillwig, M. L., Prisic, S., Coates, R. M., and Peters, R. J. (2010). Functional identification of rice syn-copalyl diphosphate synthase and its role in initiating biosynthesis of diterpenoid phytoalexin/allelopathic natural products. Plant J. 39, 309–318. doi: 10.1111/j.1365-313X.2004.02137.x
Xu, X., Yin, L., Ying, Q., Song, H., Xue, D., Lai, T., et al. (2013). High-Throughput sequencing and degradome analysis identify miRNAs and their targets involved in fruit senescence of Fragaria ananassa. PLoS One 8:e70959.
Xu, Y., Zhu, C., Xu, C., Sun, J., Grierson, D., Zhang, B., et al. (2019). Integration of metabolite profiling and transcriptome analysis reveals genes related to volatile terpenoid metabolism in finger citron (C. medica var. sarcodactylis). Molecules 24:2564. doi: 10.3390/molecules24142564
Xu, Y. H., Wang, J. W., Wang, S., Wang, J. Y., and Chen, X. Y. (2004). Characterization of GaWRKY1, a cotton transcription factor that regulates the sesquiterpene synthase gene (+)-delta-cadinene synthase-A. Plant Physiol. 135, 507–515. doi: 10.1104/pp.104.038612
Yao, S., Tang, C. P., Ke, C. Q., and Ye, Y. (2008). Abietane diterpenoids from the bark of Cryptomeria fortunei. J. Nat. Prod. 71, 1242–1246. doi: 10.1021/np8002063
Ye, Y., Wang, J., Ni, Z., Meng, X., and Xu, L. A. (2020). Small RNA and degradome sequencing reveal roles of miRNAs in strobilus development in masson pine (Pinus massoniana). Ind. Crop. Prod. 154:112724. doi: 10.1016/j.indcrop.2020.112724
Yu, Z. X., Li, J. X., Yang, C. Q., Hu, W. L., Wang, L. J., and Chen, X. Y. (2012). The jasmonate-responsive AP2/ERF transcription factors AaERF1 and AaERF2 positively regulate artemisinin biosynthesis in Artemisia annua L. Mol. Plant 5, 353–365. doi: 10.1093/mp/ssr087
Yu, Z. X., Wang, L. J., Zhao, B., Shan, C. M., Zhang, Y. H., Chen, D. F., et al. (2015). Progressive regulation of sesquiterpene biosynthesis in Arabidopsis and patchouli (Pogostemon cablin) by the miR156-targeted SPL transcription factors. Mol. Plant 8, 98–110. doi: 10.1016/j.molp.2014.11.002
Zhang, Y., Zhu, L., Xue, J., Yang, J., Hu, H., Cui, J., et al. (2021). Selection and verification of appropriate reference genes for expression normalization in Cryptomeria fortunei under abiotic stress and hormone treatments. Genes 12:791. doi: 10.3390/genes12060791
Zhao, H., Shi, X., Shen, C., Chen, C., Qu, C., Patel, G., et al. (2021). Comparative analysis of miRNA expression profiles provides insight into regulation of biosynthesis of flavonoids and terpenoids between two varieties of Toona sinensis sprouts. J. Plant Biol. doi: 10.1007/s12374-021-09321-6
Keywords: cold acclimation, transcriptome, degradome, diterpenoid biosynthesis, metabolome analysis, miRNA target
Citation: Zhang Y, Cui J, Hu H, Xue J, Yang J and Xu J (2021) Integrated Four Comparative-Omics Reveals the Mechanism of the Terpenoid Biosynthesis in Two Different Overwintering Cryptomeria fortunei Phenotypes. Front. Plant Sci. 12:740755. doi: 10.3389/fpls.2021.740755
Received: 13 July 2021; Accepted: 23 August 2021;
Published: 29 September 2021.
Edited by:
Xue-Rong Zhou, Agriculture and Food, Commonwealth Scientific and Industrial Research Organisation (CSIRO), AustraliaReviewed by:
Henrik Toft Simonsen, Technical University of Denmark, DenmarkAtsushi Fukushima, RIKEN, Japan
Copyright © 2021 Zhang, Cui, Hu, Xue, Yang and Xu. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Jin Xu, eGppbmhzaEBuamZ1LmVkdS5jbg==