In-silico based identification and functional analyses of miRNAs and their targets in Cowpea (Vigna unguiculata L.)

Cowpea (Vigna unguiculata L.) is an important leguminous plant and a good diet due to presence of carbohydrate and high protein contents. Currently, only few cowpea microRNAs (miRNAs) are reported. This study is intended to identify and functionally analyze new miRNAs and their targets in cowpea. An in-silico based homology search approach was applied and a total of 46 new miRNAs belonging to 45 families were identified and functionally annotated from the cowpea expressed sequence tags (ESTs). All these potential miRNAs are reported here for the first time in cowpea. The 46 new miRNAs were also observed with stable hairpin structures with minimum free energy, ranging from −10 to −132 kcal mol−1 with an average of −40 kcal mol−1. The length of new cowpea miRNAs are ranged from 18 to 26 nt with an average of 21 nt. The cowpea miRNA-vun-mir4414, is found as pre-miRNA cluster for the first time in cowpea. Furthermore, a set of 138 protein targets were also identified for these newly identified 46 cowpea miRNAs. These targets have significant role in various biological processes, like metabolism, transcription regulation as transcription factor, cell transport, signal transduction, growth & development and structural proteins. These findings are the significant basis to utilize and manage this important leguminous plant-cowpea for better nutritional properties and tolerance for biotic and abiotic stresses.


Introduction
MicroRNAs (miRNAs) are distinctive regulatory member of the small RNAs that regulate gene silencing at post-transcriptional level. Gene silencing by miRNAs is an important, advance and exciting area of present regulatory RNA research. They are endogenous, non-coding in nature and

Creation of single tone EST
The repeated ESTs from the same gene were eliminated and a single tone EST per miRNA was produced by using BLASTn program against the cowpea EST database with default parameters [25].

Elimination of coding sequences
The initial potential miRNA sequences of cowpea, predicted by the mature source miRNAs, were checked for protein coding. The FASTA format of initial potential sequences were subjected against protein database at NCBI using BLASTX with default parameter [26] and the protein coding sequences were removed.

Creation of hair-pen structures
The initial potential candidate cowpea miRNA sequences, confirming as non-protein coding nature, having 0-4 mismatches with the reference miRNAs and representing single tone gene were subjected to generate hair-pen or secondary structures. Publicly available Zuker folding algorithm http://www.bioinfo.rpi.edu/applications/mfold/rna/form1.cgi, known as MFOLD (version 3.6) [27] was used to predict the secondary structures. The MFOLD parameters were adjusted same as published by various researchers for the identification of miRNAs in various plant and animal species [7,8,28]. For physical scrutinizing, the hair-pen structures either showing the lowest free energy ≤−18 kcal mol −1 or less than or equal to the lowest free energy of the reference miRNAs were preferred. The Ambros et al. [29] threshold values were applied as reference to finalize the potential miRNAs in cowpea. The stem regions of the stem-loop structures were checked and confirmed for the mature sequences with either at least 16 or equal to the reference miRNAs base pairing involved in Watson-Crick or G/U base pairing between the mature miRNA and the opposite strand (miRNA*).

Convergence and phylogenetic analysis
The convergence and phylogenetic analysis was carried out for the one of conserved cowpea miRNA (vun-mir398). Simply, the vun-mir398, for its conserved behavior in different plant species was checked for convergence and phylogenetic investigation. The vun-mir398 alignment was created with Glycine max (gma), Nicotiana tabacum (nta) and Cucumis melo (cme) by the publicly accessible web logo: a sequence logo generator and ClustalW to produce cladogram tree using neighbor joining clustering method respectively. The results were saved.

Prediction of miRNAs targets
Dual schemes were used to predict the potential targets for cowpea miRNAs. In the first scheme, the newly identified cowpea miRNAs were subjected to psRNATarget (http://bioinfo3.noble.org/psRNATarget), with default parameters [30]. The cowpea miRNAs that not produced potential targets through psRNATarget, were subjected to the second scheme as described by Barozai [31]. Briefly, the cowpea mature miRNA sequences were subjected as queries through BLASTn program. The parameters were adjusted as, database: reference mRNA sequences (refseq_rnat); organism: Vigna unguiculata (taxid:4072) and Program Selection: highly similar sequences (megablast). The mRNA sequences showing ≥75% query coverage were selected and further subjected to RNA hybrid-a miRNA target prediction tool [32]. Only targets, confirming stringent seed site located at either positions 2-7 and/or 8-13 from the 5′ end of the miRNAs along with the supplementary site and having minimum free energy (MFE) ≤−20 kcal mol −1 were selected. For more stringency, these targets were subjected to the NTNU microRNA target prediction tool available at http://tare.medisin.ntnu.no/mirna_target/search#results, to confirm the RNA hybrid results. These predicted targets were further analyzed through Gene Ontology (GO) on AmiGO website.

The new cowpea miRNAs
In order to identify and characterize the potential miRNAs in cowpea, a comparative genomic approach was applied using bioinformatics tools. This is in agreement with the previous reports [8,28,31] that the homology based search by applying comparative genomics is a valid and logical approach to find interesting findings in plants at genomic level. The current study resulted a total of 46 new conserved miRNAs from the analyses of 187487 cowpea ESTs using bioinformatics tools ( . The vun-miR4414 family is observed as cluster pre-miRNA. Available miRNAs literature revealed that all these 46 miRNAs are profiled for the first time in cowpea. In the light of the empirical formula for biogenesis and expression of the miRNAs suggested by Ambros et al. [29], these miRNAs are considered as a valid candidate after justifying the criteria B, C and D. According to Ambros et al. [29] only the criterion D is enough for homologous sequences to validate as potential miRNAs in other species. The present study is in agreement with the other research groups [21,[33][34][35][36] where similarity based search by applying comparative genomics has produced novel and interesting findings in plants genomics.

Characterization of cowpea miRNAs
Characterization of newly identified candidate miRNAs is a set crucial step for their validation, as reported earlier [16,17,37]. The pre-miRNA length of the profiled cowpea miRNAs ranges from 46 to 381 nt with an average of 159 nt. The pre-miRNAs were further illustrated on the basis of their length ( Figure 1). The minimum folding free energy (MFE) of pre-miRNA is a vital and valid term of characterization. The newly identified potential cowpea pre-miRNAs have shown MFEs in range from −10 to −132 kcal mol −1 with an average of −40 kcal mol −1 as shown in Figure 2. The numbers of mismatches of mature sequences with their reference sequences were observed in a range of 0-4 with an average of three mismatches as categorized in Figure 3. These values are matched with the previously reported values in different plants [21,[37][38][39]. Mature miRNA sequences lengths were observed from 18 to 26 nt with an average of 21 nt as explained in Figure 4. These findings of mature sequences length are in agreement to prior published data in other plant species [16][17][18]36]. The 52% cowpea miRNAs sequences were found at 5′ arm, while 48% were at 3′ arm ( Figure 5(A),6). The GC content was found from 18 to 86% with an average of 42% as shown in Figure 7. Strand orientation is another important character for the generation of mature miRNAs transcripts. In this study, 24 mature miRNAs were found on minus strand while 22 were observed on plus strand of the transcripts ( Figure 8). The same results for plus and minus strand orientation of mature miRNAs are in agreement with the earlier research work [40]. The identified conserved cowpea miRNAs were also characterized on the basis of their organ of expression as presented in Figure 9. These findings are similar with the earlier reports [37] and suggesting organ dependent expression pattern of miRNAs in cowpea. The miRNA organ specific expression would be utilized to manage the organogenesis in cowpea. The secondary self-folded stem-loop structures of the cowpea pre-miRNAs are observed with at least 17 nucleotides engaged in Watson-Crick or G/U base pairing between the mature miRNA and the opposite arms (miRNAs*) in the stem region ( Figure 10). Except few where the reference miRNAs have also less base pairing and these precursors do not contain large internal loops or bulges. The mature miRNA sequences are observed in the double stranded stem region of the pre-miRNA secondary structures, as shown in Figure 5(A). Almost similar findings for various plant and animal species were reported by many researchers [16,17,20,37,41,42]. Furthermore, the newly identified cowpea miRNAs were also confirmed as non-protein coding nature by showing no significant similarity with known proteins. This validation strengthens the expressed nature for computationally identified miRNAs as non-coding RNAs. Similar results were observed in various research papers by many groups [16,43,44].

Cluster pre-miRNA gene in cowpea
In animals, a large number of miRNAs have been found in clusters and have been predicted to have similar expression profiles and functions [45]. The miRNA clusters have rarely been detected in plants. They were first reported by Jones-Rhoades and Bartel [46]. In this study, we also identified one pre-miRNA (mir4414) as cluster in cowpea having two mature miRNAs within Figure 5(B). On the basis of current available literature, this miRNA family (miR4414) was found for the first time in cowpea as a cluster.

Convergence and phylogenetic studies
The newly characterized cowpea miRNA vun-mir398, due to its conserved nature, was investigated for convergence and phylogeny. Simply, the cowpea miRNA vun-mir398 alignment and cladogram tree, using neighbour joining clustering method, were created with Glycine max (gma), Nicotiana tabacum (nta) and Cucumis melo (cme) by the publicly available Web-Logo, a sequence logo generator [47] and ClustalW, a multiple sequence alignment tool [48]. The cowpea miRNA vun-mir398 is observed in convergence with Glycine max (gma), Nicotiana tabacum (nta) and Cucumis melo (cme) as shown in Figure 11(A). The Phylogenetic cladogram tree, as illustrated in Figure 11(B), clearly showed that on the basis of sharing a more recent common ancestor the cowpea miRNA is more closely related to Glycine max (gma) than Nicotiana tabacum (nta) and Cucumis melo (cme). Zeng et al. [49] have also reported conserved nature in Euphorbiaceous plants.

The potential cowpea miRNAs targeted genes
Profiling the potential cowpea miRNAs targeted genes is a vital step for validation of the computationally identified miRNAs. A total of 138 targeted genes were predicted for the 46 potential cowpea miRNAs. The detail description is mentioned in Table 2. Different cowpea miRNAs targeting same proteins and vice versa were predicted here. This showed that one miRNA target more than one mRNAs and a single mRNA targets by many miRNAs [50]. The profiled targeted genes are categories as, 27% (37 of 138) are engaged in metabolism, 26% (36 of 138) are playing role as transcription factors, 11% (15 of 138) are involved in transport activities, 11% (15 of 138) are shown with stress related, and the rest are engaged in hypothetical protein, signal transduction, growth and development, structural proteins and diseases related. Almost all of these targets were already reported as miRNA targets in other plants [7,16,17].     Majority (27%) of the newly characterized cowpea miRNAs are observed to regulate the metabolic proteins. Such findings regarding metabolism related genes targeted by miRNAs are similar with the prior publications in plants and animals [28,43,44]. Pectin methylesterase (PME) is an important enzyme that acts on pectin, a major component of plant cell wall. PME catalyzes reactions according to the double-displacement mechanism [51]. In this study, the PME is predicted as a putative target for vun-miR1882. Thus the vun-miR1882 is a valuable resource to regulate cell wall. Another important enzyme ribulose-1,5-bisphosphate carboxylase (Rubisco) is a key enzyme in photosynthesis and photorespiration, where it catalyzes the fixation of CO 2 and O 2 , respectively. Due to its rate-limiting property in photosynthesis, it is the prime focus of improving the plant productivity [52]. The cowpea miRNA (vun-miR2657) is predicted to target this important enzyme which is the potential resource to modify Rubisco expression and ultimately plant productivity. ::::::::.:::: ::::::     The transcription factor myeloblastosis (MYB) is an important regulator of many developmental and physiological processes in plants. Ballester et al. [53], suggested that the MYB also plays a significant role in regulating the flavonoid pathway in plants. The newly identified cowpea miRNA family vun-834 is found to target the MYB transcription factors. Thus this miRNA is an important resource to fine tune the MYB regulation for the desirable traits in cowpea fruit. The transcription factor, zinc finger is believed to be involved in many biotic and abiotic stresses as responding gene to manage the plant under these stresses [54]. The same family of transcription factor is also reported to play a crucial role in plant development [55]. The newly identified cowpea miRNA families vun-miR834 and 4992 are found to target this zinc finger transcription factor family. These miRNAs are important resources to regulate the zinc finger family proteins for the betterment of cowpea under various biotic and abiotic stresses and fruit development.
Similarly 12% targeted genes by cowpea miRNAs are engaged in transport activities. ATPbinding cassette transporters comprise a highly conserved family of ATP-binding proteins that are involved in transporting of various molecules across plasma membrane. Here vun-miR1514 is identified to target ATP-binding cassette transporters. Such findings are in agreement with the other workers in the miRNA field [37,43].
Biotic and abiotic stresses like salinity, drought, temperature extremities, heavy metals, pathogen attacks, and pollution cause huge yield reductions in plants [56]. Naturally plants have various systems to protect themselves from these stresses that occur at various levels, i.e., at whole plant, tissue, cellular, sub-cellular, genetic and molecular levels [56][57][58][59][60]. Many studies suggest that plant miRNAs are involved in these stresses [9,17,61]. In this study identified miRNAs such as vun-miR1525, 2657 and 9748 also targeted heat shock proteins that expressed in response of heat stress. This suggests the role of these miRNAs during the heat stressed condition of plants. Similar findings were reported in switch grass [17].
Some miRNAs of cowpea were observed to target the protein functioning in the process of cell signal transduction. Almost similar findings were observed by many researchers in various organisms [42,43]. Protein kinases are key regulators of cell function and play crucial role in protein phosphorylation and dephosphorylation that are major signaling pathways induced by osmotic stress in higher plants. Similarly, SNF1 (sucrose non-fermenting-1) is an osmotic-stress-activated protein kinase in Arabidopsis thaliana that can significantly impact drought tolerance of Arabidopsis thaliana plants [62]. These two important proteins were targeted by cowpea miRNAs families, like vun-miR435, 2606, 2609 and 4392 respectively. Serine/threonine protein kinase (STPKs) is another protein kinase that is targeted by miRNA family (miR5241), act as sensors of environmental signals and regulate different developmental changes and also host pathogen interactions [63].
In this study, newly profiled cowpea miRNAs were also observed to target hypothetical proteins, growth and development, structural proteins and disease related proteins. Such findings were also published earlier [19,21,37].

Conclusion
The current study is resulted 46 new miRNAs and their 138 targeted genes in an important commercial plant cowpea. All these miRNAs are profiled for the first time in cowpea. These findings will serve as resources to fine tune cowpea plant at micro-molecular level. This will help us to enhance the production ability of cowpea against biotic and abiotic stress tolerance. Furthermore these miRNAs and their targets are also powerful functional genomic resources in the Kingdom plantae.