High-Throughput Sequencing of Arabidopsis microRNAs: Evidence for Frequent Birth and Death of MIRNA Genes

In plants, microRNAs (miRNAs) comprise one of two classes of small RNAs that function primarily as negative regulators at the posttranscriptional level. Several MIRNA genes in the plant kingdom are ancient, with conservation extending between angiosperms and the mosses, whereas many others are more recently evolved. Here, we use deep sequencing and computational methods to identify, profile and analyze non-conserved MIRNA genes in Arabidopsis thaliana. 48 non-conserved MIRNA families, nearly all of which were represented by single genes, were identified. Sequence similarity analyses of miRNA precursor foldback arms revealed evidence for recent evolutionary origin of 16 MIRNA loci through inverted duplication events from protein-coding gene sequences. Interestingly, these recently evolved MIRNA genes have taken distinct paths. Whereas some non-conserved miRNAs interact with and regulate target transcripts from gene families that donated parental sequences, others have drifted to the point of non-interaction with parental gene family transcripts. Some young MIRNA loci clearly originated from one gene family but form miRNAs that target transcripts in another family. We suggest that MIRNA genes are undergoing relatively frequent birth and death, with only a subset being stabilized by integration into regulatory networks.

The expanse of genetic information regulated posttranscriptionally by small RNAs is potentially large in animals and plants [11][12][13]. In humans, for example, computational and indirect experimental evidence indicates that miRNAs regulate expression of up to 1/3 of all genes [14][15][16][17]. In plants, far fewer mRNAs are directly regulated by miRNAs, although the direct and indirect consequences of miRNA-directed regulation are significant [11,12]. This is due to the roles of a large proportion of plant miRNA target transcripts that encode transcription factors required for normal growth, development, hormone response, meristem functions and stress responses [11,12,18,19]. Approximately 21 families of Arabidopsis miRNAs and their respective targets are conserved in Rice and/or Poplar [11]. Additionally, there is a growing recognition of significant numbers of miRNAs not conserved in Rice or Poplar, many of which likely arose in the recent evolutionary past [20][21][22][23][24]. In some cases, these MIRNA loci formed through inverted duplication events that yielded transcripts with self-complementary foldback structure [22]. In fact, the genomes of Arabidopsis and other plants contain a wide diversity of non-conserved, local inverted duplications that yield small RNA populations ranging from highly uniform, DCL1-dependent miRNAs (e.g. MIR163) to heterogeneous collections of bidirectional short interfering RNAs (siRNAs) formed by multiple DCLs [22,23,25,26].
Deep sequencing methods now provide a rapid way to identify and profile small RNA populations in different plants, mutants, tissues, and at different stages of development. We and others have used high-throughput pyrosequencing to analyze small RNAs across the Arabidopsis genome in wild-type and silencing-defective mutants [23,25,[27][28][29]. In this paper, we identify and analyze non-conserved and recently evolved Arabidopsis miRNAs. The data reveal a relatively large number of miRNAs that are so far unique to Arabidopsis, and suggest that many miRNAs are spawned and lost frequently during evolution.
A computational analysis to identify new MIRNA genes was done using a protocol similar to that of Xie et al. [31] ( Figure 1A). All small RNAs from the ASRP database (Set1) were used ( Figure 1A). Briefly, all loci from Set1 that yielded at least two reads (Set2) were subjected to Repeatmasker [32] and bidirectional small RNA cluster filters to eliminate siRNAs from repeat sequence classes. Small RNAs that differed at their termini by up to four nucleotide positions were consolidated and passed through a self-complementary foldback screen with settings as described [31]. Small RNAs from coding sequences and complex small RNA clusters, or that were not 20-22 nt in length, were eliminated, yielding 228 loci (Set3). These included 91 of 102 (false negative rate = 0.11) rigorously characterized MIRNA genes used as a ruledevelopment and reference set (Table 1) [31]. Failure to identify miR398a, miR399a, miR399d, miR399e, miR399f and miR447c was due to a lack of sequence reads. Known miRNAs also failed due to length (miR163), an incorrectly predicted foldback (miR164b and miR167d) or because the MIRNA gene yielded bidirectional small RNAs (miR156d and miR161). Small RNAs from known MIRNA genes were then removed from Set3, leaving 79 unique small RNA loci ( Figure 1A). Two loci that had each failed at one filter step were reclaimed manually due to their abundance and their dependence on DCL1 (http://asrp.cgrb. oregonstate.edu/db/). The 81 loci (Set4) were then subjected to a detailed foldback analysis where small RNAs from the same foldback were consolidated into a single prospective MIRNA locus. A particularly useful, although not absolute, feature was detection of miRNA* sequences, which arise from the opposite foldback arm during DCL1-mediated processing. Detection of miRNA* sequences reveals the functionality of the predicted foldback. miRNA* sequences were detected at relatively low abundance for most reference MIRNA families (Table 1).
Thirty-nine loci (Set5) beyond the reference set emerged as prospective miRNAs, including 25 for which miRNA* sequences were detected (Table 2). Eight miRNAs from Set5 were detected recently as miRNAs by Lu et al. [23]. Only one sequence, corresponding to Populus trichocarpa miR472 [24], was conserved in Poplar based on sequence/foldback similarity and target conservation. None of the new sequences were conserved in Rice. In several cases (miR830, miR833, miR851, miR861, miR862, miR863, miR864, miR865 and miR866), due to the ratio of small RNA sequences from both the 5p and 3p strands of the foldback, it was not clear which small RNA should be designated as the miRNA. In such cases, both the 5p and 3p strand sequences are designated (Table 2). In total, at least 70 MIRNA families were detected in at least one sequenced population (Tables 1 and 2, and references therein). Most miRNAs were lost or underrepresented in the dcl1-7 mutant (http://asrp.cgrb.oregonstate.edu/db/). Target RNAs were predicted from both transcript and EST databases for the 39 miRNAs in Set5 by the method of Allen et al. [22], which is based on additive, position-dependent mispair penalties. Seventy-eight or 69 of 85 validated, reference miRNA targets were predicted with a threshold score of 4 or 3.5, respectively (false negative rate = 0.08 or 0.19; Figure 1A). A total of 142 targets were predicted for the Set5 miRNAs, using the more conservative 3.5 score threshold ( Table 2). These included targets predicted for both 5p and 3p strands for eight MIRNA loci. The 3.5 score threshold was used to maintain specificity, although at the cost of sensitivity, in the target prediction algorithm.
Top-scoring predicted targets for most miRNAs from Set5, and two that were not found in this study (miR778 and miR781) [23], were tested using a standard 59RACE analysis to detect cleavage events at predicted sites opposite nucleotide 10 from the 59 end of the small RNA. The 59RACE assays were done using two gene-specific primer sets with RNA from whole seedling and inflorescence tissues. Thirteen targets for 13 miRNAs were validated, although evidence supporting miR775-and miR859-guided cleavage (At1g53290 and At3g49510, respectively) was weak ( Figure 1B, Table 2). Evidence for miR846-and miR844-guided cleavage should also be interpreted cautiously, as targets were not cleaved at the canonical position. These were considered validated because both miRNAs had SET-domain (1) e , At3g49510 (1.5), [32] MIRNAs with only predicted targets, or no predicted targets miR771   [4] miR866-5p c Top three predicted targets with a score of 3.5 or less are listed with their score in parentheses. Targets validated by 59 RACE are in bold. Remaining number of targets predicted with a score of 3.5 or less are listed in square brackets (Table S1). Dashes indicate no predicted targets with a score of 3.5 or less.  Figure 1B; http://asrp.cgrb. oregonstate.edu/db/). Additionally, targets for miR472 and miR774 were validated recently [23]. Each target-validated miRNA was represented by at least ten reads or had a predicted miRNA* sequence represented with at least one read in the database, except for miR774 and miR778 which were not sequenced here (Table 2). Twenty-two predicted miRNAs that were sequenced at least two times, or that had miRNA* sequences represented in the ASRP or MPSS databases [33], failed at the target validation step ( Table 2). The relatively high proportion of target validation failures could be due to erroneous target predictions, low-abundance targets or miRNA-guided cleavage products, or low-abundance miRNAs with limited or no activity. Additionally, for 13 miRNAs, no targets were predicted at a score threshold of 3.5 ( Table 2). While it is possible that many of these function to guide cleavage of unpredicted target RNAs, it is also possible that miRNAs exist without actual targets (see below).
Several new target families are worth noting. First, at least two targets are clearly under negative regulation by miRNAs in inflorescence tissue. Transcripts from AGL16 (At3g57230, MADSbox) and MYB12 (At2g47460, R2R3-MYB family) genes are targeted by miR824 and miR858, respectively, and are each upregulated by 2-4 fold in dcl1-7 and hen1-1 mutants (Table 2, Figure 2A). Functions for these specific transcription factors are not known. Second, the transcript for the cation/hydrogen antiporter gene CHX18 is targeted by both miR780 and miR856 (Table 2, Figure 1B). Dual targeting of the CHX18 (At5g41610) transcript may lead to secondary, 21 nt siRNAs that arise through the RDR6/DCL4-dependent pathway ( [29,[34][35][36][37][38][39]; data not shown). miR856 may also target an unrelated transcript encoding the efflux protein ZINC TRANSPORTER OF ARABIDOPSIS THALIANA1 (ZAT1). Third, miR857 was validated to target the LAC7 (At3g09220) transcript, which encodes a laccase family protein with predicted multicopper oxidase function. This represents the third miRNA family to target the laccase gene  (Tables 1 and 2). (C) Relative numbers of miRNA target family functions for conserved and nonconserved miRNAs (Tables 1 and 2). Only target classes that have been validated experimentally are included. Note that Table 2 shows many MIRNA families with weak or no predicted targets, and these are not represented in the chart. doi:10.1371/journal.pone.0000219.g002 family [18,40,41]. In fact, at least seven gene or domain families (encoding MYB, PENTATRICOPEPTIDE REPEAT (PPR), AUXIN RESPONSE FACTOR (ARF), AGO, F-box domain, kinase and LAC proteins) are now known to be targeted by multiple miRNA families [11].

Direct sequencing as a miRNA expression profiling tool
The primary purpose of this study was to analyze new MIRNA loci. Nevertheless, the data from multiple mutants and treatments allowed us to explore the suitability of direct sequencing as a miRNA profiling tool. This was done in several ways. First, the reproducibility of the pyrosequencing method was analyzed using two normalized biological replicates [Col-0 inf. 1 (35,666 reads) and Col-0 inf. 2 (42,917 reads)]. miRNA families represented in both samples by at least three reads were compared and yielded an 83% correlation ( Figure 3A). This experimental approach is essentially a random sampling from a very large population of small RNAs and hence we anticipate that this correlation would increase with greater sequencing depth.
We next compared normalized expression profiles of miRNA families in Col-0 (inflorescence tissues) and dcl1-7 (inflorescence tissues) mutant plants. As expected, based on the requirement for DCL1 in miRNA precursor processing, levels of nearly all miRNA families were decreased in dcl1-7 mutant plants ( Figure 3B). A few miRNAs were largely unaffected in dcl1-7 plants. This likely reflects the known residual activity in, and the insensitivity of some miRNAs to, the partial-loss-of-function dcl1-7 mutant allele [22].
To extend our small RNA profiling beyond different tissues and developmental stages, we generated expression profiles of miRNA families in leaves following infection by P. syringae pv. tomato (DC3000hrcC) at 0, 1 and 3 hr post-inoculation (p.i.) timepoints. This bacterium is unable to cause disease, because it is mutated in the key virulence factor delivery machine, the type III secretion system. This strain does trigger a robust basal defense response in Arabidopsis [42]. As predicted [43], miR393 was strongly upregulated (10-fold at 3 hr p.i.) ( Figure 3C). miR393 targets mRNAs encoding the auxin receptor, TRANSPORT INHIBI-TOR RESPONSE1 (TIR1), and related proteins [18,43]. Two additional miRNA families, miR160 and miR167, were significantly elevated at 3 hr p.i. by 5-fold and 6-fold, respectively ( Figure 3C). miR160 and miR167 each target mRNAs encoding members of the ARF family of transcription factors [19]. These data suggest that basal defense responses triggered by P. syringae pv. tomato (DC3000hrcC) include miRNA-mediated suppression of multiple components of auxin signaling pathways. Hence, our data extend the conclusions of Navarro et al. [43]. miR825, for which several targets were predicted but not validated (Table 2), was significantly down-regulated 3 hr p.i. (Figure 3C).
Thus, direct sequence-based profiling shows considerable promise in revealing dynamic changes in miRNA populations, although shallow sequencing depth will limit applicability of the method. While both open-ended platforms (such as direct sequencing) and closed-ended platforms (such as solid-state microarrays) can be used to profile known miRNAs, the sequence-based profiling approach affords discovery of previously unknown miRNAs. It also allows discrete measurements of complex mixtures of small RNAs that arise in heterogeneous populations from loci with diffuse boundaries [23,27,28].

Conserved vs. non-conserved miRNAs
Twenty-two and 20 MIRNA families (Tables 1 and 2) are conserved in Poplar and Rice, respectively. Conserved families in Arabidopsis were designated by having identical or related (three or fewer nucleotide substitutions) sequences, and at least one conserved target transcript, in either Poplar or Rice. A series of general and functional comparisons between conserved and nonconserved families was done. Whereas 19 of 22 conserved Arabidopsis miRNA families were represented at multiple loci, only three (MIR158, MIR447 and MIR845) of 48 non-conserved miRNAs were members of multigene families ( Figure 2B). Expansion of multigene MIRNA families has occurred through tandem and segmental duplications, as well as polyploidization events [22,44,45]. The preponderance of single gene families among the non-conserved miRNAs is consistent with recent evolutionary derivation.
Functions of conserved miRNA target genes, as a group, are less diverse than functions for non-conserved miRNA target genes (only validated or high-confidence targets were compared). This is due primarily to the high proportion of target mRNAs encoding transcription factors for conserved families (Table 1 and 2, Figure 2C). As pointed out clearly before [18,19], the vast majority of transcription factor families targeted by conserved miRNAs participate in developmental pathways, including those specifying meristem functions, organ polarity, cell division control, organ separation and hormone responses. Non-conserved miRNA target genes encode a broad range of proteins, including a limited number of transcription factors (Table 2, Figure 2C). Interestingly, several non-conserved miRNAs (miR161 and miR400) target transcripts from a clade within the large PPR family [19,41]. Another non-conserved miRNA, miR173, targets tasiRNA primary transcripts (TAS1 and TAS2) [34], which in turn yield siRNAs that also target several of the miR161-and miR400targeted PPR transcripts ( [29]; data not shown).
To what extent do the non-conserved miRNAs negatively regulate target genes? The sensitivity of target genes of conserved and non-conserved miRNAs was analyzed by transcript profiling in wild-type (Col-0 and La-er) plants and mutant plants with general miRNA deficiencies (dcl1-7 and hen1-1). To avoid biasing the analysis with miRNAs that target disproportionately high numbers of target mRNAs, such as miR161 and PPR target gene family members, the number of genes analyzed for each miRNA family was limited to two validated targets, or two predicted targets with the lowest scores. For miRNAs that have been shown to target multiple gene families, such as miR395, two targets from both gene families were analyzed. A scatterplot of fold-change in dcl1-7 and hen1-1 mutants for each target was generated. As shown previously [34], conserved miRNA target transcripts displayed a generally elevated pattern in both mutants (Figure 2A). In contrast, most target transcripts of non-conserved miRNAs were clustered around the origin, indicating that most were insensitive to either the dcl1-7 or hen1-1 mutation (Figure 2A). The exceptions that were affected at levels of 1.6-fold or greater in either mutant included transcripts from AGL16, MYB12 and a PPR gene (At1g63130), which were targeted by miR824, miR858 and miR161.1, respectively (Figure 2A). While some of these data may be skewed by tissue-specific expression patterns of miRNAs and target genes, the general trend for non-conserved miRNAs having fewer effects on target transcript levels is clear. We suggest that, in contrast to the vast majority of conserved miRNAs, a high proportion of the non-conserved miRNAs are not integrated as dominant factors within regulatory networks.

Recent evolution of non-conserved MIRNA loci
Previously, we identified a number of small RNA-generating loci with the potential to yield transcripts with self-complementary foldback potential and with extensive similarity to protein-coding gene family sequences [22]. Whereas the majority of these loci yield complex populations of siRNAs, two (MIR161 and MIR163) yield functional, non-conserved miRNAs [22]. This led to the idea that aberrant replication/recombination or transposition events from expressed gene sequences can spawn new small RNAgenerating loci with the potential to evolve into MIRNA genes that regulate members of the originating family.
This idea was tested more rigorously with the expanded set of MIRNA loci. Foldback sequences for each MIRNA locus (Tables 1   and 2) were used in FASTA searches against Arabidopsis transcript and gene databases ( Figure 4A). Nearly all foldback sequences from conserved miRNAs had hits with non-significant E-values greater than 0.05 ( Figure 4B). In contrast, 19 of 48 non-conserved miRNA foldback sequences had at least one hit with an E-value lower than 0.05 ( Figure 4B). Similarity or complementarity was detected on both 59 and 39 arms containing miRNA or miRNA* sequences.
It is important to recognize that MIRNA loci contain two regions -the miRNA and miRNA-complementary region (largely over-  lapping with the miRNA*) -with relatively high levels of complementarity or similarity to target genes. To eliminate the potential misleading influence of these sequences, which may be under selection due to the requirement for complementarity between miRNAs and their targets, on the similarity test, each foldback arm with hits (E,0.05) in the FASTA search was analyzed independently with and without the miRNA or miRNAcomplementary sequences. The top four FASTA hits ( Figure 4A) were aligned with the intact or deleted arms using a global sequence alignment method (NEEDLE). MIR163 and MIR161 were analyzed in detail previously [22] and were included here as controls. Two foldback arm sets from MIR161, with miR161.1/ miR161.1-complementary and the overlapping miR161.2/ miR161.2-complementary sequences deleted independently, were analyzed. Sixteen foldbacks, including the MIR161 and MIR163 controls, contained at least one arm with similarity or complementarity (NEEDLE score with p,0.001) to one or more genes when an alignment was done with both intact and deleted arms ( Figure 4C, Table S3). For all hits with the deleted arms, the topscore gene alignments identified sequences that were immediately flanking the region with similarity or complementarity to miRNA regions ( Figure 5). Further, in all cases in which multiple genes were hit with intact and deleted arms, the genes were closely related members of a single family. Thus, over 30% of the nonconserved MIRNA loci show evidence of common origin with specific genes or gene families.
Interestingly, the two MIR824 arms aligned best with distinct regions within one gene (AGL16). The miR824 arm was complementary to a region from exon VIII, and the miR824* arm was most similar to a duplicated region located within intron III. The two arms from MIR846 were most similar to two adjacent but distinct regions of jacalin domain-containing genes or pseudogenes (At1g57570 and At1g61230). It is likely, therefore, that some recently evolved MIRNA loci arose through juxtapositioning of sequences from two related duplicated sequences.
The protein-coding genes with extended MIRNA foldback arm similarity were analyzed in more detail. The intact foldback arms from 13 of the 16 gene-similar MIRNA loci were aligned with up to three gene sequences, and alignment quality was displayed using heat maps ( Figure 5; MIR161, MIR163 and MIR447c were not included). This revealed several sequence conservation patterns. For several MIRNA arm-gene alignments, including those containing miR778, miR780, miR824 and miR856 sequences, the aligned region containing miRNA sequences was clearly more conserved than the remaining arm segments ( Figure 5). This suggests that selection may be operating on several of the recently evolved miRNA sequences, and is indirect evidence for functional significance of the miRNA. Indeed, miR824 targets the AGL16 transcript, which is under miRNA pathway-dependent repression (Figure 2A).
For several other MIRNA arm-gene alignments, such as those containing miR862, miR447a and miR853, the miRNA-region of the alignment is similar to, or weaker than, the alignment containing flanking arm sequences ( Figure 5). In fact, among the 16 MIRNA-similar gene sets identified in this analysis (Table S3), only seven represented sets that corresponded to the best predicted or validated miRNA-target pair (Figures 4 and 6). In other words, nine of the miRNAs from loci with similarity to protein-coding genes were predicted to target transcripts from different genes (based on best target prediction scores). The target scores for eight of these sets were clearly weak or functionally implausible ( Figure 6, Table S3). Interestingly, miR447a and miR856 were validated to target transcripts in gene families distinct from those similar to the foldback sequences ( Figure 1B, Table 2; [34]). The miR447 family appears to have acquired novel target specificity and lost sequence-of-origin specificity, while miR856 may have acquired dualtargeting specificity ( Figure 6, Table 2, Table S3).
The picture that emerges from analysis of the non-conserved Arabidopsis miRNAs appears to show that new MIRNA loci are forming frequently through duplication events. These newly evolved MIRNA loci may first pass through a stage in which heterogeneous populations of siRNA-like sequences are generated, especially if the duplicated locus and resulting foldback sequence is large [22]. This, of course, assumes that the newly spawned sequence is proximal to a functional promoter. Computational evidence indicates that formation of MIR163 through inverted duplication of SAMT methyltransferase-like sequences also involved duplication of the gene-of-origin promoter [46]. Given that DCL1 has limited or insignificant activity on perfectly-paired dsRNA, acquisition of DCL1-dependence and subsequent formation of discrete small RNA products from the foldback precursor may require accumulation of drift mutations that result in foldback mispairs. From this point, we envision three evolutionary fates. The first and perhaps most common is continued sequence drift through mutation, decay of both targeting capacity and promoter (if not required for other functions) sequences, and eventual death of the locus. New MIRNA genes with neutral effects are predicted to take this route. The second fate is stabilization of the miRNAgenerating sequence with specificity for gene-of-origin, or familyof-origin, sequences. This would occur if an advantage is realized when the target gene or genes are brought under negative regulation by the miRNA, and would lead to selection of the miRNA sequence. The miR824-AGL16 regulatory pair may exemplify this evolutionary track. The third fate is chance acquisition of targeting specificity for a novel target gene or family, followed by stabilization through selection in the event of an advantage. miR856, which may target transcripts from both gene-of-origin (ZAT1) and a novel gene (CHX18), could conceivably fall into this category, although it should be noted that there are no direct data supporting a functional role for miR856. We postulate that many of the conserved MIRNA genes arose through the latter two routes, and have lost sequence relatedness to their genetic origin loci due to drift of foldback arm sequences outside of the miRNA/miRNA-complementary regions ( Figure 4).
While this paper was under review, an article from Rajagopalan et al. [47] describing results from deep sequencing of Arabidopsis small RNAs (887,000 reads) was published. They identified and named 32 new miRNAs (miR822-miR853), and another 39 candidate miRNAs that were not named. Sixteen of the named miRNAs, and eight of the candidates, were identified as miRNAs in this study (Table 2). In this paper, 16 miRNAs (miR845b, miR856-miR870) were named, eight of which were identified only here ( Table 2). Differences in tissue sampling and sequencing depth likely account for most of the differences in miRNAs identified between the two studies.

miRNA identification and target prediction
A set of computational filters based on those developed by Xie et al. [31] were used to identify new miRNAs from among sequences in the ASRP database (http://asrp.cgrb.oregonstate.edu/db/), but with the following modifications. First, only small RNAs that were represented by two or more reads were considered. Second, small RNAs arising from repeat elements [32] and bidirectional siRNA clusters [28] were removed. Third, computational assessment of foldback structure was done with sequences containing 250 nts on each side of candidate miRNAs using RNAfold, Vienna RNA package, version 1.6.1 [52].
miRNA targets were computationally predicted as described [34]. Briefly, potential targets from FASTA searches (+15/210 match/mismatch scoring ratio, -16 gap penalty and a RNA scoring matrix) were scored using a position-dependent, mispair penalty system [34]. Penalties were assessed for mismatches, bulges, and gaps (+1 per position) and G:U pairs (+0.5 per position). Penalties were doubled if the mismatch, bulge, gap, or G:U pair occurred at positions 2 to 13 relative to the 59 end of the miRNA. Only one single-nt bulge or single-nt gap was allowed. Based on a reference set of validated miRNA targets, only predicted targets with scores of four or less were considered reasonable. Conservation between Arabidopsis and Poplar (P. trichocarpa) was assessed by FASTA search, foldback analysis, and detection of similar target sequences [18,48].

miRNA target validation assays
Target validation using a 59 RACE assay was done with the GeneRacer Kit (Invitrogen, CA) as described previously [22,34,53,54]. Poly(A)+ mRNA was isolated from seedling (7 day) and inflorescence tissue (28 day, stage 1-12 flowers) of Col-0 plants, ligated to adaptor, converted to cDNA and subjected to two rounds of PCR amplification using gene-specific and adaptorspecific primers [22,34,53,54]. Amplified products were gelpurified, cloned and sequenced. Gene-specific primer sequences for miRNA targets that were successfully validated are shown in Table S4.

Foldback sequence similarity analysis
The sequences comprising foldbacks from all MIRNA loci (Tables 1  and 2) were identified using RNAfold. Foldback sequences were subjected to FASTA searches against the Arabidopsis gene and Figure 6. Targeting specificity of recently evolved MIRNAs. Two target prediction scores are shown for each of 16 miRNAs: best overall predicted target score (blue) and target scores calculated for MIRNA foldback-similar genes (grey). Left column indicates whether or not the best overall predicted target gene is in the same family as the foldback-similar gene. A dot indicates that the predicted gene is in an experimentally validated target family. Two calculations corresponding to the two major populations from the MIR161 locus (miR161.1 and miR161.2) are shown. The identities of targets are listed in Supplemental Table S3. The plot is centered on a target prediction score of 4, as this corresponds to the upper limit of a reasonable prediction. doi:10.1371/journal.pone.0000219.g006 transcript databases [22]. The 59 and 39 arms of foldbacks that had gene hits with E-values lower than 0.05 were individually aligned to the top four FASTA hits using NEEDLE [55]. Each arm was also randomly shuffled 1,000 times using SHUFFLESEQ [55] and realigned to each of the top four FASTA hits. The mean 6standard deviation of the randomized sequence scores was calculated. A Z-score was calculated for each arm-gene pair by subtracting the average score of the randomized sequence alignments from the score of the arm-gene alignment and then dividing by the standard deviation of the randomized alignments. This was repeated for arms in which the miRNA or miRNAcomplementary sequences were deleted from their respective arms. The deleted arms were aligned with gene sequences in which the target sequence was correspondingly deleted. Intact foldback arms and most-related gene sequences were also aligned and viewed using heat maps with T-COFFEE [56].