Disentangling the Potential Functions of miRNAs in the Synthesis of Terpenoids during the Development of Cinnamomum burmannii Leaves

Abstract

The essential oil of Cinnamomum burmannii (Nees and T. Nees) Blume is rich in monoterpenes and sesquiterpenes. The post-transcriptional regulatory mechanisms controlling the expression of terpenoid-related genes have not yet been clarified in C. burmannii. Here, we conducted a metabolomic analysis of the leaves of C. burmannii across four developmental stages using gas chromatography–mass spectrometry. We also identified miRNAs and their target genes involved in terpenoid biosynthesis using small RNA sequencing. A total of 135 differentially expressed metabolites were detected, including 65 terpenoids, 15 aldehydes, and 13 benzenes. A total of 876 miRNAs from 148 families were detected, among which 434 miRNAs were differentially expressed, including three known miRNAs and 431 novel miRNAs. Four miRNAs (gma-miR5368, novel_miR_377, novel_miR_111, and novel_miR_251) were predicted to regulate the expression of four differential expressed genes involved in the monoterpenoid and sesquiterpenoid synthesis. miRNAs families miR396, miR5185, and miR9408 were predicted to play diverse regulatory roles in monoterpenoid and sesquiterpenoid synthesis during the leaf development of C. burmannii. The results of our study shed new light on the roles of regulatory genes in terpenoid biosynthesis. Our findings also have implications for the further promotion of essential oil production using the leaves of C. burmannii.

1. Introduction

Cinnamomum burmannii (Nees and T. Nees) Blume is a member of the family Lauraceae that occurs in southeastern China, Indonesia, and the Philippines [1,2]. The leaves of C. burmanni are rich in terpenoids, including monoterpenoids, such as borneol, bornyl acetate, eucalyptol, limonene, pinene, and phellandrene; sesquiterpenes, such as caryophyllene, copaene, germacrene D, and phellandrene [3,4,5]; and other aromatic compounds, such as benzaldehyde, cinnamic aldehyde, and decanal [6]. Terpenes and terpenoids play important roles in attraction to pollinators and seed dispersal vectors via floral scents and fruit aromas [7], plant defense via their insecticidal and antimicrobial effects [8,9], and monoterpenes and sesquiterpenes in the essential oil of C. burmannii have antibacterial and antimicrobial properties [2,5]. For example, the extract of C. burmannii bark in Indonesia can effectively repel five common foodborne pathogenic bacteria (Bacillus cereus, Listeria monocytogenes, Staphylococcus aureus, Escherichia coli, and Salmonella anatum) [10]. The essential oil of C. burmannii leaves thus has important applications in the cosmetic and medical industries.

Terpenes and terpenoids are synthesized via the mevalonate acid (MVA) pathway and the 2-C-methyl-D-erythritol 4-phosphate (MEP) pathway in all land plants [11,12]. Previous studies have shown that mevalonate kinase (MVK) in the MVA pathway, as well as 1-deoxy-D-xylulose-5-phosphate synthase (DXS), 1-deoxy-deoxyxylulose-5-phosphate reductoisomerase (DXR), and (E)-4-hydroxy-3-methylbut-2-enyl-diphosphate synthase in the MEP pathway, are differentially expressed; these enzymes might thus play key roles in determining the chemotypes of Cinnamomum porrectum [13] and Cinnamomum camphora [14,15]. The expression of CbIDI1, which encodes isopentenyl–diphosphate isomerase (IDI), varies among vegetative organs in C. burmannii, suggesting that it plays a key role in borneol biosynthesis [16]. Terpene synthases (TPSs) have important functions in subsequent steps in the biosynthesis of monoterpenes and sesquiterpenes. The functions of seven TPS genes (CbTPS1–7) have been clarified in C. burmannii [3,17]. However, the molecular mechanisms underlying the regulation of the expression of upstream and downstream genes involved in the biosynthesis of monoterpenoids and sesquiterpenes in C. burmannii remain unclear.

microRNAs (miRNAs) are non-protein coding sequences that range from 20 to 24 nucleotides in length, and they are involved in the post-transcriptional regulation of gene expression [18]. miRNAs play key roles in regulating the growth and development of plants, as well as the synthesis of secondary metabolites [19,20]. miRNAs also positively and negatively regulate the expression of target genes [21,22]. Plant miRNA genes are often highly similar due to tandem and segmental duplication [23]. Some conserved miRNAs, such as miR396, have been shown to play key regulatory roles in Pelargonium spp. [24] and Ginkgo biloba [25]. Other non-conserved miRNAs, such as miR4995, miR5021, and miR6300, have been shown to regulate terpenoid biosynthesis in the leaves of C. camphora [26]. Few miRNAs that regulate the biosynthesis of monoterpernoids and sesquiterpenoids have been identified. In addition, our understanding of the diversity and functions of miRNAs in members of the genus Cinnamomum is limited. The roles of miRNAs and their effects on the expression of target genes during the development of C. burmanii leaves remain unclear.

Here, gas chromatography–mass spectrometry and small RNA (sRNA) sequencing were used to analyze the metabolome and to identify miRNAs, respectively, during the four stages of leaf development in C. burmannii. We then analyzed co-expression patterns among differentially expressed metabolites (DEMs) and differentially expressed miRNAs (DEmiRNAs) to identify the key miRNAs and their target genes involved in regulating the synthesis of monoterpenoids and sesquiterpenoids. Our study is the first to identify miRNAs in C. burmannii and analyze their expression patterns. Generally, the results of our study shed new light on the roles of miRNAs in terpenoid biosynthesis and will have implications for further promotion of essential oil production using the leaves of C. burmannii.

2. Materials and Methods

2.1. Collection of Leave Samples

Leaves of borneol-type C. burmannii were collected from a six-year-old plant on 16 February in the dry season of 2022. The plant samples were collected from Tianluhu Forestry Park (113°25′54″ N 23°15′33.4″ E, Guangzhou, China), where collection permission was obtained from the department of resource administration of Longdong Forest Farm. We collected C. burmannii leaves in four developmental stages: the leaf at 7 d growth after its sprout (CBTS1), at 14 d growth (CBTS2), at 21 d growth (CBTS3), and at 28 d growth (CBTS4). Three replicate leaf samples were obtained for each developmental stage for analysis of secondary metabolites and sRNA sequencing. All samples were stored at −80 °C prior to analyses.

2.2. Extraction and Detection of Metabolites via GC–MS

Each sample (500 mg) was placed in a 20-mL headspace bottle, and 20 μL of 2-octanol was added as the internal standard. Each sample was subjected to GC–MS analysis using an Agilent 7890 gas chromatograph system coupled with a 5977B mass spectrometer (Agilent Technologies, Palo Alto, CA, USA). All metabolites were further separated on a DB-Wax column (30 m × 250 μm × 0.25 μm, Agilent Technologies, Palo Alto, CA, USA), and helium was used as the carrier gas. The front inlet purge flow was 3 mL min⁻¹, and the gas flow rate through the column was 1 mL min⁻¹. The initial temperature remained at 40 °C for 4 min; it was then increased to 245 °C at a rate of 5 °C min⁻¹ for 5 min. The scanning range was 20–400 m/z. Chroma TOF4.3X software (LECO Corporation, St. Joseph, MI, USA) and the LECO-NIST database were used to identify secondary metabolites. The relative content of each compound was determined via comparison with the peak area of 2-octanol.

2.3. Differential Expression Analysis of Secondary Metabolites

Principal component analysis (PCA) was conducted in the FactoMineR package [27] in R version 3.3.2 [28]. Differentially expressed metabolites were identified with false discovery rate (FDR) > 1.50, p-value < 0.05 (Student’s t-test), and variable importance in projection value > 1.

2.4. WGCNA Analysis

We then performed a weighted correlation network analysis (WGCNA) for correlations between the contents of the top eleven terpenoids (six most abundant monoterpenoids and five sesquiterpenoids) and co-expressed genes involved in the biosynthesis of terpenoids. A trait matrix based on the relative content detected above was prepared for the six most abundant monoterpenoids and five sesquiterpenoids. Additionally, the expression information of co-expressed genes refers to Hou et al. (2023) [29], within which methods of quantifying gene expression and detecting differentially expressed genes were described across the four developmental stages of C. burmanii leaves. The correlation of co-expressed genes and these metabolites were analyzed in a weighted correlation network analysis (WGCNA) implemented in the R package WGCNA (version 1.42) [30]. A hierarchal clustering tree based on the co-expressed genes was constructed using the Dynamic Tree Cut R package [31]. Genes significantly correlated with the abundance of the 11 terpenoids were grouped into different modules and further annotated based on the GO and KEGG database, respectively. Combined with the results of DEG detection, the obtained genes from WGCNA analysis with differential expression across the four developmental stages were selected for subsequent analysis. Gene Ontology (GO, http://www.geneontology.org, accessed on 3 September 2022) and Kyoto Encyclopedia of Genes and Genomes (KEGG, http://www.genome.jp/kegg/, accessed on 3 September 2022) enrichment analyses were conducted on the identified target genes using BLAST v2.2.26 [32] with the parameters (-b 100 -v 100 -e 1e-5 -m 7 -a 2).

2.5. Construction and Sequencing of the sRNA Library

From 1 g fresh C. burmanii leaves, 2 μg of RNA was extracted using the Trizol reagent (Thermo Fisher, Shanghai, China). RNA integrity was assessed using RNA Nano 6000 Kit of Agilent Bioanalyzer 2100 system (Agilent Technologies, Palo Alto, CA, USA). The purity and concentrations of RNAs were examined by 1% agar–gel electrophoresis and NanoDrop NC 2000, respectively (Thermo Scientific, Carlsbad, CA, USA). Prior to transcriptome sequencing, adapters were ligated, cDNA was synthesized via reverse transcription, and sequences were amplified via PCR. PCR products 50 bp in length were then subjected to polyacrylamide gel electrophoresis. The AMPure XP system (Beckman Coulter, Shanghai, China) was used to purify the PCR products. A cBot Cluster Generation System was used to cluster the index-coded samples with a TruSeq PE Cluster Kit v4-cBot-HS (Illumina, Boston, MA, USA) per the manufacturer’s instructions. A total of 12 sRNA libraries were established, and they were further sequenced by the Illumina platform Hiseq 2500 (Illumina, San Diego, CA, USA) following the manufacturer’s instructions. The data of sRNA library sequencing were deposited in the SRA database under the project number PRJNA892262.

2.6. Identification of miRNAs and Target Genes

Clean reads were obtained from raw data by removing reads with adapters, reads containing poly-N, and low-quality reads. Reads less than 18 nucleotides and greater than 30 nucleotides were also removed. The clean data were compared against the Silva (http://www.arb-silva.de, accessed on 4 September 2022), GtRNAdb (http://lowelab.ucsc.edu/GtRNAdb, accessed on 4 September 2022), Rfam (http://rfam.xfam.org, accessed on 4 September 2022), and Repbase (http://www.girinst.org/repbase, accessed on 4 September 2022) databases using Bowtie v1.0.0 [33] with the parameter (-v 0). Ribosomal RNA, transfer RNA, small nuclear RNA, small nucleolar RNA, and repeats were removed. The remaining reads were mapped to the reference genome of C. burmannii (Hou et al. unpublished). Searches of the mapped reads were conducted against the miRBase v22.1 database [34], and these searches included the regions two nucleotides upstream and five nucleotides downstream of the reads. Mapped reads that were perfect matches or with one nucleotide mismatched were classified as known miRNAs. Mapped reads with one strand of miRNA precursor classified as known miRNAs were classified as novel miRNAs with a “novel_miR” prefix. MiRDeep2 v2.0.5 was used to distinguish miRNA and siRNA for each sample. Based on the distribution information of reads on the precursor sequence with regard to the characteristics of miRNA structuring (“mature”, “star”, and “loop”) using MiRDeep2 [35] with the parameters (g -1 -l 250 -b 0). The energy information of the precursor structure was estimated using Randfold v2.0 [36] with the parameter (-s 99). The prediction of new miRNA was finally realized by scoring the sequences and secondary structures using Bayesian model implemented in Randfold. MiRDeep2 was also used to classify the miRNAs into different families. TargetFinder v1.6 [37] with the parameter (-c 3) and the reference genome of C. burmanii were used to predict target genes. These target genes were further annotated based on the GO and KEGG database as the methods shown above.

2.7. Detection of DEmiRNAs and Co-Expression Analysis with Terpenoids

The quantifier.pl program installed in MiRDeep2 was used to characterize the expression of miRNAs in each sample. To quantify the expression of all known and novel miRNAs, we imported the sequencing data for each sample into the program with the parameters (-g -1 250 -b 0). sRNAs were mapped back to their precursor sequences, and the read count for each miRNA was determined. The transcripts per million (TPM) method was used to quantify the expression of miRNAs [38], where TPM = the reads mapped to total miRNAs divided by the reads mapped to a certain miRNA × 1,000,000. The DESeq2 package version 1.10. 1 in R [39,40] was used to conduct differential expression analysis between developmental stages. DEmiRNAs were identified by assuming a negative binomial distribution and using the following criteria: |log2(fold change (FC))| ≥ 0.58 and p-value ≤ 0.01. The Benjamini–Hochberg approach was used to adjust the p-values. miRNAs with their target genes involved in mono- and sesquiterpenoids were further tested with the varied contents of six monoterpenoids and five sesquiterpenoids across four developmental stages using a Pearson’s correlation method. The analysis was implemented using corrplot package installed in R [28] with the setting parameters (correlation index > 0.7, p-value < 0.05).

3. Results

3.1. Changes in Metabolites across Different Developmental Stages

GC–MS was used to characterize changes in secondary metabolites across four developmental stages of C. burmannii leaves. According to the results of PCA analysis (Figure 1A), the metabolites of CBTS1 and CBTS2 were similar, and the metabolites of CBTS1 and CBTS2 differed from those of CBTS3. The metabolites of CBTS4 were markedly different from those of CBTS1, CBTS2, and CBTS3. GC–MS revealed 135 DEMs, including 65 terpenoids, 15 aldehydes, 13 benzenes, 12 esters, 6 alkanes, 3 alcohols, 3 ketones, 3 naphthalenes, 3 phenols, and 12 others (Supplementary Material Table S1). Further analysis identified the three most abundant numbers of DEMs: 97 DEMs (CBTS1 vs. CBTS4, Figure 1B), 94 DEMs (CBTS2 vs. CBTS4), and 93 DEMs (CBTS1 vs. CBS3) (Supplementary Material Table S2). A Venn diagram revealed that 67 DEMs were shared between CBTS1 and CBTS3, between CBTS1 and CBTS4, and between CBTS2 and CBTS4 (Figure 1B). These shared DEMs comprised 37 terpenoids, 8 benzenes, and 5 aldehydes. Among the 37 terpenoids detected, six monoterpenoids (borneol, bornyl acetate, β-phellandrene, D-limonene, β-pinene, and eucalyptol) and five sesquiterpenoids (nerolidol, γ-muulrolene, β-elemene, aromandenrene, and copaene) sharply increased from CBTS1 to CBTS4 (Figure 1C). A phylogeny of co-expressed genes indicates 12 modules, of which three modules (represented by the colors plum, salmon, and pale turquoise) accounted for most co-expressed genes involved in the synthesis of the 11 terpenoids (Figure 1D). Further, GO annotation of the co-expressed genes in the three modules reveal that a total of 14 genes (salmon-6, pale turquoise-2, and plum-6) were involved in terpenoid synthesis (Figure 1F). KEGG annotation of the co-expressed genes reveals that a total of 34 co-expressed (salmon-5, pale turquoise-9, and plum-20) participated in the monoterpenoid and sesquiterpenoid biosynthesis (Figure 1G).

3.2. Characteristics and Classification of the Identified miRNAs

sRNA sequencing was used to identify miRNAs in the leaves of C. burmannii across four developmental stages. A total of 876 miRNAs were detected, and five of these miRNAs were perfect matches with mature miRNAs in the miRbase. The remaining 871 miRNAs had one strand of miRNA precursors; these were thus classified as novel miRNAs (Supplementary Material Table S3). The lengths of the five known miRNAs ranged from 19 to 24 nucleotides, and two miRNAs were 20 nucleotides in length. The lengths of the 871 novel miRNAs ranged from 18 to 25 nucleotides. Most of these miRNAs were 24 nucleotides in length (647 miRNAs, 74.28%), and only two miRNAs (0.23%) were 18 nucleotides in length (Figure 2A). The 876 identified miRNAs were further classified into 148 families (Supplementary Material Table S4), and the five largest families were miR396, miR159, miR482, miR166, and miR171 (Figure 2B). A total of 1997 target genes were predicted based on the DEmiRNAs detected (Supplementary Material Table S5). GO and KEGG analyses were conducted to annotate the functions of these target genes (Figure 2C,D). The results of GO annotation showed that target genes were involved in terpene synthesis, e.g., monoterpene biosynthetic process (8 genes, GO ID:0043693), sesquiterpene biosynthetic process (5 genes, GO:0051762). The results of KEGG annotation showed that target genes were involved in terpene synthesis, e.g., terpenoid backbone biosynthesis (4 genes, KEGG pathway ID ko00900) and sesquiterpenoid and triterpenoid biosynthesis (4 genes, ko00909).

3.3. Regulatory Networks of miRNAs and Target Genes

According to the above results, we constructed several regulatory networks of miRNAs and target genes involved in the synthesis of monoterpenoids and sesquiterpenoids (Figure 3). Three main types of miRNA–target gene regulatory networks were constructed. In the first type of network, one miRNA regulates the expression of one gene; for example, novel_miR_492, which belongs to the miR4243 family, was predicted to regulate the expression of the Cbur03G002600 gene. In the second type of network, one miRNA regulates the expression of several genes; for example, novel_miR_377, which is in the miR5185 family, was predicted to regulate the expression of the four genes Cbur04G0023450, Cbur12G012200, Cbur06G000240, and Cbur02G017610. In the third type of network, several miRNAs regulate the expression of several genes; for example, novel_miR_277, novel_miR_597, and novel_miR_851, all of which are in the miR396 family, were predicted to regulate the expression of 15 target genes, and novel_miR_231 and novel_miR_266, both of which are in the miR396 family, were predicted to regulate the expression of Cbur12G016860 and Cbur07G013130. In another sample, gma-miR5368, which is in the miR9408 family, was predicted to regulate the expression of 13 target genes, and novel_miR_111 (family unknown), novel_miR_271 (miR478 family), and novel_miR_685 (miR399 family) were predicted to regulate the expression of Cbur12G010760, Cbur09G007390, and Cbur02G035810.

3.4. DEmiRNAs and Co-Expression Analysis with Terpenoids

The TPM method was used to quantify the expression of miRNAs, and DESeq2 was used to identify DEmiRNAs between each pair of developmental stages. Among all 876 miRNAs, 434 DEmiRNAs (three known and 431 novel) were identified (Supplementary Material Table S6). Of the DEmiRNAs detected, the number of DEmiRNAs between CBTS1 and CBTS4 was the largest (252), and these comprised 95 up-regulated and 157 down-regulated DEmiRNAs (Figure 4A). Additionally, 239 and 232 DEmiRNAs were identified between CBTS2 vs. CBTS4 and CBTS1 vs. CBS3, respectively, resulting in the second and third abundant DEmiRNAs detected. A Venn diagram revealed that 88 DEmiRNAs were shared between CBTS1 and CBTS3, between CBTS1 and CBTS4, and between CBTS2 and CBTS4 (Figure 4B). In light of the target genes related to mono- and sesquiterpenoid synthesis, four DEmiRNAs (gma_miR5368, nov_miR_377, nov_miR_111, and nov_miR_251) were tested with 11 terpenoids (six monoterpenoids and five sesquiterpenoids) with significant correlation (index > 0.7, p-value < 0.05) evaluated by Pearson correlation analysis (Figure 4C). At the stages of CBTS1, the expression of gma-miR5368 was shown to possess a significantly negative correlation with eight terpenoids (β-phellandrene, D-limonene, eucalyptol, nerolidol, γ-muulrolene, β-elemene, aromandenrene, and copaene). nov_miR_111 and nov_miR_251 had a negative correlation with the content of borneol and bornyl acetate, while nov_miR_377 revealed a positive correlation with these two compounds. At the stages of CBTS2, gma-miR5368, nov_miR_377, and nov_miR_111 revealed a significantly negative correlation with all the terpenoids, while nov_miR_251 was shown to possess a significantly positive correlation with all the terpenoids. At the stages of CBTS3, gma-miR5368 and nov_miR_377 were highly positively and negatively correlated with the expression of five terpenoids (bornyl acetate, eucalyptol, nerolidol, β-elemene, and copaene), respectively. Expression of nov_miR_111 was shown to possess a significantly positive relationship with the content of borneols, and expression of nov_miR_251 was negatively correlated with the content of borneol acetate. At the stage of CBTS4, the expression of nov_miR_251 was shown to have a significantly negative correlation with all the content of six monoterpenoids while having a positive correlation with nerolidol. For nov_miR_377, however, its expression was negatively correlated with the remaining four sesquiterpenoids (γ-muulrolene, β-elemene, aromandenrene, and copaene). Expression of gma-miR5368 and nov_miR_377 was shown to possess a significantly positive and negative correlation with the content of nerolidol and copaene. Thus, correlation patterns of the four miRNAs and the content of the 11 terpenoids were shown to be varied across the four developmental stages of C. burmannii leaves.

3.5. Regulation of Monoterpenoid and Sesquiterpenoid Synthesis by miRNAs

Based on the results shown above, four miRNAs were predicted to regulate the expression of four differential expressed genes in relation to monoterpenoid and sesquiterpenoid synthesis (Figure 5). In the upstream pathway of monoterpenoid and sesquiterpenoid biosynthesis, the expression of Cbur09G007390 encoding MVK in the MVA pathway decreased from CBTS1 to CBTS4. MicroRNA gma-miR5368 was predicted to regulate the expression of Chur09G007390, with its expression increased from CBTS1 to CBTS3 but slightly decreased in CBS4. In the MEP pathway, novel_miR_377 was predicted to regulate the expression of the genes encoding enzyme DXS. The expression of Cbur04G023450 markedly increased from CBTS1 to CBTS3 but slightly decreased in CBTS4. By contrast, the expression of novel_miR_377 was high in CBS1, slightly higher in CBTS2, and markedly decreased from CBTS2 to CBTS4. The expression of Cbur02G035810 encoding HDR markedly increased from CBTS1 to CBTS4. The expression of Cbur02G009940 encoding enzyme IDI slightly increased from CBTS1 to CBTS3 but decreased in CBTS4. The opposite pattern was observed in the expression of novel_miR_111; the only exception was that a marked increase in its expression was observed in CBTS4. The expression of the gene Cbur09G005200 encoding a mono-TPS enzyme decreased from CBTS1 to CBTS4, while the expression of novel_miR_251 gradually increased from CBTS1 to CBTS4.

4. Discussion

The essential oil of C. burmannii leaves is rich in terpenoids. We detected a total of 135 DEMs across four developmental stages of d-borneol chemotype C. burmannii leaves, and these DEMs included several of the most abundant monoterpenoids, such as borneol, borneol acetate, d-limonene, eucalyptol, α-phellandrene, α-pinene, camphor, and γ-terpene, and sesquiterpenes, such as nerolidol, γ-muurolene, β-elemene, and copaene; all of these secondary metabolites have been documented in the mature leaves of C. burmannii in previous studies [4,17]. The high content of monoterpenoids and sesquiterpenes is likely associated with their roles in plant defense during the leaf development of C. burmannii, as volatile terpenes and terpenoids have insecticidal and antimicrobial activities that protect plant leaves from damage [8,9]. The essential oil of Cinnamomum zeylanicum, which is composed of monoterpenes, such as α-pinene and α-phellandrene, and sesquiterpenes, such as caryophyllene oxide and α-humulene, effectively repels adult granary weevils (Sitophilus granaries) [41]. We also found that the content of 15 aldehydes, such as hexanal, isoneral, octanal, and tridecanal, varied across the four developmental stages of C. burmannii leaves (Supplementary Material Table S1). Previous studies have shown that these aldehydes are potent invertebrate repellents [42,43]. Thus, these aldehydes might also be involved in plant defense in C. burmannii, but additional studies are needed to clarify whether these aldehydes have such a function.

The results of our study reveal that miRNAs depress a gene encoding MVK in the MVA pathway (Figure 5). MVK is a key enzyme in monoterpene biosynthesis in different chemotypes of C. camphora [14] and C. porrectum [13]. Specifically, gma-miR5368, which is in the miR9408 family, was predicted to regulate the expression of Cbur09G007390. This miRNA has been shown to be involved in the regulation of terpenoid synthesis in previous studies. A previous study shows gma-miR5368 targeted the GAPN gene, which is involved in the glycolysis/gluconeogenesis metabolic pathway in the developing seeds of Camellia oleifera [44]. Another study shows that gma-miR5368 was involved in response to chilling and cold stress in Musa itinerans [45]. Our findings indicate that gma-miR5368 has opposite effects on the expression of Cbur09G007390 (Figure 5). Additionally, the expression of gma-miR5368 was shown to possess a significantly negative correlation with three monoterpenes (β-phellandrene, D-limonene, eucalyptol) and five sesquiterpenes (nerolidol, γ-muulrolene, β-elemene, aromandenrene, and copaene) at CBTS1 (Figure 4C), and also to possess significantly negative correlation with two monoterpenes (bornyl acetate, eucalyptol) and three sesquiterpenes (nerolidol, β-elemene, and copaene). Thus, gma-miR5368, as a member of the family miR9408, is likely to play important roles in regulating the expression of key genes in the MVA pathway during leaf development in C. burmannii.

According to our findings, novel_miR_377, which is in the miR5185 family, was predicted to depress the expression of Cbur04G023450 encoding enzyme DXR in the MEP pathway (Figure 5). Enzyme DXS are key enzymes involved in the regulation of terpenoid metabolism in several Cinnamomum species, such as. C. burmannii [4,16], C. camphora [14,46,47,48], and C. porrectum [13]. Members of the miR5185 family are widespread in plants such as Allium cepa [49], Brachypodium distachyon [50], and Ribes nigrum [51]. Few studies have examined the regulatory roles of miR5185. Members of the miR5185 family have been shown to regulate the expression of genes encoding NADH dehydrogenase subunit and ubiquitin-conjugating enzyme E2 protein [49]. Additionally, members of the miR535 family are involved in the regulation of ethylene-induced fruit ripening in bananas [52,53]. Given the results in the present study, however, novel_miR_377 revealed a significantly negative correlation with all the 11 mono- and sesquiterpenoids at CBTS2, while a significant positive correlation with borneol at CBTS1 and 4 sesquiterpenoids (γ-muulrolene, β-elemene, aromandenrene, and copaene) (Figure 4C). Thus, novel_miR_377, as a member of the family miR5185, might play a role in regulating monoterpenoid and sesquiterpenoid synthesis during the development of C. burmannii leaves.

Novel_miR_251, which belongs to the miR396 family, was predicted to depress the expression of Cbur09G005200, which encodes a mono-TPS enzyme that functions downstream of the MVA and MEP pathways (Figure 5). The miR396 family is the key conserved miRNA family [54], and several miRNAs in these families have been documented in various plant species, such as Camellia sinensis [55], Juglans regia [56], Phaseolus vulgaris [57], and Pinus massoniana [22]. Members of the miR396 family regulate the biosynthesis of flavonoids [58], grain size and shoot architecture [59], and cell proliferation during leaf development [60,61]. In recent years, members of the miR396 family have been shown to play a role in regulating terpenoid synthesis. Members of the miR396 family have been shown to regulate the expression of the gene encoding DXS and terpene biosynthesis in rose-scented geranium—Pelargonium spp. [24]. Members of the miR396 family regulate secondary metabolite biosynthesis by affecting the expression of genes encoding AP2 transcription factors in G. biloba [25]. Our findings indicate that nov_miR_251 had a significantly negative correlation with the content of borneol and bornyl acetate at CBTS1, while shown to possess a significantly positive correlation with all 11 terpenoids (Figure 4C). At CBTS3, the expression of nov_miR_251 was strictly negatively correlated with the content of borneol acetate, while the expression of nov_miR_251 was shown to have significantly negative correlation with all the content of six monoterpenoids. These findings indicate that, nov_miR_251, as a member of the miR396 family, regulate the biosynthesis of monoterpenes and sesquiterpenes in C. burmannii.

5. Conclusions

To the best of our knowledge, our study is the first to clarify the regulatory roles of miRNAs in terpenoid synthesis in C. burmannii. A total of 876 microRNAs from 148 families were detected, and 434 miRNAs, including three known and 431 novel miRNAs, were differentially expressed. Four miRNAs, gma-miR5368, novel_miR_377, novel_miR_111, and novel_miR_251, were predicted to regulate the expression of four differential expressed genes involved in the monoterpenoid and sesquiterpenoid synthesis. miRNAs families miR396, miR5185, and miR9408 were predicted to play diverse regulatory roles in monoterpenoid and sesquiterpenoid synthesis during the leaf development of C. burmannii. Our findings provide new insights into the regulatory mechanisms underlying the biosynthesis of terpenoids in C. burmannii.