MicroRNA from Moringa oleifera: Identification by High Throughput Sequencing and Their Potential Contribution to Plant Medicinal Value

Moringa oleifera is a widespread plant with substantial nutritional and medicinal value. We postulated that microRNAs (miRNAs), which are endogenous, noncoding small RNAs regulating gene expression at the post-transcriptional level, might contribute to the medicinal properties of plants of this species after ingestion into human body, regulating human gene expression. However, the knowledge is scarce about miRNA in Moringa. Furthermore, in order to test the hypothesis on the pharmacological potential properties of miRNA, we conducted a high-throughput sequencing analysis using the Illumina platform. A total of 31,290,964 raw reads were produced from a library of small RNA isolated from M. oleifera seeds. We identified 94 conserved and two novel miRNAs that were validated by qRT-PCR assays. Results from qRT-PCR trials conducted on the expression of 20 Moringa miRNA showed that are conserved across multiple plant species as determined by their detection in tissue of other common crop plants. In silico analyses predicted target genes for the conserved miRNA that in turn allowed to relate the miRNAs to the regulation of physiological processes. Some of the predicted plant miRNAs have functional homology to their mammalian counterparts and regulated human genes when they were transfected into cell lines. To our knowledge, this is the first report of discovering M. oleifera miRNAs based on high-throughput sequencing and bioinformatics analysis and we provided new insight into a potential cross-species control of human gene expression. The widespread cultivation and consumption of M. oleifera, for nutritional and medicinal purposes, brings humans into close contact with products and extracts of this plant species. The potential for miRNA transfer should be evaluated as one possible mechanism of action to account for beneficial properties of this valuable species.


Introduction
Moringa oleifera Lam, a naturalized species from the monogenus family Moringaceae, is one of the best known, most widely distributed and most useful nutritional and medicinal plants [1][2][3]. Several organs of the Moringa tree are edible (e.g., pods, seeds, flowers and leaves) and are used in many countries (including many regions of Africa) for their high nutritional value [1][2]. Almost all tissues of this plant can be used in the treatment of inflammation or infectious diseases along with cardiovascular, gastrointestinal, haematological and neoplastic diseases [2][3][4]. The leaves are a source of natural antioxidants [5], vitamins A, B and C, minerals, proteins and essential amino acids [1][2]. The Italian Ministry of Health in compliance with the European Pharmaceutical Plant legislation has included Moringa seeds in the "List of Plant and Vegetal Integrators".
In 2014, Jung IL reported a tumor suppressor activity in in mammalian cells treated with cold water-soluble extract of M. oleifera leaves [6]. The author found an abnormal ribosomal RNA (rRNA) pattern and down-regulation of many genes and proteins involved in cell transformation and proliferation in mammalian cells treated with this extract. He concluded that the cold water-soluble extract of M. oleifera induced rRNA degradation. In 2015, Tian and coworkers reported a high-quality draft genome sequence of M. oleifera and compared the genome to related woody plant genomes in order to clarify the derivation of this species [7].
Plant miRNAs are a class of 18-24 nucleotide (nt) small, non-coding RNA that negatively regulate specific messenger RNA (mRNA). MiRNAs operate in a sequence-specific manner and silence specific protein-coding genes at the post-transcriptional level by targeting the 3' untranslated region (3'UTR) of mRNA [8]. This process causes mRNA cleavage and decreases protein translation [9][10][11]. In general, miRNAs are key regulators of development, stress response, growth and other important physiological processes [12][13].
In 2012, Zhang and collaborators demonstrated that osa-miR-168a and other exogenous miRNAs that are abundant in rice plants can be acquired by mice through food intake, as evidenced by their presence in sera or tissues of the mammalian. In vitro and in vivo functional studies showed that these exogenous miRNAs are able to inhibit mammalian gene expression in the liver, demonstrating the first case of cross-kingdom regulation [14]. More recently, oral administered cocktails of endogenous tumor suppressor miRNAs, exhibiting characteristics of plant miRNAs, reduced tumor burden in a mouse model of colon cancer [15]. These observations show that plant miRNAs are absorbed in the mammalian digestive tract and can target mammalian genes. Moreover, they suggest the hypothesis that engineered edible plants producing mammalian tumor suppressor miRNAs might be a new treatment modality for cancer. Such treatment might be an effective, nontoxic, and inexpensive chemo-preventive strategy for human.
Recently, Shu and collaborators presented an integrative study where comparative analysis and computational prediction have been applied to assess the cross-species transportation of miRNAs, particularly focusing on inferring the likelihood of exogenous miRNA in human circulation [16].
This work demonstrated the data-driven computational analysis is highly promising approach to study novel molecular characteristics of deliverable miRNAs allowing to bypass the complex mechanistic details.
Several miRNAs discovery methods including computational prediction, cloning strategies and others have been used [17][18][19], even though these methods demand an increased rigor in miRNA annotations [20]. High-throughput sequencing technologies have contributed markedly to the expansion of knowledge about the miRNA universe in eukaryotic cells. These technologies have revealed a number of newly evolved and species-specific miRNAs that were previously unknown [21][22].
With the effort to define medicinal factors found in this interesting plant species, we searched for understanding better genome expression and gene regulation by micro RNA (miRNA) in M. oleifera Lam. We recognize that seeds contain all information about the tissues that will develop in the adult plant; therefore, we have programmed to analyze the miRNA populations obtained from seeds. In this work, we have identified novel miRNAs and their potential target genes from M. oleifera seeds using Illumina platform technologies. We found 94 conserved miRNAs and two novel ones, and some of them were validated by qRT-PCR. Target genes were predicted by in silico and results indicate that mol-miRNAs are putatively involved in many physiological process. A selected number of miRNAs were compared to other crop species plants, resulting to be conserved across multiple species. By taking advantage of a recently developed web-application based on an algorithm that compares plant and mammalian miRNAs (http://160.80. 35.140/MirCompare), we have identified a few M. oleifera miR-NAs with functional homologies to mammalian ones. We conducted a preliminary analysis to investigate potential human gene regulation by the plant miRNA mimics. The reported analysis increases the information on plant miRNAs currently available and improves the knowledge on to the molecular mechanisms associated to nutritional and medicinal activities of this plant species.

Plant materials, sample collection and RNA extraction
Seeds collected from mature pods, before they split open and fall to the ground, harvested from M. oleifera trees grown in the District of Dschang in West Cameroon, without use of chemical fertilizer, by the Cooperative of Medical Plant Producers SOCOPOMO. Seeds were stored at -80°C until used. Seeds were germinated in a greenhouse at the Department of Biology, University of "Tor Vergata", Rome, by placing them on paper soaked in sterile water. Leaves, stems and roots were collected at one month from the beginning of germination, immediately frozen in liquid nitrogen and subsequently stored at -80°C until used. Tissues from other species (Solanum tuberosa, Olea europaea and Medicago sativa) grown under natural field conditions, were also stored at -80°C until used. Total RNA was extracted from plant tissues using the mir-Vana kit (Ambion, USA) according to the manufacturer's protocol. The assessment of RNA quality and quantity were evaluated by Agilent 2100 Bioanalyzer (Agilent Technologies, Waldbroon, Germany) and by spectrophotometry (SmartSpec Plus, Bio-Rad, USA), respectively.

Constructing and sequencing small RNA libraries
Quality control for next generation sequencing experiments was performed by Genomix4life S. r.l. (Baronissi, Salerno, Italy). Individual indexed libraries were prepared from 1 μg of purified RNA using the TruSeq SmallRNA Sample Prep Kit (Illumina, USA) according to the manufacturer's instructions. Libraries were quantified using the Agilent 2100 Bioanalyzer (Agilent Technologies, Waldbroon, Germany) and pooled so that each index-tagged sample was present in equimolar amounts with final concentrations of the pooled samples adjusted to 2 nM. The pooled samples were subject to cluster generation and sequencing using an Illumina HiSeq 2500 System (Illumina, USA) in a 1x50 single read format at a final concentration of 10 pmol. The raw sequence files generated underwent quality control analyses using FastQC [23]. This method offers a powerful means for quantitative and qualitative profiling small RNA populations in different plant species, for which limited genome information is available, such as M. oleifera.

Sequence data analysis
Several bioinformatics tools have been developed to identify conserved miRNAs and discover new ones, starting from high-throughput sequencing [24][25][26][27]. The biggest limitation of these algorithms is that genome sequencing is required to execute the analysis. The only published algorithm able to discover novel miRNAs in species without sequenced genome is miReader [28]. Jha and coworkers also presented an approach for identifying novel miRNAs in Miscanthus giganteus, whose genome has not been sequenced yet [28]. Although the M. oleifera genome was already sequenced in 2015 [7], the raw data are not yet available, therefore we proceeded with a genome-independent strategy based on Jha's work. S1 Fig summarizes the analysis workflow performed with small RNA raw data. This analysis was performed to identify a pool of RNA fragments with lengths between 18 and 24 nucleotides that would be added to the subsets of known or uncharacterized miRNAs. These data have been deposited in NCBI/GEO public database under accession number GSE70423.
Raw data filtering. Starting from a total of 31,290,964 reads the Illumina small-RNA adapter was first clipped using FASTX-Toolkit (http://hannonlab.cshl.edu/fastx_toolkit). Sequences shorter than 15nt after trimming were discarded and we obtained a total of 22,737,895 reads. Sequencing artifacts were removed and the remaining reads were trimmed on the 3p and 5p ends in order to remove low-quality bases. At the end of this procedure BLASTn [29] v 2.2.30 with parameters "-task megablast, -perc_identity 100" was used to compare the remaining reads with the Rfam database (Rfam 11.0) [30]. This stringent step removed non-coding RNA (rRNA, tRNA, snRNA, snoRNA) and degraded fragments of mRNA. Reads with lengths less than 15 nucleotides were discarded. After the alignment process, a total of 10,415,180 reads were identified as possible miRNA sequences.
Identification of known and putative novel miRNAs. The 10,415,180 reads highlighted as possible miRNAs were aligned against a specific-plant miRBASE (Release 21) with a lowredundancy rate (http://www.mirbase.org/) using BLASTn v 2.2.30 [29] with parameters "-task blastn-short, -perc_identity 100". Sequences with homology rate equal to 100% were considered conserved M. oleifera miRNAs. In the alignment process all reads with abundance values lower than 10 copies were discarded resulting in a total of 303,872 reads.
In order to identify uncharacterized putative novel miRNAs, all the discarded reads, generated from the alignment with the wide libraries of miRNAs present in the repertory of miRBase bank, have been selected by a homemade custom script and were considered as candidate novel miRNAs.
Subsequently, the miReader algorithm [28] has been used to identify a strict number of uncharacterized miRNAs without the support of a reference genome. Although, this algorithm is less qualitatively informative (i.e. prediction of secondary structures) than others algorithms, it can be considered the most effective quantitative analysis with unknown genome, as is the case of M. oleifera.
Quality controls. Quality controls were performed using FastQC software [23], before (S1 File) and after (S1 File) removing the sequencing adapter. Figures A and B in S1 File show the quality scores across all bases and sequences in more detail, highlighting the validity of sequencing and analysis processes. The Illumina small-RNA adapter inside the sequenced sample is clearly visible. ( Figure C in S1 File). In Figures A and B in S2 File is visible how sequencing quality was preserved through the filtering process, the complete removal of the small-RNA adapter (Figure C in S2 File) and the removal of all ncRNAs that do not belong to miRNAs.
Predicting conserved miRNAs across other plant species. In order to evaluate the conservation rate for M. oleifera miRNAs across all others plant species, we developed a custom script, which retrieves the sequence information for plant organisms stored inside MirBase repository and evaluates the overlapping rate in terms of sequence homology between M. oleifera and all other plant miRNAs.

Target prediction
Predicting mRNA targets in plants for both known and uncharacterized miRNAs. S2 Fig shows the process to predict gene target for known and novel-putative miRNA. Starting from the list of 74 known, psRNATarget [31] was used to predict target genes in the plant. To reduce the false-positive rate we set the maximum expectation value to 2.0 and forced the system to select only the best five predictions for each miRNA. After the target prediction analysis, the three gene sets were submitted to PlantGSEA [32] to obtain a smaller number of gene clusters that are related to plant molecular processes.
Predicting known miRNA with putative roles in human gene regulation. Several Moringa miRNAs were further analyzed to investigate their possible roles in human gene regulation. MirCompare (http://160.80. 35.140/MirCompare) is a web-application based on an algorithm that compares libraries of miRNAs belonging to organisms from plant and animal kingdoms in order to find cross-kingdom regulation. The main settings for MirCompare analysis are: 1) r-value: the best comparison rate between each miRNA couple; 2) seed-region stringency: the minimum number of matches, related to the seed region.
MirCompare was used with r-value = 0.55 and seed-region threshold = 5 to identify possible Moringa miRNAs with functional homologies to mammalian miRNAs. The combinatorial miRNA target prediction (COMIR) web tool [33], which combines four popular scoring schemes (miRanda, PITA, TargetScan and mirSVR), was used to compute the potential of a gene to be targeted by a set of miRNAs. Therefore, COMIR tool generated a list of genes that might be regulated by M. oleifera miRNAs. Subsequently, the identified genes were submitted to Gene Ontology enrichment analysis to understand the pathways, in which miRNAs might putatively play a regulating role.
Transfection. We studied biological activity of the synthetic mol-miR168a (UCGCUU GGUGCAGGUCGGGAC) that was transfected into the hepatocellular carcinoma cell line (HEP-G2) by the lipofectamin method (Hi-Fect, Qiagen, USA) according with manufacture's instruction (MiRNA mimic and inhibitor experiments protocols1, Qiagen, USA). Synthetic mol-miR168a tagged with fluorescein isothiocyanate (FITC) was used also as control of transfection. Cells were harvested 72h after transfection and characterized for the efficiency of miRNA uptake by EVOS FLoid cell imaging station (LifeTechnologies, USA) analysis, and for the effect of miRNAs on specific target genes by Western Blot assay.
Western blot analysis. Aliquots of 3x10 6 transfected cells were suspended in buffers for cytoplasmic or nuclear protein extraction, and further processed for western blot analysis as already described [34]. Primary antibodies included the rabbit monoclonal antibody SIRT1 and the goat monoclonal antibody human beta-actin (Santa Cruz biotechnology, CA, USA). Secondary antibodies included anti-goat and anti-rabbit IgG chain-specific antibodies that were conjugated to peroxidase (Calbiochem, Merck Millipore, Darmstadt, Germany). Western blot analysis for each sample was quantified by densitometry analysis (TINA software) Data analysis was performed using GraphPad Prism version 6.00 (GraphPad Software, La Jolla, CA, USA). Statistical probabilities were expressed as p <0.05 ( Ã ) or p < 0.01 ( ÃÃ ). Comparison of means of SIRT-1 protein expression in response to mol-miR168a transfection was carried out using t-Test analysis.

Experimental validation of conserved miRNAs by Quantitative Real-Time PCR (qRT-PCR)
The expression level of different miRNAs was validated using poly(A)-tailed qRT-PCR method. According to the manufacturer's protocols (Exiqon A/S, Vedbaek, Denmark), a poly-A tail was added to the mature miRNA templates (20 ng). cDNA was synthesized using a poly-T primer with a 3' degenerate anchor and a 5' universal tag at 42°C for 60 min followed by heat-inactivated for 5 min at 95°C. To provide a control for quality of the cDNA synthesis reaction and the PCR, RNA spike-in (UniSp6) was added to the sample prior to cDNA synthesis. The cDNA template was then amplified using miRNA-specific and LNA™-enhanced forward and reverse primers. SYBR1Green was used for detection. The reactions were carried out in a Rotor-Gene1Q 72-Well Rotor (Qiagen, USA) with the following amplification conditions: activation/denaturation at 95°C for 10 min followed by 40 cycles of denaturation at 95°C for 10s, annealing and extension together at 60°C for 60 sec. Finally, melting analyses were performed to confirm the absence of false-positive peaks. qRT-PCR was performed only for conserved miRNAs. All reactions were performed in triplicate for each sample. Two controls (no template control and no Reverse Transcription control) were included in all reactions. Relative expression levels of miRNAs were quantified by using the 2 -Δ ΔCt method and miR159 was used as the internal control miRNA. To determine significant differences among samples or miR-NAs we applied a One Way ANOVA analysis using GraphPad Prism version 6.00.

Results
Analyzing small RNA populations sRNAs from M. oleifera seeds with 5'-phosphate and 3'-hydroxyl groups were identified by high-throughput Illumina TruSeq smallRNA sequencing (Illumina HiSeq 1500). cDNA libraries were constructed from seeds of M. oleifera plants after removing the 5 bp adapter sequence and filtering out low quality ''n" sequences. One small RNA library was constructed with about 30 million reads. Reads cleaned from adapters, ranging from 14 to 51 nts in length, were filtered with the Rfam database [30]. Looking at the sequence distribution after the filtering process, the majority of reads were 18 to 25 nts in length. The main size groups were 21 nt and 24 nt, respectively (Fig 1). These results were consistent with previous studies in other plant species where 24 nt small RNAs were the most abundant [19]. The amounts of 21 nt and 24 nt small RNAs were approximately 18.15% and 16.06%, respectively. Identifying conserved M. oleifera-specific miRNAs Conserved families of miRNAs are present in many plant species due to their important regulatory role played. The highly conserved nature gives the opportunity to identify miRNAs in plant species for which genome sequence information is partially or fully available, or previously undefined. To identify the conserved miRNAs in M. oleifera seed plants, the massive dataset was compared to known plant mature miRNAs against miRBASE (Release 21), using BLASTn. The sequence analysis revealed the presence of 94 miRNAs belonging to 40 conserved families (Table 1). Deeper investigations conducted by miRBASE also identify secondary structures for most of the known M. oleifera miRNA sequences (S1 Table).
The most abundant miRNA families were miR166, miR156 and miR167 with eight members each. Among the other families, in miR159, seven members, and in miR396 five members were present. The remaining miRNA families had less than five members. As reported in Fig 2, the sequence length distribution had two major peaks at 21 and 22nts with high levels of reads from 20 to 22 nucleotides. Low levels of 24 nucleotide fragments could be due to the high stringent filters applied.

Predicting uncharacterized M. oleifera miRNAs
All data set was analysed by BLASTn and 21,065 sequences were identified as available for novel miRNA identification with miReader software. To avoid spurious contamination of mRNA fragments or sequencing artifacts, this filtering process was forced to discard all reads shorter than 18nts nucleotides. The length distribution of the resulted sequences is shown in Fig 2. The distribution had a peak at 18 nts decreasing as the sequence length increases. The processed set of sequences was sub-sequentially provided to miReader and, after the analysis, only two uncharacterized miRNA duplexes were predicted as novel uncharacterized miRNAs ( Table 2). This result was not surprising since only the non-conserved miRNAs are normally expressed at low level, only in specific cell-types or under specific organ development and growth conditions [35].

Identifying known miRNA conservation rates among plant species
In order to evaluate the conservation rate for the M. oleifera known miRNAs within the plant kingdom an in-house ruby script was used. In Table 3, the top 20 conserved miRNAs together with their conservation rates and their abundance in M. oleifera seeds are reported.
Moreover, we drew a plot for better understanding the relationship between miRNA conservation rate and abundance in the seed tissues of Moringa. A logarithmic scale was used to obtain a linear plot since the range of abundance values was too high. As shown in Fig 3, mol-miR166i was the most conserved miRNA throughout plant species and the most abundant miRNA in Moringa seed tissues.
Validating expression of conserved miRNAs among the organ of M. oleifera and crop plants species To validate the representative expression patterns of the miRNAs, we performed qRT-PCR analysis of eight conserved miRNAs on different organ tissues of M. oleifera and on three crop plant species. The data collected demonstrated that the expression patterns were similar between the two analytical tools (Illumina sequencing and qRT-PCR) for six of the eight miR-NAs, whereas mol-miR156f-5p showed different expression pattern as detected between the two molecular tools (Fig 4). As illustrated in Fig 5, mol-miR168d-5p and mol-miR156f-5p were more abundant than other miRNAs in seed, leaf and root tissues. Vice versa, mol-miR164c showed higher expression levels in stem tissues. However, apart from mol-miR164c other miRNAs detected by qRT-PCR had higher expression in seeds compared to the other organ tissues. Similar results were obtained by assessing conserved mol-miR168d-5p, mol-miR166j and mol-miR156t in multiple plant tissues. As shown (Fig 6A) mol-miR168d-5p was the most abundant among wellexpressed miRNAs, except for M. sativa sprout and O. europaea leaf. Moreover, when comparing mol-miR156t expression levels in different tissues, including zygotic embryo, and in seeds treated in different way, we found that it was mainly expressed in dry seed tissues rather than other tissues (Fig 6B).

Predicting target genes for conserved and uncharacterized miRNAs
The system of prediction analysis was forced to select the five best predictions for each miRNA, and from the resulting data, 48 putative target genes were identified as potentially regulated by known miRNAs (Table 4).
Despite research of targets for all putative novel miRNA, only 3 target genes were selected (Fig 7) that involve mol-miR2p-5p, the most abundant uncharacterized miRNA highlighted by the miReader algorithm.
Enrichment analysis with gene ontology terms on both target gene sets highlighted the relevance of these genes for plant organisms. As shown in Fig 8, most genes regulated by known miRNAs are involved in biological pathways localized to root, meristem and seed.  Furthermore, these miRNAs participate in fundamental processes such as the development of anatomical (root, leaf) and reproductive structure, maintenance (flower, meristem, pollen and stamen), embryonic, and post-embryonic development. Genes putatively regulated by unknown miRNAs are involved in cellular metabolism of macromolecules, in particular macromolecules affecting protein transport and biosynthesis.

Bioinformatics prediction of human gene targets for M. oleifera miRNAs
To investigate whether M. oleifera miRNAs might regulate human gene expression we used MirCompare software to search homologous human miRNAs. The most and conserved mol-miR166i resulted functional homologies with hsa-miR6503-3p ( Table 5) that is involved in regulating inflammation [36], and mol-miR393c was homologous to hsa-miR548ah-5p that is involved in immune tolerance [37]. Further, mol-mir168a showed sequence homology with hsa-miR579, a human miRNA that normally regulates TNFα expression during endotoxin tolerance [38]. COMIR software [33] predicted human genes that could be regulated by M. oleifera miR-NAs. A subset of genes relevant to leukemia and acute myeloid (Table 6), which have key roles in the regulating cellular pathways such as apoptosis, cell cycle and protein degradation, have been found. Based on these results we identified mol-miR168a homologous to hsa-miR579, a human miRNA with many target genes including SIRT1. The complete report of predicted genes ranked by COMIR score (always upper than 0.9) is given in S2 Table.  Synthetic mol-miR168a transfection and protein modulation In these preliminary studies, we tested whether an exogenous miRNA derived from M. oleifera might be functional in human cell line and we verified its ability to inhibit post transcriptional mRNA expression. In particular, transfecting Hepatoma cell line G2 with the synthetic FITC mol-miR168a was run. By 72 hours after transfection with FITC ( Fig 9A, green cells, panel b), synthetic mol-miR168a increased the number of fluorescent cells by EVOS FLoid cell imaging station analysis (about 60-70% of HEPG2 cells resulted significantly and positively transfected). The transfection of synthetic mol-miR168a, a plant miRNA that show sequences homology with hsa-miR579, determined a significant decrease of SIRT1 protein level in comparison with HF and control samples (Fig 9B and 9C).

Discussion
Medicinal plants are studied to identify new therapeutic agents and understand their mechanisms of action against a variety of human diseases. Medicinal plant compounds are known to have high biocompatibility, low toxicity, and potential biological activity [39]. In particular, M. oleifera is commonly known and used for its health benefits [2]. For centuries and in many cultures around the world, M. oleifera has been used to treat human diseases [1] and is an example of traditional medicine that is increasingly popular among African countries, due in part to their poor economic conditions. A recent collaboration between universities in Italy and Cameroon has the objective of investigating the anti-oxidant and anti-tumor properties of M. oleifera [4]. The incredible Moringa's usage as medicinal, claimed by real-life experience, is now slowly being confirmed by scientific experiences.
Here we address a new aspect of plant biology that may have substantial impact on understanding of medicinal activity. The miRNA expressed in plants may comprise an independent category of medicinal agents capable of providing beneficial effects for human consumers of valuable plants.
Recently, Shahzad and collaborators [3] reported a study of DNA markers for genetic diversity and population structure in worldwide collections of M. oleifera. However, no study on the breadth of miRNA has yet been reported for M. oleifera, therefore, the knowledge about molecular mechanisms for compounds produced by M. oleifera is still limited. Here, we describe the identification and characterization of conserved and possibly novel miRNAs from seed of M. oleifera plants grown in Cameroon. We have used bioinformatics tools to understand their potential roles in diverse biological processes and their possible role in human gene regulation.  Since the discovery of sRNAs (miRNAs and siRNAs) as regulators of gene expression in C. elegans [40] and, more recently, the discovery of plant miRNAs [41][42], the study on miRNAs has become an important and integral topic in functional genomic research.
In plants, miRNAs regulate gene expression at the post-transcriptional level by degrading or repressing translation of target mRNAs [43]. A high number of experimental and computational studies have indicated that mature miRNAs are evolutionarily conserved in plants. miRNA-mediated gene regulation has an ancient phylogenetic origin and plays an important regulatory role in physiological processes [43], many aspects of plant growth, development and environmental adaptability.
High-throughput sequencing based on Illumina technology has become a good approach for identification and expression analysis of miRNAs in several plant species, like Arabidopsis [44] and other plant species [45][46], although the presence of biases introduced during the  construction of sRNA libraries, primarily derived from the adaptor ligation steps [47], requires a carefully attention.
In this study, we identified 94 conserved miRNAs belonging to 40 families, two uncharacterized miRNA duplexes with 48 and 3 targets gene for conserved and uncharacterized miRNAs respectively. As in earlier studies [45][46], the majority of highly conserved miRNAs from Moringa-seed that were predicted by our analysis resulted to be evolutionarily conserved across plant species and to have high levels of expression. For example, miR166, miR393, miR167, miR396, miR159 and miR156 families are well conserved among other plant species [48] and have a fundamental role in plant biology. Other miRNA families, such as miR408 or miR1515 were present at lower abundance. mol-miR166 families had the highest number of reads (80,612); in particular, mol-miR166i and mol-miR166 were the most and the least abundant miRNAs in these families, respectively. Results indicate that different members of the same miRNA family have differing levels of expression.
The most and least conserved miRNA families may have evolved to play different roles in plant biology. For example, miR156, miR159, miR166 and miR160 target SPL, MYB, PHV and ARF genes respectively. These transcription factors are important to plant growth and development. The miR159 was reported to target MYB101 and MYB33 transcription factors, which are positive regulators of ABA signaling during Arabidopsis seed germination. Indeed, miR159 may play a role in seed germination [49]. MiR156 plays crucial role in the control of juvenileto-adult transition in plants by targeting the SQUAMOSA PROMOTER BINDING PROTEIN LIKE (SPL) plant-specific transcription factors. SPLs affect processes like leaf development, shoot maturation, phase change and flowering in plants [50][51][52]. A recent study showed that miR156 regulates shoot regenerative capacity and a gradual increase in miR156 leads to a decline in shoot regenerative capacity for old plants [53]. Auxin response factors were found as a predicted target of miR160, while research studies assessed that this miRNA targets ARF10, ARF16 and ARF17 and regulates various aspects of plant development in Arabidopsis [54]. MiR166/165 is an example of well-studied plant miRNAs implicated in various aspects of plant development. The miR166/165 negatively regulates its targets Class III Homeodomain Leucine-Zipper (HD-ZIPIII) transcription factors that in turn regulate the polarity establishment in leaves and vasculature and radial patterning of root. The majority of HD-ZIPIII gene family members consisting of PHABULOSA (PHB), PHAVOLUTA (PHV), REVOLUTA (REV), ATHB8, and ATHB15, are conserved in several land plants including bryophytes, lycopods and seed plants [55][56][57][58]. MiR393 is a conserved miRNA family discovered in many plants [41]. In A. thaliana, four F-box genes TIR1 (TRANSPORT INHIBITOR RESPONSE PROTEIN1), AFB1, AFB2 and AFB3 (AUXIN SIGNALING F-BOX) were validated as miR393 targets [59]. In rice,   overexpression of miR393 negatively regulates mRNAs of TIR1 and AFB2 [60]. TIR1 and AFB2 interact with IAA (INDOLE-3-ACETIC ACID) proteins, probably releasing the activities of ARFs (AUXIN RESPONSE FACTORS) and increasing resistance to auxin. The change in auxin response consequently affects diverse aspects of plant growth and development, such as flag leaf inclination, primary root growth, crown root initiation and seed development. In this study, three possible novel miRNA targets were found: mol-miR2p-5p could affect target genes involved in macromolecule metabolism, in particular cellular protein transport and biosynthesis. In Arabidopsis thaliana, the ATPTR1 gene seems to be involved in long-distance transport of di-and tri-peptides during seed germination [61] while DUF579 affects xylan biosynthesis and modulates cell wall biosynthesis [62]. Further research on possible targets in other plants may provide important evidence to facilitate the understanding of these novel miRNA functions. Additional studies into novel M. oleifera miRNAs may shed light on their roles in M. oleifera biological processes.
The final phase of our bioinformatics analysis focused on the potential for human gene regulation by the most conserved M. oleifera miRNAs. The possibility of identifying plant miR-NAs able to regulate human genome expression may be highly important in future studies on the nutritional value and medical usage of food. The combined use of MirCompare and COMIR software on the massive collection of data has identified a small number of human genes that might be regulated by M. oleifera miRNAs. These gene targets include cell-cycle regulation and signaling through the p53 pathway; genes related to some classes of cancers including leukemia, acute myeloid and lipoma. For instance, Sirtuins have important roles in cell cycle, apoptosis, metabolic regulation and inflammation. The human genome encodes seven Sirtuin isoforms SIRT1-SIRT7 with varying intracellular distribution; a number of studies  reported evidence for their roles in a spectrum of disease like cancer, diabetes, obesity and neurodegenerative diseases [63]. Recent evidence suggests that genomic stability requires cooperation of p53 and SIRT1 [64]. Our transfection experiments showed that mol-miR168a identified by MirCompare-COMIR software inhibited translation of SIRT1 mRNA in cancer cells. We focused our attention on mol-miR168a because this is the first plant miRNA involved in cross-kingdom activity (osa-miR168a shown high sequence homology with mol-miR168a) [14] and it has functional homology with hsa-miR579, a putative regulator of the SIRT1 gene.
In conclusion, we have identified a population of Moringa-specific miRNAs that could help our understanding of the regulatory role of miRNAs in this plant. Our results demonstrate that the differentially expressed miRNAs and predictions for their target genes provides a basis for further understanding M. oleifera seed miRNAs and biological processes in which they are involved. Further studies are necessary to search for more miRNAs that are novel and to validate their targets by expression analysis during seed development stage.  Table. Complete prediction report of gene target in humans, using COMIR software. (DOCX) kingdom concept. We thank Ivana Matic and Alberto Calderone for the first editing of the draft and are fully grateful to David C Pauza, University of Maryland, for his suggestion and improving of the final manuscript editing.