OPEN ACCESS ARTICLE

This Article Free via Open Access: OA OA Abstract Full Text (PDF) Supplemental Data OA All Versions of this Article: 19/6/1750    most recent tpc.107.051706v1 Alert me when this article is cited Alert me if a correction is posted Services Email this article to a friend Related articles in Plant Cell Similar articles in this journal Similar articles in ISI Web of Science Similar articles in PubMed Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via HighWire Citing Articles via ISI Web of Science (5) Citing Articles via Google Scholar Google Scholar Articles by Axtell, M. J. Articles by Bartel, D. P. Search for Related Content PubMed PubMed Citation Articles by Axtell, M. J. Articles by Bartel, D. P. Agricola Articles by Axtell, M. J. Articles by Bartel, D. P.

Common Functions for Diverse Small RNAs of Land Plants^[W],[OA]

^a Department of Biology and Huck Institutes of the Life Sciences, Pennsylvania State University, University Park, Pennsylvania 16802
^b Whitehead Institute, Howard Hughes Medical Institute, Cambridge, Massachusetts 02142
^c Department of Biology, Massachusetts Institute of Technology, Cambridge, Massachusetts 02142

¹ To whom correspondence should be addressed. E-mail mja18{at}psu.edu.

	ABSTRACT

TOP ABSTRACT INTRODUCTION RESULTS DISCUSSION METHODS REFERENCES

 
Endogenous small RNAs, including microRNAs (miRNAs) and short interfering RNAs (siRNAs), are critical components of plant gene regulation. Some abundant miRNAs involved in developmental control are conserved between anciently diverged plants, while many other less-abundant miRNAs appear to have recently emerged in the Arabidopsis thaliana lineage. Using large-scale sequencing of small RNAs, we extended the known diversity of miRNAs in basal plants to include 88 confidently annotated miRNA families in the moss Physcomitrella patens and 44 in the lycopod Selaginella moellendorffii. Cleavage of 29 targets directed by 14 distinct P. patens miRNA families and a trans-acting siRNA (ta-siRNA) was experimentally confirmed. Despite a core set of 12 miRNA families also expressed in angiosperms, weakly expressed and apparently lineage-specific miRNAs accounted for the majority of miRNA diversity in both species. Nevertheless, the molecular functions of several of these lineage-specific small RNAs matched those of angiosperms, despite dissimilarities in the small RNA sequences themselves, including small RNAs that mediated negative feedback regulation of the miRNA pathway and miR390-dependent ta-siRNAs that guided the cleavage of AUXIN RESPONSE FACTOR mRNAs. Diverse, lineage-specific, small RNAs can therefore perform common biological functions in plants.

	INTRODUCTION

TOP ABSTRACT INTRODUCTION RESULTS DISCUSSION METHODS REFERENCES

 
Most eukaryotic organisms are capable of producing small RNAs, typically between 21 and 24 nucleotides in length, that serve as sequence-specific posttranscriptional regulators. Many short interfering RNAs (siRNAs) are initially produced as short duplexes containing two-nucleotide, 3' overhangs via the enzymatic activity of endonucleases of the Dicer family acting upon long double-stranded RNA (dsRNA) (Baulcombe, 2005). One of the two siRNA strands is stably incorporated into an ARGONAUTE (AGO) protein that then functions to repress the expression of target RNAs, which are selected via base-pairing interactions with the siRNA (Hammond, 2005). This repression can take several forms, including AGO-catalyzed target cleavage (Liu et al., 2004; Meister et al., 2004) and translational repression (Doench et al., 2003). Some endogenous siRNAs also direct DNA and histone modifications to homologous genomic loci (Zilberman et al., 2003). Such siRNAs are especially numerous and diverse in flowering plants (Lu et al., 2005; Rajagopalan et al., 2006; Kasschau et al., 2007).

MicroRNAs (miRNAs) are a distinct type of small RNA that also function as negative regulators of gene expression. Primary miRNA transcripts (pri-miRNAs) adopt stem-loop RNA secondary structures from which a specific 21-nucleotide duplex is excised by a Dicer endonuclease (Bartel, 2004). Similar to siRNA duplexes, the resulting miRNA/miRNA* duplex possesses characteristic two-nucleotide 3' overhangs and is unwound coincident with the association of the mature, single-stranded miRNA with an AGO protein. The mature miRNA then serves to guide the bound AGO protein to target mRNAs based upon complementarity. In both plants and animals, the regulation of target gene expression via miRNAs is critical for numerous biological processes (Wienholds and Plasterk, 2005; Jones-Rhoades et al., 2006; Mallory and Vaucheret, 2006). In plants, miRNA-directed target cleavage can sometimes stimulate the production of dsRNA using the cleaved product as a template; the resulting dsRNA is subsequently processed by a Dicer protein to produce siRNAs referred to as trans-acting siRNAs (ta-siRNAs) because of their miRNA-like ability to repress mRNAs from loci distinct from the originating locus (Peragine et al., 2004; Vazquez et al., 2004; Allen et al., 2005; Yoshikawa et al., 2005; Axtell et al., 2006).

Although the miRNAs characterized in angiosperms have numerous distinct biological functions, there is a striking tendency for many of them to regulate developmental processes (Jones-Rhoades et al., 2006). Previous studies have shown that some miRNA–target interactions first identified in Arabidopsis thaliana are deeply conserved among land plants (Floyd and Bowman, 2004; Arazi et al., 2005; Axtell and Bartel, 2005). Interestingly, the ancient miRNA regulatory interactions are disproportionately those thought to be involved in development, suggesting that these units of posttranscriptional gene control have been indispensable during the diversification of land plants.

The sequencing of a moss (Physcomitrella patens) genome and a lycopod (Selaginella moellendorffii) genome provides an opportunity to explore the diversity of small RNA functions in nonseed plants and to more rigorously address the trajectories of small RNA evolution in three widely separated plant lineages. To date, >40 distinct miRNA families have been annotated from P. patens based upon analysis of sequenced small RNAs or bioinformatic predictions (Arazi et al., 2005; Talmor-Neiman et al., 2006a; Fattash et al., 2007). In lycopods, the expression of miR166 has been confirmed in S. moellendorffii (Floyd and Bowman, 2004), and the expression of two additional families proposed based upon microarray experiments using RNA from Selaginella uncinata (Axtell and Bartel, 2005).

In this study, we applied high-throughput small RNA sequencing to the discovery and analysis of miRNAs from the moss P. patens and the lycopod S. moellendorffii. Our results directly demonstrated the expression of a much large number of distinct miRNA families in P. patens and S. moellendorffii then has previously been reported. Only a few of these miRNAs were expressed in common between mosses, lycopods, and angiosperms; these core conserved plant miRNAs were abundantly expressed in both angiosperms and basal pants, and all, save one, appeared to be involved on the regulation of developmental processes. The majority of miRNA diversity in both lycopods and mosses appeared lineage specific and generally expressed at lower levels than were the core conserved miRNAs. The experimental isolation of 29 cleavage products directed by 14 distinct P. patens miRNAs and a ta-siRNA demonstrated that the core conserved plant miRNAs have generally retained homologous target interactions during the diversification of land plants. Surprisingly, experimentally validated targets also included several instances in which lineage-specific P. patens miRNAs and ta-siRNAs shared common functions with sequence-distinct Arabidopsis small RNAs.

	RESULTS

TOP ABSTRACT INTRODUCTION RESULTS DISCUSSION METHODS REFERENCES

 
Identification of miRNAs from P. patens and S. moellendorffii
Direct sequencing of small RNAs and computational predictions have resulted in the annotation of 44 distinct P. patens miRNA families (Arazi et al., 2005; Talmor-Neiman et al., 2006a; Fattash et al., 2007). We have recently reported a data set of 127,135 unique P. patens small RNAs from three libraries, represented by a total of 561,102 reads, which provided insights into siRNA biogenesis in plants (Axtell et al., 2006). Because of the >100-fold increase in sequencing depth over previous P. patens small RNA sequencing efforts, we reasoned that additional less-abundant miRNAs might be found by analysis of these data. In addition, a single small RNA library prepared from the above-ground tissues of S. moellendorffii was sequenced to yield 36,582 unique small RNAs, represented by 149,586 reads that matched the genome (whole-genome shotgun traces).

To produce a list of possible miRNA loci, consensus genomic sequences surrounding each small RNA were derived by assembly of the relevant whole-genome shotgun (WGS) traces, followed by RNA secondary structure prediction. As a preliminary screen, the minimum free energy structure for each locus was analyzed using MIRcheck (Jones-Rhoades and Bartel, 2004), which evaluates paired regions of a predicted RNA secondary structure for characteristic features of miRNA stem loops. The candidate miRNA loci that passed MIRcheck were then analyzed for evidence of miRNA* accumulation. A miRNA/miRNA* duplex, with two-nucleotide 3' overhangs indicative of DICER-LIKE1 (DCL1) cleavage, would be expected if a moss DCL1 ortholog excised the miRNA from the hairpin precursor as indicated by a previous study (Talmor-Neiman et al., 2006a). Indeed, the sequencing of both a miRNA and a miRNA* constitutes evidence of DCL1-like processing from a hairpin precursor, which is evidence that can supplement, or even substitute for, other criteria for confidently distinguishing miRNAs from endogenous siRNAs and other RNAs that should not be annotated as miRNAs (Rajagopalan et al., 2006; Ruby et al., 2006; Fahlgren et al., 2007).

In total, 88 distinct P. patens miRNA families, arising from 205 miRNA genes, were recovered by this analysis (Tables 1 and 3; see Supplemental Table 1 online). This set included all previously recognized P. patens miRNA families except for 13 of the 22 families uniquely reported by Fattash et al. (2007). Analysis of the deep sequencing data from P. patens also suggested refinements to a small number of previous miRNA annotations; these suggestions included reversals of miRNA/miRNA* designations, changes in the annotated positions of the dominant mature miRNA within stem-loop precursors, and consolidations of redundantly reported families (see Supplemental Table 1 online). Our resulting roster of confidently identified P. patens miRNAs (see Supplemental Table 1 online) included 12 families conserved with angiosperms and/or lycopods (Table 3). An additional six families of P. patens miRNAs annotated on the basis of computational comparisons to angiosperm miRNAs (miR167, miR172, miR395, miR414, miR418, and miR419; Fattash et al., 2007) were not observed in our sequence libraries. The expression of the remaining 76 sequence-distinct families was, to the best of our current knowledge, bryophyte specific. These included 57 novel P. patens miRNA families that have not previously been reported (Table 1).

P. patens miRNA Loci Are Often Clustered and Frequently Overlap with Protein-Coding Genes
The availability of a draft assembly and transcript annotations of the P. patens genome allowed the detailed examination of miRNA loci. Most miRNA loci characterized in angiosperms appear to arise from independent transcription units, each encoding a single miRNA-producing stem loop that does not overlap with annotated protein-coding loci (Jones-Rhoades et al., 2006). Four P. patens miR1219 loci are contained in two clusters, strongly suggesting that they are expressed from two polycistronic precursors (Talmor-Neiman et al., 2006a). We found that clustered miRNA expression was common for P. patens miRNAs (Figure 1A ; see Supplemental Table 1 online). A total of 48 of the 205 P. patens miRNA loci were arranged in closely linked clusters of two or three miRNA stem loops. Clustered miRNA stem loops were observed for both ancient miRNAs, such as miR160 and miR166, and for novel bryophyte-specific miRNAs (Figure 1A). With only a few exceptions, these clusters contained paralogs encoding mature miRNA families in the same family, rather than encoding distinct mature miRNAs unrelated in sequence as is frequently the case for animal miRNA clusters (see Supplemental Table 1 online). Similarly, recent results in maize (Zea mays) have demonstrated the expression of a polycistronic miRNA precursor for miR156 (Chuck et al., 2007), and a miRNA from the green alga C. reinhardtii is also arranged in a tandem cluster (Zhao et al., 2007).

View larger version (17K): [in this window] [in a new window] [as a PowerPoint slide]   Figure 1. Properties of miRNA Expression in P. patens. (A) Schematic depictions of representative clustered P. patens miRNAs that are potentially processed from polycistronic pri-miRNAs, indicated by black lines. In total, 48 out of the 205 P. patens miRNA loci were contained in 21 clusters containing two or three miRNA stem loops each (see Supplemental Table 1 online). Distances refer to the number of nucleotides separating the regions annotated as miRNA foldbacks. (B) Mapping of mature P. patens miRNAs with respect to annotated protein-coding genes. The numbers of mature miRNAs that overlapped various annotations of the P. patens Phypa1_1 genome assembly are indicated.

Some lineage-specific miRNA hairpins expressed in Arabidopsis possess extended complementarity to protein-coding transcripts that are frequently, though not always, the targets of such miRNAs (Allen et al., 2004; Rajagopalan et al., 2006; Fahlgren et al., 2007). This suggests that the initiation of novel miRNA–target interactions at least sometimes involves the inverted duplication of transcribed regions. Similar analyses for P. patens miRNA loci were complicated by the fact that more than half of the miRNA hairpins overlapped annotated protein-coding regions; thus, any extended homology found between such hairpins and other protein-coding transcripts would be more likely due to conserved protein domains rather than indicative of the origin of the miRNA locus. Nonetheless, of the 87 loci that did not overlap with annotated protein-coding loci, five were found to have extended complementarity with annotated transcripts (FASTA search E < 0.05): miR1026a, miR1072, miR1073, miR1220a, and miR1220b. Interestingly, none of the presumed originating transcripts of these five loci were also predicted as targets of their mature miRNAs, suggesting that these loci might have been selected to have targets distinct from their originating gene families as appears to have happened for a few Arabidopsis miRNAs (Fahlgren et al., 2007). However, the majority of P. patens miRNA loci that don't themselves overlap with annotated protein-coding transcripts did not appear to have extended complementarity to annotated transcripts.

P. patens miRNA Targets
Several targets of P. patens miRNAs and ta-siRNAs have been predicted and/or validated using EST data, WGS trace data, or directed cDNA cloning (Floyd and Bowman, 2004; Arazi et al., 2005; Axtell et al., 2006; Floyd et al., 2006; Talmor-Neiman et al., 2006a, 2006b; Fattash et al., 2007). The availability of a large set of annotated mRNAs derived from the draft genome assembly coupled with the expanded number of P. patens miRNAs discovered in this study prompted us to revisit miRNA and ta-siRNA target prediction. Targets were predicted using the methodology described by Allen et al. (2005), except that conservation of targeting in other plant genomes was not considered (see Supplemental Table 3 online). This methodology performed reasonably well, as assessed by a signal-to-noise ratio of 2.78 to 1 when compared with a set of shuffled miRNA sequences controlled for di- and trinucleotide composition.

Computational target predictions were based upon the full contingent of annotated transcripts that accompanied the Phypa1_1 genome release and thus are expected to be more comprehensive than those based upon EST or WGS sequence data. However, we recognized that this database of P. patens transcripts was probably not as accurate as those derived from the current annotations of the more thoroughly studied Arabidopsis and O. sativa genomes; thus, failure to annotate P. patens transcripts could have led to false negatives (missed predictions), while the presence of inaccurately annotated transcripts and inability to consider conservation in a related genome could have led to false positive predictions. We therefore used RLM 5'-RACE (Llave et al., 2002a) to test many of these predictions by assaying for the existence of cDNAs corresponding to miRNA-mediated cleavage products. A total of 29 targets of 14 miRNA families and a ta-siRNA were validated (Figure 2 ). As expected based upon previous observations (Floyd and Bowman, 2004; Arazi et al., 2005; Axtell and Bartel, 2005; Floyd et al., 2006; Fattash et al., 2007), most of the miRNAs conserved between Arabidopsis and P. patens had similar targets in both organisms (Figure 2A; see Supplemental Table 3 online), supporting the hypothesis that these trans-regulatory interactions have remained constant during the diversification of land plants. For instance, miR156 directed the cleavage of two SBP box genes (including Pp SBP3, which had previously been demonstrated to be a miR156 target; Arazi et al., 2005), miR166 directed the cleavage of the five class III HD-ZIP genes that had been previously been observed to possess miR166 complementary sites (Floyd et al., 2006), miR171 directed the cleavage of two GRAS domain transcription factor genes, and miR408 directed the cleavage of a plastocyanin domain–containing gene. In flowering plants, the miR159/319 family targets subsets of the MYB and TCP transcription factor gene families (Rhoades et al., 2002; Palatnik et al., 2003); none of the initially predicted P. patens miR319 targets had recognizable homology to MYB or TCP genes (see Supplemental Table 3 online). However, a gene with a cyclin domain was among the initially predicted P. patens miR319 targets and was indeed cleaved in vivo at the predicted miR319 complementary site (Figure 2A). Reducing the stringency of the search for P. patens miR319 targets allowed the identification of two homologs of MYB transcription factors that were similar to the angiosperm targets of the sequence-related miR159 (see Supplemental Table 3 online); RLM 5'-RACE demonstrated that these two genes were also cleaved by a miR319-guided activity in P. patens (Figure 2A). Thus, the initial search parameters were too stringent to identify all bona fide cleavage targets. This sensitivity trade-off when predicting plant miRNA targets with acceptable specificity without recourse to evolutionary conservation was expected. For example, the target prediction score cutoffs would have been too stringent to identify the 3' miR390 complementary site of Pp TAS3a and also would have failed to identify the 5' miR390 complementary sites of angiosperm TAS3 loci, even though these have been confirmed by direct experimentation to be functional (Axtell et al., 2006; Talmor-Neiman et al., 2006b).

View larger version (18K): [in this window] [in a new window] [as a PowerPoint slide]   Figure 2. Validated P. patens miRNA and ta-siRNA Targets. (A) Cleaved targets of conserved P. patens miRNAs. Schematics of target mRNAs along with the relative positions of the 5' residues of all sequenced uncapped cDNAs are shown. The positions of the miRNA complementary sites are indicated with red vertical lines. Gray boxes indicate the positions of regions coding for conserved protein domains, whose abbreviated names are listed below. Boxes correspond to open reading frames and horizontal black lines to untranslated regions. Shaded pixels indicate the positions of 5' residues of individual cDNA clones, which are shaded as indicated in the key below to indicate their proximity to the site predicted by miRNA-mediated cleavage between the tenth and eleventh nucleotide of complementarity. An asterisk indicates that Pp SBP3 was included as a positive control (Arazi et al., 2005). Alignments of target sites with small RNAs are provided in Supplemental Table 3 online. nt, nucleotides. (B) Cleaved targets of nonconserved P. patens miRNAs as in (A). Hatch marks in Phypa1_1 109598 indicate a region where the splicing of the observed cDNA fragments differed from the annotation.

Recurrence of miRNA-Mediated Negative Feedback Loops in Plants
The gene encoding the DCL protein chiefly responsible for miRNA biogenesis in Arabidopsis, DCL1, is itself regulated by miRNAs: miR162 is expressed from two distinct loci and can target DCL1 mRNA via a complementary site formed by the splicing of exon 12 to exon 13 (Xie et al., 2003), while miR838 is excised from a stem loop present in intron 14 of the DCL1 pre-mRNA (Rajagopalan et al., 2006). Target cleavage mediated by miR162 and excision of the miR838/miR838* duplex, both of which would lead to repression of DCL1 mRNA expression, have been hypothesized to require the activity of the DCL1 protein and thereby constitute negative feedback circuits. We found that the P. patens MIR1047 locus resided in the sense orientation relative to intron seven of a gene whose predicted protein product had strong overall similarities to DCL proteins (Figure 3A ). RLM 5'-RACE experiments detected RNA fragments terminating at positions consistent with the excision of the miR1047/miR1047* duplex from this pre-mRNA (Figure 3B). Evidence for miR1047 excision from pre-mRNA in which intron seven was not spliced was found from both poly(A)⁺ and poly(A)^– RNA samples. We also noted that approximately half of these truncated pre-mRNAs terminated at the loop-proximal cleavage position, indicating that the first DCL-mediated cleavage of the stem loop must have occurred at the loop-proximal position (Figure 3B). This was in contrast with the processing of animal pri-miRNAs (Lee et al., 2003) and certain other plant pri-miRNAs (Kurihara and Watanabe, 2004; Vaucheret et al., 2006), for which the first cleavage events can occur proximal to the base of the hairpins. Nevertheless, similar results for Arabidopsis miR838 processing (Rajagopalan et al., 2006) support the conclusion that the loop-proximal cleavage event can sometimes occur first during the biogenesis of some plant miRNAs.

View larger version (27K): [in this window] [in a new window] [as a PowerPoint slide]   Figure 3. Analogous miRNA-Mediated Feedback Loops Despite Divergent Small RNA Sequences. (A) Schematic of the Pp DCL1a pre-mRNA, which contains an expressed miRNA within intron seven. Thick rectangles represent predicted coding exons, thin rectangles represent untranslated regions, and thin lines represent predicted introns. Colored exonic regions correspond to conserved protein domains as indicated. Gray indicates the location of the MIR1047 hairpin within intron seven. (B) Cleavage at miR1047 in intron seven of the Pp DCL1a pre-mRNA. The predicted secondary structure of MIR1047 is shown, with the miRNA/miRNA* duplex highlighted in red. Fractions indicate the number of 5'-RACE products observed to terminate at the indicated positions over the total number of 5'-RACE clones sequenced from three different experiments from the indicated combinations of RNA samples and RT primers. (C) An unrooted phylogenetic reconstruction of relationships among Arabidopsis, O. sativa, and P. patens DCL proteins using C. reinhardtii DCL1 as an outgroup. The known miRNA-mediated feedback loops between At DCL1 pre-mRNA, mRNA, protein, and DCL1-dependent miRNAs are indicated in blue. The hypothesized feedback loop between Pp DCL1a pre-mRNA, protein, and ppt-MIR1047 is indicated in red. Nodes with <100% bootstrap support are noted with the percentage of supporting replicates. Accession numbers are given in Methods. (D) An unrooted phylogenetic reconstruction of relationships among Arabidopsis, O. sativa, and P. patens AGO proteins using C. reinhardtii AGO1 and AGO2 as outgroups as displayed in (B). miRNA-mediated control of Arabidopsis AGO genes is shown in blue, and miRNA-mediated control of P. patens AGO genes by ppt-miR904 is shown in red. Accession numbers are given in Methods.

The founding member of the AGO gene family, Arabidopsis AGO1, encodes a miRNA-guided endonuclease required for the trans-regulatory effects of many miRNAs (Vaucheret et al., 2004; Baumberger and Baulcombe, 2005; Qi et al., 2005). AGO1 is required for the stability and function of miR168, which targets the AGO1 mRNA; this homeostatic feedback loop is critical for proper Arabidopsis development and the normal functioning of many other Arabidopsis miRNAs (Vaucheret et al., 2004, 2006). We found that the three validated targets of the P. patens–specific miR904 were homologous to AGO genes (see Supplemental Table 3 online; Figures 2B and Figure 3D). The three P. patens AGO proteins whose mRNAs were validated targets of miR904 were clearly embedded within the AGO1 clade of land plant AGO proteins (Figure 3D), suggesting that one or more of these miR904 targets might be the major miRNA-guided endonucleases in P. patens; to reflect this similarity, we named these genes Pp AGO1a (Phypa1_1 205541), Pp AGO1b (Phypa1_1 158832), and Pp AGO1c (Phypa1_1 141045). Despite having homologous targets, no similarity between the sequences of miR168 and miR904 was detected. This suggested that the miRNA-AGO feedback loops either appeared independently in angiosperms and bryophytes by convergent evolution or that the sequences of the miRNA and AGO target complementary sites have diverged from those of a common ancestor. Both miR168 and miR904 are released from the 5' arms of their respective precursor RNAs, and their cleavage sites in both Arabidopsis and P. patens AGO mRNAs lie just upstream of the region encoding the conserved PAZ domain. These observations leave open the possibility that the miR168-AGO and miR904-AGO regulatory interactions could have descended from a common ancestor. The observation that negative feedback interactions between miRNAs, DCL genes, and AGO genes exist in both bryophytes and angiosperms suggests that these interactions are critical for control of miRNA functions in both lineages.

The phylogenetic comparison of the DCL and AGO families in Arabidopsis, O. sativa, P. patens, and C. reinhardtii provided the opportunity to consider the evolution of these families during the divergence of land plants. The DCL family appeared to have diverged and specialized early, with at least three members present in the last common ancestor of these land plants and rather modest duplication since then. We failed to observe a DCL2 homolog in P. patens, though the node supporting the placement of Pp DCL4 in the DCL4 clade was weak (Figure 3C). Whether this reflected the first appearance of DCL2 in the seed plant lineage after the divergence of the bryophytes, selective loss of DCL2 in the bryophyte lineage, or simply a failure to annotate the Pp DCL2 transcript in the Phypa1_1 draft genome was unclear. Because the genome of S. moellendorffii has not yet been assembled and annotated, we were unable to include lycopod DCL proteins in this analysis. Evolution of the AGO family appears to have been more dynamic, with much more expansion, and perhaps loss, of paralogs since the last common ancestor. None of the six identified P. patens AGO proteins fell in the AGO2/3/7-like clade (Figure 3D); whether this reflected selective loss in the bryophyte lineage, diversification of the AGO2/3/7 clade of AGO proteins in the seed plant lineage after divergence from the bryophytes, or a failure to annotate some P. patens AGO loci in the Phypa1_1 assembly was unclear. It is also possible that some of the apparent duplications and extinctions in the AGO genes were artifacts derived from errors in phylogenetic reconstruction; as for the DCL proteins, the S. moellendorffii AGO proteins were not available to assist in resolving this issue.

Pp TAS3a-d Targets Include Genes Similar to Arabidopsis ETT/ARF3 and ARF4
We previously described four ta-siRNA loci targeted by miR390 in P. patens, which we referred to as Pp TAS1-4 (Axtell et al., 2006). The existence of what we called Pp TAS1 was first suggested by Arazi et al. (2005) and was functionally characterized by Talmor-Neiman et al. (2006b), who named it TAS4. We suggest that these four P. patens ta-siRNA loci would be better referred to as Pp TAS3a-d for three reasons: All four loci possess dual miR390 complementary sites, similar to angiosperm TAS3 loci (Axtell et al., 2006; Howell et al., 2007), an unrelated Arabidopsis ta-siRNA locus has also recently been named TAS4 (Rajagopalan et al., 2006), and targets identified and validated in this study (see below) demonstrate that ta-siRNAs derived from these moss loci function to regulate AUXIN RESPONSE FACTOR (ARF) genes similar to those regulated by the angiosperm TAS3 loci.

Pp TAS3a produces abundant phased siRNAs that are dependent upon the RNA-dependent RNA polymerase Pp RDR6 (Talmor-Neiman et al., 2006b). Interestingly, Pp rdr6 mutants also display accelerated differentiation of leafy gametophores from protonemata, suggesting that the Pp TAS3 loci are important for the timing of gametophyte development (Talmor-Neiman et al., 2006b); in Arabidopsis, rdr6 mutants display an At TAS3–dependent acceleration of sporophytic vegetative phase change (Peragine et al., 2004; Adenot et al., 2006; Fahlgren et al., 2006; Garcia et al., 2006). Alignment of angiosperm TAS3 genes from multiple species allowed the identification of two conserved, nearly identical ta-siRNAs that function in Arabidopsis to regulate ETT/ARF3 and ARF4 (Allen et al., 2005; Fahlgren et al., 2006; Hunter et al., 2006), but these ta-siRNAs did not have any recognizable homology to the Pp TAS3a-d loci (Axtell et al., 2006).

Given their similar complementarity to miR390, we reasoned that Pp TAS3a-d might be paralogs derived from duplication and divergence of an ancestral locus, and perhaps the areas producing biologically relevant ta-siRNAs may have been preferentially conserved. This appears to have been the case, as alignment of the ta-siRNA–producing regions of Pp TAS3a-d revealed two 21-nucleotide blocks of strong conservation (Figure 4A ). Sense ta-siRNAs whose 5' residues were coincident with the beginning of block 1 were two to three nucleotides offset from the phasing register set by cleavage at the 3' miR390 complementary site, while antisense ta-siRNAs whose 5' ends were coincident with the beginning of block 2 were exactly coincident with the phasing register set by the 3' miR390 complementary site (Figure 4A). We found no evidence for miR390 expression in S. moellendorffii, suggesting that the miR390-TAS3 cascade of small RNA production might have been lost in that lineage.

View larger version (37K): [in this window] [in a new window] [as a PowerPoint slide]   Figure 4. Similar ARF Targets for Diverse ta-siRNAs and miRNAs. (A) T-COFFEE alignment of the regions of Pp TAS3a-d between the two miR390 complementary sites. Two 21-nucleotide blocks conserved between the four paralogs are indicated. (B) Alignments of the block 1 (+) siRNAs and the block 2 (-) RNAs with their predicted and validated targets. Lowercase letters indicate positions that cannot form either Watson-Crick or G-U pairs with the target. Arrows indicate cleavage-validated target sites. The asterisk indicates the cleavage site of the AP2 domain–containing gene Phypa1_1 129196 by Pp TAS3a-d block 1 (+)–derived siRNAs as demonstrated by Talmor-Neiman et al. (2006b). (C) An unrooted phylogenetic reconstruction of the relationships between Arabidopsis and P. patens ARF proteins. Nodes with <100% bootstrap support are noted with the percentage of supporting replicates. ARF genes known or predicted to be under small RNA–mediated regulation are indicated. Accession numbers are given in Methods.

The block 2(-) ta-siRNAs were predicted to target four distinct transcripts (see Supplemental Table 3 online), two of which encoded proteins homologous to angiosperm ARFs (Figures 4B and 4C). Cleavage products were detected for both ARF messages, thereby demonstrating that these P. patens ARF mRNAs were indeed regulated by Pp TAS3–derived ta-siRNAs (Figure 2B). Unexpectedly, these RLM 5'-RACE experiments revealed a second, longer species of uncapped mRNA for both of the P. patens ARF genes, leading us to investigate whether the longer products also corresponded to small RNA-mediated cleavage products. The putative upstream small RNA complementary sites did not appear to correspond to any ta-siRNA produced from Pp TAS3a-d nor were these two ARF genes predicted as targets of any P. patens miRNAs using established protocols (Allen et al., 2005; Rajagopalan et al., 2006). However, relaxing the stringency of the target prediction parameters revealed that the upstream cleavage sites corresponded with the bryophyte-specific miR1219 such that the observed cleavage product termini aligned exactly with the tenth nucleotide of complementarity (Figure 2B). We concluded that the two P. patens ARF genes were targeted at distinct sites by miR390-initiated, Pp TAS3a-d–derived ta-siRNAs and by the bryophyte-specific miR1219. We did not sequence any small RNAs that corresponded to the regions of the two P. patens ARF genes between the Pp TAS3a-d block 2(-) and miR1219 cleavage sites.

Phylogenetic comparison of P. patens ARF proteins to Arabidopsis ARF proteins revealed that the Pp TAS3a-d block 2(-) and miR1219-targeted ARFs were closely related to At ETT/ARF3 and At ARF4 (Figure 4C). At ETT/ARF3 and At ARF4 are targeted by miR390-phased, At TAS3–derived ta-siRNAs at two distinct complementary sites located downstream of the region coding for the ARF domain (Allen et al., 2005). Our data thus demonstrated that the arrangement of dual small RNA target sites within ARF mRNAs was shared between bryophytes and angiosperms despite divergence of both the sequence composition and of the biogenesis of the small RNAs themselves, providing a striking example of common regulatory functions being fulfilled by lineage-specific small RNAs. Moreover, the Pp TAS3a-d genes appear to transmute the expression of one miRNA (miR390) into two biologically relevant sets of related small RNAs that are capable of the coordinated regulation of a diverse set of regulatory transcripts. One set of Pp TAS3a-d ta-siRNAs regulates AP2 domain genes in concert with miR529(5'), while another set acts in concert with a lineage-specific miRNA, miR1219, to regulate the expression of ETT/ARF3 and ARF4-like genes.

Other Small RNAs in Basal Plants
The entire set of 561,102 P. patens small RNA reads that matched at least one WGS trace (Axtell et al., 2006) were categorized based upon genomic and length characteristics (Figure 5A ). miRNAs and ta-siRNAs accounted for 28 and 3% of the reads, respectively. Both subpopulations were dominated by 21-nucleotide sequences, with a considerable minority of 22-nucleotide sequences deriving from the Pp TAS3a-d ta-siRNAs (Figure 5A); these properties were consistent with those of angiosperm miRNAs and ta-siRNAs (Reinhart et al., 2002; Rajagopalan et al., 2006). Strikingly, 47% of the reads that initially matched at least one raw WGS trace did not match anywhere along the Phypa1_1 draft genome assembly. This subpopulation of reads included many short sequences of sizes inconsistent with the activities of any known Dicer endonuclease and some small RNAs between 20 and 24 nucleotides (Figure 5A). A substantial number of these reads were only sequenced once and disproportionately tended to match only a single WGS trace (data not shown), a statistically unlikely occurrence considering that shotgun coverage of the P. patens genome is 8x. This suggested that this subpopulation of reads contains some convergent sequencing errors (i.e., a sequencing error in a single small RNA and a single WGS trace that causes a spurious match). It might also contain breakdown products of RNAs and Dicer-derived small RNAs originating from genomic regions that are difficult to sample and assemble using shotgun sequencing and thus were left out of the draft genome assembly.

View larger version (32K): [in this window] [in a new window] [as a PowerPoint slide]   Figure 5. Small RNA Populations in P. patens and S. moellendorffii. (A) Categorization of 561,102 P. patens small RNA reads that matched at least one WGS trace. Reads were categorized as miRNAs (matching one or more of the 205 annotated miRNA hairpins), ta-siRNAs (from PpTAS3a-d), unknown (corresponding to 10 or fewer loci in the Phypa1_1 genome assembly), unknown-repetitive (corresponding to 11 or more loci in the Phypa1_1 genome assembly), or those not matching any loci in the Phypa1_1 genome assembly. Histograms display the length distributions of each of these small RNA populations. Reads that corresponded to the nuclear rRNAs and the chloroplast genome were not included. (B) Categorization of 149,586 S. moellendorffii small RNA reads that matched at least one WGS trace. Reads are categorized as miRNAs (matching one or more of the 58 annotated miRNA hairpins), unknown (matching 100 or fewer WGS traces), or unknown-repetitive (corresponding to 101 or more WGS traces). Histograms display the length distributions of each of these small RNA populations. Reads that corresponded to nuclear rRNA were not included.

A similar analysis of the S. moellendorffii small RNA reads found that 35% of them corresponded to miRNA hairpins, most of which were 21 nucleotides in length (Figure 5B). The remaining 65% of the reads that matched at least one S. moellendorffii WGS trace were binned into repetitive and nonrepetitive subpopulations, which comprised 11 and 54% of the total, respectively (Figure 5B). The length distributions of the S. moellendorffii unknown reads were dispersed and overall were not consistent with the small RNA sizes produced by the biochemically or genetically characterized DCL proteins from other species. However, based upon the experience with P. patens small RNAs, we expect that approximately one-half of these unknown small RNAs will not be matched to the forthcoming S. moellendorffii draft genome sequence. As with P. patens, subtracting these reads may reveal the signature of a small population of potential miRNA-independent siRNAs in S. moellendorffii.

The proportion of P. patens and S. moellendorffii potential endogenous siRNAs was lower than that observed in similar deep sequencing projects that examined Arabidopsis small RNAs (Lu et al., 2005; Rajagopalan et al., 2006; Kasschau et al., 2007), raising the possibility that endogenous siRNAs play a less prominent role in genome regulation in basal plants. However, it is also possible that basal plant endogenous siRNAs possess biochemical characteristics that render them recalcitrant to the sequencing strategy employed in this study, which is highly selective for molecules with both a 5' monophosphate and a 3' hydroxyl. This would not be without precedent. In Caenorhabditis elegans, secondary siRNAs from an exogenous trigger and a subset of endogenous siRNAs are recalcitrant to this procedure because they possess 5'-triphosphates, reflecting their likely origins as short unprimed RdRp transcripts (Ruby et al., 2006; Pak and Fire, 2007; Sijen et al., 2007). Angiosperm siRNA and miRNA/miRNA* duplexes are methylated by the HEN1 protein on the 2' hydroxyl of the 3'-most residue (Yu et al., 2005; Yang et al., 2006). Periodate treatment of P. patens total RNA demonstrated that the 3'-most residue of miR160 does not possess both free 2' and 3' hydroxyls, which was strongly suggestive of a conservation of HEN1-mediated small RNA methylation in bryophytes (Figure 6 ). Although methylation of the terminal 2' hydroxyl would not have had a major effect on sequencing frequency, other modifications of this terminal ribose might have, especially if different classes of basal plant small RNAs possess different modifications.

View larger version (31K): [in this window] [in a new window] [as a PowerPoint slide]   Figure 6. The 3'-Most Residue of P. patens miR160 Is Modified. The indicated RNA samples were oxidized with periodate, blotted, and probed for miR160. Markers (M) from top to bottom are 24, 21, and 18 nucleotides, respectively.

	DISCUSSION

TOP ABSTRACT INTRODUCTION RESULTS DISCUSSION METHODS REFERENCES

Common Functions for Plant Small RNAs, with and without Conservation of Sequence Identity
The targets of the deeply conserved P. patens miRNAs were clearly homologous with the known targets of the same miRNAs in Arabidopsis, as observed previously for deeply conserved miRNAs (Floyd and Bowman, 2004; Arazi et al., 2005; Axtell and Bartel, 2005; Floyd et al., 2006). The sequences of these targets strongly suggest that most of these deeply conserved regulatory circuits function to control transcriptional regulators that influence multicellular development and morphology. It will be fascinating to experimentally dissect how these conserved molecular modules have been recruited to control the very different developmental programs of angiosperms and bryophytes.

Analysis of the P. patens miRNA loci and their targets strongly suggests that miRNA-mediated feedback control on the miRNA pathway is also shared between bryophytes and angiosperms through the action of diverse miRNA loci generated through convergent and perhaps divergent evolutionary pathways. P. patens MIR1047 is expressed from the intron of a DCL1 homolog, and just as for Arabidopsis MIR838, processing of the stem loop by DCL1 cleaves the DCL1 pre-mRNA, implying a negative feedback loop. These two miRNAs appear to have arisen independently yet serve a common function and thus appear to constitute an example of convergent evolution for miRNA functions. Considering that large reservoirs of newly emergent and continuously evolving miRNAs have now been described in an angiosperm, a lycopod, and a moss, we expect that many additional examples of convergent miRNA functions will be found in land plants. Whether miR904 targeting of AGO1-like genes in P. patens and miR168 targeting of AGO1 in angiosperms was achieved through convergent or divergent emergence of these diverse miRNAs was less clear. Biogenesis from the same arm of the miRNA hairpins and cleavage at similar locations within the homologous AGO1 messages were both consistent with divergence from a common ancestor, but in either evolutionary scenario, our data demonstrated that analogous regulatory specificity of plant miRNAs does not always entail conservation of the mature miRNA sequences. Similar observations have been made for animal miRNAs, most notably the vertebrate miR-196 and insect miR-iab-4, which both target Ubx orthologs in their respective lineages (Yekta et al., 2004; Ronshaugen et al., 2005).

Just as for the miRNAs that appeared to mediate feedback regulation upon DCL and AGO genes in plants, the biologically relevant ta-siRNAs emanating from plant TAS3 loci differed in sequence but were found to perform analogous molecular functions in mosses and seed plants. The two P. patens ARF homologs possessed dual cleavage sites just downstream of the conserved ARF domain (Figure 2B), a location and arrangement closely mirroring that of the cleavage sites within Arabidopsis ETT/ARF3 and ARF4 (Allen et al., 2005). These data involving miR390-triggered siRNA regulation of similar angiosperm and bryophyte ARF transcripts could thus be regarded as evidence that the At TAS3 and Pp TAS3a-d genes might have descended from a common molecular ancestor. However, Arabidopsis ETT/ARF3 and ARF4 both possessed dual complementary sites to identical At TAS3–derived ta-siRNAs, whereas the P. patens ARFs possessed only one site complementary to Pp TAS3 ta-siRNAs. Cleavage at the other site in the P. patens ARF mRNAs was guided by the lineage-specific miR1219, which did not have any detectable sequence similarity to the ARF-targeting At TAS3 ta-siRNAs nor to any Pp TAS3a-d–derived ta-siRNAs. These observations add to the examples of convergent evolutionary origins for similarly functioning small RNAs, while illustrating that cleavage of a homologous target at an analogous location does not always indicate descent from a common ancestor. Moreover, the ARF-targeting ta-siRNAs are derived from the (+) strand of the angiosperm TAS3 loci, whereas they derived from the (–) strands of the moss TAS3 loci, suggesting a more complex evolutionary path for plant TAS3 loci than simple sequence divergence from a common ancestor. Data from additional plant lineages might provide important clues for distinguishing between these interesting possibilities.

The AGO1-targeting miRNAs and the miR390-phased ta-siRNAs might be the descendants of anciently derived small RNA–target interactions that have coordinately diverged from common molecular ancestors over long periods of time. This process would thus have hidden the ancestral relationships of the functionally homologous small RNA–target interactions in divergent contemporary organisms, making them appear to be nonconserved when considering primary sequence alone. In other cases, such as miRNAs excised from DCL1 pre-mRNAs, independently derived, nonconserved miRNAs appear to have converged upon an identical regulatory function. We conclude that the sequence diversity of miRNAs and ta-siRNAs expressed in various plants can sometimes obfuscate analogous regulatory functions.

	METHODS

TOP ABSTRACT INTRODUCTION RESULTS DISCUSSION METHODS REFERENCES

 
Small RNA sequencing
Sequencing of Physcomitrella patens small RNAs has been described (Axtell et al., 2006). Selaginella moellendorffii was purchased from Plant Delights Nursery. A total RNA sample derived from the above-ground parts of S. moellendorffii was isolated using the method of Chang et al. (1993). Small RNAs from this sample were sequenced as described (Axtell et al., 2006).

Annotation of miRNAs and Target Predictions
Sequenced small RNAs were filtered to remove contaminant sequences derived from rRNA or the synthetic small RNAs used as guides during the small RNA excision procedure. The remaining sequences were then compared against the relevant WGS traces available from the National Center for Biotechnology Information (NCBI) trace archive and the corresponding matches recorded. After discarding small RNAs that were repetitive (matched >1000 different WGS traces), genomic sequences surrounding small RNA hits were examined for their potential to form an RNA secondary structure typical of pri-miRNA stem loops. A 600-nucleotide window centered upon a given small RNA hit was analyzed using MIRcheck (Jones-Rhoades and Bartel, 2004). Consensus genomic sequences were derived for each potential miRNA locus that passed MIRcheck by use of BLAT (Kent, 2002) to capture redundant WGS traces followed by consensus assembly using CAP3 (Huang and Madan, 1999). All sequenced small RNAs that matched either polarity of the assembled consensus sequences of potential miRNA loci were captured and aligned. Candidate loci where two dominant sequences accumulated from only a single polarity and that, when analyzed in the context of a predicted RNA secondary structure, could form a typical miRNA/miRNA* duplex were annotated as miRNAs. In cases where no sequence corresponding to the presumed miRNA* was sequenced, a locus was still annotated as a miRNA if there was a paralogous locus from which a miRNA/miRNA* was isolated. After the P. patens draft genome Phypa1_1 assembly became available, the WGS-derived consensus miRNA loci were all successfully mapped.

The transcripts annotated from the draft P. patens genome assembly were used to predict the potential targets of P. patens miRNAs. The list of mature P. patens miRNAs (see Supplemental Table 1 online) was first collapsed to generate a nonredundant set of queries, and targets were predicted as described (Allen et al., 2005; Rajagopalan et al., 2006). To assess the signal-to-noise ratio of these predictions, a set of randomized sequences with mononucleotide compositions identical to the query and di- and trinucleotide compositions similar to the transcript database (Farh et al., 2005) was produced for each miRNA query and used to predict targets (see Supplemental Table 3 online).

Molecular Phylogenetic Analyses
P. patens proteins homologous to Arabidopsis thaliana DCL1 (At1g01040.1), AGO1 (At1g48410.1), and ETT/ARF3 (At2g33860.1) were found using<