A Transgenic Transcription Factor (TaDREB3) in Barley Affects the Expression of MicroRNAs and Other Small Non-Coding RNAs

Transcription factors (TFs), microRNAs (miRNAs), small interfering RNAs (siRNAs) and other functional non-coding small RNAs (sRNAs) are important gene regulators. Comparison of sRNA expression profiles between transgenic barley over-expressing a drought tolerant TF (TaDREB3) and non-transgenic control barley revealed many group-specific sRNAs. In addition, 42% of the shared sRNAs were differentially expressed between the two groups (|log2| >1). Furthermore, TaDREB3-derived sRNAs were only detected in transgenic barley despite the existence of homologous genes in non-transgenic barley. These results demonstrate that the TF strongly affects the expression of sRNAs and siRNAs could in turn affect the TF stability. The TF also affects size distribution and abundance of sRNAs including miRNAs. About half of the sRNAs in each group were derived from chloroplast. A sRNA derived from tRNA-His(GUG) encoded by the chloroplast genome is the most abundant sRNA, accounting for 42.2% of the total sRNAs in transgenic barley and 28.9% in non-transgenic barley. This sRNA, which targets a gene (TC245676) involved in biological processes, was only present in barley leaves but not roots. 124 and 136 miRNAs were detected in transgenic and non-transgenic barley, respectively. miR156 was the most abundant miRNA and up-regulated in transgenic barley, while miR168 was the most abundant miRNA and up-regulated in non-transgenic barley. Eight out of 20 predicted novel miRNAs were differentially expressed between the two groups. All the predicted novel miRNA targets were validated using a degradome library. Our data provide an insight into the effect of TF on the expression of sRNAs in barley.


Introduction
Transcription factors (TFs) are important regulators of gene transcription in organisms and have been classified into a number of families based on their DNA-binding domain sequences and structures. Plant-specific APETALA2 (AP2)/ethylene-responsive element-binding proteins (EREBP) are one family of TFs, which in Arabidopsis contains 145 members and shares the AP2 DNAbinding domain of 70 amino acids originally identified in the floral homeotic protein APETALA2. The AP2/EREBP family is further divided into five subfamilies, AP2, RAV, DREB, ERF and others, based on the number of the AP2 domains presented in the genes. The DREB subfamily contains 56 members classified into six groups (A-1 to A-6) [1]. All of the members contain a single AP2 domain and bind to a consensus sequence (CCGAC) of the Crepeat or dehydration response element (DRE) in the promoters of genes that are turned on in response to low temperatures and/or water deficit. The DREB TFs were among the first TFs discovered to correlate with drought stress [2]. Over-expression of DREB TFs in transgenic plants remarkably increases the ability of plants to survive under strong drought and cold conditions [2]. The DREB TFs are also responsive to salinity stress [2]. All these functions are achieved by activating expression of their target genes, which include other TFs, phospholipase C, RNA-binding proteins, sugar transport protein, desaturase, carbohydrate metabolism related proteins, osmoprotectant biosynthesis proteins, proteinase inhibitors, late embryogenesis abundant (LEA) proteins (also known as dehydrins (DHNs)) and cold-responsive (COR) proteins [3][4]. DREB TFs can be induced by drought or by other DREB members [2].
MicroRNAs (miRNAs) are another class of important gene regulators, which can cleave a specific mRNA based on reverse sequence complementarities between the miRNA and target mRNA. MiRNAs regulate a wide range of biological processes including abiotic stresses. They are single-stranded non-coding (nc) small RNAs (sRNAs) between 20-24 nucleotides (nt) in length, which are processed from hairpin precursors (pre-miRNAs) by the Dicer-Like 1 enzyme (DCL1). Pre-miRNAs are derived from primary transcripts (pri-miRNAs) transcribed from genomic DNA, which is most likely to be regulated by TFs. On the other hand, TFs are possible miRNA targets. This reciprocal regulation has been shown in animals, which forms feedback loops [5]. In animals, miRNAs and TFs can regulate the same targets to form feed-forward loops [5]. However, in plants availability of data is limited.
Small interfering RNAs (siRNAs) are the second class of small regulatory RNAs. They are similar in size to miRNAs, but are processed from perfect long double-stranded RNAs derived from viruses, transgenes or endogenous transcripts such as transposons. Endogenous siRNAs are further classified into trans-acting siRNAs (tasiRNAs), natural antisense transcript-derived siRNAs (natsiR-NAs), repeat-associated siRNAs (rasiRNAs), long siRNAs (lsiRNA) and heterochromatin siRNAs. The majority of siRNAs are 24 nt and 21 nt in length. 24-nt siRNAs mainly guide silencing at the transcriptional level, while 21-nt sRNAs mainly guide silencing at the post-transcriptional level. siRNAs can direct DNA methylation, which is particularly important in the regulation of transposable elements (TEs) that account for a large portion of the genome size and are a major driver of evolution. natsiRNAs have been demonstrated to regulate salt tolerance [6], disease resistance [7] and a ubiquitous TF, nuclear factor Y (NF-Y) [8] that contributes to drought tolerance in Arabidopsis and maize (Zea mays L.) [9].
Piwi-interacting RNA (piRNA) is a distinct class of ncsRNAs, which is found to be expressed only in animal cells [10]. Unlike miRNA and siRNA, piRNAs size between 23-29 nt, lack of sequence conservation and are processed through Dicer-independent mechanisms [10]. piRNAs function mainly in silencing TEs, but its pathway is at the forefront of defence against transposons in cells, which is achieved through association with PIWI proteins. piRNA and PIWI proteins form an active piRNA-induced silencing complex (piRISC) used to recognize and silence complementary RNA targets [10].
Barley is one of the most important cereal crops worldwide, ranking fourth in terms of production. Recently, transformation with a DREB TF (TaDREB3) from wheat (Triticum aestivum L.) driven by a duplicated 35S promoter has been shown to confer strong tolerance to drought and cold in transgenic barley (H. vulgare cv. Golden Promise) [14]. Further study revealed that the over-expression of TaDREB3 in transgenic barley leads to elevated expression of other CBF/DREB genes and LEA/ COR/DHN genes known to be responsible for the protection of cells from damage and desiccation under stress conditions [14]. In this study, we performed deep sequencing of sRNAs from a transgenic barley line expressing TaDREB3 and from nontransgenic control barley. Comparative analysis showed that many group-specific and differentially expressed sRNAs existed between the two groups, including miRNAs and other classes of sRNAs. This indicates that the over-expressed TF affects the expression of various classes of sRNAs. Our data provides a valuable resource for identifying sRNAs related to the expression of DREB TFs in barley, whose outcome would help design a new effective strategy to cope with drought stress that severely limits the crop productivity.

Drought Tolerance of Transgenic Barley Over-expressing TaDREB3
T3 generation of transgenic barley lines with TaDREB3 was evaluated for drought tolerance. Non-transgenic barley plants were used as a control. Both transgenic and non-transgenic plants were deprived of water for a week. After this treatment, distinct differences were observed between the transgenic and nontransgenic plants. The transgenic plants looked vigorous, while the non-transgenic plants were severely wilted ( Figure 1). Analysis of water content in the plants showed that the transgenic plants contained more water than the non-transgenic plants. These results are consistent with previous studies [14][15]. The most drought-tolerant transgenic barley plants with TaDREB3 and the most drought-intolerant non-transgenic barley plants were selected for further analysis of expression profiles of miRNAs and other classes of sRNAs.

Sequencing of sRNAs from Non-transgenic and Transgenic Barley
Isolation and precise quantification of sRNAs are two keys for the comparison of different sRNA expression profiles. We used the Trizol reagent (Invitrogen) to extract total RNA from leaf tissues and the Purelink miRNA Isolation Kit (Invitrogen, Carlsbad, CA, USA) to isolate sRNAs to ensure good quality sRNAs for downstream gel purification. In addition, we adjusted the RNA preparations to the same concentration for ligation with the 59 and 39 adapters and for RT-PCR. Furthermore, we ran the RT-PCR products in the same flow cell for sequencing on the Illumina Genome Analyser so that we can assume the same sequencing error rate in both groups. These measures minimised artificial differences. In the end, we obtained 6,818,851 reads of 36 bases from transgenic barley and 7,107,714 reads of 36 bases from nontransgenic barley after filtering low quality reads. Both numbers are similar to each other. Reliable endogenous non-transgenic, constitutively expressed U6 spliceosome RNA [16] was perfectly parallel in read number between the two groups, indicating that a significant bias did not exist between groups. This thus allowed for the comparison of sRNA profiles between the two groups, and analysis of the effect of TaDREB3 on the expression of sRNAs.

Size Distribution of sRNAs in Non-transgenic and Transgenic Barley
The size distribution of the adapter trimmed reads from transgenic and non-transgenic barley is specified in Tables 1 and  2, and displayed in Figure 2A. The read size with the highest expression values (read count) was 20 nt followed by 21 nt in both groups, which is consistent with the Dicer cleavage products. However, the read count of 21 nt in transgenic barley was much lower than in non-transgenic barley. By contrast, the read count of 20 nt was much higher in transgenic barley compared to nontransgenic barley despite the fact that the total read count in nontransgenic barley was slightly higher than in transgenic barley. The third most expressed length was 22 nt in non-transgenic barley and 19 nt in transgenic barley. 22-nt reads in transgenic barley were ranked fourth, while 19-nt reads in non-transgenic barley ranked fifth, in read count. The fourth abundant size of reads in non-transgenic barley was 24 nt, which was at the fifth position in transgenic barley. The 19, 22 and 24 nt were very similar in read count to each other in non-transgenic barley. However, in transgenic barley these sizes were different in read count: 19-nt reads were more abundant than 22-nt reads, whereas 22-nt reads were more abundant than 24-nt reads.
After grouping the total reads, 445,672 and 989,288 unique reads remained in transgenic and non-transgenic barley, respectively (Tables 1 and 2). The unique read number was over twice in non-transgenic barley than in transgenic barley. Figure 2B shows size distribution of the unique reads. 24 nt was the most abundant read length followed by 21 nt in both groups. This, combined with the above size distribution of the read counts, indicated that 21-nt, especially 20-nt, reads comprised many highly expressed reads (high abundance), while 24-nt reads consisted primarily of unique or low abundant reads, which is in line with previous studies. This in turn demonstrated that the sequencing result from each group was reliable and each read could thus represent its relative abundance in vivo.
Comparative analysis showed that in non-transgenic barley the number of 24-nt unique reads was more than twice of 21-nt unique reads. In contrast, in transgenic barley the numbers of 24nt and 21-nt unique reads were similar to each other. Moreover, 24-nt unique reads were more than twice in non-transgenic barley than in transgenic barley. To check the abundance of unique reads as a function of its length, we applied a threshold to the read count of 4 (limit is included). Although the numbers of reads with sizes less than 20 nt or greater than 24 nt were not obviously reduced in the two groups, the numbers of 20-24-nt reads were significantly lower, especially in non-transgenic barley ( Figure 2B), thereby bringing the unique read numbers of the two groups back into approximation. This result demonstrated that more low abundant reads existed in non-transgenic barley than in transgenic barley (this is further confirmed when setting the threshold to 10). It is to note that we cannot exclude the possibility that a portion of the low abundant reads may be due to sequencing errors. However, the sequencing error should have the same magnitude in both samples as they have been run in the same flow cell. This means that we can assume that the existence of a lower number of low copy reads (specially 24 nt long reads) in transgenic barely is a direct consequence of the TF.

sRNAs Derived from Chloroplast and Nuclear in Transgenic and Non-transgenic Barley
Previous studies showed that a number of Illumina sequencing reads map to chloroplast genomes [13], [17][18][19]. These reads were classified into chloroplast-derived sRNAs (csRNAs), which have been proposed to play roles in abiotic stress tolerance [20]. To categorize chloroplast and nuclear sRNAs, we mapped all the trimmed reads to the complete barley chloroplast genome sequence [21] by means of the Bowtie aligner. Our analysis revealed that 3322007 reads represented   (Table 3). These reads were labelled as chloroplast derived reads. Correspondingly, the unmapped reads were labelled as nuclear derived reads or more precisely, non-chloroplast reads as the nuclear genome also contains the chloroplast-derived sequences [22]. It is expected that some chloroplast sRNAs (csRNAs) may be of nuclear origin because the sRNA processing machinery most likely resides in the nucleus.
The ratio of the chloroplast reads and non-chloroplast reads is 1.29:1 in transgenic barley and 0.72:1 in non-transgenic barley ( Table 3). The chloroplast reads were more in transgenic barley and less in non-transgenic barley than the non-chloroplast reads. However, the chloroplast unique read number in both groups was much smaller than the non-chloroplast unique read number ( Table 3), indicating that the chloroplast reads are highly redundant.
The non-chloroplast unique read and total read numbers in non-transgenic barley were 2.27 and 1.38 times more, respectively,    (Table 3). In order to compare the two groups, we multiplied the read counts of transgenic barley by a factor of 1.38, which generated identical total read counts in both groups. Cross comparison showed that in transgenic barley 328,794 non-chloroplast unique reads were specific, accounting for 77.9% of the total unique reads in the plant, while in nontransgenic barley 865,500 non-chloroplast unique reads were specific, accounting for 90.3% of the total unique reads in the plant (Table 3). However, when a threshold of 10 was set, the number of group specific reads was reduced by a factor of over 1000. In contrast, the number of shared reads (detected in both groups) was only reduced by a factor of 5. 36.2% of the shared reads showed a log2-ratio of |log2| . = 1 that corresponds to a differential expression between the two groups of at least a factor of 2, i.e. twice as high in one group compared to the other (Table 3). Interestingly, the fraction of differentially expressed reads increased when setting a threshold (Table 3). These data indicate that the over-expressed TF affects the expression of both nuclear and chloroplast sRNAs in the plant.

Expression Profile of Known miRNAs in Transgenic and Non-transgenic Barley
To detect known miRNA sequences with miRanalyzer [23], all the non-chloroplast reads shorter than 17 bases were removed from our datasets. The remaining reads were first mapped (without mismatch) to known barley miRNAs in the miRBase (miRBase version 17, April 2011). Nineteen out of the 22 barley miRNAs annotated in the miRBase were detected in nontransgenic barley and 17 detected in transgenic barley ( Table 4). The two extra miRNAs in non-transgenic barley were huv-miR444a and huv-miR1126 (Table 4).
hvu-miR156 was the most abundant miRNA in transgenic barley, accounting for 52.4% of the total reads, and the second most abundant in non-transgenic barley, while hvu-miR168-5p was the most abundant in non-transgenic barley, accounting for 71.3% of the total reads, and the second most abundant miRNA in transgenic barley (Table 4). Both miRNAs together accounted for 95.2% of the total miRNA reads in non-transgenic barley and 96.6% in transgenic barley (Table 4). When threshold . = 10 was set, 14 barley miRNAs in non-transgenic barley and 13 in transgenic barley retained. The extra miRNA in non-transgenic barley was huv-miR1126.
To detect previously un-annotated mature* sequences, which are determined by calculating the theoretical mature* sequence in the secondary structure assuming a perfect 39 overhang of 2 nt, all the reads were aligned to the pre-miRNA sequences. Three novel mature* sequences were detected (Table 4, highlighted in blue). It is worth mentioning that the annotated hvu-miR168-3p sequence in the miRBase might be wrong, because firstly the annotated sequence does not form an overhang with the mature (hvu-MIR168-5p) sequence in the secondary structure, and secondly the theoretical mature* sequence coincides with the most expressed read that map to the 3p arm in both groups. Apart from the detected mature* sequences, isomiRs of many barley miRNA were also detected in the two groups (Table 5). These isomiRs have their 59 or 39 ends longer or shorter than, or one nt different from, their reference miRNA sequences (data not shown).
Next, we removed the reads mapping to barley miRNAs and aligned the remaining reads (without mismatch) to known miRNAs from other species in the miRBase. A total of 126 miRNAs in transgenic barley and 150 in non-transgenic barley were mapped, of which 25 and 49 were specific to transgenic and non-transgenic barley, respectively, while 101 were common between the two groups (Table S1). Of these miRNAs, 23 are the origin of animals, of which 14 were specific to non-transgenic barley, 7 were specific to transgenic barley and 7 were common between the two groups. miRNAs conserved between plants and animals had been the case before [24] and could be used as a   phylogenetic marker to discover the evolutionary history and pattern between plants and animals. When threshold . = 10 was applied, 65 miRNAs in non-transgenic barley and 64 in transgenic barley remained (Table S1). vvi-miR172d (or ata-miR172) is only present in non-transgenic barley. With this threshold, only one putative animal miRNA (dme-miR1-3p) remained in the two groups, indicating that most of the animal miRNAs are low abundant. Interestingly, this putative animal miRNA was upregulated in transgenic barley and had its targets found in barley (Table 6). Whether it is indeed associated with the TF is not clear. All the detected non-barley miRNAs are considered as homologous miRNAs.
After the normalization of read count, 4 miRNA families were shown to be significantly up-regulated (log2. = 1) ( Table 6), while 24 miRNA families were significantly down-regulated (log2, = 21), in transgenic barley ( Table 7). Some of these differentially expressed miRNAs were further confirmed by Northern hybridization ( Figure 3A). Many differentially expressed miRNAs target transcription factors related to plant development, morphology and flowering time (Tables 6 and 7). Other differentially expressed miRNAs target the genes involved in plant metabolic metabolisms (succinyl-CoA ligase beta-chain, potassium transporter, NADH dehydrogenase), stress response (protein kinase, resistance protein) and post-translational modifications (ribosomal protein, F-box protein) (Tables 6 and 7). Remarkably, differentially expressed miR172 targets the AP2 subfamily, a close relative of the DREB subfamily. Both AP2 and DREB subfamilies belong to the AP2 TF family.

Predicted Novel miRNAs and their Target Genes in Transgenic and Non-transgenic Barley
To identify novel miRNA candidates, the remaining unmapped reads were aligned to barley genomic or EST sequences. If some reads did not map the available barley genomic or EST sequences, these reads would then be aligned to genomic or EST sequences from the closely related species, wheat, Brachypodium, or rice. After examining protein-coding possibility, repetitive sequences and hairpin structures of mapped pre-miRNA sequences using the criteria as described [13], we identified 20 putative novel miRNAs in each group, of which 7 (Hv-miRX3, 7, 9, 10, 16, 18 and 20) were down-regulated, while 1 (Hv-miRX4) were up-regulated, in transgenic barley (Table 8). Two pairs of miRNAs (Hv-miRX7/ Hv-miRX16, and Hv-miRX14/Hv-miRX15) were each derived from the same BAC clones. One miRNA (Hv-miRX19) localised between two repetitive sequences. Most predicted novel miRNAs were low abundant (Table 8). To experimentally confirm these putative miRNAs, northern hybridization was applied. All the selected miRNAs were detected in the leaf and root tissues ( Figure 3B). It appeared that some miRNAs were expressed differently between the tissues or between the groups ( Figure 3B). These results were further confirmed by RT-PCR ( Figure 3C) with primers specific to their precursor miRNAs (Table S2), indicating that the predicted novel miRNAs are bona fide miRNAs.
To reveal the function of the predicted novel miRNAs, their targets were predicted using psRNATarget (http://bioinfo3.noble. org/psRNATarget). As is common for barley, most predicted targets were functionally unknown (Table S3). Of the functionknown targets, many encoded protein kinases, suggesting that the corresponding miRNAs may be involved in signal transduction pathways especially in developmental processes. Only a few targets encode TFs. In contrast, the predicted targets of conserved miRNAs were enriched in TFs (data not shown). Most of the predicted novel miRNAs regulated the targets through cleavage (Table S3). However, 3 miRNAs regulated their targets through both cleavage and translational inhibition (Table S3). Moreover, although barley EST sequences are very limited, 18 out of the 20 miRNAs were found to target more than one gene (Table S3). In addition, one miRNA (Hv-miRX15) had two potential cleavage sites in one target (TC226135) (Table S3).
To validate the predicted miRNA targets, which could in turn confirm the predicted novel miRNAs, we constructed a degradome library based on the miRNA-directed cleavage between the 10th and 11th nt of complementarity relative to the guiding miRNA. If a target is cleaved by a miRNA, then its 39 fragment having a 59 monophosphate and a 39 polyA tail can be ligated with an RNA adapter containing a 59 MmeI restriction site. In contrast, full-length mRNAs cannot be ligated with this adapter as they are capped at the 59 ends. After reverse transcription, second-strand  TC207028 TC237973 TC195116 TC219375  TC218364 TC195122 TC195124 TC200044  TC225150 TC210261 TC223315 TC202631  TC222091 CA032492 BF258419 EX598089  TC209643 TC203114 TC232888 BF630636  TC219515 BQ766747 TC219782 BE421183  TC233676 TC203114 BQ767847 DN156768  TC198370 TC215757 TC198750  synthesis and MmeI digestion, a 20-21-nt sequence of the mRNA plus the 59 adapter remains, to which a 39 dsDNA adapter was added and the library was amplified and then sequenced. In general, these 20221nt signatures were sufficient to unambiguously identify the transcripts of origin and the corresponding miRNAs. As can be seen, the degradome library has advantages Figure 3. Analysis of miRNAs by Northern blot hybridization and RT-PCR. Total RNA was hybridized with probes for known miRNAs (A) and putatively novel miRNAs (B) and RT-PCR with primers specific for putatively novel miRNAs (C) that are shown in over the RNA ligase-mediated 59 rapid amplification of cDNA ends (RLM 59-RACE), which requires a priori miRNA target prediction and only tests a target at a time. We obtained a total of 5,270,019 sequences from the degradome library, of which 443,276 were unique. The sizes of sequences ranged from 13 to 26 nt, but the majority was between 20221 nt. Alignment of the unique sequences with the predicted novel miRNAs revealed that all the novel miRNAs had perfectly reverse complementary sequences in the library (Table 9). While most cleavage sites took place between positions 10 and 11 with respect to the guiding miRNA, other cleavage sites occurred at positions 11-15 (Table 9). In addition, some miRNAs such as Hvu-miRX2 and Hvu-miRX5 can cleave their targets at different positions (Table 9). These data indicate that miRNA-mediated cleavage sites are diverse. Moreover, a few miRNAs can each be aligned to more than one sequence tag (Table 9), indicating that these miRNAs may have multiple targets. Because of the limited barley ESTs, we only found ESTs available for 6 sequence tags (Table 9). However, 4 of these 6 sequence tags are contained by more than one EST (Table 9). With more ESTs available, it is expected that more ESTs can be found to contain these sequence tags.

Overall ncsRNAs in Transgenic and Non-transgenic Barley
The reads not being assigned to a miRNA were aligned with the Rfam database (allowing 2 mismatches), a comprehensive collection of ncRNA families represented by multiple sequence alignments. 664 ncRNAs were mapped in non-transgenic barley and 495 mapped in transgenic barley, of which 260 ncRNAs are specific to non-transgenic barley, while 89 are specific to transgenic barley (Table S4). Among these mapped ncRNAs, 4 are homologous to piRNAs (1 from transgenic barley and 3 from non-transgenic barley), but were all low abundant (Table S4). Considering that piRNAs have never been reported in plants, these piRNA homologs may not be bona fide piRNAs.

rasiRNAs in Transgenic and Non-transgenic Barley
To analyse rasiRNAs, all the reads were aligned with the Triticeae repeat sequence database (TREP), which contains 1,716 transposable elements (TEs). A total of 3823 reads from transgenic barley and 11978 reads from non-transgenic barley mapped to the TREP database. The mapped reads accounted for 0.37% of the total reads in transgenic barley and 0.8% in non-transgenic barley. These percentages are low, considering that about 80% of the barley genome is composed of repetitive sequences primarily consisting of TEs. Altogether, 433 TE-derived sRNAs in transgenic barley and 784 in non-transgenic barley were mapped, of which 132 were significantly down-regulated, while 2 were upregulated, in transgenic barley (|log 2 | .1) (Table S5). Notably, all the TE-derived sRNAs in transgenic barley are present in nontransgenic barley.
Among the mapped TE-derived sRNAs in transgenic barley, 186 belonged to retrotransposons (Class I TEs), while 169 to DNA transposons (Class II TEs), accounting for 43% and 39%, respectively (Table S6). The unclassified TE-derived sRNAs accounted for 18% of the total mapped TE-derived sRNAs (Table  S6). In non-transgenic barley, 387 belonged to retrotransposons, 330 to DNA transposons and 67 were unclassified, accounting for 49.4%, 42.1% and 8.5%, respectively (Table S6). The precent of each TE class in one group appeared to be similar in another group, suggesting that the regulation of TE-derived sRNAs in transgenic barley may not depend on the TE type, but on the expression levels. Similar results were observed when the reads were aligned to other repetitive databases such as Repbase (Table  S7), TIGR rice repeat database (Table S8), and TIGR barley repeat database (Table S9).

natsiRNAs in Transgenic and Non-transgenic Barley
natsiRNAs are another type of siRNAs generated from doublestranded RNAs. To identify natsiRNAs, all the reads were first mapped to the reverse strand of barley genes available in the databases. The ''anti-sense'' reads were then mapped to the forward strand of all other libraries (tRNA, genes, repeats, etc.). The results showed that 14,018 natsiRNA clusters existed in nontransgenic barley, while 8,297 were present in transgenic barley. Among them 3564 were common between the two groups, of which about half were significantly up-regulated in transgenic barley (Table S10). A lot of natsiRNA originated from TEs, suggesting that TEs can generate both sense and antisense sRNAs.
Some read clusters were generated from different transcripts and each targeted more than two genes. On the other hand, the same genes can be targeted by different read clusters. The redundancy of targeted genes accounted for 82.6% of the total mapped genes in non-transgenic barley and 84.3% in transgenic barley. Despite so, about 2/3 of the mapped genes were only targeted by single reads.
The most abundant natsiRNA in the two groups was derived from a ribosomal RNA, while the second most abundant natsiRNA in the two groups originated from P450 mRNA. Both natsiRNAs were more abundant in non-transgenic barley than in transgenic barley. Interestingly, both natsiRNAs had their sense counterparts in the groups. It is not clear which type of sRNA of these two natsiRNAs is really functional in vivo. The abundance of each of the other natsiRNAs was not consistent between the two groups, suggesting the expression of natsiRNAs was affected by the TF. tsRNAs in Non-transgenic and Transgenic Barley tsRNAs are a novel class of sRNAs presented in both plants and animals [12], [25]. 55 tsRNAs were identified in non-transgenic barley, while 54 were identified in transgenic barley (Table S11). The extra tsRNA in non-transgenic barley was derived from tRNA-LeuCAG. Of the identified tsRNAs in the two groups, tsRNA derived from tRNA-Gly(TCC) was the most abundant, accounting for 83.1% of the total tRNA-derived reads in nontransgenic barley and 57.9% in transgenic barley, followed by tRNA-Ala(AGC) derived sRNA, accounting for 8.4% in nontransgenic barley and 14.1% in transgenic barley (Table S11). Most identified tsRNAs were differentially expressed between the two groups, of which 6 tsRNAs were significantly up-regulated (log2. = 1), while 20 tsRNAs were significantly down-regulated (log2, = 21), in transgenic barley (Table S11). Distinctly, the most abundant tRNA-Gly(TCC)-derived sRNA in non-transgenic barley was 4 times higher in read count than in transgenic barley (Table S11). The varied expression levels may signify that the coordination may be adopted for tsRNAs in responsible for the over-expressed TaDREB3 in the cells.

sRNAs Derived from Transgenic TaDREB3 and its Downstream Regulated Genes
To determine if there were reads derived from TaDREB3, the remaining reads were aligned to the TaDREB3 sequence. 57 reads from transgenic barley were perfectly matched to TaDREB3 (Table S12). In contrast, none of the reads from non-transgenic barley were matched despite the existence of homologous genes in the plant (our unpublished data), thereby confirming that the matched reads were derived from the TaDREB3 transcripts. Most of the matched reads localised at the 39 region of the gene (Figure 4). Only a few reads, or no read, were mapped to other regions ( Figure 4). Intriguingly, a couple of reads from transgenic barley, but not from non-transgenic barley, were perfectly matched to the antisense strand of the TaDREB3 transcripts (Table S12) and the matched regions right overlapped with the mapped sense-strand regions. These co-incidences may indicate the existence of a feedback loop in transgenic barley, in which the transcribed TaDREB3 initiated the generation of sense sRNAs through normal degradation process, then the generated sense sRNAs function as siRNAs to cleave double-stranded structures of the TaDREB3 transcripts to produce anti-sense sRNAs, and the resulting antisense sRNAs cleave TaDREB3 to release sense sRNAs.
sRNAs derived from TaDREB3-regulated genes such as HvCBF15, HvCBF16, HvCBF10B, HvCBF23, HvCor14B, HvDHN5, HvDHN8, HvDHN9 and HvDHN10 were also detected in both groups. However, the mapped read number and read count were more in non-transgenic barley than in transgenic barley. In addition, the HvDHN8 gene was only mapped by the reads from non-transgenic barley, but not from transgenic barley. Curiously, antisense reads to these genes were also detected in the two groups, but the mapped read number and read count were again more in non-transgenic barley than in transgenic barley (Table S12). qRT-PCR showed that all of the above genes were expressed at higher levels in transgenic barley than in non-transgenic barley ( Figure 5), confirming that these genes were up-regulated by TaDREB3. HvCBF16 was expressed at the lowest level in both groups ( Figure 5). Correspondingly, the detected read number and read count reverse complementary to HvCBF16 were the highest among all the genes (Table S12), suggesting that the expression of HvCBF16 might be regulated through the miRNA/siRNA pathway.

Profile of sRNAs Derived from Chloroplast in Transgenic and Non-transgenic Barley
Next, we further analysed csRNAs in transgenic and nontransgenic barley. We found that the chloroplast reads were distributed across almost the whole chloroplast genome regions in both orientations (Table S13), suggesting that a normal degradation process is dominant. However, the reads mapped to the regions encoding tRNAs, rRNAs, PetB and PsbA appeared to be more than those mapped to intragenic, intergenic or other gene regions. In addition, the predominant sizes mapped to tRNAs and rRNAs in both groups were 20 nt and 19 nt, while there was no predominant size mapped to other RNAs or regions, suggesting that the degradation process is discriminating and that tsRNAs or rRNA-derived sRNAs (rsRNAs) might be processed in a special way. A few chloroplast reads from the two groups were mapped to several TEs (Table S14). Surprisingly, a tRNA-His(GTG) derived sRNA from the chloroplast genome, which matched to the 59 end of the tRNA-His (GTG) gene, was the most abundant of all the reads sequenced (Table S15). Because no other abundant reads were mapped to the other regions of the tRNA-His (GTG) gene (Table S15, Figure 6A), we speculate that this most abundant tsRNA may be generated under a special mechanism. Interestingly, this tsRNA was only present in the leaves, but not roots ( Figure 6B), which further supports its chloroplast origin because plastids in roots do not develop chloroplasts. Curiously, this tsRNA was found to perfectly match to the antisense strand of the 50S ribosomal gene (TC245676), which was also from the chloroplast genome. Whether the tsRNA and the gene interact in vivo is worthy of further investigation.

Discussion
Transcriptional regulation in eukaryotes is known to occur through the coordinated action of multiple TFs with other transcriptional regulators such as miRNAs and other classes of sRNAs. To gain insights into this cooperative network, we used next-generation sequencing technology to deeply sequence sRNAs in transgenic barley over-expressing the TaDREB3 TF and in nontransgenic control barley. Comparison of these sRNA expression profiles between the two groups revealed that the TF strongly affects the sRNA expression.
24 miRNAs were absent in transgenic barley (Table S1), suggesting that the expression of these miRNAs might be suppressed by TaDREB3. On the other hand, 20 miRNAs were transgenic barley specific (Table S1) implicating that these miRNAs might be activated by TaDREB3. Combined together, these results indicate that TaDREB3 affects the miRNA expression in multiways. It is to note that most of these groupspecific miRNAs are low copies because when threshold . = 10 is applied, only one miRNA is different between the two groups (Supplemental Table 1). This highlights the possibility that the low abundant group-specific miRNAs may be sequenced by chance in one sample but not in the other and thus not be truly group specific. Why some miRNAs were activated, while others were suppressed, by TaDREB3 is not clear, but likely to be associated with the function of TaDREB3. It is conjectured that most of these miRNAs may not be directly regulated by TaDREB3 because the binding site for TaDREB3 was not found in the upstream region of most miRNA genes that were regulated by TaDREB3. Instead, these miRNAs might be regulated through TaDREB3-related TFs or by means of regulating the miRNA host genes (in case of intronic miRNAs). In support of this, the predicted putative targets of miRNAs were not found to be the targets of TaDREB3 (Table S3).
Differentially expressed miRNAs were observed between the two groups. The most abundant miR156 was over two-fold upregulated in transgenic barley. Coincidently, miR156 is induced by drought [26] and TaDREB3 is specially functional in drought tolerance and induced by drought [14]. This suggests that the regulation of miR156 in vivo may be true and miR156 is a good candidate for drought tolerance and for further functional analysis of the impact of TaDREB3 on the expression of miRNAs. Previous studies showed that miR156 targets SQUAMOSA promoter binding protein-like (SPL) TFs, which control flowering time, phase change, leaf initiation rate and positively regulate miR172 [27][28]. Therefore, miR156 suppresses the expression of miR172. Correspondingly, miR172 was down-regulated in transgenic barley (Table S1). Intriguingly, miR172 targets the AP2 genes, which can also control the timing of flowering [29] and enhance drought tolerance [30]. These data indicate that differentially expressed miRNAs are connected to each other in complex networks in order to coordinate with the function of TaDREB3. It is to note that miR172 can in turn positively regulate the SPL TFs, but not miR156 [28].
Analysis of the putative miRNA sequences showed that uracil and adenine were the two predominant nt at the 59 ends, accounting for 45% each, while adenine was a predominant nt at the 39 ends, accounting for 55% ( Table 8). The preference of these nt at the 59 and 39 ends is in agreement with previous findings in rice and Arabidopsis [31]. However, in the miRBase, the most dominant uracil at the 59 ends accounts for 76% in Arabidopsis  and 58% in rice, followed by adenine, accounting for 12% in Arabidopsis and 11% in rice. The 59 terminal nt of miRNAs has been implied to be an important recognition signal for the argonaute (AGO) proteins [32], which are central to the function of an active RNA-induced silencing complex (RISC) containing Dicer and associated proteins. The 59 terminal uracil is preferentially recognised by AGO1, while the 59 terminal adenine is preferentially recognised by AGO2 [32]. The size of miRNAs also has an effect on the selection of AGOs. 24-nt miRNAs are preferentially recognised by AGO4, while the other sizes of miRNAs as well as the 24-nt miRNAs with C at the 59 end are preferentially recognised by AGO5 [32].
Alignment of sRNA sequences with the TaDREB3 gene revealed that some sRNA sequences were derived from TaDREB3. We observed that (i) these sRNA sequences were only detected in transgenic barley, but not in non-transgenic barley containing the TaDREB3-homologous gene; (ii) the TaDREB3-derived sRNA sequences were concentrated on some regions rather across the whole gene; and (iii) sRNA sequences reversely complementary to the TaDREB3 gene were only detected in transgenic barley and some overlapped forward TaDREB3-derived sRNA sequences. These observations suggest that TaDREB3 may be degraded at least partially via the miRNA/siRNA pathway. If this hypothesis is true, it could explain why so many sRNA sequences are different between the two groups and would be the first direct evidence for the feedback loop model between TFs and sRNAs in plants: TaDREB3 activates the expression of sRNAs and the activated sRNAs in turn regulate the expression of TaDREB3. Nevertheless, this model needs to be confirmed experimentally.
A number of sRNA sequences were mapped to the TaDREB3regulated genes, but differ in number and read count between the two groups. In addition, only the reads from non-transgenic barley, but not transgenic barley, were mapped to the TaDREB3regulated HvDHN8 gene (Table 12). On the other hand, both sense and antisense sRNAs derived from the TaDREB3-regulated genes existed (Table S12), some of which were mapped to the same regions of these genes (data not shown). These lead us to speculate that TaDREB3 and sRNAs may both be involved in the regulation of the TaDREB3-regulated genes and that the processing of these genes into sRNA fragments may be subject to the miRNA/siRNA pathway. However, following a normal degradation pathway is also possible because the reads were distributed across the whole regions of these genes. In addition, it cannot be ruled out that other TFs may also participate in the regulation of the DREB3-regulated genes. Taken together, the regulation of these DREB3-regulated genes may be a complex process.
To our surprise, a tRNA-His(GTG)-derived sRNA accounts for approximately 1/3 of the total sRNAs sequenced, thereby making it the most abundant sRNA and making tsRNAs the most abundant sRNAs in the two groups. This phenomenon had never been reported before in any other organism. In wheat, the same tsRNA only exists in a small number (our unpublished data) despite that wheat and barley are genetically close to each other. tRNAs not only have a role in translation, but also function in other cellular processes such as reverse transcription, porphyrin biosynthesis and others [33]. We expect that tsRNAs, especially the most abundant tsRNA, could have the related functions as above. This might be true because a recent study showed that the silencing of a tsRNA (tRF-1001) in human prostate cancer cells through a siRNA strategy decreases the cell viability [11]. On the other hand, it is also possible that tsRNAs, especially the most abundant tsRNA, may play a different role in barley from their counterparts in other species. Consistent with this possibility, we found that tsRNAs including the tRNA-His(GTG)-derived sRNAs were not randomly generated in barley (Table S14). In addition, we found that the most abundant tsRNA mapped to the reverse complementary strand of a gene (TC245676) involved in sRNA binding (GO:0003723), biological processes (GO:0006412), ribosome (GO:0005840 and GO:0003735) and plastid (GO:0009536). Due to the lack of bioinformatics, we did not pursue the identification of tsRNA targets.
A large number of sRNAs were mapped to the chloroplast genome (or chloroplast-derived sequences inserted at the nuclear genome) (Table S16). Similar situation has been described in rice [17] and Arabidopsis [34][35]. This indicates that csRNAs may have conserved functions among different species. However, some csRNAs might be by-products of chloroplast RNA processing, which may not be functionally relevant. Intriguingly, many barley csRNAs were mapped to intergenic regions between the psbB, psbT, psbH, petB and petD genes (our unpublished data). In particular, one of these regions was previously identified to contain a sRNA binding site for the cleavage of the polycistron [36]. This suggests that at least part of these csRNAs may result from the miRNA/siRNA pathway. This would then reflect a more ancient application of sRNA-directed transcript processing. Exploration of the enzymatic machinery would help determine sRNA processing in the chloroplasts.
Previous studies showed that TFs activate TEs. However, we found that the TE-derived sRNAs in transgenic barley were less in both number and read count than those in non-transgenic barley. This indicates that TaDREB3 may not be responsible for the activation of TEs. On the other hand, the detected transgenic barley-specific TE-derived sRNAs suggest that TaDREB3 is responsible for the activation of some TEs. Taken together, these indicate that TaDREB3 may have a dual role in the regulation of the TE-derived sRNAs. It is likely that TaDREB3 may regulate the TE-derived sRNAs through other TFs or siRNAs, the latter of which specifically function in silencing TEs. If it is true, because part of siRNAs may be generated from endogenous doublestranded RNAs (dsRNAs), which are processed from the inverted repeats, then a possible feedback loop may exist in transgenic barley, in which the TF transcribes TEs, and siRNAs generated from the transcribed TEs repress the TF, thereby reducing the number and abundance of TEs. It is to note that siRNAs can be generated from other classes of TEs and a complicated network may be involved in the generation. This would make it difficult for siRNAs to function consistently in silencing TEs. Any subtle change in the network could change sRNA expression patterns and functions significantly.
We have mentioned above that the miRNA/siRNA pathway may be adoptedd by different classes of sRNAs. This would raise the possibility that these sRNAs might compete with each other for Argonautes for their biogenesis and functions. In mouse and yeast loss of Dicer has been shown to increase one class of sRNAs, but decrease another [37][38], thereby confirming the existence of competition. However, such competition may be intricate and vary with the physiologic state of the cell. It may be inferred that each class of sRNAs has its own accessories, which allow avoiding using Argonautes for their biogenesis at the same time in the same cells. This means that while one class of sRNAs is generated, another class of sRNAs stands or decays. This strategy could be adopted by sRNAs within the same classes for their biogenesis. This could partly explain why different abundant sRNAs exist in the same cells. However, we cannot rule out the involvement of other unknown mechanisms in the generation of these sRNAs.
In summary, we present a comprehensive analysis of the sRNA components in barley including miRNAs, rasiRNAs, natsiRNAs, tsRNA, csRNAs and other non-coding sRNAs. Because the two groups only differ in TaDREB3, any difference in the expression of sRNAs between the two groups should be related to TaDREB3. Thus, our data increases our knowledge of the relationship of sRNA expression and TaDREB3. However, fundamental studies to assign the mechanisms of biogenesis and biological functions of the TF-related sRNAs are needed. brane (Amersham Bioscience) using 206SSC. The membrane was hybridized with 32P labelled DNA oligonucleotide probe reverse complementary to predicted miRNA sequence, made with c-32P-ATP using T4 polynucleotide kinase (New England Biolabs, Beverly, MA). U6 served as a loading control and its probe was made in the same way as above. Prehybridization and hybridization were both performed at 42uC in ULTRAhyb Ultrasensitive Hybridization Buffer (Ambion). Washing was also performed at 42uC using 26SSC, 0.1% SDS, twice followed by 16SSC, 0.1% SDS once. The membrane was either imaged using a Phosphorimager (Typhoon Trio, Amerisham Bioscience) or exposed to an X-ray film.

RT-PCR and Quantitative Real-time PCR
RT-PCR was performed as described [13] and quantitative real-time PCR (qRT-PCR) with primers was carried out according to Morran et al. [14]. Briefly, leaf tissues were collected from barley plants grown under well watered and drought conditions. Total RNA was extracted from the tissues using the Trizol reagent. Reactions were performed in an RG 6000 Rotor-Gene real-time thermal cycler (Corbett Research): 3 minutes at 95uC followed by 45 cycles of 1 second (s) at 95uC, 1 s at 55uC, 30 s at 72uC (fluorescence reading acquired), and 15 s at 81uC. The transcript levels of genes and control genes encoding glyceraldehyde 3-phosphate dehydrogenase, heat shock protein 70, cyclophilin, and a-tubulin were normalised as described [45].

Supporting Information
Table S1 Conserved miRNAs mapped by reads in nontransgenic barley (GP) and transgenic barley (GP(Ta-DREB3)). The log2-ratio indicates the differential expression. Sheet 1 shows conserved miRNAs mapped by the reads with threshold . = 10. Sheet 2 shows conserved miRNAs mapped by the reads with threshold . = 1. (XLSX)    Detected repeats in the RepBase database by reads (allowing 0 mismatch) in non-transgenic barley (GP) and transgenic barley (GP(TaDREB3)). The log2ratio indicates the differential expression. Sheet 1 shows repeats detected by the reads with threshold . = 10. Sheet 2 shows summary of the repeats mapped by the reads in non-transgenic barley (GP) and transgenic barley (GP(TaDREB3)). (XLSX) Table S8 Detected repeats in the TIGR rice database by reads (allowing 0 mismatch and threshold . = 10) in non-transgenic barley (GP) and transgenic barley (GP(TaDREB3)). The log2-ratio indicates the differential expression. (XLSX) Table S9 Detected repeats in the TIGR barley database by reads (allowing 0 mismatch and threshold . = 10) in non-transgenic barley (GP) and transgenic barley (GP(TaDREB3)). The log2-ratio indicates the differential expression.

(XLSX)
Table S12 Small RNAs derived from TaDREB3 and its regulated 10 genes in non-transgenic barley (GP) and transgenic barley (GP(TaDREB3)). The log2-ratio indicates the differential expression. (XLSX) Table S13 Chloroplast genome distribution of the reads (allowing 0 mismatch) in non-transgenic barley (GP) and transgenic barley (GP(TaDREB3)). Sheet 1 shows chloroplast genome distribution of the reads in GP. Sheet 2 shows chloroplast genome distribution of the reads in GP(TaDREB3).

(XLSX)
Table S14 TE-derived small RNAs detected in the chloroplast in non-transgenic barley (GP) and transgenic barley (GP(TaDREB3)). The log2-ratio indicates the differential expression. (XLSX) Table S15 Bowtie output of the read alignments with the tRNA-His gene in non-transgenic barley (GP) and transgenic barley (GP(TaDREB3)). Sheet 1 shows bowtie output of the read alignments with the tRNA-His gene in GP. Sheet 2 shows bowtie output of the read alignments with the tRNA-His gene in GP(TaDREB3).

(XLSX)
Table S16 Summary of all reads obtained, trimmed, unique or redundant, identified as various classes of small RNAs, unclassified or unassigned from the nuclear genome or from the chloroplast genome in non-transgenic barley (GP) and transgenic barley (GP(TaDREB3)). (XLSX)