A diverse and evolutionarily fluid set of microRNAs in Arabidopsis thaliana

The most abundant reads corresponded to conserved, previously identified miRNAs

As expected, known miRNAs were the sequences most redundantly retrieved from the pool, boasting the highest read frequency of all small RNA classes and 15% of the total (Table 1). All of the miRNA families known to be conserved to poplar and rice (20 families) or just poplar (one additional family; Jones-Rhoades et al. 2006) were represented among our reads, with frequencies as high as 36,093 (miR167). Even the stress-inducible miRNAs miR395 and miR399, previously undetectable in plants grown under normal conditions (Jones-Rhoades and Bartel 2004; Fujii et al. 2005), were represented (13 and 580 reads, respectively), suggesting that some other miRNAs induced in specific conditions might also be represented by multiple reads in our data set. For a few of the miRNAs, including miR319a/b, the sequenced species differed slightly from the annotated species, suggesting refinements of the annotated species (Supplementary Database 1).

Most previously identified conserved miRNA families have multiple, paralogous loci, which combined total 92 loci (Jones-Rhoades et al. 2006). In some cases, members of the same family have slightly different sequences, which can sometimes target distinct sets of messages (Schwab et al. 2005). In other cases, paralogous loci give rise to identical mature miRNAs, raising the question of which paralogs are expressed. Mapping the miRNA* species, rare side products, or degradation fragments unique to a single paralog enabled us to confirm the expression of all but 13 of the 92 loci (Supplementary Database 1), including 19 whose transcription had not been previously confirmed, either by cloning or by mapping the 5′ end of primary transcripts (Xie et al. 2005a). The exceptions were for loci for which no reads could be uniquely mapped.

For most miRNAs, variants of the most abundant read were isolated with 5′ or 3′ heterogeneity, evidenced by missing or extra nucleotides at each terminus (Supplementary Database 1). Occasional slippage of DCL1 processing presumably gives rise to the extra bases, whereas missing nucleotides could result from slippage or end degradation. In contrast to metazoan miRNAs (Lau et al. 2001), we found that 5′ heterogeneity was common in Arabidopsis miRNA pools and only slightly less prevalent than 3′ heterogeneity, with no correlation between the extent of heterogeneity and the arm of the foldback that produced the miRNA (Supplementary Database 1). Plant miRNAs might tolerate more extensive 5′ heterogeneity because seed pairing represents a smaller portion of their targeting specificity (Mallory et al. 2004), whereas animal miRNAs truncated or extended by a single nucleotide at their 5′ end would no longer recognize many normal targets and would instead recognize many other messages (Lim et al. 2005).

The miRNA* species had ∼9% as many reads as the mature miRNAs, which was higher than the 1% observed in worms (Ruby et al. 2006). This percentage varied widely. For two of the 21 conserved families—miR395 and miR397, represented by 13 and 361 reads, respectively—no miRNA* species were observed (Supplementary Database 1). At the other extreme, miR403* was observed more frequently than mature miR403 (1643 and 66 reads, respectively). Mature miR403 directs cleavage of AGO2 mRNA (Allen et al. 2005) and is more easily detected by RNA blotting than is miR403* (H. Vaucheret, unpubl.). We infer that sequencing abundance does not always correlate with in vivo abundance, which in any event might not always predict the functional strand.

Several characteristics of Arabidopsis miRNAs and their foldbacks emerged from analysis of reads corresponding to miRNA loci that had previously been confidently identified. First, relatively few unique nonoverlapping reads mapped to authentic miRNA foldbacks (Supplementary Database 1). Although DCL1 processing on some foldbacks appeared a little sloppy, it appeared at least globally very precise in that most reads centered on the miRNA/miRNA*, even for stems that were quite extensive. Second, the miRNA* sequence was observed for most loci. Third, <1% of reads mapped to the strand antisense to that giving rise to miRNA and miRNA* (Supplementary Database 1).

Arabidopsis has many miRNAs lacking close orthologs in other sequenced plants

In addition to the 21 conserved miRNA families, another five apparently nonconserved miRNA genes (miR158, miR161, miR163, miR173, and miR447) have been confidently identified in Arabidopsis (Jones-Rhoades et al. 2006). Each was represented among our reads, with read frequency ranging from 29 for miR447 to 10,573 for miR161 (Supplementary Database 1). How many additional nonconserved miRNAs might exist in Arabidopsis? The multitude of other endogenous small RNAs, some of which derive from regions with fortuitous potential to fold into miRNA-like hairpins, has complicated miRNA identification in plants, leading to the suggestion that biogenetic requirements be confirmed using mutant backgrounds prior to annotation (Jones-Rhoades et al. 2006). High-throughput sequencing offered an alternative approach for distinguishing miRNAs from other small RNAs. Candidates from loci with a substantial number of reads deriving from the antisense strand can also be excluded because such antisense reads suggest origin from a perfect dsRNA rather than a hairpin. For the remaining candidates meeting the conventional hairpin-pairing criteria, a sequenced miRNA* species, especially one with 2-nt 3′ overhangs, provides strong evidence that a candidate originates from a DCL-processed stem-loop. As a result, demonstrating that the candidate accumulates in prescribed mutant backgrounds becomes less important, which is particularly helpful for miRNAs difficult to detect on blots.

Using these criteria, we identified 38 additional Arabidopsis miRNA families, thereby increasing by 2.5-fold the known diversity of miRNAs in Arabidopsis (Fig. 1A,B; Table 2; Supplementary Database 2; Table 2 and Supplementary Database 2 also include a 39th miRNA, miR391, which was absent in miRBase version 7.1). To receive miRNA designation, a miRNA* species (or close variant if the miRNA was sequenced at least three times) must have been observed among the reads, with the exception of miR823, which was validated using RNA blotting and the conventional set of mutants (below). The most abundant read on the foldback was deemed the miRNA, although in cases where read density was roughly equivalent for the most abundant reads from each arm of the foldback, both were together annotated to represent the new miRNA locus (using the 5′ and 3′ designations adopted in similar cases for metazoan miRNAs). Many other miRNAs might exist in Arabidopsis; another 40 candidates mapped to plausible hairpins but lacked reads representing the miRNA* species (Supplementary Table 4). Four particularly compelling candidates, each sequenced >25 times (Table 2, CandidateA–CandidateD), were carried forward in subsequent analysis, anticipating that they will eventually be validated.

Open in a separate window

Figure 1.

Newly identified Arabidopsis miRNAs. (A) Predicted secondary structures of miRNA hairpins highlighting the miRNA (red) and miRNA* species (blue). (B) The miRNA hairpins of A, shown in bracket notation with a tally of reads mapping to the hairpin and nucleotides colored as in A. (C) RNA blots demonstrating that accumulation of six detectable miRNAs depended on DCL1, not DCL2, DCL3, DCL4, RDR2, or RDR6. As a loading control, blots were stripped and reprobed for U6. (D) Sequencing frequencies of Arabidopsis miRNA families. Shown are cumulative plots for all Arabidopsis miRNA families (red squares), conserved families (violet diamonds), and all families plus the 40 sequenced candidates (gray triangles). Fourteen families, 11 of which were conserved, were sequenced at a frequency of greater than one per 1000 (dashed line).

Table 2.

Newly identified miRNAs in Arabidopsis

Open in a separate window

^aLength of mature miRNA sequence.

^bNumber of reads from seedlings (Se), rosette leaves (R), flowers (F), and siliques (Si).

^cThe sequencing of a miRNA* species is denoted “Yes” if the perfect match to the miRNA* (with the 2-nt 3′ overhangs typical of DCL products) was sequenced and was the most abundant read from that arm of the hairpin, “Yes*” if the perfect miRNA* was sequenced but was not the most abundant read from that arm of the hairpin, “Yes**” if only a close heterogeneous variant of the perfect miRNA* was sequenced, and “No” if no star species was recovered.

^dAGI codes are given for genes with top-scoring target sites for the miRNA. A complete list of predicted mRNA targets of newly identified miRNAs, along with the associated target site alignments and scores, is given in Supplementary Database 4.

^eProtein products of predicted targets in the best score class. The total number of predicted targets falling within the cutoff is given in parentheses if more than one target was predicted.

^fReported as a miRNA in Xie et al. (2005a) but was not annotated in miRBase version 7.

^gLocus was reported as a miRNA in Lu et al. (2006).

A search in plant expressed sequence tags (EST) data sets, and the Oryza sativa (rice) and Populus trichocarpa (poplar) genomes, revealed potential orthologs for only one of the newly identified miRNAs, miR828, which had recognizable orthologs in poplar and leafy spurge, each with one substitution in the mature miRNA (Supplementary Fig. 2). For all other newly identified miRNAs, potential orthologs were either absent in sequenced genomes or found only after relaxing the homology criterion to allow three point substitutions. However, most, if not all, of these candidates appeared to be false positives, because at this stringency an equivalent number of hits was found that satisfied the homology and pairing criteria but mapped to the nonhomologous arms of predicted hairpins. We concluded that most of the newly identified miRNAs do not have identifiable orthologs in the sequence databases and henceforth refer to them all as “nonconserved,” while recognizing that a few might have divergent orthologs difficult to identify with confidence, and that many might have orthologs in unsequenced species more closely related to Arabidopsis thaliana. Mirroring the search for orthologs, we found no convincing Arabidopsis paralogs of the newly identified miRNAs.

Although screening was performed on 20- to 24-nt reads, without preference for a particular length or 5′ nucleotide, 74% of the newly identified miRNA loci encoded a 21-nt miRNA, and 87% encoded a miRNA beginning with a U (Table 2). Thus, these characteristics of the conserved miRNAs (Reinhart et al. 2002) were shared by the newly identified miRNAs. Some intriguing tissue specificities were also evident (Table 2). For example, miR771 and miR839 were sequenced primarily from flowers, miR391 and miR825 appeared preferentially in rosette leaves, miR822 and miR842 were preferentially sequenced in seedlings, and miR828 was sequenced most often from siliques. For some with the most striking specificities, we speculate that expression might be at a high level within just a few specialized cells within that organ.

DCL4 processes some Arabidopsis miRNAs

Most of the newly identified miRNAs were infrequently recovered by deep sequencing, with median read frequencies of only 13, compared with 731 for the conserved families, suggesting that in plants nonconserved miRNAs are generally expressed at low levels or primarily in specific cells or growth conditions. When RNA from plants grown under normal laboratory conditions was blotted and probed for the 10 most abundant newly identified miRNAs, only eight could be detected using either DNA or LNA probes. Accumulation of six of these eight miRNAs displayed the classical biogenetic profile of DCL1-dependency, with insensitivity to defects in any of the other DCL enzymes or RDR proteins (Fig. 1C).

In contrast to the previously characterized Arabidopsis miRNAs, accumulation of two of the eight miRNAs detectable on RNA blots depended on DCL4, not DCL1 (Fig. 2A). One, miR822, was previously classified as an siRNA (ASRP1729) because it accumulates in dcl1 plants (Allen et al. 2004; Xie et al. 2004). Both miR822 and miR839 were insensitive to defects in RDR2 and RDR6, as expected for RNAs that derive from hairpins rather than dsRNA. Further supporting a hairpin precursor structure was the pattern of reads from these loci (Fig. 2B). Over 99% of reads arose from one strand, with only two of the 1892 MIR822 reads and one of the 332 MIR839 reads deriving from the antisense strand—a pattern inconsistent with a perfect dsRNA intermediate. Furthermore, the major species from each arm of the predicted foldbacks paired to each other, with 2-nt 3′ overhangs observed for the miR822:miR822* duplex. Although the cleavage precision did not appear to match that of DCL1, this preferred processing from a localized region of an RNA hairpin stem satisfied the defining feature of miRNAs. We concluded that transcripts from a few miRNA loci are processed by DCL4 rather than by DCL1. The dependency on DCL4 for miR822 and miR839 accumulation appeared even higher than that for tasiRNA accumulation; in the absence of DCL4, tasiRNA precursors are processed into 22-nt and 24-nt species by DCL2 and DCL3, respectively (Gasciolli et al. 2005; Xie et al. 2005b), whereas miR822 and miR839 species are not detectable in either dcl4-1 or dcl4-2 plants (Fig. 2A; data not shown).

Open in a separate window

Figure 2.

DCL4-dependent miRNAs in Arabidopsis. (A) RNA blots demonstrating that the accumulation of two detectable miRNAs depended on DCL4 and not on DCL1, DCL2, DCL3, RDR2, or RDR6. As a loading control, blots were stripped and reprobed for U6. Similar results were obtained with dcl4-1 and dcl4-2 alleles (data not shown). (B) Predicted secondary structures of the miRNA hairpins, with lines denoting the sequences mapping to the miRNA (top) and miRNA* (bottom) arm of each hairpin. The thickness and color of the lines correspond to the number of total reads representing each small RNA species, as indicated in the key. The two reads corresponding to the antisense of miR822 and the single read mapping to the antisense of miR839 are not depicted.

Predicted targets of newly identified miRNAs

Conserved miRNA targets can be predicted with very high confidence, whereas in single-genome analyses only the more extensively paired interactions can be predicted with reasonable confidence (Jones-Rhoades and Bartel 2004). To better predict nonconserved interactions, scoring rubrics have been developed that preferentially penalize mismatches to the 5′ and central regions of the miRNA (Allen et al. 2005; Schwab et al. 2005), which are more disruptive than those to the 3′ region of the miRNA (Mallory et al. 2004). When applying the rubric of Allen et al. (2005) in a single-genome search to predict targets of 22 unrelated miRNAs, scoring cutoffs that captured 86% of the experimentally confirmed targets of these miRNAs gave a ratio of authentic to false-positive predictions of 6.9:1, estimated by summing the number of targets predicted for the miRNAs and comparing with the average predicted for 10 shuffled cohorts. Using these score cutoffs, we applied the rubric to predict targets of the newly discovered miRNAs, achieving a lower, although still significant, estimated signal:noise ratio of 3.0:1 (Table 2; Supplementary Database 4).

One explanation for the apparently lower specificity was that for six miRNAs, the miRNA and miRNA* species were difficult to distinguish from each other, and thus both were included in the target prediction analysis, recognizing that one of the strands might contribute only false-positive predictions. Similarly, two register-shifted sequences of roughly equal abundance from the miR829 foldback were included. Another explanation might be that some of the newly identified miRNA families have fewer targets with extensive complementarity than do the previously identified families. Indeed, some might not have any biological targets, a subset of which might be “young” DCL1/DCL4 substrates whose processing will soon be lost in the course of neutral evolutionary drift unless a beneficial targeting interaction emerges first. Nonetheless, the prediction of three times as many targets as expected by chance suggested that many of the newly identified miRNAs down-regulate genes. Targets for three of the more abundant miRNAs were validated by 5′ RACE (Supplementary Fig. 3). These were CMT3, a miR823 target that encodes a CpNpG DNA cytosine methyltransferase; AGL16, a miR824 target that encodes a MADS-box transcription factor; and MYB113, a miR828 target that encodes a MYB transcription factor.

Predicted targets of the newly identified miRNAs included transcription factors in the MYB and AP2 families, which each have paralogs known to be targeted by previously identified miRNA families (Supplementary Database 4). In addition, members of transcription factor families not previously known to be regulated by Arabidopsis miRNAs, such as MADS-box, ERF (ethylene response factor), WHIRLY, and Dof (DNA binding with one finger) proteins, were among the predicted targets. F-box-containing proteins added to the list of known miRNA targets implicated in protein degradation. A PPR gene distinct from those known to be targeted by miR161, miR400, or TAS1 or TAS2 trans-acting siRNAs was also among the predictions. Other predictions extended the biological processes thought to be regulated by miRNAs. For example, nine jacalin lectins, predicted miR842 and miR846 targets, bind complex carbohydrates and are thought to be involved in initiating pathogen defense responses (Geshi and Brandt 1998). Histone variants, and epigenetic silencing machinery such as CMT3 and a bromo-adjacent homology (BAH) domain-containing protein were predicted targets, suggesting that Arabidopsis miRNAs regulate transcriptional silencing pathway components in addition to targeting miRNA biogenetic and effector proteins like DCL1 and AGO1.

Evolutionary origins of miRNA genes

Some miRNAs might have arisen from duplication of their target loci, and if so, those that were recently derived might exhibit similarity to their targets that extends beyond the mature miRNA sequence, as observed previously for miR161 and miR163 (Allen et al. 2004). Six of the newly identified miRNA loci displayed extended sequence similarity with their predicted target genes, diagnostic of common origins. Both arms of the MIR822 gene were previously observed to have an extended alignment to several DC1 domain-containing genes (Allen et al. 2005). The same pattern was seen for MIR841 and MIR826 and their predicted targets (Fig. 3A,B).

Open in a separate window

Figure 3.

Extended homology between miRNA genes and their predicted target genes, suggestive of common origin. (A) MIR841, which illustrates a pattern of extended homology (orange shading) resembling that observed previously for MIR161 and MIR163 and their respective targets (Allen et al. 2005). Segments corresponding to the mature miRNA (red) and miRNA* (blue) are indicated. The diagram (top) depicts the target gene in right–left polarity, and the alignment (bottom) depicts the target gene and miRNA* segment as their reverse complements (rc). Numbers indicate positions in the protein-coding gene (At2g38810), counting from its first annotated nucleotide. Nucleotides are shaded to indicate those shared by all (orange) or most (yellow) aligned sequences. (B) MIR826, for which extended homology suggested an evolutionary pathway whereby a later duplication creating the miRNA* arm was nested within an earlier duplication. The genomic proximity of the miRNA and target gene is shown in the top diagram. For the middle and bottom diagrams, drawing conventions are as in A. (C) MIR846, for which extended homology suggested tandem duplication within an ancestral gene whereby the duplicated regions independently gave rise to the miRNA and miRNA* segments.

A different pattern was observed for MIR842 and MIR846, suggesting an alternative pathway for miRNA gene emergence. As illustrated for MIR846, these genes appeared to derive from two regions of their predicted targets, rather than one (Fig. 3C). The simplest explanation for the dual alignment to their targets, with the miRNA arm of the hairpin aligning to one region of the target and the miRNA* arm aligning to the other region, was that a duplication within the targets preceded the duplications that gave rise to the miRNA locus.

Another interesting miRNA–target configuration involved MIR840, which was expressed from the opposite strand of its predicted target gene, AtWhirly3. This is an arrangement first observed for an Epstein-Barr Virus miRNA and its target (Pfeffer et al. 2004), but one that had not been seen in plants. AtWhirly3 encodes a homolog of potato p24, a known transcriptional regulator of plant defense and disease resistance genes. In the sense orientation, the miRNA was found within the annotated 3′ untranslated region (UTR) of a PPR mRNA, At2g02750. Although both strands encode a hairpin, our reads did not include any small RNA sequences from the AtWhirly3 strand. Either the presumptive miRNA or its star sequence could target the AtWhirly3 3′ UTR for cleavage. This implies a mechanism by which the expression of one member of a convergent gene pair influences the output of the other—a miRNA counterpart to that observed previously for a convergent gene pair that generates nat-siRNAs (Borsani et al. 2005).

Of the 44 genes for the miRNAs and candidates listed in Table 2, 35 were in regions between annotated genes, as is typical of plant miRNA genes (Reinhart et al. 2002), whereas nine overlapped protein-coding genes. One was miR840, described above. Of the remaining eight, miR837, miR838, miR848, miR852, and CandidateD overlapped introns, in the same orientation as the protein-coding host gene—an arrangement that bypasses the need to acquire an independent promoter (Baskerville and Bartel 2005). Mature miR837 also had a second match in the genome, located within the same intron that contains the miR837 stem-loop but in the antisense orientation, suggesting that miR837 might target the pre-mRNA of its host gene, an oligopeptide transporter. miR841 derived from the strand antisense to the intron of At4g13570, a gene closely related to one of its predicted targets but itself not predicted because our search was limited to spliced messages. miR777 and miR834 were localized to the 5′ UTR and 3′ UTR of genes, respectively, with their foldbacks potentially extending into annotated protein-coding regions.

A homeostatic self-regulatory mechanism for DCL1

miR838 derived from a hairpin within intron 14 of the DCL1 mRNA (Fig. 4A). The foldback potential of this intron was previously noted, and RACE mapping of the DCL1 transcript revealed a 4.0-kb fragment whose 3′ end terminates at the exon 14/15 junction, and a population of ∼2.5-kb fragments, some of which have 5′ ends falling within intron 14 (Xie et al. 2003). Because small RNAs were not detected, the fragments were attributed to aberrant splicing at intron 14 (Xie et al. 2003).

Open in a separate window

Figure 4.

An intronic hairpin positioned so as to mediate DCL1 autoregulation. (A) Intron 14 of the DCL1 primary transcript, and the predicted hairpin structure of miR838. Arrows indicate 5′ ends of RACE-mapped fragments. (B) Alternative fates of the DCL1 primary transcript, which appears to undergo either splicing to generate full-length DCL1 mRNA or processing by DCL1 itself to generate transcript fragments severed within intron 14. (C) A schematic of DCL1 post-transcriptional autoregulation. When DCL1 protein levels are high, it could compete with splicing machinery for access to intron 14, thereby supplementing miR162-mediated regulation to maintain the proper level of DCL1 mRNA.

We propose that the presence of this intronic miRNA enables a self-regulatory mechanism that helps maintain DCL1 homeostasis (Fig. 4B). When nuclear DCL1 protein levels are high, the miRNA biogenesis machinery (including DCL1 and HYL1) could compete more efficiently than the splicing machinery for the DCL1 precursor transcript. If DCL1 began to process the miRNA hairpin before the intron 14 splice sites were defined and juxtaposed during spliceosome formation, then DCL1 expression would shift toward a pool of truncated, nonfunctional DCL1 transcripts, thereby providing a regulatory feedback mechanism that supplements miR162-directed cleavage (Fig. 4C). 5′ RACE confirmed that a population of fragments had 5′ ends terminating at the ends of miR838 (Fig. 4A). The low abundance of the miRNA can be explained by the idea that four linkages must be cut to generate the miRNA:miRNA* duplex, whereas just a single cut bisects the mRNA. Perhaps very little of the duplex is fully excised, and as a result the miRNA never accumulates to sufficient levels to direct efficient target cleavage. We suggest that the processing of other intronic miRNAs might also influence the expression of their host genes—speculation bolstered by the presence of a conserved miRNA-like hairpin in the mammalian DGCR8 gene, whose protein product functions in pri-miRNA processing (Pedersen et al. 2006).

A newly identified tasiRNA gene in Arabidopsis

In a search for tasiRNA loci, we implemented a clustering algorithm that scanned the genome for phased clusters of ∼21-nt reads. This procedure found all five of the previously identified Arabidopsis tasiRNA genes (TAS1a, TAS1b, TAS1c, TAS2, and TAS3) (Supplementary Database 3) and discovered an additional locus, TAS4 (Fig. 5A), mapping between At3g25800 and At3g25790 (a MYB transcription factor).

Open in a separate window

Figure 5.

The TAS4 locus gives rise to tasiRNAs predicted to down-regulate MYB transcripts. (A) The number of reads with a 5′ terminus at each position is plotted. Bars above the axis represent sense reads; those below represent antisense reads. The miR828 complementary site is marked by a red arrow and shown below the graph, together with the fraction of 5′ RACE clones supporting the indicated cleavage site. TAS4-siR81(−), the siRNA predicted to target MYB genes, is indicated (blue bar), as is the spacing separating the phased species at each interval; spacing for the species not represented by reads is indicated in gray. (B) TAS4-siR81(−) and miR828 complementary sites in three MYB genes. Cleavage confirmed by 5′ RACE is indicated (arrows), along with the fraction of clones mapping to the cleavage site. The remaining 14 PAP2 clones mapped >20 nt from the cleavage site.

Because miRNA-directed cleavage sets the phase for and stimulates production of tasiRNAs (Allen et al. 2005; Gasciolli et al. 2005; Yoshikawa et al. 2005; Axtell et al. 2006), we searched TAS4 for miRNA complementary sites upstream of and downstream from the region that generated small RNAs. It identified a single miR828 complementary site. Cleavage at this site, validated by 5′ RACE, defined a 5′ terminus that matched that of the most proximal siRNAs arising from this locus and was in perfect register with the other predominant siRNAs (Fig. 5A). The EST mapping to this region ({"type":"entrez-nucleotide","attrs":{"text":"AU226008","term_id":"19740655"}}AU226008) corresponded to the opposite strand of the inferred primary transcript and presumably represented the RDR6-polymerized strand. Although poplar ESTs with miR828 complementary sites were found, conservation of AtTAS4 to poplar was unclear.

We also predicted three targets for TAS4-siR81(−), one of the dominant TAS4 siRNAs. The predicted targets, PAP1/MYB75/At1g56650, PAP2/MYB90/At1g66390, and MYB113/At1g66370 (Fig. 5B), encoded three MYB transcription factors that were distinct from the MYB genes targeted by miR159, and those with complementarity to miR835 and CandidateD. PAP1 and PAP2 regulate expression of anthocyanin/flavenoid and phenylpropanoid biosynthetic genes, and might also be involved in regulating leaf senescence (Borevitz et al. 2000; Pourtau et al. 2006). Intriguingly, miR828 was also predicted to down-regulate MYB113 at an independent target site (Fig. 5B), suggesting a close functional evolutionary relationship among these MYB target genes, miR828, and the TAS4 cluster. Using 5′ RACE, we identified mRNA cleavage fragments diagnostic of miR828-directed cleavage of MYB113 and tasiRNA-directed cleavage of PAP2, thereby experimentally confirming these predicted targets and demonstrating that the TAS4 locus was indeed trans-acting.

When considered as a group, the 10,469 reads from all six TAS genes were predominantly 21 nt and tended to begin with a uridine, as also observed for Arabidopsis miRNAs (Supplementary Fig. 1). MicroRNA and synthetic siRNA duplexes assemble into the silencing complex asymmetrically, such that the strand that pairs with less stability at its 5′ terminus is incorporated as the guide strand, while the other strand is degraded (Khvorova et al. 2003; Schwarz et al. 2003). Analysis of the initial 12 siRNA reads from TAS1a suggests that tasiRNAs also obey the asymmetry guidelines (Vazquez et al. 2004). The acquisition of many additional tasiRNA reads enabled us to revisit this issue. For each of the six TAS genes, each possible duplex represented by a 21mer read (including those out of register with the dominant phasing register) was considered and evaluated for which strand of the duplex yielded more reads. For 57% of the duplexes with energetically distinct terminal base-pairing, the strand with more reads was the one that appeared to be least stably paired at its 5′ terminus. Confounding this analysis, however, was the preference for a U at tasiRNA 5′ termini. If assembly or stability of the silencing complex simply preferred a U at the 5′ terminus of the guide strand, without regard for the differential pairing stabilities at the duplex ends, then there would more frequently be an A:U pair at the 5′ terminus of the guide strand than at the 5′ terminus of the passenger strand, thereby generating artifactual adherence to the pairing asymmetry guidelines. Indeed, the weak adherence vanished when repeating our analysis considering only the duplexes with an A or U at the 5′ termini of both the guide and passenger strands. Apparently, the pairing asymmetry guidelines do not apply for tasiRNAs.

Other endogenous siRNAs in Arabidopsis mapped predominantly to intergenic regions

After removing the RNAs that corresponded to the sense strand of ribosomal RNAs (rRNAs), transfer RNAs (tRNAs), small nuclear RNAs (snRNAs), and small nucleolar RNAs (snoRNAs) (Table 1; Supplemental Material), and those matching previously annotated and newly identified miRNAs and tasiRNAs, we considered the remaining RNAs that matched the nuclear genome. These included a majority (63%) of the reads and a large majority (90%) of the unique sequences (Table 1). Of the reads that did not match noncoding RNA transcripts, only 10% corresponded solely to annotated mRNAs or introns, which represented substantial depletion when considering that 49% of the sequenced genome is annotated as mRNA and intron. Of those that hit annotated mRNAs or introns, ∼46% were exclusively in the antisense orientation to a protein-coding gene. The length and 5′ nucleotide profiles of small RNAs mapping exclusively to the sense strand of genes closely resembled that of small RNAs mapping exclusively to the antisense strand of genes (Supplementary Fig. 1G,H), suggesting that a majority of the sense as well as antisense reads might be siRNAs. We considered them, together with the other small RNAs that did not match noncoding RNA transcripts, as endogenous siRNA candidates.

The candidate siRNAs included 20,720 reads that mapped to the antisense of rRNAs or to ribosomal DNA (rDNA)-like repeats but not to the mature rRNA sequences (Table 1). These are likely to include bona fide siRNAs acting by targeting the rDNA arrays for chromatin or histone modifications (Xie et al. 2004; Pontier et al. 2005; Li et al. 2006; Pontes et al. 2006). Because the fraction of the genome comprised by the rDNA arrays, as well as their copy number, was unknown, and because much of the rDNA sequence was missing from the current assembly, it was difficult to determine if this represented an enrichment.

The candidate siRNAs were mostly 24mers, the size of siRNAs associated with PolIV, heterochromatin formation, and DNA methylation (Chan et al. 2004; Xie et al. 2004; Herr et al. 2005; Kanno et al. 2005; Onodera et al. 2005; Pontier et al. 2005). As expected based on initial sequencing efforts (Tang et al. 2003), the 24mers were enriched for a 5′-terminal adenosine (Supplementary Fig. 1E,F). As for the tasiRNAs, a tendency to adhere to the pairing asymmetry guidelines for silencing complex assembly was observed only in a naïve analysis that did not consider the presumably independent preference for a particular nucleotide at the 5′ termini of the siRNA reads. When correcting for the preference for an A at the 5′ terminus of the candidate siRNAs, the strand with less stable pairing at its 5′ terminus was sequenced no more frequently than expected by chance.

Small RNAs matching annotated protein-coding genes

Some protein-coding genes had a particularly high propensity for spawning small RNAs (Supplementary Table 2). Eleven of the 20 genes most frequently hit, when read counts were normalized by the number of genomic hits and assigned only when unambiguously sense or antisense to genes, were convergently transcribed with a neighboring gene. These included an antisense gene pair, At2g16580/At2g16575, ranked 16th and 19th. Both genes have ORFs with unknown functions, with one ORF falling largely within the intron of the other. Convergent, overlapping transcription presumably generated dsRNA from which the small RNAs were derived. For the nine remaining genes in a convergent context, the 3′ termini were either uncharacterized or had nonoverlapping annotation.

RdRPs provide another mechanism for generating dsRNA. A search for reads that matched a cDNA database but failed to match the genome found 32 reads that spanned mRNA splice junctions in the antisense orientation (Supplementary Table 3). Such reads provided evidence for siRNAs generated by RdRP acting on a spliced mRNA template. The two cDNAs with the most nonoverlapping antisense hits to splice junctions encoded a TIR–NBS–LRR disease-resistance protein (At5g38850) and a basic helix–loop–helix protein (At3g23690). Both were among the top 20 genes hit by small RNAs (Supplementary Table 2). Neither had hits to introns, indicating that for these two genes the RdRP activity acted primarily on spliced templates. However, many of the genes frequently hit by small RNAs had hits to introns. Moreover, many more small RNAs matched the sense strand of mRNA splice junctions (641) than matched the antisense (Supplementary Table 3). The 20-fold difference between sense and antisense reads to splice junctions, compared with the nearly even numbers of sense and antisense reads matching protein-coding genes more generally, suggests that if RdRPs play a major role in producing siRNAs from protein-coding regions, then the templates are usually unspliced transcripts.

Candidate siRNAs derived preferentially from pericentromeric regions, with a slight preference for repeats

When candidate siRNA reads were mapped to each chromosome, plotting sequencing abundance normalized by the number of genome matches (Fig. 6), the largest peak included mostly 22mers and corresponded to a cluster of Gypsy and MuDr elements on the long arm of Chromosome 3. The next two largest peaks corresponded to two rDNA arrays with neighboring repetitive elements on Chromosome 3 and Chromosome 2. Most of the remaining small RNAs were 24mers that mapped to numerous loci dispersed throughout the genome (Fig. 6).

Open in a separate window

Figure 6.

Normalized abundance of candidate siRNAs in 0.1-Mb windows spanning the nuclear genome. Colored bars above the axis represent matches to the plus strand; colored bars below the axis represent those to the minus strand, with the colors indicating the proportion of 21mers (red), 24mers (light and dark blue), 24mers with a 5′ A (light blue), and other lengths (yellow). Below the siRNA profiles are histograms plotting the fraction of nucleotides falling within annotated protein-coding genes (black; scale, 0%–100%) and the fraction falling within repetitive elements annotated by RepeatMasker (gray; scale, 0%–100%). Centromeres are indicated by solid gray bars, and heterochromatic knobs are indicated by hashed gray bars.

Although the Arabidopsis genome is relatively compact, repetitive loci are abundant. They are most dense at and near the centromeres, and their density gradually tapers off in the 2–3 Mb on both sides of each centromere as protein-coding density increases (Fig. 6). Along the remainder of the chromosomes, repeats are present but occur at much lower density. The siRNA density did not peak at the same regions as the repeat density peaked and was instead greatest in the proximal and distal pericentromeric regions, characterized by an intermediate density of both repeats and annotated protein-coding genes (Fig. 6). A look at a diagnostic centromeric repeat class, the ∼180-base-pair (bp) repeat satellite arrays (Copenhaver et al. 1999; Nagaki et al. 2003), illustrated this result. Only 2386 reads (1055 unique sequences) matched ∼180-bp centromeric satellite repeats annotated by RepeatMasker. This was 0.43% of our reads, whereas the annotated ∼180-bp repeat represented 0.39% of the current genome assembly. Because many ∼180-bp repeats are missing from the assembly, this slight apparent enrichment was undoubtedly an overestimate; siRNAs deriving from unassembled repeats would artifactually add to the perceived density at any assembled repeats that they match. The unremarkable correspondence between candidate siRNAs and known heterochromatin was also illustrated at the heterochromatic knobs, which were rich in candidate siRNAs, but not more enriched than were the pericentromeric regions that surrounded them (Fig. 6).

The observation that siRNAs were often associated with repeats, but were not highest where repeats were most dense, raised the question of whether siRNAs derived preferentially from repeat loci. Transposons, retroelements, and low-complexity sequences identified by RepeatMasker comprised ∼15% of the current genome assembly. Of the 558,481 candidate siRNA reads, 188,502 (34%) hit these regions annotated by RepeatMasker—a modest, twofold enrichment over the 15% that would have been expected if the siRNAs derived uniformly from repetitive and nonrepetitive regions throughout the genome. The twofold enrichment was largely attributed to the depletion of both repeats and siRNA matches within annotated protein-coding genes. Of the 51% of the genome that fell between annotated protein-coding genes, ∼30% corresponded to repeats annotated by RepeatMasker. Of the 491,180 siRNAs mapping between annotated protein-coding genes, 188,502 (38%) hit regions annotated by RepeatMasker, indicating only a 1.2-fold preference for repeat regions within intergenic regions. The modest preference for repeat regions decreased further when excluding the 20,720 reads deriving from rDNA repeats. When considering only those RNAs associated with AGO4 (Qi et al. 2006), this slight enrichment increased, but not by much. Local (<100 kb) inverted-repeat regions did not appear to be overrepresented among genomic hits of intergenic siRNA candidates that fell outside of repetitive elements identified by RepeatMasker. Of the repetitive DNA detected by RepeatMasker, 80% was either class I (retrotransposon derived) or class II (DNA transposon derived). About 95% of the repeat-associated siRNAs corresponded to these two classes, in the proportion expected based on the contribution of these two classes to the genome. Representation of some of the more well-characterized transposable element families is listed in Supplementary Table 6.

Initial processing of reads

cDNA sequences were extracted from raw reads, excluding reads lacking perfect matches to the most proximal 11 nt of both adapter sequences. Unique sequences were mapped to the TAIR/NCBI genome version 6.0 (November 2005), chloroplast and mitochondrial genomes, and rRNA and tRNA sequences from TAIR (http://www.arabidopsis.org), the 5S rRNA Database (http://biobases.ibch.poznan.pl/5SData), the published chromosomal rRNA sequences (Unfried et al. 1989; Unfried and Gruendler 1990), the Arabidopsis tRNA database (http://lowelab.ucsc.edu/GtRNAdb/Athal), the Arabidopsis snoRNA database (http://bioinf.scri.sari.ac.uk/cgi-bin/plant_snorna/arabidopsis; Brown et al. 2001), miRBase 7.1, and the TIGR v2 Arabidopsis Repeat Database. Repetitive regions of the genome were mapped by RepeatMasker (http://www.repeatmasker.org) using AtRepbase (Jurka et al. 2005).

miRNA identification

Twenty-nucleotide to 24-nt sequences with more than one read, 16 or fewer hits to the genome, and no matches to annotated noncoding RNA were folded using RNAfold with 330 nt of upstream and downstream flanking sequence. For efficiency, candidate reads were clustered and only the most abundant in a set of overlapping hits was considered. Structures were evaluated using mirCheck, a script that assesses the quality of a foldback based on a battery of parameters that capture known miRNA hairpins (Jones-Rhoades and Bartel 2004). Hairpins that passed this initial filter were then manually screened. Designation as a miRNA required (1) a foldback in which the duplex region that included 25 nt centered on the most frequently sequenced read had less than eight unpaired nucleotides (summing unpaired nucleotides on both arms of the stem) and no more than three consecutive unpaired nucleotides, of which no more than two were asymmetrically bulged; (2) a sequenced miRNA* species (paired to the miRNA within the duplex with 2-nt 3′ overhangs) or for candidates with three or more reads, a slight variant of the miRNA*; and (3) a sense:antisense read ratio >0.90. In practice, all but six foldbacks that passed manual inspection and were named as miRNA loci had a ratio >99% (Supplementary Database 2). MIR824, which has >330 nt between the miRNA and the miRNA*, was found in a separate analysis of genomic regions with abundant 21-nt reads.

Phased siRNA discovery

For each unique small RNA sequence (excluding those matching miRNAs, other noncoding RNAs, or protein-coding genes) a 500-nt window, anchored at one end by that sequence, was evaluated for phased small RNAs. If three or more unique 20- to 23-nt sequences with nonoverlapping hits existed in the window, each was evaluated for phasing with any of the others in the window, allowing ±2 nt of divergence from perfect 21-nt phasing. Phased sequences were extracted from the window and the process was repeated for any remaining 20- to 23-nt sequences until two or fewer unique nonoverlapping hits were left. Each potential phase was then evaluated according to five parameters: (1) a count of all unique sequences in phase; (2) a count of the reads in the window; (3) a hits-normalized score, whereby the sum of the read frequencies of all phased sequences was divided by the sum of their genomic hits; (4) a normalized 21mer score, which divided the sum of the 21mer reads in phase by the sum of the reads in the window; and (5) a phasing score, which divided the sum of the reads in phase by the sum of the reads in the window. Cutoffs for each score were empirically adjusted to find values that captured all known tasiRNA clusters but restricted the number of false positives.

Target site prediction for miRNAs and TAS4 siRNAs

Patscan was used to search for near matches (up to six mismatches, or four mismatches and one bulged nucleotide) in TAIR version 6.0 Arabidopsis cDNA database (http://www.arabidopsis.org) to each miRNA, and target sites were scored as described (Allen et al. 2005). To assess performance, we applied this algorithm to a control set of diverse Arabidopsis miRNAs, choosing the most frequently sequenced miRNA variant to represent each known miRNA family for which mRNA targets have been experimentally validated by 5′ RACE, as listed in Jones-Rhoades et al. (2006). We also generated 10 different shuffled cohorts of these 22 miRNAs, preserving dinucleotide composition. Signal:noise ratios were calculated by comparing the total predictions for authentic miRNAs (signal) and the average for shuffled cohorts (noise). Analogously selected cohorts were also used to estimate specificity of target prediction for the newly identified miRNAs.

Hotspot identification

For candidate siRNAs, nonoverlapping 500-bp windows were ranked by scoring small RNA density as a sum of abundances of all sequences in the window, normalized by the sum of the total number of times each sequence hit the genome. Top-ranking windows were then used as seeds for extension in both directions until a 500-bp window lacking any siRNA hits was encountered.

siRNA duplex asymmetry determination

Conceptual duplexes with 2-nt 3′ overhangs were constructed by determining the reverse complement of the sequenced strand. The terminal three pairs (two nearest neighbors) on each end of the duplex were analyzed, comparing sums of the two −ΔG°₃₇ nearest-neighbor parameters for RNA duplex stability (Xia et al. 1998). More complex algorithms were also implemented and yielded similar conclusions (Supplemental Material).

RNA gel blot analysis of small RNA expression

RNA gel blots for miRNAs were performed using mutants and protocols as described previously (Vaucheret et al. 2004; Vazquez et al. 2004). Blots were probed with ³²P-end-labeled DNA [miR822, miR823, miR391(5′), miR771, miR824, miR775] or LNA [miR391(3′), miR839, miR846] oligonucleotides, each complementary to the entire length of the miRNA.

5′ RACE

5′ RACE was performed as described in Jones-Rhoades and Bartel (2004), except that RNA samples were obtained from whole siliques or seedlings, gene-specific primers (Supplemental Material) were designed to be 70–400 bases from the predicted cleavage site, and the first gene-specific amplifications for DCL1, PAP2, and MYB113 were done with the GeneRacer 5′ outer primer and were followed by two nested amplifications done with the GeneRacer 5′ nested primer.

We thank M. Axtell for many helpful discussions, G. Bell for computational advice, A. Hall and J. Jurek (Preuss Laboratory) for pilot analysis of centromeric candidate siRNAs, and M. Jones-Rhoades and G. Ruby for scripts. This work was supported by the Prix Louis D. award from the Institut de France and a grant from the NIH to D.B.