Genome-Wide Analysis of the RNA-DEPENDENT RNA POLYMERASE6/DICER-LIKE4 Pathway in Arabidopsis Reveals Dependency on miRNA- and tasiRNA-Directed Targeting^[W],[OA]

Analysis of TAS1, TAS2, and TAS3 Families

A detailed analysis of small RNAs from known TAS loci was done using populations of sequenced small RNAs from wild-type (Col-0), dcl mutant (dcl1-7, dcl2-1, dcl3-1, and dcl4-2), and rdr mutant (rdr1-1, rdr2-1, and rdr6-15) plants (Fahlgren et al., 2007; Kasschau et al., 2007; http://asrp.cgrb.oregonstate.edu/db/). tasiRNAs from TAS1, TAS2, and TAS3 families were represented by at least 836 reads (Table 1), and all were RDR6 dependent (Table 2). Only two reads were obtained for the TAS4 locus. Small RNAs from TAS loci were predominantly 21 nucleotides in length, and this size class was greatly reduced in the dcl4-2 mutant (all values normalized for sequence population size; Table 2). Additionally and as shown previously, tasiRNAs in these sequenced populations were dependent on DCL1 (see Supplemental Table 1 online; Kasschau et al., 2007).

The fact that validated tasiRNAs, or tasiRNAs with high-quality predicted targets, arise from the TAS locus positions residing up to 12 processing cycles away from the initiation cleavage sites (Figure 1) points to the importance of accurate phasing to generate small RNAs with sufficient complementarity to target sequences. Therefore, the phasing patterns of small RNAs from the miR173-guided cleavage site (TAS1 and TAS2 loci) and 3′ miR390-guided cleavage site (TAS3 loci) were analyzed. Small RNA sequence data from wild-type and rdr2-1 mutant (no effect on tasiRNA biogenesis) plants were used to develop an algorithm that converted phase position and amplitude (small RNA abundance) data into phase signals. First, a sense-antisense small RNA consolidation step was developed to unify sequence data from both strands. Small RNA pairs from sense-antisense positions that were offset by two nucleotides, with the overhang at the 3′ ends, were combined. This accounted for the possibility that either sense or antisense small RNA from a given small RNA duplex might accumulate. Second, a phase score, P, was calculated using the following formula:

in which n = number of phase cycle positions occupied by at least one small RNA read within an eight-cycle window, and k = the total number of reads for all small RNAs with consolidated start coordinates in a given phase within an eight-cycle window. Phase cycle length was set at 21 nucleotides. A positive phase signal, therefore, was limited by the power function (n − 2) to those loci with small RNAs occupying at least three cycle positions in a phased, 21-nucleotide register. Requiring occupation of three positions in an eight-cycle window limited the possibility of a positive phase score by chance. It should also be noted that the phase score amplitude is influenced by both small RNA abundance and number of positions occupied. This calculation was applied in a scrolling window format with single nucleotide advances, and scores were assigned to the fourth cycle position in phasing plots (Figure 1).

Small RNAs were detected unidirectionally downstream of the miR173-guided cleavage site from each TAS1 and TAS2 locus (Figure 1A; http://asrp.cgrb.oregonstate.edu/db/). Phase plots clearly revealed the dominant phasing event to be coincident with the miR173 target site for each locus, although both the raw read and phase transform plots revealed out-of-phase noise at each locus. Phase drift was detected at TAS1c, with small RNAs accumulating in the phase position one nucleotide forward from the initial miR173-directed phase, starting at cycle five (Figure 1A). Between cycles 1 and 5, the phase-forward:initial phase ratio of small RNA reads at TAS1c was 1:5; however, between cycles 6 and 10, the phase-forward:initial phase ratio at TAS1c was 22:1. For TAS3a, small RNAs were detected exclusively upstream of the 3′ miR390 target site (Figure 2A). The miR390-initiated phase signal, however, was slightly lower than an out-of-phase signal, which likely originated from a secondary targeting event guided by an siRNA from the 5′D2(−) position (Figure 2B; Allen et al., 2005; Axtell et al., 2006). Accurately phased, miR390-initiated RNAs were detected close to the 3′ target site at TAS3a, but phase-forward drift was detected after several cycles (Figure 2). Between cycles 1 and 6, the phase-forward:initial phase ratio of small RNA reads at TAS3a was 1:3. However, the phase-forward:initial phase ratio was 32:1 at cycles 7 and 8. This may be an adapted feature of TAS3a, as the functional tasiRNAs occur in the forward phase at the 5′D7(+) and 5′D8(+) positions (positional nomenclature from Allen et al. [2005]).

The basis for phase-forward drift at TAS3a and TAS1c could be due, at least in part, to non-21-nucleotide siRNAs resulting from misprocessing by DCL4 or alternate DCLs. Alternate-size small RNAs were detected at relatively low levels for most TAS loci, although an abundant 22-nucleotide RNA from TAS1c was detected at the fifth processing cycle from the miR173 target site (http://asrp.cgrb.oregonstate.edu/db/). Loss of DCL4 in the dcl4-2 mutant resulted in off-sized RNAs at each TAS locus (Table 2, Figure 1B). Although the abundance of total small RNAs arising from TAS1, TAS2, and TAS3 loci was only modestly affected, the ratio of 21-nucleotide to non-21-nucleotide RNAs was considerably lower in the mutant (Table 2). Importantly, this was accompanied by loss of miRNA-guided phasing (Figures 1B and 2B), and loss of functional tasiRNAs, at each locus. Therefore, although DCL4 is not required to process dsRNA formed from each TAS locus, it is clearly necessary for biogenesis of active, properly phased tasiRNAs.

A search for previously unrecognized members of TAS1, TAS2, and TAS3 families revealed two additional TAS3 family members (TAS3b and TAS3c), which contained sequences similar to ARF3- and ARF4-targeting tasiR2141 from the originally characterized TAS3a locus (Figure 2C). TAS3b was previously identified as a tasiRNA-like region by Lu et al. (2006). Unlike TAS3a, both TAS3b and TAS3c contain only a single tasiR2141-like sequence. RDR6-dependent small RNAs of predominantly 21 nucleotides were detected from both TAS3b and TAS3c (Table 2; http://asrp.cgrb.oregonstate.edu/db/), although TAS3c generated relatively few RNAs. Importantly, nucleotide sequence alignments and target validation assays revealed a functional miR390-guided cleavage site on the 3′ side of both TAS3b and TAS3c tasiRNA-generating regions (Figure 2D). The tasiR2141-like sequence was detected in the miR390-initiated phase for TAS3c, but in the +4 forward phase for TAS3b. Transcripts for each TAS3 locus were detected, and 5′ transcription start sites were mapped (see Supplemental Figure 1 online). Therefore, TAS3 tasiRNAs arise from a small multigene family, all members of which have conserved miR390-guided target sites.

Interestingly, in addition to the conserved 3′ miR390 target sites and tasiR2141-related sequences, conserved miR390 target sites on the 5′ side of TAS3a, TAS3b, and TAS3c tasiRNA-generating sequences were detected (Figure 2D). In each case, however, the site deviates from canonical plant miRNA target sites in that mispairs and G:U pairs occur at positions 9 to 11 (relative to the 5′ end of miR390; Figure 2D). Based on mutational and informatic analyses of naturally occurring target sites, mispairs and G:U pairs at these positions were predicted to have restricted miR390-guided cleavage (Jones-Rhoades and Bartel, 2004; Mallory et al., 2004; Allen et al., 2005; Schwab et al., 2005). Indeed, 5′ rapid amplification of cDNA ends (RACE) analysis failed to detect any cleavage products corresponding to an active 5′ miR390 target site, although 5′ RACE products were detected at low frequencies at several sites flanking the noncanonical target site at each locus (Figure 2D). Given the high degree of conservation, we conclude that the site is likely important, but functional in a mode that does not depend on target cleavage. Axtell et al. (2006) also detected the noncanonical 5′ miR390 site in TAS3a and showed that it functions in a noncleavage mode in vivo for TAS3-derived tasiRNA biogenesis and function.

Genome-Wide Scans for tasiRNA-Like Loci

To identify loci with tasiRNA-like features or that spawn RDR6/DCL4-dependent siRNAs, a genome-wide scan using Columbia (Col-0) and mutant small RNA sequence libraries was done. These loci were predicted to yield RDR6-dependent, 21-nucleotide RNAs in phased arrays. A set of filters to identify clusters of 21-nucleotide small RNAs that were present in either Col-0 and/or rdr2-1 mutant libraries, and absent in rdr6-15 mutant libraries, yielded overlapping sets of 50 and 46 loci in scans of genomic and transcript sequences, respectively (Figure 3; see Supplemental Table 1 online). Thirty-six of these loci were known to be targeted by at least one miRNA or tasiRNA (see Supplemental Table 1 online). These included overlapping sets of 30 miRNA targets (including seven known TAS loci) and 18 tasiRNA targets. Twenty of the 36 loci contained multiple miRNA or tasiRNA target sites. Therefore, among the 50 RDR6-dependent small RNA clusters identified, 72% arose from miRNA/tasiRNA-targeted transcripts and 40% from multiply targeted transcripts. Small RNAs (21 nucleotides) from these loci were generally lost, or accumulated to low levels, in both the dcl1-7 and dcl4-2 mutants (see Supplemental Table 1 online), indicating dependence on both miRNA-based and tasiRNA-like mechanisms.

Eighteen loci, including 17 that were targeted by miRNAs/tasiRNAs, contained at least three in-phase (21-nucleotide cycle) positions occupied by at least one small RNA within an eight-cycle window, yielding peak phase scores ranging from 1.4 to 33.1 (Table 2). Except for the known TAS1, TAS2, and TAS3 loci that came through the screen, no other non-protein-coding loci were identified. However, three noncoding loci yielded RDR6-dependent, 21-nucleotide RNAs that were phased using a more permissive algorithm in which signal required occupation of only two out of 12 cycle positions (see Supplemental Table 1 online). Two of these loci (from chromosome 4, between coordinate positions 1472767 and 1476524) were embedded in a region with duplications of ubiquitin-like genes or gene sequences (such as At4g03350; see Supplemental Table 1 online). In fact, both small RNA-generating loci contained extensive, continuous similarity to these genes and appear to have resulted from recent duplication events. The other noncoding locus was identified between opposingly oriented genes At5g27650 and At5g27660 but was represented by too few small RNAs in Col-0 plants to allow a meaningful analysis (see Supplemental Table 1 online). Lu et al. (2006) recently identified two TAS-like loci from noncoding regions (Chr1:25282658-25283382 and Chr4:13295428-13296124). These loci were not identified here because small RNAs from both loci were RDR6 independent (http://asrp.cgrb.oregonstate.edu/db/). The recently identified TAS4 locus (Rajagopalan et al., 2006) was not identified due to a lack of sequenced small RNAs.

Four gene families were represented among the 12 protein-coding loci that emerged from the phase scan (Table 2). These included a Toll/Interleukin1 receptor–nucleotide binding site–leucine-rich repeat (TIR-NBS-LRR) disease resistance gene (At5g38850; Meyers et al., 2003), the TIR1-like AFB3 gene encoding an F-box protein involved in auxin signaling (At1g12820; Dharmasiri et al., 2005), a cation/hydrogen antiporter gene ATCHX18 (At5g41610; Sze et al., 2004), and nine genes encoding closely related PPR proteins (discussed below). Transcripts from AFB3 and ATCHX18 were targeted and phase-initiated by miR393 and miR856, respectively. Interestingly, two miRNAs (miR856 and miR780) were shown to target the ATCHX18 transcript (Fahlgren et al., 2007), but phasing is initiated only at the miR856 site. Although the TIR-NBS-LRR gene generated abundant, highly phased siRNAs, we did not identify a miRNA or tasiRNA that initiates phasing. In addition to the miRNA target transcripts that yielded phased siRNAs, targets for miR168 (AGO1), miR472 (several CC-NBS-LRR genes), miR393 (AFB2), miR773 (MET2-like), and miR172 (AP2) yielded RDR6-dependent siRNAs that failed to generate a positive phase score (see Supplemental Table 1 online). siRNAs from many of these genes were also identified in other analyses (Axtell et al., 2006; Lu et al., 2006; Rajagopalan et al., 2006; Ronemus et al., 2006; see Supplemental Table 1 online).

Targeting of PPR-P Gene Transcripts and Formation of Phased siRNAs

By far, the most abundant phased, 21-nucleotide siRNA populations were identified from two clusters of PPR genes on chromosome 1 (Table 2). The Arabidopsis thaliana genome contains at least 448 PPR-related genes, which encode putative RNA binding proteins with repeats of a 35–amino acid motif (Lurin et al., 2004). Some PPR proteins are known to regulate RNA processing, stability, editing, or translation in chloroplasts or mitochondria (Lurin et al., 2004; Cushing et al., 2005; Kotera et al., 2005; Okuda et al., 2006; Shikanai, 2006). The PPR proteins are classified into subfamilies containing only tandem arrays of PPR motifs (PPR-P) or containing diversified PPR motifs interspersed with non-PPR domains (Aubourg et al., 2000; Lurin et al., 2004; Rivals et al., 2006). A few PPR-P genes from Raphanus sativus (Rfo), Zea mays (Rf-1), and Petunia hybrida (Rf) are known restorers of fertility in lines with cytoplasmic male sterility conferred by certain mitochondrial ATPase alleles (Bentolila et al., 2002; Desloire et al., 2003; Kazama and Toriyama, 2003; Wang et al., 2006; Chase, 2007).

Two regions on chromosome 1 contain 21 PPR-P genes and pseudogenes, of which 17 (81%) spawned primarily 21-nucleotide siRNAs (Figure 4A). These were dependent on RDR6, DCL4, and DCL1 but independent of RDR1, RDR2, DCL2, and DCL3 (Figure 4A, Table 2; see Supplemental Table 1 online). Phylogenetic analysis of the entire PPR-P subfamily revealed that the siRNA-generating loci belonged to a clade of 28 PPR-P genes (hereafter referred to as the PPR-P clade) that has undergone relatively recent expansion (see below). Using expression profiling arrays, the steady state levels of 23 of the 28 PPR-P clade members in inflorescence tissue was measured in each of the dcl and rdr mutants, in hen1, hst, and hyl1 mutants, and in Col-0 and Landsberg erecta (Ler) plants. The levels of transcripts for the siRNA-generating PPR-P genes were generally elevated in the mutants with miRNA or tasiRNA defects (dcl1, dcl4, rdr6, hen1, hst, and hyl1) but not in mutants with defects in other small RNA classes (dcl2, dcl3, rdr1, and rdr2) (Figure 4B). Fifty-three percent of siRNA-generated clade genes were significantly upregulated in at least one mutant involved in miRNA and tasiRNA biogenesis (23-gene set, analysis of variance, false discovery rate = 0.05). These data indicate that a high proportion of genes within the clade not only spawn siRNAs with genetic requirements similar to those of tasiRNAs, but that the PPR-P transcripts are repressed by posttranscriptional silencing mechanisms.

Computational analysis and experimental validation assays revealed 11 Arabidopsis miRNAs and tasiRNAs with potential or validated function to guide cleavage of PPR-P mRNAs (Rhoades et al., 2002; Allen et al., 2004; Sunkar and Zhu, 2004; Yoshikawa et al., 2005). These include both forms of miR161 (miR161.1 and miR161.2), miR400, four tasiRNAs from the TAS1 family, and four tasiRNAs from the TAS2 locus (Table 2; see Supplemental Figure 2 online). Remarkably, all of these small RNAs were predicted or validated to target PPR-P transcripts predominantly from the clade, with most clade members targeted by multiple small RNAs (Figure 5A, Table 2; see Supplemental Table 1 online). Five of these genes were validated here by 5′ RACE as both miR161 and TAS2 3′D6(−) targets (Table 2; see Supplemental Figure 2 online). Combined with previously validated transcript targets (Allen et al., 2004; Yoshikawa et al., 2005), eight PPR-P transcripts were validated experimentally as targets for at least one small RNA. As a group, these miRNAs and tasiRNAs were predicted to target a total of 40 PPR-P transcripts, including all 28 PPR genes in the clade. Although siRNA generation and targeting of clade PPR-P transcripts by miRNAs and tasiRNAs were correlated (see below), it was also possible that the high siRNA levels were a reflection of high accumulation levels of the target transcripts. To determine if PPR gene–derived siRNA levels were correlated with transcript abundance, the mean expression value in wild-type plants (Col-0 inflorescence tissues) and repeat-normalized siRNA reads were compared. There was no significant correlation between siRNA abundance and transcript level (Pearson correlation coefficient = 0.077; Figure 4C), suggesting that siRNA abundance was related to miRNA and/or tasiRNA targeting rather than to transcript level.

Eight of the nine PPR-P genes identified in the genome-wide phase scan contained phase initiation sites corresponding to miR161.1, miR161.2, or TAS2 3′D6(−) (Table 2, Figure 6). In fact, loci spawning the highest numbers of siRNAs appeared to correspond to genes with transcripts that were dual-targeted by both miR161 forms and TAS2 3′D6(−) (Table 2). In the case of At1g62910, which contains a duplicated coding region resulting in nearly twice as many PPR domains as most other clade proteins, two sets of miR161 and TAS2 3′D6(−) target sites and phased siRNAs were identified (Figure 6). None of the other phase patterns from the PPR-P sequences were coincident with any other predicted or validated small RNA-guided target sites, although target sites downstream of miR400 were not tested in 5′ RACE assays. Curiously, although miR400 was predicted to interact with all but one PPR transcript yielding phased siRNAs, no phased RNAs were detected in the vicinity of the miR400 target sites. In fact, each of five PPR sequences tested were negative in miR400 target validation assays, even though each miR161 and TAS2 3′D6(−) target site tested from the same genes was positive (see Supplemental Figure 2 online). miR400 target validation assays were only done using Col-0 inflorescence tissue, and it is possible that miR400-guided cleavage occurs in other tissues or under different conditions than those tested here.

The phasing patterns from miR161 and TAS2 3′D6(−) target sites were distinctive. Phased siRNAs from the TAS2 3′D6(−) target site generally occurred at the 3′ side (Figure 6); indeed, nearly all small RNAs (phased and nonphased) from the targeted PPR-P genes occurred at the 3′ side of the TAS2 3′D6(−) target site, much like tasiRNAs that originate from the 3′ side of the miR173-guided cleavage site in TAS1 and TAS2 transcripts (Figure 1). By contrast, nearly all of the phased siRNAs initiated at the miR161.1 and miR161.2 target sites occurred on the 5′ side of the initiation sites (Figure 6), which is reminiscent of tasiRNAs originating from the 5′ side of the miR390-guided cleavage site at TAS3 loci (Figure 2).

The distribution of predicted or validated target sites for the 11 miRNAs or tasiRNAs within the Arabidopsis PPR-P clade sequences was analyzed. The miR161.1/miR161.2 and miR400 target sites were separated by several hundred nucleotides in their common target transcripts (Figure 5B). However, in alignments of individual PPR domain–coding regions, it was clear that miR161 and miR400 target sequences span domain-homologous sequences (Figure 5C). In alignments of miR161 and miR400 foldback sequences, a modest similarity was detected (ClustalW score = 67), raising the possibility that their respective MIRNA genes share a common origin. The TAS1 and TAS2 tasiRNAs targeted sequences within PPR domain 3′ coding sequences, which were generally more diverse than the domain 5′ regions (Figure 5C). Although some target sequences overlap domain-homologous regions, the tasiRNA target sites represent positions across 65 nucleotides of domain sequence space (Figure 5C). This comparison suggests that targeting of PPR transcripts by miRNAs and tasiRNAs likely originated independently on multiple occasions. Interestingly, two Populus miRNAs (miR475 and miR476) also have target sites within PPR protein–coding transcripts (Lu et al., 2005). These miRNAs target PPR domain–coding sequences that are spatially distinct from the Arabidopsis small RNA target sites, again suggesting independent origination of PPR transcript targeting by distinct small RNAs.

Secondary Targeting of PPR-P Gene Transcripts

Prolific generation and amplification of RDR6/DCL4-dependent siRNAs after initial targeting by miRNAs and tasiRNAs could potentially reinforce silencing of PPR-P clade transcripts. They could also potentially trigger a cascade of reactions in which silencing spreads through the PPR gene family. To assess the potential of secondary targeting by siRNAs from the 28-member PPR-P clade, genome-wide targets were computationally predicted for antisense 21-nucleotide small RNAs arising from clade sequences. Targets were predicted based on the same methods used to predict miRNA targets (Figure 7A; Allen et al., 2005; Fahlgren et al., 2007). Secondary siRNA targets were grouped by target score and by their identity as members of the clade, nonclade PPR-P genes, or non-PPR genes. Hits resulting from perfect complementarity between siRNA and PPR gene transcripts (target score = 0) were eliminated, as these could not be assessed in terms of potential targeting in trans. Clade PPR-derived siRNAs were more likely to target members of the clade, regardless of their predicted target score (Figures 7B and 7C). The number of nonclade PPR and non-PPR hits by clade-derived siRNAs increased at higher target score levels (2.5 to 4.0; Figure 7B), as expected for increasing the allowance for mispairs and G:U pairs. This was similar to that observed for miRNA target predictions, where off-target hits were detected at higher scores using both authentic miRNAs and sequence-shuffled miRNAs (Figure 7A). The percentage of non-clade PPR hits, however, was much greater than the percentage of non-PPR hits (Figure 7C), suggesting a greater likelihood of clade-derived siRNAs interacting with other PPR transcripts rather than non-PPR transcripts. Whereas all PPR-P clade members, except for one, were potential targets of multiple clade-derived siRNAs, only 72 of the 255 PPR-P subfamily members were predicted as targets of at least one siRNA. Thirty-three nonclade PPR-P transcripts were predicted targets for two or more clade-derived siRNAs (Figure 7D). Most of these had predicted target scores of 3.5 to 4.0 (Figures 7B and 7C).

Does targeting of siRNAs from clade transcripts lead to amplification of secondary siRNAs from nonclade members? In total, 25 nonclade PPR-P loci (11%) yielded at least one siRNA. This is actually below the genome-wide level of protein-coding genes with identity to at least one small RNA (18.2%; Kasschau et al., 2007). Among the 72 and 33 nonclade PPR-P transcripts that were plausible predicted targets of at least one or two clade-derived siRNAs, respectively, only nine (12.5%) and six (18.2%) yielded siRNAs. Neither of these percentages were significantly different (P > 0.27) from the genome-wide percentage. Additionally, four of these nonclade PPR-P small RNA-generating loci (At1g43980, At4g17910, At3g60960, and At5g27300) overlapped with transposons or retroelements and produced predominantly 23- to 24-nucleotide small RNAs (Figure 7D), a feature associated with RDR2/DCL3-dependent siRNAs. Furthermore, there was no correlation between siRNA generation from nonclade PPR-P loci and predicted targeting by any of miR161, miR400, or TAS1 or TAS2 tasiRNAs (Figure 7D; data not shown). These results argue against secondary targeting by clade-derived siRNAs on nonclade target transcripts.

TAS Loci and siRNA Phasing

The tasiRNA biogenesis mechanism employs a novel use for miRNAs as site-specific RNA processing factors. Functional tasiRNAs must be in (or very close to) a 21-nucleotide register with the miRNA-guided cleavage site, otherwise the resulting siRNAs would have limited targeting potential (Vazquez et al., 2004b; Allen et al., 2005; Yoshikawa et al., 2005). Importantly, correct phasing also depends on DCL4. Although small RNAs form in the absence of DCL4, due to alternate activities provided by DCL2 and DCL3 (Dunoyer et al., 2005; Gasciolli et al., 2005; Xie et al., 2005), phasing is lost in the dcl4-2 mutant. Given the developmental phenotypes and derepression characteristics of tasiRNA target transcripts in dcl4 mutants, we infer that the off-sized, nonphased small RNAs that originate from TAS loci in dcl4-2 are largely nonfunctional. The structure of DICER from a lower eukaryote indicates that size characteristics of siRNAs are controlled by the spacing between the binding pocket of the PAZ domain, which grips the end of a substrate dsRNA with a 2-nucleotide 3′ overhang and the catalytic features at the dual active sites (Macrae et al., 2006). The 21-nucleotide siRNA-generating structure of DCL4, therefore, is a critical property of the enzyme and has been incorporated as a molecular ruler with physiologic importance. An open question concerns phase drift, which occurs after several DCL4 processing cycles of several tasiRNA precursors. Phase drift occurs at discrete points, such as processing cycle 5 for the TAS1c precursor (Figure 1), and may relate to unusual properties of the dsRNA interval interacting with DCL4 at a given cycle.

Two additional TAS3 loci (TAS3b and TAS3c) were identified (Figure 2; see Supplemental Table 1 and Supplemental Figure 1 online). These were assigned to the TAS3 family for two reasons. First, they contained sequenced or predicted tasiRNAs that were similar to the ARF3/ARF4 transcript-interacting tasiRNAs from TAS3a (Figure 2C; Allen et al., 2005; Williams et al., 2005; Hunter et al., 2006). Second, they contained conserved miR390 target sites on both the 5′ and 3′ sides of the tasiRNA sequences. Whether or not the TAS3b and TAS3c tasiRNAs are functional remains to be determined. A mutant carrying a T-DNA insertion in TAS3a displayed accelerated phase-change phenotypes similar to those seen in rdr6 and dcl4 mutants (Adenot et al., 2006), suggesting that TAS3b and TAS3c cannot substitute for TAS3a function. For each TAS3 family member, the 3′ miR390 target site, but not the 5′ target site, was cleaved (Figure 2D). The base pair interactions between miR390 and 5′ target sites contain key mispairs that prevent cleavage, and this is a highly conserved and functional feature of the TAS3 family (Axtell et al., 2006). The consequences of dual targeting of transcripts by miRNAs are discussed in greater detail below.

Genome-Wide Populations of RDR6-Dependent siRNAs

Analysis of the TAS1, TAS2, and TAS3 tasiRNA biogenesis and phasing patterns allowed development of an algorithm to detect RDR6-dependent small RNA phase signals from across the Arabidopsis genome (Figure 3). Three categories of RDR6-dependent siRNAs were detected (see Supplemental Table 1 online). In the first, tasiRNAs from six TAS loci were identified. Neither TAS3c, nor the recently recognized TAS4 (Rajagopalan et al., 2006), loci emerged from the phase scan due to low tasiRNA abundance. Second, three RDR6-dependent loci that were annotated as protein noncoding were identified, although none of these yielded a strong phase signal. Third, 38 protein-coding gene transcripts were identified as RDR6-dependent, as well as DCL1- and DCL4-dependent, generators of 21-nucleotide siRNAs. Nearly all of these are targets of at least one miRNA and/or tasiRNA, and nearly half were targeted by multiple small RNAs. Twelve of these loci yielded strong phase signals that were initiated by one or more miRNA/tasiRNA-guided cleavage events. The absence of a phase signal from the remaining 26 protein-coding genes may have been due, in part, to relatively low siRNA abundance, with more sequencing reads likely to yield more phase signals. Among the loci with a phase signal, only one (At5g38850, encoding a TIR-NBS-LRR protein) lacked a known miRNA or tasiRNA phase initiator (Table 2). These data reveal the strong linkage between endogenous deployment of the RDR6/DCL4 siRNA biogenesis system and small RNA-guided initiation.

The high proportion of dual-targeted transcripts that yield phased siRNAs through the RDR6/DCL4-dependent pathway strongly supports the two-hit trigger hypothesis proposed by Axtell et al. (2006). How efficiently a transcript gets routed through the RDR6/DCL4 siRNA biogenesis pathway may depend on how many small RNA target sites exist and on at least one of these target sites being cleaved. Small RNA-guided effector complexes may interact with factors, such as SGS3, that stabilize one of the cleavage products (Mourrain et al., 2000; Peragine et al., 2004; Yoshikawa et al., 2005). This would then lead to recruitment of RDR6 to the stabilized RNA fragment for dsRNA synthesis and subsequent DCL4-mediated processing in an end-dependent manner. The features governing which small RNA-guided cleavage fragment is stabilized and converted to dsRNA is not clear. Perhaps a cleavage product that is associated with another small RNA-AGO effector complex, such as that which associates with the 5′ fragment from miR390-guided TAS3 transcript cleavage products, marks a preferred template. The abundance of small RNA target sites, some of which do not guide cleavage, may be what directs the PPR transcripts so effectively through the pathway.

The miRNA and RDR6/DCL4-dependent siRNA systems may have coevolved as a posttranscriptional silencing mechanism to limit the activity of the RDR6/DCL4 system on most endogenous transcripts and RNAs. The need for one, and preference for multiple, small RNA–mediated initiation events may disqualify most cellular RNAs as targets. The preference for multiple targeting events also explains why most single-site miRNA target transcripts, which are the vast majority of plant miRNA targets, do not spawn siRNAs (Axtell et al., 2006). Furthermore, it may explain why viral RNAs and other targets of complex populations of siRNAs effectively amplify silencing signals (Brodersen and Voinnet, 2006).

Targeting of a Rapidly Expanding PPR-P Gene Clade

Except for TAS1 and TAS2 loci, the rapidly evolving PPR-P clade yields the most easily recognizable set of RDR6/DCL4-dependent siRNAs in Arabidopsis. Most of the PPR-P clade transcripts are targeted by multiple miRNAs and/or tasiRNAs, which efficiently drive the cleavage fragments through the RDR6/DCL4 pathway. This very likely reinforces silencing that might be directed via targeting by miR161, miR400, and TAS1/TAS2 tasiRNAs. Given the abundance of the resulting clade-derived siRNAs, it is puzzling why siRNA biogenesis does not spread outside of the clade. Despite 33 nonclade PPR-P transcripts possessing multiple, plausible target sites for sequenced, clade-derived siRNAs, the RDR6/DCL4 system was not directed to this group of transcripts. It is possible that the target prediction algorithm, which we have used successfully with miRNAs and tasiRNAs (Allen et al., 2005; Fahlgren et al., 2007), did not accurately identify prospective targets for this analysis. It is also possible that clade-derived siRNAs had the potential to trigger silencing of nonclade PPR-P transcripts, but they were in cellular concentrations that were too low to function. Functional siRNAs may have also been compartmentalized away from nonclade PPR-P target transcripts.

What functions might targeting of the PPR-P clade by the RDR6/DCL4-dependent silencing system serve? This is not immediately obvious from genetic analysis, as most of the mutant phenotypes associated with loss of RDR6 and DCL4 appear to relate specifically to loss of TAS3 tasiRNAs (Xie et al., 2005; Adenot et al., 2006; Fahlgren et al., 2006; Garcia et al., 2006; Hunter et al., 2006). Of course, physiologic functions for silencing of the clade may lay beyond those that are revealed in short-term growth and development assays. The PPR gene family is widespread and has undergone a dramatic evolutionary expansion in plants (Small and Peeters, 2000; Rivals et al., 2006). An important observation is the close phylogenetic relationship between the radish (R. sativus) fertility restorer genes (PPR-A and PPR-B) and members of the targeted PPR-P clade in Arabidopsis (see Supplemental Figure 3 online; Desloire et al., 2003). The fertility restorer genes compensate for defects in mitochondrion-encoded ATPase alleles in other species (Kazama and Toriyama, 2003; Wen et al., 2003; Akagi et al., 2004; Wang et al., 2006), although in Arabidopsis there are no known mitochondrion-based cytoplasmic male sterility systems and, therefore, no corresponding fertility restorer functions identified. A group of Populus PPR-P genes also shows relatively close phylogenetic relationship to fertility restorer genes based on amino acid sequence similarity (data not shown). Interestingly, these genes have expanded in a clade, similar to what has occurred in Arabidopsis (see Supplemental Figure 3 online). Expansion of these PPR-P clades in several plant species may reflect mechanisms to diversify alleles for interactions with multiple mitochondrial gene alleles. As the functions of the PPR-P genes are not known in Arabidopsis or Populus, the putative allele interactions may not necessarily relate to fertility or be mitochondrion specific. Remarkably, the Populus PPR-P clade members are singly or dually targeted by two miRNAs, miR475 and miR476, indicating that this clade is under posttranscriptional regulation (Lu et al., 2005). While primary targeting of PPR transcripts in Arabidopsis and Populus has occurred independently, it is unclear if Populus PPR transcripts also spawn secondary siRNAs. Regardless, silencing of the expanding clades in multiple species may serve a dosage compensation function, especially if gene expansion is providing a reservoir of PPR-P alleles to balance infrequently dispersed mitochondrial alleles. Posttranscriptional silencing, as we proposed previously (Allen et al., 2004; Fahlgren et al., 2007), may enable rapid expansion of gene families while minimizing detrimental dosage effects.

Finally, while there is compelling evidence that many MIRNA genes had evolutionary origins in duplicated segments of protein-coding genes (Allen et al., 2004; Axtell et al., 2006; Fahlgren et al., 2007), the origins of TAS gene loci are not clear. However, in view of targeting of the PPR-P clade by miRNA and RDR6/DCL4 systems, and the resulting phased siRNA patterns, we speculate that TAS loci may originate from miRNA-targeted protein-coding genes. In fact, the abundant siRNA-generating PPR-P loci look very similar to known TAS loci, except that most still appear to encode functional proteins. However, one can easily imagine loss of protein-coding potential through genetic drift and/or recombination events, but maintenance of one or more functional siRNA modules in phase with initiator small RNA cleavage sites. A nascent TAS locus would be retained if suppression of gene-of-origin family members provided an advantage.

Small RNA Libraries

Small RNA libraries were described by Kasschau et al. (2007). Small RNA sequencing for wild-type (Col-0), dcl1-7, dcl2-1, dcl3-1, dcl4-2, rdr1-1, rdr2-1, and rdr6-15 inflorescence, Col-0 and rdr6-15 seedlings, and Col-0 leaf tissues was done by picoliter-scale pyrosequencing (454 Life Sciences; Margulies et al., 2005). All small RNA sequences are available for download at the Gene Expression Omnibus (GEO) (Edgar et al., 2002; Barrett et al., 2005; http://www.ncbi.nlm.nih.gov/geo/) under accession number GSE6682 and the Arabidopsis Small RNA Project Database (http://asrp.cgrb.oregonstate.edu/db/).

Arabidopsis thaliana Mutants and Target Validation

Mutant lines for dcl1-7, dcl2-1, dcl3-1, dcl4-2, rdr1-1, rdr2-1, rdr6-15, hen1-1, hst-15, and hyl1-2 were described previously (Park et al., 2002; Reinhart et al., 2002; Peragine et al., 2004; Vazquez et al., 2004a; Xie et al., 2004, 2005; Allen et al., 2005). Targets were validated using a 5′ RACE assay with the GeneRacer kit (Invitrogen) as described previously (Llave et al., 2002; Kasschau et al., 2003; Allen et al., 2004, 2005). Total RNA was isolated from inflorescence tissue of 28-d-old, stage 1 to 12 flowers of both Col-0 and hen1-1 plants. RNA was subjected to batch and spin column purification for poly(A)⁺ enrichment with an Oligotex mRNA Midi kit (Qiagen). mRNA was subsequently ligated to the GeneRacer RNA oligo adaptor without further modification. The GeneRacer oligo(dT) primer was used to prime cDNA synthesis via reverse transcription. Two rounds of 5′ RACE reactions were performed with two gene-specific primers designed ∼200 to 500 nucleotides from the 3′ end of the predicted target sites using PCR parameters defined by the GeneRacer kit and an annealing temperature of 65°C. The 5′ RACE products were gel purified, TOPO TA cloned (Invitrogen), and sequenced.

Alignments, Phase Transformations, and Expression Profiling

Alignments between TAS3a, TAS3b, and TAS3c were performed using T-Coffee (Notredame et al., 2000). Phase transform data were produced using the following equation:

where n = number of phase cycle positions occupied by at least one small RNA read within an eight-cycle window, and k = the total number of reads for all small RNAs with consolidated start coordinates in a given phase within an eight-cycle window.

Phase cycle length was set at 21 nucleotides (see Figure 3). All microarray data were generated using Affymetrix ATH1 arrays and are available at GEO (Edgar et al., 2002; Barrett et al., 2005; http://www.ncbi.nlm.nih.gov/geo/).

Phylogenetic Analysis

PPR amino acid sequences were aligned using ClustalW (Higgins et al., 1996) with either the BLOSUM or the GONNET matrix series and with a gap opening penalty of 10 for both the pairwise and multiple alignments and a gap extension penalty of 0.1 and 0.2 for pairwise and multiple alignments, respectively. Gaps in the multiple alignments, defined as regions with at least five consecutive gap characters in all but one or two sequences, were removed. In addition, the first and last 80 positions of the multiple alignments were deleted manually to restrict the analysis to the most conserved region. Trees were produced with ClustalW using the neighbor-joining method (Saitou and Nei, 1987) and 1000 bootstrap replicates. Trees were viewed and rooted using TREEVIEW (Page, 1996).

Accession Numbers

The cDNA sequence for TAS3b was deposited in The Arabidopsis Information Resource database under accession number At5g49615. Sequence data from this article can be found in the EMBL/GenBank data libraries under accession number GSE6682 (GEO small RNA data). Col-0 (control for dcl1-7, rdr6-15, and dcl4-2), dcl1-7, rdr6-15, and dcl4-2 were from experiments described by Xie et al. (2005) (GEO accession number GSE3011, samples GSM65929, GSM65930, GSM65931, GSM65932, GSM65933, GSM65934, GSM65935, GSM65936, GSM65937, GSM65938, GSM65939, GSM65940, GSM65941, GSM65942, and GSM65943). Col-0 (control for hst-15, hyl1-2, dcl2-1, dcl3-1, rdr1-1, and rdr2-1), Ler (control for hen1-1), hst-15, hyl1-2, dcl2-1, dcl3-1, rdr1-1, and rdr2-1 were from experiments described previously (Allen et al., 2005) (GEO accession number GSE2473, samples GSM47011, GSM47012, GSM47013, GSM47014, GSM47015, GSM47016, GSM47020, GSM47021, GSM47022, GSM47028, GSM47029, GSM47030, GSM47031, GSM47032, GSM47033, GSM47034, GSM47035, GSM47036, GSM47037, GSM47038, GSM47039, GSM47040, GSM47041, GSM47042, GSM47043, GSM47044, GSM47045, GSM47046, GSM47047, and GSM47048).

Supplemental Data

The following materials are available in the online version of this article.

• Supplemental Figure 1. Mapping of Transcription Start Sites of TAS3 Family Genes by RNA Ligase-Mediated 5′ RACE.

• Supplemental Figure 2. Mapping of miR161- and TAS2 3′D6(−)-Mediated Cleavage Sites in PPR Gene Transcripts by 5′ RACE.

• Supplemental Figure 3. Phylogenetic Tree of PPR-P Clade Proteins from Arabidopsis, PPR-P–Related Proteins from Populus and R. sativus, and Outgroup PPR proteins CRP1 (Z. mays), PET309 (Saccharomyces cerevisiae), PTCD1 (Homo sapiens), and RF1A and RF1B (Oryza sativa).

• Supplemental Table 1. RDR6-Dependent 21-Nucleotide RNA-Generating Loci.