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.

Author Contributions

Conceptualization, C.H. and Q.Z.; formal analysis, C.H., B.H., P.X., Y.W. and D.L.; data curation, C.H., H.L., Y.C. and D.L.; writing—original draft preparation, C.H.; writing—review and editing, C.H., B.H. and Q.Z.; project administration, C.H.; funding acquisition, B.H. and P.X. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

The datasets presented in this study are available in online repositories. The small RNA sequencing data and full-length transcriptomic data were deposited in the public SRA database with the project numbers PRJNA892262 and PRJNA892157, respectively. The metabolome data were attached to this article.

Acknowledgments

We thank Guan-Shen Liu (Biomarker Technologies, Beijing, China) for providing us with valuable technical and analytical assistance.

Conflicts of Interest

The authors declare no conflict of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures

Figure 1. Detection of differentially expressed metabolites (DEMs) and their co-expressed genes across four developmental stages (CBTS1–4) of C. burmannii leaves. (A) PCA of 12 metabolite samples across the four developmental stages. (B) Venn diagram showing the overlapping DEMs between CBTS1 and CBTS3, between CBTS2 and CBTS4, and between CBTS1 and CBTS4. (C) Among the 67 overlapping DEMs mentioned above, the content of 11 most abundant terpenoids (six monoterpenoids and five sesquiterpenoids) significantly increased from CBTS1 to CBTS4. Symbols * and *** indicate the p values < 0.05 and <0.01, respectively using Student’s t-test analysis. (D) Hierarchical clustering of co-expressed genes shows the co-expression modules. Each branch in the phylogenetic tree corresponds to an individual gene, and the interconnected genes are grouped into twelve modules. (E) Heatmap based on gene expression patterns shows the module–metabolite relationships. The color bar shows the correlation between the twelve modules and six monoterpenes and five sesquiterpenes; blocks were colored in red and blue indicate positive and negative correlations. Numbers above represent correlation index between gene expression and abundance of terpenoids, and numbers below in brackets represent highly significant (Fisher’s exact test; p < 0.01), significant (p < 0.05), and non-significant (p ≥ 0.05) correlations. (F) Numbers of GO annotated genes in different models with the p-value < 0.05. The blocks in different color are corresponding to the modules detected in Section E. (G) Numbers of KEGG annotated genes in different models with the p-value < 0.05. The blocks in different color are corresponding to the modules detected in Section E.

Figure 1. Detection of differentially expressed metabolites (DEMs) and their co-expressed genes across four developmental stages (CBTS1–4) of C. burmannii leaves. (A) PCA of 12 metabolite samples across the four developmental stages. (B) Venn diagram showing the overlapping DEMs between CBTS1 and CBTS3, between CBTS2 and CBTS4, and between CBTS1 and CBTS4. (C) Among the 67 overlapping DEMs mentioned above, the content of 11 most abundant terpenoids (six monoterpenoids and five sesquiterpenoids) significantly increased from CBTS1 to CBTS4. Symbols * and *** indicate the p values < 0.05 and <0.01, respectively using Student’s t-test analysis. (D) Hierarchical clustering of co-expressed genes shows the co-expression modules. Each branch in the phylogenetic tree corresponds to an individual gene, and the interconnected genes are grouped into twelve modules. (E) Heatmap based on gene expression patterns shows the module–metabolite relationships. The color bar shows the correlation between the twelve modules and six monoterpenes and five sesquiterpenes; blocks were colored in red and blue indicate positive and negative correlations. Numbers above represent correlation index between gene expression and abundance of terpenoids, and numbers below in brackets represent highly significant (Fisher’s exact test; p < 0.01), significant (p < 0.05), and non-significant (p ≥ 0.05) correlations. (F) Numbers of GO annotated genes in different models with the p-value < 0.05. The blocks in different color are corresponding to the modules detected in Section E. (G) Numbers of KEGG annotated genes in different models with the p-value < 0.05. The blocks in different color are corresponding to the modules detected in Section E.

Figure 2. Structural variation in miRNAs and annotation of their target genes detected in C. burmannii leaves. (A) Length distribution of novel miRNAs detected via sRNA sequencing. (B) Numbers of miRNAs assigned to different families listed in descending order. (C) Functional annotation of 1997 target genes based on gene ontology (GO) categories. (D) Functional annotation of the target genes based on the Kyoto Encyclopedia of Genes and Genomes (KEGG) categories.

Figure 2. Structural variation in miRNAs and annotation of their target genes detected in C. burmannii leaves. (A) Length distribution of novel miRNAs detected via sRNA sequencing. (B) Numbers of miRNAs assigned to different families listed in descending order. (C) Functional annotation of 1997 target genes based on gene ontology (GO) categories. (D) Functional annotation of the target genes based on the Kyoto Encyclopedia of Genes and Genomes (KEGG) categories.

Figure 3. Regulatory network of miRNAs and target genes involved in terpenoid synthesis. The large circles in different colors indicate the various miRNAs detected, and the small circles (in orange and dark red) indicate the predicted target genes. The small dark red circles indicate the seven target genes involved in the synthesis of monoterpenoids and sesquiterpenoids, according to the results in Figure 3.

Figure 4. DEmiRNAs and co-expression analysis with terpenoids. (A) A volcano plot showing the number of DEmiRNAs between CBTS1 and CBTS4 in red (up-regulated genes) and green dots (down-regulated genes) and black dots represent the remaining non-differentially expressed miRNAs. (B) A Venn diagram showing the overlap in DEmiRNAs. (C) Correlation coefficient of miRNAs and six monoterpenoids (borneol, bornyl acetate, β-phellandrene, D-limonene, β-pinene, and eucalyptol) and five sesquiterpenoids (nerolidol, γ-muulrolene, β-elemene, aromandenrene, and copaene) evaluated by Pearson correlation analysis (p-value < 0.05) across four developmental stages (CBTS1-S4) of C. burmannii leaves.

Figure 4. DEmiRNAs and co-expression analysis with terpenoids. (A) A volcano plot showing the number of DEmiRNAs between CBTS1 and CBTS4 in red (up-regulated genes) and green dots (down-regulated genes) and black dots represent the remaining non-differentially expressed miRNAs. (B) A Venn diagram showing the overlap in DEmiRNAs. (C) Correlation coefficient of miRNAs and six monoterpenoids (borneol, bornyl acetate, β-phellandrene, D-limonene, β-pinene, and eucalyptol) and five sesquiterpenoids (nerolidol, γ-muulrolene, β-elemene, aromandenrene, and copaene) evaluated by Pearson correlation analysis (p-value < 0.05) across four developmental stages (CBTS1-S4) of C. burmannii leaves.

Figure 5. DEmiRNAs and DEGs in the MVA and MEP pathways in C. burmannii. Heatmaps showing the expression levels of the DEMiRNAs and DEGs involved in the MVA and MEP pathways across four leaf developmental stages. The symbol “⊥” indicates the regulation of gene expression by the miRNAs detected in the present study based on the target gene prediction via TargetFinder.

Figure 5. DEmiRNAs and DEGs in the MVA and MEP pathways in C. burmannii. Heatmaps showing the expression levels of the DEMiRNAs and DEGs involved in the MVA and MEP pathways across four leaf developmental stages. The symbol “⊥” indicates the regulation of gene expression by the miRNAs detected in the present study based on the target gene prediction via TargetFinder.

Supplementary Materials

The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/f14030555/s1; Table S1: Information on the differentially expressed metabolites detected in the present study; Table S2: Numbers of DEMs; Table S3: Information of novel microRNA sequences, structures, and predicted hairpin energy; Table S4: Table of miRNAs and their assigned families; Table S5: MiRNA sequences and their target genes detected in the present study; Table S6: Numbers of DEmiRNAs.

References

1. Shan, T.; Wu, C.; Shahid, H.; Zhang, C.; Wang, J.; Ding, P.; Mao, Z. Powdery-fruit disease of Cinnamomum burmannii and its influence on fruit essential oil. Int. J. Agric. Biol.; 2020; 24, pp. 1077-1083.

2. Al-Dhubiab, B.E. Pharmaceutical applications and phytochemical profile of Cinnamomum burmannii. Pharmac. Rev.; 2012; 6, pp. 125-131. [DOI: https://dx.doi.org/10.4103/0973-7847.99946]

3. Ma, R.; Su, P.; Guo, J.; Jin, B.; Ma, Q.; Zhang, H.; Chen, L.; Mao, L.; Tian, M.; Lai, C. Bornyl diphosphate synthase from Cinnamomum burmanni and its application for (+)-borneol biosynthesis in yeast. Front. Bioeng. Biotechnol.; 2021; 9, 631863. [DOI: https://dx.doi.org/10.3389/fbioe.2021.631863] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/33644023]

4. Yang, Z.; An, W.; Liu, S.; Huang, Y.; Xie, C.; Huang, S.; Zheng, X. Mining of candidate genes involved in the biosynthesis of dextrorotatory borneol in Cinnamomum burmannii by transcriptomic analysis on three chemotypes. Peer J.; 2020; 8, e9311. [DOI: https://dx.doi.org/10.7717/peerj.9311] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32566406]

5. Chen, L.; Jianyu, S.; Li, L.; Li, B.; Li, W. A new source of natural D-borneol and its characteristic. J. Med. Plants Res.; 2011; 5, pp. 3440-3447.

6. Ding, J.; Yu, X.; Ding, Z.; Cheng, B.; Yi, Y.; Yu, W.; Hayashi, N.; Komae, H. Essential oils of some Lauraceae species from the southwestern parts of China. J. Essent. Oil Res.; 1994; 6, pp. 577-585. [DOI: https://dx.doi.org/10.1080/10412905.1994.9699349]

7. Mostafa, S.; Wang, Y.; Zeng, W.; Jin, B. Floral scents and fruit aromas: Functions, compositions, biosynthesis, and regulation. Front. Plant. Sci.; 2022; 13, 860157. [DOI: https://dx.doi.org/10.3389/fpls.2022.860157]

8. Ninkuu, V.; Zhang, L.; Yan, J.; Fu, Z.; Yang, T.; Zeng, H. Biochemistry of terpenes and recent advances in plant protection. Int. J. Mol. Sci.; 2021; 22, 5710. [DOI: https://dx.doi.org/10.3390/ijms22115710]

9. Abbas, F.; Ke, Y.; Yu, R.; Yue, Y.; Amanullah, S.; Jahangir, M.M.; Fan, Y. Volatile terpenoids: Multiple functions, biosynthesis, modulation and manipulation by genetic engineering. Planta; 2017; 246, pp. 803-816. [DOI: https://dx.doi.org/10.1007/s00425-017-2749-x]

10. Shan, B.; Cai, Y.-Z.; Brooks, J.D.; Corke, H. Antibacterial properties and major bioactive components of cinnamon stick (Cinnamomum burmannii): Activity against foodborne pathogenic bacteria. J. Agr. Food Chem.; 2007; 55, pp. 5484-5490. [DOI: https://dx.doi.org/10.1021/jf070424d]

11. Vranová, E.; Coman, D.; Gruissem, W. Network analysis of the MVA and MEP pathways for isoprenoid synthesis. Annu. Rev. Plant Biol.; 2013; 64, pp. 665-700. [DOI: https://dx.doi.org/10.1146/annurev-arplant-050312-120116]

12. Tholl, D. Biosynthesis and biological functions of terpenoids in plants. Adv. Biochem. Eng./Biotechnol.; 2015; 148, pp. 63-106.

13. Qiu, F.; Wang, X.; Zheng, Y.; Wang, H.; Liu, X.; Su, X. Full-length transcriptome sequencing and different chemotype expression profile analysis of genes related to monoterpenoid biosynthesis in Cinnamomum porrectum. Int. J. Mol. Sci.; 2019; 20, 6230. [DOI: https://dx.doi.org/10.3390/ijms20246230]

14. Chen, C.; Zheng, Y.; Zhong, Y.; Wu, Y.; Li, Z.; Xu, L.-A.; Xu, M. Transcriptome analysis and identification of genes related to terpenoid biosynthesis in Cinnamomum camphora. BMC Genom.; 2018; 19, 550. [DOI: https://dx.doi.org/10.1186/s12864-018-4941-1]

15. Zhu, C.; Zhang, F.; Chen, S.; Wang, K.; Xiang, G.; Liang, X.; An, J.; Li, K.; Liu, L. Proteomics analysis and identification of proteins related to isoprenoid biosynthesis in Cinnamomum camphora (L.) Presl. Forests; 2022; 13, 1487. [DOI: https://dx.doi.org/10.3390/f13091487]

16. Li, F.; Huang, S.; Mei, Y.; Wu, B.; Hou, Z.; Zhan, P.; Hou, Z.; Huang, W.; Zhao, J.; Wang, J. et al. Genome assembly provided new insights into the Cinnamomum burmannii evolution and D-borneol biosynthesis differences between chemotypes. Ind. Crops Prod.; 2022; 186, 115181. [DOI: https://dx.doi.org/10.1016/j.indcrop.2022.115181]

17. Ma, Q.; Ma, R.; Su, P.; Jin, B.; Guo, J.; Tang, J.; Chen, T.; Zeng, W.; Lai, C.; Ling, F. Elucidation of the essential oil biosynthetic pathways in Cinnamomum burmannii through identification of six terpene synthases. Plant Sci.; 2022; 317, 111203. [DOI: https://dx.doi.org/10.1016/j.plantsci.2022.111203]

18. Axtell, M. Classification and comparison of small RNAs from plants. Annu. Rev. Plant Biol.; 2013; 64, pp. 137-159. [DOI: https://dx.doi.org/10.1146/annurev-arplant-050312-120043]

19. Wu, G. Plant microRNAs and development. J. Genet. Genom.; 2013; 40, pp. 217-230. [DOI: https://dx.doi.org/10.1016/j.jgg.2013.04.002]

20. Axtell, M.J.; Bowman, J.L. Evolution of plant microRNAs and their targets. Trends Plant Sci.; 2008; 13, pp. 343-349. [DOI: https://dx.doi.org/10.1016/j.tplants.2008.03.009]

21. Xiao, M.; Li, J.; Li, W.; Wang, Y.; Wu, F.; Xi, Y.; Zhang, L.; Ding, C.; Luo, H.; Li, Y. MicroRNAs activate gene transcription epigenetically as an enhancer trigger. RNA Biol.; 2017; 14, pp. 1326-1334. [DOI: https://dx.doi.org/10.1080/15476286.2015.1112487] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/26853707]

22. Shen, T.; Xu, M.; Qi, H.; Feng, Y.; Yang, Z.; Xu, M. Uncovering miRNA-mRNA regulatory modules in developing xylem of Pinus massoniana via small RNA and degradome sequencing. Int. J. Mol. Sci.; 2021; 22, 10154. [DOI: https://dx.doi.org/10.3390/ijms221810154]

23. Li, A.; Mao, L. Evolution of plant microRNA gene families. Cell Res.; 2007; 17, pp. 212-218. [DOI: https://dx.doi.org/10.1038/sj.cr.7310113] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/17130846]

24. Narnoliya, L.K.; Kaushal, G.; Singh, S.P. Long noncoding RNA s and miRNAs regulating terpene and tartaric acid biosynthesis in rose-scented geranium. FEBS Lett.; 2019; 593, pp. 2235-2249. [DOI: https://dx.doi.org/10.1002/1873-3468.13493] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/31210363]

25. Ye, J.; Zhang, X.; Tan, J.; Xu, F.; Cheng, S.; Chen, Z.; Zhang, W.; Liao, Y. Global identification of Ginkgo biloba microRNAs and insight into their role in metabolism regulatory network of terpene trilactones by high-throughput sequencing and degradome analysis. Ind. Crops Prod.; 2020; 148, 112289. [DOI: https://dx.doi.org/10.1016/j.indcrop.2020.112289]

26. Chen, C.; Zhong, Y.; Yu, F.; Xu, M. Deep sequencing identifies miRNAs and their target genes involved in the biosynthesis of terpenoids in Cinnamomum camphora. Ind. Crops Prod.; 2020; 145, 111853. [DOI: https://dx.doi.org/10.1016/j.indcrop.2019.111853]

27. Lê, S.; Josse, J.; Husson, F. FactoMineR: An R package for multivariate analysis. J. Stat. Softw.; 2008; 25, pp. 1-18. [DOI: https://dx.doi.org/10.18637/jss.v025.i01]

28. R Core Team. R: A Language and Environment for Statistical Computing, Version 3.2.0; R Foundation for Statistical Computing: Vienna, Austria, 2018; Available online: https://www.R-project.org (accessed on 7 September 2022).

29. Hou, C.; Zhang, Q.; Xie, P.; Lian, H.; Wang, Y.; Liang, D.; Cai, Y.; He, B. Full-length transcriptome sequencing reveals the molecular mechanism of monoterpene and sesquiterpene biosynthesis in Cinnamomum burmannii. Front. Genet.; 2023; 13, 1087495. [DOI: https://dx.doi.org/10.3389/fgene.2022.1087495]

30. Langfelder, P.; Horvath, S. WGCNA: An R package for weighted correlation network analysis. BMC Bioinform.; 2008; 9, 559. [DOI: https://dx.doi.org/10.1186/1471-2105-9-559]

31. Langfelder, P.; Zhang, B.; Horvath, S. Defining clusters from a hierarchical cluster tree: The Dynamic Tree Cut package for R. Bioinformatics; 2008; 24, pp. 719-720. [DOI: https://dx.doi.org/10.1093/bioinformatics/btm563]

32. Johnson, M.; Zaretskaya, I.; Raytselis, Y.; Merezhuk, Y.; McGinnis, S.; Madden, T. NCBI BLAST: A better web interface. Nucl. Acids Res.; 2008; 36, pp. 5-9. [DOI: https://dx.doi.org/10.1093/nar/gkn201]

33. Langmead, B.; Trapnell, C.; Pop, M.; Salzberg, S.L. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol.; 2009; 10, R25. [DOI: https://dx.doi.org/10.1186/gb-2009-10-3-r25]

34. Griffiths-Jones, S. miRBase: The microRNA sequence database. MicroRNA Protoc.; 2006; 342, pp. 129-138.

35. Friedländer, M.R.; Mackowiak, S.D.; Li, N.; Chen, W.; Rajewsky, N. miRDeep2 accurately identifies known and hundreds of novel microRNA genes in seven animal clades. Nucl. Acids Res.; 2012; 40, pp. 37-52. [DOI: https://dx.doi.org/10.1093/nar/gkr688]

36. Rehmsmeier, M.; Steffen, P.; Höchsmann, M.; Giegerich, R. Fast and effective prediction of microRNA/target duplexes. RNA; 2004; 10, pp. 1507-1517. [DOI: https://dx.doi.org/10.1261/rna.5248604]

37. Allen, E.; Xie, Z.; Gustafson, A.M.; Carrington, J.C. MicroRNA-directed phasing during trans-acting siRNA biogenesis in plants. Cell; 2005; 121, pp. 207-221. [DOI: https://dx.doi.org/10.1016/j.cell.2005.04.004]

38. Li, B.; Ruotti, V.; Stewart, R.M.; Thomson, J.A.; Dewey, C.N. RNA-Seq gene expression estimation with read mapping uncertainty. Bioinformatics; 2010; 26, pp. 493-500. [DOI: https://dx.doi.org/10.1093/bioinformatics/btp692]

39. Love, M.; Anders, S.; Huber, W. Differential analysis of count data–the DESeq2 package. Genome Biol.; 2014; 15, pp. 10-1186.

40. Love, M.I.; Huber, W.; Anders, S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. Evol.; 2014; 15, 550. [DOI: https://dx.doi.org/10.1186/s13059-014-0550-8]

41. Plata-Rueda, A.; Campos, J.M.; da Silva Rolim, G.; Martínez, L.C.; Dos Santos, M.H.; Fernandes, F.L.; Serrão, J.E.; Zanuncio, J.C. Terpenoid constituents of cinnamon and clove essential oils cause toxic effects and behavior repellency response on granary weevil, Sitophilus granarius. Ecotoxicol. Environ. Saf.; 2018; 156, pp. 263-270. [DOI: https://dx.doi.org/10.1016/j.ecoenv.2018.03.033]

42. Aldrich, J. Chemical ecology of the Heteroptera. Annu. Rev. Entomol.; 1988; 33, pp. 211-238. [DOI: https://dx.doi.org/10.1146/annurev.en.33.010188.001235]

43. Douglas, H.; Co, J.; Jones, T.; Conner, W. Heteropteran chemical repellents identified in the citrus odor of a seabird (crested auklet: Aethia cristatella): Evolutionary convergence in chemical ecology. Sci. Nat.; 2001; 88, pp. 330-332. [DOI: https://dx.doi.org/10.1007/s001140100236] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/11572013]

44. Wu, B.; Ruan, C.; Shah, A.H.; Li, D.; Li, H.; Ding, J.; Li, J.; Du, W. Identification of miRNA–mRNA regulatory modules involved in lipid metabolism and seed development in a woody oil tree (Camellia oleifera). Cells; 2022; 11, 71. [DOI: https://dx.doi.org/10.3390/cells11010071]

45. Liu, W.; Cheng, C.; Chen, F.; Ni, S.; Lin, Y.; Lai, Z. High-throughput sequencing of small RNAs revealed the diversified cold-responsive pathways during cold stress in the wild banana (Musa itinerans). BMC Plant Biol.; 2018; 18, 308. [DOI: https://dx.doi.org/10.1186/s12870-018-1483-2] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/30486778]

46. Tian, Z.; Luo, Q.; Zuo, Z. Seasonal emission of monoterpenes from four chemotypes of Cinnamomum camphora. Ind. Crops Prod.; 2021; 163, 113327. [DOI: https://dx.doi.org/10.1016/j.indcrop.2021.113327]

47. Ni, Z.; Han, X.; Chen, C.; Zhong, Y.; Xu, M.; Xu, L.A.; Yu, F. Integrating GC-MS and ssRNA-Seq analysis to identify long non-coding RNAs related to terpenoid biosynthesis in Cinnamomum camphora. Ind. Crops Prod.; 2021; 171, 113875. [DOI: https://dx.doi.org/10.1016/j.indcrop.2021.113875]

48. Yang, Z.; Xie, C.; Huang, Y.; An, W.; Liu, S.; Huang, S.; Zheng, X. Metabolism and transcriptome profiling provides insight into the genes and transcription factors involved in monoterpene biosynthesis of borneol chemotype of Cinnamomum camphora induced by mechanical damage. PeerJ; 2021; 9, e11465. [DOI: https://dx.doi.org/10.7717/peerj.11465]

49. Baghban Kohnehrouz, B.; Bastami, M.; Nayeri, S. In silico identification of novel microRNAs and targets using EST analysis in Allium cepa L. Interdiscip. Sci.-Comput. Life Sci.; 2018; 10, pp. 771-780. [DOI: https://dx.doi.org/10.1007/s12539-017-0240-9]

50. Bertolini, E.; Verelst, W.; Horner, D.S.; Gianfranceschi, L.; Piccolo, V.; Inzé, D.; Pè, M.E.; Mica, E. Addressing the role of microRNAs in reprogramming leaf growth during drought stress in Brachypodium distachyon. Mol. Plant; 2013; 6, pp. 423-443. [DOI: https://dx.doi.org/10.1093/mp/sss160]

51. Xie, X.; Li, X.; Tian, Y.; Su, M.; Zhang, J.; Han, X.; Cui, Y.; Bian, S. Identification and characterization of microRNAs and their targets from expression sequence tags of Ribes nigrum. Can. J. Plant Sci.; 2016; 96, pp. 995-1001. [DOI: https://dx.doi.org/10.1139/cjps-2016-0013]

52. Dan, M.; Huang, M.; Liao, F.; Qin, R.; Liang, X.; Zhang, E.; Huang, M.; Huang, Z.; He, Q. Identification of ethylene responsive miRNAs and their targets from newly harvested banana fruits using high-throughput sequencing. J. Agric. Food Chem.; 2018; 66, pp. 10628-10639. [DOI: https://dx.doi.org/10.1021/acs.jafc.8b01844]

53. Kong, X.; He, M.; Guo, Z.; Zhou, Y.; Chen, Z.; Qu, H.; Zhu, H. MicroRNA regulation of fruit development, quality formation and stress response. Fruit Res.; 2021; 1, 5. [DOI: https://dx.doi.org/10.48130/FruRes-2021-0005]

54. Khan, S.; Ali, A.; Saifi, M.; Saxena, P.; Ahlawat, S.; Abdin, M.Z. Identification and the potential involvement of miRNAs in the regulation of artemisinin biosynthesis in A. annua. Sci. Rep.; 2020; 10, 13614. [DOI: https://dx.doi.org/10.1038/s41598-020-69707-3]

55. Jeyaraj, A.; Zhang, X.; Hou, Y.; Shangguan, M.; Gajjeraman, P.; Li, Y.; Wei, C. Genome-wide identification of conserved and novel microRNAs in one bud and two tender leaves of tea plant (Camellia sinensis) by small RNA sequencing, microarray-based hybridization and genome survey scaffold sequences. BMC Plant Biol.; 2017; 17, 212. [DOI: https://dx.doi.org/10.1186/s12870-017-1169-1]

56. Zhao, X.; Yang, G.; Liu, X.; Yu, Z.; Peng, S. Integrated analysis of seed microRNA and mRNA transcriptome reveals important functional genes and microRNA-targets in the process of walnut (Juglans regia) seed oil accumulation. Int. J. Mol. Sci.; 2020; 21, 9093. [DOI: https://dx.doi.org/10.3390/ijms21239093]

57. Parreira, J.R.; Cappuccio, M.; Balestrazzi, A.; Fevereiro, P.; Araújo, S.d.S. MicroRNAs expression dynamics reveal post-transcriptional mechanisms regulating seed development in Phaseolus vulgaris L. Hortic. Res.; 2021; 8, 18. [DOI: https://dx.doi.org/10.1038/s41438-020-00448-0]

58. Dai, Z.; Tan, J.; Zhou, C.; Yang, X.; Yang, F.; Zhang, S.; Sun, S.; Miao, X.; Shi, Z. The OsmiR396–Os GRF 8–OsF3H-flavonoid pathway mediates resistance to the brown planthopper in rice (Oryza sativa). Plant Biotech. J.; 2019; 17, pp. 1657-1669. [DOI: https://dx.doi.org/10.1111/pbi.13091]

59. Miao, C.; Wang, D.; He, R.; Liu, S.; Zhu, J.K. Mutations in MIR 396e and MIR 396f increase grain size and modulate shoot architecture in rice. Plant Biotech. J.; 2020; 18, pp. 491-501. [DOI: https://dx.doi.org/10.1111/pbi.13214]

60. Casadevall, R.; Rodriguez, R.E.; Debernardi, J.M.; Palatnik, J.F.; Casati, P. Repression of growth regulating factors by the microRNA396 inhibits cell proliferation by UV-B radiation in Arabidopsis leaves. Plant Cell; 2013; 25, pp. 3570-3583. [DOI: https://dx.doi.org/10.1105/tpc.113.117473]

61. Rodriguez, R.E.; Mecchia, M.A.; Debernardi, J.M.; Schommer, C.; Weigel, D.; Palatnik, J.F. Control of cell proliferation in Arabidopsis thaliana by microRNA miR396. Development; 2010; 137, pp. 103-112. [DOI: https://dx.doi.org/10.1242/dev.043067]