High-throughput sequencing of Arabidopsis microRNAs: evidence for frequent birth and death of MIRNA genes

Abstract

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.

Figure 1

Identification and analysis of Arabidopsis miRNAs and targets. (A) Flowchart for the prediction of miRNAs and their targets. (B) Validation of predicted targets for 13 non-conserved miRNAs. Positions of dominant 5′ RACE products (no. 5′ ends at position/total no. 5′ ends sequenced) are indicated by vertical arrows in the expanded regions. Predicted cleavage sites are indicated by a bolded nucleotide at position ten relative to the 5′ end of the miRNA or miRNA-variant. Positions of gene-specific primers are indicated with horizontal arrows above gene diagrams.

Figure 2

Comparison of conserved and non-conserved MIRNA families. (A) Effect of dcl1-7 and hen1-1 mutations on levels of target transcripts for conserved (black) and non-conserved (red) miRNAs. Expression data are shown for two validated or high-confidence predicted targets, if available, for each family. Arrows indicate targets for miR824 (AGL16), miR858 (MYB12) and miR161.1 (At1g63130, a PPR gene). (B) Numbers of gene family members for conserved and non-conserved MIRNAs (Tables 1 and 2). (C) Relative numbers of miRNA target family functions for conserved and non-conserved 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.

Figure 3

Expression profiling of MIRNA families using high-throughput pyrosequencing. (A) Comparison of most-abundant miRNA families between biological replicates of Col-0 inflorescence (inf.) tissue. Normalized reads for each miRNA family member were consolidated. Note that MIR159 and MIR319 derived members were counted separately, even though they are frequently assigned to the same family , . (B) Fold-change of miRNAs in dcl1-7 inflorescence versus Col-0 inflorescence (left axis, bars). Total number of reads for each family is indicated (right axis, green line). As the dcl1-7 mutant contained no reads for many miRNA families, fold-change was calculated using normalized reads+1. This had the effect of dampening fold-change values for low-abundance families. (C) Fold-change of miRNA family reads in leaves at 1 hr and 3 hr post-inoculation with P. syringae (DC3000hrcC) (left axis, bars). Fold-change relative to uninoculated leaves was calculated based on normalized reads as described in panel (B). Total number of reads in the control and inoculated samples is shown (right axis, green line). Grey dashed lines indicate the p = 0.05 upper and lower thresholds.

Figure 4

Identification of MIRNA foldbacks with similarity to protein-coding genes. (A) Flowchart for identification of MIRNA foldbacks with similarity, extending beyond the miRNA target site, to protein-coding genes. (B) Arabidopsis gene or transcript hits in FASTA searches using foldback sequences for all conserved and non-conserved MIRNA loci (Tables 1 and 2). The top four hits based on E-values are shown. (C) Z-scores for the Needleman-Wunche alignment values from MIRNA foldback arms with top four gene or transcript FASTA hits. Alignments were done with intact foldback arms (I), and with foldback arms in which miRNA or miRNA-complementary sequences were deleted (D). Z-scores were derived from standard deviation values for alignments of randomized sequences. In (B) and (C), a red symbol represents an experimentally validated target, a pink symbol indicates a gene from a validated target family, and an open symbol indicates a gene that is distinct from either the validated or predicted target family.

Figure 5

Similarity between MIRNA foldback arms and protein-coding genes. Each alignment contains the coding strand for 1–3 genes, the miRNA* arm, and the miRNA arm. Orientation of the foldback arms is indicated by (+) for authentic polarity and (−) for the reverse complement polarity. Two alignments are given for MIR824 because the two arms are each most similar to distinct, duplicated regions within the AGL16 gene (At3g57230). Alignments were generated using T-COFFEE. Colors indicate alignment quality in a regional context.

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.