Regulation of terpenoid biosynthesis by miRNA in Persicaria minor induced by Fusarium oxysporum

Background Persicaria minor (kesum) is an herbaceous plant with a high level of secondary metabolite compounds, particularly terpenoids. These terpenoid compounds have well-established roles in the pharmaceutical and food industries. Although the terpenoids of P. minor have been studied thoroughly, the involvement of microRNA (miRNA) in terpenoid regulation remains poorly understood and needs to be explored. In this study, P. minor plants were inoculated with the pathogenic fungus Fusarium oxysporum for terpenoid induction. Result SPME GC-MS analysis showed the highest terpenoid accumulation on the 6th day post-inoculation (dpi) compared to the other treatment time points (0 dpi, 3 dpi, and 9 dpi). Among the increased terpenoid compounds, α-cedrene, valencene and β-bisabolene were prominent. P. minor inoculated for 6 days was selected for miRNA library construction using next generation sequencing. Differential gene expression analysis showed that 58 miRNAs belonging to 30 families had significantly altered regulation. Among these 58 differentially expressed genes (DEGs), 33 miRNAs were upregulated, whereas 25 miRNAs were downregulated. Two putative novel pre-miRNAs were identified and validated through reverse transcriptase PCR. Prediction of target transcripts potentially involved in the mevalonate pathway (MVA) was carried out by psRobot software, resulting in four miRNAs: pmi-miR530, pmi-miR6173, pmi-miR6300 and a novel miRNA, pmi-Nov_13. In addition, two miRNAs, miR396a and miR398f/g, were predicted to have their target transcripts in the non-mevalonate pathway (MEP). In addition, a novel miRNA, pmi-Nov_12, was identified to have a target gene involved in green leaf volatile (GLV) biosynthesis. RT-qPCR analysis showed that pmi-miR6173, pmi-miR6300 and pmi-nov_13 were downregulated, while miR396a and miR398f/g were upregulated. Pmi-miR530 showed upregulation at 9 dpi, and dynamic expression was observed for pmi-nov_12. Pmi-6300 and pmi-miR396a cleavage sites were detected through degradome sequence analysis. Furthermore, the relationship between miRNA metabolites and mRNA metabolites was validated using correlation analysis. Conclusion Our findings suggest that six studied miRNAs post-transcriptionally regulate terpenoid biosynthesis in P. minor. This regulatory behaviour of miRNAs has potential as a genetic tool to regulate terpenoid biosynthesis in P. minor. Electronic supplementary material The online version of this article (10.1186/s12864-019-5954-0) contains supplementary material, which is available to authorized users.


Background
Humans exploit the secondary metabolites from medicinal plants in the form of flavouring agents, perfumes, insecticides, dyes and drugs, among other products. Previous studies have shown that more than 100,000 phytochemicals have been isolated from different plant sources [1]. According to the British Nutrition Foundation, secondary metabolites (SMs) are divided into four major classes: terpenoids (volatile compounds, carotenoids and glycosides), phenolic compounds (phenolic acids, tannins and flavonoids), nitrogen-containing compounds (cyanogenic glucosides and alkaloids) and sulfur-containing compounds (thionine, defensin and lectin) [2]. The terpenoid classes consist of volatile and non-volatile compounds. Though these compounds exist in complex structures, they all are made up of isoprene units (C5). Most volatile terpenoids are monoterpenes (C10) and sesquiterpenes (C15), while non-volatile terpenoids include diterpenes (C20), triterpenes (C30) and tetraterpenes (C40) [3]. Terpenoids are widely used in pharmaceuticals, cosmetic fragrances and the food industry [4]. To date, more than 40,000 known terpenoid compounds have been identified, mostly from plants [5].
Persicaria minor (locally known as kesum) is an aromatic plant that is widespread in Southeast Asia, especially in Malaysia, Thailand, Laos, Indonesia and Vietnam. This plant has been used in dishes and medical purposes for many years. The Malaysian government enlisted P. minor in the National Agro Food Policy to ensure its adequate supply and as a platform to strengthen the agriculturalbased economy [6]. P. minor has also gained attention due to its abundance of secondary metabolites, especially flavonoids and terpenoids, which are used in the food, fragrance and pharmaceutical industries [7]. P. minor leaves produce the highest content of terpenoids compared to other organs [8]. These compounds play important roles in attracting pollinators, plant defence, and interactions with unfavourable environments. Two different pathways of terpenoid biosynthesis have been identified: the mevalonate (MVA) and non-mevalonate, or methylerythritol 4-phosphate (MEP), pathways [9]. These pathways operate in the cytoplasm and plastid, respectively [9].
Past achievements in unlocking the genes and enzymes involved in each of these pathways has opened up new possibilities for the assessment of terpenoid biosynthesis [10]. At the molecular level, molecules such as DNA, mRNA, transcription factors, and non-coding RNA play roles in regulating the genes involved in these pathways. In addition, these molecules may interact with each other, which can affect the target gene, thus influencing the production of the secondary metabolite [11]. However, the level of secondary metabolites in plants is relatively low because their production depends on plant species, environmental factors, and nutritional sources [12,13]. In addition, the type of secondary metabolites produced by the plant is stimulus-dependent. For example, in the majority of cases, terpenoid compounds are produced in high abundance under biotic stress compared to abiotic [14]. Hence, the elicitation technique was introduced to enhance secondary metabolite production in plants using a biotic elicitor [13,15]. Elicitation using fungi, especially F. oxysporum, has been reported to enhance terpenoid contents in plants [16][17][18][19].
To understand the changes at the molecular level, the relationship between miRNA and mRNA needs to be explored, because the interaction between these two RNA species may play a significant role in plant secondary metabolite production. miRNAs are a group of noncoding RNAs with small sizes (~20 bp) that act as a gene regulators by negatively regulating mRNAs via cleavage or translational inhibition [11,20]. The first miRNA, lin-4 miRNA, was discovered in C. elegans through forward genetic analysis. Later, several methods were developed to discover these tiny RNAs, such as cloning and computational methods [20]. However, these methods had their own limitations. Due to the ability of several miRNAs to target single mRNAs, using forward genetic techniques for miRNA discovery is very difficult unless the miRNA of interest has only one target [11,20]. The cloning method is also limited to highly expressed miRNAs [20]. Computational approaches could be useful for miRNA discovery but have limitations in detecting novel miRNAs in other species, since these approaches are based on homology searches. In addition, computational approaches still require experimental validation [20]. Next generation sequencing (NGS), a more recent approach, has also been introduced for miRNA study, and this method has been found to be more advanced and efficient. This approach can detect low copy number miRNAs at levels of one transcript per million [21,22]. To date, a total of 38,589 miRNAs from animals, plants, and viruses have been discovered and deposited in public miRNA databases [23]. These miRNAs have the ability to fine-tune many biological and physiological processes such as plant development and stress response by regulating a number of target genes in a time-sensitive manner [11].
Interestingly, miRNAs were reported to be involved in plant secondary metabolite regulation, such as terpenoid, phenolic and nitrogen-containing compound biosynthesis [1]. These interactions could hold potential for future genetic manipulation. In addition, using miRNA as a tool for genetic engineering could be highly beneficial because single miRNAs can target more than single mRNAs and vice versa. Authenticated in silico miRNA discovery and wet lab approaches for miRNA target validation are required prior to their application to plant genetic systems. Thus, in this study, the regulatory roles of miRNAs in the terpenoid biosynthesis pathway were explored in P. minor induced by F. oxysporum. Metabolite profiling was carried out for the F. oxysporum-inoculated plants to determine the changes in the accumulation of terpenoid compounds. In addition, small RNA libraries were constructed, which led to the characterization and analysis of miRNAs involved in the terpenoid biosynthesis pathway.

Results and discussion
Testing of F. oxysporum inoculation and terpenoid profiling in P. minor by SPME-GCMS To the best of our knowledge, no work on F. oxysporum inoculation to induce terpenoid contents has been reported to date in P. minor. Therefore, we carried out an inoculation test to determine the elicitation ability of this fungus to increase terpenoid content in P. minor. Comparison between mock-inoculated (C) and Fusariumtreated (F) plants showed different leaf morphologies over the time point of the treatment (Fig. 1). From the beginning of the treatment (0 dpi) to 3 dpi, there were no morphological changes observed in either of the C or F samples. At 6 dpi, no change was observed in the C samples, while the F samples showed wilting symptoms beginning at the tip of the leaf. At 9 dpi, the wilting spread to the middle of the leaves in the F samples. At 12 dpi, the wilting spread to the whole leaf, which led to the death of the plant. In comparison, no changes were observed in the C samples until the end point, 12 dpi. A previous report mentioned an accumulation of terpenoids when plants were attacked by pathogens [24]. Based on Fig. 1, the first symptom appeared at 6 dpi. Theoretically, the test suggested the highest terpenoid content at 6 dpi in P. minor inoculated with F. oxysporum. However, the SPME-GCMS result was carried out to reveal the terpenoid content at each time point. The initial SPME-GCMS data are provided in Additional file 1. The heatmap generated (Fig. 2) indicated a significant increase in three volatile compounds: α-cedrene (1H-3a,7-methanoazulene, 2,3,4,7,8,8a-hexahydro-3,6,8,8tetramethyl-, [3R-(3.α.,3a.beta.,7.beta.,8a.α.)]-), valencene (naphthalene, 1,2,3,5,6,7,8,8a-octahydro-1,8a-dimethyl-7-(1methylethenyl)-, [1S-(1.α.,7.α.,8a.α.)]-) and β-bisabolene on 6 dpi. Hence, samples at the 6 dpi time point were selected for small RNA library construction.
High-throughput small RNA sequencing analysis Small RNA sequencing for the C and F libraries generated approximately 12,409,685 and 42,985,084 reads, respectively. The average numbers and proportions of the different categories are presented in Table 1. Trimming of adaptor index sequences was carried out, and lowquality reads were removed to produce clean reads. Reads with lengths between 18 nt and 30 nt were used for annotation, while the rest were discarded. In total, approximately 7,724,932 and 25,726,524 unique reads were produced for the C and F libraries, respectively. The length distribution of the small RNAs in the C and F libraries is shown in Fig. 3. The most abundant small RNAs were 22 nt in length, followed by 20 nt, which belonged to the F library. Previous studies have reported that more than 60% of all plant miRNAs are 21 nt in length [25]. The clean reads from each P. minor library were annotated according to the miRBase version 21 Fig. 1 Morphological changes in P. minor leaves inoculated with F. oxysporum in a time series (0 dpi, 3 dpi, 6 dpi, 9 dpi and 12 dpi)

Differential expression of miRNAs in C and F libraries
The overall expression of miRNAs in both libraries is shown in the volcano plot (Fig. 4). The plot showed significant changes in the regulation of 58 miRNAs. Among these significantly regulated miRNAs, 31 were upregulated and 27 were downregulated. Additionally, a heatmap was generated to compare the expression of 58 significantly regulated miRNAs in the C and F libraries (Fig. 5). Further details about the significantly regulated miRNAs are documented in Table 2.    Common and specific miRNAs A Venn diagram was generated to show the common and specific miRNAs in the C and F libraries (Fig. 6). Of the 58 miRNAs responsive to F. oxysporum treatment, 42 miRNAs were common to both libraries. Two miR-NAs, pmi-miR168b and pmi-miR2111, were found to be specific to the C library, whereas 14 miRNAs were specific to the F library.

Discovery of novel miRNAs
To discover novel miRNA sequences in P. minor, all unannotated small RNAs were searched against the P. minor transcriptome data under accession number SRX669305 [26]. After searching for hairpin structures (Fig. 7) and performing MFEI calculations (Table 3), two unique sequences were identified as putative novel miR-NAs in P. minor. These putative novel miRNAs were named pmi-nov_miR12 and pmi-nov_miR13. The secondary structures of both miRNAs were validated through reverse transcriptase PCR (Fig. 8). The sizes of the pmi-nov_miR12 precursor (pre-nov_12) and pmi-nov_miR13 precursor (pre-nov_13) are~80 bp and~75 bp, respectively. The sizes of the validated precursors were slightly smaller than the predicted sizes. These smaller sizes might be due to the primers, which were designed to avoid the bulge structure in the hairpin structures. The bulge region in the hairpin structure could have reduced the efficiency of primer binding to the DNA template [27].

Prediction and functional annotation of target genes
In plants, miRNAs are involved in many biological processes by negatively regulating their target genes/transcripts [11,28]. In this study, we searched for putative target genes against the P. minor transcriptomic library software with mismatch score of 4.0 [29]. This prediction resulted in a total of 111 potential target genes in P. minor. Among them, there were 108 target genes for 58 conserved miR-NAs and 3 for two novel miRNAs ( Table 4). The majority of the miRNAs have more than one target gene. The target genes were classified according to the WEGO database into cellular component, biological process and molecular function categories (Fig. 9).

Analysis of target transcripts involved in terpenoid pathway
The terpenoid pathway consists of two distinct pathways, MVA and MEP. From the target prediction (Table 4), a total of seven miRNA targets were discovered that seemed to be directly involved in these terpenoid biosynthesis pathways. In the MVA pathway, the target involved were diphosphomevalonate decarboxylase (MVD), targeted by pmi-miR530; sesquiterpene synthase and farnesyl diphosphate synthase (FDS), targeted by pmi-miR6173; 3-hydroxy-3-methylglutaryl-coenzyme A reductase (HMGR), targeted by pmi-miR6300; and mevalonate kinase (MVK), targeted by pmi-nov_13. In the MEP or non-mevalonate pathway, the identified targets were 1-deoxy-d-xylulose-5-phosphate synthase (DXS) and 1-deoxy-d-xylulose-5-phosphate reductoisomerase (DXR), which were targeted by pmi-miR396a and pmi-miR398f/g, respectively. In addition, two more targets, peroxidase and alcohol dehydrogenase (ADH), were targeted by pmi-miR396a and pmi-nov_12, respectively. Both of these targets were included in this analysis, because peroxidase is involved in the early signalling of secondary metabolite biosynthesis, and ADH is an enzyme involved in green leaf volatile (GLV) biosynthesis. GLV compounds have been reported to have similar  physiological and functional properties to those of terpenoids [30,31].

Expression profiles of miRNAs and their targets by RT-qPCR
Real-time quantitative polymerase chain reaction (RT-qPCR) was performed to experimentally validate the expression of five conserved miRNAs (pmi-miR396a, pmi-miR398f/g, pmi-miR530, pmi-miR6173, and pmi-miR6300) and two novel miRNAs (pmi-nov_12 and pmi-nov_13) and their target transcripts. The expression profiles are shown in Fig. 10. Pmi-miR396a showed upregulation to 2.5-fold at 9 dpi. Its target, peroxidase57, was upregulated at 3 dpi and then downregulated at 6 dpi and 9 dpi. Peroxidase is an important enzyme that acts as a scavenger for reactive oxygen species (ROS) [32]. ROS are produced by plants as an early response to stress. Simultaneously, plants require a mechanism to prevent ROS from damaging the cell. Hence, high production of peroxidase during early inoculation (3 dpi) may help neutralize the excessive ROS inside the plant cell. In addition, peroxidase may also be involved in terpenoid biosynthesis. In cucumber inoculated with red spider mites (Tetranychus urticae Koch), increased production of (E, E)-α-farnesene resulted [33]. In A. thaliana, miR396 is encoded by two different gene loci, which lead to the biogenesis of two miR396 members, miR396a and miR396b, targeting the GRF transcription factor, which regulates the number of cells in the leaf [34,35]. Overexpression of miR396a and miR396b resulted in reduced leaf size [34]. In addition, recent studies have revealed that the application of a target mimic approach in miR396a contributes to plant defence against fungal pathogens. The study suggested that the low activity of miR396a induces the plant defence mechanism through the accumulation of hydrogen peroxidase (H 2 O 2 ) and callus formation [36]. In rice, a similar approach was used, leading to decreased accumulation of miR396a and increased GRF transcription factor activity. That kind of interaction activates auxin biosynthesis and ARF, thus resulting in greater yield [37]. In addition to peroxidase, another target of pmi-miR396a in P. minor was DXS. Negative correlations between these two were observed at 3 dpi. DXS is involved in the MEP pathway, which produces monoterpene and diterpene as the main products [9,10]. The descending pattern of DXS expression corresponds with the increase in pmi-miR396a, especially at 3 dpi and 9 dpi. For the target transcript peroxidase57, the target was not strongly regulated with high expression of pmi-miR396a at 3 dpi. This kind of regulation was also demonstrated for miR398 in A. thaliana [38]. In A. thaliana, miR398 targeted two types of copper superoxide dismutase (CSD), CSD1 and CSD2. The upregulation of miR398 consistently leads to the downregulation of CSD1, whereas CSD2 did not show any correlation [38]. In P. minor, the expression profile of pmi-miR398f/g showed an increasing pattern from 3 dpi to 9 dpi, and the maximum was observed at 9 dpi. The upregulation of pmi-miR398f/g led to the downregulation of the DXR target transcript. Similar to DXS, DXR is involved in the MEP pathway of terpenoid biosynthesis [10]. The miR398 family is thought to be involved in various plant stress responses by participating in the oxidative burst process [39,40]. In A. thaliana, the expression of miR398 was downregulated under salinity, oxidative stress, Pseudomonas syringae infection and phosphate deficiency. This downregulation leads to the accumulation of a target transcript, CSD, which acts as a ROS scavenger [39]. In P. minor, pmi-miR398f/g was upregulated upon treatment with F. oxysporum. The expression pattern seems to contrast with that of miR398 in A. thaliana. However, in Medicago truncatula, miR398a/b,  Adjusted minimum folding energy f Minimum folding energy index which targets copper-containing proteins involved in copper homeostasis, was upregulated under drought conditions [41].
Pmi-miR6300 showed a decreasing pattern from 3 dpi until 9 dpi, resulting in the accumulation of the target transcript, HMGR. The highest accumulation of HMGR was observed at 9 dpi, when RT-qPCR showed an upregulation of more than three-fold. HMGR is a vital enzyme involved in sesquiterpene biosynthesis in the MVA pathway [10]. A very limited amount is known about the role of the miR6300 family in plants. In barley, miR6300 was downregulated under drought treatment [42]. In addition, miR6300 has also been discovered in chickpea [43]. However, the target transcripts of miR6300 in both the mentioned studies are still unknown.
Pmi-miR530 targeted MVD, which is involved in the MVA pathway. No significant change was observed in pmi-miR539 expression until 6 dpi. On the other hand, the target transcript MVD showed a gradual increase from 3 dpi until 6 dpi. At 9 dpi, pmi-miR530 was drastically upregulated (exceeding 12-fold), while MVD was repressed. The spike in pmi-miR530 expression was supported by a previous study in chickpea. In chickpea inoculated with F. oxysporum, there was a 17-fold upregulation of miR530 compared to that in control plants. In addition, miR530 in chickpea targeted zinc knuckle protein, which is involved in regulating plant development [44].
Two target transcripts, FDP and sesquiterpene synthase, which are involved in the MVA pathway, were targeted by pmi-miR6173. RT-qPCR analysis showed a decreasing pattern in pmi-miR6173 expression. The expression patterns for both targets were similar. The expression of target transcripts increased at 3 dpi. However, the expression declined at 6 dpi, and no significant change was observed from 6 dpi to 9 dpi. A decreasing expression pattern of miR6173 was also observed in the herbs Sedum alfredii and Medicago sativa. In S. alfredii, miR6173 was downregulated under cadmium treatment. Target prediction showed that miR6173 in S. alfredii targeted a number of targets, such as calcium-binding EF-hand family protein, ATP synthase subunit α, and aspartic protease. However, further work is needed to discover the role of miR6173 in S. alfredii [45]. In M. sativa, miR6173 was among the miR-NAs downregulated under drought treatment. miR6173 in M. sativa targeted splicing factor 3A subunit 2, which plays an important role in mRNA splicing [46].
The single gene involved in GLV biosynthesis, ADH, was targeted by the novel miRNA pmi-nov_12. Our results showed that pmi-nov_12 was downregulated from 3 dpi to 6 dpi and then upregulated at 9 dpi, where a 1.2-fold increase was observed. The expression pattern of pmi-nov_12 was downregulated and upregulated at a series of time points, showing a dynamic expression pattern. This type of expression was previously reported in Populus tremula. In P. tremula treated with salinity and abscisic acid, miR398 was downregulated in the first 48 h and then showed high expression at 72 h. In P. minor, although dynamic regulation occurs in pmi-nov_12, a negative correlation can still be observed between pmi-nov_12 and ADH. ADH was reported to be involved in GLV biosynthesis by catalysing the conversion between alcohol and aldehyde via the hydroperoxide lyase pathway [47]. Another novel miRNA in P. minor, pmi-nov_ 13, was downregulated in response to F. oxysporum treatment, and its target transcript (MVK) was upregulated, especially at 6 dpi. MVK is involved in the MVA pathway of terpenoid biosynthesis. Notably, pmi-nov_12 and pmi-nov_13 are unique miRNAs in P. minor that have never been reported in another plant species   before. However, for pmi-miR6173, pmi-nov_12 and pmi-nov_13, the expression of the target genes was comparable to miRNA expression even with decreasing metabolite content at the end of the treatment (9 dpi). Hence, we suggest that these miRNAs may inhibit mRNA translation, which is another mode of miRNA action [48,49].
miRNAs as post-transcriptional regulator in P. minor In plants, miRNAs play roles in various biological processes, such as plant development, signal transduction, stress response, and secondary metabolite regulation [50]. In this study, miRNAs were discovered that regulate both the MVA and MEP pathways in terpenoid biosynthesis in P. minor (Fig. 11). Four miRNAs, pmi-miR6300, pmi-nov_12, pmi-miR530, and pmi-miR6173, were found to be involved in the MVA pathway, whereas two miRNAs, pmi-miR396a and pmi-miR398f/g, were found to be involved in the MEP pathway. Most of the miRNAs regulate upstream genes in the terpenoid biosynthesis pathway.
In the MVA pathway, among the four responsive miR-NAs, three (pmi-miR6300, pmi-6173 and pmi-nov_13) were downregulated, whereas the fourth miRNA (pmi-miR530) showed a sharp increase at 9 dpi. The Pmi-miR6300-regulated gene encoded HMGR, a rate-limiting enzyme in the MVA pathway and the first rate-limiting enzyme recognized in the terpenoid biosynthesis pathway [51]. In A. thaliana, the HMGR gene has undergone a duplication process that led to HMG1 and HMG2. HMG1 is expressed throughout the plant, whereas the expression of HMG2 was observed in only the meristem and floral parts of A. thaliana plants [51]. In addition, the expression of HMGR is affected by internal factors such as plant development and phytohormone [52]. HMGR is also affected by external stimuli, such as light, wounding, elicitor treatment and pathogen invasion [52,53]. The mutant hmg1 in A. thaliana exhibited infertility, premature senescence, and dwarfness [54]. In the MVA pathway, HMGR catalyses the rate-limiting steps that directly affect downstream products. In A. thaliana, the hmg1 mutant showed a 65% reduction in triterpene compound accumulation compared to the wild type [55]. Moreover, elicitor treatments such as fungal inoculation and wounding also increased the expression level of HMGR and subsequently led to higher accumulation of sesquiterpene compounds    [53,56]. In P. minor, terpenoids are produced in small quantities under normal conditions. RT-qPCR results showed that pmi-miR6300 was upregulated at 0 dpi, which led to the suppression of its target, HMGR. Indeed, after treatment, pmi-miR6300 was downregulated and led to the high accumulation of HMGR. In X. strumarium, HMGR was targeted by miR1134 and miR5021 [50]. However, the study of both miRNA and HMGR in X. strumarium is still at the prediction level and requires further experimental validation.
In the MVA pathway, the HMGR enzyme catalyses the conversion of HMG-CoA to mevalonate, which is later converted into mevalonate-5-phosphate through the enzyme MVK. Until now, no clear interaction between miRNA and MVK has been reported in plants, while in mice, miR122 targets MVK, which is involved in cholesterol biosynthesis [1,57]. In this study, we investigated the interaction and reported that pmi-nov_13 targets MVK in P. minor. The second miRNA, pmi-miR530, was found to regulate the MVD enzyme, which catalyses mevalonate diphosphate into IPP. Pmi-miR530 was increased drastically at 9 dpi, which led to the downregulation of its target. As for pmi-nov_13 and MVK, no miRNA in plants has previously been reported to target MVD. In mice, miR124 is involved in hypocholesterolaemia by targeting MVD [58]. In plants, information about the roles of MVK and MVD is quite limited compared to yeast [59]. However, the latest study of MVK in Gingko biloba discovered that overexpression of MVK led to the accumulation of the terpene trilactone [60]. In P. minor, the accumulation of MVK at 6 dpi may contribute to sesquiterpene biosynthesis, although the rate-limiting enzyme had its highest expression at 9 dpi. In the MEP pathway, two miRNAs (pmi-miR396a and miR398f/g) were discovered to be involved by targeting DXS and DXR. DXS is involved in the upstream reaction in the MEP pathway by catalysing the condensation process between pyruvate and D-glyceraldehyde 3phosphate to produce DXP [61]. The discovery of DXS was first reported in E. coli and then in A. thaliana via mutation of the cla gene [61,62]. Transgenic A. thaliana overexpressing the DXS gene resulted in an abundance of terpenoids compared to wild plants, and DXS was categorized as a rate-limiting enzyme in the MEP pathway [61]. Further study of DXS in A. thaliana showed that DXS exists as the paralogues DSX1 and DXS2. Of these two paralogues, only DXS1 was found to be involved in terpenoid biosynthesis [63].
In A. thaliana, the majority of the DXS genes were expressed in photosynthetic organs, such as leaves, and floral parts, while homologues to DXS2 in tomato were expressed in roots and trichomes [64]. In the second step of the MEP pathway, the DXR enzyme converts DXP to MEP via rearrangement of the molecular structure followed by a reduction process by NADPH. This step is reversible. Unlike DXS, the role of DXR as a ratelimiting enzyme is still unclear, and its function may depend on plant species, organ and development stage [65]. However, similar to DXS, DXR genes are distributed in different plant parts and stimulated by light response [65]. The mutant dxr exhibited an albino phenotype, defects in gibberellin and abscisic acid biosynthesis and improper formation of trichomes and stomatal closure [66].
Interestingly, in P. minor, the majority of miRNAs involved in the MVA pathway were downregulated except for pmi-miR530, which was upregulated at 9 dpi. The decreasing pattern exhibited by pmi-miR6300, pmi-nov_ 13 and pmi-miR6173 resulted in the accumulation of their target transcripts, which most likely regulate sesquiterpene production. These findings were supported by the metabolite profile in Fig. 2, which shows prominent production of α-cedrene, valencene, and βbisabolene after treatment with F. oxysporum. Meanwhile, the increasing pattern shown by pmi-miR396a and pmi-miR398f/g resulted in the suppression of DXS and DXR in the MEP pathway. The differences in miRNA expression involved in regulating MVA and MEP are probably due to the effect of F. oxysporum, which stimulates terpenoid production in the MVA rather than the MEP pathway. Downstream products in the MVA pathway (sesquiterpene), which play crucial roles in plant defence mechanisms, may influence the selection of the terpenoid pathway in the plant system. In addition, the interaction of two regulators, miRNAs and transcription factors, can also affect terpenoid biosynthesis. In A. thaliana and P. cablin, expression of miR156 reduced the production of sesquiterpene by suppression of the target, SPL [67]. In P. minor, pmi-miR156b/c was predicted to target the SPL gene. However, there is no evidence that this interaction would affect sesquiterpene production, because a previous report showed that the interaction between miR156 and SPL is conserved only in regulating floral development [11,68].

miRNA cleavage site determination by degradome analysis
Degradome sequencing was carried out to determine the miRNA cleavage site. The sequences were deposited under accession number SRX3921398. The statistical data are shown in Table 5. There were 17,532,759 degradome sequencing reads. Filtering and removal of adaptors revealed 15,493,710 clean reads. Mapping to the P. minor transcriptome resulted in 3,079,840 sequence tags. However, to increase the hit search, another P. minor degradome library, with the accession number SRX3921610, was used together with this degradome library. The results showed that out of the 7 studied miRNAs involved in the terpenoid biosynthesis pathway, cleavage sites for pmi-miR396a and pmi-miR6300 and their targets, peroxid-ase57 and HMGR, were detected between positions 10 and 11 of the miRNA sequences (Fig. 12), and this result was also supported by previous findings [28,29].

Spearman correlation analysis
Overall, correlation analysis revealed that pmi-miR530 and pmi-miR398f/g had positive correlations with metabolite content, whereas the remaining miRNAs (pmi-396a, pmi-miR6300, pmi-miR6173 and pmi-nov_13) showed negative correlations (Fig. 13a). Pmi-nov_12 was excluded from this analysis because it was not involved in the terpenoid biosynthetic pathway. For the target genes, the majority of them displayed positive correlations with metabolite content, except for HMGR (Fig. 13b). The low abundance of αcedrene, valencene and β-bisabolene at 9 dpi may have

Conclusion
In this study, a total of 58 conserved and two novel miR-NAs responsive to F. oxysporum treatment in P. minor were identified. Five of the 58 conserved miRNAs and a novel miRNA were found to be directly involved with the MVA and MEP pathways. However, some of the miRNAs, such as pmi_miR6173, pmi-nov_12 and pmi-nov_13, might need further confirmation at the protein level against their own targets. This work provides the framework for further exploration of miRNA-mediated regulatory mechanisms in terpenoid biosynthesis in P. minor. Moreover, our study also found some miRNAs targeting mRNAs encoding transcription factors, which suggested the role of miRNAs in regulating various biological processes, including plant growth and development. The above seven studied miRNAs could be utilized to regulate secondary metabolite biosynthesis by manipulating the MVA and MEP pathways through the RNAi mechanism.

Plant and fungus cultures
P. minor plants were cultured for approximately 6 weeks on Murashige and Skoog media in a controlled chamber room at Kompleks Rumah Tumbuhan, Universiti  Metabolite profiling in P. minor was carried out at a series of time points, beginning with 3 days postinoculation (dpi), 6 dpi, and 9 dpi, and using SPME-GCMS. In addition, P. minor plants treated with sterile distilled water were used as a control (0 dpi). Every time point had three biological replicates. The SPME method was adapted for the collection of volatiles. P. minor leaves were harvested (approximately 1 g) in liquid nitrogen and ground into small pieces. The ground tissues were then promptly transferred into labelled SPME vials to avoid evaporation of volatile compounds. The vials were closed with screw caps and then heated at 65°C for 15 min in a water bath. Volatile compounds were collected by introducing an SPME needle via a septum cap. SPME fibre (100 μM polydimethylsiloxane, PDMS), which was located inside the SPME needle, was injected into the vials to absorb the volatiles. Alternatively, polyacrylate fiber also could be used for volatiles detection [70]. The equilibrium time for SPME fibre in the collection vial was set for 30 min at 55°C-60°C. GC-MS analysis was performed on an Agilent 7890A gas chromatograph (GC) directly coupled to the mass spectrophotometer (MS) of an Agilent 5975C inert MSD with a triple axis detector. The column used was a non-polar column, HP-5MS (30 m length × 0.25 mm) diameter and film thickness 0.25 μM. Helium was used as the carrier gas, with a flow rate of 1.3 mL/min. A splitless injection was set at 50°C hold for 3 min, increased to 100°C at a rate of 20°C/ minutes, and held at 250°C for 3 min. The peaks were identified by searching the NIST/EPA/NIH mass spectral library (version 11), and the results were combined in a GC-MS chromatogram.
The data were normalized according to the sum from each sample. One-way ANOVA was carried out followed by Fisher's LSD test. The P-value was set at P < 0.05.

Small RNA library construction
To generate high-quality small RNA libraries, the crucial step is to obtain high-quality total RNA. Low-quality starting material can reduce the quantity of small RNA [72]. Total RNA was extracted, and the quality (purity and concentration) was measured using a Nanodrop spectrophotometer (ND-1000) and Qubit, respectively. The integrity of the extracted RNA was determined by Bioanalyzer analysis (Agilent 2100) using an RNA 6000 chip. RNA samples with RNA integrity number (RIN) values over 7.0 were selected for small RNA library construction. All quality control steps were enclosed in Additional file 2. Mock-inoculated (C) and Fusariumtreated (F) small RNA libraries were constructed with two biological replicates, C1 and C2, and F1 and F2, respectively. For the F library, the time point at which P. minor emitted volatiles at the maximum level was selected for small RNA library construction. Mockinoculated P. minor were prepared by adding sterile distilled water to the plant and used as a control. RNA from both the C and F samples was extracted using PureLink® Plant RNA Reagent (Thermo Fisher Scientific, USA) using the recommended protocol. Briefly, approximately 0.1 g of leaf sample from P. minor was ground and transferred to 1.5 mL. A volume of 0.5 mL of Pure-Link® Plant RNA Reagent was added and vortexed. The samples were incubated at room temperature for 5 min. Next, the samples were centrifuged at 12,000×g for 2 min at 4°C. Approximately 450 μL of supernatant was transferred into new tubes, and 0.1 mL of 5 M NaCl and 0.3 mL of chloroform were added and mixed to the supernatant. The mixtures were centrifuged at 12,000×g for 2 min at 4°C, and approximately 350 μL of supernatant was transferred to new tubes, mixed, and incubated for 10 min at room temperature. The tubes were then centrifuged at 12,000×g for 2 min at 4°C, and a small pellet was observed at the bottom of each tube. The supernatant was discarded, and 75% ethanol was added to wash the pellet. Then, the tubes were centrifuged again at 12,000×g for 1 min at 4°C. The ethanol was discarded, and the pellet was diluted with 20 μL of nuclease-free water. Then, DNase treatment was carried out for each RNA sample. The concentration of RNA was determined using a Nanodrop spectrophotometer (Thermo Scientific, USA). The integrity of the RNA was assessed by running a 1% agarose gel and Bioanalyzer analysis (Agilent 2100 Bioanalyzer, USA). A small RNA library was constructed using the NEBNext Small RNA Library preparation kit as reported by a previous study [22]. Then, the libraries were sent to Universiti Malaya, Malaysia for sequencing using an Illumina HiSeq 2500™ in Rapid Run mode.

Bioinformatic analysis
Raw data were analysed using CLC Genomics Workbench version 8 [69]. The adaptors were removed, and the low-quality reads were filtered. Only clean reads with lengths between 18 nt and 30 nt were selected for further analysis. The annotation process was carried out by mapping the clean reads to the miRNA and non-coding-RNA databases miRBase version 21 and Rfam version 12 [23,73]. A maximum of two mismatches were allowed [74]. Then, each library was 11normalized to transcripts per million (TPM); normalized expression = (actual miRNA count/total count of clean reads) × 1,000,000 [75]. Differential gene expression was calculated using Baggerley's test [76]. Fold change was calculated using |log2 fold change| ≥ 2, and the P-value was set at P < 0.05 [77].

Novel miRNA prediction
Novel miRNA prediction was carried out based on a previous study [78]. For secondary structure validation, forward and reverse primers were designed using the PrimerQuest Tool by Integrated DNA Technologies (https://sg.idtdna.com/), and the parameters were set to PCR 2 Primer ( Table 6). The primer was synthesized at First Base (Malaysia). The PCR products were measured and visualized using 4% agarose gel.

Prediction of target transcript and gene ontology
PsRobot software was deployed to predict miRNA targets with a score of 4.0 to obtain more targets [29,79]. The miRNA libraries generated in this study and the transcriptomic library of P. minor (http://www.ncbi.nlm. nih.gov/sra/SRX669305) were used in this target prediction. Then, the target genes were classified through gene ontology analysis using WEGO software (http://wego. genomics.org.cn/) [80]. The target genes involved in the terpenoid pathway were simply determined by the gene annotations in the transcriptomic library [81]. The genes were mapped with the KEGG database (http://www.genome.jp/kegg/) to identify the involvement of these genes in the terpenoid pathway [82].
Validation and expression profile using real-time quantitative PCR (RT-qPCR) To validate the existence and expression of those miR-NAs and target genes involved in the terpenoid pathway under F. oxysporum inoculation, RT-qPCR was carried out using Maxima SYBR Green qPCR Master Mix (Thermo Fisher, USA) at a series of time points: 0 dpi, 3 dpi, 6 dpi, and 9 dpi. For miRNA, the miRNA sequence was used as the forward primer, and the universal primer     from Qiagen was used as the reverse primer (Table 7). For the target genes, forward and reverse primers were designed using the PrimerQuest Tool by Integrated DNA Technologies (https://sg.idtdna.com/), and the parameters were set to PCR 2 Primer intercalating dye (Table 8). 5.8 s rRNA and tubulin were used as reference genes for the miRNAs and target genes, respectively. Relative gene expression was calculated using the Livak method [83].

Degradome sequencing
To determine the miRNA cleavage site, RNA samples from C and F were pooled together and sent to BGI (China) for degradome sequencing. After retrieving the raw data, low-quality reads were discarded, and adaptors were removed using Skewer software [84]. Once again, psRobot was deployed to map the degradome sequence against the reference gene (P. minor transcriptome) and to determine the cleavage site [29].