Identification of microRNAs and relative target genes in Moringa oleifera leaf and callus

MicroRNAs, a class of small, non-coding RNAs, play important roles in plant growth, development and stress response by negatively regulating gene expression. Moringa oleifera Lam. plant has many medical and nutritional uses; however, little attention has been dedicated to its potential for the bio production of active compounds. In this study, 431 conserved and 392 novel microRNA families were identified and 9 novel small RNA libraries constructed from leaf, and cold stress treated callus, using high-throughput sequencing technology. Based on the M. oleifera genome, the microRNA repertoire of the seed was re-evaluated. qRT-PCR analysis confirmed the expression pattern of 11 conserved microRNAs in all groups. MicroRNA159 was found to be the most abundant conserved microRNA in leaf and callus, while microRNA393 was most abundantly expressed in the seed. The majority of predicted microRNA target genes were transcriptional factors involved in plant reproduction, growth/development and abiotic/biotic stress response. In conclusion, this is the first comprehensive analysis of microRNAs in M. oleifera leaf and callus which represents an important addition to the existing M. oleifera seed microRNA database and allows for possible exploitation of plant microRNAs induced with abiotic stress, as a tool for bio-enrichment with pharmacologically important phytochemicals.

Prediction of miRNA targets. In order to predict the putative genes regulated by miRNAs identified in our samples, we took advantage of psRNATarget 28 . The cDNA library used is the TAIR v10 of A. thaliana, while the scoring scheme follows the V2, released in 2017 by the authors.
Gene Ontology (GO) enrichment analysis. To better understand the biological functions underlying target genes, we performed enrichment analyses on the three categories of GO (biological processes, molecular functions, cellular compartments). GO annotation of the predicted targets was performed by uploading the list of highlighted genes in the PANTHER classification system 29 . Experimental validation of miRNAs using quantitative real-time PCR analysis. Total RNA was extracted from M. oleifera leaves, callus and cold stress treated callus using the NucleoSpin miRNA kit (Macherey-Nagel GmbH&Co., Germany), according to the manufacturer's protocol. Quality and quantity of the total RNA were evaluated by Agilent 2100 Bioanalyzer (Agilent Technologies, CA, USA) and by spectrophotometry (NanoDrop 2000, ThermoFischer Scientific, USA), respectively. cDNA synthesis was performed with a miScript II RT Kit (Qiagen, USA) from 50 ng of sample. A quantitative real-time PCR analysis of mol-miR159a, mol-miR167a, mol-miR168a, mol-miR156a, mol-miR395a, mol-miR162a, mol-miR166i, mol-miR160h, mol-miR398b, mol-miR396a was performed using miScript SYBR Green PCR Kit according to the manufacturer's protocol (Qiagen, USA) and the reaction was performed using a Rotor-Gene Q (Qiagen, USA) machine. The amplification conditions were: activation/denaturation at 95 °C for 10 min followed by 40 cycles of denaturation at 94 °C for 10 sec, annealing at 60 °C for 30 sec and extension at 70 °C for 30 sec. All reactions were performed in triplicate for each sample and 5S rRNA was used as the internal control gene. Relative expression levels of miRNAs were quantified by using the 2 −ΔΔCt method as widely described in our previous work 30 .

Results
Construction and sequencing of small RNA libraries. To identify M. oleifera miRNAs, total RNA was extracted from leaf, callus (Non Treated Callus; NTC) and cold stress treated callus (Treated Callus; TC) in order to construct nine small RNA libraries. Next, the libraries were sequenced using Illumina NextSeq. 500 System sequencing platform (Illumina, CA, USA) and analysed by bioinformatics approach. A total number of 71,900,401 raw reads was obtained from leaf (3 replicates): 72,962,350 from NTC (3 replicates) and 69,876,775 from TC (3 replicates). After removing the adaptors, low quality reads and contaminants, 63,698,781 clean reads from leaf, 71,445,523 from NTC and 68,889,789 from TC were obtained. As already reported in our previous work 20 the number of raw reads for seed tissue before and after cleaning process, was 31,290,964 and 22,737,895 -respectively. The clean and unique reads were subsequently subjected to size distribution analysis as shown in Fig. 1. The filtering analysis reported the length of clean and unique reads of small RNAs to be between 21 to 24 nucleotides (nt). The 24 nt was the most abundant size in seed and leaf (23% and 8.7%, respectively) followed by 21, 23 and 22 nt in size. In both NTC and TC, the most frequent small RNAs were 24 nt in length (8% and 12%, respectively), followed by small RNAs of 23, 21 and 22 nt in length, in TC and of 21, 22 and 23 nt in length in NTC. www.nature.com/scientificreports www.nature.com/scientificreports/ Known miRNAs in M. oleifera. As for the length distribution of mature miRNAs is concerned (high and low confidence miRNAs), 21 nt is the most frequent class with the majority reads in all samples, followed by 20, 22, 19 and 18 nt. The majority of miRNAs from leaf belongs to 21 nt category (61%, 2.9*10^5 total number of reads). The abundance of reads supporting this length is notably higher also in seed (55.1%) with a great number of reads (3.2*10^5). NTC and TC presented a highly frequent 21 nt miRNAs size (63.8% and 61.7%, respectively) with a high number of reads (1.3*10^5 and 1.2*10^5, respectively) ( Fig. 2A).
All the predictions obtained by the analysis with mirDeep2 software 25 were classified as high-confidence. Table 1 shows the most abundant high-confidence miRNAs (over 100 reads) in all experimental conditions. For each of them, a detailed set of sequence (mature, star, precursor) and structural information is reported in Supplementary Table 1. As for the expression of high-confidence miRNA isoforms, the most abundant in NTC is mol-miR398c_1 with 4.3*10^4 reads (mean value of triplicate). In TC, mol-miR398c_3 is the most frequent isoform with 3.1*10^4 reads (mean value of triplicate) while in leaf and seed, mol-miR166u_5 and mol-miR166u_2 with 1.6*10^5 and 1.3*10^5 reads (mean value of triplicate), respectively, are the most abundantly expressed isoforms (Supplementary Table 1).   TGACAGAAGAGAGTGAGCAC  399  693  589  2,427   156f 3  TGACAGAAGAGAGTGAGCA  0  0  331  0   156t 1-3  TTGACAGAAGAAAGAGAGCAC  0  0  99  126   157  157a 1-8  TTGACAGAAGATAGAGAGCAC  1,797 674  64, Continued www.nature.com/scientificreports www.nature.com/scientificreports/ A subsequent secondary comparison between the sRNA sequences of ten libraries and known mature miRNA sequences from other plants was performed. The outcome of this analysis generated a total of 300 low-confidence miRNAs that share their mature sequence with miRBase (Release 21) 27 records. All of them belong to a total of 56 miRNA families (Table 1, Supplementary Table 2). In most families, more than one precursor was identified. The most abundant isoforms resulted to be mol-miR159_1 in NTC, TC and leaf with 2.7*10^4; 3.1*10^4 and 5.1*10^4 reads (mean value of triplicate), respectively, and mol-miR393b_1 in seed with 9.7*10^3 reads (mean value of triplicate) (Supplementary Table 2).
In high and low confidence miRNA families, mol-miR166 showed a higher number of isoforms in seed (27 isoforms, 2.7*10^5 total number of reads), in leaf (23 isoforms, 1.5*10^5 total number of reads) and in TC (22 isoforms, 1.0*10^4 total number of reads), while in NTC, miRNA166 families possessed 16 isoforms with high abundance (2.6*10^4 total number of reads). Thirty-six percent of miRNA families has only one condition-specific isoform: among these, mol-miR6300 has the higher number of reads in leaf (8,719 reads), callus (7,374 and 5,057 reads in NTC and TC, respectively) and seed (4,809 reads) while mol-miR398 showed high expression only in seed (3,343 reads) (Fig. 2B).
As for low-confidence miRNAs, we assessed the correlation between conservation rate and abundance in all experimental conditions. mol-miR159_1 was the most abundantly expressed in M. oleifera leaf, NTC and TC (5.1*10^4, 2.7*10^4 and 3.1*10^4 reads respectively) with a low conservation rate among plants (rate 30). Mol-miR156a_2 was found to be the most conserved miRNA among plants (rate 165), although uncommon in all experimental conditions except for the leaf (104 reads). In addition, mol-miR395a_2, mol-miR171a_1, mol-miR164a_1, mol-miR172a_1, mol-miR167a_2 and mol-mir169b/l_1 showed a high conservation rate among plants and a relatively low abundance, except for mol-miR164a_1 in the seed. The remaining miRNAs, with the exception for mol-miR396b_1, showed a conservation rate below 50 with different abundance degree. (Fig. 3).
Novel miRNAs in M. oleifera. A total of 392 novel microRNAs were predicted by identifying all the potential precursors (pre-miRNAs) and modelling their secondary structures. All the discovered sequences ranged in length between 18 and 25 nt (Supplementary Table 3). The sequences of most of the novel miRNAs were 21 nt long (102 sequences) while the second most abundant class was 18 nt long (69 seq). The lengths of novel miRNA precursors ranged from 37 to 105 nt, where the 43 nt was most abundant (17 seq). Leaf has the highest number of novel miRNAs (213 seq). NTC, TC and seed follow with 121, 85  www.nature.com/scientificreports www.nature.com/scientificreports/ most abundant miRNA isoform was found in the leaf with 21 nt in length (mol-miR-n31/n211/n239; 7.8*10^5 reads), followed by seed (mol-miR-n323; 2.0*10^5 reads), TC (mol-miR-n217; 1.8*10^4 reads) and NTC (mol-miR-n290; 1.4*10^5 reads) with 20 nt in length.
Identification of "cold-sensitive" miRNAs. Cold stress has been shown to have an important influence on plant growth and development by modulation of gene expression. Some studies document that chilling stress increases the amount/slows down decomposition of polyphenols in some plants (Nicotiana tabacum L., Vitis vinifera L.) 31,32 . Moreover, chilling stress response induced increased expression of selected miRNAs in wild tomato (Solanum habrochaites) 17 , Brassica rapa L. 33 and Astragalus membranaceus 34 . We performed a comparative analysis between the known miRNAs identified in M. oleifera and the ones published in the aforementioned studies. As shown in Fig. 4A, 25 miRNAs are conserved among the selected species. By looking at their global abundance within different conditions (Fig. 4B), it is not possible to determine whether the subset in its entirety is sensitive to the cold stress. However, miR-166e-3p in Brassica rapa L., mol-miR319_1 and mol-miR396h_1 in M. oleifera, ame-miR396-5 in A. membranaceus are upregulated after cold stress treatment.
Identification of condition-specific and backbone miRNAs. In addition to analysis of the gene regulation mediated by miRNAs, the analysis of differential internal regulation of miRNA expression in different tissues and experimental conditions was performed. In plants, this fine-tuning mechanism is based on the regulation of miRNA's decoy time 35 . As consequence, miRNA patterns in plants are modulated according to the tissue and the environmental conditions. On the other hand, we also defined a set of miRNAs whose expression is stable across different experimental conditions (backbone). In order to evaluate the modulation of miRNA patterns among tissues and cold-stress, we overlaid the high-confidence, low-confidence and novel miRNAs identified in each condition, separately (Fig. 5). A detailed list of tissue-specific and backbone miRNAs is also reported in the Supplementary Table 4. Among 56 miRNA families belonging to all experimental conditions, mol-miR166 and mol-miR159 families are the largest families identified with ten and eight members, respectively. A total of five different isoforms belongs to the mol-miR156 and mol-miR396 families. Several miRNAs appear to be tissue-specific: 135 unique miRNAs were found in the seed; among these, mol-miR166 (8 members), mol-miR167 (8 members), mol-miR156 (6 members) and mol-miR399 (6 members) were identified as well as 46 novel miRNAs.
Prediction of miRNA targets. Plant miRNAs play important roles in diverse biological processes by cleaving target mRNAs or supressing their translation. In order to understand the biological functions of M. oleifera miRNAs, PNRD (plant non-coding RNA database) 36 was used as a reference for the prediction of putative target genes for tissue-specific and backbone miRNAs. The analysis revealed 7,677 unambiguous target genes for 54 miRNA families conserved across all the experimental conditions (Supplementary Table 5).
The majority of target genes were transcription factors; these include key regulators of abiotic and biotic stress response, plant reproduction and plant growth/development gene such as squamosa promoter-binding proteins (SPL or SBP) 37  www.nature.com/scientificreports www.nature.com/scientificreports/ Some of the identified miRNA targets are genes directly involved in resistance to biotic and abiotic stresses, including protein degradation (ubiquitin carboxyl-terminal hydrolase 1) and oxidoreductase activity (ALDH22a1).

RNA binding family protein, Methyltransferases and DEAD/DEAH box RNA helicase family protein), plant development (Pentatricopeptide repeat (PPR) superfamily protein, RING/U-box superfamily protein), translation process (Ribosomal protein) and oxidoreductase activity (2-oxoglutarate and Fe(II)-dependent oxygenase superfamily protein, Galactose oxidase).
Gene Ontology (GO) enrichment analysis. The enrichment analysis of miRNA target genes represents a fundamental procedure for evaluating the miRNA-mediated biological modulations that occur in Moringa leaves, seed, TC and NTC. The terms belonging to GO biological processes (BP), molecular functions (MF), and cellular compartments (CC) have been enriched and the results for the top 20 terms plotted in a bar plot (Supplementary Figures 1_NTC, TC, LEAF, SEED). Analyses performed on targets for backbone miRNAs, identified certain biological functions and molecular processes that are crucial for the proper survival of the plant. On the other hand, the enrichment analysis of targets for condition-specific miRNAs highlights biological modulations that are essential in a specific, biological context. Biological processes. In the leaf, the majority of miRNAs seems to be responsible for fine regulation of genes involved in leaf symmetry 52,53 ; other genes regulate the development of different plant organs, such as anther 54 , floral whorl 55 , stamen 56 and so on. In callus, most miRNAs are involved in positive heterochronic gene regulation of development, especially in stress conditions; other genes are involved in the biosynthesis and metabolic processes of beta-D-glucan 57 (callus) and the lignin/phenylpropanoid catabolic processes 58 . In the seed, most frequently expressed genes regulate the maintenance of DNA methylation, negative regulation of intracellular signal transduction and specification of axis polarity 59 .

Cellular compartments.
In leaf and seed, most of the regulations take place in the late endosome, vacuole and cell projection. In callus, most of the genes are localized in the 1.3-beta-D-glucan synthase complex (NTC) and mainly in the chloroplast endopeptidase Clp complex (TC). The presence of several genes located in the chloroplast could be correlated with the enrichment of stress-response genes. In fact, most of the transcriptional activities in the nucleus are regulated in part by signals derived from plastids. This process is named "retrograde signaling" and it is well amplified during responses to chemical, physical and biological stress 60 .

Molecular functions.
Despite most of the enriched terms do not provide a detailed description of the involved molecular function, they clearly highlight the fine-tuning role of microRNAs in plants. Indeed, most of the targets are involved in xenobiotic transmembrane transporting ATPase activity 61 (leaf), oxidoreductase and ATP-dependent peptidase activity 62 (TC), superoxide dismutase copper chaperone activity 63 (NTC), phosphorelay sensor kinase activity and xenobiotic transmembrane transporting ATPase activity 64 (seed).

Validation of conserved miRNAs by qRT-PCR.
To confirm the expression pattern of the M. oleifera miRNAs, 11 conserved miRNAs with different expression profiles were randomly selected for real-time PCR analysis. The data collected demonstrate that the expression patterns are similar between the two analytical tools (Illumina sequencing and qRT-PCR) for six of the eleven miRNAs, whereas other show different expression pattern as detected between the two molecular tools.
As illustrated in Fig. 6, mol-miR396a and mol-miR159a are more abundant than other miRNAs in leaf and callus (both NTC and TC). Mol-miR162a is more abundantly expressed in leaf than in callus (both NTC and TC), while mol-miR398b and mol-miR168a are more expressed in calluses. On the other hand, mol-miR166i and mol-miR160h show similar levels of expression in all tissues.

Discussion
In the past decade, M. oleifera has gained growing attention for its nutraceutical and pharmacological functions as well as potential benefits for human health; several studies have demonstrated the great number of bioactive compounds contained in leaves, seeds, pods and flowers 1,2,4 . The large number of these bioactive compounds (such as secondary metabolites) might explain the pharmacological properties of this "miracle tree". Recent studies have confirmed the high content of polyphenols (such as flavonoids and phenolic acids) in M. oleifera leaves 65 . These secondary metabolites are synthesized in the plant for its specialized need, in particular, ecological conditions; unlike primary metabolites, secondary metabolites are directly involved in the defence mechanism against environmental injuries. In human, M. oleifera polyphenols (such as flavonoids) are involved in protection against chronic diseases associated with oxidative stress, including cardiovascular disease and cancer 4 . Plant metabolic engineering has been used for enhancing biosynthesis of these pharmacological important phytochemicals [6][7][8][9] . One of the ways to increase polyphenol biosynthesis is by modulating levels of miRNAs -the ultimate regulators of biosynthesis and accumulation of secondary metabolites in plants [66][67][68] .
In our previous study 20 we reported a high throughput sequencing of miRNAs from M. oleifera seed; some of them were involved in regulation of human genes when transfected into Hepatoma cell lines (HepG2) 69 . Being the part of the plant that contains the highest number of bioactive compounds, the leaf was chosen as model tissue for high-throughput sequencing in this study. Moreover, callus cells, given their undifferentiated nature, represent an optimal model for the identification and quantification of miRNA types present in M. oleifera in its basal condition. Low temperature was chosen as appropriate stress condition in the experimental procedure, considering that, due to natural geographical distribution of M. oleifera, this plant hardly ever comes into contact with this type of adverse condition. Indeed, cold stress has been shown to have an important influence on plant growth and development by modulation of gene expression. Despite a preliminary, explorative analysis did not highlight www.nature.com/scientificreports www.nature.com/scientificreports/ the presence of a "cold-sensitive" signature of microRNAs, specific isoforms were up-regulated in different plants. The application of the cold stress to a callus for further induction of somaclonal variations, could open, in turn, the possibility to develop a M. oleifera cell line better adapted to a temperate climate. For these reasons, M. oleifera calluses both untreated and exposed to abiotic stress seemed to be an interesting candidate models for this study.
The recent availability of data for the M. oleifera genome, allowed to re-evaluate our previous predictions from seeds and to identify a set of high-confidence and novel miRNAs from cold-stressed and non-stressed callus and leaf. A subsequent sequence alignment (Blast) compared to all plant-related mature miRNAs annotated in miRbase, allowed for identification of an additional set of low-confidence miRNAs. In the present study, a total of 131 known high-confidence miRNA isoforms (21 miRNA families), 300 known low-confidence miRNA isoforms (56 miRNA families) and 392 novel miRNA isoforms were identified in M. oleifera leaf, seed and callus with or without abiotic stress in the ten libraries. The majority of the identified miRNAs were 21 nt in length (60%; mean value of all experimental conditions), which is the canonical size for miRNAs generated from DCL1 processing. This result is similar to that of the conserved miRNAs predicted in B. oleracea 70 and other plant species such as soybean 71 , maize 72 , switchgrass 73 and Chinese cabbage 74,75 . In our study, the number of reads of the conserved 21 nt miRNAs is very high, especially for the leaf (4.6*10^5 reads). The result was consistent with our previous reports in which sequence length 21nt has the highest number of miRNA reads in M. oleifera seeds.
Mol-miR166 and mol-miRNA156 families contain larger number of isomiRs in the seed, leaves and in the treated callus; these microRNAs are involved in many biological processes including leaf development 76 , apical dominance 77 , floral transition and development 78 . MiRNA166, together with miRNA165, belong to another important class of miRNAs involved in Shoot Apical Meristem (SAM) maintenance. These two miRNAs share the same function of targeting/repressing class III Homeodomain-Leucine Zipper (HD-ZIP III) expression. These miRNAs play an important role in meristem maintenance, adaxial identity of leaves, lateral root growth and procambium identity.
The mol-miR156 family was a largely conserved miRNA family, followed by mol-miR395, mol-miR164 and mol-miR172. In our analysis, mol-miR156 family shows higher expression in leaf compared to other experimental conditions; these miRNAs belong to the largest families highly conserved in all land plants 79 . Previous studies demonstrated their role in SPL pathway regulating plant development, flowering and plastochron length.
However, there is a large number of miRNAs present in only few other species. For instance, mol-miR6300 was discovered in soybean and subsequently identified in Finger millet and in Chickpea, as well as M. oleifera seed (as reported in our previous work 20 ). In flowers and leaves of Chickpea, this miRNA was found abundantly expressed; however, as predicted by other authors, its annotation in miRBase is ambiguous and needs experimental confirmation 80 . In our analysis, all the experimental conditions (especially for leaf) showed an up-regulation www.nature.com/scientificreports www.nature.com/scientificreports/ of this miRNA compared to other miRNAs. It would be interesting to investigate in more detail the biological function of this miRNA considering that it is highly expressed in many parts of this medicinal plant.
In this study, mol-miR159, mol-miR393 and mol-miR396 were found highly expressed in all experimental conditions. These results are in line with different studies, comprising various plant species, in which miRNA159 was shown to be very frequently detectable miRNA involved in fundamental plant biology roles, such as plant growth, and development 81 . In A. thaliana, miRNA159 has been shown to regulate anther and silique development by targeting MYB33. Mutation of the miRNA159-binding site on MYB33 mRNA resulted in pleiotropic defects including severely impaired fertility, stunted anthers, small siliques, and small seeds. Interestingly, miRNA159 accumulation was up-regulated by gibberellin (GA) application and GA-deficient mutants showed low miRNA159 accumulation. Treating these mutants with GA was itself sufficient to increase the accumulation of miRNA159 to wild-type levels and above, demonstrating the interplay between miRNAs and hormones in plant development 82 .
It is noteworthy to mention that miRNA159 is involved in cross-kingdom regulation of mammalian gene expression; Chin et al. found that plant miRNA159 could be detected in human sera and its levels were inversely correlated with breast cancer incidence and progression targeting human Transcription Factor 7 (TCF7) gene 83 .
Gene targets were computationally predicted for both conserved and novel miRNAs in order to elucidate the biological functions of miRNAs in M. oleifera leaf, seed and callus. As demonstrated by other studies, mol-miR159, very common in all experimental conditions, was predicted to target MYB transcription factors; these target genes have been reported to play an important role in abscisic acid (ABA) signalling during A. thaliana seed germination 84 . Moreover, mol-miR159 and its target gene (GAMYB) were involved in modulation of grapevine floral development in response to gibberellin (GA) treatments and this interaction has important implications for the molecular breeding of high-quality seedless grapevine berry. Mol-miR393, more common in seed than in leaf and callus, was shown to be involved in regulating various aspects of plant growth and development; for instance, mol-miR393 targets the Transport Inhibitor Response 1 (TIR1) and Auxin Signaling F-Box (AFB) genes, involved in flower development thus contributing to the maize grain filling rate by regulating maize growth, development and environment stress response 85 . According to our analysis, this miRNA regulates enzymes involved mainly in Hydrolase and Methyltransferase activity including transcriptional regulation. Growth-Regulating Factors (GRF) genes are targeted and regulated by mol-miR396, very common in M. oleifera seed and leaf. In A. thaliana leaves in phase of development, expression of miRNA396 in the distal region of the leaf blade restricts GRF activity, thereby confining cell proliferation to the proximal leaf blade. As the leaf matures, miRNA396 level increases in the developed leaf leading to decreasing levels of GRF, stopping leaf blade growth. This conserved miRNA and its targets produced similar result in our analysis. Our results showed a very high conservation rate for mol-miR156 and a moderate abundance in M. oleifera leaf. This miRNA suppresses SPL and represses flowering in many plant species. As reported by various authors, the miRNA156 family is one of the largest miRNA gene families and members of the family are highly conserved in all land plants.
The predicted potential targets of M. oleifera novel miRNAs were involved in different molecular functions such as transcription factor activity, Auxin signalling and Ethylene responses. Many of these target genes are involved in plant growth and development, or in stress responses. Furthermore, our predicted analysis found a novel miRNA (mol-miR-n111) overexpressed in treated callus respect to non-treated callus that regulate many target genes, including MYB activity; previous reports have indicated that this transcription factors is also a target of miRNA159 86,87 In our analysis, mol-miR159 was up-regulated in all experimental conditions; in plants, MYBs are involved in different processes including primary and secondary metabolism, cell fate and identity, development, and responses to biotic and abiotic stresses including cold tolerance. The novel mol-miR-n323, commonly expressed in seed, was predicted to regulate different target genes, including Histone Deacetylase 5 (HDA5). In A. thaliana, this protein is present in a protein complex involved in the regulation of flowering time 88 ; it's regulated by mol-miR482, uncommon in seed. In eukaryotes, epigenetic mechanisms (such as acetylation) play a crucial role in the regulation of gene expression; the HDACs protein complex has been associated with various developmental processes such as flowering, floral organ identity, seed development and circadian clocks.
In order to identify conserved and novel miRNAs involved in the biosynthesis of secondary metabolites in plants, we performed an overlaid analysis to identified specific miRNAs. Mol-miR408, more common in treated callus respect to others experimental conditions, are involved in the biosynthesis of benzyliso-quinoline alkaloids (BIA). This alkaloid is synthesized mainly by an agronomically and economically important medicinal plant, the Opium poppy (Papaver somniferum L.), used for main morphine alkaloid production 67 . In Argemone mexicana L. (Papaveraceae), BIA is used to treat different disorders, given its antimicrobial, antiparasitic, antimalarial, pesticide, cytotoxic and neurological properties 89 . Another important miRNA implicated in the biosynthesis of the secondary metabolites in plants is mol-miR156. This miRNA, targeting SPL9 gene, is up-regulated in M. oleifera leaf and seed and seems to be involved in the Anthocyanin and Sesquiterpenoid/triterpenoid biosynthesis. Anthocyanins are flavonoids with antioxidant properties and can therefore can be used potentially as dietary nutraceuticals for human health. Data from several epidemiological studies have reported an inverse correlation between anthocyanin intake and risk of cardiovascular disease (CVD) or CVD-related mortality 90 . Different authors have shown a beneficial effect of anthocyanin rich foods on gut microbiota; the result of in vitro and in vivo studies highlighted a significant proliferative effect on Bifidobacterium spp., known for their wide use in probiotics and for the treatment of irritable bowel syndrome, as well as inhibition of Clostridium histolyticum, pathogenic in humans 91 . Mol-miR396b_1, upregulated in M. oleifera seed and leaf, is well known for being able to regulate the biosynthesis of Flavonol glycoside by target Kaempferol 3-O-beta-D-galactosyltransferase 92 . Dietary flavonoids isolated from different medicinal plants have received an increased attention due to their considerable benefits in the prevention and management of modern diseases such as cancers, diabetes, and cardiovascular diseases 93 . www.nature.com/scientificreports www.nature.com/scientificreports/ In conclusion, this is the first comprehensive identification of conserved and novel miRNAs in M. oleifera leaf, seed and callus (stressed and non-stressed). This dataset represents an important supplement to the existing M. oleifera miRNA database and shall be useful in understanding the possible role for enhancing biosynthesis of pharmacologically important phytochemicals by plant miRNAs during biotic and abiotic stresses. Further studies are necessary in order to elucidate this complex regulatory network potentially able to improve human health in socially neglected populations by a "miracle" tree with high nutritional and medicinal value.

Data availability
The sequences from the small RNA library have been deposited in the Gene Expression Omnibus (GEO) database, the accession number is GSE119247 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE119247).