Identification and Characterization of Wilt and Salt Stress-Responsive MicroRNAs in Chickpea through High-Throughput Sequencing

Chickpea (Cicer arietinum) is the second most widely grown legume worldwide and is the most important pulse crop in the Indian subcontinent. Chickpea productivity is adversely affected by a large number of biotic and abiotic stresses. MicroRNAs (miRNAs) have been implicated in the regulation of plant responses to several biotic and abiotic stresses. This study is the first attempt to identify chickpea miRNAs that are associated with biotic and abiotic stresses. The wilt infection that is caused by the fungus Fusarium oxysporum f.sp. ciceris is one of the major diseases severely affecting chickpea yields. Of late, increasing soil salinization has become a major problem in realizing these potential yields. Three chickpea libraries using fungal-infected, salt-treated and untreated seedlings were constructed and sequenced using next-generation sequencing technology. A total of 12,135,571 unique reads were obtained. In addition to 122 conserved miRNAs belonging to 25 different families, 59 novel miRNAs along with their star sequences were identified. Four legume-specific miRNAs, including miR5213, miR5232, miR2111 and miR2118, were found in all of the libraries. Poly(A)-based qRT-PCR (Quantitative real-time PCR) was used to validate eleven conserved and five novel miRNAs. miR530 was highly up regulated in response to fungal infection, which targets genes encoding zinc knuckle- and microtubule-associated proteins. Many miRNAs responded in a similar fashion under both biotic and abiotic stresses, indicating the existence of cross talk between the pathways that are involved in regulating these stresses. The potential target genes for the conserved and novel miRNAs were predicted based on sequence homologies. miR166 targets a HD-ZIPIII transcription factor and was validated by 5′ RLM-RACE. This study has identified several conserved and novel miRNAs in the chickpea that are associated with gene regulation following exposure to wilt and salt stress.


Introduction
MicroRNAs (miRNAs) are small, endogenous, non-coding RNAs that are present in animals, plants and some viruses. These RNAs participate in the regulation of target genes by binding to complementary mRNAs, resulting in either their cleavage or translational repression. miRNAs are involved in diverse processes in different organisms, including developmental timing in worms, cell death and fat metabolism in flies, hematopoiesis in mammals and leaf development, floral patterning and environmental stress responses in plants [1].
MIRNA genes are transcribed as independent transcriptional units by RNA polymerase II enzymes to generate primary miRNAs (pri-miRNAs). pri-miRNAs form imperfect folded structures that are processed by Dicer-like1 nuclease (a member of the RNase III endonuclease family) to precursor miRNAs (pre-miRNAs). The secondary structures of these precursors are well conserved in plants. The pre-miRNA contains a miRNA-star miRNA (miRNA*) intermediate duplex from which the miRNA* eventually is degraded. However, recent studies have revealed the higher accumulation of miRNA* under certain conditions in plants, indicating the probable role of miRNAs in modulating plant growth and development [2]. Mature miRNAs are 19 to 24 nucleotides (nt) in length and interact with an RNA-induced silencing complex (RISC) to cleave specific target mRNAs or inhibit their translation. This complementarity plays an important role in determining the fate of the mRNA. When the complementarity between the miRNA and mRNA is perfect or near perfect, the mRNA is cleaved; however, if there are many mismatches between them, translational repression occurs. There are also instances in which miRNAs and mRNAs with perfect complementarities lead to the repression of translation and not to the usual cleavage.
The first identified miRNAs were lin4 and let7 in Caenorhabditis elegans, which is a model nematode [3,4]. The first plant miRNAs were identified in Arabidopsis [5] and later in other plants. Currently, 7,321 mature miRNAs have been reported in 72 plant species (miRBase version 20) [6]. Among dicots, the maximum number of miRNAs occurs in the legume family (1,460), followed by Brassicaceae (863). Although the legume family has the best representation in terms of the number of miRNAs, chickpea is a notable omission from the list.
Chickpea (Cicer arietinum) is the world's second most widely grown legume and is cultivated in more than 40 countries. The Indian subcontinent is the principal chickpea-producing and consuming region, contributing almost 70% of the world's total production [7]. Chickpea seeds are a rich source of protein and starch for the human population and the records of chickpea cultivation date back to 6,000 BC. Globally, chickpea is grown on 11.5 million hectares (ha) to produce 10.4 million tons with an average yield of approximately 0.9 t/ha, which is far below its yield potential of 6 t/ha under optimal growth conditions [7]. The disparity between the actual and potential yields can be explained by large numbers of biotic and abiotic stresses that adversely affect its productivity. Among the biotic stresses, wilt infection that is caused by the fungus Fusarium oxysporum f.sp. ciceris is a major concern. Abiotic stress conditions, such as terminal drought and salt stress, also lead to major losses. ICC4958 is a drought tolerant chickpea cultivar and gets affected at terminal drought, which occurs at the pod filling and seed-developing stage of the crop [8,9]. However, recent studies on salinity tolerance and ion accumulation in chickpea have revealed it as a highly sensitive crop to salinity when compared to other species in cropping systems [10,11]. Thus, salinity is another major constraint in chickpea yield. A better understanding of genes and their interactions with the environment can play a very important and determinant role in tackling these stress conditions. The recently available transcriptome and genome sequences that have been reported by independent groups are important resources that will facilitate the attainment of these goals in the chickpea [12,13,14]. Hu et al. (2013) identified 28 potential miRNA candidates belonging to 20 families from 16 ESTs and 12 GSSs in the chickpea using a comparative genome-based computational analysis [15]. A total of 664 miRNA targets were predicted, including genes encoding transcription factors (TFs) in addition to those that function in the stress response, signal transduction, methylation and a variety of other metabolic processes. These findings lay the foundation for the elucidation of miRNA function in the development and stress responses of the chickpea. miRNAs have been discovered primarily using direct cloning and bioinformatic approaches. All of the miRNAs in plants have been identified via the cloning of small RNAs or a computational approach, in which the homologs of known miRNAs are searched.
We have generated small RNA libraries corresponding with the control conditions, Fusarium wilt infection and salt stress, which were sequenced using the Illumina sequencing platform to identify miRNAs in the chickpea. This study is the first report in which small RNA libraries have been constructed and sequenced to identify miRNAs in the chickpea.

Sequence analyses
Three separate small RNA libraries that were constructed from the total RNA of the control, Fusarium wilt-infected and saltstressed plants were subjected to Illumina Solexa sequencing. This sequencing generated 29,170,463 raw reads, which after processing by UEA sRNA workbench 2.4-Plant version sequence file preprocessing tool (http://srna-tools.cmp.uea.ac.uk/), produced approximately 12,135,571 total unique reads (Table 1). After removing the adaptor sequences, filtering the low-quality tags and eliminating the t/rRNA sequences, the putative small RNA population accounted for approximately 88.5%, 79.1% and 78.4% in the control, wilt-infected and salt-stressed libraries, respectively ( Figure S1). The majority of small RNAs (approximately 50%) from the control and salt-stressed libraries were 24 nt in length (Figure 1), which is similar to other plant species, such as Arabidopsis thaliana, Solanum lycopersicum and Medicago truncatula [16,17,18]. Notably, in the wilt-infected library, small RNAs that were 20 nt in length accounted for 20% of the population, but when unique reads were analyzed, the small RNA distribution revealed that 24 nt was the major size class. Similar patterns have been reported in cucumber and soybean [19,20]. In soybean, the unique and redundant sequence classes possessed 24 nt and 23 nt long small RNAs, respectively, as the most abundant sequences. For the differential expression analysis, the total numbers of miRNA reads in each given sample were normalized as transcripts per million, and the fold changes between the treated and control samples were calculated. Out of 122 conserved miRNAs, 44 were upregulated in response to wilt stress, but in the case of salt stress, the number of down regulated miRNAs was greater than that which was observed in response to wilt stress. However, the differential expression of novel miRNAs under both of these stresses showed relatively similar patterns, with approximately 60% of the miRNAs being down regulated under either wilt and/or salt stress (Figure 2a, b).

Identification of conserved miRNAs in chickpea
The unique reads that were obtained from the miRCat analysis tool (UEA small RNA workbench) were mapped to the miRNAs that were available in miRBase version 18 (http://www.mirbase. org/) [21,22,23]. The small RNA sequences that matched the known miRNAs from the miRBase database were identified as conserved miRNAs in the chickpea. The sequence analyses revealed the presence of 122 miRNAs belonging to 25 conserved families. The most abundant family was miR156 with 14 members. Among the others, miR171 (12 members), miR169 and miR172 (9 members each), miR166 and miR167 (8 members each), miR319 and miR399 (6 members each) and miR396 (5 members) were present. The remaining miRNA families had less than five members, with some families, such as miR530, miR162, miR5232 and miR408, being represented by only one member ( Table 2; Table S1). In a recent report of the chickpea genome, the sequences of 20 unique miRNA families were reported, of which MIR169_2 and MIR159 were the most abundant [14]. In our study, out of 25 conserved miRNAs families, 16 possessed miRNA* sequences, thus providing additional evidence in support of the authenticity of the miRNAs. However, no miRNA* sequences were obtained for miR2111, miR162, miR164, miR390, miR394, miR397, miR530, miR408 or miR5213. Detailed information regarding the precursor structures of the conserved miRNAs is provided in Table S2.

Identification of legume-specific miRNAs in chickpea library
We obtained four legume-specific miRNAs, including miR2111, miR2118, miR5213 and miR5232, in our libraries that were previously reported in another legume, Medicago [24,25]. To date, miR5232 has only been reported in Medicago in a study involving miRNA regulation during arbuscular mycorrhizal symbiosis [25]. Accordingly, miR5232 may be a legume-specific miRNA that is involved in the biotic stress response. The multiple sequence alignment of the mature miRNAs in addition to the precursor sequences of these four legume-specific miRNAs revealed that they were most closely similar to Medicago and consequently has been conserved throughout evolution (Figure 3a, 3b). However, in recent studies, sequences that are similar to miR2118 have been reported in other non-leguminous plant systems, such as tomato and rice [26]. Apart from Fabaceae, miR2118 family members are most abundant in the Rutaceae and Solanaceae plant families [27]. Even the nomenclature of the miR2118 family is inconsistent in miRBase: the miR2118-like sequences have been disparately named miR482 (sly-miR482), miR5300 (as in tomato) and miR2809. The variation in the miR2118 sequence is species specific. Thus, miR2118 sequences in the chickpea are more similar with mtr-miR2118a [24] in comparison with other plant systems.

Identification of novel miRNAs in chickpea
We identified 59 novel miRNAs using the miRCat module of the UEA sRNA workbench, which aligned the pooled reads from all three of the libraries to the chickpea genome (NCBI Genome: PRJNA175619) [13], the ESTs database from NCBI and transcriptome data from the chickpea transcriptome database [28], and applied prediction criteria for plant miRNAs [29] ( Table 3; Table S3). The low abundance of novel miRNAs in our data supports the earlier notion of the lower expression levels of novel miRNAs compared with those of conserved miRNAs [30]. The precursor miRNA candidates were then tested using RandFold with a cutoff of 0.1. The minimum free energy that was required to form the predicted hairpin structure for the precursor was in the range of 297.2 to 226.03 Kcal/mol, which is similar to the values that were reported for the precursors of other plant species (Table S4). The secondary structures of the precursors of five validated novel miRNAs were evaluated using the Mfold software ( Figure 4) [31]. The data analysis revealed the presence of miRNA* sequences for all of the 59 novel miRNAs of the chickpea. The miRNA* supports the release of the miRNA duplex from the predicted hairpin structure [32]; therefore, the presence of miRNA* sequences further supports the identity of these small RNA sequences in our libraries as novel miRNAs.

Expression patterns of known and novel miRNAs in chickpea
Total RNA from the tissues of control, wilt-infected and saltstressed plants were used to validate the miRNAs. The poly (A)RNA of these three samples was reverse transcribed into cDNA for the validation of the expression of eleven conserved and five        novel miRNAs using qRT-PCR. The expression levels of the chickpea miRNAs under wilt stress were significantly altered compared with those of the control conditions. In contrast, those that were observed under salt stress did not greatly change. Among the validated conserved miRNAs, miR530 was upregulated seventeen-fold during wilt stress, suggesting that it is an important candidate miRNA that is involved in the plant wilt stress response. miR156_1 and miR156_10 were slightly upregulated under both the wilt and salt stresses. miR2118, which is one of the legume-specific miRNAs, was also upregulated by approximately 0.5-fold during wilt stress compared with the control seedlings (Figure 5a, 5b). Conversely, no significant expression patterns were detected with respect to the known miRNAs in response to salt stress in the chickpea. The expression analysis of the novel miRNAs revealed that three out of five (car-miR008, car-miR011 and car-miR015) were approximately three fold upregulated on average during salt stress (novel chickpea miRNAs have been designated as ''car-miRNA'' throughout the manuscript, in which ''car'' is an abbreviation for Cicer arietinum). However, the expression patterns that were observed during wilt stress revealed limited information because little significant changes occurred.

Prediction and validation of miRNA targets in chickpea
The putative miRNA targets in chickpea were predicted using the psRNATarget program [33]. The predicted target genes (approximately 358 different transcripts) were extensively involved in different biological processes involving a large number of gene families. Some of these genes encoded TFs, DNA replication proteins and those that are involved in cellular metabolism in addition to a variety of stress response-associated proteins. The target prediction analysis revealed the involvement of some of the miRNAs in regulating metabolic processes through the target genes. In chickpea, miR159 is involved in the metabolism of amino acids, fatty acids and lipids. One of the target genes of miR159 encodes acyltransferase, which is essential for ester biosynthesis. The chickpea miR156 and miR166 target genes encode squamosa promoter-binding protein and homeoboxleucine zipper protein, respectively, as previously reported [34,35]. Table 4 describes details of the target genes of validated miRNAs; a complete list is provided as supporting information (Table S5; Table S6). The most widely targeted class of genes is the protein kinases, which are associated with plant defense mechanisms via cell signaling-related processes. The novel car-miRNA008 targets the chalcone synthase (CHS) gene. Chalcone, which is an intermediate in flavonoid biosynthesis, is involved in natural defense mechanisms. CHS expression is also involved in salicylic acid defense pathways. car-miR2118 and car-miR5213 target two defense-response chickpea genes encoding Toll/ Interleukin-1 receptor-nucleotide binding site-leucine-rich repeats (TIR-NBS-LRR). Members of the TIR-NBS-LRR gene family are genuine targets for miR2118 [24]. Additionally, miR5213 suppresses defense-response genes in Medicago. The cleavage of such transcripts as mediated by miR5213 is notably conserved in AM symbiosis-capable plants, such as Medicago truncatula, Glycine max, Lotus japonicus, Populus trichocarpa and Cicer arietinum, but not in plants for which this symbiosis is not observed, such as Arabidopsis thaliana [25].
Members of the miR166/165 family target HD-ZIP III TF genes by cleaving the mRNA at complementary base pairs in leguminous plants [34,36,37]. These results are similar to those of earlier predicted reports involving other plant systems. The target gene of miR166 was experimentally validated by modified 59RLM-RACE [38,39]. All of the positive clones were sequenced, Table 2. Cont. and cleavage was observed at the 17 th and 18 th positions of the mRNA by the 59 end of miR166 (Figure 6), unlike the previously reported miRNA-target recognition parameters [40]. Although our results are not in agreement with previous studies, such as those involving the soybean, in which miR166 target validation by 59RACE and degradome sequencing confirmed cleavage at the  10 th and 11 th positions [41], there have been reports of the miRNA (belonging to different families)-mediated cleavage of target mRNA, thus defying the recognition rule. A total of 18 miRNA/target pairs of Pinus taeda possessed non-conventional cleavage sites, such as pta-miR951:AW065026, which is cleaved at the 16 th and 17 th positions [42]. Similar results have been reported in other plant species, such as mtr-miR397:AC135467 [24], ath-miR168:AGO1 [43], pvu-miR171:gi62704692 [44] and ath-miR398a:CSD1 [39]. Thus, it appears that the sequence of the target gene and the miRNA sequence determine the cleavage site apart from the conventional complimentary region-based target cleavage. Therefore, it is quite possible that chickpea has a different cleavage site for miR166 (pair miR166:TC04758) compared with other plant species.

Analyses of GO terms and KEGG pathways
The GO terms of the target genes were annotated according to their biological processes, molecular functions or involvement as cellular components. The enzyme mapping of the annotated sequences was performed using direct GO for the enzyme mapping and the Kyoto Encyclopedia of Genes and Genomes (KEGG) for the definitions of the KEGG orthologs. The miRNAtargeted genes belonged to various biological processes, cellular components and molecular functions as depicted in Figure 7. The maximum numbers of target genes were involved in biological processes, including both metabolic and cellular processes.
However, the target genes that were involved in binding were the most abundant (80%) within the molecular functions category.

Discussion
In this study, high-throughput deep sequencing was used to gain in-depth knowledge of gene regulation by miRNAs in the chickpea under biotic and abiotic stresses. Salt stress is one of the major constraints to increasing chickpea productivity. Soil salinity levels affect germination in plants. Under salt stress conditions, chickpea plants show high levels of anthocyanin pigmentation in their foliage and reduced growth rates [45]. Among the biotic stresses, Fusarium wilt is one of the major soil/seed-borne diseases severely affecting chickpea growth. Its causative agent is Fusarium oxysporum f.sp. ciceris, which is a fungal pathogen.
Most of the miRNAs that were obtained in our library have a preference for the 59-U as has been reported in other plants, which is in accordance with the defined structures of the mature miRNAs [1,46]. The lengths of the chickpea precursors ranged from 61 to 220 nt, which were similar to those of the soybean (55 to 239 nt) and peanut (75 to 343 nt) [30,35]. The calculation of the minimum free energy (MFE) values further added credence to these predicted hairpin structures as putative miRNA precursors. The chickpea precursors had minimum free energy values ranging from 297.2 to 226.03 Kcal/mol with an average of 250.1419 Kcal/mol, which was similar to the 250.01 Kcal/ mol that was observed in Arachis hypogaea and the reported value  (Table S2) [30]. Greater increases in miRNA expression were observed following wilt stress compared with salt stress, suggesting the significant role of small RNAs in the response to pathogen attack. The total number of miRNAs was greater in the wilt stress library than in the salt stress library. Four legume-specific miRNAs were identified in the chickpea libraries, including miR2111, miR2118, miR5232 and miR5213, which were previously reported in Medicago. The sequence conservation among the different legumes and the precursor sequence similarity of these four chickpea miRNAs further substantiate their accurate identification in this study. car-miR5232 cleaves only two transcripts encoding an ATPase E1-E2 type and an expressed protein of unknown function, in concordance with a similar study in Medicago, in which miR5232 targets were experimentally confirmed by degradome sequencing [25].
The significance of miRNA* in authenticating the presence of miRNA has previously been established. A comparison of chickpea miRNA* and mature miRNA data revealed that they vary in abundance in response to the different stress treatments, which has also been previously reported [47,48]. Our target search analysis indicated that miRNA* act upon different transcripts than do their parental miRNAs (data not shown), which has been observed in plants, animals and humans [25,49,50]. For example, miR393 and its miRNA* counterpart regulated the expression of genes belonging to entirely different protein families; i.e., TIR1 and SNARE, respectively [51,52].

Expression patterns during biotic stress
This study is the first attempt to identify miRNAs that are associated with fungal attack in the chickpea. Alterations in the expression of genes that are involved in defense during pathogen attack have been previously reported. These genes are regulated by small RNAs. miR393 was the first miRNA whose role in pathogen attack was demonstrated [51]. Eleven conserved and five novel miRNAs were analyzed in the chickpea under wilt and salt stress. Interestingly, miR530 was significantly upregulated during wilt stress. This observation suggests that its target genes are expressed at lower levels, which included the zinc knuckle proteins and microtubule-associated proteins. Zinc knuckle proteins are involved in the regulation of morning-specific growth in Arabidopsis [53]. The target of miR530 varies in different plants under different conditions and tissues. In Populus trichocarpa, this miRNA targets zinc knuckle (CCHC type) family proteins along with a homeobox TF [54], whereas in soybean, it targets genes that encode the CONSTANS interacting protein and nuclear transcription factor Y [55]. In Eugenia uniflora, miR530 targets wall-associated receptor kinase-like 14, S-acyltransferase tip-1 and a protein of unknown function in rice [56,57]. In a recent study in maize plants that were resistant to the fungus Exserohilum turcicum, miR530 was identified as a novel miRNA and was predicted to target genes that are involved in kinase activities in addition to DNA-binding TFs [58]. Based on the significant upregulation of miR530 in response to Fusarium infection and its unique target genes in the chickpea, it appears to be involved in the response to pathogen attack.
The three legume-specific miRNAs (miR2111, miR2118 and miR5213) play critical roles during pathogen attack. In the chickpea, miR2111 targets a Kelch repeat-containing F-box protein. F-box proteins are responsible for the controlled ubiquitin-dependent degradation of cellular regulatory proteins and are involved in defense responses, auxin responses and floral organ development [59,60,61]. Targets of F-box proteins are central regulators of key cellular events and include G1 cyclins and inhibitors of cyclin-dependent kinases [62]. It appears that miR2111 and F-box proteins act together to regulate the defense response in chickpea following biotic stress. Other than F-box proteins, miR2111 also targets TIR domain-containing NBS-LRR disease resistance proteins. miR2118 and miR5213 also target the same class of R genes. TIR, which is an F-box protein, is a receptor for the plant hormone auxin [63,64,65,51], and LRR consists of tandem Kelch repeats [66]. Interestingly, the chickpea miR2118 was upregulated in response to wilt infection and down regulated following salt stress. miR2118 has also been shown to be suppressed after Verticillium fungal attack in cotton [67]. Fusarium wilt leads to symptoms that are similar to those of Verticillium wilt, whose common host plant is cotton. miR2118 functions through three novel target transcripts encoding TIR-NBS-LRR disease resistance proteins, but its functional regulation remains unclear. In the soybean, miR2118 targets the protein family that is associated with disease resistance in addition to zinc finger proteins [55] and replication termination factor 2 in response to biotic (Asian soybean rust) and abiotic (water deficiency) stresses [35].
Other miRNAs also target disease resistance genes. For example, novel car-miRNA023 target proteins are involved in disease resistance. The highly conserved miRNA171 family targets more than 20 genes that are involved in different processes and Table 4. Predicted target genes of miRNAs in chickpea. miRNA family Target  Putative Functions of Predicted Targets   Conserved miRNAs   mir156_1  TC12891, TC03863, TC05745,TC07318,  Squamosa promoter-binding TF family protein, SCP1-like   TC03684, TC29077, TC07318, TC15422,  small phosphatase   TC04572   mir156_10  TC29077, TC15422, TC12891, TC03863,  SCP1-like small phosphatase, squamosa promoter-TC05745, TC03684, TC07318, TC19303,    pathways in the chickpea. One particular member, miR171_7, targets a disease resistance-responsive dirigent-like protein (DIR). The conspicuous involvement of disease resistance genes in the response to pathogen attack has been previously established. ESTs encoding dirigent proteins were identified in the SSH library of a chickpea that was infected with Fusarium wilt [68]. Dirigent proteins impart disease resistance through their involvement in lignification during biotic stress. Similar studies have been reported involving Gossypium barbadense that was infected with Verticillium fungus, in which two DIR genes were isolated from the SSH library [69].
Many of the genes that are targeted by miRNAs are involved in disease resistance and growth-related processes. Therefore, it can be surmised that these miRNAs are involved in the regulation of plant development and pathogen growth by acting both as positive and negative regulators, depending on their target genes.

Expression patterns during abiotic stress
Our library allowed for the identification of a large number of conserved salt-responsive miRNAs, including miR390, miR172, miR171, miR169, miR408, miR159, miR396, miR2111, miR5213, miR397, miR393, miR162, miR168, miR166, miR167, miR156, miR530, miR399, miR160, miR319, miR164, miR398, miR2118 and miR394. Among these miRNAs, miR156, miR396 and miR319 were upregulated in response to salt stress, which was confirmed using qRT-PCR. Our results agreed with a previous study involving Arabidopsis, in which 10 salt-responsive miRNAs (miR156, miR165, miR319, miR393, miR396, miR167, miR168, miR171, miR152 and miR394) were reported to be involved in the high salinity stress response [70]. In the chickpea, the transcript levels of miR156 family members were elevated in response to salt stress compared with those of miR166 and others as has been reported in previous studies. Some of the miRNAs that are regulated under salt stress in other plant systems were not found in our library. This phenomenon may be due to different stages or stress conditions; i.e., particular treatment methods or species-specific responses.
Previous studies have demonstrated that miR169 family members are associated with high salt stress [71]. From our target prediction analysis, miR169-targeted genes belong to the nuclear TF family, which contains a CCAAT-binding complex. This CCAAT-binding complex is a eukaryotic promoter element that is evolutionary conserved [72]. Recent studies have demonstrated that these proteins play significant roles in abiotic stress-response pathways [39,73]. The genes that are targeted by miR169 function in transcriptional regulation, suggesting their significant involvement in the salt stress response.
In this study, the salt-responsive miRNA miR390 explicitly targeted protein kinases and the CZF1 TF. The CZF TF is associated with intracellular signal transduction, is involved in the negative regulation of programmed cell death and responds to fungal attack via plant defense mechanisms. CZF1 contains a zinc finger with a CCCH-type domain and has been reported in Arabidopsis thaliana to be salt-inducible. A parallel study in upland cotton reported that the LZF TF acted in response to salt stress [74], and its network of protein-protein interactions was deduced. The chickpea miR396 exhibited higher expression levels under salt stress and was also reported to be salt-responsive in rice. Additionally, transgenic lines over expressing osa-mir396c showed reduced tolerances to salt and alkali stresses compared with wild type plants [75].
In our analysis of the miRNA expression data under both biotic and abiotic stresses, few were upregulated under both types of stresses. miR396 and a member of the miR156 family were upregulated in response to both the wilt and salt stresses at levels of approximately 1.5-fold, indicating the relative similarity between fungal infection-and salinity stress-responses in the chickpea, which was stated in a previous report, in which the chickpea responded to fungal infection (Ascochyta blight) more similarly to high salinity stress than to drought or cold stresses [76]. Additionally, cross talk exists between the stress-signaling pathways that involve several kinases and TFs that are important targeting candidates for several miRNAs under wilt and salt stresses [77,78]. Our data indicate that miR172, miR319, miR171, miR390 and miR396 have serine/threonine protein kinases and MAPK protein kinases as their target genes, which involve signaling pathways. It can be presumed that together, these miRNAs might mediate defense mechanisms under stress conditions via transcriptional regulators. miRNAs also target genes that are directly or indirectly involved in the defense against various stresses. For example, car-miR08 targets a chalcone synthase gene, which is an intermediate in flavonoid biosynthesis. Flavonoids are secondary metabolites that serve variable functions, including those involving pigmentation, UV protection and antifungal defense. Therefore, it can be conferred that these miRNAs come into play during stress management in plants by targeting the genes that are involved either directly or indirectly.
The explicit role of miRNAs in regulating defense mechanisms by the complementary binding of target genes is evident through exhaustive literature reviews. This study will aid in the elucidation of the stress response mechanisms that are utilized by the chickpea. Further, there is limited available knowledge describing comprehensive studies of miRNA expression in the chickpea in response to particular stresses.

Plant materials and stress treatments
The chickpea cultivar ICC4958 was used throughout the study. ICC4958 is a Fusarium wilt-resistant and salt-sensitive chickpea cultivar [79,45]. The plants of the ICC4958 cultivar were grown on a 16-h day/8-h night photoperiod cycle at 2562uC. Fourteenday-old seedlings were subjected to the wilt and salt stresses separately. The stress treatments were performed as follows: for wilt stress, two-week-old plants that were grown under hydroponic conditions were exposed to a toxin that was isolated from the fungus Fusarium oxysporum f.sp. ciceris for one day. For salt stress, the roots of two-week-old seedlings were immersed in a 150 mM NaCl solution for 12 h. All of the tissues (control, wilt-stressed and salt-stressed) were harvested at their respective time points, snapfrozen in liquid nitrogen and maintained at 280uC for further analyses.

Small RNA library preparation and sequencing
Total RNA was isolated using the TRIzol reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer's protocol. For the construction of the small RNA library, low molecular weight (LMW) RNA was enriched by the LiCl method. Equal amounts of RNA were pooled from the root and shoot tissues for each group to generate a LMW RNA library. The RNA was run on a 15% polyacrylamide gel, and the 20 to 30 nt small RNA fraction was extracted and eluted. A preadenylated adaptor was ligated to the 59 end of the small RNAs with T4 ligase. The ligation product was eluted, and subsequently, 39 end adaptor ligation was performed [80] followed by RT-PCR. The PCR products were checked for quality and quantified using a Bioanalyzer (Agilent, Germany). The samples were then se-quenced using the Illumina Genome Analyzer IIx (Illumina Inc., USA).

Computational sequence analysis for identification of miRNAs
The total reads were trimmed and filtered using the UEA small RNA workbench 2.4-Plant version sequence file pre-processing tool (http://srna-tools.cmp.uea.ac.uk/) [81]. The unique tags were generated following a series of processing steps, which included adaptor trimming (using the adaptor removal tool), the elimination of low-quality sequences and the removal of contaminated and other non-coding RNAs, including tRNAs, rRNAs, etc. The UEA sRNA toolkit-Plant version filter pipeline (http://srna-tools.cmp. uea.ac.uk/) was used to exclude the low-complexity and lowquality sequences and eliminate the t/r RNA population by mapping them to plant t/r RNAs from the "Rfam" database, Arabidopsis tRNAs from ''The Genomic tRNA Database'' and plant t/rRNA sequences from the ''EMBL'' release 95. Then, the miRCat pipeline (miRNA categorization) was used to predict novel miRNAs and their precursors using default parameters [82]. The secondary structures of the small RNA sequences were folded using RNAfold (http://rna.tbi.univie.ac.at/cgi-bin/RNAfold.cgi) to predict potential miRNA precursors. The small RNA sequences that had characteristic hairpin structures, together with additional minimal folding free-energy indices (MFEI) [83,84], were considered to be candidate miRNAs by miRCat. The small RNA sequences that matched the following criteria were considered to be valid miRNA precursors: i) no more than 3 consecutive mismatches between the miRNA and miRNA*; ii) at least 17 of the 25 nt surrounding the miRNA must be involved in base pairing; iii) the hairpin must be at least 75 nt in length; and iv) at least 50% of the bases in the hairpin should be paired. The folding structures of the precursors of the new miRNA with the miRNA* were carried out using the UEA sRNA toolkit-RNA hairpin folding and annotation tool, which uses the Vienna Package to obtain the secondary structure of a precursor sequence, highlighting the miRNA/miRNA* sequences on the hairpin structure [85]. The data discussed in this publication has been deposited in Gene Expression Omnibus [86] repository under the accession number GSE57857 (http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc =GSE57857). miRNA validation by poly(A) tail assay-based quantitative real-time PCR (qRT-PCR) The predicted chickpea miRNAs were validated by performing poly(A)-tailed RT-PCR on sixteen miRNAs, including eleven conserved and five novel miRNAs. The total RNAs from the treated and control samples were extracted using the TRIzol reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer's instructions. A 1-mg aliquot of this RNA was used for the poly(A) tailing using the Poly(A) Tailing Kit (Ambion, USA) according to the manufacturer's instructions and then purified using the RNeasyMinElute Cleanup Kit (QIAGENGmBH, Germany). The poly(A) RNA (2 mg) was then reverse-transcribed into cDNA that was primed by a standard poly(T) anchor adaptor using an RTQ primer. For the RT-PCR reaction, the conditions were as follows: 65uC for 10 min, 4uC for 2 min, 50uC for 60 min and 70uC for 15 min. Three biological replicates per sample were used for the analyses.
The poly(T) cDNA was diluted 10-fold and used to perform qRT-PCR using KAPA FAST SYBR Green chemistry (Kapa Biosystems, USA). For the qRT-PCR, the sequences of the specific miRNAs that were validated served as the forward primer and RTQ uni-primer, having an adaptor sequence as the reverse primer (Table S7). The 5S rRNA was used as the reference gene for all of the reactions. Three biological replicates were used per sample in addition to three technical replicates along with a notemplate control and no-RT enzyme control. The data were analyzed using the 2[-Delta DeltaC(T)] method [87] and reported as the means 6 standard errors (SE) of three biological replicates.

Prediction and validation of chickpea miRNA target genes
The chickpea transcript dataset, which was downloaded from the chickpea transcriptome database (CTDB), was used to determine the potential target mRNA candidates for the miRNAs using the psRNATarget program with default parameters (http:// plantgrn.noble.org/psRNATarget/). To reduce the false-positive prediction rate, the cut-off threshold was set at 0 to 3.0 points. Thus, all of the sequences with #3.0 points were considered to be miRNA targets. The functional annotations of the predicted target transcripts were performed using the NCBI nucleic acid and protein databases. Based on the predicted data, miRNA166 was validated using modified 59 RACE. For this validation, the FirstChoice RLM-RACE Kit (Ambion, USA) was used with minor modifications, and the cDNA amplification was carried out using 1 mg of total RNA. A single PCR fragment was cloned into the pGEM-T Easy Vector (Promega, USA) and sequenced to identify the 59end of the target gene.

Analyses of GO terms and KEGG pathways
The GO terms of the target genes were annotated according to their biological processes, molecular functions or involvement as cellular components using Blast2GO [88]. The enzyme mapping of the annotated sequences was performed directly using the GO terms, and the Kyoto Encyclopedia of Genes and Genomes (KEGG) was used to define the KEGG orthologs. Figure S1 Elimination summary of the reads.

(TIF)
Table S1 Conserved miRNAs that were identified in the three libraries with their detailed information. (XLSX)