The impact of microRNAs on transcriptional heterogeneity and gene co-expression across single embryonic stem cells

Abstract

MicroRNAs act posttranscriptionally to suppress multiple target genes within a cell population. To what extent this multi-target suppression occurs in individual cells and how it impacts transcriptional heterogeneity and gene co-expression remains unknown. Here we used single-cell sequencing combined with introduction of individual microRNAs. miR-294 and let-7c were introduced into otherwise microRNA-deficient Dgcr8 knockout mouse embryonic stem cells. Both microRNAs induce suppression and correlated expression of their respective gene targets. The two microRNAs had opposing effects on transcriptional heterogeneity within the cell population, with let-7c increasing and miR-294 decreasing the heterogeneity between cells. Furthermore, let-7c promotes, whereas miR-294 suppresses, the phasing of cell cycle genes. These results show at the individual cell level how a microRNA simultaneously has impacts on its many targets and how that in turn can influence a population of cells. The findings have important implications in the understanding of how microRNAs influence the co-expression of genes and pathways, and thus ultimately cell fate.

Figure 1. Single-cell sequencing of mESC transfected with either miR-294 or let-7c.

(a) Scheme of the bioinformatics and statistical pipeline for single-cell sequence analysis. Briefly, (i) reads are aligned with TopHat software and then (ii) cells with low coverage or with no evidence of miRNA transfection are removed and (iii) finally non-expressed genes are filtered out. Remaining genes are used (iv) for downstream analysis. Details provided in material and methods section. (b) ES distribution of miR-294 target genes across miRNA transfected, Dgcr8^−/− or WT cells. (c) PCA analysis on filtered and normalized data showing individual cells colour coded by condition. The PCA analysis separates cells according to their condition. (d) Average expression of 16 pluripotency/differentiation markers in each condition.

Figure 2. Differential gene expression of mESC transfected with either miR294 or let-7c.

(a) Volcano plots showing expression fold change of mean across cells in each condition on x axis and -log₁₀(FDR) on y axis. Individual cells within a condition were treated as repeats for differential expression analysis. Significantly differentially expressed miRNA targets (FDR <10%) for miR-294 and let-7c identified in previous population-based array experiments are highlighted as black triangles and P-values for the enrichment analysis is shown in upper left corner. miR-294 and let-7c targets are similarly highlighted in Dgcr8^−/− versus WT mESC comparison. Conversely to let-7c, miR-294 is broadly expressed in mESC and its targets are significantly upregulated in miRNA-deficient Dgcr8^−/− cells versus WT mESC. (b) Identification of discriminant processes between miR-294 and let-7c-transfected cells. Pathways are sorted according to their relevance by mean of the frequency they were selected by the RFE method to correctly assign a cell to its category. Only pathways selected at least in 10% of 1,000 trials are shown.

Figure 3. Cell-cell correlations within Dgcr8^−/− and Dgcr8^−/− transfected with either miR294 or let-7c.

(a) Density plots of distances (|1–|PCC|) between pairs of cells within each condition. Top plot is for all genes, whereas bottom plots is for CellNet pluripotency regulators as defined in ref. . Top plot inset shows PCC distribution between pairs of cells within each condition. It is worth noting that miR-294 increases homogeneity when considering all genes or pluripotency genes. (b) Heatmaps showing all within condition cell–cell PCC using all genes. Cells are ordered by hierarchical clustering and subpopulation of cells identified with Dynamic tree cut method (see Methods). Above each heatmaps the percentage of cells in each of the cell cycle phases is reported for each identified sub-population of cells as defined by Cyclone. (c) Inter-cluster distance distribution between clusters (that is, subpopulations) of cells shown in b. Inter-cluster distance for a cell x was defined as its average distance from all the other cells, expect the ones in the same subpopulation as cell x. Let-7c inter-cluster distance is significantly higher than Dgcr8^−/− inter-cluster distance (two tailed Wilcoxon test P=1.94e−8), meaning that let-7c is producing more distinct and better separated sub-populations of cell compared to the ones identified in Dgcr8^−/− cells. (d) Identification of pathways driving the formation of subpopultions of cells within let-7c or Dgcr8 knockdown conditions. Pathways are sorted according the corrected P-value (that is, FDR) returned by the ANOVA. Only top five significant pathways are shown. (e) Percentage of cells in each of the cell cycle phases for miRNA-transfected and Dgcr8^−/−-deficient cells. Exact Fisher's test is used to the compare number of cells in each cell cycle stage between miRNA-transfected versus Dgcr8-knockdown cells. (f) Average expression of six pluripotency and ten differentiation markers across the three identified subpopulations of let-7c-transfected cells.

Figure 4. Gene co-expression across cells in Dgcr8^−/− and Dgcr8^−/−transfected with either miR294 or let7c.

(a) RMI for a set 33 of high-confidence miR294 targets in miR294-transfected (blue), let-7c-transfected (red) and Dgcr8^−/− cells (yellow). Histograms show permutation test for random sets of genes expressed at the same levels as the targets under each condition. Green dotted line represents RMI of the targets themselves with associated p-values. (b) Same as a, except for a set of 41 high-confidence let-7c targets. (c) RMI distribution for each of the three conditions using larger sets of miRNA targets (defined as down with addition of miRNA and predicted by Targetscan, miRanda-miRSVR or previous population array data). As power of RMI is reduced by the larger size of these target sets, a distribution of RMIs was computed by randomly extracting a subset of 10 genes with replacement 10,000 times from the corresponding list of targets. (d) Spearman's correlations among 36 cell cycle-regulated transcripts in miRNAs transfected and Dgcr8^−/− cells show an increase of cell cycle-dependent transcription in let-7c-transfected cells. Genes are grouped by cell cycle phases (squares) and ordered in the same way across conditions. Above each heatmap the average expression of each gene and its standard deviation is reported. (e) Differential RMI among 36 cell cycle-regulated transcripts (Supplementary Data set 7) in miRNAs transfected versus Dgcr8^−/− cells show an increase of cell-cycle-dependent transcription in let-7c transfected cells. Genes are grouped by cell cycle phases and their RMI value was compared in miRNA transfected cells versus Dgcr8^−/− cells. The number of genes used for each cell cycle phase is reported in the upper part of the plot. Significant changes of RMI determined by permutation test are indicated with asterisks (let-7c G1/S P=0.026, G2 P=0.04, G2/M P=3e−3; miR-294 G2/M P=0.016). (f) Differential RMI among hallmark gene sets from MSigDb in miRNAs transfected versus Dgcr8^−/− cells. For each hallmark gene set we computed RMI across the three conditions (miR-294, let-7c and Dgcr8^−/−) and then computed the RMI fold-change (miR-294 versus Dgcr8^−/− and let-7c versus Dgcr8^−/−) and its significance determined by permutation test. Only genes set with significant changes are shown.