Skip to main content

Parallel analysis of RNA ends enhances global investigation of microRNAs and target RNAs of Brachypodium distachyon

Abstract

Background

The wild grass Brachypodium distachyon has emerged as a model system for temperate grasses and biofuel plants. However, the global analysis of miRNAs, molecules known to be key for eukaryotic gene regulation, has been limited in B. distachyon to studies examining a few samples or that rely on computational predictions. Similarly an in-depth global analysis of miRNA-mediated target cleavage using parallel analysis of RNA ends (PARE) data is lacking in B. distachyon.

Results

B. distachyon small RNAs were cloned and deeply sequenced from 17 libraries that represent different tissues and stresses. Using a computational pipeline, we identified 116 miRNAs including not only conserved miRNAs that have not been reported in B. distachyon, but also non-conserved miRNAs that were not found in other plants. To investigate miRNA-mediated cleavage function, four PARE libraries were constructed from key tissues and sequenced to a total depth of approximately 70 million sequences. The roughly 5 million distinct genome-matched sequences that resulted represent an extensive dataset for analyzing small RNA-guided cleavage events. Analysis of the PARE and miRNA data provided experimental evidence for miRNA-mediated cleavage of 264 sites in predicted miRNA targets. In addition, PARE analysis revealed that differentially expressed miRNAs in the same family guide specific target RNA cleavage in a correspondingly tissue-preferential manner.

Conclusions

B. distachyon miRNAs and target RNAs were experimentally identified and analyzed. Knowledge gained from this study should provide insights into the roles of miRNAs and the regulation of their targets in B. distachyon and related plants.

Background

The wild grass Brachypodium distachyon (Brachypodium) belongs to the Pooideae family of grasses and has emerged as a model system for temperate grasses including agronomically important crops such as wheat (Triticum monococcum, T. durum and T. aestivum), barley (Hordeum vulgare) and oat (Avena sativa), as well as for biofuel plants such as switchgrass (Panicum virgatum) and Miscanthus (Miscanthus x giganteus) [1, 2]. It has many attractive features that resemble those of other models like Arabidopsis (Arabidopsis thaliana) and rice, including a relatively small and sequenced diploid genome, a well-developed transformation system and a growing bank of available Transfer-DNA insertional mutants [3–6]. Furthermore, as a wild grass, Brachypodium has never been subjected to human selection, which has drained the genetic diversity of domesticated cereal grasses. For these reasons, Brachypodium is a promising model system among temperate grasses for the identification of traits that may have been lost during domestication. In addition, phylogenetic surveys indicate that some important cereal crops, such as wheat and barley, are much more closely related to Brachypodium than they are to rice, sorghum or maize [7], and therefore studies involving Brachypodium are easily translatable and more applicable to these crops than studies of other model systems.

Brachypodium is also proposed as a model system for biofuel crops for understanding topics such as the characteristics of the grass cell wall involved in the process of saccharification [8]. However, despite its increasing importance as a model system, many areas of fundamental research in Brachypodium are underdeveloped, such as that related to microRNAs (miRNAs) and their targets.

miRNAs are critical regulators of the expression of many specific target genes involved in development, stress and other processes [9–11]. MIRNA genes are transcribed by RNA polymerase II, and processed by a Dicer-like (DCL) protein from stem-loop structures. Mature miRNAs are loaded onto Argonaute (AGO) proteins to form RNA-induced silencing complexes (RISCs) and perform target RNA regulation by base pairing [12, 13]. Although it is reported that some plant miRNAs inhibit target mRNA translation [14–17] or direct DNA methylation of target genes [18], the majority of plant miRNAs, once loaded, guide RISCs to cleave their target RNAs, which are subsequently degraded. In some cases, miRNA-mediated cleavage of a target RNA can trigger the production of secondary small interfering RNAs (siRNAs), typically comprising trans-acting siRNAs (tasiRNAs) and/or phased siRNAs (phasiRNAs) [19–21]. The secondary siRNAs generated from a limited number of primary targets can silence large gene families such as Pentatricopeptide Repeat (PPR) genes in Arabidopsis or Nucleotide Binding Site-Leucine-Rich Repeat (NBS-LRR) genes in other dicot plants [21–24]. Other than the tasiRNAs generated by the 21-nucleotide miR390 [25], most of the other tasiRNAs (or phasiRNAs) are generated by 22-nucleotide-long miRNAs [26–29]. In rice and Brachypodium, 22-nucleotide miR2118 and miR2275 can trigger 21-nucleotide and 24-nucleotide phasiRNA generation, respectively, from hundreds of non-coding loci of unknown function [30–32].

In diverse plants and samples, numerous miRNAs have been identified and profiled using deep sequencing and bioinformatic tools [33]. miRNA prediction from deep-sequencing data relies on the unique features of plant miRNAs among the other types of small RNAs [21, 34–36]. Currently, more than 5,100 miRNA precursors from 67 plant species have been annotated and deposited in miRBase [37]. In Brachypodium, 136 mature miRNAs from 135 miRNA precursors have been identified from different tissues and drought- or cold stress-treated tissues [38–42]. However, the number of annotated miRNAs for Brachypodium is relatively low compared to other model plants, such as the 338 miRNAs for Arabidopsis and 708 miRNAs for rice. Annotation of miRNAs from the other temperate grasses and biofuel plants is also very limited. Moreover, there have been rising concerns about incorrect plant miRNA annotation. These misannotations were mostly due to computational analysis without experimental validation. For example, by adopting the recently updated plant miRNA criteria [36], re-evaluation of about 400 miRBase-annotated rice miRNAs via a large-scale analysis revealed that about 150 had characteristics more typical of siRNAs [34]. Thus, for the continued development of Brachypodium as a model plant, it is a worthwhile challenge to identify with high confidence a more complete set of miRNAs.

Similarly, to understand the function of miRNAs, the identification of their mRNA targets, preferably with experimental support, is essential. Computational prediction of target RNAs relies on the feature of complementary base pairing between miRNA and target RNA pairs. Parallel analysis of RNA ends (PARE) is a technique for RNA degradome analysis that provides high-throughput experimental evidence of miRNA-mediated target cleavage [43–45]. Using PARE, target RNA cleavage products can be cloned and deeply sequenced. This has been a powerful approach for global target identification and validation in Arabidopsis, legumes, rice and a few other plants [21, 43–52]. The miRNA-mediated target cleavage is an endonucleolytic site-specific cleavage that cuts within the target RNA strand of the miRNA-target RNA duplex between nucleotides 10 and 11 relative to the 5′ end of the miRNA sequence. The single-nucleotide resolution of PARE data can therefore provide detailed information about site-specific cleavages guided by the individual miRNAs in the samples examined. Of temperate grasses and biofuel plants, only barley has been studied with a focus on global target identification using PARE data [52]. Although the generation and analysis of PARE data from Brachypodium was recently reported for the study of general mRNA decay [53], miRNA-mediated cleavage was only briefly presented.

In this study, we investigated Brachypodium miRNAs by deeply sequencing 17 small RNA libraries from various tissues and environmentally stress-treated plants. Using a computational pipeline, we identified not only conserved miRNAs that have not been reported in Brachypodium, but also non-conserved miRNAs that were not found in other plants. In addition, we identified a new regulation of miR5200 under submergence stress. Finally, we identified many miRNA target cleavages using four PARE libraries from key tissues. The knowledge gained from this study should be beneficial for understanding the role of miRNAs and their target regulation in Brachypodium and related temperate grasses and biofuel plants.

Results

Massive sequencing of Brachypodiumsmall RNAs

To enhance the knowledge and understanding of the small RNA population in Brachypodium, a very large set of small RNAs was cloned and deeply sequenced with Illumina technology. In total, 17 libraries were constructed and sequenced to achieve a broad coverage of small RNAs in various tissues and stresses (Table 1). The plant material is described in Materials and methods. Briefly, eight tissue libraries are from root, seedling, stem, leaf and stem tissues, and two biological replicates of leaf and panicle tissues. Nine stress libraries include drought, salt, cold, heat, submergence treatments and their control, phosphate-starvation and its control, and a pool of abiotic, biotic, and mechanical stresses. After trimming the adapter sequences, a total of 94,101,750 reads were obtained and aligned against the Brachypodium genome [31] using Bowtie. About 80%, or 75,048,757, matched perfectly to the Brachypodium genome. This resulted in an average abundance of approximately 4.4 million and approximately 1.2 million distinct genome-matched small RNA sequences per library. To our knowledge, the population sequenced in this paper is the largest such effort to date for Brachypodium and greatly enhances the small RNA resources for this plant.

Table 1 Summary statistics of small RNA libraries

Identification of BrachypodiummicroRNAs and microRNA precursors

Although many miRNAs have been identified from Brachypodium, the large set of small RNAs in this study provides a tremendous opportunity for identification of new miRNAs. In a recent release of miRBase (Release 19), a total of 135 miRNA precursors have been annotated for Brachypodium[37]. Since miRNA identification for Brachypodium is not saturated yet compared to other well-studied plants, such as Arabidopsis and rice, we applied a computational pipeline adapted from our recent reports [21, 34] to identify more miRNAs and their precursors. Out of approximately 12 million genome-matched distinct small RNA sequences, approximately 37,500 sequences passed the first filter, which consists of cut-offs for abundance, size and hits to the genome (Figure 1). These sequences were then evaluated for miRNA:miRNA-star (miRNA*) pairing, and 3,711 miRNA candidates were extracted from 5,627 putative miRNA precursor loci. In a third filter, strand bias and abundance bias cut-offs were used to distinguish the miRNA precursors from siRNA loci, resulting in 154 sequences from 199 loci. Finally, stem-loop structures of miRNA precursor candidates were visually inspected to ensure accuracy and 81 miRNAs generated from 121 precursors were identified (Additional file 1: Table S1). Excluding 52 miRNAs that were previously annotated for Brachypodium or similar to known plant miRNAs, 29 miRNAs from 30 precursors were newly identified in Brachypodium (Table 2).

Figure 1
figure 1

Pipeline for the prediction of miRNAs. miRNAs were identified by a series of filters. Black arrows are used to indicate the small RNAs that passed the criteria. The dashed arrow indicates small RNAs that failed to pass the criteria. Numbers of candidate miRNAs and precursors after each filter are indicated. Gray lines indicate where data from two pipelines has been combined. aStrict criteria are ≤4 mismatches + bulges in the miRNA:miRNA* pairing and a bulge size of ≤1 nucleotide. bLess strict criteria are defined as ≤8 mismatches + bulges in the miRNA:miRNA* pairing. cNon-redundant miRNAs from miRBase Release 18. dUp to two nucleotide mismatches from miRBase Release 19 were allowed. eConserved miRNAs include manually identified miR444 precursors and miR2218 precursors. miRNA, microRNA; nt, nucleotide; TP2M, transcripts per 2 million.

Table 2 Newly identified microRNAs

Because a pipeline approach designed to identify new miRNAs and precursors must be stringent enough to minimize the misannotation of false miRNAs, by design it will invariably miss some of the conserved miRNAs. Therefore, a homology pipeline with less strict criteria and manual analysis was developed to look for known miRNAs that did not follow all of the stringent criteria detailed in the first pipeline. Specifically, small RNA sequences that failed the first pipeline but matched previously annotated plant miRNAs were re-evaluated (Figure 1). This pipeline allowed less abundance and less strict stem-loop structures, and it had a lower cut-off value for abundance bias. Additionally, members of the miR444 family (Additional file 2: Figure S1) and miR2118 were added manually. This miRNA family has an intron in the precursor, which interferes with the folding filters in the pipelines. Together with conserved miRNAs identified from the previous, stringent pipeline, in total 75 conserved miRNAs generated from 115 precursors were identified. In total, we identified 116 mature miRNAs generated from 161 precursors in Brachypodium.

miR162 may not be functional in Brachypodium

We observed that miR162, one of the well-conserved miRNAs in plants, had quite low abundance in Brachypodium. In our data, miR162 was sequenced just once in a small RNA library from submergence-treated seedlings out of 17 libraries. We validated the lack of miR162 expression using an RNA blot (Figure 2a). While miR162 was highly expressed in rice and Arabidopsis, Brachypodium miR162 was not detected. One miR162 precursor was found in the Brachypodium genome (Figure 2b). Compared to the miR162 found in rice and Arabidopsis, the adenine (A) typically found at the 20th nucleotide, was replaced by a guanine (G) in Bdi-miR162 (Figure 2c). To see the impact of this substitution, we compared the miRNA:miRNA* pairing structure of authentic Bdi-miR162 with a hypothetical Bdi-miR162 with an A at the 20th nucleotide position (Bdi-miR16220A in Figure 2d). Substitution of the G by an A at the 20th nucleotide position of the hypothetical Bdi-miR162 (Bdi-miR16220A) increased the stability of the miRNA:miRNA* pairing. There was no bulge and the structure was similar to the stable structure of the Osa-miR162a:miR162* pairing in rice. Furthermore, there was no evidence of cleavage in its predicted target, the Brachypodium DCL1 homolog, in our PARE data. These data imply that the altered sequence of Bdi-MIR162 may affect miRNA processing and cause the reduced level of Bdi-miR162, resulting in little or no DCL1 cleavage.

Figure 2
figure 2

Characteristics of Bdi-miR162. (a) RNA blot of miR162 for Brachypodium (Bdi), rice (Osa) and Arabidopsis (Ath). Mixed probes of P32-end labeled Bdi-miR162 and Osa-miR162a were used to detect mature miR162 from Brachypodium panicle, rice panicle and Arabidopsis flower tissues. miR166 was used as a control. (b) Predicted secondary structure of the Bdi-MIR162 precursor. The location of Bdi-miR162 is indicated with a black line. (c) Sequence alignment of miR162 for Brachypodium, rice and Arabidopsis. Nucleotides that differ are shown in red. (d) Comparison of base pairings in the folded MIR162 precursor of Brachypodium and rice. Base pairing in a mutated Bdi-MIR162 precursor is shown between them. Position 20 of Bdi-miR162 is highlighted in yellow. miRNA is shown in red and miRNA* in blue. miRNA, microRNA.

Alternative folding generates microRNAs by stem-sharing

Clustering of miRNA genes is less common in plants than in animals [54]. About 10% to 15% of annotated miRNAs in rice, Arabidopsis and poplar are found in tandem clusters at a distance of less than 3 kb [55]. Most are polycistronic miRNA precursors that encode identical mature miRNAs. We also identified several miRNA precursors that are located in tandem in the Brachypodium genome. Interestingly, we noticed that two independent miRNA precursors overlapped each other in some miRNA clusters including Bdi-MIR9483ab, Bdi-MIR2118ab and Bdi-MIR166hj (Figure 3a,b,c). In these cases, one arm of a stem is shared by two adjacent miRNA precursors and serves as a common miRNA* strand. In these tandem-overlapping hairpins, the second hairpin structures are more stable and favored than the first. Thus, if the RNA secondary structure of the entire locus is predicted, the second miRNA precursor is preferred and the first one loses the hairpin structure. Our small RNA data indicate that both overlapping miRNA precursors may generate mature miRNAs. It is not clear whether Bdi-MIR9483a and Bdi-MIR9483b are equally functional since they encode identical mature miRNAs (Figure 3a). However, in the cases of the Bdi-MIR2118ab and Bdi-MIR166hj clusters, it is apparent that both of the precursors form stem-loop structures and are processed to generate their own mature miRNAs, due to distinct sequences of mature miRNAs and their expression. Two sets of mature miRNA were expressed at different levels in Bdi-MIR2118ab and Bdi-MIR166hj loci (Figure 3b,c). We hypothesized a model that alternative RNA folding of overlapped miRNA precursors generates miRNAs by sharing one arm of a stem (Figure 3d). In this model, the nascent transcript preferentially folds into the first stem-loop structure and is processed to generate the first mature miRNA. After transcription, both of the two possible RNA folding structures may fold to generate both of the respective mature miRNAs.

Figure 3
figure 3

Alternative RNA folding generates miRNAs by sharing one arm of a stem. Examples of alternative RNA folding in precursors of Bdi-miR9483 (a), Bdi-miR2118 (b) and Bdi-miR166 (c). Shown are schematic diagrams of adjacent miRNA precursors that share one arm of a stem. miRNA placement along the precursor is indicated with arrows, the middle arrow representing the shared arm of two miRNA precursors and the abundance of small RNAs matching each position is plotted below. Purple and red lines and filled circles indicate two mature miRNAs corresponding to each precursor and the blue line and filled circle indicates the miRNA*. Predicted secondary structures from each precursor, as well as that from the full sequence, are shown in the right panel. (d) Model for alternative RNA folding generating miRNAs by sharing one arm of a stem. The nascent transcript preferentially processed by DCL activity generates the first miRNAs (purple). After transcription, two possible RNA folding structures are responsible for generating the two mature miRNAs (purple and red). DCL, Dicer-like protein; miRNA, microRNA; TP2M, transcripts per 2 million.

Validation of stress-responsive microRNAs

To identify stress-regulated miRNAs, the sequencing frequency of miRNAs was compared between each stress-treated sample and its control library. We first examined how well-known and well-conserved examples of stress-regulated miRNAs are represented in our Brachypodium small RNA libraries. For example, miR399 and miR827, which are known to be induced under phosphate starvation conditions in Arabidopsis[56–60] and rice [34, 58, 61, 62], were highly represented in the phosphate-starved shoot library compared to its control library by about 9- and 4-fold, respectively. The upregulation in both cases was further validated by using the splinted ligation-mediated miRNA detection method (Figure 4a). Plants grown in a phosphate-depleted Murashige and Skoog (MS) agar media for 3 weeks accumulated both more miR399 and more miR827 than plants grown in a MS agar media.

Figure 4
figure 4

Validation of stress-responsive miRNAs. (a) Validation of phosphate-starvation inducible Bdi-miR399 and Bdi-miR827 by splinted ligation-mediated miRNA detection. Bdi-miR168 was used as a control. (b) Expression pattern of Bdi-miR5200 during submergence treatment and subsequent recovery. Three-week-old seedlings were subjected to submergence stress by completely immersing the plants in water, and subsequent recovery upon removal from the water. Shoots were collected before submergence treatment (0 d), after one-day (1 d) or two-day (2 d) treatments, and one (3 d) or three (5 d) days after recovery. Then 20 μg of total RNA were used to detect Bdi-FTL1 using an RNA blot. Bdi-miR168 and rRNAs were used as controls. miRNA, microRNA; Recov., recovery; rRNA, ribosomal RNA; Submer., submergence.

Next, we investigated miRNA expression using stress treatment and found that miR5200 was submergence-inducible; it was 80 times more represented in the submergence library compared to the control. To validate this regulation experimentally, the expression of miR5200 in response to submergence treatment and subsequent recovery was examined. Accumulation of miR5200 increased under submergence conditions and slowly decreased as the plants recovered from the submergence (Figure 4b). An FT (Flowering time locus T)-like gene, which is a predicted target of miR5200, was downregulated during the submergence treatment while miR5200 was upregulated. In addition, Bdi-FTL1 mRNA abundance increased during the recovery from the submergence. This result implies that miR5200 may play a role in regulating the flowering time under submergence conditions.

Identification of cleaved microRNA targets on a global scale

To investigate the targets of Brachypodium miRNAs, our strategy was to use the PARE approach to sequence deeply RNA decay intermediates that have a 5′ monophosphate and a 3′ poly(A) tail. This is the structure of the downstream fragments resulting from miRNA-guided cleavage. In Arabidopsis, PARE captured the vast majority of these cleavage products, which had been previously validated in independent experiments [44]. For Brachypodium, PARE libraries were constructed from RNA from four key organs. As shown in Table 3, these libraries were sequenced to a total depth of nearly 70 million reads. Of these, nearly 61 million matched the genome resulting in a non-redundant set of more than 5 million PARE sequences. These libraries represent a rich resource of PARE (or RNA degradome) data, which previously had been lacking for Brachypodium.

Table 3 Summary statistics of PARE libraries

To identify miRNA target cleavage events, we first generated a set of predicted target cleavage sites for 116 Brachypodium miRNAs and then examined all the distinct PARE sequences for those that precisely matched these predictions. Two prediction programs using different alignment algorithms were used to predict a large set of potential targets. Specifically, psRNATarget and CleaveLand (with the native and a modified scoring system) were used, as described in Materials and methods. These programs identified a total of 2,141 predicted miRNA targets in transcripts corresponding to 1,968 genes (loci) as indicated in Table 4. As shown in Additional file 2: Figure S2, psRNATarget and the two versions of CleaveLand each generated unique targets and thus a larger set of combined potential target cleavage sites then any one alone. When the total set of predicted targets cleavage sites were matched to the PARE data, 262 of the sites corresponding to 246 distinct coding regions (CDSs) had PARE sequences starting precisely at the predicted site, that is, between nucleotides 10 and 11 (from the 5′ end of the miRNA) of the miRNA-target RNA duplex, in at least one of the four tissue libraries (Table 4). The individual targets are listed in Additional file 1: Table S3, which also includes which programs predicted each target and the corresponding target scores.

Table 4 microRNA targets with PARE sequences at the expected site

A computational approach for characterizing microRNA target cleavage sites supported by PARE data

Since this study is the first in-depth analysis of cleaved Brachypodium miRNA targets, we sought to examine each target that had a PARE sequence precisely at the expected cleavage site for characteristics associated with the cleavage of well-studied miRNA targets in Arabidopsis. Hence, we developed a pipeline to evaluate cleaved miRNA targets for three parameters related to the prominence of the cleavage products, which is discussed in more detail below: abundance, rank and peak percentage. Our approach was first to establish cut-off criteria for each of these parameters using 137 Arabidopsis miRNA targets that had been experimentally validated for miRNA-guided cleavage, independent of PARE (Additional file 1: Table S4). The second step was to use these prominent criteria to characterize miRNA target cleavage events in Brachypodium.

We analyzed two sets of previously published Arabidopsis PARE data that were generated from the inflorescence tissues of Col-0 wild type (WT) and an xrn4 mutant [44]. XRN4 degrades the downstream fragments resulting from miRNA-mediated cleavages so the lack of this enzyme in the xrn4 mutant improved the sensitivity of PARE in Arabidopsis[44, 63]. The first parameter evaluated was the abundance of the PARE sequence at the target cleavage site. Often, miRNA-mediated cleavage at the precise cleavage site within a target mRNA occurs at a higher frequency than non-specific mRNA degradation events, and therefore the abundance of the PARE sequence at that site is higher than background. Figure 5a shows that most PARE sequences in the Arabidopsis samples had low abundance, while only approximately 2% to 7% of PARE sequences were more abundant than 10 TP10M. However, 73% to 82% of PARE sequences matching miRNA-mediated cleavage sites were at least 10 TP10M. Thus, a minimum abundance of 10 TP10M was chosen as the cut-off criterion for this parameter. The second parameter evaluated was the abundance rank within a given target mRNA (Figure 5b). Between 80% and 84% of the Arabidopsis miRNA-guided cleavages were either the most abundant or the second-most abundant (within the top two) PARE sequences in the annotated mRNA. An abundance rank of occurrence within the top two most abundant sequences within the mRNA was set as the cut-off for the second criterion of the prominence of the miRNA-guided target cleavage. The third and final parameter evaluated the pattern of cleavage, specifically, whether the abundance of the PARE sequence at the cleavage site was a high ‘peak’ above other degradation intermediates along the mRNA or was of similar abundance as the other PARE sequences from that RNA (Figure 5c). The peak percentage value was calculated as 100 times the sum of the top two PARE sequences from the target RNA, divided by the sum of all PARE abundances from that RNA. A value for ‘background’ degradation was first calculated from all annotated mRNAs in the WT library with PARE data. About 50% of all WT annotated transcripts had a peak percentage of less than 10%. For the validated miRNA targets with PARE cleavages, a minimum cut-off value of twice background, or a peak percentage of at least 20%, captured 61% to 74% of the known miRNA target cleavages in the WT and xrn4 backgrounds, respectively, and was therefore chosen as the cut-off criterion for this parameter.

Figure 5
figure 5

Examination of Arabidopsis PARE libraries using independently validated miRNA targets to establish criteria cut-offs. Comparisons of validated Arabidopsis miRNA targets and total genome-matched sequences in WT and xrn4 are shown. Histograms in (a) and (b) and lines in (c) show the percentages of PARE sequences matching miRNA target cleavage sites or all the PARE sequences from WT and xrn4 PARE libraries. Gray areas in (a), (b) and (c) indicate the selected cut-offs. (a) Distribution of PARE sequences by normalized abundance. (b) Distribution of PARE sequences by abundance rank in an annotated transcript for miRNA targets. (c) Distribution of PARE sequences by the peak percentage in annotated transcripts for miRNA targets or all genes. (d) A flowchart for evaluating miRNA target cleavage sites based on the three prominence criteria. amiRNA targets with PARE sequences at the precise canonical location. (e) Distribution of validated miRNA targets in WT and xrn4 according to level. miRNA, microRNA; PARE, parallel analysis of RNA ends; TP10M, transcripts per 10 million; WT wild type.

Using these three criteria, a pipeline was established to evaluate predicted miRNA targets that have PARE sequence data matching the precise predicted cleavage site, and to quantify the number of prominence criteria each fulfilled (Figure 5d). Targets meeting none of the criteria are designated Level 1, and targets meeting all the prominence criteria are designated Level 4. When the 137 previously validated Arabidopsis miRNA targets were evaluated with these criteria, 110 and 115 targets were identified (Levels 1 to 4) for WT and xrn4, respectively (Figure 5e). Although the number of known targets detected at any level for WT and xrn4 was similar, the number of Level 4 targets was greater (63 versus 49) and the number of Level 1 targets was fewer (13 versus 21) for xrn4, consistent with the increased accumulation of miRNA cleavage products in the mutant [44, 64]. All of the Arabidopsis miRNA targets examined were previously experimentally validated and most fulfill multiple criteria. However, it is apparent that all levels contain bona fide Arabidopsis miRNA targets (Figure 5e).

PARE data support precise cleavage guided by most BrachypodiummicroRNAs, with higher prominence levels evident among annotated conserved microRNAs

Next, we sought to use the same pipeline outlined in Figure 5d to further characterize the precise target cleavage sites detected with PARE for Brachypodium. Specifically, 262 sites in 246 predicted miRNA targets from Table 4 that had PARE data indicating precise target cleavage were examined. The level designation for each and the individual criteria passed are listed in Additional file 1: Table S3. Data for the complete set of previously annotated and newly identified miRNAs from Figure 1 are compiled in Figure 6. Even without the benefit of an xrn4 mutant, which is not available for Brachypodium, a high proportion (81%) of the 116 Brachypodium miRNAs had PARE sequences at the precise site(s) predicted for target cleavage. Among annotated miRNAs, 85% had PARE cleavage support whereas 70% did among those newly identified in this study. When the individual cleavage sites for the multiple targets predicted for each miRNA are considered, 42 sites were Level 4 targets, and 38 were level 2 or 3. The Level 1 cleavage sites predominated for the Brachypodium miRNA targets with precise PARE data, particularly those of the newly identified miRNAs. Level 1 did not predominate when PARE levels were compared for the Arabidopsis target sites used to establish the cut-offs (Figure 5e). However, the target cleavage sites examined for Arabidopsis were biased toward those that were easier to detect because only those that had prior experimental validation were examined by PARE.

Figure 6
figure 6

Target prominence level distribution among precise cleavage sites with PARE data. The number of distinct cleavage sites of each of the levels for (a) annotated miRNAs and (b) newly identified miRNAs are shown. miRNAs are listed by families in ascending order as in the supplemental tables. Non-conserved miRNAs are indicated with # in (a). miRNA, microRNA; PARE, parallel analysis of RNA ends.

An important observation made above for the level designations of Arabidopsis miRNA targets was also observed for those of Brachypodium; all prominence levels contain valid miRNA targets. For Brachypodium, this is indicated by the examples in the degradation plots (D-plots) in Figure 7, which indicate the location and abundance of the PARE sequence at the precise cleavage site as well as all PARE sequences matching these target RNAs. Although the PARE sequences at the Level 1 and Level 2 miRNA cleavage sites are not as prominent as those in the Level 3 and Level 4 examples, they are homologous to experimentally validated targets of conserved miRNAs in other plant systems. Additional D-plots for Level 1 to 4 targets are shown in Additional file 2: Figure S3. These data indicate that the level numbers represent a prominence score that is useful for characterizing miRNA target cleavage, and even a Level 1 cleavage at the precise predicted site for a miRNA should be considered quite strongly as evidence for cleavage.

Figure 7
figure 7

Brachypodium PARE data demonstrate miRNA-guided target cleavage. Examples of degradation plots for miRNA targets for each prominence level. The abundance of all the PARE sequences matching miRNA targets are plotted with blank diamonds. The PARE sequence corresponding to the miRNA-mediated cleavage is represented by a black diamond. (a) Bradi1g52240, encoding a GIBBERELLIC-ACID INSENSITIVE (GAI), REPRESSOR of GAI (RGA), and SCARECROW (SCR) (GRAS) transcription factor, is a Level 4 target of Bdi-miR171bcef. (b) Bradi3g05510, encoding a SQUAMOSA PROMOTER BINDING PROTEIN-LIKE (SPL) transcription factor is a Level 3 target of Bdi-miR156b-ij. (c) Bradi4g39330, encoding a laccase, is a Level 2 target of Bdi-miR397ac. (d) Bradi3g51077, encoding argonaute 1a, is a Level 1 target of Bdi-miR168. nt, nucleotide; miRNA, microRNA; PARE, parallel analysis of RNA ends; TP10M, transcripts per 10 million.

Analysis of PARE data for identifying microRNA cleavage-initiating phased small RNAs

A small group of miRNAs initiate the production of phased siRNAs from their target loci. The cleavage guided by the miRNA defines the beginning of the phase register, which can have a 21- or 24-nucleotide periodicity. Examples of miRNAs initiating 21-nucleotide phased loci are miR173 [19], miR390 [25] and miR2118 [30], whereas miR2275 [30] is the only miRNA currently known that initiates 24-nucleotide phased loci. It has previously been reported that Brachypodium contains hundreds of 21- and 24-nucleotide phased loci in its panicles [31]. Since Brachypodium produces miR390, miR2118 and miR2275, all known for initiating phased loci, these miRNAs may initiate some or all of these previously reported loci.

PARE data captures the downstream products of miRNA-mediated cleavage and thus could help identify the location of the cleavages initiating these phased loci and the identity of the miRNAs responsible. Using the same method as we used previously [31], over 1,100 phased loci were identified in leaf, panicle, stem and root small RNA libraries (as detailed in Materials and methods). The PARE sequences corresponding to each locus +/-5 phase cycles (105- or 120-nucleotide for 21- or 24-nucleotide periods, respectively) were collected. The 13 nucleotides upstream of the PARE sequence start was added to the first 11 nucleotides of the PARE sequence to create the miRNA-complementary site ±1. These miRNA-complementary sites were compared to Bdi-miR390, Bdi-miR2118 and Bdi-miR2275, and scored such that 1 was added for each mismatch and 0.5 was added for each G:U pair. Duplexes with a score of 5 or less were considered targets of the respective miRNA. PARE sequences were associated with 945 phased loci. The miRNAs initiating 437 of the phased loci were identified (Table 5). For all but six of the loci, the PARE sequence generated by that respective miRNA was within ±1 of the dominant phase register. An example of the PARE and small RNA sequence distribution from representative loci initiated by miR2118 is shown in Additional file 2: Figure S4. Identifying the miRNAs initiating approximately 500 phased loci showed that PARE data can be a valuable tool for analysis of these loci and their cleavage sites.

Table 5 PARE data identifies miRNA-initiated phased loci

Only about half of the phased loci were accounted for by miR390, miR2118 and miR2275. Since it has been reported that 22-nucleotide miRNAs can potentially initiate phased loci [27, 28], the remaining seven 22-nucleotide miRNAs from Additional file 1: Table S2 were examined for their potential to target the phased loci. PARE data indicates that Bdi-miR5163b-3p targets one locus, Bradi4g10171, encoding an NB-ARC (nucleotide binding human APAF-1, plant R protein, Caenorhabditis elegans CED-4) domain-containing R protein (Table 5, Figure 8). Phasing at this locus was found in the leaf, panicle, stem and root. Figure 8a shows the complementarity between Bdi-miR5163b-3p and Bradi4g10171. The PARE and small RNA were not evenly distributed across the gene, but rather concentrated after the PARE sequence indicative of Bdi-miR5163b-3p-mediated cleavage (Figure 8b,c). The position with the maximum phase score coincides with the exact position of Bdi-miR5163b-3p cleavage (Figure 8d), and nearly all of the small RNAs within the phased region precisely coincide with the phase register (Figure 8e). Taken together, this shows that Bdi-miR5163b-3p is a novel miRNA that initiates phasing of its target locus.

Figure 8
figure 8

Evidence that Bdi-miR5163b-3p-guided cleavage induces phased siRNAs from an NB-ARC mRNA. (a) The sequence alignment of the Bradi4g10171, an NB-ARC gene, and Bdi-Bdi-miR5163b-3p. Lines indicate a perfect match. The circle indicates a G:U pair. The black arrow indicates the cleavage site with the TP10M PARE abundance listed above. (b,c) PARE and small RNA distributions for Bradi4g10171. The region shown is 9771971 to 9776753 on chromosome 4 and the gray boxes on the x-axis are the exons. Gray diamonds and circles represent PARE and small RNA sequences, respectively. The filled black diamond at the position indicated by the arrow is the PARE site generated by Bdi-miR5163b-3p cleavage. (d) The phasing score profile for the phased region. (e) The small RNA abundance profile for the phased region. PARE, parallel analysis of RNA ends; siRNA, small interfering RNA; TP10M, transcripts per 10 million; TP2M, transcripts per 2 million.

It has been proposed that phasiRNAs are often generated from protein-coding loci by miRNA-triggered cleavage and function to suppress gene families, such as those of the PPR, NB-LRR and MYB genes in Arabidopsis, Medicago, apple, and peach [21, 23, 65–67]. The aforementioned NB-ARC target of Bdi-miR5163b-3p is a member of the NB-LRR family of disease resistance genes as are several other predicted NB-ARC targets of this miRNA. Although we were not able to obtain cleavage evidence for the others from PARE data, all the NB-ARC targets of Bdi-miR5163b-3p are phasiRNA loci. It is notable that Bdi-miR5163b-3p does not share sequence similarity with miR472/482, miR2118, miR1507 and miR2109/miR5213, which are known to target NB domains of NB-LRR and trigger phasiRNA generation [21, 22, 24, 67]. In addition, Bdi-miR5163b-3p targets the LRR domain instead of the NB domain. Bdi-miR5163b-3p has not been identified in other plants yet. These results imply that Bdi-miR5163b-3p may be a Brachypodium-specific miRNA that regulates NB-ARC genes and plays a role in disease resistance.

PARE identifies differential cleavage by microRNA family members

Many miRNA family members have slight differences in sequences and differential expression patterns, which could cause their specificity for target regulation [34]. We investigated differential regulation of targets by miRNA family members using PARE data. We observed that two miRNA families, miR166 and miR156/529, have members with both sequence variation and differential regulation patterns. In Brachypodium, ten MIR genes encode miR166. While six MIR166 genes generate identical guide sequences, MIR166h, MIR166j and MIR166f give rise to miRNA variants with slight sequence variations. Interestingly, the three miR166 sequences showed differential expression patterns in diverse tissues (Figure 9a). The most abundant sequence, miR166a-eg, was ubiquitously expressed except in mature leaf tissue of six-week-old plants. However, miR166f and miR166h exhibited tissue-preferential expression in roots and panicles, respectively. The expression pattern of miR166j was not distinguishable due to sequence similarity with miR166h. Given the clustered structure of MIR166h and MIR166j (Figure 3c), it is likely that miR166j expression is similar to miR166h. To address whether the differentially expressed miR166 family members also differentially regulate their target genes, we examined the cleavage events in a target gene encoding a homeobox domain-leucine zipper (HD-ZIP) transcription factor. As shown in Figure 9a, the target gene has two major PARE sequences corresponding to cleavage sites by miR166a-eg, miR166h, miR166j and miR166f. Since miR166f has a 5′-end two-nucleotide offset relative to the other miR166 family members, the miR166f-mediated cleavage site is also shifted relative to the cleavage site mediated by other family members. In root tissue, two different PARE sequences were mapped to the cleavage sites of both the root-preferential miR166f and the remaining miR166 family members. However, the panicle library showed only one prominent cleavage site by miR166a-eg, miR166h and miR166j, while miR166f-mediated target cleavage was detected at only a basal level. This implies that miR166 members may regulate HD-ZIP gene expression in a tissue-preferential manner.

Figure 9
figure 9

miRNA family members show different patterns of tissue-preferential expression and target cleavage. (a) miR166 family-mediated cleavage of Bradi1g13910, which encodes an HD-ZIP III transcription factor. (b) miR156 family-mediated cleavage of Bradi4g33770, which encodes a SPL transcription factor. The schematic diagrams of genomic structures for miRNA target genes are shown. The gray boxes indicate exons. The sequences of miRNA target sites and miRNA family members were aligned. Lines indicate a perfect match. Circles indicate G:U pairs. PARE sequence distribution for target sites are shown with histograms. The expression patterns of miRNA family members were examined for diverse tissues using the splinted ligation-mediated miRNA detection method. Bdi-miR168 was used for the control. bp, base pair; HD-ZIP, homeobox domain-leucine zipper; SPL, Squamosa Promoter Binding Protein-Like; miRNA, microRNA; PARE, parallel analysis of RNA ends; TP10M, transcripts per 10 million.

It has been reported that the regulation of OsSPL14 expression by Osa-miR156 and Osa-miR529 may play a crucial role in both tillering and panicle-branching in rice [34, 68, 69]. To explore this regulation in Brachypodium, we first examined the expression pattern of Bdi-miR156 and Bdi-miR529. Bdi-miR156 was highly expressed in young shoot tissues, and downregulated in mature leaf tissues (Figure 9b). In addition, Bdi-miR156 expression in panicle tissue was much lower than that in the other tissues, whereas Bdi-miR529 was preferentially expressed in panicle tissue. These expression patterns are similar to those in rice [34, 70]. Next, we addressed the effect of these two miRNAs on a target gene encoding BdiSPL14L (Brachypodium distachyon Squamosa Promoter Binding Protein 14-Like, an OsSPL14 homolog in Brachypodium), using PARE data. In leaf tissue of three-week-old plants, a prominent PARE sequence was mapped to a cleavage site by miR156, while the PARE sequence from the site of miR529-mediated cleavage had a very low abundance. In contrast, the PARE sequence corresponding to the site of cleavage mediated by panicle-preferential miR529 was more abundant than that of miR156 in the panicle library. This result implies that cleavage of BdiSPL14L during vegetative growth is predominantly miR156-guided, whereas miR529 and miR156 regulate the same gene during reproductive development. In conclusion, the two examples of miRNA families shown to be differentially expressed and to guide differential target cleavage highlight the potential importance, in general, of distinct miRNA family members for post-transcriptional control.

Discussion

In this study, we have provided a better understanding of Brachypodium miRNAs in several ways. The deep sequencing of 17 small RNA libraries as well as four PARE libraries from various tissues and environmentally stress-treated plants allowed for the identification of many new miRNAs and the characterization of their roles in target cleavage. The miRNAs identified include both conserved miRNAs that have not been reported in Brachypodium, as well as non-conserved miRNAs that have not been found in other plants and thus Thus, these miRNAs may function in gene regulatory mechanisms specific to Brachypodium or those shared with related species. Our analysis uncovered a total of 81 new miRNA precursors of conserved and non-conserved miRNAs. These confirmed 11 conserved and five non-conserved miRNA precursors that were recently reported [42]. In addition, evidence that miR162-mediated DCL1 regulation may not function in Brachypodium was presented. We also identified a new mode of miR5200 regulation in response to submergence stress. More than 260 targets of new and known miRNAs with PARE sequences at the precise cleavage site were identified and characterized with a set of prominence criteria. Validated miRNA targets were found at all prominence levels and even low prominence was found to be characteristic of some targets. Combining PARE data with the small RNA data identified the miRNAs responsible for initiating approximately 500 phased loci, including one of the novel miRNAs, Bdi-miR5163b-3p. PARE data also supported differential target cleavage in various tissues by miRNA variants.

Null function of miR162-mediated DCL1 regulation in Brachypodium

It has been proposed that miR162-mediated negative feedback regulation of DCL1, encoding a protein responsible for miRNA biogenesis, is a widely conserved mechanism in flowering plants [71, 72]. In Arabidopsis, DCL1 is targeted not only by miR162, encoded at two distinct loci, but also by miR838, generated from a stem-loop in intron 14 of the DCL1 pre-mRNA [73]. Both miRNAs repress DCL1 mRNA expression, while DCL1 is required to process those miRNAs. This negative feedback mechanism may be required for a sophisticated homeostatic regulation of miRNA processing and subsequent target regulation in flowering plants. Even though miR838 seems to be an Arabidopsis-specific miRNA, miR162 is widely conserved and identified in most flowering plant species, including Arabidopsis, poplar, Medicago, tomato, rice, maize and others. However, miR162 expression was found to be extremely low in wheat [39, 74], barley [75] and sorghum [76]. Our data, as well as a recent publication [39] indicate extremely low expression of miR162 in Brachypodium despite the presence of a miR162 precursor. Among the 17 libraries reported here, miR162 was only sequenced once and was beneath detection when examined by RNA blot (Figure 2). We also could not find any PARE sequence corresponding to a Bdi-miR162 target site in a Brachypodium gene that is homologous to DCL1 (Bradi1g77087.1).

We identified a one-nucleotide difference in miR162 compared to the conserved sequence and hypothesized that this change may affect the stability of miRNA:miRNA* pairing, resulting in inefficient processing of the duplex. We cannot rule out other possibilities like reduced stability, less efficient nuclear export or inefficient AGO1 loading of the miRNA:miRNA* duplex. However, the lack of PARE sequences corresponding to the MIR162 locus supports the inefficient processing of the precursor by DCL1. We also cannot rule out the possibility that the MIR162 promoter is not active in the tissues analyzed in Brachypodium and its relatives. It has been reported that altered miRNA:miRNA* base pairing can affect miRNA processing in viruses, animals and plants [77, 78]. In animals, those mutations are associated with cancer and other diseases [79]. In Arabidopsis, it has been reported that a natural variation in miR164a* from the C24 ecotype affects efficient processing of miR164a, resulting in altered leaf serration and an accessory bud formation phenotype in the C24 ecotype [80]. The functional effect of low expression of miR162 in Brachypodium is still unclear. In the moss Physcomitrella patens, Ppt-miR1047, embedded in intron 7 of PpDCL1a, targets a host gene, PpDCL1 mRNA [72]. Despite the sequence and position differences between Ath-miR838 and Ppt-miR1047, they appear to have analogous functions in the negative feedback regulation of DCL1 pre-mRNAs. This finding suggests that the miRNA-DCL1 feedback loop is critical for control of miRNA function in land plants. Therefore, it is likely that Brachypodium and its relatives may have some other mechanism to regulate DCL1 gene expression, or have evolved to control DCL1 function without a miRNA-DCL1 feedback loop.

miR5200 regulates gene expression of a flowering time regulator under submergence conditions

We have identified a submergence-inducible miR5200 that is predicted to target two genes encoding FT-like genes in Brachypodium. Interestingly, Bdi-FTL1 was downregulated during submergence when miR5200 is induced (Figure 4b). It has been proposed that flowering is inhibited by submergence both in Arabidopsis and rice [81]. Ectopic expression of the rice SUBMERGENCE1A (SUB1A) gene in Arabidopsis led to flowering inhibition by reducing the level of CONSTANS (CO) and FT mRNAs. In addition, ectopic over-expression or submergence-induced expression of SUB1A in rice affected the abundance of mRNAs encoding Hd1 and Hd3a, which encode two conserved flowering regulators of CO and FT, respectively. It is possible that Brachypodium also contains a similar regulatory mechanism where submergence represses flowering time by downregulation of Bdi-FTL1. In addition, miR5200 may also play a crucial role in downregulation of Bdi-FTL1. However, miR5200 was not identified in either rice or Arabidopsis. Our miRNA conservation analysis (Additional file 2: Figure S5) showed that wheat and barley express miR5200 or related small RNAs at low levels. Thus, miR5200 in wheat and barley also could be upregulated under submergence and target FT-like genes.

Known and potential outcomes of modulating miR156/529 expression

Recently, it has been shown that over-expression of miR156 in rice, maize, switchgrass and Brachypodium caused morphological alterations including an increased number of branches, resulting in a favorable phenotype for improved biomass production [70, 82–84]. However, over-expression of miR156 can also result in other defects like dwarfism and male sterility, likely because miR156 represses many SPL genes. Increased expression of OsSPL14 due to sequence alteration of the miR156 target site leads to less tillers and more panicle branching, resulting in increased grain yields in rice [68, 69]. Additionally, miR529, which is in the same miRNA family as miR156, has been proposed to play a crucial role in panicle branching [34]. Sequence variations and differential expression patterns between miR156 and miR529 could cause their distinct function in SPL gene repression. In our data, Bdi-miR156 and Bdi-miR529 also showed a similar differential regulation of their target gene, Bdi-SPL14L, when they were expressed (Figure 9b). In contrast to miR156, the effects of over-expressing miR529 have not been examined yet for transgenic plants. Since Bdi-miR529 and Bdi-miR156 share several target genes and may differ for others, their over-expression should be significant and, depending on the outcome, may suggest strategies to improve food or bioenergy crops.

Power of PARE to identify cleaved targets of BrachypodiummicroRNAs

The PARE data we examined in this study represents a major advance because, prior to this work, the cleavage of miRNA targets had only been demonstrated by examining a limited set of Brachypodium miRNA targets [39, 53]. The distinct cleavage of Bdi-SPL14L guided by Bdi-miR156 and Bdi-miR529 as discussed above, was precisely pinpointed at sites four nucleotides apart using PARE (Figure 9b), reminiscent of similar observations from rice PARE data [34]. Cleavage guided by Bdi-miR166 family members in an HD-ZIP III mRNA at sites two nucleotides apart was also validated with PARE (Figure 9a). In total, PARE evidence for precise cleavage at more than 263 miRNA cleavage sites was obtained. The targets with PARE cleavage evidence corresponded to 70% of the new miRNAs and 85% of the annotated miRNAs.

We found that using two distinct target prediction programs, CleaveLand and psRNATarget, was a more powerful way to predict more miRNA targets for cleavage validation with PARE than either algorithm alone (Additional file 2: Figure S2). Accordingly, this larger pool enhanced the number of cleaved targets precisely identified with PARE. We developed a characterization tool to examine the prominence of the PARE sequence at the precise cleavage site of these targets using criteria established based on Arabidopsis miRNA targets with independent experimental validation. This tool proved useful in several ways. First, it identified the PARE-validated targets that had the most prominent target cleavage data as exemplified by the D-plots for Level 3 or 4 targets (Figure 7 and Additional file 2: Figure S3). Second, since experimentally confirmed miRNA targets are present at all levels, by knowing which criteria a particular Level 2 or 3 target has passed (as listed in Additional file 1: Table S3), potential explanations and future experiments can be identified. For instance, a target with a D-plot exhibiting many decay intermediates towards the 3′ end or throughout the RNA may explain why the precise target may not be in the top two most abundant PARE sequences or pass the peak percentage criteria. Alternatively, the PARE sequence at the cleavage site may be of low abundance. This could be due to low accumulation of the miRNA or the target RNA in the samples examined, or because the miRNA largely inhibits translation and cleavage is less prominent. Cleavage can be detected for translationally regulated miRNA targets as is the case for miR172 and miR398 targets [16, 85]. Moreover, for some targets, a low prominence level may be a conserved characteristic of the target RNA. This is clearly the case for AGO1. Cleavage of AGO1 mRNA is guided by miR168 and AGO1 is a Level 2 target in Arabidopsis and a Level 1 target in rice and Brachypodium (Figure 7d and Additional file 2: Figure S6). As PARE libraries are sequenced from diverse samples and closely related species, some of the Brachypodium targets of new or annotated miRNAs may be classified in higher levels or found to be characteristically Level 1. The PARE data contributed by this work is extensive but likely will be augmented with additional Brachypodium PARE data in the future.

Conclusions

The knowledge gained from the identification of new miRNAs and cleavage-validated miRNA targets enhances our understanding of regulatory pathways in Brachypodium and should inform similar studies in related plants. The novel regulation of miRNA expression detected under submergence conditions is one example. Another is the differential regulation of closely related miRNA family members, which correlates with differential target cleavage detected with PARE. The PARE data contributed by this study also contains RNA decay intermediates derived from mechanisms independent of miRNAs and thus should prove useful for studying other components of the RNA degradome.

Materials and methods

Plant growth and stress treatment

Brachypodium (Line Bd21) was used for all libraries. For the BDI small RNA and PARE libraries, seeds were germinated in soil and the plants were grown in a growth chamber at a constant 20°C, under a 20 h light, 4 h dark cycle. Leaf, stem and panicle tissues were harvested from six-week-old plants. Leaf tissue was taken from the top of the blades of the first and second leaves, not including the leaf sheath. Stem tissue was harvested from just above the first top node to just below the second node, including both nodes as well as the leaf sheath of the second leaf. Panicle tissue was harvested from plants with immature, unfertilized flowers from the tip of the spikelet to just above, but not including, the panicle node at approximately 10 h into the subjective day. To collect root tissue, seeds were husked, surface-sterilized with a 50% bleach solution, rinsed with sterile water and germinated on half-strength Murashige and Skoog (MS) medium containing 0.2% (w/v) Phytagel in sterile plastic boxes. Root tissue was harvested from three-week-old plants from the tip of the root to just below, but not including the seed.

For the OBD set of small RNA libraries from tissues, Brachypodium plants were grown with supplemental light for a 16 h photoperiod at 25°C. For OBD02, leaf and stem tissues were collected at days 7 to 8, 14 to 15, and 27 to 28 and pooled. Panicle tissue was collected at days 28 to 35 at four times: 0700 (dawn), 1300, 1900 and 0100 hours. For OBD03, whole seedlings including roots were pooled. Two samples were grown in soil under a 16 h photoperiod and tissues were collected from three- and four-day-old plants. Two additional samples included five-day-old etiolated seedlings and light-grown seedlings grown on moistened filter paper.

Stress treatments and controls except phosphate starvation and its control for BDI small RNA libraries were conducted using three-week-old plants planted in soil under normal growth conditions with a 20 h photoperiod as above. Stresses began 5 h to 6 h after the lights came on. For heat stress, plants were moved to 40°C for 12 h. For cold stress, plants were moved to 4°C in a cold room for 12 h. For drought stress, plants were removed from the soil and rinsed, blotted dry with paper towels, and air-dried for 12 h. For salt stress, the base of the plant pot was submerged in 300 mM NaCl solution for 12 h to a depth of approximately 1 inch, without covering the soil. For submergence stress, plants were fully submerged in 20°C water for 24 h. For phosphate starvation stress, seeds were planted in MS media without KH2PO4, and control plants were planted in MS media; plants were harvested at three weeks, 5 to 6 h after lights on. Control plants were sampled at the same time as the stress treatments.

For OBD04, tissues from plants treated with abiotic, biotic and mechanical stresses were pooled. For abiotic stress, leaves and stems were collected from plants subjected to heat (37°C), cold (4°C), salt (200 mM NaCl) and osmotic (300 mM mannitol) treatments. For biotic stress, leaves were collected from plants with the following treatments. For the bacterial treatment, Pseudomonas fluorescens and Pseudomonas syringae pv. tomato DC3000 (Pto DC3000) were delivered to 28–29-day-old Brachypodium leaves under vacuum at 1 × 107 cfu/mL. Leaves of four-week-old plants were placed in a 1.5 ml microcentrifuge tube containing 20 μg/mL victorin. Tissue was harvested 24 h after treatment with victorin. Mechanical stress was applied by crushing the leaves of 29–30-day-old plants with tweezers. Unless otherwise indicated, libraries were made from a pool of samples collected every 4 h after the stress began, that is, 4 h + 8 h + 12 h + 16 h + 20 h + 24 h.

All collected tissues were immediately frozen in liquid nitrogen and stored at -80°C until RNA extraction.

Small RNA library construction and sequencing

The BDI set of small RNA libraries were constructed as described [31, 86] with minor modifications. Briefly, total RNAs were isolated using the Trizol reagent (Molecular Research Center), and low molecular weight RNA was purified from total RNA by PEG8000/NaCl precipitation. Following denaturing urea PAGE purification of the 20- to 30-nucleotide components, small RNAs were ligated first with the 5′ RNA adapter (5′ GUUCAGAGUUCUACAGUCCGACGAUC 3′) and then the 3′ RNA adapter (5′ PUCGUAUGCCGUCUUCUGCUUG-idT 3′). In each step, the ligated products were purified with PAGE gel. After the first strands of cDNAs were synthesized using RT primer (5′ CAAGCAGAAGACGGCATACGA 3′) and purified with PAGE gel, the final bands were amplified using 18 cycles of PCR using PCR primers (forward PCR primer: 5′ AATGATACGGCGACCACCGACAGGTTCAGAGTTCTACAGTCCGA 3′, reverse PCR primer: 5′ CAAGCAGAAGACGGCATACGA 3′) and purified with PAGE gel for sequencing. The OBD set of small RNA libraries was constructed as described [31]. Libraries were sequenced on an Illumina GAIIx instrument.

Computational analysis of small RNA sequencing data and microRNA identification

The raw sequencing data was trimmed by removing adapter sequences, and mapped to the Brachypodium genome 8X assembly [87] using Bowtie [88]. The two pipelines used for the computational analysis of the small RNA libraries and the identification of miRNAs is outlined in Figure 1. For the stringent criteria pipeline, small RNA sequences 20 to 24 nucleotides in length, with an abundance of ≥10 TP2M in at least one library, and ≤20 total genome hits, were evaluated for potential pairing of miRNA and miR* using a modified version of miREAP optimized for plant miRNAs [21, 34]. Strand bias was calculated as the total abundance of sequences matching the sense strand divided by the total abundance of all sequences matching both strands. Abundance bias was calculated as the sum of the abundances of the top two sequences for the precursor divided by the total abundance of all matching sequences. Sequences with a strand bias of ≥0.9 and an abundance bias of ≥0.7 were selected and the stem-loop structure was predicted using UNAFold [89]. Precursor sequences with ≤4 mismatches and bulges in the miRNA:miRNA* pairing, and a bulge size of ≤1 nucleotide were selected for final manual analysis. Each precursor structure was inspected manually to give visual confirmation of the secondary structure, that the predicted miRNA sequence was the most abundant sequence variant at the precursor and that additional miRNA sequences originate from the same precursor. A custom Perl script for this pipeline is available as Supplemental Dataset 1 in Jeong et al. [34].

For the sequence homology pipeline, annotated miRNAs from 47 flowering plant species in miRBase Release 18 were mapped to the Brachypodium genome 8X using Bowtie [88] with zero mismatches. miRBase sequences 20 to 24 nucleotides in size with ≤20 genome hits were then matched to the small RNA libraries to identify conserved miRNAs. Further conserved candidates were selected from the small RNA sequences that had an abundance ≥1 TP2M in at least one library and did not pass the stringent criteria pipeline. If multiple sequences mapped to the same genome position, the most abundant sequence was chosen. A putative precursor was selected from sequences 150-nucleotide upstream and 150-nucleotide downstream from the 5′ end of the miRNA. For sequences within 100 nucleotides of each other, the most abundant miRNA was selected for choosing the precursor location. Stem-loop structures of the putative precursors were predicted using UNAFold [89], and precursor sizes and positions were adjusted manually to achieve the best folding structure. Sequences with ≤8 mismatches plus bulges within the miRNA/miRNA* pairing were then evaluated for strand bias and abundance bias with cut-offs of 0.9 and 0.4, respectively. Finally, precursors were evaluated manually for additional miRNA candidates.

Identified miRNAs were matched to miRBase Release 18 with up to two mismatches allowed for naming. Sequences that did not match any known miRNA were considered novel; those sequences only found in Brachypodium were considered non-conserved and annotated, and those sequences matching at least one additional species were considered conserved miRNAs.

Additionally, miRNAs with known family members that did not pass any of the pipelines were manually searched for in miRBase Release 19 using the criteria of up to three mismatches.

microRNA detection using RNA blots

Low molecular weight RNA was isolated from 150 μg of total RNA by PEG/NaCl precipitation and separated on 15% denaturing PAGE gels, transferred to Hybond N+ nylon membranes (Amersham Biosciences), and fixed by UV cross-linking. The DNA oligonucleotide probes listed in Additional file 1: Table S4 were purified and hybridized [34].

microRNA RNA detection with splinted ligation-mediated miRNA detection

miRNAs were detected using the miRNAtect-It miRNA labeling and detection kit (Affymetrix) as in [90–92]. Specific bridge oligonucleotides shown in Additional file 1: Table S5 were ligated to both 5 to 20 μg of total RNA and a P32-labeled detection oligo. Products were separated on a 15% denaturing PAGE gel and visualized using a Typhoon Phosphoimager.

mRNA detection using RNA blots

To detect mRNA, 20 μg of total RNA was loaded on 1.5% formaldehyde agarose gels, transferred to Hybond N+ nylon membranes, and probed with P32-labeled complementary DNA sequences generated by PCR using the primers listed in Additional file 1: Table S5 and the Church and Gilbert hybridization buffer [93].

PARE library construction and sequencing

PARE libraries were constructed as described [94] with minor modifications. Briefly, total RNA was extracted with Trizol Reagent (Molecular Research Center) and PolyA RNA was isolated. A 5′ RNA adapter was added to ligation-competent RNA molecules containing a 5′ monophosphate, and the resulting products were reverse-transcribed into cDNA using an oligo(dT) primer with a 3′ adapter sequence, and amplified by PCR. In each step the ligated products were purified with PAGE gel. The resulting product was digested with MmeI to capture 20- or 21-bp fragments. MmeI cleavage products were ligated with a double-stranded DNA 3′ adapter with degenerate overhangs and amplified by PCR. Final bands were purified with PAGE gel and used for high-throughput sequencing. Libraries were sequenced on an Illumina GAIIx instrument.

Identifying microRNA targets

Putative miRNA targets were identified using three target prediction schemes. Using the CleaveLand v2.0 suite of software [95], miRNAs were matched to Brachypodium coding sequence CDSv1.2 [96] using SHRiMP V2.2.0 [87], and miRNA targets with a target score ≤7 were identified using the Perl script described previously [97]. A second, modified version of CleaveLand (CleaveLandM) was run, again using SHRiMP V2.2.0, to match to the coding sequence, but modified to use the scoring system detailed in [98]. Targets were identified with a score ≤4. Finally, targets were predicted using psRNATarget [99] with the following parameters: the maximum expectation was 5.0, length for complementarity scoring (hspsize) was set to the length of the miRNA, and the range of central pairing was 10 to 11 [100].

Characterizing microRNA targets in PARE libraries

PARE sequences were trimmed of adapters and down to 20 nucleotides. PARE sequences were then matched to the Brachypodium genome 8X assembly and CDS (MIPSV1.2) using Bowtie. Sequences matching the genome more than 20 times, as well as regions of single-, double- or triple-nucleotide sequences >12 nucleotides in length were removed. PARE libraries were normalized to TP10M by multiplying by a factor of 10,000,000 divided by the number of genome-matched sequences Using custom-written Perl scripts, predicted miRNA targets were matched to the PARE data from each of the four tissue libraries in Table 3, and PARE sequences that paired with the miRNA at the precise canonical cleavage site were identified. For each PARE sequence passing these filters, the data from the library in which the precise PARE sequence at the target cleavage site had the highest abundance was used for further characterization using the prominence criteria outlined in the text.

Previously published Arabidopsis inflorescence wild-type (Col-0) and xrn4 mutant PARE libraries [44] were exactly matched to the TAIR9 genome and cDNAs. Sequences matching the genome more than 20 times, sequences matching rRNAs, tRNAs, small nuclear RNAs and small nucleolar RNAs, and sequences with single-, double- or triple-nucleotide repeats >12 nucleotides in length were removed. PARE parameters, including abundance, peak percentage and abundance rank in the transcript, were established for the 137 previously reported Arabidopsis miRNA targets listed in Additional file 1: Table S4, and criteria for the Brachypodium PARE libraries was established from parameters of the Arabidopsis miRNA targets. Abundance is the abundance of the PARE sequence at the exact miRNA cleavage site. The peak percentage was calculated as the abundance of the PARE sequence exactly matching the miRNA cleavage site divided by the sum of all PARE sequences matching the transcript times 100. The abundance rank was determined by listing all PARE sequences matching the transcript in order from largest to smallest.

Phasing analysis

The 21- and 24-nucleotide phased loci were identified using Perl scripts as previously described [31] with a minimum phasing score ≥25. The small RNA libraries BDI02, BDI04, BDI05, BDI06, BDI08, OBD01, OBD02 and OBD03 were used for this analysis. To set the boundaries of the phased loci, the libraries were combined and the 5′ end was set to the position where the phase score first was ≥20 in any library. The 3′ boundary was set at the last phase window with a score ≥20. The PARE sequences within five phase periods of these loci were pulled from the BDI25 PARE library and scored for complementarity to Bdi-miR390, Bdi-miR2118, Bdi-2275 and the 22-nucleotide small RNAs listed in Additional file 1: Table S2. Each mismatch between the target and miRNA had a penalty of 1 and each G:U pairing had a penalty of 0.5. Those with a score of 5 or less were considered targets of the respective miRNA.

Accession numbers

Newly generated small RNA data have been deposited at the National Center for Biotechnology Information Gene Expression Omnibus (NCBI GEO) under accession number GSE52441 (GSM1266839 to GSM1266857). Two panicle small RNA libraries were previously published [31] and registered under GEO numbers GSM506620 and GSM506621. The data are also available online [101, 102].

Abbreviations

AGO:

ARGONAUTE protein

bp:

base pair

CDS:

coding sequence

DCL:

Dicer-like protein

D-plot:

degradation plot

kb:

kilobase

HD-ZIP:

homeobox domain-leucine zipper

LRR:

leucine-rich repeat

miRNA:

microRNA

MS medium:

Murashige and Skoog medium

NB:

nucleotide binding

PAGE:

polyacrylamide gel electrophoresis

PARE:

parallel analysis of RNA ends

PCR:

polymerase chain reaction

phasiRNA:

phased siRNA

RISC:

RNA-induced silencing complex

RT:

reverse transcriptase

siRNA:

small interfering RNA

tasiRNA:

trans-acting siRNA

TP10M:

transcripts per 10 million

TP2M:

transcripts per 2 million

WT:

wild type.

References

  1. Draper J, Mur LA, Jenkins G, Ghosh-Biswas GC, Bablak P, Hasterok R, Routledge AP: Brachypodium distachyon. A new model system for functional genomics in grasses. Plant Physiol. 2001, 127: 1539-1555. 10.1104/pp.010196.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  2. Mur LA, Allainguillaume J, Catalán P, Hasterok R, Jenkins G, Lesniewska K, Thomas I, Vogel J: Exploiting the Brachypodium tool box in cereal and grass research. New Phytol. 2011, 191: 334-347. 10.1111/j.1469-8137.2011.03748.x.

    Article  PubMed  Google Scholar 

  3. Christiansen P, Andersen CH, Didion T, Folling M, Nielsen KK: A rapid and efficient transformation protocol for the grass Brachypodium distachyon. Plant Cell Rep. 2005, 23: 751-758. 10.1007/s00299-004-0889-5.

    Article  PubMed  CAS  Google Scholar 

  4. Vain P, Worland B, Thole V, McKenzie N, Alves SC, Opanowicz M, Fish LJ, Bevan MW, Snape JW: Agrobacterium-mediated transformation of the temperate grass Brachypodium distachyon (genotype Bd21) for T-DNA insertional mutagenesis. Plant Biotechnol J. 2008, 6: 236-245. 10.1111/j.1467-7652.2007.00308.x.

    Article  PubMed  CAS  Google Scholar 

  5. Vain P, Thole V, Worland B, Opanowicz M, Bush MS, Doonan JH: A T-DNA mutation in the RNA helicase eIF4A confers a dose-dependent dwarfing phenotype in Brachypodium distachyon. Plant J. 2011, 66: 929-940. 10.1111/j.1365-313X.2011.04555.x.

    Article  PubMed  CAS  Google Scholar 

  6. Thole V, Peraldi A, Worland B, Nicholson P, Doonan JH, Vain P: T-DNA mutagenesis in Brachypodium distachyon. J Exp Bot. 2012, 63: 567-576. 10.1093/jxb/err333.

    Article  PubMed  CAS  Google Scholar 

  7. Kellogg EA: Evolutionary history of the grasses. Plant Physiol. 2001, 125: 1198-1205. 10.1104/pp.125.3.1198.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  8. Gomez LD, Bristow JK, Statham ER, McQueen-Mason SJ: Analysis of saccharification in Brachypodium distachyon stems under mild conditions of hydrolysis. Biotechnol Biofuels. 2008, 1: 15-10.1186/1754-6834-1-15.

    Article  PubMed  PubMed Central  Google Scholar 

  9. Chen X: Small RNAs in development – insights from plants. Curr Opin Genet Dev. 2012, 22: 361-367. 10.1016/j.gde.2012.04.004.

    Article  PubMed  PubMed Central  Google Scholar 

  10. Jin H: Endogenous small RNAs and antibacterial immunity in plants. FEBS Lett. 2008, 582: 2679-2684. 10.1016/j.febslet.2008.06.053.

    Article  PubMed  CAS  Google Scholar 

  11. Sunkar R, Li YF, Jagadeeswaran G: Functions of microRNAs in plant stress responses. Trends Plant Sci. 2012, 17: 196-203. 10.1016/j.tplants.2012.01.010.

    Article  PubMed  CAS  Google Scholar 

  12. Bartel DP: MicroRNAs: target recognition and regulatory functions. Cell. 2009, 136: 215-233. 10.1016/j.cell.2009.01.002.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  13. Voinnet O: Origin, biogenesis, and activity of plant microRNAs. Cell. 2009, 136: 669-687. 10.1016/j.cell.2009.01.046.

    Article  PubMed  CAS  Google Scholar 

  14. Brodersen P, Sakvarelidze-Achard L, Bruun-Rasmussen M, Dunoyer P, Yamamoto YY, Sieburth L, Voinnet O: Widespread translational inhibition by plant miRNAs and siRNAs. Science. 2008, 320: 1185-1190. 10.1126/science.1159151.

    Article  PubMed  CAS  Google Scholar 

  15. Lanet E, Delannoy E, Sormani R, Floris M, Brodersen P, Crété P, Voinnet O, Robaglia C: Biochemical evidence for translational repression by Arabidopsis microRNAs. Plant Cell. 2009, 21: 1762-1768. 10.1105/tpc.108.063412.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  16. Dugas DV, Bartel B: Sucrose induction of Arabidopsis miR398 represses two Cu/Zn superoxide dismutases. Plant Mol Biol. 2008, 67: 403-417. 10.1007/s11103-008-9329-1.

    Article  PubMed  CAS  Google Scholar 

  17. Yu B, Wang H: Translational inhibition by microRNAs in plants. Prog Mol Subcell Biol. 2010, 50: 41-57. 10.1007/978-3-642-03103-8_3.

    Article  PubMed  CAS  Google Scholar 

  18. Wu L, Zhou H, Zhang Q, Zhang J, Ni F, Liu C, Qi Y: DNA methylation mediated by a microRNA pathway. Mol Cell. 2010, 14: 465-475.

    Article  Google Scholar 

  19. Allen E, Xie Z, Gustafson AM, Carrington JC: microRNA-directed phasing during trans-acting siRNA biogenesis in plants. Cell. 2005, 121: 207-221. 10.1016/j.cell.2005.04.004.

    Article  PubMed  CAS  Google Scholar 

  20. Yoshikawa M, Peragine A, Park MY, Poethig RS: A pathway for the biogenesis of trans-acting siRNAs in Arabidopsis. Genes Dev. 2005, 19: 2164-2175. 10.1101/gad.1352605.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  21. Zhai J, Jeong DH, De Paoli E, Park S, Rosen BD, Li Y, González AJ, Yan Z, Kitto SL, Grusak MA, Jackson SA, Stacey G, Cook DR, Green PJ, Sherrier DJ, Meyers BC: MicroRNAs as master regulators of the plant NB-LRR defense gene family via the production of phased, trans-acting siRNAs. Genes Dev. 2011, 25: 2540-2553. 10.1101/gad.177527.111.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  22. Shivaprasad PV, Chen HM, Patel K, Bond DM, Santos BA, Baulcombe DC: A microRNA superfamily regulates nucleotide binding site-leucine-rich repeats and other mRNAs. Plant Cell. 2012, 24: 859-874. 10.1105/tpc.111.095380.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  23. Howell MD, Fahlgren N, Chapman EJ, Cumbie JS, Sullivan CM, Givan SA, Kasschau KD, Carrington JC: Genome-wide analysis of the RNA-DEPENDENT RNA POLYMERASE6/DICER-LIKE4 pathway in Arabidopsis reveals dependency on miRNA- and tasiRNA-directed targeting. Plant Cell. 2007, 19: 926-942. 10.1105/tpc.107.050062.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  24. Li F, Pignatta D, Bendix C, Brunkard JO, Cohn MM, Tung J, Sun H, Kumar P, Baker B: MicroRNA regulation of plant innate immune receptors. Proc Natl Acad Sci USA. 2012, 109: 1790-1795. 10.1073/pnas.1118282109.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  25. Axtell MJ, Jan C, Rajagopalan R, Bartel DP: A two-hit trigger for siRNA biogenesis in plants. Cell. 2006, 127: 565-577. 10.1016/j.cell.2006.09.032.

    Article  PubMed  CAS  Google Scholar 

  26. Cuperus JT, Fahlgren N, Carrington JC: Evolution and functional diversification of miRNA genes. Plant Cell. 2011, 23: 431-442. 10.1105/tpc.110.082784.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  27. Chen HM, Chen LT, Patel K, Li YH, Baulcombe DC, Wu SH: 22-Nucleotide RNAs trigger secondary siRNA biogenesis in plants. Proc Natl Acad Sci USA. 2010, 107: 15269-15274. 10.1073/pnas.1001738107.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  28. Cuperus JT, Carbonell A, Fahlgren N, Garcia-Ruiz H, Burke RT, Takeda A, Sullivan CM, Gilbert SD, Montgomery TA, Carrington JC: Unique functionality of 22-nt miRNAs in triggering RDR6-dependent siRNA biogenesis from target transcripts in Arabidopsis. Nat Struct Mol Biol. 2010, 17: 997-1003. 10.1038/nsmb.1866.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  29. Manavella PA, Koenig D, Weigel D: Plant secondary siRNA production determined by microRNA-duplex structure. Proc Natl Acad Sci USA. 2012, 109: 2461-2466. 10.1073/pnas.1200169109.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  30. Johnson C, Kasprzewska A, Tennessen K, Fernandes J, Nan GL, Walbot V, Sundaresan V, Vance V, Bowman LH: Clusters and superclusters of phased small RNAs in the developing inflorescence of rice. Genome Res. 2009, 19: 1429-1440. 10.1101/gr.089854.108.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  31. The International Brachypodium Initiative: Genome sequencing and analysis of the model grass Brachypodium distachyon. Nature. 2010, 463: 763-768. 10.1038/nature08747.

    Article  Google Scholar 

  32. Song X, Li P, Zhai J, Zhou M, Ma L, Liu B, Jeong DH, Nakano M, Cao S, Liu C, Chu C, Wang XJ, Green PJ, Meyers BC, Cao X: Roles of DCL4 and DCL3b in rice phased small RNA biogenesis. Plant J. 2012, 69: 462-474. 10.1111/j.1365-313X.2011.04805.x.

    Article  PubMed  CAS  Google Scholar 

  33. Meyers BC, Souret FF, Lu C, Green PJ: Sweating the small stuff: microRNA discovery in plants. Curr Opin Biotechnol. 2006, 17: 139-146. 10.1016/j.copbio.2006.01.008.

    Article  PubMed  CAS  Google Scholar 

  34. Jeong DH, Park S, Zhai J, Gurazada SG, De Paoli E, Meyers BC, Green PJ: Massive analysis of rice small RNAs: mechanistic implications of regulated microRNAs and variants for differential target RNA cleavage. Plant Cell. 2011, 23: 4185-4207. 10.1105/tpc.111.089045.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  35. Friedlander MR, Chen W, Adamidi C, Maaskola J, Einspanier R, Knespel S, Rajewsky N: Discovering microRNAs from deep sequencing data using miRDeep. Nat Biotechnol. 2008, 26: 407-415. 10.1038/nbt1394.

    Article  PubMed  Google Scholar 

  36. Meyers BC, Axtell MJ, Bartel B, Bartel DP, Baulcombe D, Bowman JL, Cao X, Carrington JC, Chen X, Green PJ, Griffiths-Jones S, Jacobsen SE, Mallory AC, Martienssen RA, Poethig RS, Qi Y, Vaucheret H, Voinnet O, Watanabe Y, Weigel D, Zhu JK: Criteria for annotation of plant MicroRNAs. Plant Cell. 2008, 20: 3186-3190. 10.1105/tpc.108.064311.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  37. Kozomara A, Griffiths-Jones S: miRBase: integrating microRNA annotation and deep-sequencing data. Nucleic Acids Res. 2011, 39: D152-D157. 10.1093/nar/gkq1027.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

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

    Article  PubMed  CAS  Google Scholar 

  39. Wei B, Cai T, Zhang R, Li A, Huo N, Li S, Gu YQ, Vogel J, Jia J, Qi Y, Mao L: Novel microRNAs uncovered by deep sequencing of small RNA transcriptomes in bread wheat (Triticum aestivum L.) and Brachypodium distachyon (L.) Beauv. Funct Integr Genomics. 2009, 9: 499-511. 10.1007/s10142-009-0128-9.

    Article  PubMed  CAS  Google Scholar 

  40. Zhang J, Xu Y, Huan Q, Chong K: Deep sequencing of Brachypodium small RNAs at the global genome level identifies microRNAs involved in cold stress response. BMC Genomics. 2009, 10: 449-10.1186/1471-2164-10-449.

    Article  PubMed  PubMed Central  Google Scholar 

  41. Baev V, Milev I, Naydenov M, Apostolova E, Minkov G, Minkov I, Yahubyan G: Implementation of a de novo genome-wide computational approach for updating Brachypodium miRNAs. Genomics. 2011, 97: 282-293. 10.1016/j.ygeno.2011.02.008.

    Article  PubMed  CAS  Google Scholar 

  42. Bertolini E, Verelst W, Horner DS, Gianfranceschi L, Piccolo V, Inze D, Pe ME, Mica E: Addressing the role of microRNAs in reprogramming leaf growth during drought stress in Brachypodium distachyon. Mol Plant. 2013, 6: 423-443. 10.1093/mp/sss160.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  43. Addo-Quaye C, Eshoo TW, Bartel DP, Axtell MJ: Endogenous siRNA and miRNA targets identified by sequencing of the Arabidopsis degradome. Curr Biol. 2008, 18: 758-762. 10.1016/j.cub.2008.04.042.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  44. German MA, Pillay M, Jeong DH, Hetawal A, Luo S, Janardhanan P, Kannan V, Rymarquis LA, Nobuta K, German R, De Paoli E, Lu C, Schroth G, Meyers BC, Green PJ: Global identification of microRNA-target RNA pairs by parallel analysis of RNA ends. Nat Biotechnol. 2008, 26: 941-946. 10.1038/nbt1417.

    Article  PubMed  CAS  Google Scholar 

  45. Gregory BD, O’Malley RC, Lister R, Urich MA, Tonti-Filippini J, Chen H, Millar AH, Ecker JR: A link between RNA metabolism and silencing affecting Arabidopsis development. Dev Cell. 2008, 14: 854-866. 10.1016/j.devcel.2008.04.005.

    Article  PubMed  CAS  Google Scholar 

  46. Li YF, Zheng Y, Addo-Quaye C, Zhang L, Saini A, Jagadeeswaran G, Axtell MJ, Zhang W, Sunkar R: Transcriptome-wide identification of microRNA targets in rice. Plant J. 2010, 62: 742-759. 10.1111/j.1365-313X.2010.04187.x.

    Article  PubMed  CAS  Google Scholar 

  47. Zhou M, Gu L, Li P, Song X, Chen Z, Wei L, Cao X: Degradome sequencing reveals endogenous small RNA targets in rice (Oryza sativa L. ssp. indica). Frontiers Biol. 2010, 5: 67-90. 10.1007/s11515-010-0007-8.

    Article  CAS  Google Scholar 

  48. Addo-Quaye C, Snyder JA, Park YB, Li YF, Sunkar R, Axtell MJ: Sliced microRNA targets and precise loop-first processing of MIR319 hairpins revealed by analysis of the Physcomitrella patens degradome. RNA. 2009, 15: 2112-2121. 10.1261/rna.1774909.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  49. Song QX, Liu YF, Hu XY, Zhang WK, Ma B, Chen SY, Zhang JS: Identification of miRNAs and their target genes in developing soybean seeds by deep sequencing. BMC Plant Biol. 2011, 11: 5-10.1186/1471-2229-11-5.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  50. Pantaleo V, Szittya G, Moxon S, Miozzi L, Moulton V, Dalmay T, Burgyan J: Identification of grapevine microRNAs and their targets using high-throughput sequencing and degradome analysis. Plant J. 2010, 62: 960-976.

    PubMed  CAS  Google Scholar 

  51. Zhang JZ, Ai XY, Guo WW, Peng SA, Deng XX, Hu CG: Identification of miRNAs and their target genes using deep sequencing and degradome analysis in trifoliate orange [Poncirus trifoliata L. Raf] [corrected]. Mol Biotechnol. 2012, 51: 44-57. 10.1007/s12033-011-9439-x.

    Article  PubMed  CAS  Google Scholar 

  52. Curaba J, Spriggs A, Taylor J, Li Z, Helliwell C: miRNA regulation in the early development of barley seed. BMC Plant Biol. 2012, 12: 120-10.1186/1471-2229-12-120.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  53. Zhang J, Mao Z, Chong K: A global profiling of uncapped mRNAs under cold stress reveals specific decay patterns and endonucleolytic cleavages in Brachypodium distachyon. Genome Biol. 2013, 14: R92-10.1186/gb-2013-14-8-r92.

    Article  PubMed  PubMed Central  Google Scholar 

  54. Axtell MJ, Westholm JO, Lai EC: Vive la difference: biogenesis and evolution of microRNAs in plants and animals. Genome Biol. 2011, 12: 221-10.1186/gb-2011-12-4-221.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  55. Merchan F, Boualem A, Crespi M, Frugier F: Plant polycistronic precursors containing non-homologous microRNAs target transcripts encoding functionally related proteins. Genome Biol. 2009, 10: R136-10.1186/gb-2009-10-12-r136.

    Article  PubMed  PubMed Central  Google Scholar 

  56. Fujii H, Chiou TJ, Lin SI, Aung K, Zhu JK: A miRNA involved in phosphate-starvation response in Arabidopsis. Curr Biol. 2005, 15: 2038-2043. 10.1016/j.cub.2005.10.016.

    Article  PubMed  CAS  Google Scholar 

  57. Aung K, Lin S-I, Wu C-C, Huang Y-T, Su C-l, Chiou T-J: pho2, a phosphate overaccumulator, is caused by a nonsense mutation in a microRNA399 target gene. Plant Physiol. 2006, 141: 1000-1011. 10.1104/pp.106.078063.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  58. Bari R, Datt Pant B, Stitt M, Scheible W-R: PHO2, MicroRNA399, and PHR1 define a phosphate-signaling pathway in plants. Plant Physiol. 2006, 141: 988-999. 10.1104/pp.106.079707.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  59. Chiou T-J, Aung K, Lin S-I, Wu C-C, Chiang S-F, Su C-L: Regulation of phosphate homeostasis by microRNA in Arabidopsis. Plant Cell. 2006, 18: 412-421. 10.1105/tpc.105.038943.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  60. Hsieh L-C, Lin S-I, Shih AC-C, Chen J-W, Lin W-Y, Tseng C-Y, Li W-H, Chiou T-J: Uncovering small RNA-mediated responses to phosphate deficiency in Arabidopsis by deep sequencing 10.1104/pp. 109.147280. Plant Physiol. 2009, 151: 2120-2132. 10.1104/pp.109.147280.

    Article  PubMed  PubMed Central  Google Scholar 

  61. Wang C, Ying S, Huang H, Li K, Wu P, Shou H: Involvement of OsSPX1 in phosphate homeostasis in rice. Plant J. 2009, 57: 895-904. 10.1111/j.1365-313X.2008.03734.x.

    Article  PubMed  CAS  Google Scholar 

  62. Hu B, Zhu C, Li F, Tang J, Wang Y, Lin A, Liu L, Che R, Chu C: LEAF TIP NECROSIS1 plays a pivotal role in the regulation of multiple phosphate starvation responses in rice. Plant Physiol. 2011, 156: 1101-1115. 10.1104/pp.110.170209.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  63. Souret FF, Kastenmayer JP, Green PJ: AtXRN4 degrades mRNA in Arabidopsis and its substrates include selected miRNA targets. Mol Cell. 2004, 15: 173-183. 10.1016/j.molcel.2004.06.006.

    Article  PubMed  CAS  Google Scholar 

  64. Rymarquis LA, Souret FF, Green PJ: Evidence that XRN4, an Arabidopsis homolog of exoribonuclease XRN1, preferentially impacts transcripts with certain sequences or in particular functional categories. RNA. 2011, 17: 501-511. 10.1261/rna.2467911.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  65. Xia R, Zhu H, An YQ, Beers EP, Liu Z: Apple miRNAs and tasiRNAs with novel regulatory networks. Genome Biol. 2012, 13: R47-10.1186/gb-2012-13-6-r47.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  66. Zhu H, Xia R, Zhao B, An YQ, Dardick CD, Callahan AM, Liu Z: Unique expression, processing regulation, and regulatory network of peach (Prunus persica) miRNAs. BMC Plant Biol. 2012, 12: 149-10.1186/1471-2229-12-149.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  67. Fei Q, Xia R, Meyers BC: Phased, secondary, small interfering RNAs in posttranscriptional regulatory networks. Plant Cell. 2013, 25: 2400-2415. 10.1105/tpc.113.114652.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  68. Jiao Y, Wang Y, Xue D, Wang J, Yan M, Liu G, Dong G, Zeng D, Lu Z, Zhu X, Quian Q, Li J: Regulation of OsSPL14 by OsmiR156 defines ideal plant architecture in rice. Nat Genet. 2010, 42: 541-544. 10.1038/ng.591.

    Article  PubMed  CAS  Google Scholar 

  69. Miura K, Ikeda M, Matsubara A, Song XJ, Ito M, Asano K, Matsuoka M, Kitano H, Ashikari M: OsSPL14 promotes panicle branching and higher grain productivity in rice. Nat Genet. 2010, 42: 545-549. 10.1038/ng.592.

    Article  PubMed  CAS  Google Scholar 

  70. Xie K, Wu C, Xiong L: Genomic organization, differential expression, and interaction of SQUAMOSA promoter-binding-like transcription factors and microRNA156 in rice. Plant Physiol. 2006, 142: 280-293. 10.1104/pp.106.084475.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  71. Xie Z, Kasschau KD, Carrington JC: Negative feedback regulation of dicer-like1 in Arabidopsis by microRNA-guided mRNA degradation. Curr Biol. 2003, 13: 784-789. 10.1016/S0960-9822(03)00281-1.

    Article  PubMed  CAS  Google Scholar 

  72. Axtell MJ, Bowman JL: Evolution of plant microRNAs and their targets. Trends Plant Sci. 2008, 13: 343-349. 10.1016/j.tplants.2008.03.009.

    Article  PubMed  CAS  Google Scholar 

  73. Rajagopalan R, Vaucheret H, Trejo J, Bartel DP: A diverse and evolutionarily fluid set of microRNAs in Arabidopsis thaliana. Genes Dev. 2006, 20: 3407-3425. 10.1101/gad.1476406.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  74. Yao Y, Guo G, Ni Z, Sunkar R, Du J, Zhu JK, Sun Q: Cloning and characterization of microRNAs from wheat (Triticum aestivum L.). Genome Biol. 2007, 8: R96-10.1186/gb-2007-8-6-r96.

    Article  PubMed  PubMed Central  Google Scholar 

  75. Schreiber AW, Shi BJ, Huang CY, Langridge P, Baumann U: Discovery of barley miRNAs through deep sequencing of short reads. BMC Genomics. 2011, 12: 129-10.1186/1471-2164-12-129.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  76. Zhang L, Zheng Y, Jagadeeswaran G, Li Y, Gowdu K, Sunkar R: Identification and temporal expression analysis of conserved and novel microRNAs in sorghum. Genomics. 2011, 98: 460-468. 10.1016/j.ygeno.2011.08.005.

    Article  PubMed  CAS  Google Scholar 

  77. Duan R, Pak C, Jin P: Single nucleotide polymorphism associated with mature miR-125a alters the processing of pri-miRNA. Hum Mol Genet. 2007, 16: 1124-1131. 10.1093/hmg/ddm062.

    Article  PubMed  CAS  Google Scholar 

  78. Sun G, Yan J, Noltner K, Feng J, Li H, Sarkis DA, Sommer SS, Rossi JJ: SNPs in human miRNA genes affect biogenesis and function. RNA. 2009, 15: 1640-1651. 10.1261/rna.1560209.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  79. Ryan BM, Robles AI, Harris CC: Genetic variation in microRNA networks: the implications for cancer research. Nat Rev Cancer. 2010, 10: 389-402. 10.1038/nrc2867.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  80. Todesco M, Balasubramanian S, Cao J, Ott F, Sureshkumar S, Schneeberger K, Meyer RC, Altmann T, Weigel D: Natural variation in biogenesis efficiency of individual Arabidopsis thaliana microRNAs. Curr Biol. 2012, 22: 166-170. 10.1016/j.cub.2011.11.060.

    Article  PubMed  CAS  Google Scholar 

  81. Pena-Castro JM, van Zanten M, Lee SC, Patel MR, Voesenek LA, Fukao T, Bailey-Serres J: Expression of rice SUB1A and SUB1C transcription factors in Arabidopsis uncovers flowering inhibition as a submergence tolerance mechanism. Plant J. 2011, 67: 434-446. 10.1111/j.1365-313X.2011.04605.x.

    Article  PubMed  CAS  Google Scholar 

  82. Chuck G, Cigan AM, Saeteurn K, Hake S: The heterochronic maize mutant Corngrass1 results from overexpression of a tandem microRNA. Nat Genet. 2007, 39: 544-549. 10.1038/ng2001.

    Article  PubMed  CAS  Google Scholar 

  83. Chuck GS, Tobias C, Sun L, Kraemer F, Li C, Dibble D, Arora R, Bragg JN, Vogel JP, Singh S, Simmons BA, Pauly M, Hake S: Overexpression of the maize Corngrass1 microRNA prevents flowering, improves digestibility, and increases starch content of switchgrass. Proc Natl Acad Sci USA. 2011, 108: 17550-17555. 10.1073/pnas.1113971108.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  84. Fu C, Sunkar R, Zhou C, Shen H, Zhang JY, Matts J, Wolf J, Mann DG, Stewart CN, Tang Y, Wang ZY: Overexpression of miR156 in switchgrass (Panicum virgatum L.) results in various morphological alterations and leads to improved biomass production. Plant Biotechnol J. 2012, 10: 443-452. 10.1111/j.1467-7652.2011.00677.x.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  85. Schwab R, Palatnik JF, Riester M, Schommer C, Schmid M, Weigel D: Specific effects of microRNAs on the plant transcriptome. Dev Cell. 2005, 8: 517-527. 10.1016/j.devcel.2005.01.018.

    Article  PubMed  CAS  Google Scholar 

  86. Lu C, Meyers BC, Green PJ: Construction of small RNA cDNA libraries for deep sequencing. Methods. 2007, 43: 110-117. 10.1016/j.ymeth.2007.05.002.

    Article  PubMed  Google Scholar 

  87. Brachypodium genome 8X assembly. [ftp://ftp.jgi-psf.org/pub/compgen/phytozome/v9.0/Bdistachyon/]

  88. Langmead B, Trapnell C, Pop M, Salzberg SL: Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol. 2009, 10: R25-10.1186/gb-2009-10-3-r25.

    Article  PubMed  PubMed Central  Google Scholar 

  89. Markham NR, Zuker M: UNAFold: software for nucleic acid folding and hybridization. Methods Mol Biol. 2008, 453: 3-31. 10.1007/978-1-60327-429-6_1.

    Article  PubMed  CAS  Google Scholar 

  90. Maroney PA, Chamnongpol S, Souret F, Nilsen TW: A rapid, quantitative assay for direct detection of microRNAs and other small RNAs using splinted ligation. RNA. 2007, 13: 930-936. 10.1261/rna.518107.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  91. Maroney PA, Chamnongpol S, Souret F, Nilsen TW: Direct detection of small RNAs using splinted ligation. Nat Protoc. 2008, 3: 279-287. 10.1038/nprot.2007.530.

    Article  PubMed  CAS  Google Scholar 

  92. Jeong DH, Green PJ: Methods for validation of miRNA sequence variants and the cleavage of their targets. Methods. 2012, 58: 135-143. 10.1016/j.ymeth.2012.08.005.

    Article  PubMed  CAS  Google Scholar 

  93. Church G, Gilbert W: Genomic sequencing. Proc Natl Acad Sci USA. 1984, 81: 1991-1995. 10.1073/pnas.81.7.1991.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  94. German MA, Luo S, Schroth G, Meyers BC, Green PJ: Construction of parallel analysis of RNA ends (PARE) libraries for the study of cleaved miRNA targets and the RNA degradome. Nat Protocols. 2009, 4: 356-362. 10.1038/nprot.2009.8.

    Article  PubMed  CAS  Google Scholar 

  95. CleaveLand v2.0. [http://axtell-lab-psu.weebly.com/cleaveland.html]

  96. Brachypodium coding sequence CDSv1.2. [http://compbio.cs.toronto.edu/shrimp/]

  97. Addo-Quaye C, Miller W, Axtell MJ: CleaveLand: a pipeline for using degradome data to find cleaved small RNA targets. Bioinformatics. 2009, 25: 130-131. 10.1093/bioinformatics/btn604.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  98. Allen E, Xie Z, Gustafson AM, Sung GH, Spatafora JW, Carrington JC: Evolution of microRNA genes by inverted duplication of target gene sequences in Arabidopsis thaliana. Nat Genet. 2004, 36: 1282-1290. 10.1038/ng1478.

    Article  PubMed  CAS  Google Scholar 

  99. psRNATarget. [http://plantgrn.noble.org/psRNATarget/]

  100. Dai X, Zhao PX: psRNATarget: a plant small RNA target analysis server. Nucleic Acids Res. 2011, 39: W155-W159. 10.1093/nar/gkr319.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  101. Brachypodium small RNA data. [http://mpss.udel.edu/brachy_sRNA2]

  102. Brachypodium PARE data. [http://mpss.udel.edu/brachy_pare2]

Download references

Acknowledgements

We thank Yang Yen (South Dakota State University) for helping to establish the conditions for plant growth, Emanuele De Paoli for helpful discussions, Bruce Kingham and Greg May (National Center for Genome Resources) for assistance with Illumina sequencing, and Mayumi Nakano for her work on the web interface. This research was supported primarily by the Office of Biological and Environmental Research of the Department of Energy under Plant Feedstock Genomics for Bioenergy award number DE-FG02-07ER64450 to PJG, with additional support from the National Science Foundation, award number MCB-1231726 and the National Institutes of Health, award number AI043288 to JCC. NF was supported, in part, by a USDA-NIFA postdoctoral fellowship.

Author information

Authors and Affiliations

Authors

Corresponding author

Correspondence to Pamela J Green.

Additional information

Competing interests

The authors declare that they have no competing interests.

Authors’ contributions

D-HJ, SAS, SP, MG, MAG and MA performed the experiments. D-HJ, SAS, LAR, JZ, BCM and PJG analyzed the data. DFG provided materials. NF, SEF, TCM and JCC developed and provided data for the OBD libraries. D-HJ, SAS, LAR and PJG wrote the manuscript. All authors read and approved the final manuscript.

Electronic supplementary material

13059_2013_3327_MOESM1_ESM.docx

Additional file 1: Table S1: Precursors of Brachypodium miRNAs. Table S2. Conserved and non-conserved annotated miRNAs. Table S3. miRNA targets with PARE sequences at the predicted cleavage sites. Table S4. Arabidopsis miRNA targets. Table S5. Oligomers used in this study. (DOCX 227 KB)

13059_2013_3327_MOESM2_ESM.pptx

Additional file 2: Figure S1: Predicted secondary structure of Bdi-MIR444 precursors. Figure S2. Target prediction programs identified distinct subsets of targets. Figure S3. Representative D-plots for miRNA targets at each stringency level. Figure S4. Bdi-miR2118 cleavage induces phasing. Figure S5. Bdi-miR5200 in diverse plant species. Figure S6. D-plots of AGO1 for Arabidopsis, rice and Brachypodium.(PPTX 817 KB)

Authors’ original submitted files for images

Rights and permissions

Reprints and permissions

About this article

Cite this article

Jeong, DH., Schmidt, S.A., Rymarquis, L.A. et al. Parallel analysis of RNA ends enhances global investigation of microRNAs and target RNAs of Brachypodium distachyon. Genome Biol 14, R145 (2013). https://doi.org/10.1186/gb-2013-14-12-r145

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/gb-2013-14-12-r145

Keywords