MicroRNA discovery and profiling in human embryonic stem cells by deep sequencing of small RNA libraries

Abstract

We used massively parallel pyrosequencing to discover and characterize microRNAs (miRNAs) expressed in human embryonic stem cells (hESC). Sequencing of small RNA cDNA libraries derived from undifferentiated hESC and from isogenic differentiating cultures yielded a total of 425,505 high-quality sequence reads. A custom data analysis pipeline delineated expression profiles for 191 previously annotated miRNAs, 13 novel miRNAs, and 56 candidate miRNAs. Further characterization of a subset of the novel miRNAs in Dicer-knockdown hESC demonstrated Dicer-dependent expression, providing additional validation of our results. A set of 14 miRNAs (9 known and 5 novel) was noted to be expressed in undifferentiated hESC and then strongly downregulated with differentiation. Functional annotation analysis of predicted targets of these miRNAs and comparison with a null model using non-hESC-expressed miRNAs identified statistically enriched functional categories, including chromatin remodeling and lineage-specific differentiation annotations. Finally, integration of our data with genome-wide chromatin immunoprecipitation data on OCT4, SOX2, and NANOG binding sites implicates these transcription factors in the regulation of nine of the novel/candidate miRNAs identified here. Comparison of our results with those of recent deep sequencing studies in mouse and human ESC shows that most of the novel/candidate miRNAs found here were not identified in the other studies. The data indicate that hESC express a larger complement of miRNAs than previously appreciated, and they provide a resource for additional studies of miRNA regulation of hESC physiology. Disclosure of potential conflicts of interest is found at the end of this article.

Figure 1. Small RNA library generation and data analysis pipeline

(A) The small RNA library generation and sequencing scheme is shown. Small RNAs were isolated from undifferentiated H1 hESC and isogenic spontaneously differentiating cultures. Following 3’ and 5’ linker ligation, RT-PCR was performed to generate two independent cDNA libraries of small RNAs that were then used as templates for massively parallel pyrosequencing (454 sequencing). (B) The flow chart describes the data analysis pipeline. “Seqs” represent nonredundant sequences derived after collapsing multiple reads of identical sequence. The columns flanking the middle column indicate the number of sequences and reads remaining at each step of the data analysis. At the end of the pipeline, 189 sequences in the Undiff-hESC dataset and 121 in the Diff-hESC dataset met our criteria for canonical hairpin-derived sequences. Loci numbers are higher because some canonical sequences map to more than one locus in the genome. ¹The initial step of the data analysis was removal of sequences corresponding to 18 nt and 24 nt RNA markers that had been spiked into the total RNA prior to gel electrophoresis. ²Percentages of total reads from Undiff-hESC and Diff-hESC datasets that were classified into the designated categories and filtered out at each step are listed in the middle boxes.

Figure 2. Global view of known miRNAs detected in hESC sequencing datasets

The percentage of total reads for a given miRNA in Undiff-hESC or Diff-hESC reflects its relative abundance in each cell population. The 100 most abundant miRNAs are shown arranged in order of decreasing abundance in Diff-hESC (as the lower abundance miRNAs would not be visible on the graph with the same scale). Selected miRNAs present at high abundance in Undiff-hESC are identified by name. A full list of all the known miRNAs identified in this study and their relative abundance in each dataset is available in Supplemental Table 1.

Figure 3. Differential expression of known miRNAs between Undiff-hESC and Diff-hESC

(A) Expression ratios (percent of total reads in Undiff-hESC divided by percent of total reads in Diff-hESC) are shown for all known miRNAs that were detected in both Undiff-hESC and Diff-hESC datasets. Specific data pertaining to the 10 most differentially expressed miRNAs at both ends of the spectrum are displayed in the inset. (B) The absolute number of reads obtained for miRNAs that were solely detected in either Undiff-hESC or in Diff-hESC is shown.

Figure 4. qRT-PCR (TaqMan) assays of selected miRNAs found by deep sequencing to be over-expressed in H1 Undiff-hESC relative to H1 Diff-hESC

Values on the y-axis (Relative Quantification) represent the relative expression of a given miRNA in Undiff-hESC relative to Diff-hESC as measured by qRT-PCR. (A) Results of qRT-PCR in H1 hESC were consistent with those from deep sequencing in 13 out of 14 cases. The one miRNA for which over-expression in H1 Undiff-hESC was not confirmed is indicated by an asterisk. (B) The 13 miRNAs confirmed to be over-expressed in H1 Undiff-hESC showed the same pattern of over-expression in BG01 hESC.

Figure 5. Dicer-dependent expression of novel miRNAs

(A) Dicer mRNA expression (measured by qRT-PCR) in vector control hESC vs. Dicer knockdown hESC is shown. The Relative Quantification method (RQ) was used and Dicer expression in vector control hESC is arbitrarily set to 100. As shown, Dicer transcript levels are reduced by 84% in the Dicer knockdown hESC compared to vector control hESC. (B) Custom TaqMan® Small RNA Assays were used to measure the expression of the three novel miRNAs indicated in both Dicer-knockdown and vector control H1 hESC. As negative controls, three snoRNAs (which are not expected to undergo Dicer processing) were measured using TaqMan qRT-PCR assays of similar design. The degree of expression of each small RNA (snoRNA or miRNA) in the Dicer-knockdown cells was compared to that in the vector control cells, and expressed as a fold-change relative to vector control. The expression of all three novel miRNAs was diminished significantly in Dicer-knockdown cells, whereas the snoRNAs did not show such a reduction and appeared in fact to show modestly elevated expression under Dicer-knockdown conditions. Three known miRNAs, serving as positive controls, showed the expected decrease in expression in Dicer knockdown hESC relative to vector control hESC.

Figure 6. Undiff-hESC-associated miRNAs

The set of five novel (A) and nine known (B) miRNAs that were at least 4-fold over-expressed in Undiff-hESC, or represented by at least 10 reads in the Undiff-hESC (in cases where the miRNA was not detected in Diff-hESC), is shown. P-values are calculated using Fisher's exact test. For each novel miRNA, the putative hairpin secondary structure generated by RNAshapes is shown. The red letters in the hairpin indicate the mature miRNA sequence obtained in this study. Putative secondary structures for all 13 of the novel miRNAs discovered in this study are given in Supplemental Table 4.