Multiple cell and population-level interactions with mouse embryonic stem cell heterogeneity

Spatial and temporal accuracy of gene expression programmes is central to cell choices during differentiation. As cells grow and divide, they dilute and turnover their contents and are exposed to intrinsic and extrinsic sources of stochasticity. For cells to differentiate, gene expression programmes must be resistant to these effects, yet reliably integrate appropriate autonomous and external signals. In recent decades, predominant approaches to investigate cell choices have been molecular, with mechanistic understanding emerging from insight into regulation, molecular interactions and effects of specific regulators. Gene regulation by higher scales of organisation – cells and tissues – has been comparatively neglected, with data largely taken from ensemble measures of gene expression from dead cells. These methods lose cell context and cell diversity and cannot monitor the emergence and maintenance of gene expression differences between cells. However, advancements in live imaging and image analysis technologies now permit a more detailed investigation into these different levels of regulation.

To study the dynamics of gene expression states in cell lineages, we investigate the heterogeneity in pluripotency factor expression in mouse embryonic stem cells (mESCs). Expression of proteins such as Nanog, Rex1 (Zfp42) and Stella (Dppa3) is highly heterogeneous in mESCs (Chambers et al., 2007; Hayashi et al., 2008; Toyooka et al., 2008). For Nanog, expression is bimodal, with high and low local maxima (Chambers et al., 2007). Nanog expression relates to phenotypic behaviour, with low-expressing cells showing a tendency to differentiate and high-expressing cells tending towards self-renewal (Chambers et al., 2007; Abranches et al., 2014). Treatment of ESCs with 2i inhibitors (Ying et al., 2008) favours self-renewal, and shifts Nanog expression towards a unimodal high distribution. In culture containing serum and LIF, cells can fluctuate between high and low states (Chambers et al., 2007; Kalmar et al., 2009; Miyanari and Torres-Padilla, 2012; Abranches et al., 2014; Singer et al., 2014) making it a potentially excellent culture system for understanding the mechanisms of how gene expression differences arise between cells. Despite several studies on ‘spontaneous' fluctuations of Nanog, triggers for the spontaneous switching are not known, necessitating a more comprehensive investigation of the regulatory influences governing expression.

We propose that key regulation of pluripotency factor expression will be identifiable in the dynamic behaviour of cells and their niche. Cell cycle dynamics are intimately associated with cell fate choice in many systems (Budirahardja and Gonczy, 2009). Is ESC cell cycle behaviour a determinant of gene expression? In mESCs, treatments prolonging cell cycles do not perceptibly alter the expression of pluripotency genes such as Nanog (Li et al., 2012; Li and Kirschner, 2014). However, although early embryonic cell cycles can be highly synchronous, many eukaryotic cycles are highly heterogeneous (Brooks, 1981; Di Talia et al., 2007; Muramoto and Chubb, 2008) and, with different signalling associated with different cycle stages, cycle variability potentially provides a driver of gene expression heterogeneity. The heterogeneity of the ESC cycle has not been determined. Other sources of heterogeneity come from cell history and environment. How does past behaviour of a cell influence future gene expression choices? Different cells have different neighbours and so potentially experience different signals and mechanical triggers. Standard ensemble or static measures of gene expression do not register dynamic cell properties such as cell cycle behaviour, cell history and environmental dynamics, and perturbation experiments often confound analysis due to the complexity of molecular interactions regulating most cellular processes.

To determine the contributions of cell and population-level processes to pluripotency factor gene expression, we investigated the regulation of Nanog expression using high-content imaging of multiple generations of unperturbed mESCs. Our large-scale data approach reveals the complexity of interactions underlying Nanog expression dynamics. We identify interactions between Nanog reporter expression, cell cycle and cell density, and reveal how expression is confined into an attractor state. We address how coupling between cellular processes is modulated during the transition to the pluripotent ground state. Finally, we introduce a new technique to distinguish cell-autonomous and non-autonomous regulation of cellular choices without experimental perturbation. Our approaches are generally applicable to understanding the regulation of gene expression decisions and cell behaviour in development.

To image fluctuations in pluripotency factor gene expression along cell lineages, we used TNGA cells (Chambers et al., 2007), which have GFP inserted directly after the Nanog translational start codon. We chose a stable GFP reporter, which is ideal for observation of long-term fluctuations of gene expression within complete cell cycles and along cell lineages, appropriate for a gene expressed over 2 days and multiple cell cycles in the early mouse embryo (Chambers et al., 2003). A destabilised GFP or direct transcriptional reporter would provide reduced signal-to-noise ratios and require potentially damaging illumination, features unsuitable for quantitative long-term imaging. To facilitate cell tracking, we expressed H2B-mRFP to label nuclei (Fig. 1A). Nuclei were tracked to generate large data arrays of coordinates for mother, daughter and granddaughter cells. Coordinates were used to extract the GFP intensity per unit volume at each time point. An example lineage is shown in Fig. 1A, with the mother cell indicated by a white arrow, its daughters with yellow arrows and granddaughters with blue. We used large data sets, typically 400-800+ cell lineages per generation per condition. We captured three independent experiments, each with five to seven imaging fields of view for two complete generations. We then captured three further independent pairwise experiments, each with six or seven fields of view, comparing daughter lineages in LIF with daughters in LIF/2i.

To determine the relationship between cell cycle dynamics and pluripotency, we first characterised the basic properties of timing and heterogeneity of ESC cell cycles in LIF. Cycle time was highly heterogeneous within cell populations (Fig. 1B). Median cycle durations were 11-13 h for daughters and 12-14 h for granddaughters (Fig. 1B). The first cycle for cells after addition of 2i-containing media had a longer duration than the corresponding controls (Fig. 1C; Kolmogorov–Smirnov test, P<0.0016; see supplementary material Appendix S1 for statistics). Increased cycle time was also observed after up to five passages in 2i, relative to controls of similar age (supplementary material Fig. S1D; P<5.8×10⁻⁵), indicating a sustained reduction of doubling rate.

Variability in cell cycle durations can be used to infer general principles of cell cycle regulation (Brooks, 1981), with several mammalian tissue culture lines showing exponential distributions for differences in division time between sister cells, indicating a control step in the cycle crossed at random (the transition probability model). To test if the transition probability model applies to ESCs, we attempted different fits for frequency plots of interdivision time for sister cells, using a single exponential and a more complex function (Murphy et al., 1984) based upon the Eyring-Stover survival theory (Wullstein et al., 1980) with an environmentally sensitive parameter τ in the exponential term (Fig. 1D,E; supplementary material Appendix S1). We carried out a chi-squared goodness-of-fit test, which is independent of heavily weighted bins, for Eyring-Stover and exponential models. The test rejected the exponential in most experiments. By contrast, the Eyring-Stover fit was retained in the majority of cases (supplementary material Appendix S1). For non-related cells, neither model fits the data. Together, these results indicate that the transition probability model is a poor descriptor of the ESC cycle and, unlike other mammalian cell lines, there is likely to be more than one critical step controlling transition through the cell cycle. The more reliable Eyring-Stover fit implies a model with environmental regulation of a rate-limiting step. Environmental influences include growth factor signalling, which can be perturbed in mouse ESCs using 2i. Variability in cycle time (characterised by the coefficient of variation, CV) was slightly increased after multiple passages in 2i (2i/LIF, CV=0.29; LIF, CV=0.23). Difficulties in tracking late-passage 2i-treated cells precluded acquisition of data sets of suitable scale for curve fitting. Instead, we considered cells during their first cycle after 2i treatment. These displayed a slightly extended interdivision time (0.43 h; Fig. 1E). However, this extended interdivision time between sisters is smaller than the overall mean change in cycle times between LIF and 2i/LIF culture (2.2 h), implying that multiple cell cycle transitions are sensitive to 2i.

Cell cycle slowing is coupled to differentiation in many systems (Budirahardja and Gonczy, 2009), so it was surprising that 2i, which reverses differentiation, extended cycles. Studies of cancer stem cell models identified a slow-cycling stem cell state (Sharma et al., 2010). Is pluripotency also a slow-cycling state? To test this in unperturbed cells, we compared cycle time and median Nanog reporter expression (within a cell cycle) for both daughter (Fig. 2A-C) and granddaughter (supplementary material Fig. S1E), in standard serum/LIF culture. Relationships between variables are described using Pearson correlation coefficients, which measure the direction and strength of linear dependent relationships between different measurements.

Cycle times were correlated with Nanog reporter expression. Although low expression occurred in both short and long cell cycles, highest levels tended to be in longer cycles. The correlation was weak (r=0.13) but significant (P=0.0018) for daughter cells from three independent experiments. A similar correlation was observed for granddaughters (r=0.14, P=0.0004; supplementary material Fig. S1E) and between cycle duration and rate of change of Nanog reporter expression (r=0.14; Fig. 2B; supplementary material Table S1). Correlations were not cycle phase dependent, as Nanog reporter levels from the first 5 h of cycles gave similar correlation values to complete cycles (Fig. 2C), implying Nanog reporter expression is not correlated with cycle time specifically because longer cycles have more time to accumulate GFP. Correlations between reporter expression and cycle duration also occurred following treatment with 2i (supplementary material Table S1). In previous experiments in mESCs, artificial extension of G1 did not alter Nanog levels (Li et al., 2012) and serum level modulation showed a similar resistance of Nanog and Oct4 (Pou5f1) to loss of growth potential (Li and Kirschner, 2014). Together, these data suggest that the extended cell cycle effects we observed are a feature rather than a driver of enhanced pluripotency.

Interactions between expression and cell cycle do not occur for all pluripotency regulators. The Rex1 transcription factor is also heterogeneously expressed in mESCs. We tested whether a reporter cell line (OCRG9) with GFP inserted into the Rex1 coding sequence (Toyooka et al., 2008) would reveal connections between expression and cycle time. The observed correlation value was −0.005 (Fig. 2D; n=434, three repeats). Coherence of Nanog and Rex1 expression has been observed (Toyooka et al., 2008; Singer et al., 2014), although coherence was only partial, so our observations of differences between Rex1 and Nanog in cell cycle coupling might reflect gene-specific differences.

A recent report using TNGA cells suggested disparities between GFP and endogenous Nanog expression (Faddah et al., 2013), with poor correspondence for the GFP-negative population, although the GFP-positive population represented Nanog protein very well. We observed a moderate overall correlation between GFP and Nanog protein levels (r=0.40, n=349 cells; supplementary material Fig. S1F); however, in agreement with Faddah et al., the GFP^– population was poorly correlated with Nanog protein levels. Independently considering the GFP⁺ population elevated the correlation with Nanog protein levels (r=0.56, P=1×10⁻¹⁹, n=226 cells). Given the measurement noise inherent in comparing two different fluorescent channels, this correlation reflects a lower bound estimate. These data, together with the Faddah et al. study, indicate that the GFP⁺ population is a good measure of Nanog protein levels. To test the effect of the GFP^– population on our live imaging data, we repeated the analysis with GFP^– cells screened from the data. The correlations between expression and cycle time were 0.18 (P=0.0001) for whole cycles and 0.19 (P=5×10⁻⁵) for the first 5 h of cycles. After 2i addition, correlations were 0.22 (P=3×10⁻⁷) for full cycles and 0.13 (P=0.0004) for the first 5 h. The unchanged correlation values imply that the interactions between GFP and cell cycle in the TNGA cells are not a consequence of the GFP^– population.

It is not clear how Nanog heterogeneity relates to cell lineage and cell cycle stage and how expression dynamics relate to current expression state. To address these issues, we mapped expression and cycle times onto cell lineage data (Fig. 3A) to identify sources of stability and change.

The intergenerational relationships between cells in lineages are displayed as correlation heatmaps in Fig. 3B,C (see also supplementary material Fig. S2). Correlations between different members of lineages for cell cycle duration are shown in Fig. 3B (highest correlations in red, lowest in blue). More closely related cells had more similar cycle times, with daughter pairs and granddaughter pairs both showing strong correlations (r=0.69±0.07 and r=0.66±0.12, respectively) and dilution of this similarity down lineages. Restricting analysis to the GFP⁺ daughters gave a similar correlation (r=0.74±0.04). Environmental regulation is not clearly apparent here, as granddaughter cousins show lower correlations than granddaughter sisters, although existing at roughly the same time and place. Strongly correlated daughter cycle times were also observed for OCRG9 cells (r=0.69±0.05). Fig. 3C shows a similar analysis for the Nanog reporter. All cells within a lineage are very strongly correlated, indicating reporter expression fluctuates very slowly. As with cell cycles, closely related members of a lineage were more correlated than more distantly related cells; however, fluctuations within a lineage over two complete cycles were small, with mother-daughter pairs having high correlation values (r=0.77±0.03) and daughter and granddaughter pairs showing yet higher correlations in reporter expression (r=0.91±0.01 and r=0.86±0.05, respectively). Restricting analysis to GFP⁺ daughters also shows a very high correlation (r=0.83±0.03). Some cells fluctuated more rapidly, in agreement with earlier observations (Kalmar et al., 2009; Abranches et al., 2014), although this behaviour was rare. With a mother-daughter correlation of 0.77, lineage correlations would become undetectable after six or seven cell cycles (3 days). Repopulation of full Nanog heterogeneity by purified high reporter cells was previously shown to be complete within 6 days (Chambers et al., 2007), suggesting that fluctuation dynamics are not enhanced in fractionated populations. Strongly correlated expression in sisters was also observed for Rex1-GFP (r=0.76±0.09). The timescales indicated by these high correlations between related cells are longer than the range of fluctuation times of ∼2 cell cycles reported for cultured human cells (Sigal et al., 2006) and Dictyostelium (Muramoto et al., 2010) and are in line with recent work using different Nanog reporters in culture (Singer et al., 2014) and early mouse embryos (Xenopoulos et al., 2015).

To gain insight into the origins of the strong correlations in cell behaviour within cell lineages, we considered a simple model, using the observed correlation values between mothers and daughters. The model generates two daughters from one mother using linear combination of mother data and a Gaussian random variable along with the intergenerational correlation values known experimentally for Nanog (r^N=0.77) and cycle time (r^C=0.6) separately. Sampling pairs of values generated a correlation for Nanog between simulated sister pairs of 0.59 (supplementary material Appendix S1) and a cycle correlation of 0.36. Experimentally, these correlations are higher, with r^N=0.91 and r^C=0.69 for daughter-daughter pairs. These differences between model and data are consistent with a role of the cell environment in stabilising gene expression between generations. Alternative models are: (1) a latent property of the mother, such as reporter RNA, is inherited to both daughters where it is revealed as an enhanced correlation; or (2) mother expression at the point of division may deviate from the median, but will be closer to that of the daughters. Correlations in Nanog and cycle behaviour between sisters were not significantly affected by 2i treatment. In side-by-side experiments, r^N=0.84 for both 2i/LIF and LIF alone, and r^C=0.75 versus 0.68, respectively. These data indicate that the processes repressed by 2i, involving MAP kinase and GSK3 signalling, are not required for intergenerational stability.

Sister cells are highly correlated in expression of the Nanog reporter, but correlations fade along cell lineages. At what time in the cell cycle do these differences appear? Fluctuations of reporter expression were measured within individual cycles (Fig. 3D). The difference in expression between sisters was small and relatively stable in the first half of cycles but increased more steeply in the second half. These data suggest the first half of the cycle is dominated by maternally expressed protein and RNA, which when diluted out reveals the dynamic behaviour of each daughter. The analysis in Fig. 3D is insensitive to fluctuations of both sisters in a correlated manner. To investigate the extent to which sister fluctuations are correlated we used a bivariate mean squared displacement (MSD) analysis on Nanog intensity values, decomposing sister time series into summation (D1+D2) and difference (D1–D2) components. In the case of independent fluctuations, the summation and difference MSDs are equal. Differences between summation and difference MSDs reflect the degree to which sister intensity fluctuations are linked. Fig. 3E shows that the summation and difference components are not equal, with the difference component showing a significantly lower trajectory, indicating sister cell fluctuations are not independent. The slight curvature of the MSD plots suggests confinement, perhaps indicating a restriction on divergence between cells.

How is the directionality of expression fluctuations altered by transition to the pluripotent ground state? Do all cells induce pluripotency gene expression in 2i? Alternatively, is the transition dominated by selection, with either death or slower cell cycling of low expressers? To distinguish between these possibilities, we measured the change in reporter levels in raw difference plots (Fig. 3F) showing median changes in GFP for all individual cells. Large changes can be observed in a small percentage of cells, with potential for both up and down transitions. The percentage of down transitions was reduced in the first cell cycle after addition of 2i (Fig. 3G). These data imply the transition into ground-state pluripotency is an induction rather than a selection. Supporting this view, the positive correlation between cell cycle duration and Nanog reporter (with or without 2i) implies no growth advantage in increasing Nanog. Furthermore, cell death counts in 2i/LIF (15 cells with 822 daughter lineages) were no greater than in LIF alone (24 cell deaths with 754 lineages), implying no widespread purging of sections of the population by 2i.

How does the directionality of fluctuations relate to cell state? We observed that fluctuations in Nanog reporter exhibit a clear directionality that depends upon the level at the start of a cell cycle (Fig. 3H). We divided cell cycle data into equal quadrants, with quadrant 1 representing the beginning of the cycle and quadrant 4 the end. The change in expression from quadrant 1 to quadrant 4 shows a negative slope when plotted against starting expression (gradient=−0.20; −0.18 for GFP⁺ cells; Fig. 3H), indicating that high-expressing cells tend to reduce expression by the end of the cycle, whereas low-expressing cells tend to increase expression. This supports the view of ESC culture as an epigenetic ‘attractor' state (Huang et al., 2005; Huang, 2011). After 2i treatment, the tendency of high-expressing cells to lower their expression was reduced (gradient=0.008; −0.077 for GFP⁺). During the transition to the ground state – a new attractor – the population will be out of equilibrium and not revert to the initial mean. During this transition, Nanog reporter expression initially decreases in many cells (Fig. 3G), sometimes substantially, implying heterogeneity in the response to dedifferentiation cues.

The simple model described above suggested Nanog and cell cycle regulation by environmental factors. To investigate any local signalling effects, we compared the difference in behaviour between cells as a function of the distance between them. Fig. 4A,B show the difference in Nanog reporter between cells as a function of distance for related (red) and unrelated (blue) daughter cell pairs at the beginning (A) and end (B) of cell cycles (see also supplementary material Fig. S3A,B). There was no relationship between intercellular distance and GFP (r=0.006). The same comparison is shown in Fig. 4C,D for cycle durations, revealing no evidence for intracellular distance as a determinant of the difference in cycle duration (r=−0.023). These data also indicate that daughters that move apart quickly are no more or less similar than those remaining in close proximity. Correlations between intercellular distance and reporter levels/cycle time were also absent in cells treated with 2i. Together, these data suggest no strong environmental determinants differentiating gene expression and cycle behaviour over the length scale of a field of view (193.5 µm²).

To investigate the possibility of environmental effects over greater length scales, we compared cell behaviours between individual imaging fields of view. We compared the field of view (FOV) average cycle time against FOV average GFP intensity (Fig. 4E). The correlation between cycle time and expression was higher in FOV average data (r=0.60) than single-cell data (r=0.13, see above) for daughters (three independent experiments). For the corresponding granddaughters, the correlation value was 0.63. The repeat test (daughters only, three independent experiments) gave r=0.37. After bootstrapping the data by randomising values between fields (see supplementary material Appendix S1), the probability that the correlation value of 0.6 between Nanog/cycle duration would occur randomly was 0.018 and 0.12 for r=0.37. So the increased correlation observed between field-averaged cell cycle and Nanog reporter expression may constitute a weak effect. A recent study (Kumar et al., 2014) found that individual ESC colonies had homogeneous expression of pluripotency markers, including Nanog, which the authors interpreted as inheritance of expression states over multiple generations. The differences in magnitude of the effects in our data and the Kumar et al. study might be due to culture conditions. In our serum/LIF culture, most cells grow as a rough monolayer, although compact vertically projecting colonies are occasionally observed. A mathematical model for culture progression from single founder cells, constrained by the inheritance values measured in our study (Fig. 4F), suggests high local correlations in cell behaviour would not arise from inheritance. Although the model used high mother-daughter correlation values (r=0.77 for Nanog), the simulated FOV cell cycle-Nanog correlation at our culture densities declined to the level of the single-cell correlation by the time the simulated culture was at the cell density used experimentally (Fig. 4G). We infer that any strong local correlations in cell behaviour would be derived, in part, from local signalling, rather than purely inherited behaviour. We commenced imaging when there were 10-30 cells in a FOV. Based on the cell cycle times we measured in this study, it is unlikely that one founder cell in a FOV could generate 10-30 cells between the time of plating and the time of imaging (18 h), so our model might overestimate the inheritance component.

Cell movements are an integral part of early mouse development (Plusa et al., 2008), sorting cells out and introducing them to novel stimuli (Xenopoulos et al., 2012), implying that cell motility might regulate the heterogeneous behaviours of ESCs. We investigated motility of TNGA cells using MSD analysis of distance moved. MSD plots of TNGA cells indicate a more complex model of translocation than random walk diffusion alone. Fig. 5A shows MSD as a function of lag between time points. The fit is non-linear and slightly upwardly curving, indicating random walk with flow, perhaps resulting from cells moving into available space, with resistance to motility from increasing cell density at source. Similar trajectories were observed with LIF and LIF/2i. Although the plot suggests enhanced motility of cells in LIF compared with cells entering 2i, suggesting an effect of enhanced local pockets of cell density in 2i, no strong change in diffusion coefficient was observed with data pooled from three independent experiments (1.13 µm²/min for LIF and 0.98 µm²/min for LIF/2i). A slight flattening of the trajectory was observed in some MSD plots, perhaps an effect of the limited size of a FOV. Motility showed weak but significant negative correlations with both cell cycle and Nanog reporter (supplementary material Table S2).

An alternative metric to describe local environment is cell density. Density was calculated by measuring the level of nuclear red fluorescence within 50 pixel diameter circles centred on cell centroids. We measured correlations for each cell between density and reporter expression at each time point (supplementary material Fig. S3C,D). Peak correlation values were slightly positive, but the spread was high and apparent weak correlations non-significant. However, if we use the median cell density from whole cycles, this showed significant positive correlations with both cell cycle duration and Nanog reporter expression (Fig. 5B; supplementary material Table S2) for daughters (r^C=0.25, r^N=0.34) and granddaughters (r^C=0.18, r^N=0.28). Anti-correlations between density and motility were observed in both LIF and 2i/LIF and are likely to reflect obstructions to cell migration. Correlation between Nanog reporter and density was clear in GFP⁺ cells (r=0.26, P=1×10⁻⁸) and persisted into 2i, although the link between cell cycle and density was lost in 2i. Similar to Nanog reporter cells, density was correlated both with cycle duration (r=0.12, P=0.01) and Rex1 reporter expression in OCRG9 cells (r=0.22, P=6×10⁻⁶; Fig. 5C). The implied role of density in the regulation of ESC behaviour may parallel the anecdotal image of the ESC colony with a dense 3D mass of pluripotent cells surrounded by the flatter and more polarised differentiating cells. Although in our serum/LIF cultures the structures that the cells form are generally monolayer-like, considerable heterogeneity in cell aggregate morphology within a culture does exist.

We have developed a high-content imaging and analysis platform for parallel measurement of multiple dynamic cellular and population features of mouse ESCs, together with gene expression, using large data sets. Our analysis revealed that mESC cell cycles are highly variable in duration. Analysis of this variability indicated that cycles are regulated at multiple transition points, unlike other standard cell lines. High Nanog reporter expression is associated with longer cell cycles, and 2i, which drives pluripotency, increased both cycle duration and variability. Fluctuations in cycle duration and gene expression were slow, with closely related cells retaining very similar cycle times and expression. The expression state of the cell is a strong indicator of its future state, although high expressers tend to reduce their expression and low expressers increase their expression. The cell environment also interacts with Nanog expression and cycle behaviour. Local intercellular signalling interactions are not strong over short timescales; however, cell density emerges as a recurrent feature, for both cycle behaviour and pluripotency factor gene expression. The link between density and cycles might be a contact inhibition phenomenon. For density and gene expression, there may be similarity with the early embryo, with Nanog becoming restricted to the inner cell mass, then epiblast (Chambers et al., 2003), which will perceive a greater number of cell-cell contacts than prospective extra-embryonic tissue. Parameterised models underpin the importance of cell-cell coupling in the long-term stability of gene expression states. Analysis of expression of the Rex1 pluripotency factor indicated partially overlapping features with Nanog regulation.

Comparing our results with previous studies of Nanog dynamics identifies apparent differences, which can be explained by the enhanced scale of our data set and the different approaches used. In an early study using TNGA cells (Kalmar et al., 2009), Nanog showed fast switching between states. Although we observed some large fluctuations over timescales of cell cycles, these were infrequent. In this early study cells were flow-sorted before imaging, providing a different population context. A recent study using a destabilised reporter also observed fast fluctuations (Abranches et al., 2014). Stable GFP reporters reveal a time-integrated view of transcription, showing the combined behaviour of several bursts of a destabilised reporter. A recent study imaging endogenous pluripotency factor levels using antibodies and single-molecule RNA fluorescent in situ hybridisation (Kumar et al., 2014) suggests very low heterogeneity in closely related cells and implies long-term stable transcriptional behaviour, which parallels the stability of Nanog reporter expression in pre-implantation mouse embryos (Xenopoulos et al., 2015) and a recent in vitro study (Singer et al., 2014). Recent work also showed no obvious differences in fluctuation range or rate between LIF and LIF/2i (Abranches et al., 2014; Singer et al., 2014). Our study shows an increase in fluctuation range in 2i, which is likely to be because we measured the transition to 2i, not the 2i steady state. If the fluctuation range did not alter after 2i treatment, other mechanisms, such as selection, would be required to generate a uniform high Nanog state. We saw no such evidence of selection, with slower cell cycles and no increase in cell death after 2i treatment. The Abranches et al. study observed no bias in mitotic division time related to Nanog reporter level. Our data are consistent with this, although we observed a clear anti-correlation between reporter levels at the beginning of the cell cycle and at the end, an observation made clearer by the scale of our data set and the clarity of a stable reporter generating a time-integrated signal.

ESC cell cycle control appears more complex than in other mammalian cell culture models. Initial studies on interdivision times between sister cells revealed a single rate-limiting transition (Brooks, 1981) for several mammalian cell lines. Subsequent work showed the variability between unrelated cells in a population could be explained by two rate-limiting steps (Brooks et al., 1980) or a more complex environmentally regulated step (Murphy et al., 1984). Our data indicate that a single rate-limiting step (the transition probability model) is not sufficient to explain interdivision times of ESC sisters, and that a more complex environmentally regulated model does not fit data from non-related cells in the population. In addition, the increase in interdivision times between sisters in 2i is small compared with the overall increase in cycle time observed in 2i. Together, these data are consistent with a model in which ESC cycle progression is actively regulated at multiple phases.

Average cycle durations were increased in high Nanog reporter cells, although high variability in cycle duration was observed in all conditions. Cycle times were further increased after 2i, although longer cell cycles were still associated with cells with higher reporter expression. Differentiation is usually associated with slowing of the cycle, so observing slower cycles for a less differentiated state was initially surprising. Slow-cycling stem cell states were previously inferred in cancer biology, although differences in cycle times (Sharma et al., 2010) are more extreme. Previous studies did not observe changes in Nanog expression caused by disruption of growth potential or G1 (Li et al., 2012; Li and Kirschner, 2014). Together, these data indicate that longer cell cycles are a feature, rather than a cause, of the pluripotent state. Frequent use of system-wide regulators such as cell cycle kinases and associated networks might not be compatible with cells remaining in the metastable attractor state proposed for stem cells (Huang, 2011). Consistent with the attractor view, cells expressing high levels of the reporter tend to decrease reporter expression. One might view this as an epigenetic barrier, such as the side of one of Waddington's valleys (Waddington, 1957) or the wall of an attractor (Huang, 2011) driving reversion to a local mean. The persistence of a significant proportion of the cell population showing overall down transitions following 2i treatment suggests a probabilistic search through the new attractor landscape, not switch-like behaviour. An alternative explanation is that the different behaviours after the dedifferentiation stimulus reflect pre-existing heterogeneities in state, inferred as a source of differential induced pluripotent stem cell (iPSC) reprogramming potential (Pour et al., 2015).

We have identified a variety of cellular and population-level features coupled with Nanog fluctuations. However, most features are highly heterogeneous – our data indicate tendencies, not determinism, which raises the conjecture that ‘stemness' is unlikely to be explained, or derived, by a ‘magic bullet'. It will be interesting to see how these conclusions are borne out in other developmental contexts. Our approaches concern the central problem in developmental biology of how cells become different, and these methods are therefore expected to be generically applicable to understanding development in a wide range of systems, and ultimately provide the basis for large-scale dynamic imaging screens to identify the molecular regulators of the interactions and phenomena that have been revealed.

For imaging Nanog fluctuations, we used TNGA cells (from Austin Smith) (Chambers et al., 2007). To image Rex1 expression, we used OCRG9 cells (from Hitoshi Niwa) (Toyooka et al., 2008). To facilitate cell tracking, cells were stably transformed with a plasmid expressing H2B-mRFP from a PGK promoter. Selection of clones used an IRES-hygromycin cassette downstream of H2B-mRFP. Cells were cultured in Glasgow Minimal Essential Medium (GMEM, Gibco) with 10% fetal bovine serum (FBS) and LIF, 1 mM sodium pyruvate, non-essential amino acids (Lonza), 2 mM L-glutamate and 7.7 ppm 2-mercaptoethanol on gelatin-coated culture dishes.

For imaging, cells were plated into 8-well µslides (Ibidi) for 10-20% confluence at imaging onset, allowing up to two complete cycles to be readily tracked without problems inherent in low culture densities. At higher starting densities it was rare to obtain two complete cell generations owing to increased cycle duration and death. Plating was carried out 18 h prior to imaging. 2 h before imaging, medium was replaced. Further experiments compared cells cultured in conventional serum/LIF culture with serum/LIF and 2i (obtained from Philip Cohen, College of Life Sciences, University of Dundee). 4 h before imaging, medium was replaced with fresh conventional or 2i medium (Ying et al., 2008). For long-term 2i treatment, H2B-RFP TNGA cells were co-cultured at a 50:50 ratio with parental TNGA cells for several passages prior to imaging, to facilitate cell tracking. We used a widefield fluorescence system designed for fast imaging of photosensitive samples (Stevense et al., 2010; Corrigan and Chubb, 2014). Images were captured using a GFP/mCherry filter set (Chroma 59022), 40×1.30 NA objective, UV (GG420, Schott) and neutral density (Chroma) filters to attenuate illumination. Bleedthrough from red into the GFP channel was corrected for post-imaging (see below). For the OCRG9 cells, which express Oct4-CFP in addition to GFP, CFP bleedthrough was less than 10% of signal, so of negligible effect on measured correlations. Fifty-two z-slices were acquired with 0.78 μm step size and 50-150 ms exposure per channel. Stacks were collected every 15 min for up to 72 h and 12-14 fields were collected in parallel using a motorised xy stage. Environmental control of CO₂, temperature and humidity was controlled with a perspex chamber (Digital Pixel) in a temperature-controlled room. Three independent repeats were carried out for two generation studies, three for comparisons of LIF with LIF/2i and Rex1 studies, and two repeats were used for late-passage 2i studies.

Movies were deconvolved using Volocity (PerkinElmer) with calculated point spread functions (PSFs). FOVs processed without deconvolution or deconvolved using measured PSFs gave similar correlation values. A graphical user interface was developed in MATLAB (Mathworks) to record cell tracks. Cell positions were recorded using a mouse click on visually determined nuclear centroids. Lineages of all cells initially in a FOV were tracked and coordinates stored in MATLAB arrays. Cells lost/dying during tracking were excluded as they could not contribute to cycle durations. Manual tracking is required to accurately follow lineages for two complete cell generations. Tracking was performed by multiple individuals, with cross-checking for reproducibility. Recorded coordinates were used to calculate reporter, H2B-mRFP and background intensity in boxes of 5×5×3 voxels centred on mouse click coordinates. Five frames either side of mitosis were removed from fluorescence data due to mitotic shape convolution effects. Compensation for bleedthrough used a custom-built function, which also subtracted the background. For FOV calculations, we also used another method for background correction with equivalent results. A description of this, and other mathematical treatment of the data, can be found in supplementary material Appendix S1. Unless otherwise described, P-values on correlation coefficients were calculated using the MATLAB function ‘corrcoef’ (applying Fisher's z-transform followed by a t-test). In box plots, we used the MATLAB default boxplot function. Horizontal red lines indicate the median, and the top and bottom of the box, respectively, denote the upper quartile (q3) and lower quartile (q1). Points are considered as outliers if their value is greater than q3+1.5×(q3–q1) or less than q1−1.5×(q3–q1). Whiskers extend to the most extreme value that is not labelled as an outlier. The threshold between GFP+ and GFP– cells was determined as the position of the minimum of the intensity histogram between the GFP– peak and the GFP+ peak.

We thank Marios Stavridis for advice throughout the project.