Strand selection determinants of human miRNAs

Expression level differences between two miRNAs originating from the same precursor cannot result from differential transcription efficiency; therefore, they must be caused by differential stability due to biased strand incorporation into the RISC. To test whether particular sequence features may explain this bias, we compared the sequence composition of the 33 miRNA pairs that exhibit an expression difference greater than 1000-fold between high- and low-expressed duplex strands to all 847 annotated human miRNAs, as well as the 48 new miRNAs identified in this study (see Additional file 1: Table S1). We find three characteristic sequence features: (i) a U-bias at the 5'-end of the highly expressed strand (One-sided Binomial test, p = 7.9*10^-08); (ii) a C-bias at the 5'-end of the low expressed strand (One-sided Binomial test, p = 4.2*10^-12); (iii) an excess of purines (A and G) in the highly expressed strand and, consequently, an excess of pyrimidines (U and C) in the low expressed strand (One-sided Wilcoxon test, p < 0.02 in all comparisons) (see Table ​Table1,1, Figures ​Figures22 and ​and33).

Open in a separate window

Figure 2

Sequence features of human miRNA strand selection: 5' end nucleotide. Sequence logos showing nucleotide composition of high-expressed strand (A) and low-expressed strand (B) from 33 miRNA pairs with large strand selection bias, and from 103 miRNA pairs with little bias: high-expressed (C) and low-expressed (D).

Open in a separate window

Figure 3

Sequence features of human miRNA strand selection: purine/pyrimidine bias. (A) Purine content of the high-expressed strand from 33 miRNA pairs with large strand selection bias (blue) and the high-expressed strand from all other expressed miRNA pairs (yellow). (B) Pyrimidine content of the low-expressed strand from 103 miRNA pairs with little strand selection bias (red) and the low-expressed strand from all other expressed miRNA pairs (yellow). Purine (C) and pyrimidine (D) content of high- and low-expressed strands from 103 miRNA pairs with little strand selection bias. Colour labels as on panels (A) and (B).

In contrast, we find no specific sequence features characteristic of the miRNA pairs originating from the 103 precursor duplexes with an expression difference between the strands less than 10-fold (One-sided Wilcoxon test, p > 0.5 in all comparisons) (see Figures ​Figures22 and ​and3).3). Further, when the entire range of the expression differences within miRNA pairs is considered, rather than only the two extreme groups, we find a gradual increase in the 5'-nucleotide preference and the purine/pyrimidine content bias as the expression difference within pairs increases (see Additional file 1: Figure S3). This result was not caused by low expression of the miRNA pairs, as restricting the analysis to 50 miRNA pairs with each miRNA represented by more than 10 sequence reads did not affect the result.

To test reproducibility of this finding, we examined miRNA sequences in a technical replicate of the pooled cerebral cortical sample, as well as in a cortical sample from a single adult human male (see Methods). In both cases, we can identify the same sequence features characterising high- and low-expressed strands from precursors with a large strand selection bias, but not strands from precursors showing little bias (see Additional file 1: Figures S4 and S5). Further, all three features can be reproduced in another four samples representing other tissues and species: in the human cell line (HeLa) [26], human embryonic stem cells [27], human embryoid bodies [27], and mouse embryonic stem cells [28] (see Methods and Additional file 1: Table S4). Thus, all sequence features we find associated with miRNA strand selection in human brain are common among mammalian tissues.

Although reproducible, the observed sequence features might potentially be explained by factors other than strand selection bias within miRNA pairs. First, it is conceivable that Illumina sequencing itself is biased towards RNAs with a 5'-terminal U and a high purine content, and discriminates against RNAs with a 5'-terminal C and a high pyrimidine content. After examining all 3,650 unique sequences that can be mapped to the precursor region of known human miRNAs, we find neither a purine/pyrimidine bias, nor a terminal U/C bias (two sided Binomial test, p > 0.5) (see Additional file 1: Figure S6).

Second, the identified sequence features could reflect a common origin for the 33 examined miRNA pairs within the miRNA families. To assess this possibility, we retained only one randomly chosen miRNA pair to represent a family. After doing so, we still observe all three sequence features among the remaining 20 miRNA pairs (p < 0.05 in all tests).

Third, the observed sequence features could be associated with the high expression level of miRNAs, rather than with strand selection bias. Indeed, one of the miRNAs within a pair must be highly expressed to produce a greater than 1000-fold expression level difference. To assess this bias, we analyzed the sequence composition of 10 miRNA pairs with an expression level higher than 925 copies for at least one of the two miRNAs, and an expression difference within the pairs less than 10. Similarly, we sub-sampled 10 of the 33 miRNA pairs with an expression difference greater than 1,000-fold, while retaining a similar copy number distribution (see Additional file 1: Table S3). We find that while all three sequence features are still clearly observable for the 10 miRNA pairs with a large strand selection bias (p < 0.01 in all tests), none of the features are present among the 10 high-expressed miRNA pairs with little strand selection bias (p > 0.4 in all tests). Taken together, these results indicate that the identified sequence features are likely to be associated with strand selection bias within the human miRNA pairs.

Strand selection determinants of fly miRNAs

To test whether the same sequence features determine expression differences within miRNA pairs in another organism, we used published data containing 2,010,618 Roche/454 sequence reads corresponding to 702,945 unique sequences collected from fruit fly (Drosophila melanogaster) S2 cells and from three tissues at 6 developmental stages [29]. After mapping these reads to the fly genome, we find sequences corresponding to 138 known fly miRNAs with a copy number equal or greater than two (see Methods). In addition, we identify 112 novel counterparts of known miRNAs (see Figure ​Figure1B).1B). Thus, in flies as well as in humans, most miRNA precursors produce two expressed mature miRNA sequences.

We then determined whether specific sequence features determine strand selection bias in fly miRNA pairs. Since the fly sequence data has a smaller number of reads, we used expression differences greater than 30 as the strand selection bias cutoff for pairs, resulting in 30 miRNA pairs -- comparable to the number used in the human data analysis. Despite the smaller magnitude of expression differences within these miRNA pairs, we find a significant preference for the 5'-terminal U in the high-expressed duplex strand when compared to all expressed miRNAs (One-sided Binomial test, p = 0.0011). As in humans, the preference for a terminal uridine increases as the expression difference in miRNA pairs increases (from 76% to 85%). In contrast to humans, however, we find a preference for both C and G at the 5'-end of the low-expressed strand (One-sided Binomial test, p = 0.012 and p = 0.0078, respectively). Further, we find no detectable bias in the purine/pyrimidine composition between the high- and the low-expressed strands within 30 miRNA pairs (Wilcoxon test, p > 0.5) (see Table ​Table1).1). Analyzing an additional D. melanogaster dataset containing sequences from 14 small RNA libraries and processed using Illumina platform [30], we obtain almost identical results (see Additional file 1: Table S5). Thus, sequence elements associated with preferential strand selection in humans and flies are not the same.

Strand selection bias and miRNA processing

Can these additional sequence features found in humans facilitate asymmetric strand incorporation into an RISC? The presence of specific nucleotides on the 5'-end of both selected and excluded miRNAs is the first feature. In the plant Arabidopsis thaliana, the 5'-terminal nucleotides determine the sorting of various small RNAs into different Argonaute protein complexes [31]. It is intriguing to speculate that despite the vast differences between RNAi systems in humans and plants, human Argonaute proteins may have retained the ability to discriminate 5'-terminal nucleotides, even though for a different purpose.

Another feature is the preference for purine nucleotides in the high-expressed strand. We speculate that purine residues may facilitate strand loading into the RISC through sequence-independent interactions. Intriguingly, the PAZ domain of Argonaute proteins, which is involved in the initial steps of RNA recognition and binding, contains a large number of invariant aromatic residues involved in RNA binding [32-34]. It is therefore conceivable that the stacking and hydrophobic interactions between these residues and the purine-rich strand of the miRNA precursor duplex may contribute to strand selection.

Finally, as strand selection specificity can result in enormous expression level differences within some miRNA pairs, which in extreme cases can exceed one million copies, it is tempting to speculate that some of the sequence features we find may actively target the excluded miRNA strand for degradation. It is important to mention, however, that sequence features identified here show significant associated with strand selection bias, but do not explain all strand selection preferences found among miRNAs. Notably, recent observation that strand selection preferences may vary among tissues implies that strand selection factors other that sequence features might exist [19].

Strand selection bias and miRNA-based regulation

Next we asked why strand selection bias exists for some miRNA pairs and, conversely, why some miRNA pairs are not biased. In our dataset, 15 of the 20 highest expressed miRNAs have a greater than 1000-fold expression difference between high- and low-expressed miRNA, while 19 showed a greater than 100-fold difference. Similarly, among the 100 most highly expressed miRNAs, 33 and 70 have a greater than 1000-fold and 100-fold expression difference, respectively. Thus, a large proportion of highly expressed miRNA precursors tends to generate only one highly expressed mature miRNA. This is intuitively understandable, since in the absence of strand selection bias each precursor would generate a pair of highly expressed mature miRNAs, collectively affecting the expression of the two groups of target genes. This condition must impose a strong constraint on the expression of an miRNA pair, as two different groups of target genes would have to be regulated in a coordinated fashion among various tissues and developmental stages.

If strand selection bias serves to limit the influence of a second mature miRNA on gene regulation, we expect to find fewer target genes for the low-expressed miRNA strand as compared to high-expressed one. Using the TargetScanS algorithm to identify evolutionarily conserved target sites [35], we indeed find that the high-expressed miRNAs have significantly more potential target genes than their low-expressed counterparts (One-sided Wilcoxon test, p = 0.0000013). This result confirms the notion that strand selection bias limits the regulatory role of a second mature miRNA that can be generated from a highly expressed precursor. It also indicates that the strand selection bias observed in the examined dataset is likely to persist throughout various tissues and various evolutionary lineages of vertebrates.

In contrast to the 33 miRNA pairs exhibiting a large strand selection bias, we find no significant difference in the target gene numbers between the high- and low-expressed miRNAs originating from the 103 pairs showing little bias (One-sided Wilcoxon test, p = 0.09). Further, the numbers of target genes for these miRNAs were significantly lower on average than the corresponding numbers for the high-expressed miRNAs from the 33 biased pairs (One-sided Wilcoxon test, p = 0.00001). Thus, expression of a single miRNA per precursor may be more compatible with the regulation of a large number of target genes due to lesser constraints imposed by the second miRNA. In fact, many of these miRNAs, such as the let-7 family, mir-29 family, mir-26 family, mir-140-3p, mir-30d, mir-7, mir-191, mir-21, mir-100, mir-15a, and mir-16, play an important role in development and differentiation [36-41].

Mapping of sequence reads to the human genome

For each of the three brain sequencing datasets, as well as for hESCs, hEBs and HeLa cells dataset, to remove the adapter sequence at the 3'-end of the sequence reads, all unique sequences were trimmed using the custom trimming procedure. Specifically, we first filtered out the low compositional complexity sequences using the DUST algorithm [45]. The remaining sequences were trimmed by matching the adapter sequence to the 3'-end allowing 3 mismatches if the length of the match was greater than 10 nucleotides and allowing 2 mismatches if the length of the match was between 4 and 10 nucleotides. The trimmed sequences were mapped to the human genome (hg18) [46] using the Illumina-supplied algorithm, ELAND. Only sequences perfectly matching the genome and with a length ranging from 18 to 28 nucleotides were retained. Using these parameters, approximately 60% of total sequence reads of each sample could be mapped. Of these sequences, we expect less than 1% to map incorrectly, as determined by a mapping of scrambled reads with the same length and nucleotide content distribution.

Quantification of known human miRNAs

Human miRNA annotations, including the sequences and genomic locations of miRNA precursors and mature miRNA sequences, were downloaded from miRBase version 11 [47-49]. The genomic location of each mature miRNA was then derived using these data. All mapped sequences with a copy number equal to or greater than two were used for the following quantification procedure. First, all sequences mapping within three nucleotides upstream or downstream of the annotated 5'-position of the mature miRNAs were retained. Then, for each mature miRNA, the sequence with a maximal copy number was designated as the reference sequence. Finally, the expression level of each miRNA was calculated as a sum of the copy number of the reference sequence and the sequences mapping at the same 5'-end position as the reference sequence. The expression measurements obtained using this method and using the annotated mature miRNA sequence as a reference were similar, but the former method resulted in a slightly better correlation between experimental replicates and was used throughout the analysis (data not shown). The only exception was for the analysis of 5' and 3' miRNA sequence heterogeneity, when all sequences mapping within five nucleotides of the annotated 5'-end were retained.

Novel microRNA detection

For the miRNA precursors with one annotated miRNA, small sequences mapping to the opposite arm of the precursor hairpin were analysed. The sequence with the maximal copy number was considered as a novel miRNA candidate. A further criterion required the existence of at least 14 basepairs between an annotated miRNA and a novel miRNA candidate within the precursor hairpin. The quantification process for novel miRNAs was the same as for known miRNAs.

Mapping fly sequencing data to the Drosophila melanogaster genome

Roche/454 and Solexa sequences of the Drosophila melanogaster small RNA libraries were downloaded from GEO [GSE7448] [29] and [GSE11624] [50], and mapped to the Drosophila melanogaster genome (dm3, BDGP Release 5 from UCSC) using SOAP [51], respectively. Only sequences perfectly matching the genome and with a length ranging from 16 to 28 nt were retained. The minimal length cutoff was determined by mapping scrambled sequences to ensure that the proportion of falsely mapped reads was less than 1%. Detection of novel fly miRNAs and miRNA quantification were performed exactly as for human miRNAs.

Mapping mouse sequencing data to mouse genome

Solexa sequences of mouse ESC small RNA library were downloaded from GEO [GSM314552] [28] and mapped to the mouse genome (mm9) using SOAP. Only sequences perfectly matching the genome and with a length ranging from 18 to 28 nt were retained. Detection of novel mouse miRNAs and miRNA quantification were performed exactly as for human miRNAs.

miRNA target prediction

Potential miRNA target sites were determined using the TargetScanS algorithm [6,7,52]. TargetScanS predicts conserved targets of miRNAs by searching the 3'-UTR region of all annotated human genes for the presence of 7-mer and 8-mer sequences perfectly matching the miRNA seed region and conserved among 4 vertebrate species (human, mouse, rat, and dog).

Sequence logo and statistical analyses

Sequence logos were generated using WebLogo [53]. Statistical analysis of the sequence features was performed using R [54]. Specifically, in the 5'-nucleotide prevalence analysis, the expected probability of the Binomial test was calculated based on the 5'-nucleotide base composition of all 847 annotated human miRNAs, plus 48 new miRNAs. For the purine/pyrimidine bias analysis, the nucleotide content of high- and low-expressed miRNAs was compared to the nucleotide content of high- and low-expressed miRNAs in all 253 expressed miRNA pairs and to the expressed miRNA pairs, excluding the tested group, using one-sided Wilcoxon test. Specifically, the one-sided Wilcoxon test results were: for 33 biased miRNA pairs (see Additional file 1: Table S1), when using all 253 pairs as one control, p = 0.013, p = 0.0011; using all 253 pairs and excluding the tested group (33 biased miRNAs) as one control, p = 0.0011, p = 0.00023; for 103 low-bias miRNA pairs (see Additional file 1: Table S2), using all 253 pairs as one control, p = 0.55, p = 0.56; and using all 253 pairs and excluding the tested group (low-bias miRNA pairs) as one control, p = 0.65, p = 0.64. (see Figures ​Figures2,2, ​,33 and Additional file 1: Figure S7)

For the 10 biased and 10 low-bias high-expressed miRNA pairs (see Additional file 1: Table S3), the purine/pyrimidine bias analysis was performed by comparing the purine content and pyrimidine content of the high- and low-expressed miRNAs within the pairs, due to the small sample size. The one-sided Wilcoxon test results were: for purine enrichment of 10 biased miRNA pairs, p = 0.00015; for 10 low-bias miRNA pairs, p = 0.44.

For HeLa data, the statistical results were as follows: preference for U and C at the 5'-ends of high- and low-expressed miRNAs, respectively (One-sided Binomial test, p = 8.6+10^-7 and p = 7.2+10^-8) and purine/pyrimidine bias (One sided Wilcoxon test, p < 0.05). (see Additional file 1: Table S4)

For hESCs data, the statistical results were as follows: preference for U and C at the 5'-ends of high- and low-expressed miRNAs, respectively (One-sided Binomial test, p = 1.8+10^-4 and p = 3.2+10^-7) and purine/pyrimidine bias (One sided Wilcoxon test, p < 0.05). (see Additional file 1: Table S4)

For hEBs data, the statistical results were as follows: preference for U and C at the 5'-ends of high- and low-expressed miRNAs, respectively (One-sided Binomial test, p = 1.7+10^-5 and p = 4.3+10^-8) and purine/pyrimidine bias (One sided Wilcoxon test, p < 0.05). (see Additional file 1: Table S4)

For Mouse ESC cells data, the statistical results were as follows: preference for U and C at the 5'-ends of high- and low-expressed miRNAs, respectively (One-sided Binomial test, p = 0.008185 and p = 3.664e+10^-5) and purine/pyrimidine bias (One sided Wilcoxon test, p < 0.05). (see Additional file 1: Table S4)