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 (trans-acting 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 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.

Materials and Methods

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(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.

Bioinformatics Analysis of sRNA Sequences

Adaptor and low-quality sequences were trimmed as suggested by Sunkar et al. (2008). Sunkar et al. (2008) by using the “sequence file pre-processing tool” from UEA sRNA workbench V2.5.0-Plant version (http://srna-tools.cmp.uea.ac.uk/) (Stocks et al., 2012). High quality trimmed sequences (reads with no “N,” no more than 6 bases with quality score <13) with a length of 16–30 nt were further subjected to remove, if mapped with plant t/rRNAs from “Rfam” (excepted miRNAs), Arabidopsis tRNAs from “The Genomic tRNA Database,” and plant t/rRNA sequences from “EMBL release 95” by using the “filter tool” from UEA sRNA workbench V2.5.0-Plant version.

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₂ (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₂ 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/) (Wu et al., 2012) 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/) (Wu et al., 2012) 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/) (Chen et al., 2012). 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.

Results

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⁻¹ (average of −54.29 kcal mol⁻¹). The chromosomal distributions of predicted miRNAs are listed in Figure 1.

FIGURE 1

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.

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) 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).

FIGURE 2

Figure 2. The Venn diagram illustrates the numbers of common and unique differentially expressed miRNAs induced by drought stress. Up- and Down-regulated miRNAs (≥2 or ≤ −2 fold changes) under the control and drought stages of drought-susceptible and tolerant genotypes of sorghum are given in (A) with the confidence of P ≤ 0.05. The numbers of extended drought-responsive miRNAs with the confidence of ≥5 reads at least in one genotype are illustrated in (B).

TABLE 1

Table 1. Drought stress-responsive miRNAs (fold change >2 “at-lest” in one genotype) and their putative targets are listed.

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 stress-related 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.

FIGURE 3

Figure 3. Gene ontology (GO) analysis of miRNAs target genes identified in sorghum genotypes drought-susceptible and drought-tolerant. (A) Blue bars indicate the enrichment of miRNA targets in GO terms. Green bars indicate the percentage of total annotated sorghum gene mapping to GO terms. (B) Further dissection of “response to stimulus” exposes various stress responsive target genes, including drought (7 genes) and heat (6 genes) highlighted with red and yellow color, respectively.

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 novel-sbi-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.

FIGURE 4

Figure 4. Negative correlation between drought-responsive miRNAs and their predicted target genes; where (A,B) represents miRNA up, target gene down, and vice versa, respectively. The miRNA expressions obtained by next-generation sequencing is represented by red color. The color blue represents the expression level of target genes obtained from transcriptome analyses of sorghum.

Identification of Sorghum Tass

Phase-initiators direct the cleavage of trans-acting siRNA (tasi-RNA or TAS; a newly identified class of 21 nt short siRNAs) and later initiate the production of tasi-RNA clusters. Initially, we identified 135 unique clusters comprising 155 unique tasi-RNA sequences. The predicted tasi-RNAs were found to be conserved in plants. Further investigation revealed that six phase-initiators (e.g., miR390, miR5567, miR5568, miR6220, miR6225, and miR6230) directing 31 unique TASs lead to 64 unique tasi-RNAs (Table S7 in Supplementary Material). The length of the phase-initiator, a key determinant for triggering tasi-RNA biogenesis, is usually 22 nt (except 21 nt long miR390) in Arabidopsis (Chen et al., 2010). A similar exception was observed for sorghum miR390 (21 nt long). The length of other phase-initiators varied from 21 to 24 nt (e.g., 21/24 nt, miR5567; 21 nt, miR5568; 23/24 nt, miR6220, and 23/24 nt, miR6225). The variable length of the phase-initiators suggests the cleavage of double-stranded RNAs (dsRNAs) by multiple Dicer-like (DCL) proteins, thereby generating siRNA classes with different sizes (Axtell, 2009). Biogenesis of tasi-RNAs is dependent on miRNA triggers and requires either dual miRNA target sites (known as “two-hit” model) or single-target site (known as “one-hit” model) in the non-coding RNA precursor (Axtell et al., 2006). While most of the target genes were cleaved by only one miRNA at a single recognition site, we identified two target sites for five TASs (e.g., sbiTAS_miR5567/6225a-d and sbiTAS_miR5568/6220a) (Figure S2 in Supplementary Material). The length of each TAS in this study is restricted to 251 bp and the corresponding P-values varying from 4.47E-08 to 8.20E-05. The targets for tasi-RNAs were predicted by degradome supported psRobot server with default parameters. Predicted target genes for unique 64 tasi-RNAs were found to be involved in plant development and stress responses (Table S8 in Supplementary Material).

Discussion

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).

FIGURE 5

Figure 5. Distribution of sRNA reads produced from control and drought stress libraries of sorghum; where, (A,B) represents reads distribution in control and stress library of drought-susceptible genotype, respectively. (C,D) Represent reads distribution in control and stress library of drought-tolerant genotype, respectively. The highest sRNA read abundance was found for 24 nt in all 4 libraries.

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).

TABLE 2

Table 2. Conserved^β miRNA families in monocot and dicot species.

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 (Zhou et al., 2010) 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 (Zhou et al., 2010; 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 plants (Table 1; 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 novel-sbi-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 (Yang et al., 2009; Golldack et al., 2011). Moreover, miR5385 was up-regulated by drought in the tolerant cultivar, but down-regulated 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 novel-sbi-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 heat-shock 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), novel-sbi-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 novel-sbi-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 miR390-derived 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.

FIGURE 6

Figure 6. Diagrammatic representation of miR390 directed TAS3a (A) and TAS3b (B) genes in sorghum. Predicted miR30-TAS binding sites are shown in red; all tasi-RNAs are shown in other than red color with an underline; tasi-RNAs have been named by using a standardized nomenclature in which the register is given a number D1, D2… starting at the miRNA target site [17]; the processing direction is noted by a 5′ or 3′ prefix; the orientation is indicated by adding the suffix [+] for the positive (original transcript) strand, or [−] for the negative (RDR generated and denoted by #) strand.

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 genotype-specific 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 drought-responsive 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.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We thank to Dr. Monika Dalal, Indian Institute of Millets Research, Hyderabad (presently working as a Senior Scientist in NRC on Plant Biotechnology, New Delhi), for providing germplasm of sorghum. SM gratefully acknowledges the Department of Science and Technology (DST) for DST-INSPIRE fellowship grant.

Supplementary Material

The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fpls.2015.00506

Abbreviations

AGO, ARGONAUTE; ARF, auxin response factor; DCL, dicer-like 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.

References

Allen, E., Xie, Z., Gustafson, A. M., and Carrington, J. C. (2005). microRNA-directed phasing during trans-acting siRNA biogenesis in plants. Cell 121, 207–221. doi: 10.1016/j.cell.2005.04.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Ambros, V., Bartel, B., Bartel, D. P., Burge, C. B., Carrington, J. C., Chen, X., et al. (2003). A uniform system for microRNA annotation. RNA 9, 277–279. doi: 10.1261/rna.2183803

PubMed Abstract | CrossRef Full Text | Google Scholar

Anders, S. (2010). Analysing RNA-Seq Data with the DESeq Package. Heidelberg: European Molecular Biology Laboratory.

Google Scholar

Audic, S., and Claverie, J. M. (1997). The significance of digital gene expression profiles. Genome Res. 7, 986–995.

PubMed Abstract | Google Scholar

Axtell, M. J. (2009). The small RNAs of Physcomitrella patens: Expression, function and evolution. Ann. Plant Rev. 36, 113–142. doi: 10.1002/9781444316070.ch5

CrossRef Full Text | Google Scholar

Axtell, M. J., Jan, C., Rajagopalan, R., and Bartel, D. P. (2006). A two-hit trigger for siRNA biogenesis in plants. Cell 127, 565–577. doi: 10.1016/j.cell.2006.09.032

PubMed Abstract | CrossRef Full Text | Google Scholar

Barrera-Figueroa, B. E., Gao, L., Diop, N. N., Wu, Z., Ehlers, J. D., Roberts, P. A., et al. (2011). Identification and comparative analysis of drought-associated microRNAs in two cowpea genotypes. BMC Plant Biol. 11:127. doi: 10.1186/1471-2229-11-127

PubMed Abstract | CrossRef Full Text | Google Scholar

Bartel, B., and Bartel, D. P. (2003). MicroRNAs: at the root of plant development. Plant Physiol. 132, 709–717. doi: 10.1104/pp.103.023630

CrossRef Full Text | Google Scholar

Bonnet, E., Wuyts, J., Rouzé, P., and Van de Peer, Y. (2004). Evidence that microRNA precursors, unlike other non-coding RNAs, have lower folding free energies than random sequences. Bioinformatics 20, 2911–2917. doi: 10.1093/bioinformatics/bth374

PubMed Abstract | CrossRef Full Text | Google Scholar

Buchanan, C. D., Lim, S., Salzman, R. A., Kagiampakis, I., Morishige, D. T., Weers, B. D., et al. (2005). Sorghum bicolor's transcriptome response to dehydration, high salinity and ABA. Plant Mol. Biol. 58, 699–720. doi: 10.1007/s11103-005-7876-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Calviño, M., Bruggmann, R., and Messing, J. (2011). Characterization of the small RNA component of the transcriptome from grain and sweet sorghum stems. BMC Genomics 12:356. doi: 10.1186/1471-2164-12-356

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, H. M., Chen, L. T., Patel, K., Li, Y. H., Baulcombe, D. C., and Wu, S. H. (2010). 22-Nucleotide RNAs trigger secondary siRNA biogenesis in plants. Proc. Natl. Acad. Sci. U.S.A. 107, 15269–15274. doi: 10.1073/pnas.1001738107

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, T. W., Gan, R. C., Wu, T. H., Huang, P. J., Lee, C. Y., Chen, Y. Y., et al. (2012). FastAnnotator-an efficient transcript annotation web tool. BMC Genomics 13:S9. doi: 10.1186/1471-2164-13-S7-S9

PubMed Abstract | CrossRef Full Text | Google Scholar

Chinnusamy, V., Ohta, M., Kanrar, S., Lee, B. H., Hong, X., Agarwal, M., et al. (2003). ICE1: a regulator of cold-induced transcriptome and freezing tolerance in Arabidopsis. Genes Dev. 17, 1043–1054. doi: 10.1101/gad.1077503

PubMed Abstract | CrossRef Full Text | Google Scholar

Dai, X., and Zhao, P. X. (2011). psRNATarget: a plant small RNA target analysis server. Nucleic Acids Res. 39, W155–W159. doi: 10.1093/nar/gkr319

PubMed Abstract | CrossRef Full Text | Google Scholar

Du, Z., Zhou, X., Ling, Y., Zhang, Z., and Su, Z. (2010). agriGO: a GO analysis toolkit for the agricultural community. Nucleic Acids Res. 38, W64–W70. doi: 10.1093/nar/gkq310

PubMed Abstract | CrossRef Full Text | Google Scholar

Eldem, V. Çelikkol Akçay, U., Ozhuner, E., Bakır, Y., Uranbey, S., and Unver, T. (2012). Genome-wide identification of miRNAs responsive to drought in peach (Prunus persica) by high-throughput deep sequencing. PLoS ONE 7:e50298. doi: 10.1371/journal.pone.0050298

PubMed Abstract | CrossRef Full Text | Google Scholar

Feng-Xi, Y., and Di-Qiu, Y. (2009). Overexpression of Arabidopsis MiR396 enhances drought tolerance in transgenic tobacco plants. Acta Botanica Yunnanica 31, 421–426. doi: 10.3724/SP.J.1143.2009.00421

CrossRef Full Text

Ferreira, T. H., Gentile, A., Vilela, R. D., Costa, G. G., Dias, L. I., Endres, L., et al. (2012). microRNAs associated with drought response in the bioenergy crop sugarcane (Saccharum spp.). PLoS ONE 7:e46703. doi: 10.1371/journal.pone.0046703

PubMed Abstract | CrossRef Full Text | Google Scholar

Golldack, D., Lüking, I., and Yang, O. (2011). Plant tolerance to drought and salinity: stress regulating transcription factors and their functional significance in the cellular transcriptional network. Plant Cell Rep. 30, 1383–1391. doi: 10.1007/s00299-011-1068-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Grabherr, M. G., Haas, B. J., Yassour, M., Levin, J. Z., Thompson, D. A., Amit, I., et al. (2011). Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat. Biotechnol. 29, 644–652. doi: 10.1038/nbt.1883

PubMed Abstract | CrossRef Full Text | Google Scholar

Guilfoyle, T. J., and Hagen, G. (2007). Auxin response factors. Curr. Opin. Plant Boil. 10, 453–460. doi: 10.1016/j.pbi.2007.08.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Guo, H. S., Xie, Q., Fei, J. F., and Chua, N. H. (2005). MicroRNA directs mRNA cleavage of the transcription factor NAC1 to downregulate auxin signals for Arabidopsis lateral root development. Plant Cell Online 17, 1376–1386. doi: 10.1105/tpc.105.030841

PubMed Abstract | CrossRef Full Text | Google Scholar

Gutierrez, L., Mongelard, G., Floková, K., Pacurar, D. I., Novák, O., Staswick, P., et al. (2012). Auxin controls Arabidopsis adventitious root initiation by regulating jasmonic acid homeostasis. Plant Cell Online 24, 2515–2527. doi: 10.1105/tpc.112.099119

PubMed Abstract | CrossRef Full Text | Google Scholar

Hsieh, L. C., Lin, S. I., Shih, A. C., Chen, J. W., Lin, W. Y., Tseng, C. Y., et al. (2009). Uncovering small RNA-mediated responses to phosphate deficiency in Arabidopsis by deep sequencing. Plant Physiol. 151, 2120–2132. doi: 10.1104/pp.109.147280

PubMed Abstract | CrossRef Full Text | Google Scholar

Jeong, D. H., Park, S., Zhai, J., Gurazada, S. G., De Paoli, E., Meyers, B. C., et al. (2011). Massive analysis of rice small RNAs: mechanistic implications of regulated microRNAs and variants for differential target RNA cleavage. Plant Cell Online 23, 4185–4207. doi: 10.1105/tpc.111.089045

PubMed Abstract | CrossRef Full Text | Google Scholar

Jia, Y., Gu, H., Wang, X., Chen, Q., Shi, S., Zhang, J., et al. (2012). Molecular cloning and characterization of an F-box family gene CarF-box1 from chickpea (Cicer arietinum L.). Mol. Boil. Rep. 39, 2337–2345. doi: 10.1007/s11033-011-0984-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Jiangfeng, D. U., Yongjun, W. U., Xiaofeng, F., Junxia, C., Liang, Z., and Shiheng, T. (2010). Prediction of sorghum miRNAs and their targets with computational methods. Chinese Sci. Bull. 55, 1263–1270. doi: 10.1007/s11434-010-0035-4

CrossRef Full Text | Google Scholar

Jogeswar, G., Pallela, R., Jakka, N. M., Reddy, P. S., Venkateswara, R. J., Sreenivasulu, N., et al. (2006). Antioxidative response in different sorghum species under short-term salinity stress. Acta Physiol. Plant. 28, 465–475. doi: 10.1007/BF02706630

CrossRef Full Text | Google Scholar

Kantar, M., Unver, T., and Budak, H. (2010). Regulation of barley miRNAs upon dehydration stress correlated with target gene expression. Funct. Integr. Genom. 10, 493–507. doi: 10.1007/s10142-010-0181-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Katiyar, A., Smita, S., Chinnusamy, V., Pandey, D. M., and Bansal, K. (2012). Identification of miRNAs in sorghum by using bioinformatics approach. Plant Signal. Behave. 7, 246–259. doi: 10.4161/psb.18914

PubMed Abstract | CrossRef Full Text | Google Scholar

Kikuchi, K., Ueguchi-Tanaka, M., Yoshida, K. T., Nagato, Y., Matsusoka, M., and Hirano, H. Y. (2000). Molecular analysis of the NAC gene family in rice. Mol. Gen. Genet. 262, 1047–1051. doi: 10.1007/PL00008647

PubMed Abstract | CrossRef Full Text | Google Scholar

Kozomara, A., and Griffiths-Jones, S. (2013). miRBase: annotating high confidence microRNAs using deep sequencing data. Nucleic Acids Res. 42, D68–D73. doi: 10.1093/nar/gkt1181

PubMed Abstract | CrossRef Full Text | Google Scholar

Kulcheski, F. R., de Oliveira, L. F., Molina, L. G., Almerão, M. P., Rodrigues, F. A., Marcolino, J., et al. (2011). Identification of novel soybean microRNAs involved in abiotic and biotic stresses. BMC Genomics 12:307. doi: 10.1186/1471-2164-12-307

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, B., and Dewey, C. N. (2011). RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics 12:323. doi: 10.1186/1471-2105-12-323

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, F., Guo, S., Zhao, Y., Chen, D., Chong, K., and Xu, Y. (2010a). Overexpression of a homopeptide repeat-containing bHLH protein gene (OrbHLH001) from Dongxiang Wild Rice confers freezing and salt tolerance in transgenic Arabidopsis. Plant Cell Rep. 29, 977–986. doi: 10.1007/s00299-010-0883-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, W. X., Oono, Y., Zhu, J., He, X. J., Wu, J. M., Iida, K., et al. (2008). The Arabidopsis NFYA5 transcription factor is regulated transcriptionally and posttranscriptionally to promote drought resistance. Plant Cell Online 20, 2238–2251. doi: 10.1105/tpc.108.059444

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, Y. F., Zheng, Y., Addo-Quaye, C., Zhang, L., Saini, A., Jagadeeswaran, G., et al. (2010b). Transcriptome−wide identification of microRNA targets in rice. Plant J. 62, 742–759. doi: 10.1111/j.1365-313X.2010.04187.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, D., Song, Y., Chen, Z., and Yu, D. (2009). Ectopic expression of miR396 suppresses GRF target gene expression and alters leaf growth in Arabidopsis. Physiol. Plant. 136, 223–236. doi: 10.1111/j.1399-3054.2009.01229.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Long, T. A., Tsukagoshi, H., Busch, W., Lahner, B., Salt, D. E., and Benfey, P. N. (2010). The bHLH transcription factor POPEYE regulates response to iron deficiency in Arabidopsis roots. Plant Cell Online 22, 2219–2236. doi: 10.1105/tpc.110.074096

PubMed Abstract | CrossRef Full Text | Google Scholar

Lu, S., Sun, Y. H., Amerson, H., and Chiang, V. L. (2007). MicroRNAs in loblolly pine (Pinus taeda L.) and their association with fusiform rust gall development. Plant J. 51, 1077–1098. doi: 10.1111/j.1365-313X.2007.03208.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Ludlow, M., and Muchow, R. (1990). A critical evaluation of traits for improving crop yields in water-limited environments. Adv. Agron. 43, 107–153. doi: 10.1016/S0065-2113(08)60477-0

CrossRef Full Text | Google Scholar

Luo, Q. J., Mittal, A., Jia, F., and Rock, C. D. (2012). An autoregulatory feedback loop involving PAP1 and TAS4 in response to sugars in Arabidopsis. Plant Mol. Biol. 80, 117–129. doi: 10.1007/s11103-011-9778-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Meyers, B. C., Axtell, M. J., Bartel, B., Bartel, D. P., Baulcombe, D., Bowman, J. L., et al. (2008). Criteria for annotation of plant MicroRNAs. Plant Cell 20, 3186–3190. doi: 10.1105/tpc.108.064311

PubMed Abstract | CrossRef Full Text | Google Scholar

Moldovan, D., Spriggs, A., Yang, J., Pogson, B. J., Dennis, E. S., and Wilson, I. W. (2009). Hypoxia-responsive microRNAs and trans-acting small interfering RNAs in Arabidopsis. J. Exp. Bot. 61, 165–177. doi: 10.1093/jxb/erp296

PubMed Abstract | CrossRef Full Text | Google Scholar

Mukri, G., Biradar, B. D., Uppar, D. S., and Salimath, P. M. (2010). Influence of different temperature regimes on seed setting behavior and productivity traits in rabi sorghum. Electr. J. Plant Breed. 1, 1042–1048.

Google Scholar

Ooka, H., Satoh, K., Doi, K., Nagata, T., Otomo, Y., Murakami, K., et al. (2003). Comprehensive analysis of NAC family genes in Oryza sativa and Arabidopsis thaliana. DNA Res. 10, 239–247. doi: 10.1093/dnares/10.6.239

PubMed Abstract | CrossRef Full Text | Google Scholar

Pantaleo, V., Szittya, G., Moxon, S., Miozzi, L., Moulton, V., Dalmay, T., et al. (2010). Identification of grapevine microRNAs and their targets using high−throughput sequencing and degradome analysis. Plant J. 62, 960–976. doi: 10.1111/j.0960-7412.2010.04208.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Patel, R. K., and Jain, M. (2012). NGS QC Toolkit: a toolkit for quality control of next generation sequencing data. PLoS ONE 7:e30619. doi: 10.1371/journal.pone.0030619

PubMed Abstract | CrossRef Full Text | Google Scholar

Paterson, A. H., Bowers, J. E., Bruggmann, R., Dubchak, I., Grimwood, J., Gundlach, H., et al. (2009). The Sorghum bicolor genome and the diversification of grasses. Nature 457, 551–556. doi: 10.1038/nature07723

PubMed Abstract | CrossRef Full Text | Google Scholar

Raman, S., Greb, T., Peaucelle, A., Blein, T., Laufs, P., and Theres, K. (2008). Interplay of miR164, CUP−SHAPED COTYLEDON genes and LATERAL SUPPRESSOR controls axillary meristem formation in Arabidopsis thaliana. Plant J. 55, 65–76. doi: 10.1111/j.1365-313X.2008.03483.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Remy, E., Cabrito, T. R., Baster, P., Batista, R. A., Teixeira, M. C., Friml, J., et al. (2013). A major facilitator superfamily transporter plays a dual role in polar auxin transport and drought stress tolerance in Arabidopsis. Plant Cell Online 25, 901–926. doi: 10.1105/tpc.113.110353

PubMed Abstract | CrossRef Full Text | Google Scholar

Rhoades, M. W., Reinhart, B. J., Lim, L. P., Burge, C. B., Bartel, B., and Bartel, D. P. (2002). Prediction of plant microRNA targets Cell 110, 513–520. doi: 10.1016/S0092-8674(02)00863-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Rostamza, M., Richards, R. A., and Watt, M. (2013). Response of millet and sorghum to a varying water supply around the primary and nodal roots. Ann. Bot. 112, 439–446. doi: 10.1093/aob/mct099

PubMed Abstract | CrossRef Full Text | Google Scholar

Song, Q. X., Liu, Y. F., Hu, X. Y., Zhang, W. K., Ma, B., Chen, S. Y., et al. (2011). Identification of miRNAs and their target genes in developing soybean seeds by deep sequencing. BMC Plant Biol. 11:5. doi: 10.1186/1471-2229-11-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Stocks, M. B., Moxon, S., Mapleson, D., Woolfenden, H. C., Mohorianu, I., Folkes, L., et al. (2012). The UEA sRNA workbench: a suite of tools for analysing and visualizing next generation sequencing microRNA and small RNA datasets. Bioinformatics 28, 2059–2061. doi: 10.1093/bioinformatics/bts311

PubMed Abstract | CrossRef Full Text | Google Scholar

Sun, S. J., Guo, S. Q., Yang, X., Bao, Y. M., Tang, H. J., Sun, H., et al. (2010). Functional analysis of a novel Cys2/His2-type zinc finger protein involved in salt tolerance in rice. J. Exp. Bot. 61, 2807–2818. doi: 10.1093/jxb/erq120

PubMed Abstract | CrossRef Full Text | Google Scholar

Sunkar, R., Zhou, X., Zheng, Y., Zhang, W., and Zhu, J. K. (2008). Identification of novel and candidate miRNAs in rice by high throughput sequencing. BMC Plant Biol. 8:25. doi: 10.1186/1471-2229-8-25

PubMed Abstract | CrossRef Full Text | Google Scholar

Turner, M., Nizampatnam, N. R., Baron, M., Coppin, S., Damodaran, S., Adhikari, S., et al. (2013). Ectopic expression of miR160 results in auxin hypersensitivity, cytokinin hyposensitivity, and inhibition of symbiotic nodule development in soybean. Plant Physiol. 162, 2042–2055. doi: 10.1104/pp.113.220699

PubMed Abstract | CrossRef Full Text | Google Scholar

Unver, T., and Budak, H. (2009). Conserved microRNAs and their targets in model grass species Brachypodium distachyon. Planta 230, 659–669. doi: 10.1007/s00425-009-0974-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, W., Vinocur, B., Shoseyov, O., and Altman, A. (2004). Role of plant heat-shock proteins and molecular chaperones in the abiotic stress response. Trends Plant Sci. 9, 244–252. doi: 10.1016/j.tplants.2004.03.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, H. J., Ma, Y. K., Chen, T., Wang, M., and Wang, X. J. (2012). PsRobot: a web-based plant small RNA meta-analysis toolbox. Nucleic Acids Res. 40, W22–W28. doi: 10.1093/nar/gks554

PubMed Abstract | CrossRef Full Text | Google Scholar

Xie, F., Frazier, T. P., and Zhang, B. (2011). Identification, characterization and expression analysis of MicroRNAs and their targets in the potato (Solanum tuberosum). Gene 473, 8–22. doi: 10.1016/j.gene.2010.09.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, O., Popova, O. V., Süthoff, U., Lüking, I., Dietz, K. J., and Golldack, D. (2009). The Arabidopsis basic leucine zipper transcription factor AtbZIP24 regulates complex transcriptional networks involved in abiotic stress resistance. Gene 436, 45–55. doi: 10.1016/j.gene.2009.02.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, B., Pan, X., Cannon, C. H., Cobb, G. P., and Anderson, T. A. (2006). Conservation and divergence of plant microRNA genes. Plant J. 46, 243–259. doi: 10.1111/j.1365-313X.2006.02697.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, L., Zheng, Y., Jagadeeswaran, G., Li, Y., Gowdu, K., and Sunkar, R. (2011a). Identification and temporal expression analysis of conserved and novel microRNAs in Sorghum. Genomics 98, 460–468. doi: 10.1016/j.ygeno.2011.08.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, W., Gao, S., Zhou, X., Xia, J., Chellappan, P., Zhou, X., et al. (2010). Multiple distinct small RNAs originate from the same microRNA precursors. Genome Biol. 11:R81. doi: 10.1186/gb-2010-11-8-r81

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, X., Zou, Z., Gong, P., Zhang, J., Ziaf, K., Li, H., et al. (2011b). Over-expression of microRNA169 confers enhanced drought tolerance to tomato. Biotechnol. Lett. 33, 403–409. doi: 10.1007/s10529-010-0436-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, L., Liu, Y., Liu, Z., Kong, D., Duan, M., and Luo, L. (2010). Genome-wide identification and analysis of drought-responsive microRNAs in Oryza sativa. J. Exp. Bot. 61, 4157–4168. doi: 10.1093/jxb/erq237

PubMed Abstract | CrossRef Full Text | Google Scholar