Differential expression of microRNAs and other small RNAs in barley between water and drought conditions.

Drought is a major constraint to crop production, and microRNAs (miRNAs) play an important role in plant drought tolerance. Analysis of miRNAs and other classes of small RNAs (sRNAs) in barley grown under water and drought conditions reveals that drought selectively regulates expression of miRNAs and other classes of sRNAs. Low-expressed miRNAs and all repeat-associated siRNAs (rasiRNAs) tended towards down-regulation, while tRNA-derived sRNAs (tsRNAs) had the tendency to be up-regulated, under drought. Antisense sRNAs (putative siRNAs) did not have such a tendency under drought. In drought-tolerant transgenic barley overexpressing DREB transcription factor, most of the low-expressed miRNAs were also down-regulated. In contrast, tsRNAs, rasiRNAs and other classes of sRNAs were not consistently expressed between the drought-treated and transgenic plants. The differential expression of miRNAs and siRNAs was further confirmed by Northern hybridization and quantitative real-time PCR (qRT-PCR). Targets of the drought-regulated miRNAs and siRNAs were predicted, identified by degradome libraries and confirmed by qRT-PCR. Their functions are diverse, but most are involved in transcriptional regulation. Our data provide insight into the expression profiles of miRNAs and other sRNAs, and their relationship under drought, thereby helping understand how miRNAs and sRNAs respond to drought stress in cereal crops.

Introduction

Drought is a major constraint to crop production and occurs every year in many places in the world. Due to global warming, improving crops tolerant to drought is now becoming a primary objective of plant breeding. Barley (Hordeum vulgare L.) is an important cereal crop, ranking fourth among all grains in terms of quantity produced and area cultivated. Compared to its close wheat relative, barley is more tolerant to drought. In addition, barley has a short growing season with a high degree of natural and easily inducible variation. Therefore, barley would be an excellent model plant for investigating the genetic basis of drought tolerance.

MicroRNAs (miRNAs) are single-stranded small RNAs (sRNAs) with 20–24 nucleotides (nt) in length and encode no protein. miRNAs are generated from hairpin precursors (pre-miRNAs) that are formed from miRNA primary transcripts (pri-miRNAs), which are transcribed from genomic DNA. Other classes of sRNAs, including natural antisense transcript-derived siRNAs (natsiRNAs), repeat-associated siRNAs (rasiRNAs), long siRNAs (lsiRNA), heterochromatin siRNAs, secondary siRNAs, tiny noncoding (nc) RNAs (tncRNAs), 21U-RNAs, scan RNAs (scnRNAs), promoter/termini-associated sRNAs (PASRs/TASRs), transcription initiation RNAs (tiRNAs), transcription start site-associated RNAs (TSSa RNAs), splice site RNAs (spliRNAs) and sRNAs derived from rRNAs, snoRNAs, tRNAs and chloroplasts (Hackenberg et al., 2012; references therein), are also about 20–24 nt in size, but generated from long linear double-stranded RNAs rather than hairpin structures. Some dsRNAs can be generated from single-stranded sense transcripts by RDR6, a member of the RNA-dependent RNA polymerase (RdRP) family. In addition, secondary siRNAs can be further subdivided into ‘phased’ or ‘trans-acting’ (tasiRNAs) and natsiRNAs, which can act in cis (cis-natsiRNAs) or in trans (trans-natsiRNAs) (Axtell, 2013). miRNAs are predominantly involved in targeted mRNA degradation or translational repression (Huntzinger and Izaurralde, 2011), while siRNAs primarily mediate transcriptional silencing of genome loci crucial in chromatin remodelling and the maintenance of genome integrity or heterochromatic state (Elbashir et al., 2001; Pontes et al., 2008). The other classes of sRNAs are functionally diverse or unknown. Nevertheless, all of these classes of sRNAs are important in the regulation of biological processes (Jones-Rhoades et al., 2006).

There is extensive evidence indicating that miRNAs are regulated by drought in various plant species, including rice (Zhao et al., 2007; Zhou et al., 2010), populous (Li et al., 2011; Lu et al., 2008; Shuai et al., 2013), Arabidopsis (Li et al., 2008, 2011; Liu et al., 2008; Sunkar and Zhu, 2004), wheat (Kantar et al., 2011; Yao et al., 2010), maize (Li et al., 2011; Xu et al., 2011), soybean (Kulcheski et al., 2011), Medicago truncatula (Trindade et al., 2010; Wang et al., 2011), Phaseolus vulgaris (Arenas-Huertero et al., 2009), cassava (Ballen-Taborda et al., 2013), cowpea (Barrera-Figueroa et al., 2011), Prunus persica (Eldem et al., 2012), potato (Zhang et al., 2014), switchgrass (Xie et al., 2014), cotton (Wang et al., 2013), sugarcane (Ferreira et al., 2012) and tobacco (Frazier et al., 2011). In return, miRNAs turn on many genes in response to drought stress. Recent studies showed that overexpression of drought-regulated miRNAs leads to transgenic plants tolerant to drought (Li et al., 2008; Zhang et al., 2011). miRNAs thus have great potential to be used as a tool for improving drought tolerance in barley and other cereal crops. However, so far, little information is available on drought-regulated miRNAs or other sRNAs in barley, especially identified through genome-wide high-throughput sequencing. In barley, only dehydration-regulated miRNAs were identified using a computer-based approach and publicly available barley ESTs (Kantar et al., 2010). On the other hand, genome-wide identification and analysis of drought-responsive miRNAs and other sRNAs are very helpful in identifying regulatory networks between miRNAs and other sRNAs in plants. In this work, we analysed expression profiles of both miRNAs and other sRNAs in barley under water and drought conditions. We found that a lot of miRNAs and other sRNAs were regulated by drought. While low-expressed miRNAs and rasiRNAs were down-regulated, tsRNAs were up-regulated under drought conditions in barley, suggesting that different classes of sRNAs respond differently to drought, some of which are likely to be transcriptionally regulated by each other. Our data provide a valuable resource of miRNAs and other sRNAs for potential use in improving drought tolerance in barley and other cereal crops.

Results

Deep sequencing of sRNAs from Golden Promise barley under water and drought conditions

To identify drought-regulated miRNAs and other sRNAs in barley which might serve as future targets for improving drought tolerance of cereal crops, we deeply sequenced sRNAs in Golden Promise (GP) barley grown under water and drought conditions. Under drought conditions, GP [designated as GP(−w)] showed severe wilt symptoms, while under water conditions GP [designated as GP(w)] revealed no symptoms (Figure1). A total of 7 113 852 sequencing reads from the GP(w) sample and a total of 7 185 002 sequencing reads from the GP(−w) sample were obtained (Table 1). After removal of the adaptor sequences, low quality and short reads <16 nt, 6 671 598 and 6 664 821 sequencing reads remained in the GP(w) and GP(−w) samples, respectively (Table 1). The total number of reads is highly similar in both samples, which allows a meaningful comparison (especially of the number of unique reads which strongly depends on the total number of reads). After collapsing the reads into unique reads, 1 207 726 and 773 114 sequences were obtained in the GP(w) and GP(−w) samples, respectively (Table 1).

Figure 1

Phenotype of GP barley under water and drought conditions. The plants were watered for 3 weeks after germination. After that one plant was stopped for watering for 5 days, at which time the photograph was taken.

Table 1

The sequencing yield and outcome of processing

sRNA library	Raw reads	After filtration	Unique reads (UR)	Chloroplast read count (UC)	Chloroplast UR	Nuclear UC	Nuclear UR
GP(w)	7 113 852	6 671 598	1 207 726	2 752 340	64 836	3 755 898	1 091 816
GP(−w)	7 185 002	6 664 821	773 114	3 596 877	55 474	2 853 614	667 832

Cross-comparisons showed that 191 277 sequences are shared among both samples, while 1 282 989 were specific to the GP(w) and 836 578 to the GP(−w) sample. Among the shared reads, 16 827 were differentially expressed between the two samples (|log2| > 1), of which 5809 were up-regulated (log2 ≥ +1), while 11 018 were down-regulated (log2 ≤ −1) under drought conditions.

Length distribution of sRNA reads from water and drought-treated barley

The analysis of the read lengths of 18–25 nt gives valuable hints on the different sRNA species present in the samples. 20-nt sRNAs are the most frequent in GP(w) and GP(−w), followed by 21 nt (Figure S1A). Compared to GP(w), more 20-nt sRNAs are present in GP(−w). 21-nt sRNAs are the other way around. At a unique read level, the 24-nt peak is dominant, followed by 21 nt in both samples (Figure S1B). 24-nt reads have been proposed to be often generated from transposable elements (TEs) (Hackenberg et al., 2012). The 21-nt peak contains less unique reads but reads with higher read count in GP(w) compared to GP(−w). Shorter reads than 21 nt are more frequent in GP(−w) than in GP(w) in both read count and unique read number.

Classification of sRNAs from water and drought-treated barley

The mapping statistics to the barley nuclear and chloroplast genomes (International Barley Genome Sequencing Consortium, 2012; Saski et al., 2007) showed several differences (Table 1). We observed both higher read count and more unique reads mapped to the nuclear genome compared to the chloroplast in GP(w). On the contrary, we observed higher number of unique reads from the chloroplast in GP(−w) (Table 1). sRNAs from the GP(w) and GP(−w) samples were further grouped by mapping the reads to specific databases such as the RFAM (Gardner et al., 2009, release 10.0), miRBase version 19 (Griffiths-Jones et al., 2008), the TIGR repeat database (Ouyang and Buell, 2004), the TREP repeat database (Wicker ), RepBase (Jurka ) and tRNA databases. We observed that 79% and 87% of the reads from the GP(w) and GP(−w) samples, respectively, could be assigned to any of the database reference sequences (Table S1). The remaining unmapped reads may be due to unavailable genome or EST sequences, or due to sequencing errors. In Table S1, the higher mean read count in GP(−w) compared to GP(w) indicates that a higher number of low copy reads must exist in GP(w). The influence of sequencing quality can be ruled out as mapped reads have been used for the calculation.

About 12% in GP(w) and 8.3% in GP(−w) were mapped to miRNAs, 0.73% in GP(w) and 0.47% in GP(−w) were mapped to repetitive elements, 5% in GP(w) and 6.9% in GP(−w) were mapped to tRNAs, and 0.83% in GP(w) and 0.53% in GP(−w) were mapped to ncRNAs from the RFAM database (Table S1). Approximate 9.5% of the total sRNAs in GP(w) and 6.7% in GP(−w) were mapped to un-annotated genomic regions (Table S1). The antisense strands of the above genomic features were also detected. 39.6% of the total sRNAs in GP(w) and 51.6% in GP(−w) were mapped to the antisense strands of barley genes (from the HVGI database) whose proportion is the highest among all the mapped sRNAs (Table S1). Likewise, sRNAs originating from the chloroplast genome were also mapped to the sense and antisense strands in both samples. Significantly, a chloroplast-derived sRNA in both samples, which mapped to a chloroplast-encoded trnH-GUG gene, is the most frequent, accounting for 28.9% of the total sRNAs in GP(w) and 36.7% of the total sRNAs in GP(−w). In general, miRNAs, sRNAs derived from mRNAs, ncRNAs (RFAM) and those mapped to the genome but not classified are more frequent in GP(w), while antisense sRNAs (designated all antisense sRNAs as putative siRNAs) and sRNAs derived from tRNAs (tsRNAs) tend to be more frequent in GP(−w). For an unbiased comparison, the sRNA read counts are normalized to reads per million (RPM), calculated by dividing the read count by the total reads in the dataset and then multiplying 106. This result is shown in Figure2.

Figure 2

RPM values of different classes of sRNAs. RPM-reads per million; miRNA homologs-miRNAs from non-barley species in the miRBase; HVGI-barley EST database; Unclassified (barley genome)-mapped genome sequences are not annotated.

Length distribution of sRNAs in water and drought-treated barley

Figure S2 shows length distribution of sRNAs (18–25 nt). The patterns are similar between water and drought conditions. Only a subset of sRNAs derived from barley genes (HVGI) are smearly distributed. This proportion of sRNAs may be generated by degradation pathways rather than from a regulated process. Reads mapped to barley genes and chloroplast-derived reads have a dominant size at 20 nt. In contrast, the dominant size of miRNAs is 21 nt followed by 20 nt, which is consistent with DCL cleavage products. However, the dominant 21 nt vary in amount between water and drought conditions. While 74.6% of all miRNA reads have a length of 21 nt under water conditions, under drought conditions this percentage drops to 61.2% (Figure S3). This bias is in fact due to miR156. Firstly, miR156 is stronger expressed under drought conditions. Secondly, the by far most frequent miR156-derived read is 20 nt long and not 21 nt as the annotated sequence in miRBase. Note that the canonical miRBase sequence is also present but with an over two magnitudes lower read count compared to the most frequent read. In contrast to the above sRNAs, ncRNAs (from the RFAM database), repeat-associated sRNAs (in both sense and antisense directions) and the ‘genome mapped’ sRNAs, all have a dominant size of 24 nt (Figure S2). Of all the sRNAs, miRNAs and the reads that map the antisense direction of transcripts (designated as siRNAs) have the sharpest distributions. In the case of miRNAs, 99.1% and 98.8% of the reads are between 20 and 21 nt long in the GP(w) and GP(−w) samples, respectively. This confirms that these reads were processed, but not degradation products.

Differential expression of miRNAs between water and drought conditions

A total of 33 barley miRNAs from GP(w) and 31 from GP(−w) were detected using miRBase (Table S2). All the barley miRNAs identified in GP(−w) exist in GP(w). The two miRNAs missed in GP(−w) are hvu-miR399 and hvu-miR6177, which are both of low copies in GP(w). Among the 31 common miRNAs, 13 were significantly down-regulated (log2 ≤ −1), while only one (hvu-miR5049b) was significantly up-regulated (log2 ≥ +1) under drought conditions (Table S2). However, most of those miRNAs have low copy numbers. When a threshold of 10 RPM is adopted, only five, which were all down-regulated miRNAs, remain. It is worth noting that hvu-miR5048a is listed in the table as the reads do not map to this miRNA when using the strict criterion for known miRNAs, but when allowing one mismatch, this miRNA can be detected.

After removal of the barley miRNA-mapping reads, the remaining reads were further mapped (without mismatch) to miRNAs from other species in miRBase. A total of 416 putative homologous miRNAs from each of the GP(w) and GP(−w) samples were identified. However, when a threshold of 10 RPM was applied, only 74 putative homologous miRNAs from each sample remain (Table S3), of which only one (gma-miR6300) was significantly up-regulated under drought conditions (Table S3). In contrast, under the same conditions, 20 putative homologous miRNAs belonging to 10 miRNA families were significantly down-regulated (log2 ≤ −1) (Table S3). This expression profile is similar to that of the previously known barley miRNAs under the same drought conditions.

Novel miRNA identification and their differential expression between water and drought conditions

To identify novel miRNAs, we combined previously described criteria (Hackenberg et al., 2011) and alignment patterns of pre-miRNAs to the barley sRNA datasets. We hereby identified three novel miRNAs, designated as hvu-miRX33, hvu-miRX34 and hvu-miRX35, respectively. However, their miRNA* was not found in the datasets. hvu-miRX33 and hvu-miRX34 are of 22 nt, while hvu-miRX35 is of 21 nt (Table 2). All of these three miRNAs have higher negative minimal free energy index (MFEI) than other different type of RNAs, which is >0.85 (Zhang et al., 2006). hvu-miRX34 and hvu-miRX35 start with U at their 5′ ends, which is consistent with the preference of most miRNAs described previously (Kuang et al., 2009; Zhang et al., 2009). Another miRNA, hvu-miRX33, starts with G at its 5′ end. hvu-miRX33 was only expressed under water conditions. The other two miRNAs were down-regulated under drought conditions. hvu-MIRX35 gene (Accession No. dbj|AK252755.1) was found to contain another miRNA discovered previously (Schreiber et al., 2011). Homologous sequences of hvu-miRX33 and hvu-miRX35 are also present in wheat, but hvu-miRX34 only exists in barley.

Table 2

Identified novel miRNAs in water and drought-treated barley

miRNA	miRNA sequence	Read count (RC) GP(w)	RCGP(−w)	Precursor sequence
hvu-miRX33	CGGUAGGGCUGUAUGAUGGCGA	159	0	UCGCCAUCAUACGCCCAACCGUGCAUUUGAUAUGCAUAUAUAUGCAUCACGAGCCACGGUAGGGCUGUAUGAUGGCGA
hvu-miRX34	UGAGAAGGUAGAUCAUAAUAGC	279	2	CUCAGAUGAGAAGGUAGAUCAUAAUAGCUAAAAGAUAGUAUUCUCCGCAUCUCAUCUAAG
hvu-miRX35	UAAUCUUCUGGAAAUAUGCUU	26	12	AUAUCUGUCGGAUCUCUCUUAAUUUUGUAAUCUUCUGGAAAUAUGCUUAGGUGUAAU

Bolded are mature miRNA sequences.

Taken together, 9% of the reads from water conditions and 8% of the reads from drought conditions are mapped to miRNAs. The lower relative frequency of miRNAs under drought conditions may indicate that drought suppresses miRNA expression in barley.

Differential expression of ncRNAs between water and drought conditions

We further mapped all reads to the RFAM database and identified 185 ncRNAs in each of the GP(−w) and GP(w) samples (Table S4). Of these ncRNAs, five were found to be significantly up-regulated (log2 ≥ +1), while 21 were significantly down-regulated (log2 ≤ −1), under drought. snoRNAs, which primarily guide chemical modifications of rRNAs, tRNAs and small nuclear RNAs (snRNA), were dominantly regulated by drought. When using a threshold of 10 RPM, only 12 drought down-regulated ncRNAs and one (derived from RNaseP_arch) drought up-regulated ncRNA remained. Figure S4 shows the absolute read counts, relative fractions and unique read fractions for the different analysed sRNA species as a function of the regulation class (up-, down- or non-regulated). Drastic differences can be seen between the different RNA species; while all tsRNAs seem to be regulated, the vast majority of chloroplast-derived reads are non-regulated.

Differential expression of tsRNAs between water and drought conditions

Mapping all reads to tRNA databases showed that all 61 codons of tRNAs derived sRNAs in each sample, of which sRNAs derived from 19 codons of tRNAs were up-regulated (log2 ≥ +1), while sRNAs derived from two codons of tRNAs were down-regulated (log2 ≤ −1) under drought conditions (Table S5). When a threshold of 10 RPM was applied, only up-regulated tsRNAs remained at an anticodon level under drought conditions (log2 ≥ +1). This expression pattern contrasts to that of miRNAs or ncRNAs. Whether there is a correlation among tsRNAs, ncRNAs and miRNAs under drought conditions is interesting.

Differential expression of repeat-associated sRNAs between water and drought conditions

Alignment of all reads to the TREP repeat database, which contains 1716 TEs, revealed that 811 TEs derived sRNAs in GP(−w), while 1028 TEs derived sRNAs in GP(w) (Table S6). Of these TE-derived sRNAs, 771 were common between the two samples, 257 were specific to GP(w), while 40 were specific to GP(−w). In addition, a high proportion of sRNAs were found to be derived from the antisense strand of TEs. For example, 421 and 652 TEs derived antisense sRNAs in GP(−w) and GP(w), respectively. These antisense reads are likely to be transcribed from the antisense strand of the genome and thus may not be siRNA/piwiRNA-like sRNAs that can regulate the repeat expression. However, we cannot exclude the possibility that some of the antisense reads might originate from a repeat that inserted in antisense into any other transcribed element. In this case, these antisense reads may be able to function in the suppression of repeat expression. Nevertheless, all of these hypotheses need to be confirmed experimentally.

Among the identified antisense TE-derived sRNAs, 393 were common between the two samples, and 259 and 28 were specific to GP(w) and GP(−w), respectively. Of the common TE-derived sRNAs, only 12 sense and 9 antisense TE-derived sRNAs were significantly up-regulated (log2 ≥ +1) under drought conditions. In contrast, 406 sense and 200 antisense TE-derived sRNAs were significantly down-regulated (log2 ≤ −1) under the same conditions. When a threshold of 10 RPM was adopted, only 31 TE-derived sRNAs were significantly down-regulated, and no TE-derived sRNA was significantly up-regulated under drought conditions.

Table 3 shows the expression values as a function of different TE types. In general, there are no TE-specific differences between GP(w) and GP(−w), but a general down-regulation of all TE types under drought conditions, which is different from the regulation of tsRNAs under the same conditions. The majority of the mapped sense TE-derived sRNAs belong to retrotransposons (Class I TEs), accounting for 76% of the total mapped TE-derived sRNAs in GP(−w) and 73.5% in GP(w). TE-derived sRNAs belonging to DNA transposons (Class II TEs) account for 19.9% and 21.9% in GP(−w) and GP(w), respectively. There are 4.1% of TE-derived sRNAs from GP(−w) and 4.7% from GP(w) not being classified. Antisense TE-derived sRNAs gave a similar result as above (Table 4). Similar results were also observed from alignment of the reads to other repeat databases such as Repbase, the TIGR rice repeat database and the TIGR barley repeat database (data not shown).

Table 3

Expression of TREP repetitive elements

Name	GP(w) RC	GP(w) RPM	GP(−w) RC	GP(−w) RPM	log2[GP(−w)/GP(w)]
Gypsy	4509	466.6	2574	258.6	−0.85
Harbinger	378	41.6	160	20.7	−1.01
Mutator	146	23.1	69	10.5	−1.13
Unknown	960	128.4	460	60.3	−1.09
Mariner	713	92.1	358	39.6	−1.22
CACTA	1724	242.4	893	116.9	−1.05
Copia	5152	510.2	2930	294.8	−0.79

The table depicts two different levels, more specific (top) and more general (bottom).

Table 4

Reads that map in antisense direction to TREP repeats

Name	GP(w) RC	GP(w) RPM	GP(−w) RC	GP(−w) RPM	log2[GP(−w)/GP(w)]
Gypsy	865	115.5	469	57.6	−1.00
Unknown	194	27.8	88	14.0	−0.99
CACTA	451	64.1	234	31.4	−1.03
Copia	2810	282.6	1658	171.0	−0.72

Differential expression of chloroplast-derived sRNAs between water and drought conditions

Putative chloroplast-derived reads are strongly present in both GP(−w) and GP(w) (Table 1). Further analysis revealed that these sRNAs were distributed across the whole chloroplast genome. The comparison of chloroplast-derived sRNAs from the GP(−w) and GP(w) samples showed that the majority of chloroplast-derived sRNAs were not regulated by drought. However, tsRNAs were found again to tend to be up-regulated by drought.

Differential expression of antisense sRNAs between water and drought conditions

Antisense sRNAs (putative siRNAs) were identified via two steps: firstly, all the reads were mapped to the reverse strand of barley genes available in the databases, and then, the ‘antisense’ reads were mapped to the forward strand of all other libraries (tRNA, genes, repeats, etc.) in order to detect the RNA species from which the siRNAs are generated. In this way, a total of 5747 siRNA were identified to be common between water and drought conditions, and 914 siRNAs were found to be differentially expressed (|log2| ≥ 1) (Table S7). Among them, 172 siRNAs were up-regulated (log2 ≥ +1), while 642 were down-regulated (log2 ≤ −1) under drought conditions (Table S7). The up-regulated siRNAs were attributed to antisense tsRNAs, which, like sense tsRNAs, were also up-regulated under drought conditions.

Co-regulated miRNAs and other sRNAs in drought-treated and DREB3 transcription factor overexpressing transgenic barley

DREB transcription factors are among the first transcription factors discovered to be regulated under drought conditions and to enhance drought tolerance (Morran et al., 2011). We previously identified that overexpression of a DREB3 transcription factor from wheat (TaDREB3) in barley strongly affects the expression of sRNAs (including miRNAs) in the plants (Hackenberg et al., 2012). To identify truly drought-responsive miRNAs and siRNAs in barley under drought stress, we compared three sRNA datasets from the transgenic barley, GP(−w) and GP(w) samples. All of these three sRNA datasets were obtained under uniform conditions, for example the plants were grown at the same time in the same size of pots and same soil under the same glasshouse, and the same concentration of sRNAs was used for sequencing, which was performed at the same time in the same flowcell and generated a similar number of raw sequencing reads. Thus, these three datasets are comparable. For simplicity, we define consistent drought regulation as up-regulation in both GP(−w) and GP(TaDREB3) or down-regulation in both GP(−w) and GP(TaDREB3) when compared to GP(w). Any other combinations like up-regulation in GP(TaDREB3) but down-regulation in GP(−w) are considered as inconsistent regulation. As expected, both sample-specific and shared sRNA species were present in the datasets. Sample-specific sRNA species are generally of low read count and can be of several origins including sequencing errors and sample fluctuations. More importantly, we found many regulated (|log2| ≥ 1) sRNA species, both consistently and inconsistently regulated. The consistently regulated sRNAs are shown in Table S8.

Figure S5 shows that the proportion of consistent and inconsistently regulated sRNAs fluctuates strongly as a function of RNA species. TE-derived sRNAs are 100% consistently regulated, which are all down-regulated compared to GP(w). All other sRNA species show both consistent and inconsistent regulation. Among them, tsRNAs are the most variable class of sRNAs, which have a proportion of 65% tsRNAs that are inconsistently regul-ated between the transgenic and GP(−w) samples. As tsRNAs tended to be up-regulated under drought condition, this result confirms that the expression of tsRNAs is drought-regulated because the transgenic barley was not treated by drought. On the other hand, it also confirms our previous result that TaDREB3 strongly affects the expression of sRNAs in plant. Whether the transcription factor overcomes drought stress by regulating some tsRNA expression is interesting. All the significantly co-regulated sRNAs among the transgenic GP(−w) and GP(w) samples were listed in Table 5. miRNAs account for the most, followed by tsRNAs and other ncRNAs. Most of the co-regulated miRNAs and ncRNAs are down-regulated, while all the co-regulated tsRNAs are up-regulated compared to their corresponding RNAs in GP(w) (Table 5).

Table 5

Significantly co-regulated elements from all classes of sRNAs

Name	GP(w) RC	GP(w) RPM	GP(−w) RC	GP(−w) RPM	GP(TaDREB3) RC	GP(TaDREB3) RPM	log2[GP(−w)/GP(w)]	log2[GP(TaDREB3)/GP(w)]
hvu microRNA
hvu-miR5049a	546	103.2	102	17.6	116	23.4	−2.55	−2.14
hvu-miR5052	11	2.1	2	0.3	3	0.6	−2.59	−1.78
hvu-miR5053	23	4.3	6	1.0	5	1.0	−2.07	−2.11
hvu-miR1130	32	6.0	9	1.6	8	1.6	−1.96	−1.91
hvu-miR6191	9	1.7	2	0.3	3	0.6	−2.30	−1.49
hvu-miR1120	371	70.1	167	28.8	86	17.3	−1.28	−2.02
hvu-miR5048a hvu-miR5048b	4534	857.1	1722	296.8	2103	424.0	−1.53	−1.02
hvu-miR444b	64	12.1	29	5.0	30	6.0	−1.28	−1.00
Homologous microRNAs
osa-miR827a	1618	305.9	736	126.9	165	33.3	−1.27	−3.20
bdi-miR5054	1325	250.5	226	39.0	610	123.0	−2.68	−1.03
ath-miR169b	1834	346.7	516	88.9	667	134.5	−1.96	−1.37
ptc-miR169s	1665	314.8	486	83.8	613	123.6	−1.91	−1.35
gma-miR169v	1642	310.4	486	83.8	613	123.6	−1.89	−1.33
rgl-miR5139	302	57.1	106	18.3	109	22.0	−1.64	−1.38
hvu-miR5048a hvu-miR5048b	235	44.4	75	12.9	92	18.5	−1.78	−1.26
gma-miR393h	9009	1703.1	5139	885.8	1654	333.4	−0.94	−2.35
osa-miR5072	243	45.9	114	19.7	72	14.5	−1.23	−1.66
bdi-miR5064	274	51.8	148	25.5	72	14.5	−1.02	−1.84
tRNA (anticodon)
ValAAC	37 385	7067.6	113 732	19604.7	78 418	15808.7	1.47	1.16
RFAM
C0719	495	93.6	182	31.4	138	27.8	−1.58	−1.75
MIR812	3445	651.3	1878	323.7	1046	210.9	−1.01	−1.63
RNaseP_bact_b	148	28.0	79	13.6	60	12.1	−1.04	−1.21
tRNA (chloroplast)
trnS-GCU	5078	960.0	21 892	3773.7	11 765	2371.8	1.97	1.30
trnS-UGA	2100	397.0	6176	1064.6	4822	972.1	1.42	1.29
trnfM-CAU	2633	497.8	6432	1108.7	6454	1301.1	1.16	1.39

Examination of co-regulated antisense sRNAs among the transgenic GP(−w) and GP(w) samples showed that only small portion of antisense sRNAs were up-regulated in both the transgenic and GP(−w) samples. Most antisense sRNAs were down-regulated in both samples (Table S7). Further analysis showed that the co-regulated antisense sRNAs were generated from all sources of elements including tRNAs, repeats, nuclear and chloroplast, but the majority were derived from nuclear-encoded genes (Table S7).

Experimental validation of drought-regulated miRNAs and siRNAs

To experimentally confirm whether the above co-regulated miRNAs are indeed regulated by drought, Northern hybridization was applied using total RNAs isolated from leaf and root tissues of GP treated by water and drought conditions, respectively. Of 14 selected co-regulated miRNAs (11 conserved and three novel miRNAs, Table S9), three (hvu-miR444b, hvu-miR5049a and hvu-miRX35) were not detected in any of the tissues or conditions (data not shown), which may be due to low expression levels and/or low sensitivity of Northern hybridization. Of the detected miRNAs, most did not give significant change in expression level between water and drought conditions, for example hvu-miR168-5p and hvu-miRX34 (Figure3). Only hvu-miR159b, hvu-miR166a, ath-miR172a, osa-miR393a and hvu-miR5048 were obviously differentially expressed between water and drought conditions. Both hvu-miR159b and hvu-miR166a were up-regulated in leaf tissues, but down-regulated in root tissues, under drought conditions. ath-miR172a and hvu-miR5048 were down-regulated in both tissues under drought conditions. osa-miR393a was scarce in root tissues under both water and drought conditions, but appeared to be up-regulated under drought conditions. hvu-miR156, hvu-miR5048 and hvu-miRX33 were generally expressed less in leaf tissues than in root tissues under both water and drought conditions. hvu-miR159b and hvu-miR166a were expressed higher in root tissues only under water conditions. Under drought conditions, these miRNAs were expressed slightly higher in leaf tissues than in root tissues. ath-miR172a and osa-miR393a were expressed higher under water conditions than under drought conditions in both leaf and root tissues. Intriguingly, hvu-miR5052 was expressed higher in leaf tissues under water conditions, but under drought conditions, this miRNA was reversely expressed. hvu-miRX33 was up-regulated, while ath-miR172, osa-miR393a, hvu-miR5052 and hvu-miR5048 were down-regulated, under drought conditions.

Figure 3

Northern analysis of miRNAs in leaf and root tissues of GP and WI4330 barley under water and drought conditions.

Expression of the selected drought-regulated miRNAs was also examined in another barley cultivar, WI4330, by Northern hybridization. WI4330 is a breeding line and not highly tolerant to drought according to our unpublished data. The result showed that most of the miRNAs were expressed similarly as in GP (Figure3), indicating that miRNA expression tends to be conserved among barley cultivars. Only ath-miR169b was expressed slightly differently between the two cultivars. This miRNA acted to be down-regulated in both leaf and root tissues under drought conditions in WI4330.

The expression levels of the drought-regulated miRNAs were further quantitated by qRT-PCR. Of the selected 14 miRNAs (Table S9), only hvu-miR5052 was not regulated in either leaf or root tissues (Figure S6A). All other tested miRNAs were regulated to some degree by drought. hvu-miR156 and hvu-miR159b are expressed at high levels and are up-regulated under drought conditions in both tissues. Moderately expressed hvu-miR166a was up-regulated in leaves, but down-regulated in roots under drought conditions. ath-miR169b, ath-miR172a, hvu-miR444b, hvu-miR5048, hvu-miR5049a, hvu-miRX33 and hvu-miRX34 had a similar expression pattern as hvu-miR166a. In contrast, hvu-miR168-5p was only drought up-regulated in leaves. osa-miR393a and hvu-miRX35 appeared to be only expressed in leaf tissues, but not in root tissues. However, while osa-miR393a was up-regulated by drought conditions, hvu-miRX35 was down-regulated under the same conditions. Overall, the qRT-PCR result is consistent with that from the Northern hybridization or from the deep sequencing data.

Co-regulated siRNAs were also examined by Northern hybridization. Twelve co-regulated siRNAs designated as hvu-siRNA1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11 and 12, respectively, were selected (Table S9). siRNA4, siRNA7 and siRNA8 were derived from tRNAs, and siRNA1, siRNA2, siRNA3, siRNA10 and siRNA12 were derived from barley genes, while the rest were derived from unknown elements. It is notable that most of these siRNAs especially siRNA2 exist in multiple copies (Table S8). Only siRNA11 exists in a single copy, and siRNA3 and siRNA8 exist in two copies (Table S8). Northern hybridization showed that all siRNAs except for siRNA4, which was not detectable, were detected in at least one tissue and/or condition. siRNA10 was up-regulated while siRNA8 and siRNA9 were down-regulated in leaf and root tissues under drought conditions (Figure4). siRNA6 and siRNA12 were unchanged in expression level in both tissues under drought conditions. Unlike siRNA10, siRNA1 and siRNA2 were only up-regulated in leaf tissues under drought conditions. In root tissues, these two siRNAs as well as siRNA7 were almost undetectable. siRNA11 was not detectable in leaf tissues. siRNA3 was down-regulated in leaf tissues, but up-regulated in root tissues, under drought conditions.

Figure 4

Northern analysis of siRNAs in leaf and root tissues of GP and WI4330 barley under water and drought conditions.

Northern analysis of WI4330 showed that the expression pattern of above siRNAs was the same as in GP under the same conditions (Figure4). This indicates that the expression of siRNAs is also conserved in barley. We tried to use qRT-PCR to quantitate the expression levels of siRNAs, but failed for all siRNAs except siRNA11 which was present only in root tissues, being up-regulated by drought (Figure S6B). The failure of qRT-PCR for these siRNAs may be due to that they exist as clusters in the cells (Table S8) and hence cannot be amplified into products with a uniform size.

Prediction and validation of targets of drought-regulated miRNAs and siRNAs

To analyse the molecular function of drought-regulated miRNAs, their targets were predicted using psRNAtarget (http://plantgrn.noble.org/psRNATarget/). For most drought-regulated miRNAs, target genes can be predicted, but only some target genes are functionally known (Table S10). The known functions are mostly related to stress tolerance and developments, which emphasizes the importance of drought-regulated miRNAs. Notably, a considerable number of targets are transcription factors. Some predicted targets also encode housekeeping enzymes such as carboxylase and ribonuclease. In addition, all miRNAs have more than one target. Intriguingly, some genes are targeted by more than one miRNA as previously reported for other miRNAs’ targets. Gene Ontology (GO) analysis showed that three GO terms (ADP binding, defence response and response to stress) were significantly enriched (Table S10), and likely, these genes play important roles in the plant for response of drought stress.

To verify these predicted targets, a degradome library was constructed using leaves and roots of GP plants treated by water and drought. The advantage of this degradome library is that it requires no priori miRNA target prediction and can test all targets at once. Searching this library for miRNA targets showed that only nine predicted target genes could be verified (yellow highlight in Table S10). This small number may result from the limited barley EST sequences that were used to predict or map targets of the miRNAs. In addition, it could also be possible that the miRNAs whose cleavage products were not found in the library may function in certain particular tissues or developmental stage or via the transcriptional or translational repression mechanisms according to previous studies (Shamimuzzaman and Vodkin, 2012). Furthermore, it cannot be ruled out that most of the predicted targets may not be true. Some degradome reads were aligned to noncanonical (10–11) positions, suggesting that miRNA-mediated cleavage sites may not be restricted to positions 10–11.

To further confirm the predicted targets verified by the degradome library, qRT-PCR was applied on the basis of inverse relationship between miRNA and target, that is increasing miRNA level will reduce its target level or vice versa. However, this standard does not apply for the miRNAs that function in inhibiting transcription or translation of their targets. Eight genes targeted by five miRNAs were chosen (Table S10). We did not choose miR156's target like Squamosa promoter-binding-like transcription factor, because it has been well defined previously. Among the chosen targets, two are functionally unknown, while the others encode various functions (Table S10). The qRT-PCR result showed that most of these targets have a reverse expression level to their miRNAs in both leaf and root tissues (compare Figures S6A and S7). However, osa-miR393a shows a non-inverse relationship to one of its target genes, suggesting that this target may not be a bona fide target or it is regulated via a noncleavage mechanism.

Targets of drought-regulated siRNAs were also predicted using psRNAtarget. Unlike the miRNAs above, most siRNAs found no predicted targets (Table S10). For the predicted targets, most of them are not characterized. For the characterized predicted targets, their functions are quite diverse. One target was found to contain a BTB/POZ domain and is a regulator of hormone-responsive gene expression in barley (Woodger et al., 2004), suggesting that siRNAs may be involved in signalling pathways. In view of the fact that siRNAs function in the same way as miRNAs, we used the existing barley degradome library to verify the predicted targets of siRNAs (Table S11). The results showed that many targets had their cleavage products present in the library. Some targets even had their cleavage products occurring at different positions as did some of miRNA's targets. Unexpectedly, targets of the only single-copy siRNA11 in the cells were not found in the degradome library. We used qRT-PCR to examine two predicted targets (Table S10, Accession No: DN160496 and TC252151) of siRNA11, but failed as both targets were expressed at levels below the detection limits of qRT-PCR (data not shown). We did not pursue confirmation by quantifying the expression levels of predicted targets of the other siRNAs, because these siRNAs could not generate products with a uniform size by qRT-PCR.

Discussion

In this study, we used deep sequencing technology to analyse miRNAs and other sRNAs in barley under water and drought conditions. Seventy-seven per cent and 87% of sRNAs from the GP(w) and GP(−w) samples, respectively, were assigned to known entities from barley databases. The difference of the mapped sRNAs between the two samples is mainly attributed to tsRNAs, especially the most abundant trnH-GUG-derived sRNAs from the chloroplast genome, which were more prevalent in the GP(−w) sample than in the GP(w) sample. Another major contributor is siRNAs derived from antisense genes. Other sRNAs also contribute to the difference but not significantly. About 50% of rasiRNAs and those that mapped repeats [retro (transposons)] in an antisense orientation were differentially regulated between water and drought conditions. All repeat classes were down-regulated. Only a small fraction of barley-specific miRNAs are regulated by drought, most of which were down-regulated under drought stress. Overall, siRNAs and tsRNAs were up-regulated, and other sRNAs were down-regulated by drought (Figure5). The distributions of sRNAs within each sample are also different. While miRNAs account for 8.3%–12%, TE-derived sRNAs only account for 0.47%–0.73%, of the total sRNAs. This phenomenon has been observed in Arabidopsis where miRNAs and TE-derived sRNAs account for 14.4% and 16.9% of the total sRNAs, respectively, while intergenic region-derived sRNAs account for 49.1% of the total sRNAs (Lu et al., 2012). Intergenic region-derived sRNAs in barley only account for <9.5%. These results indicate that different sRNAs respond differently to drought stress, and sRNA generation is condition dependent. However, we found that the length distribution of sRNAs was not significantly affected by drought, despite that some proportion of sRNAs may be generated by degradation pathways. This suggests that sRNA length may be controlled independently of drought stress.

Figure 5

Distribution of overall sRNAs (siRNAs inclusive) of up-, down- and nonregulation under drought stress.

All tsRNAs were up-regulated under drought stress. In animals, tsRNAs also increase under stresses (Saikia et al., 2012), suggesting that a common mechanism might exist between plant and animal for tsRNA generation under stresses. Curiously, tsRNAs derived from the chloroplast were also up-regulated under drought stress. This suggests that (i) part of the chloroplast tsRNAs might be from the nucleus and/or (ii) tsRNAs in barley are genome-wide regulated (regulation is not limited to some particular organelles). The latter raises the interesting question on how two different sources of tsRNAs are coordinately regulated in plants. Perhaps, the processing of tsRNAs in the genome and chloroplast shares the same mechanism. In animals, tsRNAs have been found in the RISC and play a critical role in many biological processes (Lee et al., 2009). We believe that drought-regulated tsRNAs would have roles in drought tolerance. However, the molecular mechanism of tsRNA regulation by drought desires and merits further investigation.

Intriguingly, tsRNAs and repeat classes are drought-regulated oppositely. In addition, some cis-natsiRNAs were found to overlap with some trans-natsiRNAs, thereby forming double-stranded RNAs against each function. Furthermore, cleavage products of some siRNAs’ predicted targets were detected. These, combined with the fact that miRNAs can trigger siRNAs, lead us to speculate that regulatory networks may exist among sRNAs. The sRNAs may regulate or be regulated directly or indirectly via transcriptional silencing, post-transcriptional silencing and/or translational inhibition. It is likely that post-transcriptional silencing could be a dominant mechanism for all functional sRNAs because of the detection of many siRNAs’ cleavage products. In contrast, transcriptional silencing could mainly be conducted by rasiRNAs and miRNAs through methylation and chromatin remodelling to control the transcription of repeat sequences (i.e. transposons and retrotransposons) and genes that derive siRNAs, respectively. Translational inhibition would involve a third factor like transcription factors for controlling the transcription of genes that derive siRNAs. The role of such regulatory networks is to fine-tune the expression of miRNAs and siRNAs in order to coordinate with the expression of other genes in response to drought stress. However, we need to mention here that we could not find inverse expression of the natsiRNA pairs by analysing our sRNA datasets, suggesting that the functional mechanism of natsiRNAs remains elusive in barley.

All up-regulated miRNAs were of low copies, suggesting that these miRNAs may be subject to weak or no selection according to previous studies (Liang and Li, 2009). This would in turn mean that they tend to turn over quickly in evolution (Liang and Li, 2009). Nevertheless, low-expressed miRNAs are generally thought to have limited biological importance, and their effect on potential target mRNAs would, in most cases, be weak or negligible. However, we cannot rule out the possibility that these low-expressed miRNAs may participate in processes that can tolerate functional variation under drought stress. Apart from low expression, some miRNAs such as osa-miR393a were found to be significantly differentially expressed between leaf and root tissues. It is surprising that miR169 was not significantly regulated by drought in barley. This result is in agreement with that from Trindade et al. (2010) but not with that from Wang et al. (2011) in M. truncatula. Wang et al. (2011) showed that miR169 is down-regulated in M. truncatula under drought conditions. The discrepancy has been proposed to result from different drought treatments. In Arabidopsis, miR169 is also down-regulated under drought (Li et al., 2008), but in rice, this miRNA is up-regulated by drought (Zhao et al., 2007; Zhou et al., 2010). The up-regulation of miR169 in rice may be attributed by its upstream dehydration-responsive element (DRE), which is associated with abiotic stresses including drought stress (Zhao et al., 2007). The consequence of the unchanged or slightly up-regulated miR169 in barley under drought is yet unknown, but presumably, it would have little impact on plant response to drought stress. However, the ultimate proof of this hypothesis will come from its transgenic experiments overexpressing and knocking-down miR169.

Our previous study showed that overexpression of TaDREB3 in barley significantly affects the expression of various classes of sRNAs including miRNAs (Hackenberg et al., 2012). Comparing sRNAs from drought-treated barley and those from the transgenic barley reveals that many sRNAs including miRNAs were co-regulated. In other words, these sRNAs were consistently up- or down-regulated in both the drought-treated and transgenic plants. This provides an opportunity to identify true drought-responsive sRNAs, which might be used for improvement of plant drought tolerance in the future. However, further analysis finds that most of these co-regulated sRNAs are derived from protein-coding genes and hence unlikely to be generated via the miRNA/siRNA pathway. Therefore, these sRNAs may not play an important role in drought tolerance especially if those genes turn out to be not related to drought response. Why are some sRNAs co-regulated in plants? One possible explanation is that the expression of the protein-coding genes, from which the sRNAs are derived, may be controlled by other drought or TaDREB3-regulated genes. In Arabidopsis, the roles of different classes of sRNAs have been proposed: siRNAs maintain genome stability, whereas miRNAs mediate gene expression diversity (Ha et al., 2009). It is interesting to see if these properties are conserved in other species.

Experimental procedures

Plant materials, RNA isolation and sRNA sequencing

Golden Promise plants were grown in 6-inch pots in coco-peat soil in a glasshouse at 22–23 °C day/16 °C night with a 12-h day/night light cycle under water condition. Three weeks after germination, leaf and root tissues were harvested and further harvested after water was withheld for 5 days. Total RNA was extracted using TRIzol reagent and used for both Northern hybridization, reverse transcription PCR (RT-PCR) and quantitative real-time PCR (qRT-PCR). 18–30-nt sRNAs isolated from leaf tissues via 15% polyacrylamide gel were sequenced using the 36-base Illumina platform.

Processing of sRNA reads and alignment to reference sequences

The raw sequencing reads were processed by trimming reads at the first base that has a Phred score quality ≤2, eliminating the adapter sequence, removing reads with a strong compositional bias (frequency of the most base ≥0.9), and collapsing identical sequences into unique sequences, thereby assigning a read count/unique read (the number of times in a given RNA molecule was sequenced).

All reads were mapped to the nuclear and chloroplast genome sequences without mismatches and labelled as nuclear and chloroplast upon their inferred origin. Known miRNAs were detected by mapping all reads to miRBase without mismatches. Other classes of sRNAs were detected by mapping all reads to TIGR repeat database, Triticeae repeat sequence database (TREP) and all RepBase repeats, RFAM, HVGI and NCBI without mismatches. The reads mapped to a given library were not used for next library mapping. Antisense sRNAs (putative siRNAs) were detected by mapping all reads to the antisense strand of the elements.

Prediction of novel miRNAs and targets

Prediction of novel miRNAs was performed with miRDeep (Friedlander et al., 2008). Pre-miRNA stem–loop structures were examined with MFOLD (Zuker, 2003). Only pre-miRNAs meeting a minimum free energy lower than −16 kcal/mol were retained for next analysis. miRNA targets were predicted using psTarget (http://plantgrn.noble.org/psRNATarget). Gene annotation was obtained from the same source.

Construction of degradome library

Degradome library was constructed according to Addo-Quaye et al. (2008). Poly(A) RNAs isolated using the Oligotex Kit (Qiagen, Santa Clarita, CA) were ligated with a 5′ RNA adaptor containing a MmeI restriction site using T4 RNA ligase. Following reverse transcription, second-strand synthesis, MmeI digestion, ligation of a 3′ dsDNA adaptor and gel-purification, the cDNA products were amplified by PCR and sequenced with the Illumina HiSeq platform.

Differential expression analysis of sRNAs and their predicted targets by Northern hybridization and qRT-PCR

For Northern hybridization, 50 μg total RNA was electrophoresed on 15% polyacrylamide gel containing 7 m urea and transferred to Hybond-N membrane using 20 × SSC (3 m NaCl and 0.3 m C6H9Na3O9). The membrane was hybridized with 32P-labelled oligonucleotide probe reverse complementary to miRNA sequence. A U6 snRNA-specific probe served as a loading control. Prehybridization and hybridization were performed at 37 °C in 125 mm Na2HPO4 (pH 7.2), 250 mm NaCl, 7% SDS and 50% formamide. Washing was performed at 37 °C using 2 × SSC (twice) and 1 × SSC (once). The membrane was subjected to autoradiography at −80 °C for about 7 days.

For qRT-PCR, total RNA samples were polyadenlyated and reverse-transcribed using the NCode™ VILO™ miRNA cDNA Synthesis Kit (Invitrogen, Carlsbad, CA) containing poly(A) polymerase, ATP, SuperScript™III RT and a universal RT primer. cDNA was amplified with a miRNA- or target-specific forward primer and the RT primer using the following condition: 3 min at 95 °C followed by 45 cycles of 1 s at 95 °C, 1 s at 55 °C, 30 s at 72 °C (fluorescence reading acquired) and 15 s at 81 °C. Normalization was performed using three biological replicates and four control genes encoding glyceraldehyde 3-Pi dehydrogenase, heat-shock protein 70, cyclophilin and α-tubulin (Burton et al., 2004).