MicroRNAome Profile of Euphorbia kansui in Response to Methyl Jasmonate

miRNAs play vital regulatory roles in different plant developmental stages and in plant response to biotic and abiotic stresses. However, information is limited on the miRNA regulatory mechanism to methyl jasmonate (MeJA). In this study, we used the microRNAome profile to illustrate the relevant regulatory mechanisms of Euphorbia kansui in response to methyl jasmonate (MeJA) through Illumina RNA-Seq. As a result, we identified 875 miRNAs corresponding to 11,277 target mRNAs, among them, 168 known miRNA families representing 6019 target mRNAs sequences were obtained. 452 miRNA-mRNA pairs presented an anti-correlationship (Cor < −0.50 and p-value of correlation ≤ 0.05). The miRNA with a fold change ≥ 2 and a p (p-Value) < 0.05 in pairwise comparison were identified as significant differentially expressed miRNAs (DEMs). The DEMs in MeJA treatment of 0, 24, 36 and 48 h were compared by using Short Time Expression Miner (STEM) cluster and 4 significant gene profiles (p-value ≤ 0.02) were identified. Through the kyoto encyclopedia of genes and genomes (KEGG) pathway and gene ontology (GO) enrichment analysis on all miRNA targets, we identified 33 mRNAs in terpenoid biosynthesis, which were regulated by miRNAs under MeJA treatment, so the miRNA maybe involved in the response of E. kansui plant to exogenous MeJA and the results would provide very useful information on illustrating the regulatory mechanism of E. kansui and also provide an overall view of the miRNAs response to MeJA stress of a non-model plant.


Introduction
Euphorbia kansui Liou is a classical traditional Chinese medicine (TCM), which dried root has long been used to treat asthma, edema, ascites [1]. There are many secondary metabolites in E. kansui, diterpense and triterpenes are the main compounds, which show anti-virus, skin irritating and modulation of multidrug resistance effects [2][3][4]. The extracts in folk are widely used in anti-tumor, anti-cancer, anti-fertility, anti-virus, anti-inflammatory activity, regulating the immune system and other activities [5][6][7][8].
Jasmonic acid (JA) and its derivatives referred to as jasmonates, are key signaling molecules and play roles in many biological processes, such as growth inhibition, senescence, wound response, plant defense and the secondary mechanism [9,10]. The JA system is an important component in complex plant hormone signaling systems [11], which can elicit three major classes of plant secondary metabolites, terpenoids, alkaloids and phenyipropanoids [9,12]. Methyl jasmonate (MeJA) is the most important kind of JA system, which has also been used to elicit defense responses in many species through enhancing the secondary metabolites production [13][14][15][16], such as volatile terpenoids from Amomum villosum [17], triterpene in Euphorbia pekinensis [18] and tropane alkaloids in Hyoscyamus niger [19].
MicroRNAs (miRNAs) are short non-coding RNAs with a length of [19][20][21][22][23][24][25] nucleotides. Many studies have demonstrated that miRNAs play crucial roles in a wide range of biological processes and stress responses [20][21][22][23]. In plant, miRNAs are post-transcriptional regulators of gene expression, which are related to growth, development and stress, and in different stages of plant development the miRNAs expression differs [24][25][26][27]. Biswas et al. [28] reported that miRNAs play an important role in different processes in a post-transcriptional manner by down-regulating target gene products. Shriram et al. [29] also reported that when plants are stressed, miRNAs downregulated their target mRNAs and acculated the accumulation and function of positive regulators. Some studies have found that miRNAs are involved in the drought stress of wheat [30,31] or dehydration stress in Triticeae [32]. A large-scale multiomics analysis showed that miRNAs involved the potential regulatory in wheat stem sawfly [33]. Shen et al. reported that some miRNAs of Oryza sativa responded to a drought, salt and the cold [34].
At present, miRNA have attracted more and more attention, and in different plants where they are more conserved and novel miRNAs have been identified [35,36]. Some miRNAs were regulated to respond to stress conditions, including MeJA, salicylic acid (SA), abscisic acid (ABA), melatonin and rice stripe virus (RSV), which suggests that miRNAs direct post-transcriptional regulation of their respective target genes to replay the stress [34,[37][38][39][40].
Clarifying the regulatory relationships between miRNAs and their corresponding mRNA targets, many biological questions have been further understood [40][41][42]. Expression level analysis of miRNAs is particularly important in exploring their biological functions. Using Illumina RNA-Seq, a second-generation sequencing-based technology, it is possible to measure miRNA expression levels in tissues of interest. In this study, we used Illumina RNA-Seq to examine the microRNA transcriptome of the control group and MeJA treatment groups at different handling times of Euphorbia kansui. Through analysis, the microRNAome profile was used to illustrate the relevant regulatory mechanisms of E. kansui in response to MeJA. The results showed that miRNAs involved in the response of E. kansui plant to exogenous MeJA, which will provide very useful information on illustrating the regulatory mechanism of E. kansui and also provide an overall view of miRNAs response to MeJA stress of a non-model plant.

Differentially Expressed miRNAs (DEMs) Respond to Methyl Jasmonate (MeJA) Treatment
We analyzed differentially expressed miRNAs (DEMs) to study the regulation of miRNA in MeJA treatment. From Figure 2a,b, we could obtain the number of differentially expressed total miRNA and known miRNA respectively.

Differentially Expressed miRNAs (DEMs) Respond to Methyl Jasmonate (MeJA) Treatment
We analyzed differentially expressed miRNAs (DEMs) to study the regulation of miRNA in MeJA treatment. From Figure 2a,b, we could obtain the number of differentially expressed total miRNA and known miRNA respectively.   The expression patterns of these DEMs in response to MeJA treatment of 0, 24, 36 and 48 h were compared by using a Short Time Expression Miner (STEM) cluster of 4 significant gene profiles (p-value ≤ 0.02): STEM profile 5 (150 genes), profile 14 (38 genes), profile 19 (10 genes) and profile 3 (30 genes) were identified ( Figure 3). It could be found that in profile 5 and 3, after MeJA treatment 24 h, the miRNA expression levels decreased to the lowest. While in profile 14, miRNA expression levels increased to the highest at MeJA treatment 24 h, and in profile 19, miRNA expression levels increased directly to the highest at MeJA treatment 48 h.

Target RNA Prediction and Enrichment Analysis
In total, 875 conserved miRNAs, corresponding to 11,277 target RNA in 8 samples were identified; among them, 168 known miRNA families representing 6019 target RNA were obtained. The details of the target RNA are shown in Table S3. In known miRNA, the most abundant miRNA was the MIR169 family, accounting for 358 of the total miRNA. MIR156 (308), MIR166 (258), MIR171 (237) and MIR395 (232), family also had high abundance of expression (Table S4). In addition, we identified 707 novel miRNA, corresponding to 5258 RNA sequences.
In general, the target mRNAs were negatively regulated by their mature miRNA with either translational repression or mRNA degradation. According to correlation between miRNA-mRNA, an association of heatmap was built in Figure 4. 452 miRNA-mRNA pairs presented an anti-correlationship between the expression level of miRNAs and mRNAs (Cor < −0.50 and p-value of correlation ≤0.05) (Table S5).

Target RNA Prediction and Enrichment Analysis
In total, 875 conserved miRNAs, corresponding to 11,277 target RNA in 8 samples were identified; among them, 168 known miRNA families representing 6019 target RNA were obtained. The details of the target RNA are shown in Table S3. In known miRNA, the most abundant miRNA was the MIR169 family, accounting for 358 of the total miRNA. MIR156 (308), MIR166 (258), MIR171 (237) and MIR395 (232), family also had high abundance of expression (Table S4). In addition, we identified 707 novel miRNA, corresponding to 5258 RNA sequences.
In general, the target mRNAs were negatively regulated by their mature miRNA with either translational repression or mRNA degradation. According to correlation between miRNA-mRNA, an association of heatmap was built in Figure 4. 452 miRNA-mRNA pairs presented an anti-correlationship between the expression level of miRNAs and mRNAs (Cor < −0.50 and p-value of correlation ≤0.05) (Table S5).
terms, including 146 DEGs in cellular process and 160 DEGs in metabolic process. 363 DEGs (T1-VS-T2) in biological process were enriched in 699 GO terms, including 256 DEGs in the cellular process and 249 DEGs in the metabolic process. 273 DEGs (T1-VS-T3) in the biological process were enriched in 687 GO terms, including 272 DEGs in the cellular process and 267 DEGs in the metabolic process. 137 DEGs (T2-VS-T3) in the biological process were enriched in 396 GO terms, including 91 DEGs in the cellular process and 91 DEGs in the metabolic process.

Verification of DEM Profiles Associated with MeJA Treatment with Real-time Quantitative PCR Detecting System (QPCR)
In order to experimentally verify the DEMs profile obtained by sequence data associated with MeJA treatment, Real-time Quantitative PCR Detecting System (QPCR) analysis was conducted respectively. A total of seven miRNA in the terpenoid biosynthesis pathway, including MIR1446-x, MIR535-y, MIR845-y, MIR394-y, MIR5563-x, novel-m0022-5p, novel-m0346-3p were selected to analyze (two biological repeats per miRNA), the QPCR verification results of seven miRNA in the terpenoid biosynthesis calculation as shown in Table S9 and miRNA expression patterns can be seen from Figure 6. The results suggested that the expression levels of the miRNAs verified were mostly consistent with those transcriptome data (Table S8)

miRNA Involved in the Response of Euphorbia kansui Plant to Exogenous MeJA
The methyl jasmonate (MeJA) has been identified as a vital cellular regulator that mediates the secondary metabolism containing biosynthesis of flavonoids, alkaloids and terpenoids and so on. Many studies have showed that miRNAs are involved the potential regulatory when plants are stressed by drought, dehydration, salt and cold [30][31][32][33][34]. In the present study, in order to find MeJA-mediated miRNAs, we used the E. kansui leaves after MeJA mock-treatment 24, 36 and 48 h for miRNA analysis to profile their transcriptional alterations in response to MeJA elicitation. 875 conserved miRNAs corresponding to 11,277 target RNAs were identified, which contained 168 known miRNAs and 707 novel miRNAs. Among miRNAs identified, differentially expressed miRNAs in response to MeJA were obtained. From expressed profiles and the correlation between DEGs and DEMs, 452 miRNA-mRNA pairs presented an anti-correlationship between the expression level of miRNAs and mRNAs were found. As a result, miRNA maybe involved in the response of the E. kansui plant to exogenous MeJA, and the expression patterns of miRNA regulate various metabolic pathways.

miRNA Involved in the Regulation of Terpenoid Biosynthesis after Exogenous MeJA Treatment in Euphorbia kansui
MicroRNAs (miRNAs) are endogenous and short non-coding small RNAs, which participate in biological and metabolic processes in plants through negatively regulating their target mRNA. The negative interaction between the miRNA and its target mRNA regulates its expression through mRNA cleavage, translational repression, or epigenetic modifications [43,44].

miRNA Involved in the Response of Euphorbia kansui Plant to Exogenous MeJA
The methyl jasmonate (MeJA) has been identified as a vital cellular regulator that mediates the secondary metabolism containing biosynthesis of flavonoids, alkaloids and terpenoids and so on. Many studies have showed that miRNAs are involved the potential regulatory when plants are stressed by drought, dehydration, salt and cold [30][31][32][33][34]. In the present study, in order to find MeJA-mediated miRNAs, we used the E. kansui leaves after MeJA mock-treatment 24, 36 and 48 h for miRNA analysis to profile their transcriptional alterations in response to MeJA elicitation. 875 conserved miRNAs corresponding to 11,277 target RNAs were identified, which contained 168 known miRNAs and 707 novel miRNAs. Among miRNAs identified, differentially expressed miRNAs in response to MeJA were obtained. From expressed profiles and the correlation between DEGs and DEMs, 452 miRNA-mRNA pairs presented an anti-correlationship between the expression level of miRNAs and mRNAs were found. As a result, miRNA maybe involved in the response of the E. kansui plant to exogenous MeJA, and the expression patterns of miRNA regulate various metabolic pathways.

miRNA Involved in the Regulation of Terpenoid Biosynthesis after Exogenous MeJA Treatment in Euphorbia kansui
MicroRNAs (miRNAs) are endogenous and short non-coding small RNAs, which participate in biological and metabolic processes in plants through negatively regulating their target mRNA. The negative interaction between the miRNA and its target mRNA regulates its expression through mRNA cleavage, translational repression, or epigenetic modifications [43,44].
Due to the lack of an available genome sequence for E. kansui, we found 424 DEMs with KEGG pathway annotation. The pathways involved in terpenoid biosynthesis including terpenoid backbone biosynthesis pathway, sesquiterpenoid and triterpenoid biosynthesis pathway, steroid biosynthesis pathway, ubiquinone and other terpenoid-quinone biosynthesis pathway, and diterpenoid biosynthesis pathway were annotated in different comparisons of control and three MeJA treatments by KEGG analysis (Table 1).
We identified 33 mRNAs in terpenoid biosynthesis, which were regulated by miRNAs under MeJA treatment (Table S8). At present, information on miRNA involvement in the biosynthesis of the secondary metabolites in plants is limited. Some reports were focused on the involvement of some special miRNA in the regulation of a plant's secondary metabolism [45][46][47][48][49]. Shen et al. [50] found that auxin (indole acetic acid, IAA) repressed the expression of key terpenoid indole alkaloid pathway genes in Catharanthus roseus seedlings. Hazra et al. [51] reported that some non-responsive genes of reactive oxygen species were controlled by MeJA through the down regulation of five secondary metabolite biosynthesis specific miRNAs. From our results, we found the high abundance of expression was focused on MIR156, MIR166, MIR169, MIR171 and MIR395 family. It was previously reported that MIR156, MIR166, MIR169, MIR171 were all connected with the development of seedling. MIR395 family participated in the development of Euphorbiaceae and might be involved in metabolic processes [52]. These MIR families were up-regulated or down-regulated under the inducing of cold, drought, tension, compression, dehydration, salinity, treatment of phytohormone abscisic acid and ultraviolet-B (UV-B) radiation and other abiotic and biotic stresses [53][54][55][56][57]. In the present studies, we found some miRNAs in a different pathway of terpenoid biosynthesis, which were induced by MeJA (Table S8). We selected seven miRNA's in the terpenoid biosynthesis pathway to verify the QPCR analysis. From Figure 6, we observed the verification results of seven miRNAs, which were mostly in agreement with the miRNA-seq data.
MIR1446-x targets unigene063000 coding PRPOL, which takes part in ubiquinone and other terpenoid-quinone biosynthesis, MIR845-y targets unigene014530 coding DHCR24, which is an important enzyme in steroid biosynthesis, novel-m0022-5p targets Unigene003418 coding MHGR, which is an important enzyme in terpenoid backbone biosynthesis, MIR394-y targets unigene069994 coding VTE4, which is an important enzyme in ubiquinone and other terpenoid-quinone biosynthesis, MIR5563-x targets unigene069994 which is a sesquiterpene synthase) (Table S10). All these indicated that these DEMs were involved in the regulation of terpenoid biosynthesis by their target mRNA at an early stage of MeJA treatment.

miRNA Involved in Other Pathways
From the miRNA targets KEGG enrichment analysis, we found that some miRNA involved in other pathways, such as "Glutathione metabolism" and "Ether lipid metabolism" (Table S10) were affected by MeJA obviously.
Glutathione (GSH) is a low molecular weight tripeptide, which is playing a vital role in metabolism and cell function, e.g., providing homeostasis, involvement in the ascorbate-glutathione cycle and the detoxification of xenobiotics, and GSH is also important to detoxify the exogenous substances, isolate heavy metal and take part in other stress [58]. In our study, we found Unigene028077, which was predicted to be glutathione reductase, cytosolic-like (GR), was an important component of the antioxidant machinery to respond to abiotic stress in plant and also play an important role in tolerance to abiotic stress [59]. GR was regulated by novel-m0076-3p; novel-m0247-5, novel-m0415-3p, novel-m0438-5p, novel-m0439-5p and novel-m0639-5p.
Miao et al. reported that miRNA-34a target genes were enriched in many pathways, including ether lipid metabolism, alpha-linolenic and so on [60]. From the KEGG enrichment, we found that the miRNA6483-y and novel-0375-3p target Unigene027979 in ether lipid metabolism were affected by MeJA. Moreover, Unigene027979 was predicted to be choline/ethanol aminephosphotransferase 1 (Jatropha curcas), which is the key enzyme in the final step of the synthesis of PtdCho (phosphatidylcholine) via the Kennedy pathway.

Plant Material and MeJA Treatment
Healthy Euphorbia kansui plants were grown at the Botanical Garden of Northwest University in Shaanxi Province (Shaanxi, China). The 2 to 3 year old plants at the same development stage were assigned to four groups (Two biological repeats per group), sprayed with 20 µM MeJA (in Milli-Q water). The MeJA solution was sprayed as a fine mist to completely wet the adaxial side of each leaves. Leaves from the group treated with MeJA 24 (T1), 36 (T2) and 48 (T3) h were collected separately and directly into liquid nitrogen. All materials were stored at −80 • C until analysis.

MicroRNA Sequencing
Total RNA from four groups (two biological repeats per group) was separately isolated from leaves using RNAprep Pure Plant Kit (TIANGEN, DP441, Beijing, China), according to the manufacturer's instructions. The RNA molecules in a size range of 18-30 nt were separated and purified by polyacrylamide gel electrophoresis (PAGE). The 5 adapters were ligated to the RNAs as well, then the 3 adapters were added and the 36-44 nt RNAs were enriched. The ligation products were reverse transcripted by PCR amplification and the 140-160 bp size PCR products were enriched to generate a cDNA library and sequenced using Illumina HiSeqSE50 by Sagene Biotech Co., Ltd. (Guangzhou, China). E. kansui leaves of four treatment groups (two biological replicates for each treatment) were sequenced and have been deposited into the NCBI (National Center for Biotechnology Information) Sequence Read Archive database (accession number: SRP126443).

Sequencing Data Processing and Analysis
Reads obtained from the sequencing machines included dirty reads containing adapters or low quality bases which would have affected the following assembly and analysis. To get clean tags, raw reads were further filtered through after the low quality reads were removed, which contained more than one low quality (Q-value ≤ 20) base or contained the unknown nucleotides (N), reads without 3 adapters, reads containing 5 adapters, and then removed the reads containing 3 and 5 adapters but no small RNA fragment between them. From the reads containing polyA in small RNA fragments and reads shorter than 18 nt (not include adapters), we obtained the clean tags.
All of the clean tags were aligned with small RNAs in the GeneBank database (Release 209.0), Rfam database (11.0) and reference genome to identify and remove rRNA, scRNA, snoRNA, snRNA and tRNA. The tags that mapped to exons or introns might have been fragments from mRNA degradation and the tags that mapped to repeat sequences, all of these tags were removed.
All of the clean tags were then searched against the miRBase database (Release 21) to identify known miRNAs of E. kansui, then those species the miRNAs alignment with other species was a dependable way to identify the known miRNAs. All of the unannotated tags were aligned with reference transcriptome. According to their transcriptome positions and hairpin structures predicted by software Mireap_v0.2, the novel miRNA candidates were identified.

miRNA Expression Profiles
Total miRNA consists of existing miRNAs, miRNA, known miRNA and novel miRNA, miRNA expression level was calculated based on their expression in each sample and normalized to transcripts per million (TPM) according to the following formula. TPM = Actual miRNA counts/Total counts of clean tags × 106 In addition, the expression of existing miRNA, known miRNA and novel miRNA was also analyzed individually. Meantime, the analysis of miRNA families was used to identify if the miRNAs existed in other species. The analysis result was marked as "+" or "−" which referred to existing or non-existing, respectively.
The mRNA-sequence on the four groups described in 4.1 (two biological repeats per group) were also analyzed using Illumina HiSeqTM 4000 (Gene Denovo Biotechnology Co., Guangzhou, China) and the mRNA transcriptome of Euphorbia kansui in response to MeJA were obtained. The sequences have been deposited into the NCBI Sequence Read Archive database (Accession number: SRP126436).
The correlation analysis between miRNA-mRNA profiles were performed using the method of correlation coefficient. Correlation < −0.50 is the negative RNA-MIR profile association in response to MeJA treatment.

Differentially Expressed miRNAs Analysis and Target Prediction
Differentially expressed miRNAs (DEMs) between the control and three treatments were identified with a fold change ≥ 2 and p-value < 0.05 in a comparison as significant DEMs. The calculation formula was as follows. Based on the sequences of the known miRNAs and novel miRNAs, the candidate target genes were predicted using patmatch softwares (v1.2, https://www.arabidopsis.org/cgi-bin/ patmatch/nph-patmatch.pl).
miRNA target genes function categorization and pathway enrichment were performed based on the GO and KEGG databases. DEGs were subjected to enrichment analysis of KEGG pathways. The calculation formula was as follows (N represents the number of all genes with pathway annotation, n represents the number of DEGs in N, M is the number of all genes that are annotated to the specific pathway terms; m represents the number of DEGs in M). This analysis could recognize the most biochemical metabolic pathway and signal transduction pathways for differentially expressed genes (DEGs).
DEGs were also subjected to enrichment analysis of GO functions. The calculation formula is as above (N represents the number of all genes with GO annotation, n represents the number of DEGs in N, M is the number of all genes that are annotated to the specific GO terms; m represents the number of DEGs in M). This analysis could recognize the main biological functions of DEGs.

Target Gene Prediction
Based on the sequence results of all miRNAs, including the existing miRNAs, known miRNAs and novel miRNAs, we used the software patmatch (v1.2, https://www.arabidopsis.org/cgi-bin/ patmatch/nph-patmatch.pl) to predict the target genes. The candidate target mRNAs were predicted as follows: (1) No more than 4 mismatches between sRNA and target (G-U bases count as 0\.5 mismatches) and no more than two adjacent mismatches in the miRNA/target duplex; (2) no adjacent mismatches in positions 2\−12 of the miRNA/target duplex (5\ of miRNA) and no mismatches in positions 10-11 of miRNA/target duplex; (3) No more than 2\.5 mismatches in positions 1\−12 of the of the miRNA/target duplex (5\ of miRNA);. (4) Minimum free energy (MFE) of the miRNA/target duplex should be \>\= 74% of the MFE of the miRNA bound to its perfect complement.

Quantitative PCR (QPCR) Analysis and Statistical Analysis
In order to verify the miRNA, total RNA was separately isolated from leaves as mentioned above, then All-in-One™ miRNA First-Strand cDNA Synthesis Kit (GeneCopoeia QP013, Rockville, MD, USA) was used for reverse transcription, according to the manufacturer's instructions. The reaction mixture was gently mixed, after brief centrifugation and incubation for 60 min at 37 • C. Then at 85 • C incubated for 5 min (termination of the reverse transcription reaction). Expression levels of miRNA were investigated by Quantitative PCR (QPCR) (ABI Stepone Plus, Boston, MA, USA). The QPCR analysis was performed in a 20 µL reaction system using GoTaq qPCR Master Mix (A6002, Promega, Fitchburg, WI, USA), which contained 10 µL 2 ×qPCR Master Mix, 0.5 µL Forward primer (10 µM), 0.5 µL Reverse primer (10 µM), 1.0 µL cDNA template (the reverse transcription product is diluted 5 times) and 8.0 µL ddH 2 O. The small RNA-specific primers for QPCR are listed in Table 2. The QPCR conditions consisted of an initial denaturation at 95 • C for 10 min, followed by 40 cycles of denaturation at 95 • C for 15 s, annealing at 60 • C for 20 s, and extension at 72 • C for 20 s, then dissociation at 95 • C for 15 s, 60 • C for 1 min and 95 • C for 15 s. Experiments were performed in triplicate. The comparative C t method was used to analyze the data, firstly the 2 −∆∆Ct was calculated and then the standard (STD) was calculated for analysis.

Conclusions
The results of methyl jasmonate responsive microRNA transcriptome on Euphorbia kansui showed that miRNA may be involved in the response of E. kansui plant to exogenous MeJA. KEGG pathway analysis, GO enrichment on RNA-Seq data and QPCR verification showed that miRNAs might play roles on terpenoid biosynthesis through their target RNAs in E. kansui response to MeJA. In addition, based on the results of KEGG enrichment analysis on miRNA targets, miRNA might also be involved in glutathione metabolism and ether lipid metabolism. However, further studies would be needed in the future.