Genome wide in-silico miRNA and target network prediction from stress responsive Horsegram (Macrotyloma uniflorum) accessions.

Horsegram (Macrotyloma uniflorum (Lam.) Verdc.) is a drought hardy food and fodder legume of Indo-African continents with diverse germplasm sources demonstrating alternating mechanisms depicting contrasting adaptations to different climatic zones. Tissue specific expression of genes contributes substantially to location specific adaptations. Regulatory networks of such adaptive genes are elucidated for downstream translational research. MicroRNAs are small endogenous regulatory RNAs which alters the gene expression profiles at a particular time and type of tissue. Identification of such small regulatory RNAs in low moisture stress hardy crops can help in cross species transfer and validation confirming stress tolerance ability. This study outlined prediction of conserved miRNAs from transcriptome shotgun assembled sequences and EST sequences of horsegram. We could validate eight out of 15 of the identified miRNAs to demonstrate their role in deficit moisture stress tolerance mechanism of horsegram variety Paiyur1 with their target networks. The putative mumiRs were related to other food legumes indicating the presence of gene regulatory networks. Differential miRNA expression among drought specific tissues indicted the probable energy conservation mechanism. Targets were identified for functional characterization and regulatory network was constructed to find out the probable pathways of post-transcriptional regulation. The functional network revealed mechanism of biotic and abiotic stress tolerance, energy conservation and photoperiod responsiveness.

Introduction

MicroRNAs (miRNAs/ miRs) are 19–24 nucleotides long endogenous molecular bigwigs in post-transcriptional gene regulatory networks[1]. These miR genes are capped, polyadenylated and spliced like other RNA polymerase II transcripts. The mature miRs are located in a hairpin structure within the primary transcript (pri-transcript) and are preprocessed by at least two RNase mediated steps to mature miRNA. Pri-miRNAs range from 50–100 nucleotides that coil into hairpin loop structures encompassing paired stems and unpaired loops[2]. Intensive investigation over decades has enriched the miRNA’s knowhow in biogenesis and explicit regulatory machinery. Plant miRNA precursors are less conserved, whereas mature miRNAs are conserved at higher magnitude in comparison to animal miRNAs[3]. Hence, exploring conserved miRNAs makes sense to identify the potential miRNAs from flora with targeted traits.

Macrotyloma uniflorum (Lam.) Verdc. (Horsegram) originated from South-Western India is a multipurpose pulse crop mainly cultivated for providing nutritional security in the form of food, livestock supplement and green manure in India and Africa[4-6]. Being a diploid (2n = 20, 22 or 24], short duration (120–180 days to maturity) plant species and adapted to grow on wide range of agro-climatic conditions, horsegram can be weighed as appropriate model for moisture stress tolerance genes/QTLs investigation challenging undernourishment in drought prone regions[6,7]. Further, it can be envisaged as nutraceuticals, forage crop[8] and anti-calcifying inhibitors[9]. Recently, substantial improvement in accumulation of EST[5,10] and transcriptome data[11] has compelled its position as a future crop with versatile utility. The present study details the comprehensive computational approach[12-14] to predict miRNAs fromavailable ESTs and Transcriptome Shotgun Assembly (TSA) sequences of horsegram based on miRNA homolog search. Potential pathways contributing to drought tolerance is studied as inherent trait of horsegram, which could have common or divergent gene regulatory networks[6,10,15].

Flora without whole genome sequences are having alternate resource sequence sources in public databases such as Genome Survey Sequences (GSS), Expressed Sequence Tags (EST) and Bacterial Artificial Chromosome (BAC) sequences rendering plentiful resources to mine conserved microRNAs[16-18]. The database miRBase hosts 8746 reported miRNAs belonging to four phyla of plant species enlisted both mature and precursor miRNAs (https://www.miRbase.org/) (Release 21, June 2014). Plant miRs share functional similarity with small interfering RNAs (siRNAs) in guided target cleavage as microRNA targets sites of coding mRNA sequences[3,19]. Currently, the data generated from next generation sequencing (NGS) studies have been employed in miRNAs prediction and their impact on multiple traits in related species. The in-silico homology based prediction of miRNA is advantageous and orthologs of previously reported miRNAs could be deciphered with their evolutionary significance among species[20-22]. In recent years, the multi-faceted functionalities of miRs in plants are being understood effectively despite the fact that miRNAs are less studied in plants.

This is the first report in horsegram with comprehensive analysis of conservation and phylogeny of miRs. Previously, a comprehensive study was performed in horsegram from ESTs in predicting eight novel miRs[11]. Here we report differential expression of identified miRs through a stringent in-silico schema elucidating conserved miRNAs, their characterization, validation with their target annotation and networking.

Materials and methods

Query and reference datasets

Totalling 27,997 shotgun assembled contigs of transcriptome were collected from TSA sequence set: GANR01000001-GANR01027997 from European Nucleotide Archive (www.ebi.ac.uk/ena), SSH-Mu library of moisture stressed cDNA of Macrotyloma (https://www.ncbi.nlm.nih.gov/nucest NCBI dbEST ID 75866463 from) and 1008 ESTs were downloaded (www.ncbi.nlm.nih.gov/dbEST/) from publically available NCBI EST database to represent query sequences for miRNA homolog search. About 8496 miRBase mature Viridiplantae miRNAs were used (https://www.miRbase.org) (Released 21: June, 2014) database[23] in the present investigation and clustered by CD- HIT-EST, with threshold value of 100[24]. Out of 3777 clustered sequences, only non-redundant miRNAs were used as reference miRNAs for finding the homologs in M. uniflorum (Lam.) Verdc. candidate sequences to create a local nucleotide sequence database. To elucidate likely role of miRNAs involved in drought stress tolerance trait, the predicted microRNAs were screened for target genes listed in DroughtdB database (https://pgsb.helmholtz-muenchen.de/droughtdb/).

Bioinformatics tools employed

For conserved miRNA prediction of horsegram candidate sequences from Viridiplantae miRNAs reported in miRBase database, NCBI BLAST version 2.2.27[25] an alignment tool was used (https://ftp.ncbi.nih.gov/blast/executables/blast+). The representative sequences of all plant miRNAs were obtained after clustering with CD-HIT-EST with threshold value of 100[24]. MFOLD online tool[26] was (https://unafold.rna.albany.edu/?q=mfold) used for pre-miRNAs secondary structures prediction. The miRNA targets were deduced using psRNAtarget server[27]. Online circoletto[28] tool was used to illustrate the miRNAs and its target genes by circos plot[29]. The phylogenetic analysis of predicted miRNAs with their closely related miRNA families was performed with MEGA7[30].

Computational prediction of potential miRNA homologs

The prediction of miRNA homologs was performed using TSA and EST sequences of M. uniflorum available at NCBI and EBI databases. The plant miRNAs were obtained from NCBI was clustered and the representative sequences generated were employed as reference miRNA database using BLAST-2.2.27+. The query sequence consisting of TSA and EST of horsegram were subjected to nucleotide blast (blastn) against reference sequences of plant miRNAs with the following set parameters: (i) word match size-7; (ii) length of mature miRNA sequence ≥ 18 nt without gap; (iii) mismatch range- 0–2 (iv) e-value-0.1. The filtered sequences were utilized to find out the coding and non-coding candidate sequences by performing blastx against non-redundant protein database. The protein coding regions were removed whereas the non-coding regions were exploited for secondary structure prediction and validation.

Secondary structure prediction and validation

With maximum stringency only non-coding sequences were subjected to Zuker’s folding algorithm based secondary structure prediction in MFold[26]. A sliding window of 100nt from ~ 80nt upstream and ~ 80nt nucleotide downstream of the mature miRNA were set to find the precursors[33] with parameters set as: (a) Folding temperature- 37 °C, (b) Ionic conditions of 1 M NaCl without divalent ions, (c) linear RNA sequence, (d) percent sub-optimality number of 5 and (e) maximum interior/bulge loop size-30. As reported elsewhere, to validate the structures of pre-miRs adjusted minimal folding free energy (AMFE) and minimal folding free energy index (MFEI) were calculated[34] (Table 1).

Table 1

Characteristics of predicted precursor miRNAs of horsegram.

TSA ID	MFP	LT	LT	LP	AU%	GC%	A	C	G	U/T	MFE	AMFE	FEI
GANR01006328	miR482	355	355	119	57.99	42.01	38	28	22	31	− 42.01	− 31.34	0.74
GANR01007318	miR482	460	460	130	39.24	40.76	38	21	32	39	− 40.76	− 35.61	0.87
GANR01008465	miR156	609	609	100	50	50	31	19	19	31	−50	−39	0.78
GANR01009673	miR2673	775	775	118	53.39	46.61	28	17	38	35	− 46.61	− 26.77	0.57
GANR01016903	miR5653	2143	2143	119	52.11	47.89	36	25	32	26	− 47.89	− 28.49	0.59
GANR01019897	miR156	670	670	107	56.08	43.92	30	26	21	30	− 43.92	− 26.82	0.61
GANR01022354	miR1507	950	950	91	54.45	45.05	31	23	18	19	− 45.05	− 36.04	0.8
GANR01022920	miR168	515	515	122	41.81	58.19	19	36	35	32	− 58.19	− 48.85	0.83
GANR01023649	miR2118, miR482	644	644	106	60.38	39.62	36	19	23	28	− 39.62	29.05	0.7
GANR01024080	miR390	576	576	102	56.87	43.13	24	24	20	29	− 43.13	44.41	1.03

MFP, microRNA family/families present in stem-loop structure; LT, length of TSA; LP, length of precursor miRNA; MFE -minimal folding free energy (−Kcal/mol); AMFE, adjusted minimal folding (−Kcal/mol); MFEI, minimal folding free energy index.

The secondary structures were screened for[35]: (i) Less than three nucleotide substitutions within predicted mature miRNA to reference miR, (ii) The candidate sequence must fold into an appropriate and proper stem-loop hairpin secondary structure, (iii) The localization of mature miRNA must be in one arm of the stem-loop structure, (iv) Less than six mismatches between mature miRNA and its corresponding star sequence (miRNA*) and (v) The secondary structure of putative pre-miR must contain high negative MFE and MFEI levels. Putative horsegram miRNAs were designated based on the standard nomenclature system[23,36].

Phylogeny of horsegram miRNA

To enumerate the conserved nature and its phylogenetic relationship of miRs and its precursors were aligned with reported plant miRNAs using BLASTn with e-value, maximum mismatch and hits number set as 10, 3 and 5 respectively. The homologous precursor miRNAs were aligned to predicted miRNAs with Phylogeny.fr web tool (https://www.phylogeny.fr/index.cgi)[37] by integrating multiple sequence alignment tools, MUSCLE[38] and GBlocks[39] to refine the alignment; phylogenetic tree construction by PhyML[40] and Tree rendering by TreeDyn[41]. The alignment output is used in MEGA7[30] for interpretation of molecular clock to estimate interfamily evolution.

Phylogenetic analysis

The analysis involved 146 nucleotide sequences. All positions containing gaps and missing data were eliminated. Evolutionary analyses were also conducted in MEGA7[30].

Functional annotation of miRNA targets

As annotations are not available for horsegram, the closely related soybean was used as a reference to annotate miR targets of the predicted miRNAs. Putative identified query miRs were hit against horsegram mRNA sequences using psRNATarget tool[27]. The parameters set target prediction were (i) 9–11 nucleotide mismatch range for translational inhibition, (ii) maximum 4 mismatches without gap at complementary site, (iii) multiplicity of target sites -2, (iv) maximum expectation value -3 and (v) number of hits -10. The homologs of putative miRNAs and their targets were represented in Glycine max genome as circos plot to demonstrate multiple targets of identified putative miRNAs in horsegram[28]. In addition, to elucidate the role of identified miRs in drought stress response, specific miRNA target genes among drought tolerance contributing genes reported in Droughtdb (https://pgsb.helmholtz-muenchen.de/droughtdb/) were confirmed. The genes reported for horsegram are insufficient, thus, A. thaliana was used as reference organism to determine gene targets for the new miRNAs predicted. The identified putative miRNAs were used as query against the A. thaliana DFCI gene index (AGI) release 15 and A. thaliana TAIR10, cDNA, removed miRNA gene (release date 14th December 2010) using psRNATarget tool. The parameters set for the target gene prediction were as follows: (i) range of central mismatch for translational inhibition 9–11 nucleotide; (ii) a maximum of 64 mismatches without gaps at complementary site; (iii) multiplicity of target sites 2; (iv) maximum expectation value 3 and (v) Number of hits -5. Similarly, to understand the specific possibilities of putative horsegram miRNAs involved in drought stress, the drought genes were downloaded from Droughtdb in order to search for drought miRNA target genes. The visual representation of miRNA and its targets was represented by drawing circus plot to show the multiple targets of horsegram microRNAs predicted.

Functional annotation and metabolic pathway analysis

Mercator[42] AgriGO[43] and B2g[44] were used to determine the functional roles and to classify Gene Ontology (GO) terms into molecular functions, biological processes and cellular components. Further, corresponding pathways were mined through KASS server, which assigns KEGG Orthology (KO) terms and employs KEGG pathways for miRNA target genes. The simplified schema of workflow adopted in present investigation for miRNA prediction is given in Fig. 1.

Figure 1

Pipeline for prediction of miRNAs and its in-silico validation. From non-coding seqeuences and secondary structure analyses miRs were predicted in-silico as depicted in this pipeline.

Quantitative expression analysis of identified miRNAs

Stress responsive horsegram germplasm were identified from the germplasm core developed as described earlier (6). Selected stress responsive variety Paiyur 1 was raised in glass house under control and moisture stress conditions. Plants were maintained at field condition for control and at temporary wilting point for stress condition. Leaf tissues of control and stress plants were collected and fixed in RNALater and stored at − 80 °C till RNA isolation.

High quality RNA samples were extracted from 100 mg (wet weight) leaf tissue using Plant RNA Easy mini spin coloumn kit following the manufacturer’s guidelines (Qiagen). Quantity and quality of the RNA was quantified with the advanced Qiaexpert.

RNA samples of one microgram each were 3′ polyadenylated using Poly A RNA polymerase (Sigma Aldrich). Forward primers sequences of all 15 miRs and a common Poly T reverse primers were designed (Supplementary material 1).

Quantitative rtPCR using One Step PrimeScript™ RT-PCR Kit with SYBR green (Takara) was performed for all 15 identified miRs. qPCR capturing was conducted using the Illumina Eco RT PCR machine.

For qPCR cycling conditions set were 50 °C for 10 min for cDNA sysnthesis; 95 °C for 3 min for polymerase activation followed by 45° cycles of 95 °C 15 s, 60 °C for 30 s capturing with SYBR green filter at end of each cycle.

Results

miRNAs in M.uniflorum

Significant miRNA homologs within reported 8496 miRNAs were identified by executing nucleotide blast (BLASTn) with 27,997 TSA contigs, SSH-Mu library sequences of moisture stressed horsegram cDNA (NCBI dbEST ID 75866463) and 1008 EST sequences as query resulted in 16 ESTs and 6303 TSA hits (E value 0.1) for further analysis. Stringent filtering (mismatch < 3 and length > 17) was incorporated to narrow down the hits to 5807 (8 ESTs and 5799 TSA) sequences for putative candidate sequences identification with non-coding regions. Based on coding potential, 214 noncoding sequences (7 ESTs and 207 TSA) were filtered out to predict ten distinct pre-miRs coding for 15 conserved mature miRs (Tables 1 and 2) clustering into nine different miR-families (Figs. 2 and 3).

Table 2

Details of the predicted horsegram miRNAs (tentative names given to putative predicted miRNAs. Submitted to MiRBase).

Horsegram miRNA	miRNA homolog	Mature miRNA sequence	LM	NM	Loc	Strand
mun-miR482d-3p	gma-miR482d-3p	gguaugggagguguagggaaga	22	0	3′	−
mun-miR482b-5p	gma-miR482b-5p	uuccuucccaauccccccaua	21	0	5′	−
mun-miR482a-5p	gma-miR482a-5p	agaauuugugggaaugggcuga	22	0	5′	+
mun-miR482-5p	pvu-miR482-5p	ggaaugggcugauugggaagca	22	0	5′	+
mun-miR482a-3p	gma-miR482a-3p	ucuucccaauuccgcccauuccua	24	0	3′	+
mun-miR156r	gma-miR156r	ugcucucuaucuucugucag	20	0	5′	−
mun-miR2673a	mtr-miR2673a	ucuuccucuuccucuucc	18	0	3′	+
mun-miR5653	ath-miR5653	ugaguugaguugaguugag	19	1	3′	+
mun-miR156k	osa-miR156k	gacagaagagagagagcaca	20	0	5′	+
mun-miR1507a	gma-miR1507a	agacgauguauggaaugaga	20	0	3′	−
mun-miR168a-5p	ath-miR168a-5p	ucgcuuggugcaggucgggaa	21	0	5′	+
mun-miR2118	pvu-miR2118	aggaauggguggaaucggcaa	21	0	3′	−
mun-miR482a	sly-miR482a	aggaauggguggaauuggaaa	21	2	3′	−
mun-miR390b-5p	gma-miR390b-5p	gugcuaucccuccugagcuu	20	0	5′	−
mun-miR390a-5p	ath-miR390a-5p	gcgcuaucccuccugagcuu	20	1	5′	−

LM, Length of the mature miRNA.

Figure 2

Secondary structures of predicted miRNA precursors of horsegram using MFOLD tool[26].

Figure 3

Identified miRNAs with (a) length variation and (b) family size in horsegram. The differentially expressed miRs were found to have more members in horsegram indicating network of ncRNAs playing key role in regulating stress tolerance mechanism.

The predicted miRNAs were analyzed for various structural features to distinguish from other small RNAs such as tRNAs, rRNAs and mRNAs (Table 1 and Fig. 2). Most crucial characteristic feature of stable secondary structure is its minimal folding energy (MFE) ranging from 37.3 to 59.6 (−kcal/mol) for ten in-silico predicted premiRNAs. The MFE of the precursor miRNA was retained low to achieve thermodynamic stability[45]. Owing to their sequence length polymorphism, pre-miRNAs characterization was based only on MFEI (Minimum free energy index) following previous reports[22,46]. The MFEI ranges from 0.57 to 1.03 for ten in-silico predicted pre-miRNAs (Fig. 2). The (A + U) % of precursor horsegram miRNAs ranges from 39.24–60.38 satisfies the criteria included by Zhang et al.[18]. Nucleotide distribution (A = 31.6%, U = 32.31%, G = 26% and C = 23.8%) of the predicted pre-miR are heterogeneous as given Tables 1 and 3. On an average pre-miRs were of 111.4 bp length and mature miRs ranged from 18 to 24 (Fig. 3a) with mean of 20.73 (Table 3). Of predicted miRs, miR482 family had more members (six) representing its presence as one of the supreme molecule in the consortia of horsegram regulatory miRs. Similarly, miR156 and miR390 were found to have two members each whereas, other miRs have singlets from each family (Fig. 3b). Although, no miRs were identified from ESTs and SSH library validating them as coding sequences, the TSA sequences showed miRNA frequency of 1 in 1937 contigs (Table 4). The results of the predicted miRs and its targets in horsegram were statistically analyzed and summarized in Tables 3 and 4.

Table 3

Summary statistics for precursor miRNAs of M. uniflorum.

Parameter	Mean	Standard Deviation	Minimum	Maximum
LP	111.4	12.02	91	130
LM	20.733	1.437	18	24
AU%	52.232	6.86	39.24	60.38
GC%	45.718	25.82	39.62	58.19
A	31.6		19	38
U/T	32.31		58.1	19
G	26		18	38
C	23.8		17	36
MFE	− 31.06		− 37.3	− 59.6
AMFE	28.37		− 31.34	− 48.85
MFEI	0.752		0.57	1.03

LP, Length of the pre-cursor miRNA; LM, Length of the mature miRNA.

Table 4

Summary of the outcomes of the bioinformatics approach adopted to identify miRNAs in horsegram.

Database accessed	NCBI dbESTEBI
Candidate number of ESTs	1050
Candidate number of TSA	27,997
Total number of candidate sequences	29,047
Candidate number of miRNAs	8496
Number of contigs containing potential precursors	10
Number of microRNA families	8
Number of microRNAs predicted	15
Frequency of horsegram miRNA	1 miRNA per 1937(approx) sequences

Evolutionary relationship with soybean

Soybean (G. max L.), the closely related model legume crop with complete genome information cracked has been used to establish a comparative syntenic map of predicted horsegram miRNAs to determine the putative regions of homologous miR genes. From mapping results, except mun-miR-482a-3p and miR2673a, rest of the putative horsegram miRNAs have orthologs widespread in G. max genome. miR156r and miR156k had more number of orthologs. The mapped miRs are depicted in soybean genome (Fig. 4). Comparative genomics enabled us to infer miR gene function, which could enumerate future research focus contributed by horsegram miRNAs in related species and in distant crops as well.

Figure 4

Synteny mapping of putative horsegram miRs with soybean genome. Synteny map explains the conserved sequences across species and their cross species transferability.

Conservation and phylogenetic analysis

The pre-miR homologs were identified by performing blast of horsegram pre-miRNA against the miRBase database. The hits were filtered to retrieve only miRNAs of same family. In this study, high degree of conservation was exemplified by comparison of mun-miR168 to other plant precursor miRNAs (Figs. 5 and 6). It is evident that horsegram miRNAs share similarity with related legume species (Fig. 7). It is conclusively apparent that, horsegram miRNAs seem to have evolved at different rates in different time period similar to other plants.

Figure 5

Pre-miR sequence conservation blocks for miR168 family of horsegram with related plants. Non-coding RNAs are conserved across species with specific role in development, metabolism and energy conservation.

Figure 6

Phylogenetic analysis were done using MEGA7[30] for predicted miRNA precursors of horsegram with their closely related plant miRs. Nine clusters of miR precursors clustered with related plant species.

Figure 7

Phylogenetic tree of mun-miR families were constructed using PhyML[40] and TreeDyn[41]. miR family 1507 is having related species miR in a separate cluster and their role in stress tolerance is also unique.

Putative target genes, their functions and networks

Inferring the function of miRNA targets is crucial to significantly substantiate the functional role of miRNAs in gene expression and regulation. For predicted 15 horsegram miRs, 39 target genes were identified and classified into different groups based on their functional annotation. As revealed by Mercator annotation results (Fig. 8) of predicted miR targets, diverse processes ranging from RNA transcription regulation, protein posttranscriptional modifications, development, signalling, biotic and abiotic stress tolerance and glycolysis were being regulated by the identified 15 miRs. In addition, target genes also regulate cell wall degradation, hormone biosynthesis and redox components like ascorbate and glutathione synthesis (Table 5). The annotation of target genes using the Mercator tool categorized them into12 broad biochemical processes (Fig. 9). Differential expression of drought responsive miRNAs and gene networks Transcriptome analysis of Illumina sequence data from eight samples representing shoot and root tissues of contrasting horsegram genotypes M-191 (drought sensitive) and M- 249 (drought tolerant) (NCBI Bioproject PRJNA216977) was used to decipher the horsegram miRNA expression in root and shoot samples each under normal and drought stress condition respectively[11]. The miRNA abundance estimated from blast output indicates the behavioral miRNA expression (Fig. 10). The clustering simplifies the differential gene expression levels of horsegram miRNA and samples based on the similar expression. Network predicted and depicted as target- target interaction (Fig. 11) designates machineries of energy conservation confirming earlier hypothesis of structural compaction and energy conservation to survive under stress conditions[14,15]. This hypothesis may best fit and appropriate for other crops as well. To extend this hypothesis we identified target homologs in other plants also (Fig. 12).

Figure 8

Mercator annotation results. Functional groups of identified transcripts based on annotation results drawn using Mercator[42]. Predicted miR targets ranging from RNA transcription regulation, protein post transcriptional modifications, development, signalling, biotic and abiotic stress tolerance and glycolysis were being regulated by the identified 15 miRs.

Table 5

Functional annotation of horsegram miRNA targets.

miRNA_Acc	Target_Acc	Target_Description
mun-miR482d-3p	AT2G28380.1	DRB2/dsRNA-binding protein 2
	AT2G36470.1	Plant protein of unknown function (DUF868)
	AT4G06560.1	Transposable element gene
	AT3G30713.1	Transposable element gene
	AT5G65700.1	BAM1/Leucine-rich receptor-like protein kinase family protein
	AT5G65700.2	BAM1/Leucine-rich receptor-like protein kinase family protein
	AT4G06587.1	Transposable element gene
mun-miR482b-5p	AT1G78270.1	AtUGT85A4, UGT85A4/UDP-glucosyl transferase 85A4
	AT2G32700.5	LUH/LEUNIG_homolog
	AT2G32700.2	LUH/LEUNIG_homolog
	AT2G32700.4	LUH/LEUNIG_homolog
	AT2G32700.6	LUH/LEUNIG_homolog
	AT2G32700.1	LUH/LEUNIG_homolog
	AT2G32700.7	LUH/LEUNIG_homolog
	AT1G48550.1	Vacuolar protein sorting-associated protein 26
	AT1G48550.2	Vacuolar protein sorting-associated protein 26
	AT5G28650.1	WRKY74, ATWRKY74/WRKY DNA-binding protein 74
	AT3G58640.2	Mitogen activated protein kinase kinase kinase-related
	AT3G58640.1	Mitogen activated protein kinase kinase kinase-related
mun-miR482a-3p	AT1G72050.2	TFIIIA/transcription factor IIIA
	AT1G72050.1	TFIIIA/transcription factor IIIA
	AT4G03080.1	BSL1/BRI1 suppressor 1 (BSU1)-like 1
	AT5G12000.1	Protein kinase protein with adenine nucleotide alpha hydrolases-like domain
mun-miR482a-3p	AT1G22930.2	T-complex protein 11
mun-miR156r	AT2G16000.1	Transposable element gene
	AT3G45775.1	Transposable element gene
	AT2G43370.1	RNA-binding (RRM/RBD/RNP motifs) family protein
	AT5G50570.1	SPL13A, SPL13/Squamosa promoter-binding protein-like (SBP domain) transcription factor family protein
	AT5G50570.2	SPL13A, SPL13/Squamosa promoter-binding protein-like (SBP domain) transcription factor family protein
	AT2G42200.1	SPL9, AtSPL9/squamosa promoter binding protein-like 9
	AT5G50670.1	SPL13B, SPL13/Squamosa promoter-binding protein-like (SBP domain) transcription factor family protein
	AT3G57920.1	SPL15/squamosa promoter binding protein-like 15
	AT1G27370.1	SPL10/squamosa promoter binding protein-like 10
	AT1G69170.1	Squamosa promoter-binding protein-like (SBP domain) transcription factor family protein
	AT1G27370.2	SPL10/squamosa promoter binding protein-like 10
	AT5G43270.1	SPL2/squamosa promoter binding protein-like 2
	AT1G27370.4	SPL10/squamosa promoter binding protein-like 10
mun-miR1507a	AT2G15970.2	COR413-PM1 | cold regulated 413 plasma membrane 1
	AT4G10465.1	Heavy metal transport/detoxification superfamily protein
	AT2G15970.1	COR413-PM1, WCOR413, WCOR413-LIKE, ATCOR413-PM1, FL3-5A3, ATCYP19 |/cold regulated 413 plasma membrane 1
	AT4G06678.1	Transposable element gene
	AT5G35142.1	Transposable element gene
	AT5G44925.1	Transposable element gene
	AT3G21250.1	ATMRP6, MRP6, ABCC8 | multidrug resistance-associated protein 6 |
	AT3G21250.2	MRP6, ABCC8 | multidrug resistance-associated protein 6 |
	AT3G54510.2	Early-responsive to dehydration stress protein (ERD4)
	AT3G54510.1	Early-responsive to dehydration stress protein (ERD4)
mun-miR168a-5p	AT1G48410.3	AGO1/Stabilizer of iron transporter SufD/Polynucleotidyl transferase
	AT1G48410.1	AGO1/Stabilizer of iron transporter SufD/Polynucleotidyl transferase
	AT1G48410.2	AGO1/Stabilizer of iron transporter SufD/Polynucleotidyl transferase
mun-miR2118	AT2G21230.2	Basic-leucine zipper (bZIP) transcription factor family protein
	AT2G21230.1	Basic-leucine zipper (bZIP) transcription factor family protein
	AT2G21230.3	Basic-leucine zipper (bZIP) transcription factor family protein
	AT5G45050.2	TTR1, ATWRKY16, WRKY16/Disease resistance protein (TIR-NBS-LRR class)
	AT5G45050.1	TTR1, ATWRKY16, WRKY16/Disease resistance protein (TIR-NBS-LRR class)
	AT3G05680.1	EMB2016/embryo defective 2016
	AT3G05680.2	EMB2016/embryo defective 2016 |
	AT1G50120.1	Unknown function/Expressed
mun-miR482a	AT5G15850.1	COL1, ATCOL1/CONSTANS-like 1
	AT1G74600.1	Pentatricopeptide (PPR) repeat-containing protein
	AT5G15840.1	CO, FG/B-box type zinc finger protein with CCT domain
	AT5G15840.2	CO, FG/B-box type zinc finger protein with CCT domain
	AT2G21230.2	Basic-leucine zipper (bZIP) transcription factor family protein
	AT2G21230.1	Basic-leucine zipper (bZIP) transcription factor family protein
	AT2G21230.3	Basic-leucine zipper (bZIP) transcription factor family protein
	AT1G18300.1	atnudt4, NUDT4/nudix hydrolase homolog 4
mun-miR390b-5p	AT1G58050.1	RNA helicase family protein
	AT5G67610.1	Uncharacterized conserved protein (DUF2215)
	AT5G67610.2	Uncharacterized conserved protein (DUF2215)

Figure 9

Target groups based on functional groups identification from B2g[44]. (a) Cell cycle, cell division and basic cell process were the key functions predicted among the miR target sequences. (b) Cell intergrity and intact cell membrane were being targeted indicating the activation of degradation pathways under stress conditions.

Figure 10

Expression heatmap of putative miRNAs differential expression in horsegram. Expression analysis clearly indicated two clusters of miRs. There is enough variation among the

source tissue and gene families. The results of in-silico and real time quantitative expression pattern of miR are similar confirming prediction efficiency of the pipeline.

Figure 11

Predicted miRNA target genes function network. The function based protein target network drawn using sting server[82] elucidates the role of miRs identified in cell development and maintenance of minimal functions during stress conditions. Most of the reductions in functions identified were correlated to energy conservation mechanism.

Figure 12

Predicted miR target whole genome wide homologs across species were identified using OrthoVenn2[83]. Like the conserved miRs across species, the target range is also conserved in different plant species indicating existence of similar mechanism.

Discussion

Comprehensive studies on plant miRNAs endorse stage specificity and multitude of targets[47-49]. The non-coding miRNAs play a regulatory role in its target protein coding gene expression[50]. The miRNA multiplexes with RNA induced silencing complex (RISC) guiding the repression or cleavage of its target messenger RNA by seed nuclei base-pairing[3]. Complementarity between miRNAs and their target genes are high to regulate developmental processes, metabolism and stress responses[51,52]. Hence, there is enormous necessity to identify and validate miRNAs for further downstream applications in plants[22]. In horsegram like deficit moisture stress tolerant crop, major genes interaction network influence in expression of a tolerant phenotype. From contrasting stress responsive genotype data, and qPCR of identified miRNAs from a single stress responsive genotype under irrigated and deficit stress (Fig. 13), we could identify differential expression of miR genes. These predicted 15 miRs were clustered in to 11 different families and are conserved. To validated identified miRNAs qPCR was performed as reported earlier (53). Of the total eight miRs validated, two clusters were formed. Mun-miR 482 was found in both the clusters; whereas, mun-miR1507 was found in one cluster which distinguish the tested variety Paiyur1 for its ability to react for stress (Fig. 14). Of the 15 identified mun-miRs, miR 156, miR 171 and miR 390 were found to be differentially expressed in germinating seeds of halophyte Reaumuria soongorica under salt stress conditions (54). Legume genomes are rich in SNPs. Saturated map of SNPs were reported in legumes like vigna (55) and In pigeonpea (56). Of which highest haplotype desnsity of 0.7380 was reported for serine threonine kinase coding disease resistance gene (56) which was identified as an important miR target in the present ivestigation. Variations in these miR sequences are expected among germplasm to change their expression. This could be a reason for non expression of other genes. The conserved nature of maximum plant miRNAs bolstered the miRNA search in different plant species utilizing available genomic resources like ESTs, GSS and mRNA sequences[34,57-59] to a great extent and we followed the trend in identifying them in a drought tolerant crop. Conserved miRNAs in ginger, garlic, coffee and tea were identified[60-63]. For predicted pre-miRNAs, MFEI resolution for length variation extends from 0.57 to 1.03 in horsegram. Pre-miRs illustrate the absence of large internal loops/bulges and at least 20 nucleotides for Wobble base pairing (G/U base pairings) or Watson–Crick base pairing between the miRNA and the star sequence[62,64]. Estimated (A + U) % of pre-miRs range satisfies the predetermined criteria[18]. The putative genomic region of these non-coding sequences were determined and intergenic sequences were utilized for secondary structure prediction with the stringent filtering criteria to attain potential precursor miRNA[26]. Resultant secondary structures were inspected for precursor miRNA and the positioning of the mature miRNA within its stem-loop structure manually. The sequences with suitable secondary structure[66] were characterized for structure attributes and its possible existence as potential precursor of the predicted miRs. The transcription of miRNAs from sense and antisense strands of genes were already reported[67]. Our results stand by the possibility of the miRNA in both sense and antisense strands (Table 2) similar results were already reported in potato, tobacco and B. rapa[16,68,69]. The frequency of miRNA using expressed sequence tags is 1 in 1000 confirming previous reports[62]. The interaction of miRNAs with its target genes can enumerate evolutionary role of microRNAs[48,70,71]. Utilization of miRNAs in RNA interference (RNAi) mediated gene regulation has emerged as an important tool in novel traits engineering either by over expression of a miRNA or target genes by synthetic miRNAs. Gene silencing/knockout may help in understanding miRNA in plant responses to stress resulting increased productivity with improved nutritional value[72]. Cellular functions were not known for eight of the target proteins identified in the present investigation. Biotic/Abiotic stresses are critically affecting growth and development of plants. There are reports on contrasting mechanisms confirming grades of deficit moisture stress condition to control among different horsegram germplasm sources[6,15]. Among the identified miRNAs mun-miR156r, mun-miR156k, mun-miR482a, mun-miR390b-5p and mun-miR482a-5p were noticed to be involved in RNA processing, protein synthesis and modifications as well as plant development by altering the expression of their respective targets, whereas the genes encoding abiotic stress, biotic stress (NBS-LLR resistant class) and signalling associated proteins were targeted by mun-miR482a-5p, mun-miR482d-3p and mun-miR1507a respectively. Thus, there is concurrent substantiation that microRNAs can play pivotal role in crop improvement[73] and the miRNAs predicted from this investigation can be useful for future research. Closely related homologs of predicted putative miRs exhibited high degree of conservation as in mum-miR482 with related plant mature miRNAs (Fig. 6). The precursor sequences from miR482, miR390 and miR150 represented a strong candidate promulgating its necessary importance at post-transcriptional gene regulation in horsegram. The conserved nature of precursor and mature miRNAs has been reported in various plant groups from earlier studies[16,34,73,74]. The potential of miRNAs to bind corresponding target mRNA are comparable with their complementarity to degrade target mRNA. Therefore, to infer the contribution of microRNAs in cellular functions and regulatory gene networks, the miRNA target gene prediction is a crucial step[64]. The majority of the predicted miRNA targets in horsegram depict energy conservation mechanism which play important role in survival during adverse conditions. Target gene prediction reveals the regulation of multiple genes by single miRNA with different levels of regulation[22,47,75]. Also, the genes targeted belong to more than one gene family, which shows the multitude of miRNA functions in various metabolic processes (Table 4). We could spot 39 target genes of predicted miRs in horsegram and its homologs in other plant species. This indicates that, drought tolerance phenotype can be manipulated by adjusting the expression of identified miRs and their targets.

Figure 13

Control and stress treatments of horsegram plants. Paiyur1 plants were grown under glass house conditions. C. control plant pot was maintained with frequent recharging of soil moisture by irrigation; S. Plants grown under severe soil moisture deficit condition. Day time the plants suffered temporary wilting and regains during evening hours.

Figure 14

Quantitative rtPCR based differential expression of identified mun-miRs. Heatmaps were generated using Heatmapper[84]. Quantitative expression of miRs formed two clusters. Expression of miR1507 and miR390 families defines Paiyur1 phenotypic response to stress conditions. Down-regulation of other family miRs could be correlated to decline in cellular functions under stress conditions.

Conclusions

Great concern is being publicized recently to unravel various mechanisms involved in drought tolerance with emphasis to changing climatic conditions. Even though, there are extensive inquiries on miRNAs discovery and their functional prediction were completed, some of the non-model plants with considerable traits were not subjected to systematic investigation to elucidate contrasting mechanisms and their key players. Eight novel miRNAs were[76] predicted from horsegram. However, these non-validated miRNAs were predicted with low stringency and standard nomenclature was not followed. Inappropriately, the drought hardy horsegram miRNAs still remain unknown and lacks extensive validation. In this study, 15 conserved miRNAs belonging to nine different families were identified from EST and TSA sequences and validated in this first report with its tissue specific differential expression (Fig. 10). Mun-miR482a-5p and mun-miR482d-3p are miR482 family miRNAs reported to regulate stress response in soybean and apple[77,78]. The same trend is observed in our results confirming the earlier report. Mun-miR390b-5p is a miR390b family miRNA earlier found for protein degradation and post transcriptional modifications in Arabidopsis, maize and soybean[79-81]. Upregulation of miR390b and downregulation of 1507 family miRNA was observed (Fig. 14) in qPCR of Paiyur1 variety which is a stress responsive accession. The present investigation indirectly links stress tolerance to energy conservation as indicated by the network rather than direct response to stress conditions (Fig. 15). Hence, we confirm the interplay of miRNA-stress response-NBS-LLR class R-protein response[77] in energy conservation to survive the stress conditions. Additionally, the study reveals that identified miRNA regulated target genes have differential biological functions including cell wall degradation, hormone synthesis and synthesis of redox component like ascorbate and glutathione (Table 5). These findings can further help in translational research of related legumes and may result in selection of better germplasm for higher productivity and nutritional security. The network predicted from this investigation confirms the earlier hypotheses: structural compaction to overcome stress tolerance[14,15]. Altogether, the outcomes of the present investigation deliver clarity that, horsegram is a drought adapted crop and can be considered as a model crop for drought tolerance research. Thus, the predicted horsegram miRNAs may unravel the unique tolerance capability associated to its metabolic pathways and the present workflow represents simple and straightforward approach for the prediction and characterization of miRNAs in those plants for which genomes are yet to be splintered.

Figure 15

Control and stress treatments of horsegram plants and perspective conclusion of the predicted stress tolerance mechanism. Illustration was drawn with wet lab results and glass house grown plant sample analyses. Under stress conditions AAO expression defines the tolerance with declined total sugar and chlorophyll contents. Pith autolysis and shrunken stomata was observed under soil moisture deficit stress.

Supplementary Table.