Bimodal Genes Vary Coherently

We next used the smFISH data to determine whether the bimodal genes were correlated, which would suggest their control by a single pair of distinct cell states, or varied independently, which would suggest a multiplicity of states. The data revealed a cluster of bimodal genes that correlated with one another. Rex1, Nanog, and Esrrb displayed the strongest correlations (r = ~0.7, Figures 2C and S2B), while genes with strong overlap between modes, such as Tcl1, Lifr, Sox2, and Tet1, displayed somewhat weaker, but still significant, correlations (r = ~0.5, Figures 2C and S2B), beyond those observed between bimodal and nonbimodal genes (e.g., r = 0.2 for Rex1 and Oct4). Thus, a cell in the high- or low-expression state of one bimodal gene is likely to be in the corresponding expression state of others. Some correlations were negative: expression of the de novo methyltransferase Dnmt3b was reduced in the Rex1-high state (r = −0.46, Figure 2C). Note that cell cycle effects did not explain these correlations (Figure S2C). Together, these data suggest that bimodal genes appear to be broadly coregulated in two distinct states.

Long-tailed genes exhibited more complex relationships. Those with very large variation (CV > 1.5) were correlated with the expression state of the bimodal gene cluster, but not with one another (Figures 2C, 2D, and S2B). For example, genes like Dppa3, Tbx3, and Prdm14 burst predominantly in the Rex1-high state, but even in this state, most cells showed no transcripts of these genes (p < 0.001, see Supplemental Information for statistical analysis). Thus, it appears that these genes are expressed in infrequent, stochastic bursts that occur mainly in one of the two cellular states.

Interestingly, expression of Socs3, a negative regulator and direct target of Stat signaling (Auernhammer et al., 1999), appeared conditional on expression of its bimodally expressed receptor Lifr (note absence of Socs3 expression in low-Lifr cells in Figure S2B). Analysis of additional regulators not measured here could in principle reveal additional states or more complex distributions. Overall, however, the multidimensional mRNA distributions measured here are consistent with a simple picture based on two primary states and stochastic bursting.

The Two Primary States Exhibit Distinct DNA Methylation Profiles

These data contained an intriguing relationship between three factors involved in DNA methylation: the de novo methyltransferase Dnmt3b, the hydroxylase Tet1, which has been implicated in removing methylation (Ficz et al., 2011; Ito et al., 2010; Koh et al., 2011; Wu et al., 2011), and Prdm14, which represses expression of Dnmt3b (Grabole et al., 2013; Leitch et al., 2013; Ma et al., 2011; Yamaji et al., 2013). While Rex1 was anticorrelated with Dnmt3b expression and positively correlated with Tet1 (Figure 3A), Prdm14 showed a long-tailed distribution conditioned on the Rex1-high state (Figure 2D). Based on these relationships and the observation that methylation increases during early development (Meissner et al., 2008), we hypothesized that the Rex1-low state might exhibit increased methylation compared to the Rex1-high state.

Open in a separate window

Figure 3

The Two Rex1 States Are Differentially Methylated

(A) smFISH measurements show that Rex1 bimodality is correlated with Tet1 and anticorrelated with Dnmt3b expression.

(B) Locus-specific bisulfite sequencing of the Dazl promoter. Methylation levels shown are in the Rex1-high (top), Rex1-low (middle), and Rex1-low-to-high reverting (bottom) populations.

(C) Global levels of 5mC measured by quantitative ELISA in the Rex1-high, -low, and -low-to-high reverting cells. Data shown are mean ± SD from two independent experiments. , p < 0.05; , p < 0.01; by two-sample t test.

(D) Histogram of promoter methylation shows bimodality in the Rex1-high (top) and -low (bottom) states, as quantified by RRBS.

(E) Scatter plot of promoter methylation between Rex1-high and -low states. Each point is the methylation fraction of a single gene promoter, color-coded by the number of CpGs in that promoter. Divergence from the diagonal implies differential methylation between states. Inset: Single CpGs in the promoter of the specific gene labeled, color coded by distance from TSS; see Figure S3C for additional genes.

To test this hypothesis, we sorted Rex1-high and -low cells using the Rex1-dGFP reporter line and performed locus-specific bisulfite sequencing at known targets of methylation Dazl, Mael, and Sycp3 (Borgel et al., 2010) (Figures 3B, S3A, and S3B). These promoters exhibited 2–3 times greater methylation in Rex1-low cells compared to Rex1-high cells, indicating the two states are differentially methylated in at least some genes. In contrast, Rex1-low cells that subsequently reverted to the Rex1-high state recovered the methylation levels of Rex1-high cells, indicating that methylation is reversible. Similarly, quantitative ELISA analysis demonstrated both differential methylation and reversibility in global methylation levels (Figure 3C).

We next asked more generally which genes exhibited differential promoter methylation. We again sorted Rex1-high and -low cells and assayed DNA methylation by reduced representation bisulfite sequencing (RRBS), analyzing regions 2 kb upstream to 500 bp downstream of each ESC-expressed mRNA transcriptional start site (Meissner et al., 2008; Marks et al., 2012). The distributions of methylation levels across genes were bimodal in both Rex1 states, with the more highly methylated peak shifted to even higher methylation levels in Rex1-low cells (Figure 3D). By analyzing the shift in methylation on a gene-by-gene basis, we found that the increase in methylation in Rex1-low cells occurred predominantly through increased methylation of the promoters that were more highly methylated in Rex1-high cells (Figures 3E and S3C). Thus, the change in promoter methylation occurs in a specific subset of promoters. Furthermore, the overall methylation level of a gene was related to the number of CpGs in its promoter, such that the larger the CpG content of a promoter, the lower its methylation in both states. Not all gene promoters were well covered by RRBS. However, among those that were, several key ESC regulators including Esrrb, Tet1, and Tcl1 all showed increased levels of methylation in the Rex1-low state. Figure 3E (inset) and Figure S3C show methylation levels of individual CpGs. These results provide a view of the change in promoter methylation that occurs during transitions between the Rex1-high and -low states.

Bursty Transcription Generates Dynamic Fluctuations in Individual Cells

Evidently, cells populate two transcriptional states, each characterized by distinct methylation profiles. To understand the dynamic changes in gene expression that occur as individual cells switch between these states, we turned to time-lapse microscopy. We analyzed transcriptional reporter cell lines for Nanog and Oct4, each containing a histone 2B (H2B)-tagged fluorophore expressed under the control of the corresponding promoter (Figures S4A and S4B; see also Supplemental Information). We imaged the reporter cell lines for ~50 hr periods with 15 min intervals between frames and segmented and tracked individual cells over time in the resulting image sequences. For each cell lineage, we quantified the instantaneous reporter production rate, defined as the rate of increase of total fluorescent protein in the cell, corrected for the partitioning of fluorescent protein into daughter cells during cell division (Supplemental Information). The H2B-fluorescent protein degradation rate is negligible under these conditions (Figure S4C), enabling us to use the reporter production rate as a measure of instantaneous mRNA level. An advantage of this approach is that it provides relatively strong fluorescence signals per cell, but still enables high time-resolution analysis (Locke and Elowitz, 2009). Consistent with static smFISH distributions, the production rate distributions of the Nanog and Oct4 fluorescent reporters were bimodal and unimodal, respectively (Figure 4A).

Open in a separate window

Figure 4

Movies Reveal Transcriptional Bursting and State-Switching Dynamics in Individual Cells

(A) Distribution of Nanog and Oct4 production rates from representative movies in serum + LIF, and Gaussian fits to the components. Production rates were extracted from a total of 376 and 103 tracked cell cycles for Nanog and Oct4, respectively.

(B) Production rate distributions of individual cell lineage trees, each consisting of closely related cells descending from a single cell. Lineage trees are color-coded by the state they spend the majority of time in.

(C) Example single lineage traces exhibiting step-like changes in production rates within a state.

(D) Cell cycle phase distribution of steps within the Nanog-high state. Step occurrences are normalized by the frequencies of each cell cycle phase observed in the tracked data.

(E) Representative trace showing apparent steps from simulations under the bursty transcription model, using parameters estimated from mRNA distribution for the Nanog-high state (see Supplemental Information; see Figure S4E for simulation of Oct4 dynamics).

(F) Example traces of individual cells switching between Nanog-low and Nanog-high states.

(G) Empirical transition rates (mean ± SD) between the two Nanog states (N^Hi, Nanog-high; N^Lo, Nanog-low).

Dynamically, cells remained in either one of two distinct Nanog expression states for multiple cell cycles (Figure 4B). During these periods, expression levels varied over the full range of Nanog expression levels within each state, with no evidence for persistent substates. However, closer examination revealed fluctuations within a single state, which typically occurred in discrete steps separated by periods of steady expression (Figure 4C). Using a computational step detection algorithm (Figure S4D and Supplemental Information), we found that Nanog and Oct4 reporters exhibited 0.38 and 0.29 steps per cell cycle, respectively. These steps occurred in a cell cycle phase-dependent manner (Figure 4D), with down-steps clustered around cell division events and up-steps more broadly distributed across cell cycle phases.

Could these step-like dynamics arise simply from transcriptional bursting? To address this question, we simulated single-cell mRNA and protein traces using the bursty transcription model, with parameters determined from the NB fits of the static mRNA distributions (Supplemental Information). These simulations generated dynamic traces resembling those observed experimentally (Figures 4E and S4E). In the simulations, mRNA half-life and burst frequency determine the characteristics of detectable steps (Figure S4F); in general, steps were most prominent at low burst frequencies and short mRNA half-lives and became difficult to discriminate at high burst frequencies and long mRNA half-lives.

Step-like dynamics appear to be a natural consequence of stochastic expression, with up-steps reflecting burst-like production of mRNA and down-steps resulting from ~2-fold reduction in mRNA copy number at cell division (Figure S4G). This interpretation is consistent with the observed clustering of down-steps around cell division events and a more uniform cell cycle distribution of up-steps (Figure 4D). Because large bursts can effectively cancel out mRNA dilution at cell division events, they may appear underrepresented near cell division events. Note that most cell cycles showed no up-steps, suggesting that they are not due to increased gene dosage after chromosome replication (Brewster et al., 2014; Rosenfeld et al., 2005).

Dynamic Transitions between Cellular States

We next asked how cells transition dynamically between states. Previous work has relied on cell sorting, which can distort the signaling environments. By contrast, movies enable direct observation of switching events within a mixed cell population. Since the Nanog reporter production rate fluctuates even within a single state, we used a hidden Markov model (HMM) to classify each cell as either Nanog-high or Nanog-low at every point in time (Supplemental Information). We trained the HMM using time-series data of Nanog reporter production rates, sampled at fixed intervals across all tracked cell cycles, and used it to identify switching events and estimate switching frequencies.

Transitions from the Nanog-low to the Nanog-high state or vice versa occurred at a rate of 2.3 ± 0.25, or 7.9 ± 1.2, transitions per 100 cell cycles (mean ± SD), respectively (Figures 4F and 4G). These events did not correlate between sister cells (Table S1), consistent with independent, stochastic events. Interstate switching on average showed bigger and longer-lasting fold changes than intrastate steps in production rates (Figure S4H). Together, these results suggest that gene expression dynamics are dominated by a combination of step-like changes due to bursty transcription on shorter timescales and abrupt, apparently stochastic interstate switching events on longer timescales.

“2i” Inhibitors Modulate Bursty Transcription and State-Switching Dynamics

We next asked how gene expression heterogeneity and dynamics change in response to key perturbations. Dual inhibition of MEK and GSK3β, known as “2i,” were shown to enhance pluripotency and reduce Rex1 and Nanog heterogeneity (Marks et al., 2012; Wray et al., 2010). However, it has remained unclear how 2i affects the distribution of other heterogeneously expressed regulatory genes and what impact it has on dynamic fluctuations in gene expression.

We found that addition of 2i to serum + LIF media reduced variability in the mRNA levels of most genes (Figure 5A). In principle, this could reflect the elimination of one cellular state and/or changes in burst parameters. In 2i, the bimodal genes from Figure 2A became unimodal, suggesting that 2i suppresses one of the two cellular states (Figures 5A and S5A). In the case of Nanog, the remaining state increased its mean transcript level by ~1.5-fold, to what we term Nanog-SH (super high). Tet1, Sox2, and Tcl1 also became unimodal, but displayed an overall decrease in absolute expression. With long-tailed genes, we found that mean Dppa3 expression decreased slightly, while Prdm14 and Tbx3 became more homogenously expressed, exhibiting an increase in mean expression by ~300% and ~1,000%, respectively. These changes reflect the fact that nearly all cells were now observed to express Prdm14 and Tbx3. Thus, 2i appeared to reduce variability in most genes, either by eliminating bimodality or by increasing their burst frequency.

Open in a separate window

Figure 5

2i and DNA Methylation Modulate Bursty Transcription and State-Switching Dynamics

(A) Comparison of mRNA distributions and CV between cells grown in serum + LIF and 2i + serum + LIF. Top: For each gene, the CV in serum + LIF is plotted on the left, and the CV for 2i + serum + LIF is plotted on the right. Dnmt3b in 2i + serum + LIF is represented in gray to reflect its marginal case of poor quality of fit in both bimodal and long-tailed models. Bottom: The left half of each violin represents the mRNA distribution in serum + LIF, while the right represents 2i + serum + LIF. Each gene’s distributions are normalized by a value corresponding to the larger 95th percentile between the two treatments.

(B) Distribution of Nanog production rates from movies in 2i + serum + LIF.

(C) Empirical transition rates between the two Nanog states in the presence of 2i (N^Lo, Nanog-low; N^SH, Nanog-SH).

(D) Mixing time in each condition is estimated from autocorrelation, A(τ), of production rate ranks shown in Figure S5D, right panels. Red, Nanog-high in serum + LIF; purple, Nanog-SH in 2i + serum + LIF. Error bars: standard deviation, bootstrap method.

(E) Comparison of transcriptional heterogeneity between Dnmt TKO (black line) and the parental line (blue bars) as measured by smFISH for Rex1, Nanog, Esrrb, and SDHA. Note that for Rex1/Nanog/Essrb, there are fewer “off” cells in the leftmost bins for the TKO than WT.

(F) Rex1-dGFP distribution as measured by flow cytometry grown in serum + LIF with 5-aza or DMSO (carrier control). Time points were taken after 2, 4, and 6 days.

(G) Cells were grown in 2i + serum + LIF and subsequently replated into serum + LIF with 5-aza or DMSO (carrier control). Time points were taken after 2, 4, and 6 days. GFP levels were measured by flow cytometry.

Recently, it was shown that 2i-treated cells exhibit differentiation propensity similar to sorted Rex1-high subpopulations in embryoid body formation, suggesting they may represent similar functional states (Marks et al., 2012). We used the time-lapse movie system to compare the dynamic behavior of 2i-treated cells to that of cells in the Rex1-high subpopulation. Consistent with mRNA measurements, 2i shifted most cells into a Nanog-SH state (96% of total), characterized by ~3-fold higher median production rate compared to the Nanog-high state in serum + LIF (Figure 5B). Only a small fraction of cells showed expression overlapping with the Nanog-low state in serum + LIF at the beginning of the movies (after 6 days in treatment). Moreover, these cells switched to the Nanog-SH state at a >40-fold higher rate than the Nanog-low to Nanog-high switching rate measured in serum + LIF, with no reverse transitions observed (Figures 5C and S5B). These observations suggest that 2i increases the Nanog-low to Nanog-high switching rate and reduces or eliminates Nanog-high to Nanog-low transitions (Figure 5C).

What effect, if any, does 2i have on the dynamics of gene expression noise? Static distributions suggested that 2i increased Nanog burst frequency ~45% from 0.39 to 0.55 burst/hr using the Nanog mRNA half-life previously estimated (Table S1 in Sharova et al., 2009) and assuming no change between conditions. To analyze the effects on dynamics, we computed the “mixing time,” previously introduced to quantify the mean timescale over which a cell maintains a given expression level relative to the rest of the cell population (Sigal et al., 2006). Simulations of the bursty gene expression model showed that higher burst frequencies lead to faster mixing times, while burst size has little effect (Figure S5C). Together with smFISH measurements, this model predicted that Nanog mixing times should be faster in 2i. Qualitatively consistent with this prediction, the mixing time of Nanog production rate was reduced from 8.5 hr in Nanog-high in serum + LIF media to 1.7 in Nanog-SH cells in 2i-containing media (Figures 5D and S5D).

Together, these results indicate that 2i impacts ESC heterogeneity at three levels: First, it reduces gene expression variation in many, but not all, genes. Second, it eliminates one cell state by increasing the rate of transitions out of the Nanog-low state and inhibiting the reverse transition. Third, as shown for Nanog, 2i increases burst frequency and reduces mixing times, effectively speeding up the intrastate equilibration rate within the cell population.

smFISH Hybridization, Imaging, and Analysis

The RNA FISH protocol from Raj et al., 2008 was adapted for fixed cells in suspension. See Supplemental Experimental Procedures for details. Semiautomated dot detection and segmentation were performed using custom Matlab software. A Laplacian-of-Gaussian (LoG) Kernel was used to score potential dots across all cells. The distribution of these scores across all potential dots was thresholded by Otsu’s method to discriminate between true dots and background dots (see Figures S1A–S1D). Please see Table S2 for the smFISH probe sequences used in this study.

mRNA Distribution Fitting

The Negative Binomial Distribution is defined as

where n = number of transcripts per cell, p₀ = probability of transition from the “on” promoter state to the “off” promoter state, and r = number of bursting events per mRNA half-life. The average burst size is computed as b = (1 − p₀)/p₀. Using this model, individual mRNA distributions were fit using maximum likelihood estimation. To discriminate between unimodal and bimodal fits, two tests were used to ensure that the improvement of the fit is counterbalanced by the additional degrees of freedom from the added parameters. To be considered bimodal, a distribution was required to pass both Akaike Information Criteria (AIC) and the log-likelihood ratio test (p < 0.05).

Fluorescence Time-Lapse Microscopy and Data analysis

Reporter cells were mixed with unlabeled parental cells at 1:9 ratio and plated at a total density of 20,000 cells/cm² on glass-bottom plates (MatTek) coated with human laminin-511 (BioLamina) >12 hr before imaging. Images were acquired every 15 min for ~50 hr with daily medium change. Cells were segmented and tracked from the acquired images using our own Matlab code (see Supplemental Information for image analysis methods).

2i Perturbation and Analysis

For 2i treatment we supplemented serum + LIF media with MEK inhibitor PD0325901 at 1 μM and GSK3 inhibitor CH99021 at 3 μM. Cells grown in serum + LIF media were treated with 2i for 6 days before harvesting for smFISH assay and imaging for movies.

We thank Jordi Garcia-Ojalvo, Xiling Shen, Georg Seelig, Arjun Raj, and David Sprinzak for helpful comments on the manuscript; the Kathrin Plath Lab, the Austin Smith Lab, and RIKEN for kindly providing reporter and knockout cell lines; and the Caltech FACS Facility for assistance with cell sorting. This work was supported by the National Institutes of Health grants R01HD075605A, R01GM086793A, and P50GM068763; the Weston Havens Foundation; Human Frontiers Science Program; the Packard Foundation; a Wellcome Trust Investigators Grant to M.A.S.; and a KAUST, APART, and CIRM Fellowship to J.T. This work is funded by the Gordon and Betty Moore Foundation through Grant GBMF2809 to the Caltech Programmable Molecular Technology Initiative.