Application of massively parallel sequencing to microRNA profiling and discovery in human embryonic stem cells

Abstract

MicroRNAs (miRNAs) are emerging as important, albeit poorly characterized, regulators of biological processes. Key to further elucidation of their roles is the generation of more complete lists of their numbers and expression changes in different cell states. Here, we report a new method for surveying the expression of small RNAs, including microRNAs, using Illumina sequencing technology. We also present a set of methods for annotating sequences deriving from known miRNAs, identifying variability in mature miRNA sequences, and identifying sequences belonging to previously unidentified miRNA genes. Application of this approach to RNA from human embryonic stem cells obtained before and after their differentiation into embryoid bodies revealed the sequences and expression levels of 334 known plus 104 novel miRNA genes. One hundred seventy-one known and 23 novel microRNA sequences exhibited significant expression differences between these two developmental states. Owing to the increased number of sequence reads, these libraries represent the deepest miRNA sampling to date, spanning nearly six orders of magnitude of expression. The predicted targets of those miRNAs enriched in either sample shared common features. Included among the high-ranked predicted gene targets are those implicated in differentiation, cell cycle control, programmed cell death, and transcriptional regulation.

Figure 1.

Distributions of sequence counts for the different classes of small RNAs. (A) The box plot (left) shows the relative expression levels of sequences in each of the eight major classes of small RNAs (log₁₀ transformation) in the hESC small RNA library. miRNAs were the most highly expressed class (mean sequence count = 253, median sequence count = 9). The most abundant miRNA in this library was miR-103, which had 91,398 instances of the most common isomiR in this library and nearly 120,000 in the matched library (EB). The highest log-transformed count between the two libraries (right) for all miRNAs identified as differentially expressed (black) is roughly normal (mean = 2.32, median = 2.17), representing a tag count of 1743 and 151, respectively. The miRNAs detected in at least one of the libraries but not significantly differentially expressed are shown in light gray for comparison. There is a slight enrichment of miRNAs with lower absolute expression in this group (mean = 1.35, median = 1.24), suggesting miRNAs with higher absolute expression levels are more likely to be identified as differentially expressed. (B) Total counts for the eight classes in the box plot are summarized. They are represented as a fraction of the total sequences that had at least one perfect alignment to the human reference genome (1,631,559 total).

Figure 2.

The repertoire of isomiRs and 3′ modifications of hsa-miR-130a. A diverse variety of isomiRs were observed for many of the known and novel miRNAs. Sequences representing the miR* were also commonly observed (highlighted in purple). The reference miRNA sequence from miRBase was, in many cases, not the most frequently observed isomiR (shown in blue). In this example, the most abundant miR* did not correspond to the current miRBase entry, although within the predicted pre-miRNA structure it derives from the correct 2-nt offset position relative to the reference miRNA (far right). Although variation at the 3′ end is generally much more common than at the 5′ end, target preference of the 5′ isomiRs may differ. Sequences with evidence of 3′ additions of nucleotides (red) were common, with certain miRNAs more heavily modified than others (summarized in Supplemental Table 4). The predicted structure of the pre-miRNA is represented at the bottom in dot-bracket notation, and as a graphic (far right). Experimentally inferred Drosha and Dicer1 cleavage positions are indicated in dark and light blue arrows, respectively. Large arrows represent cleavage sites for the most abundant isomiR whereas small arrows indicate those for the other isomiRs.

Figure 3.

Clustering of over-represented Gene Ontology (GO) classes in predicted targets of differential microRNAs. Shown are heat map representations of GO (Ashburner et al. 2000) terms over-represented among predicted cooperative targets (Y-axis) of hESC-enriched miRNAs (A) and EB-enriched miRNAs (B). All genes with statistically over-represented GO annotations were included (P < 0.01, X-axis) as identified by GoStat (Beissbarth and Speed 2004). GO terms common to both sets of genes were those involved in transcriptional regulation, differentiation, and development. Those GO terms unique to hESC-enriched miRNA targets were associated with programmed cell death, response to stress, and cell motility while those unique to EB-enriched miRNA targets describe various aspects of cell proliferation regulation as well as nucleotide and nucleic acid metabolism.