Dynamics of embryonic stem cell differentiation inferred from single-cell transcriptomics show a series of transitions through discrete cell states

Abstract
Introduction
Results
Discussion
Materials and methods
Data availability
References
Article and author information
Metrics

Abstract

The complexity of gene regulatory networks that lead multipotent cells to acquire different cell fates makes a quantitative understanding of differentiation challenging. Using a statistical framework to analyze single-cell transcriptomics data, we infer the gene expression dynamics of early mouse embryonic stem (mES) cell differentiation, uncovering discrete transitions across nine cell states. We validate the predicted transitions across discrete states using flow cytometry. Moreover, using live-cell microscopy, we show that individual cells undergo abrupt transitions from a naïve to primed pluripotent state. Using the inferred discrete cell states to build a probabilistic model for the underlying gene regulatory network, we further predict and experimentally verify that these states have unique response to perturbations, thus defining them functionally. Our study provides a framework to infer the dynamics of differentiation from single cell transcriptomics data and to build predictive models of the gene regulatory networks that drive the sequence of cell fate decisions during development.

https://doi.org/10.7554/eLife.20487.001

Introduction

During differentiation, cells repeatedly choose between alternative fates in order to give rise to a multitude of distinct cell types. A major challenge in developmental biology is to uncover the dynamics of gene expression and the underlying gene regulatory networks that lead cells to their different fates. Given the complexity of gene regulatory networks, with their large number of components and even larger number of potential interactions between those components, building detailed predictive mathematical models is challenging. The lack of sufficient data requires a large number of assumptions to be made in order to constrain all the parameters in such models (Karr et al., 2012).

Our accompanying work on extracting cell states and the sequence of cell state transitions from gene expression data (Furchtgott et al., 2016) suggested that following the dynamics of a key set of genes was sufficient to trace these transitions, and in several instances the set of key genes that we discovered were also functionally important for lineage decisions. We asked whether we could similarly determine the suitable parameters to quantitatively describe cell state transitions during early mammalian germ layer development and build predictive mathematical models of the underlying regulatory network.

Early differentiation of pluripotent mouse embryonic stem (mES) cells, which are derived from the inner cell mass of the peri-implantation stage embryo (see pictorial summary in Figure 1—figure supplement 1A), recapitulate various aspects of in vivo germ layer differentiation (Evans and Kaufman, 1981; Keller, 2005). During this stage, both mES cells and cells in vivo express key pluripotency factors, such as Nanog, Sox2, Oct4, Klf4, Jarid2, and Esrrb, which mutually activate one another to form a pluripotency circuit (Kim et al., 2008; Young, 2011; Zhou et al., 2007). Following implantation, naïve pluripotent ES cells of the inner cell mass downregulate Klf4 and upregulate Otx2, Dnmt3a, and Dnmt3b, as they transition into ‘primed’ pluripotent cells found in the epiblast (Buecker et al., 2014; Nichols and Smith, 2009). Over the next few days of differentiation, TGF-beta signaling factors, with the aid of WNT/beta-catenin signaling, promote and inhibit the differentiation of pluripotent cells into mesendodermal (characterized by genes such as Brachyury (T), FoxA2, Mixl1 and Gsc) and ectodermal (characterized by Eras, Sez6, Stmn3, and Stmn4) cell fates, respectively (Gadue et al., 2006; Hart et al., 2002; Li et al., 2015; Lindsley et al., 2006; Tada et al., 2005; Watabe and Miyazono, 2009). Mesendodermal progenitors further differentiate into mesoderm and definitive endoderm progenitors. Mesoderm cells are usually distinguished by the expression of Gata4 and Eomes, and endoderm cells by Sox17 and FoxA2, although in mouse these genes are shared between both lineages, with differences only in their timing and level of expression (Arnold and Robertson, 2009; Kanai-Azuma et al., 2002; Kim and Ong, 2012; Lumelsky et al., 2001; Rojas et al., 2005). Along the ectodermal lineage, BMP signaling pushes ectodermal cells toward epidermis, while in the absence of BMP signaling, ectodermal cells acquire a neural fate (Wilson and Hemmati-Brivanlou, 1995). Epidermal cells are characterized by Keratins, whereas neural cells express Sox1 and Pax6 (Koch and Roop, 2004; Pevny et al., 1998; Sansom et al., 2009; Streit and Stern, 1999). The cells at the physical border between epidermal and neural cells give rise to neural crest cells (expressing Sox10, Msx2, Snai1 and Slug) in response to WNT and BMP signaling, which are often described as a fourth germ layer because of the diverse range of tissues to which they give rise (Gans and Northcutt, 1983; Knecht and Bronner-Fraser, 2002; Nicole and Chaya, 1991). Despite the detailed understanding of early embryonic development revealed by decades of work in genetics and developmental biology, a quantitative understanding of how the underlying gene regulatory network leads cells through a series of cell fate decisions has remained elusive.

We use single-cell RNA-seq to determine how gene expression patterns change as mouse embryonic stem cells differentiate into different germ-layer progenitors. We employ a Bayesian framework (Furchtgott et al., 2016) to simultaneously infer cell states, the sequence of transitions between these states, and the key sets of genes whose expression patterns provide a parameter space in which the cell states and cell state transitions are inferred. Our computational analysis, together with experimental validation using flow cytometry and live-cell imaging of a new Otx2 reporter mES cell line, suggest that cells reside in discrete states and rapidly transition from one state to another.

Using the inferred gene expression dynamics and by requiring models to replicate the existence of the observed discrete cell states, we extract probability distributions of the parameters of a model gene regulatory network. Intriguingly, requiring the model to have discrete cell states leads to the prediction that each cell state has a distinct response to perturbations by signals and changing transcription factor expression levels. We experimentally verify three distinct categories of predictions, each testing whether cells exhibit such state-dependent behavior in response to a different type of perturbation. The experimental results conclude that whether (i) Sox2 overexpression represses Oct4, (ii) Snai1 overexpression represses Oct4, and (iii) LIF and BMP promote pluripotency or differentiation into neural crest, all depend on cell state. Finally, we discuss the biological implications of our results.

Results

Acquiring single-cell transcriptomics data during early differentiation

We differentiated populations of mES cells by exposing them to one of four combinations of signaling factors and small molecules to perturb key paracrine signaling pathways involved in early mammalian patterning (Power and Tam, 1993; Tam et al., 2006): FGF, WNT, and/or TGF-beta signaling for up to five days (Figure 1A; see also Figure 1—figure supplement 1B, Materials and methods). Although cells in each population were differentiated in a monolayer culture and therefore exposed to nearly uniform conditions, we observed significant heterogeneity in the expression – as measured by immunofluorescence – of various known early germ layer marker genes (such as T, Pax6, Slug, FoxA2, and Gata4) in each population, suggesting a diversity of cell types under the same signaling conditions (Figure 1B). Further, undifferentiated pluripotent cells persisted in differentiating populations (Figure 1—source data 1, Figure 2—source data 1). Therefore, to capture the cell-to-cell variability within differentiating populations, we collected and transcriptionally profiled single cells every 24 hr over the course of five days of differentiation (Figure 1—source data 1) using a modified version of CEL-seq (Hashimshony et al., 2012). We obtained gene expression data from a total of 288 cells (Figure 1—figure supplement 1C–J; Materials and methods) with a median of 508,939 mapped reads, 48,475 transcripts and 7032 genes detected per cell. We then randomly subsampled 20,000 reads from each cell to eliminate any technical biases that may have resulted from differences in read numbers across cells (Figure 1—figure supplement 1K).

Figure 1 with 1 supplement see all

Download asset Open asset

Single-Cell Gene Expression Profiling of mESCs during early germ layer differentiation.

(A) Mouse embryonic stem cells (mESCs) were exposed to various differentiation conditions to perturb FGF, WNT, and TGF-beta signaling for up to five days of differentiation. Single cells, collected every 24 hr during differentiation, were transcriptionally profiled using CEL-Seq. (See also Figure 1—figure supplement 1B and Figure 1—source data 1). (B) Images of immunostained mESCs undergoing differentiation show cell-to-cell variability in their expression of known germ layer marker genes. (Scale bar = 100 μm).

https://doi.org/10.7554/eLife.20487.002

Figure 1—source data 1 Differentiation conditions and duration of single cells sorted into seven 96-well plates.: https://doi.org/10.7554/eLife.20487.003
Download elife-20487-fig1-data1-v2.docx

Bayesian statistical approach discovers appropriate coordinate systems to infer cell states and state transitions

One of the challenges in analyzing single-cell gene expression data is the high dimensionality of the data set and the concomitant sparsity of the data (the number of data points divided by the dimensionality is small) (Advani and Ganguli, 2016). Conventional analysis of single-cell gene expression data relies on multi-gene or multi-cell correlation estimates, such as PCA (Seurat) (Satija et al., 2015), ICA (Monocle) (Trapnell et al., 2014) and WGCNA (Li et al., 2016; Saadatpour et al., 2014) to reduce the dimensionality of expression data. However, discovering cell types and their lineage relationships using these methods has been challenging (Furchtgott et al., 2016).

In the accompanying paper, Furchtgott et al. develop a Bayesian framework that simultaneously infers (i) cell cluster identities of the cells, ${C} \equiv {c_{1}, c_{2}, \dots, c_{N}}$ ,, (ii) the sets of transitions ${T}$ between these clusters, (iii) the key sets of marker genes ${α_{i}}$ that define each cell cluster and (iv) the sets of transition genes ${β_{i}}$ that define the transitions between clusters, from single-cell gene expression data ${g_{i}}$ , by means of an iterative algorithm to determine the maximum likelihood estimates of these variables (Furchtgott et al., 2016).

Here, we employed this Bayesian framework to discover cell types and infer their lineage relationships for early mouse germ layer differentiation. We started by clustering the single-cell gene expression data for the 288 cells into 12 seed clusters ${c_{1}^{0}, c_{2}^{0}, \dots, c_{12}^{0}}$ using Seurat (Satija et al., 2015) as well as k-means (Figure 2—figure supplement 2A, B and C), restricting the analysis to transcription factors (2672 total) because of their functional role in orchestrating global gene expression (Spitz and Furlong, 2012). Seurat identifies cell clusters by performing density-based clustering on a two dimensional t-distributed Stochastic Neighbor Embedding (t-SNE) map of the gene expression data (Van der Maaten and Hinton, 2008). These clusters ${C}^{0} = {c_{1}^{0}, c_{2}^{0}, \dots, c_{12}^{0}}$ , ranging in size from 14 to 47 single cells, served as a seed for the iterative algorithm (described below).

We next considered every possible group of 3 clusters (e.g., $c_{1}^{0}, c_{2}^{0}$ and $c_{3}^{0}$ ) from a total of $^{12} C_{3} = 220$ such combinations. For each triplet of clusters, we first determined the probability that each gene $i$ was a marker gene ( $α_{i} = 1$ ), a transition gene ( $β_{i} = 1$ ) or neither ( $α_{i}, β_{i} = 0$ ) based on the distribution of their expression patterns in cells of each cluster, where ${g_{i}}$ is the single-cell gene expression data of the $i -$ th gene. Marker and transition genes are defined as follows (Figure 2—figure supplement 1A, Materials and methods,; Furchtgott et al., 2016): (i) A marker gene $i (α_{i} = 1)$ has a distribution of expression levels that is highest in one cluster, and well separated from the distribution of its expression levels in the other two clusters. Marker genes distinguish one of the clusters from the other two. (ii) A transition gene j $(β_{j} = 1)$ has a distribution of expression levels that is lowest in one cluster, and well separated from the distribution of its expression levels in the other two clusters. Each such transition gene establishes relative relationships between the three clusters (Furchtgott et al., 2016). (iii) Genes that are neither marker ( $α = 0$ ) nor transition genes ( $β = 0$ ) do not follow constraints (i) and (ii) on expression level distributions. Computing the probability of each gene being a marker gene, a transition gene, or neither allowed us to determine the most likely set of transitions $T$ between each triplet of clusters. Each gene’s contribution to the posterior probability $T$ is weighted by the odds ratio that the gene is a transition gene (Figure 2—figure supplement 1B). For example, for clusters $c_{1}^{0}, c_{2}^{0}$ and $c_{3}^{0}$ , a gene whose expression is lower in $c_{2}^{0}$ casts a vote against $c_{2}^{0}$ being the intermediate state (i.e., against the transition $T = c_{2}^{0}$ , where $c_{2}^{0}$ is intermediate, Figure 2—figure supplement 1B right) that is weighted by its odds of being a transition gene for those three clusters (Figure 2—figure supplement 1B, left). This Bayesian framework led to a summation of these weighted votes to determine the most likely set of transitions between each set of three clusters and concomitantly the most likely marker and transition genes corresponding to these clusters and transitions (Figure 2—figure supplement 1B, right).

For the seed cluster set ${C}^{0}$ , we determined 179 sets of transitions between clusters and identified 1035 transcription factors that were high probability marker or transition genes for at least one of the identified transitions. For a gene to be defined as a marker or transition gene, we used a probability cutoff of 0.5. Moreover, we used a probability cutoff of 0.6 for a triplet of clusters to count as a transition event. We next re-clustered the single cells in the gene expression space defined by these 1035 marker or transition genes using Seurat, to obtain a new cluster set ${C}^{1} = {c_{1}^{0}, c_{2}^{0}, \dots, c_{10}^{0}}$ consisting of 10 clusters. In this process, cells changed cluster identities, and certain clusters merged (Figure 2A, Figure 2—figure supplement 1C).

Figure 2 with 3 supplements see all

Download asset Open asset

Iterative Bayesian algorithm converges upon a set of cell clusters and local transitions that together define a multi-potent lineage tree.

(A) Iterative determination of the most likely sets of transitions ${T}$ and re-clustering of cells in the resulting subspace of transition and marker genes, starting from a seed set of cluster identities ${C}^{0}$ . With each iteration, the cluster identities as well as the total number of clusters change, as shown by the Seurat t-SNE maps (each dot represents a cell, colored based on its cluster identity). The inferred sets of transitions between clusters at each iteration are represented as a lineage tree (each circle represents a cell cluster). After five iterations, the algorithm converged upon a set of 9 clusters (shown in box). (See also Figure 2—figure supplement 2). (B) Left: Top ten genes (x-axis) with highest probability of being marker genes for clusters C₁ (yellow), C₂ (light red) and C₃ (light green) plotted against their probability of being marker genes. Right: Cell-cell correlation matrix computed using these 30 marker genes for the 108 cells belonging to clusters C₁, C₂ and C₃ shows three clear blocks of high correlation along the diagonal. (C) Left: Top ten genes (x-axis) with highest probability of being transitioned genes for clusters C₁, C₂ and C₃, plotted against their probability of being transitioned genes (y-axis). The transition genes belong to one of two classes, those that show high expression in cells belonging to C₁ and C₂ but low expression in C₃ (red), and those expressed at high levels in cells in clusters C₁ and C₃ but low levels in C₂ (green). The cell-cell correlation matrix computed using these 20 transition genes shows that the 29 cells belonging to cluster C₁have intermediate levels of correlation with cells in both C₂ and C₃, whereas the 46 cells in C₂ show low correlation levels with the 33 cells in C₃. (D) The global cell-cell correlation matrix computed for all 288 cells using the 889 genes used for the final iteration of clustering shows a barely detectable structure. (E) The inferred clusters and their lineage relationships can be represented in a three-dimensional coordinate system where the x- and y- axes are the normalized log expression level of the two classes of transition genes (genes in Figure 2B, left) and the z-axis measures the normalized log expression level of the marker genes for cluster C₁ (Figure 2A left in yellow). Each dot represents a single cell, and cells are colored based on their cluster identity.

https://doi.org/10.7554/eLife.20487.005

Figure 2—source data 1 Plate and well id’s of cells belonging to each cluster.: https://doi.org/10.7554/eLife.20487.006
Download elife-20487-fig2-data1-v2.docx
Figure 2—source data 2 Triplet probabilities of final tree.: https://doi.org/10.7554/eLife.20487.007
Download elife-20487-fig2-data2-v2.docx

By iteratively determining the most likely sets of transitions and the most likely marker and transition genes, and by re-clustering the cells within the subspace of these genes, the algorithm converged (i.e., the number of genes of the re-clustering subspace became less than 10% of the total number of transcription factors) upon the most likely set of cell clusters (Figure 2—source data 1), the sets of transitions between these cell clusters (Figure 2—source data 2), as well as a set of 889 genes categorized as marker or transition genes for at least one set of transitions after five iterations (Figure 2A; Figure 2—figure supplement 2D).

The final cluster set consists of 9 cell clusters ranging in size between 14 and 57 cells; every cell was mapped to a cluster, and we observed mixing of cells from different experimental conditions to the same cluster as well as cells from the same experimental conditions being assigned to different clusters (Figure 1—source data 1; Figure 2—source data 1; Figure 3—figure supplement 1A). We combined the local sets of transitions between different triplets of clusters (Figure 2—source data 2) in order to infer the most parsimonious lineage tree between the clusters (Figure 2A) (Furchtgott et al., 2016). Importantly, we obtained identical final clusters starting with different seed cluster sets using k-means clustering with the gap statistic, as well as with different threshold probability parameter values for defining transition and marker genes, showing that our results were robust to the choice of seed clusters, threshold probability value and clustering method (Figure 2—figure supplement 2A, B and C; Figure 2—figure supplement 3A and B; Materials and methods). The cluster identities as well as their lineage relationships were unchanged when the analysis was repeated with a subset of cells; in which either an entire cluster was removed or a random set of half (144) cells were removed (Figure 2—figure supplement 3C and D). Further, we found that the clustering configuration does not change depending on whether the analysis is restricted to only transcription factors or includes all genes (Figure 2—figure supplement 3E). However, using all genes resulted in greater error rates along the topology of the inferred lineage tree compared to when only transcription factors were used (Furchtgott et al., 2016; Figure 2—figure supplement 3F).

The inferred lineage relationships between the final clusters could be visualized in the subspace of inferred marker and transition genes. We illustrate this first for the three clusters C₁, C₂, and C₃. We identified three classes of marker genes, each consisting of high-probability marker genes specific to one of the three clusters (Figure 2B). Each gene class is denoted by its highest probability member gene in curly brackets (e.g., {Otx2}). When the cell-cell Pearson correlation matrix between all 288 cells was determined using the 889 genes used for the final iteration of clustering and lineage determination, the matrix showed a barely detectable structure of nine blocks (with very low contrast) along the diagonal with marginally higher correlation levels, each corresponding to a cell cluster (Figure 2D). As expected, the low level of contrast observed in Figure 2D improves dramatically when the same correlation measures are taken across cells in a triplet, using marker or transition genes for this triplet; illustrating the locally defined nature of marker and transition genes (Figure 2B, right; Figure 2C, right; Figure 2—figure supplement 2E). The same matrix computed using high-probability marker genes for clusters C₁, C₂, and C₃ (Figure 2B, left) showed three distinct blocks of high correlation along the diagonal, each corresponding to a different cluster (Figure 2B, right). Similarly, when the cell-cell correlations were measured using the two classes of inferred transition genes (Figure 2C, left), each consisting of high-probability transition genes present in C₁ and downregulated either in C₂ or in C₃, the correlation matrix showed intermediate correlation levels between C₁ and either C₂ or C₃, and low correlation levels between C₂ and C₃ (Figure 2C, right). The distribution functions of the expression levels of these transition genes in each of the three different clusters (C₁, C₂ and C₃) led to the inference that clusters C₂ and C₃ are connected via cluster C₁ with a probability of 0.83 (Figure 2E, Figure 2—figure supplement 1A and B).

We visualized the gene expression changes that characterize transitions from one cell cluster to another by plotting the cells in C₁, C₂ and C₃ in a three-dimensional gene expression subspace (Figure 2E), using as axes the mean normalized expression levels of the two transition gene classes down-regulated in C₂ or C₃ (in red and green in Figure 2C) and of the marker gene class specific to C₁ (Figure 2B in orange). These axes constitute a low-dimensional coordinate system for the inferred set of transitions between C₁, C₂ and C₃.

Similarly, the inferred transitions across all sets of three clusters (Figure 2—source data 2) together form a lineage tree (Figure 3A) that spans all nine identified cell clusters, which can be visualized in gene expression space through a series of local transition and marker gene classes (Figure 3C; Figure 3—source data 1). We next investigated the gene expression variability among cells within each cluster by performing principal component analysis (PCA) on the transcription factor gene expression for cells within each cluster. Importantly, we found that for all clusters, no principal component is statistically significant (compared to randomizations of the data; Figure 3B, Figure 3—figure supplement 1B), validating that within each inferred cluster, the cells have the same identity within the resolution of our data.

Figure 3 with 1 supplement see all

Download asset Open asset

Cells transition from one discrete state to another during differentiation.

(A) Computationally inferred cell clusters and sequence of transitions are shown in the appropriate subspace of gene expression. Each dot represents a single cell, and cells are colored based on their cluster identity. For a linear transition sequence of cell states (such as from C₀ to C₁), the transitions are represented in a two dimensional plot with the axes defined by the normalized mean log of the unique reads of genes that are most differentially regulated in the two states, while for lineage bifurcations between alternative daughter cell states, the plots are shown in three dimensions, where the x and y axes are normalized mean log unique reads of the associated set of transition genes, and the z axes are the normalized mean log unique reads of the marker genes associated with the inferred progenitor state. Labeled in parenthesis next to each cluster are the abbreviated names of the putative corresponding cell types found in vivo (Epi: epiblast; bi_Ec: bi-potent ectoderm; ME: mesendoderm; NE: neural ectoderm; NC: neural crest; M: mesoderm; DE: definitive endoderm). (B) Top: Plot of the variances of the first ten principal components of the gene expression of cells in cluster C₀. The red line is the maximum principal component variance over 1000 randomizations of the data, showing that no principal component is statistically significant. Bottom: variances of the first principal component of each cluster, normalized by the maximum principal component variance of the randomized gene expression data for the corresponding cluster. (C) A list of high probability genes that belong to the various marker and transition gene classes that define the axes of the plots in Figure 3A, each represented by one gene in curly brackets. The curly brackets contain the gene name with the highest probability for that class, and other high probability genes (as in Figure 2A and B) are listed in the table. While some of the genes are used only once, others such as *Otx2* and *Oct4* are repeatedly reused in different subspaces to describe the transition. (D) Flow cytometry analysis of cell populations sampled every 24 hr during differentiation and immunostained for nine genes (two shown at a time for each density contour plot): Klf4, Otx2, Oct4, Sox2, Slug, Pax6, FoxA2, Gata4 (each taken from a different gene class shown in Figure 3C), and T recapitulate the predicted structure and temporal ordering of transitions through discrete cell states. Axes represent the log of gene expression, normalized by the range between the minimum and maximum across each gene. Plots in pink and green represent C₂ and C₃ lineages following the split from C₁, respectively. (E) Live cell microscopy of *Otx2* reporter (mCitrine) cell line to infer the dynamics of cell state transition from C₀ to C₁. Sample images (shown) at t = 0, 6, 12, 18, and 24 hr of differentiation. Cells were terminated at approximately 25 hr into differentiation and immunostained for Nanog (ES marker gene, Figure 2—figure supplement 2A), which shows an anti-correlation between Otx2 and Nanog expression levels. (Scale bar = 100 μm) (F) Top: Time series (x-axis) traces of single-cell Otx2 (y-axis) expression dynamics taken every 15 min show that the duration of transition from Otx2-low (C₀) to Otx2-high (C₁) is approximately 4 hr, which is well within the time frame of one cell cycle (~10 hr). The end-point (t = 25 hr) Otx2 levels show a clear separation between high and low (histogram of ~200 cells shown to the right in gray), indicating that some cells have made the transition from C₀ to C₁ while others not. Each trace is colored by its relative end-point Nanog immunofluorescence intensity level. Otx2 levels are normalized by the mean level at t = 0. Bottom: Histogram (y-axis = log (cell count)) of residence durations of ~400 cells in the Otx2-low C₀ state, showing that transition times vary across multiple cell cycle lengths (time lapse length = 48 hr). Inset bar shows mean as well as upper (white) and lower quartiles of the transition durations of cells.

https://doi.org/10.7554/eLife.20487.011

Figure 3—source data 1 Probabilities of membership in marker and transition gene classes in final tree. Listed are, for direct triplets along the lineage tree, the genes with the highest probabilities of belonging to transition gene and marker gene classes, and their associated probabilities. Genes belonging to the classes shown in curly brackets in Figure 3C (probability greater than 0.5%) are shown.: https://doi.org/10.7554/eLife.20487.012
Download elife-20487-fig3-data1-v2.xlsx

The inferred dynamics of differentiation can therefore be visualized in a low-dimensional subspace of gene expression, showing that differentiation occurs through a sequence of discrete cell state transitions.

Correspondence of cell states discovered ab initio from single-cell data to known in vivo cell types

Inspection of the genes that make up the local transition and marker gene classes (Figure 3C; Figure 3—source data 1) allowed us to match clusters to embryonic cell types found in vivo that show similar gene expression.

Cluster C₀ is characterized by the high expression of pluripotency genes Oct4, Sox2, Sall1, Etv5, Jarid2, Esrrb, Klf4 and Klf5, whereas cluster C₁ has lower Jarid2, Esrrb, Klf4 and Klf5, and higher Otx2, Bptf, Cbx1 and Dnmt3a/b expression compared to cluster C₀, suggesting that clusters C₀ and C₁ correspond to naïve ES and primed epiblast pluripotent cell types, respectively (Borgel et al., 2010; Goller et al., 2008; Kim et al., 2001; Nichols and Smith, 2009; Tesar et al., 2007; Zhou et al., 2007).

Clusters C₂ and C₃, which branch out from C₁, show differential expression of pluripotency genes relative to C₁; Bptf and Cbx1 are downregulated in both C₂ and C₃, Oct4, Etv5 and Dnmt3a are downregulated in cluster C₃ but maintained in C₂, and Sox2, Otx2 and Dnmt3b are downregulated in cluster C₂ but maintained in cluster C₃. Cluster C₂ is further characterized by a high expression level of primitive streak markers Mixl1 and T (Hart et al., 2002; Tada et al., 2005), whereas cluster C₃ is characterized by Sez6, Stmn3 and Stmn4, which have recently been shown to characterize the previously elusive mammalian bi-potent ectoderm progenitor population (Li et al., 2015). Together, these patterns strongly suggest that clusters C₂ and C₃ represent mesendoderm and bi-potent ectoderm progenitor cell types, respectively.

The bi-potent ectoderm progenitor-like cluster C₃ is then followed by a lineage split into clusters C₅ and C₆. While Stmn4 is downregulated in both C₅ and C₆ compared to C₃, Sez6 is downregulated in only C₅, and Stmn3 as well as neural progenitor marker Pax6 are downregulated in C₆ but maintained in C₅. Cluster C₅ is further characterized by Smarce1 and Zic2, and cluster C₆ by Slug and Msx2, suggesting that C₅ and C₆ may be related to neural progenitor and neural crest cells, respectively (Brown and Brown, 2009; Nicole and Chaya, 1991; Vogel-Ciernia and Wood, 2014).

Cluster C₄, although similar in its expression level of Mixl1 and T to cluster C₂, shows higher expression of other primitive streak genes such as FoxA2 and Tcf3 (Merrill et al., 2004) and lower expression of Etv5. Cluster C₄ is then followed by a bifurcation between clusters C₇ and C₈. Cluster C₇ shows high expression levels of Gata4 and Snai1, indicative of its relation to mesoderm, and cluster C₈ is characterized by high FoxA2 compared to clusters C₄ and C₈, suggestive of its relation to definitive endoderm (Kim and Ong, 2012; Rojas et al., 2005). We predict that cluster C₄ represents a primed bi-potent mesendoderm cell type relative to cluster C₂ (Nakanishi et al., 2009).

Together, these results suggest that the cell clusters and sets of transitions computationally inferred from single-cell transcriptomics data correspond to known in vivo cell types and their lineage relationships.

Differentiation occurs through a series of discrete cell state transitions

The fact that gene expression in each cell cluster does not vary significantly – as measured by the relative sizes of the largest eigenvalues of the PC components of the gene expression data (or percent variance explained thereby) versus that of the same data randomly shuffled (Figure 3B; Figure 3—figure supplement 1B) – allows for genes to be sorted into a few gene classes that show highly correlated expression patterns across clusters (Figure 3C). This suggests that one can validate the inferred sequence of cell state transitions and its gene expression dynamics by measuring the expression of one gene from each class in differentiating cells over time.

In order to confirm the gene expression dynamics over the inferred sequence of cell state transitions, we assessed populations of cells for their expression levels of key transition and marker genes (each taken from a different gene class) via immunostaining and flow cytometry. We sampled mES cell populations every 24 hr during differentiation and immunostained each for Klf4, Otx2, Oct4, Sox2, Pax6, Slug, FoxA2, Gata4 and T. (Although T is not assigned to a specific gene class, it is highly expressed in the mesendoderm-like states C₂ and C₄, and it thus allows us to distinguish C₂ from the earlier epiblast-like state C₁.) The flow cytometry density contour plots shown (Figure 3D) are characterized by high-density peaks which are separated from one another by regions of low density, mirroring the discreteness of the cell states inferred from single-cell transcriptomics data. The relative locations of these high-density peaks and the time at which they appear and disappear recapitulate the inferred gene expression dynamics of the cell state transitions of the lineage tree.

During the first two days of differentiation, all cell populations downregulated Klf4 and upregulated Otx2, as shown in the first row of density contour plots in Figure 3D. This is consistent with the first observed state transition in our inferred lineage tree from the naïve ES C₀ state to the primed epiblast-like state C₁. On day three of differentiation (third column of plots in Figure 3D), Sox2 and Oct4 are asymmetrically downregulated relative to the preceding population, as is seen in mesendoderm-like state C₂ and bi-potent ectoderm-like state C₃ relative to the epiblast-like state C₁. Sox2-high, Oct4-low cells on day three are either high for Pax6 or for Slug, consistent with comparisons between the neural ectoderm-like state C₅ and neural crest-like C₆. On day four, the Pax6-high and Slug-high populations become proportionally larger as the Pax6/Slug-low population shrinks, supporting the inferred temporal ordering that C₅ and C₆ arise from the bi-potent ectoderm-like state C₃. Oct4-high, Sox2-low cells on day three of differentiation are high for T, but show two discrete levels of FoxA2, mirroring the difference between the two mesendoderm-like states C₂ (FoxA2-low) and C₄ (FoxA2-high). Further, we found that Etv5, a gene whose expression dynamics had hitherto not been implicated with early mesendodermal differentiation in mammals, was significantly downregulated from C₂ to C₄, as predicted from the single-cell gene expression data (Figure 3—figure supplement 1C). Finally, at days four and five, we observe FoxA2-high, Gata4-low and FoxA2-low, Gata4-high cell populations, which correspond to the primed mesendoderm and definitive endoderm-like states C₄ and C₈ and the mesoderm-like state C₇, respectively. We thus confirmed that differentiating cell populations recapitulate the gene expression dynamics of cell state transitions inferred from single-cell data (Figure 3A).

The observation that the majority of randomly sampled cells are found to belong to one of nine discrete cell states (both transcriptionally and at the protein level) suggests that cell state transitions occur within a relatively short timeframe compared to the amount of time cells spend within each state. We tested this hypothesis on the first cell state transition from the naïve ES C₀ state to the primed epiblast-like state C₁ (Figure 3A). To do so, we generated an Otx2-mCitrine fusion protein reporter mES cell line (Materials and methods) and observed the single-cell-resolution dynamics of Otx2 expression for up to two days (Figure 3E and F).

In agreement with our hypothesis, we observed that Otx2 levels, at the end of 24 hr of differentiation, show a bimodal distribution (Figure 3F, top), and cells tend to occupy either an Otx2-low state (corresponding to ES state C₀) or an Otx2-high state (corresponding to epiblast-like state C₁). We find that cells transition from an Otx2-low to an Otx2-high state well within the duration of a single cell cycle (mean transition duration of 4.52 hr compared to the cell-cycle length of approximately 10 hr). In contrast, cells tend to stay in either Otx2-low or -high states for up to multiple cell cycles, with a large amount of cell-to-cell variability in the residence duration (Figure 3F, bottom). Together with our results from the analysis of single-cell transcriptomics data, these observations show that cells reside in discrete states in gene expression space and correspondingly undergo abrupt state transitions.

A probabilistic model that replicates the observed discrete cell states predicts state-dependent interpretation of perturbations

Our analysis of single-cell gene expression data suggested a lineage tree composed of discrete cell states, and identified genes associated with individual cell states and transitions between them. While we predict the existence of discrete cell states based on their gene expression pattern, finding unique physiological properties that can define and distinguish their existence functionally would lend even greater support to this prediction. We therefore next sought to find properties of cell states that distinguished them functionally from one another. In order to do so, we built a predictive and testable quantitative model of the underlying gene regulatory network based on the expression patterns of the marker and transition genes.

From the 889 genes that were categorized as either marker or transition genes for all the high probability triplets, we first chose genes involved only in the triplets that fall directly along the inferred lineage tree. That is, we removed genes that were categorized as transition or marker genes for triplets consisting of ‘indirect’ lineage relationships, where at least one cell state is skipped between two cell states connected through the lineage tree. For instance, we did not consider the genes categorized as marker or transition genes only in the triplet C_0, C₁ and C₅, because C₃ is skipped between C₁ and C₅.

Since some transition genes inferred from our Bayesian analysis are re-used to infer multiple local state transitions (Figure 3C ,e.g., Oct4, Otx2), we classified transcription factors based on their distinct binarized patterns of expression across all nine cell states, with genes showing the same patterns belonging to the same gene module (Materials and methods, Figure 4—figure supplement 1A, Figure 4—source data 1). Hence, we categorized the 321 marker and transition genes involving ‘direct’ triplets along the tree into 26 gene modules, each of which showed distinct patterns of expression across the cell states. Further, because our goal was to test whether different cell states were functionally distinct (i.e., respond differently to the same signals and gene expression changes), we also noted the expression pattern of signaling factor genes belonging to FGF, WNT, LIF and BMP signaling pathways along the lineage tree (Figure 3A). These signaling factor genes constituting each of these modules were selected based on GO categories, leading to a total of 29 gene modules (Materials and methods). We denote each gene module by a representative gene in square brackets; for example, the gene module that uniquely characterizes the ES state C₀ is denoted as [Klf4] (Figure 4—source data 1 and 2).

Owing to the large number of gene modules, and consequently even larger number of potential interactions between these modules, even the simplest mathematical model would consist of hundreds of parameters. However, for most of these parameters, direct experimental measurements are not available. In order to overcome this challenge, we exploited recent developments based on renormalization group approaches to determine which parameters are relevant for the observed data (Machta et al., 2013). We adapted the seminal model of artificial neural networks, known as the Hopfield model (Fard et al., 2016; Hopfield, 1984; Maetschke and Ragan, 2014), to construct an effective gene regulatory network between the 29 gene modules. By construction, we required that this mathematical model produces the nine cell states seen in Figure 3A. We considered a network that contains direct interactions, in which each module j exerts a drive on module i, which is equal to an interaction strength $J_{i j}$ (positive or negative) multiplied by the concentration of module j. The total drive on module i is the sum of the drives from the different modules. Given our observation of discrete cell states, we further considered that the total drive on module i affects expression in a highly non-linear manner, with high gene expression for drives that exceed a critical drive $ϕ_{0}$ , and low gene expression otherwise (Figure 4—figure supplement 1B). For simplicity, we assumed that the expression of every gene module exhibits a non-linear, step-function response, when subjected to the same drive; thereby reducing the number of parameters of the model. Indeed there are numerous genes that manifest sigmoidal-like response in expression, in the presence of internal and external stimuli (Lebrecht et al., 2005; Segal and Widom, 2009). Thus the effective dynamics of expression levels $m_{i}$ of each module $i$ are given by the non-linear equation:

\frac{d m_{i}}{d t} = H (\sum_{j} J_{i j} m_{j} - ϕ_{0}) - \frac{m_{i}}{τ_{i}}

where $H$ is the Heaviside step function and $τ_{i}$ is the effective lifetime of module $i$ (Materials and methods).

We determined the set of interactions $J_{i j}$ that are consistent with the observed cell states (C₀-C₈, Figure 3A) being stable fixed points of the network. If state ${\vec{m}}^{α} = {m_{1}^{α}, \dots, m_{29}^{α}}$ with expression level $m_{i}^{α}$ in module $i$ is a stable fixed point of the network, then the interactions $J_{i j}$ must be such that the total drive on each module that is expressed in ${\vec{m}}^{α}$ is greater than the critical drive, and the total drive on each module that is not expressed in ${\vec{m}}^{α}$ is less than the critical drive:

m_{i}^{α} = 1 \Rightarrow \sum_{j} J_{i j} m_{j}^{α} \geq ϕ_{0}

m_{i}^{α} = 0 \Rightarrow \sum_{j} J_{i j} m_{j}^{α} < ϕ_{0}

Thus, for each stable state, we have 29 constraints on the possible values of $J_{i j}$ , one for each module. Given that we have nine cell states, there are 29*9 = 261 inequalities that constrain the values of the 29² = 841 different parameters, $J_{i j}$ . The problem is therefore underdetermined even for our simplified model of the underlying network, and there are an infinite number of solutions that would allow for the observed cell states to be stable.

By using a linear programming method to obtain an ensemble of 10,000 sets of $J_{i j}$ interactions (Materials and methods), each satisfying the constraint that all nine cell states are stable fixed points, we estimated the probability distribution for the 841 parameters of the model (Figure 4A and Figure 4—figure supplement 2), giving us a probabilistic model of the underlying network. We further assumed that all the possible 10,000 sets of $J_{i j}$ interactions that reproduced the nine stable cell states were equally likely, since we did not have any experimental evidence to distinguish between them.

Figure 4 with 3 supplements see all

Download asset Open asset

Quantitative modeling of the network underlying germ layer differentiation.

(A) The inferred gene regulatory network from 10,000 sampled solutions that stabilize each of the nine cell states. Each circle represents a gene module. Mean positive and negative interactions between the modules are shown in red and green, respectively, and their thickness and transparency are proportional to the absolute magnitude of the mean and the coefficient of variation (c.v.), respectively. The colored circles represent the gene modules expressed uniquely in only one of the cell states (color code matched with Figure 3A for each state). (B, C) Subsets of the network consisting of gene modules that are expressed in (and stabilize) the naïve ES C₀ state (B) and epiblast-like C₁ (C) state. As cells transition from C₀ to C₁, expression of [*Klf4*], [*Apex1*], [*Ets2*], [*Atf2*] modules is downregulated (shown in gray) while [*Hes6*] and [*Otx2*] modules are upregulated, leading to changes in the effective interactions between gene modules that are common to both C₀ and C₁ states, such as [*Sox2*] and [*Oct4*]. (D) [*Sox2*] overexpression (x-axis) plotted against the probability of [*Oct4*] downregulation (y-axis) computed over 10,000 models (Materials and methods). In the C₁ state (solid line), [*Oct4*] is downregulated in an increasing fraction of models following [*Sox2*] overexpression, while in C₀, [*Oct4*] is stable in ~96% of the models (dotted line). In order to obtain the error bars for this and subsequent predictions, we randomly sampled three subsets of 3333 from the 10,000 models. For each set we computed the mean and standard error of the proportion of models that show downregulation of *Oct4* in response to *Sox2* overexpression. (**E, F**) Subsets of the model consisting of gene modules that are expressed in the epiblast-like C₁ (E) and mesendoderm-like C₂ (F) states, and their interactions with [*Snai1*], which is not normally expressed in C₁ or C₂. As cells transition from the C₁ to C₂ state, [*Hes6*], [*Sox2*], [*Otx2*], [*Churc1*] are downregulated (shown in gray), while [*Hmga1*], [T], [*Atf2*], [*Hes1*], [*Ets2*], [*Apex1*], [*Brd7*], [*Hmgn2*] and [*Smarce1*] are upregulated, leading to changes in the effective interactions between [*Snai1*] and modules that are common to both C₁ and C₂, such as [*Oct4*]. (G) The probability of [*Oct4*] being downregulated (y-axis) as a function of [*Snai1*] overexpression (x-axis). In the C₁ state (solid line), the over expression of [*Snai1*] has no effect on [*Oct4*] levels in ~94.5% of the 10,000 models whereas in the C₂ state (dotted line), the overexpression of [*Snai1*] leads to [*Oct4*] downregulation in up to 19% of the models. (H) The C₃ state shows a downregulation of [*Oct4*] and [BMP], and upregulation of [*Tead1*], [*Apex1*], [*Pax6*], [*Smarce1*], [*Ets2*], [*Atf2*], [*Hes1*], [*Fhl1*], [Hmgn2] modules relative to C₁. (I) Cells in different states are predicted to respond differently to morphogens. Plot showing the percentage of models (y-axis) where states C₁ and C₃ (x-axis) transition to C₆ (characterized by unique marker gene module [*Msx2*]), in response to [LIF]+[BMP]. C₁ cells remain stable in response to [LIF]+[BMP] signaling in >98% of the models whereas C₃ cells are destabilized and move to the C₆ state in ~11% of the models.

https://doi.org/10.7554/eLife.20487.014

Figure 4—source data 1

Gene modules used for modeling the network.

* The [BMP] and [Aes] modules have the same binary pattern. ** The [FGF] and [WNT] modules have the same binary pattern.: https://doi.org/10.7554/eLife.20487.015
Download elife-20487-fig4-data1-v2.docx
Figure 4—source data 2 Binary expression profiles of the gene modules used for modeling the network in the 9 cell clusters.: https://doi.org/10.7554/eLife.20487.016
Download elife-20487-fig4-data2-v2.docx

We used this probabilistic model to make testable predictions as to how different cell states respond to perturbations: to see if different cell states are defined not only by their distinct transcriptional profiles, but also functionally distinct in their phenotypic responses to the same perturbations. There are a vast number of testable predictions that one could extract from our gene regulatory network model. However, given the low throughput nature of perturbation experiments, we selected three distinct probabilistic predictions, each probing different aspects of the model gene regulatory network.

First, we considered changes in the effective interaction between two gene modules as a function of cell state (i.e., how the expression level changes of one gene module affect the expression of another gene module differs across cell states due to the difference sets of gene modules present in each state). To this end, we looked at two classes of gene module pairs: (i) gene modules that are co-expressed in two mother-daughter cell states and (ii) gene modules that are never co-expressed in any cell state.

Gene modules [Sox2] and [Oct4] are highly expressed in both the ES cluster C₀ and the epiblast-like C₁ cluster, after which they are asymmetrically downregulated in the mesendoderm-like C₂ and ectoderm-like C_3. We find that for 67.5% of the 10,000 sampled solutions, [Sox2] and [Oct4] have mutually inhibitory interactions (i.e., negative coupling constants). Although both [Sox2] and [Oct4] are present together in the C₀ and C₁ states, their effective interactions are altered in different ways in each cell state by the presence of other gene modules. As cells transition from state C₀ to C₁, they downregulate gene modules [Klf4], [Atf2], [Apex1] and [Ets2], and upregulate [Hes6] and [Otx2], among others (Figure 4B and C), leading to changes in the effective interaction strength between [Sox2] and [Oct4]. By incrementally increasing [Sox2] levels relative to its base value and assessing the fraction of models that show [Oct4] downregulation, we found that [Oct4] levels are predicted to be more stable to [Sox2] overexpression in state C₀ than in C₁ (Figure 4D), thus distinguishing C₀ and C₁ functionally (Geula et al., 2015).

On the other hand, [Snai1] and [Oct4] are not expressed together in any of the nine cell states. We investigated the predicted effects of [Snai1] overexpression on [Oct4] in the epiblast-like state C₁ and mesendoderm-like state C₂, both of which normally express [Oct4] but not [Snai1]. Although [Snai1] has a negative interaction with [Oct4] in 79.2% of the models, the modules expressed in C₁ exert a greater positive drive on [Oct4] (Figure 4E and F) than those expressed in C₂. This leads to the prediction that [Oct4] is less sensitive to [Snai1] overexpression in state C₁ compared to C₂ (Figure 4G).

We next considered the effect of morphogen signals in different states. Specifically, we considered the LIF, BMP, WNT and FGF signaling pathways, which are known to play a significant role in patterning the early embryo, as well as are central to our in vitro differentiation process (Materials and methods). We grouped signaling genes by their respective pathways (defined by GO categories) and assigned each group to a module based on its average expression pattern across the nine cell states. Because WNT and FGF modules show no changes in expression across all cell states (most likely due to the large number of genes that fall into the relevant GO categories), we focused on investigating the effects of LIF and BMP signaling on cells in the epiblast-like C₁ and in the bi-potent ectoderm-like state C₃ (Figure 4H). Given an initial state C₁ or C₃, we calculated the probabilities that cells either remain in the same state or move to a different state in response to [LIF] and [BMP] (Materials and methods). Our simulations found that for ~98% of the models, cells that are initially in state C₁ either remained stabilized in C₁ or moved to state C₀ in response to [LIF] and [BMP] addition. However, in response to the same perturbation, the vast majority of cells in the C₃ state either transitioned to the neural crest-like state C₆(11.2%) or stayed in the C3 state (86.1%) (Figure 4I).

To summarize, we predict that [Oct4] expression is less sensitive to [Sox2] overexpression in state C₀ than in C₁; [Oct4] expression is less sensitive to [Snai1] overexpression in state C₁ compared to C₂; and cells in state C₃, but not in C₁, can transition to state C₆ following [LIF]+[BMP] exposure.

Importantly, we further noted that the model predictions were robust to changes in the probability cutoff for the genes we considered: although the number of gene modules changed (27 modules for a cut off of 0.7 and 24 for 0.9), we found that the models made the same qualitative predictions (Figure 4—figure supplement 3).

Thus, by categorizing genes into different modules by their expression patterns across the observed cell states, these modules provide a starting point for modeling the gene regulatory network responsible for cell fate decisions, allowing us to make predictions for how the network gives rise to distinct phenotypic responses to the same perturbation across different cell states.

Interpretation of Sox2, Snai1, and LIF+BMP are cell state dependent

We next experimentally tested the qualitative aspects of the model’s predictions of state-dependence in cells’ responses to perturbations. We first tested how cells’ Oct4 levels respond to Sox2 overexpression in the naïve ES and epiblast-like states C₀ and C₁. We transiently transfected cells with a plasmid containing a Tet-inducible bi-directional promoter, flanked by the open reading frames of Sox2 and mCerulean, which we used as a fluorescent reporter of induction (Figure 5—figure supplement 1A). We induced overexpression in cells either in the undifferentiated C₀ state or the epiblast-like C₁ state, which correspond to Day 0 and Day 2 of differentiation, respectively (Figure 3D, Figure 5—figure supplement 1D). As a control, we used identical populations that were transfected with a plasmid containing only mCerulean under the inducible promoter. In such experiments, we typically saw mCerulean fluorescence appear approximately three hours into induction and persist for about three to four days after transfection. We therefore induced overexpression for 24 hr to minimize the effect of plasmid loss but still allow for several cell cycles to occur during induction. Following induction, we fixed and immunostained the cells for Oct4, and analyzed the results via flow cytometry. In agreement with our predictions (Figure 4D), we found that Sox2 overexpression correlates ( $R = - 0.3258, p = 1.48 \times 10^{- 13}$ ) with downregulation of Oct4 in the epiblast-like state C₁ (significant relative to control, $p = 5.72 \times 10^{- 31}$ ; see also Figure 5—figure supplement 1C), whereas this effect was not observed in undifferentiated cells (state C₀) (Figure 5A and B).

Figure 5 with 1 supplement see all

Download asset Open asset

Experimental validation shows that interpretation of *Sox2*, *Snai1*, and LIF+BMP is cell state dependent.

(A) Comparison of the effects of *Sox2* overexpression (x-axis) on Oct4 levels (y-axis) in the naïve ES state C₀ (left) and epiblast-like C₁ state shows negative correlation between *Sox2* overexpression and Oct4 levels in the C₁ state, but not in C₀. Plots showing mCerulean (marker) -only overexpression in C₀ or C₁ are indistinguishable from *Sox2* overexpression in C₀ (Figure 5—figure supplement 1C). (B) Fraction of Oct4-high cells (y-axis; defined as greater than 2σ below the mean log of Oct4 of non-transfected control cells) plotted against binned *Sox2* overexpression level confirms model prediction (Figure 4D) that *Sox2* overexpression leads to downregulation of Oct4 in C₁ but not C₀. (C) Comparison of the effects of *Snai1* and mCerulean-only (left) overexpression on Oct4 levels (x-axis) in the epiblast-like C₁ and mesendoderm-like C₂ states (y-axis; T-low and –high, respectively) shows downregulation of Oct4 in response to *Snai1* overexpression in the C₁ state but not in C₂. (D) Fraction of Oct4-high cells in *Snai1* overexpressing cells, normalized by this fraction in mCerulean overexpressing control cells (y-axis), plotted against binned *Snai1* overexpression level (x-axis) confirms the prediction (Figure 4G) that *Snai1* overexpression leads to greater downregulation of Oct4 in C₂ compared to C₁. (E) Live cell images of Oct4-mCitrine cells at t = 0, 6, 12, 18, 24 hr of LIF+BMP exposure. At t = 0, cells are either in state C₁ (Oct4-high) or C₃ (Oct4-low) (Figure 5—figure supplement 1E). (Scale bar = 100 μm) Cells were fixed at t = 24 hr and immunostained for *Msx2*. (F) Time series (x-axis) traces of single-cell Oct4 expression (y-axis) taken every 15 min from live cells. Each trace is colored by its relative end-point *Msx2* immunofluorescence intensity level. (G) The initial Oct4 reporter (mCitrine) intensity (y-axis) and final Msx2 immunofluorescence (x-axis) are negatively correlated. Each dot represents a single cell. Histogram of Oct4 reporter intensity at t = 0 levels shown in gray. Based on this histogram, we defined a range of threshold values for determining Oct4-high and –low (shown in overlapping region of orange and green along y-axis). (H) Plot showing fraction of Msx2-high (y-axis; as defined by greater than 2σ above background) confirms prediction (Figure 4I) that Msx2 is upregulated with a greater probability in the C₃ state compared to C₁ (x-axis) in response to LIF+BMP exposure.

https://doi.org/10.7554/eLife.20487.020

We then tested the effects of Snai1 overexpression on Oct4 in the epiblast-like state C₁ and mesendoderm-like state C₂, using the same experimental framework as described above. On day three of differentiation, cell populations either contain a mixture of C₁, C₂ and (minimally) C₄ cell states, or a combination of C₁, C₃ and C₅ (or C₆), depending on the signaling conditions (Figure 3D). Using the signaling conditions that yield the former set of cell states (C₁, C₂ and C₄₎, we transfected cells at 2.5 days into differentiation, and drove overexpression of Snai1 12 hr later in a population consisting primarily of cells in C₁ and C₂ states (Figure 5—figure supplement 1D). After 24 hr of Snai1 overexpression and further differentiation, we fixed and immunostained the cells for T to distinguish cells in C₁ (T-low) and C₂ (T-high) states. We also immunostained the cells for Oct4 to distinguish the C₁ state from other T-low states that arise during the last 24 hr of differentiation following the initiation of induction. We found that the fraction of C₁ cells within the transfected population was significantly reduced relative to control ( $p = 1.98 \times 10^{- 13}$ ), suggesting that cells in this state had downregulated Oct4 levels in response to Snai1 overexpression. On the other hand, the fraction of C₂ cells within the transfected population and their Oct4 levels were maintained relative to control, in agreement with our predictions (Figure 4G; Figure 5C and D).

Finally, we tested whether cells in epiblast-like C₁ and bi-potent ectoderm-like C₃ states respond differently to LIF+BMP signaling, as predicted by our model. In order to investigate the relationship between a cell’s initial state and its final state in response to LIF+BMP exposure, we needed to assess cells’ initial states non-invasively. We found that 2.5 days into differentiation, we could obtain populations that consist primarily of cells in epiblast-like state C₁ and bi-potent ectoderm-like state C₃ (Figure 5—figure supplement 1E), which have high and low expression of Oct4, respectively. We therefore utilized an Oct4-mCitrine mES cell line that we had previously engineered (Thomson et al., 2011) to distinguish cells in C₁ and C₃ states after 2.5 days of differentiation. At this point, 1200 U/mL LIF and 25 ng/mL BMP4 were added to the media, after which we followed individual cells’ Oct4 expression dynamics for approximately 24 hr via live-cell microscopy, followed by fixing and immunostaining for Msx2, a unique marker gene for the neural crest-like cell state C₆ (Figure 5E and F). As predicted by the model (Figure 4I), only cells that had low Oct4 levels (and were therefore in the bi-potent ectoderm-like state C₃) prior to LIF+BMP exposure showed upregulation of Msx2 in response to LIF+BMP ( $R = - 0.5056, p = 0.0044,$ Figure 5G and H). Together, these results show that the inferred cell states reflect phenotypic discreteness in cells’ responses to perturbations, and that the gene expression changes that define these responses mirror those predicted by our model gene regulatory network.

Discussion

By using learned sparse patterns of gene expression from established experimental systems (Furchtgott et al., 2016), we can analyze single-cell transcriptomics data to uncover the gene expression dynamics of differentiation. This method naturally identifies a small set of transcription factors whose expression profiles are multimodal across neighboring cell states. Given that transcription factors are key orchestrators of gene expression and therefore cell fate decisions (Spitz and Furlong, 2012), multimodal distributions of the expression levels of even a small set of transcription factors can define cell states in a population of cells.

While cell states can be characterized by the gene expression patterns of key sets of genes, these states can only be fully validated by demonstrating distinct physiological properties. To discover distinct properties of the cell states in early mES cell differentiation, we built probabilistic models of the underlying network. Requiring these models to have discrete cell states leads to the prediction that each cell state has a distinct response to perturbations by signals and changing levels of gene expression. Thus, the cell states we discovered can be functionally defined by their responses to perturbation. Our experimental tests show, as predicted by the model network, that Oct4 is either downregulated or unaffected by overexpression of Sox2 or Snai1, depending on the cell state. Previous studies have already shown that Sox2 and Oct4, along with Klf4, constitute part of a positive feedback loop that stabilizes the pluripotent ground state (Kim et al., 2008; Young, 2011). It is also known that in undifferentiated cells, Snai1 overexpression leads to downregulation of Oct4 expression and, subsequently, to exit of pluripotency (Galvagni et al., 2015). However, our results demonstrate that these interactions are state-dependent by showing that the effective positive interactions between Sox2 and Oct4 become destabilized as Klf4 levels drop and cells transition to a primed, epiblast-like pluripotent state. Similarly, the negative interaction exerted by Snai1 on Oct4 becomes attenuated in the presence of early primitive streak genes such as T. We also predict and show that LIF+BMP exposure pushes bi-potent ectoderm-like cells toward an Msx2-positive neural crest-like state, but this effect is not seen in epiblast-like cells. These results are further supported by the fact that both LIF and BMP signaling pathways can be used to keep cells in the pluripotent cell state (Chambers, 2004; Tam et al., 2006; Ying and Smith, 2003), and that BMP signaling plays a significant role in the differentiation of neural crest cells (Knecht and Bronner-Fraser, 2002). Together, these findings signify that the inferred cell states directly reflect differences in cells’ responses to perturbations and show that these cell states can also be defined by their unique responses to perturbations.

Comprehensive interrogation of gene expression through RNA sequencing is impossible without the termination of cells, providing only static snapshots of gene expression during differentiation. Despite this and the complexity of the underlying network, we discover that both cell states and the sequence of cell state transitions can be accurately determined by monitoring the levels of just a few transition or marker genes. Monitoring the expression dynamics of these key genes in live cells using microscopy will allow us in the future to continuously track the cell-fate decisions of individual cells. The inferred gene modules therefore represent the ‘order parameters’ by which cell-state transition dynamics can be directly measured. Live cell microscopy experiments will also allow us to measure, in conjunction with cell state transition dynamics, changes in individual cells’ spatial environment, movement, lineage history, and cell cycle dynamics in order to address fundamental biological questions as to how these factors affect cell fate decisions. Finally, our results suggest that cell-to-cell heterogeneity within differentiating populations arises largely as a consequence of cells’ variability in their timing of cell state transitions. Our inferred cell clusters show mixing of cells from different time points (Figure 1—source data 1, Figure 2—source data 1), suggesting that the observed states themselves do not change over time and that at the population level, differentiation occurs as a change in the proportions of cells in various cell states rather than through changes in the cell states themselves (Figure 3D). Since cells interpret perturbations differently even in consecutive states (Figure 5), this suggests that heterogeneity arising from timing variability is further amplified in response to signal addition or fluctuations in gene expression level. These findings emphasize the importance of understanding how the timing of cell state transitions is controlled during development.

Share this article

Cite this article

Single-Cell Gene Expression Profiling of mESCs during early germ layer differentiation.

Figure 1—source data 1

Iterative Bayesian algorithm converges upon a set of cell clusters and local transitions that together define a multi-potent lineage tree.

Figure 2—source data 1

Figure 2—source data 2

Cells transition from one discrete state to another during differentiation.

Figure 3—source data 1

Quantitative modeling of the network underlying germ layer differentiation.

Figure 4—source data 1

Figure 4—source data 2

Experimental validation shows that interpretation of Sox2, Snai1, and LIF+BMP is cell state dependent.

Author details

Sumin Jang

Contribution

Contributed equally with

For correspondence

Competing interests

Sandeep Choubey

Contribution

Contributed equally with

Competing interests

Leon Furchtgott

Contribution

Competing interests

Ling-Nan Zou

Contribution

Competing interests

Adele Doyle

Contribution

Competing interests

Vilas Menon

Contribution

Competing interests

Ethan B Loew

Contribution

Competing interests

Anne-Rachel Krostag

Contribution

Competing interests

Refugio A Martinez

Contribution

Competing interests

Linda Madisen

Contribution

Competing interests

Boaz P Levi

Contribution

Competing interests

Sharad Ramanathan

Contribution

For correspondence

Competing interests

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

Categories and tags

Research organism

Further reading