Identification of novel drought-responsive microRNAs and trans-acting siRNAs from Sorghum bicolor (L.) Moench by high-throughput sequencing analysis

Small non-coding RNAs (sRNAs) namely microRNAs (miRNAs) and trans-acting small interfering RNAs (tasi-RNAs) play a crucial role in post-transcriptional regulation of gene expression and thus the control plant development and stress responses. In order to identify drought-responsive miRNAs and tasi-RNAs in sorghum, we constructed small RNA libraries from a drought tolerant (M35-1) and susceptible (C43) sorghum genotypes grown under control and drought stress conditions, and sequenced by Illumina Genome Analyzer IIx. Ninety seven conserved and 526 novel miRNAs representing 472 unique miRNA families were identified from sorghum. Ninety-six unique miRNAs were found to be regulated by drought stress, of which 32 were up- and 49 were down-regulated (fold change ≥ 2 or ≤ −2) at least in one genotype, while the remaining 15 miRNAs showed contrasting drought-regulated expression pattern between genotypes. A maximum of 17 and 18 miRNAs was differentially regulated under drought stress condition in the sensitive and tolerant genotypes, respectively. These results suggest that genotype dependent stress responsive regulation of miRNAs may contribute, at least in part, to the differential drought tolerance of sorghum genotypes. We also identified two miR390-directed TAS3 gene homologs and the auxin response factors as tasi-RNA targets. We predicted more than 1300 unique target genes for the novel and conserved miRNAs. These target genes were predicted to be involved in different cellular, metabolic, response to stimulus, biological regulation, and developmental processes. Genome-wide identification of stress-responsive miRNAs, tasi-RNAs and their targets identified in this study will be useful in unraveling the molecular mechanisms underlying drought stress responses and genetic improvement of biomass production and stress tolerance in sorghum.


Introduction
Sorghum (Sorghum bicolor (L.) Moench) is the 4th most important nutritious cereal crop of the world after wheat, rice and maize. Sorghum is widely grown for food, feed, fiber and fuel and serve as a major staple food crop of millions of people, especially in semi-arid tropics including Africa, China, India, Mexico and USA. Drought is one of the major agronomic problems contributing to severe yield losses in sorghum worldwide. Further significant reduction in average yield has been reported due to the drought condition at flowering and grain filling phases, ultimately leading to a reduction in the grains number and size in sorghum. However, few genotypes of sorghum are relatively drought tolerant and known for its adaptation to water-limited environment (Ludlow and Muchow, 1990). Therefore, sorghum may serve as a model crop system for understanding the physiological and molecular mechanisms underlying drought tolerance. Dwindling fresh water resources and projected increase in incidences of drought under global climate change scenario warrant development of climate resilient sorghum varieties and hybrids. This necessitates a thorough understanding of the mechanisms of the physiological and molecular levels that contribute to drought tolerance. The interplay between drought and changes in gene expression has been studied only a few instances in sorghum (Buchanan et al., 2005). Therefore, molecular analysis of gene expression in sorghum needs further studies.
Several classes of small RNAs (sRNAs) regulate gene expression in plants. These include miRNAs (microRNAs), nat-miRNAs (natural antisense miRNAs), hc-siRNAs (heterochromatic small interfering RNAs), tasi-RNAs (transacting small interfering RNAs), nat-siRNAs (natural antisense small interfering RNAs), and ra-siRNAs (repeat-associated small interfering RNAs). These sRNAs are involved in growth and differentiation, genome stability, gene expression and defense in plants. miRNAs are coded by MIR (MICRORNA) genes in plants. MIR genes are transcribed by RNA Pol II. The 5 ′ methyl capped and 3 ′ polyadenylated pri-miRNA transcripts forms stem-loop secondary structures, and these processed into mature miRNAs of typically ∼21 nt length by DCL1-SE-HYL1 complex. The biogenesis of tasi-RNAs is triggered by the miRNA that targets the TAS (tasi-RNA producing locus) gene transcripts (Allen et al., 2005). Based on the similarity of the mature miRNA sequence, the miRNA (MIR) genes coding for identical or nearly identical mature miRNAs are grouped in to same family. For naming the miRNA in the same family usually zero to two mismatches are considered. Different MIR genes code for same family of miRNAs, the loci is named with the same number and Abbreviations: AGO, ARGONAUTE; ARF, auxin response factor; DCL, dicerlike protein; DCL1, endoribonuclease Dicer-like 1; DS, drought susceptible; dsRNAs, double-stranded RNAs; DT, drought tolerant; GO, gene ontology; MFEI, minimal folding free energy index; miRNA, microRNA; NGS, next-generation sequencing; PTGS, post-transcriptional gene silencing; RdRP, RNA-dependent RNA polymerase; RISC, RNA-induced silencing complex; RWC, relative water content; siRNA, small interfering RNA; sRNA, small RNA; tasi-RNA, trans-acting siRNAs; TFs, transcription factors; TPM, transcript per million genome-matched reads; UTR, untranslated region. a sequential alphabetical suffix (for more details please refer to Meyers et al., 2008).
Nearby 7000 plant miRNAs have been deposited in the miRBase 21, including 427 from Arabidopsis, 713 from rice, 321 from maize, 241 from sorghum and 525 from brachypodium. Compared with the number of miRNAs identified in rice, miRNAs reported in sorghum are very less. This suggests the potential for identification of many more miRNAs. Sorghum crop is highly tolerant to drought and heat stress, and the extension of MIR gene family members, including miR169 family, was suggested as one of the possible reasons for adaptation of sorghum to abiotic stresses (Paterson et al., 2009). Recently, small RNA sequencing approaches has been used to identify spatial and temporal expression pattern of miRNAs in sweet sorghum (Calviño et al., 2011;Zhang et al., 2011a). However, no drought stress regulated miRNAs and tasi-RNAs have been identified so far either through experimental or computational approaches in sorghum. To date, only four TAS families (TAS1, TAS2, TAS3, and TAS4) targeted by three miRNAs miR173, miR828, and miR390 have been identified in Arabidopsis. For the identification of genotype-specific and drought-responsive miRNAs, we constructed a small RNA library using RNA samples of two genotypes namely M35-1 (drought tolerant/DT) and C43 (drought susceptible/DS) grown under control and drought stress conditions. The miRNAs and tasi-RNAs were identified by sequencing the small RNA libraries using the Solexa deep sequencing technology combined with bioinformatics analysis. The comprehensive profile of miRNAs and tasi-RNAs, and their regulatory cascades in sorghum identified in this study will provide a basic platform for genetic improvement of sorghum.

Plant Materials and Stress Treatment
Sorghum [Sorghum bicolor (L.) Moench] seeds of two genotypes i.e., M35-1 (drought tolerant/DT) (Jogeswar et al., 2006) and C43 (drought susceptible/DS) (Mukri et al., 2010) were obtained from the Indian Institute of Millets Research (formerly known as Directorate of Sorghum Research), Hyderabad, India. Plants were grown at 28-32 • C day/night temperature with 12/12 h light/dark period in the National Phytotron Facility of Indian Agricultural Research Institute, New Delhi, India. Plants of both genotypes were irrigated alternately with water and Hoagland's solution at 3 days interval. Thirty days after sowing, drought stress was imposed by withholding irrigation to one set of plants (drought stressed) until the leaf relative water content reached to about 60-65% in both the genotypes. Thus the level of stress experienced by plants measured as RWC was similar in both the genotypes. The fully opened uppermost leaves from untreated (control) and treated (drought stressed) seedlings were harvested, frozen in liquid nitrogen and stored at −80 • C in the same day for the construction of small RNA libraries.

Small RNA Library Construction and Sequencing
Total RNA was isolated from the leaves using the TRIzol reagent (Invitrogen, USA), according to the manufacturer's protocol. The RNA quality was examined using gel electrophoresis Frontiers in Plant Science | www.frontiersin.org (28S:18S > 1.5) and Bioanalyzer (Agilent 2100, RIN = 8.0). Small RNA sequencing libraries were prepared using Illumina Truseq small RNA library preparation kit following manufacturer's protocol. Briefly, one microgram of total RNA was ligated with 3 ′ Illumina RNA adapter using Truncated T4RNA Ligase (NEB), followed by 5 ′ RNA adapter ligation and RT-PCR. The PCR amplicons were separated in a 6% PAGE along with a custom size ladder to select small RNA fraction. The final library was quantified using Agilent Bioanalyser DNA1000 chip and normalized to 10 nM. Seven pM concentration of the small RNA library was used for cluster generation and sequencing analysis (by Sandor Pvt. Ltd., Hyderabad, India) using the Illumina Genome Analyzer IIx according to the manufacturer's protocol (Illumina Inc., USA). The Illumina FASTQ files generated from this study have been submitted to the EMBL-EBI ArrayExpress (https://www.ebi.ac.uk/arrayexpress) with the accession number E-MTAB-1630.

Prediction of Conserved and Novel miRNA Members in Sorghum
The unique reads were submitted to the miRCat pipeline (UEA sRNA workbench V2.5.0-Plant version). The miRCat (miRNA Categorization) was run with the following parameters: The minimum sRNA abundance: 1 read; the minimum sRNA size: 16 nt; the maximum sRNA size: 30 nt, the minimum length of hairpin: 60 nt; the maximum number genome hits: 16. The 100 nt flanking regions of aligned reads were extracted from the genome and folded using RNAfold (http://rna.tbi.univie.ac. at/cgi-bin/RNAfold.cgi). The miRCat trimmed and analyzed the resulting secondary structure to verify the characteristic miRNA as per the plant miRNAs annotation criteria (Meyers et al., 2008) and executed the following additional checks: (1) The miRNA and miRNA * are derived from opposite stem-arms and form a duplex with two nucleotide 3 ′ overhangs. (2) The number of mismatches between miRNA and miRNA * should be less than or equal to four. (3) The frequency of asymmetric bulges is one or less and size of the bulge is no more than 2 nucleotides within the miRNA/ miRNA * duplex. The maximum number of occurrences (reads) for a particular miRNA family was denoted as miRNA abundance. The folding structure of precursors was examined using RNA folding/annotation tool (UEA sRNA workbench V2.5.0-Plant version) that uses the Vienna Package to obtain the secondary structure of a precursor sequence and highlighting up miRNA/miRNA * sequences on hairpin structure.

Differential Expression Analysis of miRNAs
To compare abundance of miRNAs in control and treatment library, the count of each miRNA was normalized to transcripts per million (TPM) following appropriate statistical method (Audic and Claverie, 1997). (i) Normalization criterion: TPM = (actual miRNA count/total count of clean reads) * 1,000,000. Afterwards, the fold-change and P-value of the normalized expression were calculated by using the following formula, (ii) Fold-change criterion: Fold change = log 2 (miRNA TPM in the treatment library/miRNA TPM in the control library), (iii) P-value: The P-value of precursor miRNA candidate was tested using randfold (using a cutoff of 0.1) tool (Bonnet et al., 2004) and integrated on the UEA sRNA workbench. The Poisson distribution model was used for estimating the statistical significance of miRNA expression changes under control and treatment conditions. Upregulation of any miRNA expression levels was considered a positive value, while negative values indicate down-regulation. To identify drought-responsive miRNAs, several standards (Eldem et al., 2012) were followed as: (1) normalized count was at least 1 TPM in either control or stress library; (2) P-value = 0.01 as the threshold; (3) log 2 ratio of the normalized count under stress or control libraries was >1 or < −1. Unique (active form of the miRNA) or identical mature miRNAs, generated from two or more homologous miRNA genes were only considered for differential expression analysis.

Tasi-RNA Analysis
To detect phased small RNA clusters corresponding to tasi-RNAs, the "The UEA Small RNA Workbench V2.5.0" (http:// srna-workbench.cmp.uea.ac.uk/tools/ta-si-prediction/) (Stocks et al., 2012) was used. Small RNAs (sRNAs) that are not identical to the genome were rejected. A minimum abundance 2 and p-value threshold of 0.0001 was used to detect statistically significant clusters. Potential phase-initiators for these TASs were predicted by using psRobot program (http://omicslab.genetics. ac.cn/psRobot/)  and psRNATarget program (http://plantgrn.noble.org/psRNATarget/) (Dai and Zhao, 2011) with default parameters and an assumption that the 10th nucleotide position on the sRNA serving as the phase-initiator corresponds to the cleavage start position of its targeted TAS. The overlapping TAS on chromosome was merged as a single cluster. The identical tasi-RNAs within TAS gene were removed and retain only unique tasi-RNAs for further analysis.

Target Gene Prediction
The potential targets of sorghum miRNAs were predicted using the psRobot; plant small RNA analysis toolbox (http://omicslab. genetics.ac.cn/psRobot/)  with strict parameters and psRNATarget (http://plantgrn.noble.org/psRNATarget/) (Dai and Zhao, 2011) with default parameters. For psRobot program, we used the following parameters for the target prediction-Penalty score threshold: 3.0; Five prime boundary of essential sequence: 1; Three prime boundary of essential sequence: 31; Maximal number of permitted gaps: 0; Position after which with gaps permitted: 1. For psRNATarget program, we used the following standards for the target prediction-Maximum expectation: 3.0; Length for complementarity scoring (hspsize): 20; Target accessibility (UPE): 25.0; Flanking length around the target site for target accessibility analysis: 17; Range of central mismatch leading to translational inhibition: 9-11 nt. All the mutual and unique targets were accepted and listed in this article. Further, all the targets regulated by sorghum miRNAs identified in this study were subjected to AgriGO toolkit (http:// bioinfo.cau.edu.cn/agriGO/) (Du et al., 2010) to investigate gene ontology enrichment. The singular enrichment analysis (SEA) was performed to find enriched GO terms within annotated miRNA targets.

RNA-Seq Library Construction, Sequencing, and Data Analysis
One microgram of total RNA was used for the preparation of RNA-Seq library using Illumina TruSeq mRNA library preparation kit following the manufacturer's protocol. In short, poly-A RNA was isolated from total RNA and chemically fragmented. First and second strand synthesis were followed by end repair, and adenosines were added to the 3 ′ ends. Adapters were ligated to the cDNA and 200 ± 25 bp fragments were gel purified and enriched by PCR. The libraries generated were quantitated using an Agilent Bioanalyzer (DNA 1000 kit; Agilent Technologies, Santa Clara, CA) and a 2 × 101 cycle paired end sequencing (sequenced by Sandor Pvt. Ltd., Hyderabad, India) was performed using an Illumina HiScanSQ sequencer (Illumina Inc.). Initially, raw reads were processed by the NGSQC toolkit (http://www.nipgr.res.in/ngsqctoolkit.html) and high quality reads were subjected to de-novo assembly using Trinity assembler (Patel and Jain, 2012). Assembled transcripts were quantified by standard pipeline (Trinity→RSEM→R→DESeq), and those transcripts were removed, which has zero FPKM in all four samples (Anders, 2010;Grabherr et al., 2011;Li and Dewey, 2011). These transcripts were further processed by the transdecoder tool to retrieve the full length coding sequence and subsequent annotated by FastAnnotator (http://fastannotator. cgu.edu.tw/) . Pathway enrichment analysis was performed for the predicted transcripts by KEGG Automatic Annotation Server (KAAS; www.genome.jp/tools/kaas/) for the classification of spatial and temporal governed pathways. The Illumina FASTQ files generated from this study have been submitted to the EMBL-EBI ArrayExpress (https://www.ebi.ac. uk/arrayexpress) with the accession number E-MTAB-3571.

Analysis of Small RNA Population in Sorghum
To identify drought stress regulated miRNAs from sorghum, four libraries of small RNAs from two genotypes namely M35-1 (drought tolerant or DT) and C43 (drought susceptible or DS) grown under irrigated and drought stress conditions were constructed and sequenced independently. The relative water content (RWC) of drought stressed seedlings were about 65% in both genotypes ( Figure S1 in Supplementary Material), indicating similar level of stress imposed to these genotypes. The Solexa (Illumina) sequencing of small RNA libraries from M35-1 and C43 led to the generation of 19,821,595 and 18,890,943 primary reads, respectively, under well irrigated condition, whereas 22,878,203 and 21,740,738 primary reads, respectively, were generated under drought stress conditions. Initially, the raw reads were filtered for possible adaptor contaminations and subsequently mapped to the sorghum genome. For each library, more than 90% reads were mapped to the genome, suggesting slight contamination in library construction and Illumina sequencing. The remaining raw reads (after adaptor removal) were further filtered for poor quality sequences, invalid sequences and sequences smaller than 16 nt and larger than 30 nt and as a result, a total of 9,472,887 clean reads were obtained which represented 975,457 unique reads under irrigated condition. Similarly, a sum of 7,296,968 clean reads with 944,141 unique reads were also determined from small RNA libraries of drought stress treated seedlings (Table S1 in Supplementary Material). The unique reads (after filtering out for t/rRNA matches) from each library that perfectly mapped (with no mismatch) to the sorghum genome represented the small RNA (sRNA) population.

Identification of Conserved and Novel miRNAs
To predict miRNAs, all unique reads (gained after filtering) were submitted to the miRCat pipeline (UEA sRNA workbench V2.5.0-Plant version) and were mapped (with no mismatch) against the reference genome of sorghum. The 100 nt flanking regions of aligned reads were extracted from the genome and folded using RNAfold. The resulting secondary structures (potential precursor) were then trimmed and analyzed by miRCat. To recover the putative miRNA candidate, the potential precursors were subjected to a series of stringent criteria suggested for the annotation of plant miRNAs (Ambros et al., 2003;Meyers et al., 2008;Kozomara and Griffiths-Jones, 2013). The probable miRNA candidates that perfectly matched (miRNAs with ≤3 mismatches) with mature miRNAs of sorghum in the miRBase 21 were acknowledged as known miRNAs. The sequences that matched with miRBase entries of other plant species were designated as conserved miRNAs. Finally, the sequences that showed no homology to any previously known and conserved plant miRNAs were denoted as novel miRNA in sorghum. In this study, 97 miRNAs belonging to 26 miRNA families were found identical with the known miRNAs of plant species in miRBase 21. Among these, 32 miRNAs were perfectly matched with mature miRNAs of sorghum and the remaining 65 miRNAs were highly conserved in other plant species and hence are referred to as known and conserved miRNAs, respectively (Table S2 in Supplementary Material). In addition, a total of 526 miRNAs derived from 518 loci identified in this study did not show sufficient homology with any of previously reported plant miRNAs in miRBase 21, and hence classified as novel miRNAs (Table S2 in Supplementary Material). The identification of miRNA * candidates provided further support to consider them as true miRNAs. In this study, 47 and 130 miRNA * were predicted from 97 conserved and 130 novel miRNAs, respectively. A close observation revealed that 541, 499, and 274 reads of miR167h * , miR169e * , and miR167d * accumulated, respectively. Likewise, 117, 63, and 47 reads of novel miRNAs namely novel-sbi-miR-383 * , novel-sbi-miR-77b * , and novel-sbi-miR-51 * accumulated, respectively (Table S2 in Supplementary Material). Furthermore, minimal folding free energy index (MFEI), a sufficient criterion to distinguish miRNA from coding mRNAs and non-coding RNA, i.e., tRNAs and rRNAs, defines that a candidate RNA sequence is more likely to be a miRNA when the MFEI is greater than 0.85 (Zhang et al., 2006). In this study, 96.91% conserved (except for miR166f, miR396f, and miR399b) and 58.56% novel pre-miRNAs of sorghum had an MFEI = 0.85. Additionally, nucleotide bias analysis showed that uracil was the most prominent nucleotide at the beginning in 84 (86.60%) conserved and 177 (33.65%) novel miRNAs. The pre-miRNAs of novel and conserved miRNAs bore a canonical stem-loop structure with free energies ranging from −16 to −131.8 kcal·mol −1 (average of −54.29 kcal mol −1 ). The chromosomal distributions of predicted miRNAs are listed in Figure 1.

Drought-responsive and Genotype-specific miRNA
We compared the normalized count of each miRNAs in a stressed library against the control library to compute the stress regulated expression of miRNAs. The number of mature miRNAs per million clean reads (known as transcripts per million or TPM) FIGURE 1 | Distribution of predicted miRNAs. Chromosomal distribution (A) showed that the maximum numbers of miRNAs were predicted from chromosome 1 and 4. Distribution of conserved miRNAs (B) showed that miR169, miR166, and miR167 were the most abundant miRNAs in sorghum genome.
was calculated to identify the normalized expression level. A total of 96 unique miRNAs (fold change = 2), belongs to 8 known and 88 novel families, showed differential expression under drought stress as compared to the control, in at least one genotype (Figure 2; Table 1; Table S3 in Supplementary Material). Out of the 96 drought stress regulated miRNAs, 23 and 9 miRNAs were up-and down-regulated, respectively, under the drought stress, in both genotypes. Forty four miRNAs were upregulated in drought tolerant M35-1, while these miRNAs were downregulated in drought sensitive C43 genotype under drought stress. In contrast, 19 miRNAs were downregulated in drought tolerant M35-1, while these miRNAs were upregulated in drought sensitive C43 genotype under drought stress. The novel-sbi-miR-259 was undetectable in M35-1, while it was downregulated in C43 (Figure 2; Table 1; Table S3 in Supplementary Material).

Targets of Known and Novel miRNAs
Recently, 64 target genes for 17 miRNAs (Jiangfeng et al., 2010), 125 target genes for 42 miRNAs (Zhang et al., 2011a) and 72 potential target genes for 31 miRNAs (Katiyar et al., 2012) were reported in sorghum. Here we identified more than 1300 unique potential targets for 49 conserved and 383 novel miRNA families in sorghum (Table S4 in Supplementary Material). For the remaining 96 new miRNAs, no target could be identified, which might due to the stringency of target prediction used in this study. The further BLAST analysis identified miRNA targets homologous to conserved target genes of several plant species. The putative target genes were considered to be a key factor in a wide range of biological processes. Inconsistent with previous reports (Rhoades et al., 2002;Bartel and Bartel, 2003;Song et al., 2011), many of the miRNA target genes predicted in this study encode transcription factors belonging to SPB, Zinc finger, WRKY, WD-40, NAC, MYB, HSFs, GRAS, ARFs, and bHLH families. Several targets for novel miRNAs identified in this study were genes with unknown function. Further functional analysis of miRNA-target gene pair will contribute to our understanding of the role of these miRNAs in sorghum.

Go Enrichment Analysis of miRNA Target Genes
Widely used standard for functional annotation and enrichment analysis through gene ontology (GO) was carried out for more than 1300 unique potential targets of 472 unique sorghum miRNAs using AgriGO database (Du et al., 2010) with default parameters. Significantly high percentages of these targets were found to be involved in cellular, metabolic, response to stimulus, biological regulation, and developmental processes (Figure 3). Further dissection of "response to stimulus" led to the identification of seven genes with over-represented GO term "response to water deprivation" (GO: 0009414), one with "response to water stimulus" (GO: 0009415), five with "response to heat" (GO: 0009408) and two with "response to heat acclimation" (GO: 0010286) ( Figure S3 in Supplementary Material; Table S5 in Supplementary Material). To explore the modulated biological processes between tolerant and sensitive genotypes under drought stress, differentially expressed genes regulated by miRNAs specific to drought tolerant and sensitive genotypes were analyzed separately using GO term enrichment analysis. Gene ontology enrichment analysis was performed using an FDR-adjusted p-value of ≤0.05 to select pathways that were statistically enriched with miRNA targets. The significantly enriched GO terms, including response to abiotic stimulus, response to stress, response to stimulus, response to external stimulus, response to organic substances, response to starvation, response to hormone stimulus, response to gravity, response to endogenous stimulus, response to chemical stimulus, and response to auxin stimulus were found among the up-or down-regulated genes in at least one genotype of sorghum (Table S5 in Supplementary Material). Interestingly, three stressrelated biological process GO terms varied significantly between tolerant and sensitive genotypes under drought stress. We noticed that GO terms "response to inorganic substance" (GO: 0010035, FDR p-value = 0.038) and "response to stress" (GO: 0006950, FDR p-value = 0.024) were significantly enriched only in sensitive genotype, whereas "response to abiotic stimulus" (GO: 0009628, FDR p-value = 0.045) only in tolerant genotype under drought stress. The results indicate that the inherent preparedness and responsiveness of tolerant cultivar toward drought stress is much higher as compared with sensitive cultivar.

Digital Expression Analysis of miRNAs and their Targets
To investigate the relationship between predicted miRNAs and their targets, expression levels were calculated for randomly selected 20 miRNAs (eight conserved and 12 novels) and their corresponding target genes from small RNA and RNA-Seq sequencing experiments, respectively. A negative correlation (miRNA up, target gene down, and vice versa) was observed between the expression levels of several up-or down-regulated miRNAs and their targets (Figure 4). For instance, the induced expression of miR156b, miR396b-c, miR396d-e, miR396f, miR5385, novel-sbi-miR-46, novel-sbi-miR-48, novel-sbi-miR-141, novel-sbi-miR-176, novel-sbi-miR-224, and a novel-sbi-miR-335 resulted in enhanced level of accumulation of their target genes under drought stress in DT-genotype. Conversely, drought-induced down-regulation of miR156a, miR319a-b, miR529, novel-sbi-miR-111, novel-sbi-miR-120a-b, novel-sbi-miR-227a-c, novel-sbi-miR-268, novel-sbi-miR-376, and novelsbi-miR-350 led to up-regulation of their target gene. The lists of differentially expressed target genes under drought stress are listed in Table S6 in Supplementary Material.

High-throughput Sequencing of Sorghum Small RNAs
In the past few years, high-throughput sequencing is being used to identify miRNAs in several crop plants, including barley, brassica, cowpea, cucumber, maize, peanut, rice, sorghum, soybean, and wheat. In sorghum, 17 (Du et al., 2010), 29 (Zhang et al., 2011a) and 31 (Katiyar et al., 2012) new miRNAs were identified recently. Compared with the number of miRNAs that have been identified in other cereal crops such as rice, only limited miRNAs have been discovered in sorghum. Additionally, drought-regulated miRNAs have not been identified either through experimental or computational methods in sorghum.
In the present study, we constructed and sequenced sRNA libraries from seedlings of M35-1 (drought tolerant) and C43 (drought susceptible) genotypes grown under irrigated and drought stress conditions. As a sequencing throughput, the small RNAs population of size 17-29 nt have been found. After eliminating the adapter sequences, the highest read abundance was found for 21-24 nt small RNAs, which is consistent with the size of the DCL enzyme products (Figure 5). The read abundance for each miRNA families highly differed. The miR166, miR167, miR156, and miR399 are the largest miRNA families with 11, 10, 8, and 7 members, respectively, in sorghum. The miR160 and miR396 families have 6 members each (Figure 1; Table S2 in Supplementary Material). Likewise, miR156g * and miR398 * were found to accumulate at high levels under drought condition in DS and DT genotypes of sorghum, suggesting that these two miRNA * might function in stress response mechanism irrespective of genotypes. The miR166f * , miR167g * , miR169e * , miR398 * , and novel-sbi-miR-383 * accumulated at high levels exclusively in DS genotype, whereas, miR166g * , miR167h * , and miR169h * accumulated at high levels exclusively in DT genotype, suggesting the function of these miRNA * in a genotype specific manner. We also observed that the majority of the miRNA transcripts showed higher abundance in tolerant genotype as compared to susceptible genotype under drought treatment (Table S4 in Supplementary Material).

Monocot-and Dicot-abundant miRNAs
Previous studies have shown that some miRNAs are conserved across the plant kingdom, while other miRNAs have been reported in either monocot or dicot plants. In this study, 49 conserved miRNAs belonging to 26 miRNA families were identified from previously deposited miRNAs in miRBase 21 ( Table 2). Among these, 30 (61.22%) miRNA families (denoted with the symbol "&" in Table 2) were found to be highly conserved between dicot and monocot plants. The remaining 19 (38.78%) miRNA families (denoted with the symbol "#" in Table 2) were found conserved in monocot plants only. In this study, miR156g-h was identified first time in monocot, which was previously reported only in dicot plant, suggesting the conserved nature of this family in monocot and dicots. Likewise, variants of selected miRNA members, namely miR156b, miR169b, miR396a, and miR399a were found in sorghum as reported previously in monocots. However, other variants of the same family were found to be conserved in dicot plants. This suggests that these miRNAs might have evolved after the divergence of monocots and dicots. In the present study, we found monocot abundant miR156g-h, miR166c-i, miR167f-j, miR169b, miR171a-b, miR172a-c, miR319a-b, miR395, miR396a, miR437, miR529, miR2118a-b, miR2118d, miR2118e, miR2275, miR5385, and miR6221 (denoted with the symbol "@" in Table 2) which were not previously reported in sorghum. Similarly, miRNA families, namely miR437, miR156g-h, miR2118d, miR5385, and miR6221 (denoted with the symbol "$" in Table 2) were identified first time in monocots. In addition, we found that the expression levels of conserved miRNAs were higher than that of unique miRNAs. For example, conserved miRNAs such as miR156c-f, miR156g-h, miR160b-f, miR166c-i, miR166j-k, miR167d-e, miR167f-j, miR396a, and miR398 exhibited high levels of read abundant (more than 750 TPM) in all four libraries. These observations support the earlier conclusion that phylogenetically conserved miRNAs are highly abundant (Sunkar et al., 2008).

Response of Known and Novel miRNAs to Drought Stress
We compared the small RNA expression profiles of drought tolerant and sensitive genotypes of sorghum and identified 96 drought-responsive miRNAs with more than twofold change at-least in one genotype ( Table 1; Table S3 in Supplementary Material). Of the 96 drought regulated miRNAs, 63 miRNAs showed opposite regulation in drought tolerant and drought sensitive genotypes, i.e., 44 of them were upregulated in M35-1 but downregulated in C43, while 19 of them were downregulated in M35-1 but upregulated in C43 (Figure 2; Table 1, Table S3 in Supplementary Material). Among the conserved miRNAs, miR160a, miR169d-l, 396b-c, 396d-e, miR529, miR2118e, miR2275, and miR5385 were found as drought stress responsive in sorghum. The miR160 is known to regulate ARF10/ARF16/ARF17 repressors family and overexpression of miR160 confers auxin hypersensitivity (Turner et al., 2013). In Arabidopsis, ARF17, a target of miR160, negatively regulate the expression of auxin inducible Gretchen Hagen3 (GH3) genes, encoding acyl-acid-amido synthetase which fine-tuning adventitious root initiation in the Arabidopsis (Gutierrez et al., 2012). In cereals, adventitious (nodal/crown) roots, contributes significantly to soil moisture uptake under drought stress (Rostamza et al., 2013). In the present study, we also found that miR160a targets several ARF family members, including ARF10, ARF16, and ARF17. Thus, miR160 upregulation in M35-1 and downregulation in C43, might be contributing to better tolerance of M35-1. We also observed that the downregulation of miR396b-c and miR396d-e in the DS genotype of sorghum as reported previously in rice  and cowpea (Barrera-Figueroa et al., 2011), but were upregulated in the DT genotype of sorghum as reported previously in Arabidopsis (Liu et al., 2009) and tobacco (Feng-Xi and Di-Qiu, 2009) under drought-stress. Previous studies showed that miR396a-overexpressing transgenic plants were more drought tolerant than wild-type plants (Liu et al., 2009;Yang et al., 2009). Thus, upregulation of miR396 in DT genotype appears to be contributing to drought stress tolerance of M35-1. Likewise, downregulation of miR529 in both genotypes of sorghum was observed as previously shown in rice Jeong et al., 2011). In addition to the conserved miRNA family members, we found drought regulation of several members of novel miRNA families. These miRNAs might be involved in lineage-or species-specific stress response pathways and functions. In general, many of the miRNA families consist of more than one miRNA gene, which may have identical or diverse mature miRNA sequences. These homologous miRNA genes may functionally diverge from each other during the evolutionary process. For instance, in our study, miR156 and miR164 families showed clear evidence for functional diversification. While one member of miR156 family known as miR156b was induced by drought stress, while another member miR156a was significantly down-regulated in drought sensitive genotype. On the other hand, all members of the miR396 family were up-regulated by drought stress in the tolerant genotype, but not in the sensitive cultivar. Detailed analyses revealed that the majority of the sorghum miRNAs was expressed in a genotype specific manner.

MicroRNAs Regulated Targets in Response to Drought Stress
In general, plants respond to environmental stresses by regulating target genes through up-or down-regulation of miRNAs to cope with these stresses. In the present study, 108 targets were predicted for 49 known miRNAs using bioinformatics analysis, and many of these were conserved in plant species, indicating broad conservation of the known miRNA regulatory roles in    Table S4 in Supplementary Material). In addition to the well-documented conserved targets, a few of the known miRNAs, including miR2118c, miR529, and miR399a were found to target additional genes in sorghum that have not been previously reported. Selected miRNAs, such as miR529 (targeting genes coding for SBP, DCD, cellulase, protease-related, and ubiquitin), and miR398 (targeting genes coding for Cu/Zn SOD, selenium-binding protein, and cytochrome C) showed multiple target genes, indicating that these miRNAs have diverse roles. Furthermore, a single gene may be regulated by several miRNAs. For example, squamosa promoter binding protein (SBP) gene is targeted by miR156a, miR156b, miR156c-f, miR156g-h, miR529, novel-sbi-miR-119, novel-sbi-miR-383, novel-sbi-miR-329, and a novel-sbi-miR-350 in sorghum. In addition, several conserved members of identical miRNA families were found to possess conserved target genes. For instance, all members from the sbi-miR156 and sbi-miR160 family targeted to SBP and auxin response factor (ARF), respectively. Similarly, all members of the sbi-miR399 and sbi-miR164 families targeted phosphate transporter (PHT), and NAC domain containing protein, respectively. Similar results were observed previously in several plant species and these miRNA target genes have been found to be involved in plant growth and/or responses to environmental changes (Unver and Budak, 2009;Xie et al., 2011). The putative targets of the drought regulated miRNAs offered important clues on the drought response in sorghum. For instance, sbi-miR164 (targeted to NAC transcription factor) was found to be downregulated by drought stress at-lest in one genotype of sorghum. This is consistent with previous reports in Arabidopsis and sugarcane (Guo et al., 2005;Raman et al., 2008;Ferreira et al., 2012). NAC transcription factors regulate the development, growth and stress responses, including cold, drought and pathogen attack (Kikuchi et al., 2000;Ooka et al., 2003). Similarly, several droughts regulated miRNAs such as novel-sbi-miR-4, novel-sbi-miR-41, novel-sbi-miR-87, novel-sbi-miR-391, and novel-sbi-miR-412, which target ARFs, were upregulated in DT M35-1 but down regulated in DS C43 sorghum genotype under drought stress. This is consistent with earlier studies that showed stress adaptation by ARFs regulated auxin-mediated mechanisms (Guilfoyle and Hagen, 2007). Proper regulation of auxin transport is critical for drought tolerance (Remy et al., 2013). These results suggest that the predicted targets such as NAC and ARFs play important roles in the drought response in tolerant genotype of sorghum. The suppression of novel-sbi-miR-335 in the sensitive genotype under drought stress could be one of the factors associated with susceptibility of sorghum cultivar C43. The induction of novel-sbi-miR-335 in tolerant genotype under drought stress might contribute to drought-tolerance of the M35-1 cultivar. Moreover, the up-regulation of novel-sbi-miR-389, responsive to dehydration stress in drought tolerant M35-1 cultivar was consistent with earlier observations in barley (Kantar et al., 2010) and soybean (Kulcheski et al., 2011). Some miRNAs such as novel-sbi-miR-360a-c showed up-regulation in drought sensitive C43, but showed down-regulation in drought tolerant M35-1 indicating different adaptive mechanisms by respective genotype. Interestingly, the novel-sbi-miR-266 and a novel-sbi-miR-339 miRNA family showed decreased expression levels under drought stress in both the genotypes. These miRNAs may control a genotype-independent common drought response mechanism. All these miRNA members putatively target a zinc finger protein. Downregulation of novel-sbi-miR-112 (targeted to "basic region/leucine zipper protein-60" or bZIP60) and novelsbi-miR-254 (targeted to "homeobox-leucine zipper protein-17" or HD-ZIP17) in response to drought may increase the abundance of bZIPs, and contribute to drought tolerance in the DS genotype of sorghum. The importance of bZIPs in stress tolerance of plants was also reported previously Golldack et al., 2011). Moreover, miR5385 was up-regulated by drought in the tolerant cultivar, but downregulated in sensitive cultivar. The predicted target for this miRNA is a basic-helix loop-helix (bHLH) transcription factor. Several reports elucidated the role of bHLH in response to abiotic stresses, e.g., freezing (Chinnusamy et al., 2003), iron deficiency (Long et al., 2010) and salt stress (Li et al., 2010a). Abiotic stresses lead to accumulation of excess concentrations of reactive oxygen species (ROS), resulting in oxidative damage to the cell. Peroxidases help restrict ROS build up and thus oxidative damage to the cells. We predicted peroxidise family as target for six miRNA families, namely novel-sbi-miR-30, novel-sbi-miR-106, novel-sbi-miR-177, novel-sbi-miR-204, novel-sbi-miR-272, and a novel-sbi-miR-217 in sorghum. In tolerant genotype, the expression level of novel-sbi-miR-272 was lower than that of sensitive cultivar during drought stress. The increase of novelsbi-miR-272 in the sensitive genotype C43 under drought stress may lead to a decrease in the transcript levels of peroxides and thus could be one of the factors associated with vulnerability of this genotype. Several droughts-responsive miRNA families such as novel-sbi-miR-105a-b, novel-sbi-miR-180a-c, and novel-sbi-miR-416 targeted to Kelch repeat-containing F-box protein in sorghum. These proteins are known to be involved in response to biotic and abiotic stresses (Sun et al., 2010;Jia et al., 2012). In sorghum, we observed similar expression patterns for these miRNAs irrespective of genotypes. Three miRNA families (e.g., Novel-sbi-miR-213, miR169c, and a miR169d-l) were down-regulated by drought stress in sensitive genotype C43. These miRNAs targeted to nuclear factor Y (NFY) transcription factor. Li et al. (2008) reported that NFYA5 transcript is strongly induced by drought stress in an abscisic acid (ABA)dependent manner, whereas miR169 was down-regulated by drought stress through an ABA-dependent pathway. Transgenic plants overexpressing miR169 and NFYA5 knockout plants exhibited hypersensitivity to drought stress. Thus, NFYA5 is important for drought tolerance in plants (Li et al., 2008). In contrast, tomato plants over-expressing miR169c exhibited better tolerance to drought due to reduced stomata opening (Zhang et al., 2011b). Thus, drought tolerance may likely involve divergent mechanisms in the different plants. Plants respond to heat or drought stress by the induction of the synthesis of heatshock proteins (HSPs). HSPs play a crucial role in protecting plants under diverse abiotic stresses (Wang et al., 2004). Four miRNA families (e.g., miR396d-e; novel-sbi-miR-26; novel-sbi-miR-85a-k, and a novel-sbi-miR-336) that targeted to HSPs/HSFs were found to be down-regulated during drought stress in the drought sensitive genotype of sorghum. In contrast, miR396d-e and novel-sbi-miR-336 were found to be up-regulated under drought stress in the drought tolerant genotype of sorghum. In addition, GO analysis exposed the response of several target genes in heat (GO: 0009408), water deprivation (GO: 0009414), water stimulus (GO: 0009415) and heat acclimation (GO: 0010286). The miR156b (targets IRREGULAR XYLEM 1), miR396b-c (targets mitochondrial HSP70-2); novel-sbi-miR-27 (targets HSP21), novel-sbi-miR-55 (targets HOMEBOX 7), novel-sbi-miR-171 (targets early responsive to dehydration 15), novel-sbi-miR-254 (targets nitrate transmembrane transporter), novel-sbi-miR-268 (targets HSP101), novel-sbi-miR-272 (targets HSP1; WRKY; Heat-intolerant 1 and Cytochrome P450), novelsbi-miR-277 (targets F-box family protein), novel-sbi-miR-329 (targets abscisic acid responsive elements-binding factor 2), and novel-sbi-miR-389 (DEHYDRIN XERO1), were observed as a regulator of drought and heat responsive proteins (Table S5 in Supplementary Material). The novel-sbi-miR-268 and a novelsbi-miR-272 was found to be up-regulated in sensitive genotype, whereas, down-regulated in the tolerant genotype of sorghum. Likewise, miR396b-c and novel-sbi-miR-254 was found to be down-regulated in sensitive, but up-regulated in the tolerant genotype of sorghum. These results suggested that the two genotypes studied here have differential molecular mechanisms to respond to drought stress. In addition to conserved targets, genotype-specific miRNA regulated gene targets were also observed ( Table S4 in Supplementary Material). Although many newly evolved miRNAs that may exhibit weak expression, imperfect processing and lack of targets are believed to serve no biological function, many of them have been shown to target and regulate specific genes or gene families in sorghum (Li et al., 2010b;Pantaleo et al., 2010;Zhang et al., 2010).

Conserved miR390-TAS3 Pathway in Sorghum
Tasi-RNAs were originally discovered in Arabidopsis, where four TAS families (e.g., TAS1/2, TAS3, and TAS4) code for TAS transcripts which are targeted by three miRNA families, namely miR173, miR390, and miR828. The tasi-RNA prediction approach includes the detection of phased 21 nt sRNAs, characteristic of tasi-RNA loci and the assessment of statistical significance in term of P-value. Using this approach, various research groups have demonstrated the existence of tasi-RNAs in Arabidopsis, rice, grapevine, brassica, apple and soybean. To examine whether there are any homologous genes of predicted TASs in sorghum, we aligned 31 TASs with non-redundant (nr) nucleotide database. Surprisingly, we did not find homologs of 29 TASs. However, two miR390 targeted TASs (Table S7 in Supplementary Material) showed significant similarity with TAS3 gene of Sorghum bicolor and Saccharum officinarum. Further, these two TASs were aligned with Arabidopsis TAS genes, and found that sorghum TAS is similar to that of AtTAS3 gene, hence named these as sorghum TAS3 genes. Thus, the miR390-TAS3 pathway is conserved in sorghum also in addition to the previously reported plants such as Arabidopsis, rice, grapevine, brassica, apple and soybean. Conversely, we did not find miR828-TAS4 and miR173-TAS1/TAS2 pathway in sorghum. In consistency with previous observations that TAS1, TAS2, and TAS4 families do not occur in monocots (Axtell et al., 2006), in our study, we also did not find these tasiRNAs. miR828-TAS4 pathway is widely represented in dicot species (Luo et al., 2012), while the miR173-TAS1/TAS2 pathway has been found to be restricted in Arabidopsis. We mapped TAS loci with coding or non-coding regions, and found that 11 TASs did not have homology to coding regions, 18 TASs showed partial homology to coding regions, and two TASs showed significant similarity (>50% query coverage and sequence identity) with coding regions. In addition, miR5567 and miR6225 target TASs were found to be highly conserved ( Figure S2 in Supplementary Material). We also observed that two putative TAS3 genes have similar sequences, but different predicted structures. Previous reports disclosed the connection of Arabidopsis tasi-RNAs (TAS1, TAS2, and TAS3) with hypoxia stress (Moldovan et al., 2009). In addition, TAS4-derived siRNAs regulate the biosynthesis of anthocyanins in response to phosphorous deficiency (Hsieh et al., 2009;Luo et al., 2012). Likewise tasi-RNAs from Pinus taeda was predicted to target the transcripts of two disease resistance proteins of pine, suggesting its role in the response to pathogens (Lu et al., 2007). However, no drought-stress related tasi-RNAs have been reported in plants. In the present study, we noticed that most of miR390derived tasi-RNA sequences were present in both DS and DT libraries under control and stress conditions. We assume that neither miR390 nor tasi-RNAs from TAS3 genes were affected by drought stress. Furthermore, we noticed that all tasi-RNAs derived from miR5567 targeted TASs, and miR6230 targeted TASs were found only in the control library of drought sensitive genotype. Targets of miRNAs and tasi-RNAs play important roles in many pathways involved in the development or responses to the environment. Our result revealed that 5 ′ cleavage product of miR390-TAS3 produces two tasi-RNAs [e.g., SbiTAS3a-siR390-5 ′ D4 (+); sbiTAS3b-siR390-5 ′ D3 (+)] that target transcripts of auxin response factors (Figure 6). This is consistent with previous reports in Arabidopsis (Allen et al., 2005). Similarly, miR5568-TAS derived tasi-RNA [e.g., SbiTAS-siR5568c-3 ′ D2 (+)], found in drought-stressed library of drought sensitive genotype, targets to universal stress protein (USP). Genes encoding the universal stress protein domain confer enhanced ability of plants to tolerate stress. The other tasi-RNAs are novel, which will help further broaden the scope of tasi-RNA mediated gene regulation.

Conclusions
Small RNA-seq in combination with bioinformatic analysis identified 97 conserved and 526 novel miRNAs in sorghum. In addition, we discovered several novels and conserved miRNAs regulated by drought stress. Common as well as genotypespecific miRNA expression patterns were discovered from DS and DT-genotypes of sorghum elucidated the underlying molecular mechanisms and diverse physiological pathways. The targets for predicting miRNAs were found to be involved in cellular, metabolic, response to stimulus, biological regulation, and developmental processes. Additionally, the discovery of two TAS3 genes (orthologs of Arabidopsis TAS3 genes) targeted by miR390 revealed conserved miR390-TAS3 pathway in sorghum. The prediction of miR390-TAS3 derived tasiARF suggested their role in the auxin signaling pathway. The outcomes of this study provide valuable information for further functional characterization of miRNAs in response to drought stress in sorghum. Our findings laid the groundwork for functional characterization of droughtresponsive miRNAs and manipulating miRNAs or their targets for improving biomass production and stress tolerance in sorghum.

Author Contributions
AK initiated the research, performed the downstream analyses, interpreted the results and drafted the manuscript. SS helped in computational analyses and data management. SM conducted the plant stress treatment and collected the tissues. VC interpreted the results and designed the experiments. DP and KB conceived the idea of the study and drafting the manuscript. All authors have read and approved the manuscript for publication.