Nanog-dependent feedback loops regulate murine embryonic stem cell heterogeneity

Abstract

A number of key regulators of mouse embryonic stem (ES) cell identity, including the transcription factor Nanog, show strong expression fluctuations at the single-cell level. The molecular basis for these fluctuations is unknown. Here we used a genetic complementation strategy to investigate expression changes during transient periods of Nanog downregulation. Employing an integrated approach that includes high-throughput single-cell transcriptional profiling and mathematical modelling, we found that early molecular changes subsequent to Nanog loss are stochastic and reversible. However, analysis also revealed that Nanog loss severely compromises the self-sustaining feedback structure of the ES cell regulatory network. Consequently, these nascent changes soon become consolidated to committed fate decisions in the prolonged absence of Nanog. Consistent with this, we found that exogenous regulation of Nanog-dependent feedback control mechanisms produced a more homogeneous ES cell population. Taken together our results indicate that Nanog-dependent feedback loops have a role in controlling both ES cell fate decisions and population variability.

Figure 1. Quantifying the molecular effects of Nanog fluctuations

(a) The lentiviral vector construct to conditionally regulate Nanog expression levels: dLTR, deleted long-terminal repeat; FLAP, sequence element that improves transduction efficiency; rtTA, a TetOn tetracycline (doxycycline)-controlled transcriptional activator; WRE, woodchuck hepatitis virus post-transcriptional regulatory element. (b) Flow-cytometric comparison of the distribution of Nanog expression levels in wild-type Nanog GFP⁵⁴ and NanogR ES cells. In both cases, GFP levels reflect Nanog levels. (c) Experimental design. Scale bar 100 μm. (d) Effect of Nanog downregulation and rescue on protein expression levels in the ES cell TRN as measured by western blot. Full scans are given in Supplementary Fig. S1 (e) Decomposition of the extended ES cell TRN after Nanog depletion. Colours and grayscale denote relative expression levels measured by qPCR.

Figure 2. Transcriptome changes during periods of transient Nanog depletion

(a) Heat map of significant gene expression changes. (b) Mean fold expression changes for pluripotency (elements of the extended ES cell regulatory network as detailed in Ref. and updated in Supplementary Table S1), cell cycle and lineage associated gene sets. Stars indicate significance by 2 sample t-test with p-values: * < 0.05, ** < 0.01, *** < 0.0001; arrows indicate earliest time at which a significant expression change was observed (p < 0.05). Error bars show ± one standard error, n = 3. (c) Machine-learning classification of genome-wide expression patterns of 1032 pluripotent and somatic cell samples (left panel) and Nanog downregulation time course samples (right panel). The large ranges in the Pluripotency Scores and Lineage Scores in the left hand panel reflect the wide variety of cell types used to construct the classifier. The Nanog depletion time series represent the earliest stages of ES cell differentiation, and therefore, naturally show high Pluripotency Scores and low Lineage Scores. Nevertheless, a movement away from the pluripotent state is clearly detected using these classifiers; thus, highlighting their sensitivity and range. (d) Similarity matrix of samples in classification space. (e) PCA of the Nanog depletion time-course data (first two components are shown).

Figure 3. Expression changes and promoter occupancy by Nanog

(a) Hierarchical clustering of expression changes of pluripotency, lineage and cell cycle associated genes (see Supplementary Table S1). The blue bar shows the number of times each gene has been reported as a target of Nanog in six recent papers that examined promoter occupancy^,–. The green bar shows the category to which the genes belong. Pluripotency associated genes are frequently high confidence targets of Nanog. (b) Expression patterns for high confidence direct targets of Nanog. Genes were selected as high confidence Nanog targets if they were identified in at least three of six recent papers that examined Nanog target gene promoter occupancy^,–. (c) Mean fold changes for high confidence direct targets of Nanog. Expression patterns are not uniform, so mean fold changes are shown separately for those genes that were upregulated and downregulated during the time-course. Error bars show ± one standard error, n = 3. (d) The total number of Nanog target genes that changed significantly after Nanog depletion along with the number of other ES cell TRN members that also regulate that gene. Most commonly, Nanog regulates expression in concert with 1 to 5 other transcription factors.

Figure 4. Gene expression changes are regulated in a highly combinatorial manner

Rescue efficiency (see online Methods) plotted against rescue time. Genes are grouped by: (a) the total number of factors in the ES TRN which directly regulate their expression and (b) the most significant regulatory combinations (odds-ratio > 1, Fisher exact test p-value < 0.05, and 3 or more target genes). Stars highlight combinations that correspond to feedback loops in the ES cell TRN. A few combinations (Nanog + 6, Nanog + 7, All 9 factors, in panel a; DZ, DKY, NO, NDKOS, NDKORSZ, in panel b) have negative rescue efficiency at the later time-points, indicating that Nanog reintroduction resulted in further movement away (rather than back toward) to the initial state. Promoter occupancy data is from Ref.

Figure 5. Single cell gene expression patterns

(a) Heat maps of single cell expression profiles. Highly expressed genes are in red; absent (not expressed) genes are in light blue. The coloured sidebar identifies the following classes of genes: housekeeping genes (black); pluripotency associated genes (purple); cell cycle associated genes (yellow); lineage associated genes (brown). (b) Flow-cytometric analysis of Nanog expression levels during transient Nanog downregulation. (c) SVM classification of single cell expression profiles plotted in the first two principal components. Training datasets (0hrs and 36hrs) are in the left panel and the 24hrs and 36hrsR datasets are in the right panel. In both panels, the black line separates the pluripotent and lineage primed classes.

Figure 6. Feedback in the ES cell TRN

(a) The feedback-rich wild-type TRN. (b) The feedback-depleted NanogR TRN (+dox, Nanog active). (c) The total number of feedback loops that each transcription factor participates in is shown for the wild-type ES cell TRN (red) and NanogR TRN (blue). (d) Feedback centrality (as detailed in Refs. ,) for the wild-type ES cell TRN (red) and NanogR TRN (blue).

Figure 7. Cell-cell variability in wild-type and NanogR populations

(a) Heat maps of single cell expression profiles in wild-type CCE and NanogR ES cell populations. Highly expressed genes are in red; absent (not expressed) genes are in light blue. (b) SVM classification of NanogR and wild-type CCE mouse ES cells plotted in the first two principal components. The black line separates the NanogR and wild-type classes. (c) Violin plots of single cell relative expression levels in NanogR (blue) and wild-type CCE ES cells (red) for each factor in the extended ES cell TRN. Note, that these plots show expression variation but not covariance. In order to quantify overall (multivariate) variability, the median dispersion of the populations was calculated (see online Methods). Expression of Nac1 is not shown since it was not detected in sufficient numbers of cells in either NanogR or CCE ES cell populations to estimate its distribution. (d) The distance to mediancentre may be used as a test statistic to assess significant differences in overall (multivariate) variability in NanogR and wild-type CCE cells. The star indicates that the NanogR cells are significantly less variable as a population than the wild-type CCE ES cells (p ≤ 0.05 by multivariate analogue of Levene’s test). Error bars show ± one standard error, n = 66 (NanogR) and n = 77 (wild-type).