Comparative miRNAome Analysis Revealed Numerous Conserved and Novel Drought Responsive miRNAs in Cotton (Gossypium spp.)

Negative regulations of gene expression by endogenous, non-coding miRNAs have been shown to play important role in abiotic stress responses in plants. However, limited knowledge is available on water stress responsive miRNAs in cotton. In this study, differentially expressed drought responsive miRNAs and their targets were identified under natural field conditions th rough high-throughput small RNA sequencing by comparing leaf samples of drought tolerant Gossypium hirsutum cv. KC3 and drought susceptible G. barbadense cv. Suvin. Totally four small RNA libraries were constructed and sequenced by employing ion proton technology. Altogether, there were 39 407 089 reads with a mean read length of 22 nt. In KC3, there were 5 138 unique miRNA reads that were differentially expressed with at least two folds under water stressed conditions. In contrast, Suvin have shown 8 469 unique miRNA reads that were differentially expressed with minimum of two folds under water stressed conditions. Comparison of miRNAs expressed under water stressed conditions between KC3 and Suvin, have resulted 7 494 miRNA reads and interestingly majority of them were down regulated with at least two folds. Besides identify ing large numbers of novel miRNAs, several abiotic stress responsive conserved miRNAs were also noticed. Of particular interest were miR750, miR2, miR14, miR276, miR279, miRbantam and miR5176 that were highly down regulated in KC3 under water stress conditions. Strikingly, miR2 and miR-bantam were previously shown to target pre-apoptotic genes in biological systems. Further, the identified miRNAs were also targeting different classes of dehydrogenases, protein kinases and transcription factors. Our results revealed for the first time that there were large numbers of water stress related miRNAs that might be sequentially and/or complexly involved in gene regulation that confers drought tolerance in cotton under field conditions and they have enormous potential in elucidating the molecular mechanism of miRNA based gene regulation and more importantly in genetic improvement of drought tolerance in cotton. Further, this is the first report on experimentally identifying miRNAs in G. barbadense.


Introduction
Diverse types of ribonucleic acids (RNA) with various functions such as coding, noncoding/non-messenger, structural, enzymatic and regulatory functions were reported so far. Among them, RNAs with regulatory activities are termed as micro-RNAs (miRNAs). They are about 22-24 nucleotide (nt) in length and found in all the living organisms including crop plants (Ambros, 2004). As negative regulators of plant gene expression, miRNAs can function by means of minimum two ways: 1) miRNAs base pair with mRNA targets by precise or nearly precise complementarity and cleave / destruct the target mRNA through the RNA interference machinery and 2) they inhibit protein synthesis through an unknown mechanism (Tang et al., 2003). Some known miRNA genic loci form clusters in the genome and these miRNA clusters and the miRNAs in a given cluster are often expressed simultaneously (Wei et al., 2013). Confirmations for the role of miRNA in cell growth and proliferation, developmental timing and phase switch such as vegetative to reproductive growth and response to biotic and abiotic stresses are increasing endlessly and now it is largely deliberated that their regulatory role is more persistent than what was assumed earlier.
Cotton (Gossypium spp.,) is a highly preferred natural textile fiber because of its economic value and comfort and cotton seed oil is a significant food source for human and livestock. Gossypium consists of at least 45 diploid and 5 allotetraploid species; however, only four species viz., G. hirsutum, G. barbadense, G. arboreum and G. herbaceum (described in the order of highest to lowest global acreage), have been domesticated for their abundant seed trichomes and provide the foundation for the global textile industry. India has the largest area dedicated to cotton production worldwide; however, the cotton fiber production is largely limited by water stress since more than 60% of this area is under rainfed (Boopathi and Pathmanaban, 2012). Genetic improvement of cotton with drought resistance will be a viable and cost effective option to mitigate such losses. Though significant advancement has been documented, the progress in this direction is not up to the expected level due to lack of knowledge on drought resistance mechanisms. Although progress has been made in unravelling the complex stress response mechanisms, particularly in the identification of water stress responsive protein-coding genes, the key factors such as miRNAs that regulate the expression of these genes remain unclear. Such gap severely hampers the precise manipulation of drought resistance in cotton. Even though there is no direct proof that miRNAs are the crucial regulatory elements in development of drought tolerance in cotton, we postulate that differential expression of miRNAs might be involved in such complex process. As a first step towards the understanding of their regulatory mechanisms and networks of target genes in drought tolerance in cotton, differential expression of miRNAs between drought tolerant (Gossypium hirsutum cv. KC3) and drought susceptible (G. barbadense cv. Suvin) cotton lines were compared in the present study under water stressed and irrigated conditions in the field using a deep sequencing approach.
In general, miRNAs are being identified in plants by forward genetics (whereby miRNAs are isolated from an organism showing abnormal phenotypic characteristics), reverse genetics (whereby a specific gene is knocked out or over expressed), cDNA sequencing and computational prediction approaches (Ambros, 2004). In the early days, the identification of cotton miRNAs depended largely on computational approaches based on the conserved characteristics of miRNAs (Qiu et al., 2007;Zhang et al., 2007;Barozai et al., 2007), which was accomplished by searching the cotton genomic sequences and/or ESTs that were homologous to known mature miRNAs or pre -miRNAs. Later, after complete sequencing of G. raimondii  and G. arboreum (Li et al., 2014) it was proposed that there were totally 348 and 431 miRNA sequences, respectively in these two species.
The first experimental validation and exploration of cotton miRNAs by cloning and sequencing identified only three cotton miRNAs (Abdurakhmonov et al., 2008). However, recently deep sequencing has emerged as an effective approach for large-scale miRNA discovery. Small RNA libraries prepared from cotton ovules and seedlings have been analysed by high-throughput sequencing (Kwak et al., 2009;Ruan et al., 2009). Pang et al. (2009) reported 27 conserved and four novel miRNA families from G. hirsutum by sequencing the small RNAs of leaves and ovules. Along with data from microarrays and Northern blots, they examined the changing pattern of certain miRNAs and their targets. Similar comparative miRNAome analysis revealed seven fiber initiation-related miRNAs expressed in cotton ovules besides 36 novel miRNAs and two conserved miRNAs . Targets of these miRNAs were experimentally validated and they described complex regulatory networks that coordinate fiber initiation responses. Yin et al. (2012) reported genome-wide profiling of miRNAs and other small non-coding RNAs in the Verticillium dahliae inoculated cotton roots. Small RNAs and their targets were also identified during cotton somatic embryogenesis through high-throughput small RNA and degradome sequencing .  have identified 93 miRNAs, including 63 novel miRNAs that are involved in fiber initiation and seed development. Similarly, 257 novel miRNAs were reported by Xue et al. (2013) in elongating cotton fiber cells. In another comparative miRNA expression study, Wei et al. (2013) identified sixteen conserved miRNA families during anther development in genetic male sterile and wild type cotton.  has studied genetic variation of miRNAs and their target genes and they were genetically mapped in cotton.
Despite of these progresses, only 378 mature miRNA sequences from cotton (Gossypium spp.), were described in miRBase (http://www.mirbase.org/cgi-bin/browse.pl; Release 21; June, 2014). As several thousands of miRNAs were identified in plants, this number is extremely lesser than the expected number of miRNAs in cotton and they are remain to be recognized. Further, among the reported miRNAs, only very few that were specific to drought and other abiotic stress has been reported in cotton using computational (Boopathi and Pathmanaban, 2012) or experimental (Yin et al., 2011;Wang et al., 2013) approaches, even though such abiotic stress specific miRNAs have been expressively established in Arabidopsis, rice and other crops (Sunkar et al., 2008). Identification of miRNAs in cotton, in response to water stress, certainly have important implications for gene regulation under drought and such knowledge can impressively enhance the drought tolerance breeding efficiency besides contributing significantly to the goal of having a complete profile of miRNAs in cotton. Besides, knowledge on miRNA-guided stress regulatory networks would provide new tools for the genetic improvement of other abiotic stress tolerances. In addition, all the miRNA studies employed G. hirsutum, G. arboreum, G. herbaceum and G. raimondii. However, miRNA expressions in G. barbadense, which has superior fiber quality traits that suit to modern textile mills, remain to be explored. Hence, in view of the importance of miRNAs in regulation of water stress resistance in cotton, this study was conducted with the following objective: to capture and characterize the differentially expressed miRNAs in drought tolerant and susceptible cotton cultivars using next generation sequencing platform.

Small RNA libraries
Totally four different small RNA libraries were constructed from the leaf tissues of KC3 and Suvin collected from two extreme water regimes and intensely sequenced using ion proton TM technology. After the deletion of low quality reads and several kinds of contaminant tags 9 826 785 clean reads from leaf tissues of KC3-C were obtained (Table 1). Similarly, clean reads of 8 672 895 from KC3-WS; 8 627 325 from Suvin-C and 8 017 500 from Suvin-WS were obtained (Table 1) and they were used for analysis of differentially expressed small RNAs. Altogether, there were 39 407 089 reads found in all the four small RNA libraries with a mean read length of 22 nt. The sequence reads of the each library was normalized to RPKM and each reads that were differentially expressed were recorded in RPKM (Table 2 and Figure 1). Small RNA reads with a size of 27 nt were found to be present in significantly greater number in all the four comparisons followed by 28 or 29 nt reads ( Figure 1). Small RNAs in the size range of 24-33 nt that were differentially expressed in drought susceptible cotton cultivar, Suvin, under two different water regimes were significantly higher than the other three comparisons. However, predominantly the size range of 20-23 nt were differentially expressed between KC3-WS and Suvin-WS ( Figure 1). Totally 8 469 differentially expressed distinct reads were found between Suvin-C and Suvin-WS followed by 7 494 differentially expressed distinct reads between KC3-WS and Suvin-WS (Table 2). ESTs from NCBI and the genome sequence of G. raimondii , were used as reference sequences to predict the miRNA precursors since no allopolyploid cotton genome sequence was available at the time of this analysis. Nevertheless, to a maximum of 15.35% of the small RNA reads obtained from this study were mapped with G. raimondii genome sequences ( Table 2). The distinct reads that were identified as rRNA, snRNA, snoRNA, tRNA, tmRNA during the analysis (number of reads in each small RNA class is given in Table 2) were separated and the remaining reads (which constitutes the majority of the reads; Table 2) were used further.
1.2Differentially expressed conserved miRNA families RPKM reads of each miRNA were used to notify their abundance and differential expression in two different contrasting cotton cultivars under two different water regimes. Simple blast analysis with matured miRNA   Table 1). Interestingly, several miRNA families were found to exhibit more than 10 fold change in their expression level in KC3-WS (for example, miR397 and miR398 were more than 10 fold up regulated and miR160, miR164, miR166, miR390, miR399, miR479, miR2111 and miR2796 were more than 10 fold down regulated; Supplementary Table 1). In contrast, some of the miRNA family had lowest fold change (miR535 was in 2.031 fold up regulation and miR160 was in 2.060 fold down regulation).
When comparing the miRNAs that were differentially expressed between the miRNA libraries prepared from drought tolerant, KC3 and drought susceptible, Suvin under water stress in the natural field conditions, totally 142 conserved miRNAs that belongs to 29 conserved families were noticed (Supplementary  Table 1). There were 14 miRNAs belonging to 11 miRNA families were up regulated, 44 miRNAs belonging to 18 miRNA families were down regulated and seven miRNA families were both up and down regulated (Supplementary Table 1). Remarkably, pronounced fold change (between 2.022 and 3473.119) under down regulated miRNAs was identified in this situation. In the same way, extensive variation in the fold change (between 2.024 and 38.453) under up regulated miRNAs was also noticed. The highest induction was observed for miR399 (38.453 folds) followed by miR166 (9.173 folds) and miR160 (7.819 folds) in Suvin under irrigated conditions and highest down regulation (perhaps, the highest down regulated miRNA that was found in this report) was seen for miR2118 (which recorded 3473.119 folds) followed by miR167 ( (miR14)]. It was further noticed that miRNA families, hme-bantam, miR2, miR14, miR156, miR160, miR166, miR279, miR2111, miR3027, and miR3050 were immensely downregulated and one member of miR169 family was significantly up-regulated in Suvin-WS (Supplementary Table 1).

Identification of novel miRNAs
In addition to conserved miRNAs, deep sequencing had also helped to identify new miRNAs in cotton that have not yet reported in miRBase. Since, the sequencing of the G. hirsutum and G. barbadense genome is incomplete and information on complete genome-wide cotton small RNA population is unknown, accurate identification of non-conserved miRNA in cotton is a difficult task. In order to ensure the truthful identification of miRNA and considering of the complicated nature of plant small RNA, a series of strict filtrations was used to enhance reliability of the result presented here. ESTs from NCBI and the genome sequence of G. raimondii  were used as reference sequences to predict several numbers of miRNA precursors. However, with the availability of genomic sequencing data from G. hirsutum and G. barbadense in near future, there is a greater potential to find several additional miRNAs of Gossypium with this same data sets.
Though several hundreds of new miRNA loci were identified as eligible candidates, following a BLASTn search and hairpin structure prediction, only 59 putative unique miRNAs were detected in all the four small RNA libraries (Table 3). The RNA structures were predicted by mFold software and manually checked according to the criteria of Meyers et al. (2008). The lowest minimum free energy (MFE) among all hairpin structures of the novel miRNAs precursors was −43.9 kcal/mol (Table 3). All precursors of novel miRNAs had regular hairpin structures and two of these (Ghi-miR(contig-4796) and Ghi-miR(contig-13785) are presented in Figure 3. GC content of the identified new miRNA precursor sequences was ranged from 23-71 percent (Table 3).  (Table 4). Contig_16568 accumulated abundantly with more than 17 RPKM. Contig_16528, Contig_1094, Contig_15570 and Contig_5490 were also detected in relatively high levels with more than 10 RPKM in these two libraries.
Interestingly, the majority of the novel miRNAs was detected in more than 3.5 RPKM (Table 4).
Cotton Genomics and Genetics 2016, Vol.7, No.2, 1-23 http://cgg.biopublisher.ca 8 Figure 3 Secondary hairpin structures predicted for putative novel miRNAs. The highlighted nucleotides are the sequence of putative novel miRNA found in this study and the remaining are the sequences of G. raimondii.

Targets of miRNAs identified in this study
It is known that miRNAs mediate plant abiotic stress responses through regulating their target genes and thus serving as key players in the gene regulation networks. The targets of miRNAs were predicted by the web tool psRNATarget (Dai and Zhao, 2011) using the Gossypium (cotton) DFCI Gene Index (CGI) Release 11 (http://compbio.dfci.harvard.edu/index.html) as the sequence library for the targets search. The predicted targets of some of the drought-responsive miRNAs, identified in this study, were found to be involved in diverse cellular processes including development, transcription, protein degradation, detoxification, nutrient status, and crossadaptation to various abiotic stresses (data not shown). Besides, miRNAs reported here have also shown (elsewhere) to target several genes that are associated with abiotic stress resistance.
sequence alone is sufficient for the identification of the physiologically relevant target genes (e.g., hme-bantam by Brennecke et al., 2003).
This study has helped to identify several drought responsive candidate miRNAs that has potential in revealing the role of small RNA in drought tolerance in crop plants and cotton in particular. Among the several differentially expressed conserved miRNAs, those that were found between KC3-WS and Suvin-WS has more pronounced role in deciphering the molecular mechanism of drought resistance in cotton. The drought responsive miRNAs identified in this study were also found to be associated with drought or other abiotic stresses elsewhere and they were described hereunder.
The conserved miRNA, hme-bantam with 23 nt length has found to be highly down regulated (35.018 fold change in RPKM reads; Supplementary Table 1) in the drought susceptible parent, Suvin. Bantam microRNA in Heliconius melpomene was predicted independently by two groups using computational approaches (Brennecke et al., 2003;Lai et al., 2003). Northern blotting confirmed that hme-bantam with a 23 nt sequence was the most commonly expressed. Interestingly, bantam microRNA simultaneously stimulates cell proliferation and prevents apoptosis and targets the 3' UTR of the pro-apoptotic gene hid in Drosophila. In addition to bantam miRNA, miR2 family have also been shown to regulate the expression of the proapoptotic genes such as reaper, grim and skl and to limit apoptosis in Drosophila (Zhang and Cohan, 2013). Two members of miR2 was identified in this study and both were highly down regulated (159.578 and 34.575 folds; Supplementary T able 1) in the drought susceptible, Suvin. Drosophila melanogaster presents perhaps the most important model organism for understanding the basic principles and molecular mechanisms of animal development. Similarly, Drosophila genetics has played an important role in understanding the functional roles of several animal and plant miRNAs. Interestingly, high level expression miR2 in root tissues were reported in Mangrove (Khraiwesh et al., 2013) and more specifically differential expression of miR2 under drought was noticed in maize (Xu et al., 2010). Another report in cucumber has also shown that miR2 targeted a gene coding a tetratricopeptide repeat-containing protein (Martinez et al., 2011) that usually interacts with heat shock proteins (Yin et al., 2011). Thus, it can be strongly hypothesised that vigorous up regulation of hme-bantam and miR2 in the drought tolerant KC3 accelerates antiapoptotic activity in the cells during water stress and provide useful leads to promote research that confirms that these miRNAs endorses drought tolerance in cotton. Similarly several other conserved miRNAs identified in this study were found to be linked to water or other abiotic stresses and cellular processes.
For example, miR14 was 76.243 fold down regulated in Suvin under water stress in this study (Supplementary Table 1). The targets for mir14 include cell-death effectors other than hid (Xu et al., 2003) and mutants of miR14 were shown to be sensitive to high salinity (Leung and Sharp, 2010). In foxtail millet, miR14 has reported to target GRAS transcription factor that has been involved in development and other processes (Yi et al., 2013). This study has identified that there were nine members of small RNAs (maximum of 4.058 fold down regulated) under water stressed conditions in Suvin and all of them have shown > 90% percent identify with miR156 family. Squamosa Promoter Binding like (SPL) proteins are negatively regulated by the miR156, whose cellular levels are higher in younger plants and progressively decrease as the plant ages and especially during flowering Yamaguchi et al., 2009). Among the miRNA-target that were studied in cotton, miR156-SPL2, shown significant regulation relationship in roots and leaves under salinity stress . Up regulation of miR156 in rice shoot and root tissues in response to heat (Sailaja et al., 2014) and drought (Zhou et al., 2010) stresses was also reported as it found in KC3-WS in this study. Further, high-throughput sequencing and expression analysis of microRNAs and their targets under salt stress in Caragana intermedia has shown that the miR156 family exhibited the highest abundance among all of the miRNA families .
More than 8 fold up regulation of miR159 was noticed in KC3-WS when compared with Suvin-WS (Supplementary Table 1). Ding et al. (2013) reviewed that miR159 was induced by ABA and drought treatments and miR159 play a key role in ABA response by directing the degradation of MYB mRNAs. Bertolinia et al. (2013) addressed the role of miR159 in reprogramming leaf growth during drought stress in Brachypodium Cotton Genomics and Genetics 2016, Vol.7, No.2, 1-23 http://cgg.biopublisher.ca distachyon. Shuai et al. (2013) shown that the expression patterns and RT-qPCR results of drought-responsive miRNAs in Populus were consistent and both indicating that miR159 members were up-regulated after drought treatment. Similar kind of up regulation of miR159 under drought was also reported in legumes (Mantri et al., 2013), finger millet (Nageshbabu et al., 2013) and in cotton .
Eight members of up regulated miR162 with maximum elevated expression of 15.071 folds (Supplementary Table  1) was found in water stressed KC3 leaves in this study. Similar kind of up regulation of miR162 under drought was reported in cotton  and Brachypodium (Bertolinia et al., 2013). Up regulation of miR162 under salt stress was has also been reported in roots of Vigna unguiculata (Paul et al., 2011). In cotton, Ruan et al. (2009) have shown that miR162 target zinc finger protein and miR162 was also found to be key regulator of dicer-like protein (DCL) that are essential for miRNA or small RNA biogenesis (Liu et al., 2009 , 2013) which is shown to be involved in ABA and abiotic stress signalling in several plant species (Nicolas et al., 2014). Clear role for up and down regulation of miR166 was not yet deciphered; hence, additional studies on miR166 under drought stress should shed more information on precise molecular details on stress resistance mechanisms in plants.
Down regulation of 5 members of miR167 was noticed in Suvin-WS (Supplementary Table 1). Microarray analysis showed that miR167 was up regulated by drought stress in Arabidopsis and miR167 has been reported to play major roles in drought and ABA response in plants (Liu et al., 2008). Similar to more than 5 fold down regulation of miR169 in KC3-WS, down regulation of miR169 has also been noticed in drought resistant, soybean (Ni et al., 2013). It has experimentally validated that GmNFYA3 (a gene encoding Nuclear factor Y (NF-Y), which is a heterotrimeric transcription factor), was target gene of miR169 and demonstrated that GmNFYA3 is a positive regulator of plant tolerance to drought stress (Ni et al., 2013). In addition, Xu et al. (2014) evidently demonstrated that miR169 family members might be involved in flowering time regulation under abiotic stresses in plants including drought. Plants often flower earlier under water stress conditions and early flowering is considered as one of the key drought resistance component traits in plants (Komashita et al., 2008). The stress-activated transition to flowering enhances the chance of a plant population surviving under the threatening water stress.
All the four members of miR171 were up regulated in KC3-WS (Supplementary Table 1). Induced expression of miR171 was shown in barley under dehydration stress and the identified targets for miR171 was scarecrow like 6 (SCRL6). It encodes a transcription factor that is involved in diverse plant developmental processes such as leaf or root (Kantar et al., 2010).
Elevated up regulation of miR276 (57.625 fold; Supplementary Table 1) was found in this study in KC3-WS and to the best of our knowledge, there is no report in plant that specifies the role of miR276 under water stress. However, it has been reported in insects and shown to have key role in larval development and age dependent behavioural changes (Hori et al., 2010;Kozomara and Griffiths-Jones, 2011).
Yet another miRNA family that has high level of induction in KC3-WS were two members of miR279 (maximum of 58.068 fold changes; Supplementary Table 1). Up regulated expression of miR279 was mainly found in roots and it is expressed at lower levels in leaves, and scarcely expressed in stems and booting panicles in the common wild rice, O. rufipogon . However, their role in drought tolerance in plant is yet to be determined at molecular level.
In our previous study in cotton, we have computationally identified miR319 and established through in silico analysis that miR319 might be involved in abiotic stress resistance (Boopathi and Pathmanaban, 2012). Interestingly, miR319 was found to be induced 4.939 folds in KC3-WS (Supplementary Table 1). Transgenic plants overexpressing miR319 displayed morphological changes and exhibited enhanced drought and salt tolerance associated with increased leaf wax content and water retention but reduced sodium uptake (Zhou et al., 2013). As one of the first experimentally characterized and most conserved miRNA families, the miR319 targets Teosinte Branched / Cycloidea / Proliferating cell factors (TCP) genes, which encode plant -specific transcription factors sharing a conserved TCP domain with a basic helix-loop-helix structure. The TCP family is known to be largely involved in plant development such as the control of cell proliferation in leaf morphogenesis Zhou et al., 2013).
It is known that miR393 was commonly up regulated during drought stress in Arabidopsis, rice and sugarcane (Ding et al., 2013). The target of miR393 encodes TIR1 (transport inhibitor response 1), an auxin receptor in Arabidopsis. The TIR1 enzyme is a positive regulator of auxin signalling by promoting the degradation of Aux/IAA proteins through ubiquitination. Thus, increased levels of miR393 would down regulate auxin signalling and may reduce plant growth under drought stress (Ding et al., 2013). It was also shown in this study that miR393 was down regulated more than two folds in KC3-WS; Supplementary Table 1) which might favour the growth of KC3 under drought. Similarly, more than two fold down regulation of miR394 was noticed in KC3-WS. Downregulated miR394 was noticed in Populus, whose predicted targets are annotated as dehydration-responsive protein and F-box proteins, which were reported to play significant roles in the abiotic stress-response pathway (Shuai et al., 2013).
High induction of 36 members of miR396 (maximum of 11.456 folds) was displayed in KC3-WS in this study (Supplementary Table 1). miR396 family was found to have maximum number of members in this study and interestingly, all were up regulated in KC3-WS. Up regulation of miR396 under water stress was also noticed in Arabidopsis and tobacco (Liu et al., 2009;Yang and Yu, 2009). miR396 has been shown to target six growth regulating factor (GRF) transcription factors with roles in the coordination of cell division and differentiation during leaf development in Arabidopsis (Wang et al., 2011). Transgenic Arabidopsis plants overexpressing miR396 displayed narrow-leaf phenotypes because of reduction in cell number, achieved through repression of the expression of GRF genes and were more tolerant to drought than wild-type plants, likely because of the lower stomatal density (Liu et al., 2009). Overexpression of miR396 also conferred increased tolerance to drought stress in tobacco (Yang and Yu, 2009). These observations suggest that miR396 plays important roles not only in leaf development but also in drought tolerance in plants (Ding et al., 2013).
Three members of miR399 were up regulated (maximum of 7.150 fold) in KC3-WS (Supplementary Table 1). miR399 was found as one of the largest family in plants and it appears to target four sites in the 5′ UTR of a putative ubiquitin conjugating enzyme E2 (UBC; Sunkar and Zhu, 2004). It was found that there is no substantial regulation of miR399 by salt, drought, or cold stress. Low K or low N did not induce the expression of miR399. But miR399 was highly induced (whereas the target UBC mRNA was reduced) by low-phosphate (Pi) stress (Lu and Huang, 2008;Mantri et al., 2013).
Interestingly, 15 members of miR482 were up regulated (to the maximum tune of 15.957 folds) in KC3 under water limited conditions (however, one member of miR482 was down regulated (2.085 fold); Supplementary  Table 1). Up regulation of miR482 under drought was reported in finger millet (Nageshbabu et al., 2013). Four UDP-glucosyltransferases (UDPGs) were found to be targeted by miR482 (Shuai et al., 2013). The UDPGs are enzymes that attach a sugar molecule to a specific acceptor in plants. As in Arabidopsis, the UDPG is a key regulator of stress adaption through auxin IBA (Tognetti et al., 2010) and plays a role in fine-tuning nitrogen assimilation in cassava (Kannangara et al., 2011). miR482 has also shown to target resistance gene receptor kinases (Nageshbabu et al., 2013;Shuai et al., 2013) and involved in disease resistance in cotton  and soybean nodulation (Nageshbabu et al., 2013). The ability of drought stress to increase and decrease the Cotton Genomics and Genetics 2016, Vol.7, No.2, 1-23 http://cgg.biopublisher.ca 13 abundance of same members of miR482 highlights the need for more in-depth and detailed characterization of stress-responsive miRNAs in plants, particularly in cotton.
Up regulation of all the 10 members of miR535 family were reported in this study in KC3-WS (Supplementary  Table 1). Similarly, miR535 was the most abundant miRNA in drought-exposed peach (Eldem et al., 2012). It has been suggested that stress-regulated miRNAs might be regulated by stress-related transcription factors. Interestingly, the promoters of miR535 in Arabidopsis contained all of the stress-related cis-elements investigated, including the drought-responsive element, low temperature-induced element, MYC-binding site, and MYBbinding site (Abdel-Ghany and Pilon, 2008).
The most abundant conserved miRNA family that has highly induced in drought tolerant cotton cultivar KC3-WS is miR750: there was 196.813 orders of magnitude of changes when compared with drought susceptible, Suvin. It has been shown that expression of miR750 was the most abundant miRNA after infecting silkworm with Bombyx mori cytoplasmic polyhedrosis virus and has six target genes including kinases . Nevertheless, the role of miR750 in plant abiotic stress resistance is yet to be determined. Similarly, high accumulation of several conserved miRNA families such as miR3027, miR3050, miR5176, miR7496, miR7499 and miR7500 (20.390, 15.514, 33.688, 6.660, 11.525 and 7.535 folds, respectively) were found in KC3-WS. However, their role in drought or abiotic stress resistance in crop plants have not yet been reported. Thus, this study has identified several drought responsive highly induced conserved miRNAs in cotton and they provided strong lead to specifically characterize such miRNA roles in drought tolerance in plants, particularly in cotton.
Two members of miR2118 were up regulated in KC3-WS with maximum of 7.806 folds (Supplementary Table 1). Similar kind of accumulation of miR2118 upon drought and ABA treatments in legumes were reported and showed that they may target regulatory elements for cellular processes (Arenas-Huertero et al., 2009). In our previous study in cotton, we have computationally identified miR2118 and established through in silico analysis that miR2118 might be involved in abiotic stress resistance (Boopathi and Pathmanaban, 2012). In addition, significant level of expression of miR2118 was reported in drought tolerant cowpea genotype (Barrera-Figueroa et al., 2011) and it was also up regulated in response to drought stress in Medicago truncatula (Wang et al., 2011). Further, it has been shown that miR2118 targets TIR-NBS-LRR domain protein, which usually expressed in response to multiple abiotic stresses such as drought, cold, salinity and ABA (Wang et al., 2011;Mantri et al., 2013).
More than six fold up regulation of two members of miR2949 and two fold induction of miR3476 were displayed in KC3-WS (Supplementary Table 1). Though the role of miR2949 and miR3476 in drought tolerance has not been reported in plants, high accumulation of miR2949 and miR3476 in cotton was predicted as positive regulators of the process of fiber initiation  and somatic embryogenesis .
The increased rates of miRNA detection afforded by deep-sequencing technologies provide challenges to the level of confidence required to annotate a sequence as a miRNA. A typical RNA deep-sequencing experiment will identify millions of short sequences. Increased coverage results in detection of sequences of ever-lower abundance. It therefore becomes more and more challenging to distinguish true miRNAs from fragments of other transcripts, other short RNAs and spurious transcription. Guidelines for microRNA annotation were established in 2003(Ambros et al., 2003, requiring evidence of expression of a ~22 nt sequence (for example, cloning, sequencing or northern blot), together with evidence for a microRNA precursor structure (predicted stem-loop flanking the mature sequence). Updated annotation criteria were recently suggested to distinguish microRNAs from other classes of short RNAs in plants (Meyers et al., 2008). These standards have proved extremely powerful in maintaining a clean data set of miRNA sequences for the community. The eukaryotic genome also contains millions of predicted hairpins, so a flanking stem-loop structure should be considered necessary but not sufficient to annotate a sequence as a miRNA. If poorly analysed, a single data set thus has the potential to generate a large number of suspicious annotations, thus overwhelming the real miRNAs. However, correct interpretation of RNA deep-sequencing data provides several additional signals to help distinguish miRNAs from other sequences (Kozomara and Griffiths-Jones, 2011).
The lowest minimum free energy (MFE) among all hairpin structures of the novel miRNAs precursors was −43.9 kcal/mol (Table 4), which is significantly lower than the threshold of −30 kcal/mol reported in a previous study [25]. Among the 36 novel miRNAs identified in this report, none of the novel miRNA had homologs in miRBase Release 19. Hence, we assumed that they were Gossypium-specific, or restricted to closely related species. Further, there is a chance to get more novel or conserved miRNAs if new genomic resources are developed in near future as discussed above.
In conclusion, our previous computation study (Boopathi and Pathmanaban, 2012) has shown that conserved miRNAs were shown to be linked with abiotic stress resistance. Here we provided experimental proof for those and other miRNAs for their roles in drought tolerance in cotton. During water stress responsive transcript profiling, it is important to target the different tissues and their precise stage besides the dynamics of the stress (such as time, duration and intensity), for isolation of RNA. Thus, it will be more useful if the ultimate economically important tissues, the flower, boll and fiber, are used during sampling instead of leaves that was used in this study. Further, instead of using genotypes with different genetic backgrounds, near-isogenic lines, which differ only in the target trait, would be ideal genetic material which ensures that differentially expressed genes are linked to the drought resistance and not to the genetic background. Our laboratory has already started to work on these lines. The number of new miRNAs identified in this study is relatively low when compared with other crops. This is mainly due to lack of genome sequence in G. hirsutum and G. barbadense. Release of their genome sequences in the near future may help to identify more number of novel miRNAs in cotton in response to drought stress using this data.

Plant materials
Two contrasting accessions of Gossypium, which are considered as ideal genetic materials for studying drought tolerance mechanism in cotton, were used in this study. One was short duration upland cotton, G. hirsutum L. var. KC 3 (hereafter referred as KC3) that has good performance under rainfed tract of Tamil Nadu, India and provide assured minimum yield with relatively good quality fiber even under severe water stress. Another line was the long duration G. barbadense L. var. Suvin (hereafter referred as Suvin) which provide poor yield with inferior fiber qualities under drought stress . These two lines have also shown significant dissimilarity at molecular level (Thiagu et al., 2011) and hence they were selected for development of different types of mapping population at this laboratory. Therefore, the differentially expressed miRNAs identified in this study could also be genetically mapped, since there is no genome sequence currently available in G. hirsutum or G. Barbadense.

Experimental conditions in the field
Field trial was conducted under upland condition at TNAU, Coimbatore, India during the regular cotton cultivating season, kharif (July -December, 2012). Both the cotton lines were grown in a single row by following 75 x 40 cm spacing between and within rows, respectively with two replications and evaluated under two water regimes: irrigated (control or non-stress) and water limited (stress). A buffer area of 3.0 m wide and 0.75 m deep along the length of the experimental plot divided the control and stress plots. All the standard agronomic practices were equally applied to both the treatments. Both the lines were surface irrigated as per the regular cultivation practices, except for water limited regime. Water stress was imposed by withholding irrigation when the plants were in peak vegetative phase to stress plots i.e., from 30 days after sowing (DAS) and there were continuous rain free period of 20 days. The leaf samples were collected from both the accessions and treatments on 50 th DAS with two replications and stored in Shrimpex lysis buffer (Shrimpex Biotech Services Private Limited, Chennai, India).

Small RNA library construction and sequencing
Cotton Genomics and Genetics 2016, Vol.7, No.2, 1-23 http://cgg.biopublisher.ca Small RNA library construction and sequencing were done at Shrimpex Biotech Services Private Limited, Chennai, India. Briefly, the miRNA was isolated according to the manufacturer's protocol and concentration was measured using Nanodrop Spectrophotometer. To avoid DNA contamination, the miRNA was treated with DNase @ 0.5 unit per μg of miRNA (the reaction mixture contains, 0.5 unit of DNase, 10X DNase buffer and the miRNA template to a final volume of 30 μL) and incubated at 37°C for 30 minutes. The reaction was stopped by incubating at 65°C for 10 minutes with 1 μL of 50 mM EDTA. The miRNA was purified from DNase using Shrimpex spin column concentration module according to the manufacturer's protocol. The quality of the obtained miRNA before and after DNase treatment was analyzed using Agilent 2100 small RNA kit.
The miRNA library was prepared using Ion Total-RNA Seq kit v2 and the miRNAs were enriched using magnetic bead clean up module for small RNA enrichment. In the next step, Ion RNA-Seq adapter was ligated to the enriched miRNA to facilitate reverse transcription of first strand cDNA synthesis and further purification using magnetic bead clean up module. The purified cDNA were amplified using a barcoded forward primer and a common reverse primer and quantified using Agilent Bioanalyzer 2100 DNA 1000 kit. Smear analysis was performed to determine the size distribution range of the four samples to confirm that the size of the four libraries was in the range of 94~114 bp.
The molar concentration of each library was determined and diluted to the same molar concentration. Equal volumes of each barcoded library of the replicated samples were mixed to prepare a pool of barcoded libraries and were amplified using emulsion PCR (ePCR). In ePCR, the clonal bead populations are generated in water-in-oil microreactors. After ePCR, the templates were denatured and bead enrichment step was performed to separate the beads with templates from non-template beads. The enriched beads were then loaded onto the Ion P1 Chip and sequenced.

Sequence analysis for miRNA identification
Sequences were trimmed to eliminate low-quality reads and other contaminants from the adaptor tags. All low quality reads were removed such as reads with 3' and 5' adapter contaminants, those without insert tags, and those with poly-A sequences. The remaining high-quality sequences were trimmed of their adapter sequences. All highquality sequences were considered significant and number of reads in each nucleotide class was classified. Each read were annotated through BLAST searching for identification of rRNA, scRNA, snoRNA, snRNA and tRNA using 'Rfam' (ftp://ftp.sanger.ac.uk/pub/databases/Rfam/) and those matched reads were avoided for further analysis.
Initially, a simple blast search of resulted small RNA sequences identified in both the replications (in the size range of 20 -33 nt) with miRBase was done to categorize their family classification based on conservation (version 19 (August, 2012), http://microrna.sanger.ac.uk). Sequences that were identical or related to known miRNA sequences (with maximum of single nucleotide substitution) were considered as known or conserved miRNAs and others were considered as putative novel miRNAs.

Identification of novel miRNAs
Considering the complicated structure of miRNA gene and its biogenesis, a series of strict filters were used during the new miRNA identification process to enhance reliability of the result presented here. Currently, no allopolyploid cotton genome sequence was available; therefore ESTs from NCBI and the genome sequence of G. raimondii , the closest living relative of the progenitor D-genome donor of allotetropolyploid cottons (G. hirsutum and G. barbadense) were used as reference sequences to predict the miRNA precursors. The 20-33 nt reads that were assumed as putative novel miRNA were mapped to cotton nucleotide sequences available at NCBI and D genome sequence (http://cgp.genomics.org.cn/page/species/index.jsp). Approximately 100 nucleotides flanking both side of the query sequence were obtained and the characteristic secondary hairpin structure of each miRNA precursors was predicted using the miRNA prediction software, mFold 3.2 (Zuker, 2003). During this process the following characteristics were considered to affirm novel cotton miRNAs: i) formation of flawless secondary structure (including formation of a perfect stem-loop structure, the read sequence of this study sat at one arm of the stem and asymmetric bulges was minimal in size (one base at most) and infrequent (typically one or less) within the miRNA/miRNA* duplex) and ii) the minimum free energy should be lower than the threshold of−30 kcal/mol (Bonnet et al., 2004).

Differentially expressed miRNAs
Differential expression analysis of both conserved and novel miRNAs between the two contrasting genotypes under two different water regimes were identified as described below. The miRNA reads of the each libraries were assembled de nova into contigs and normalized to reads per kilobase per million (RPKM; Mortazavi et al., 2008) and RPKM was used as a statistical value for the analysis of differences in expression among the samples. The increased or decreased level of expression (up-or down-regulation) of miRNA in each sample under given treatment were represented in RPKM fold change (as described in Rada-Iglesias et al., 2011) and minimum of two fold change in RPKM were considered further. Such analysis was helped to identify major (thus insignificant miRNAs were not included in this report) and differentially expressed miRNAs in leaf samples between i) KC3 irrigated control (hereafter mentioned as KC3-C) and KC3 water stressed (hereafter mentioned as KC3-WS) ii) KC3-WS and Suvin water stressed (hereafter mentioned as Suvin-WS) iii) Suvin irrigated control (hereafter mentioned as Suvin-C) and Suvin-WS and iv) KC3-C and Suvin-C.

Identification of miRNA's targets
The targets of miRNAs was predicted by the web tool psRNATarget (Dai and Zhao, 2011) using the Gossypium (cotton) DFCI Gene Index (CGI) Release 11 (http://compbio.dfci.harvard.edu/index.html) as the sequence library for the targets search. No more than three mismatches between miRNAs and targets were allowed, especially no mismatches were allowed within the maximum expectation region.