In embryonic development, cells differentiate through stereotypical sequences of intermediate states to generate particular mature fates. By contrast, driving differentiation by ectopically expressing terminal transcription factors (direct programming) can generate similar fates by alternative routes. How differentiation in direct programming relates to embryonic differentiation is unclear. We applied single-cell RNA sequencing to compare two motor neuron differentiation protocols: a standard protocol approximating the embryonic lineage, and a direct programming method. Both initially undergo similar early neural commitment. Later, the direct programming path diverges into a novel transitional state rather than following the expected embryonic spinal intermediates. The novel state in direct programming has specific and uncharacteristic gene expression. It forms a loop in gene expression space that converges separately onto the same final motor neuron state as the standard path. Despite their different developmental histories, motor neurons from both protocols structurally, functionally, and transcriptionally resemble motor neurons isolated from embryos.https://doi.org/10.7554/eLife.26945.001
Embryonic development proceeds through defined intermediate states, such as germ layer intermediates, and lineage-specific progenitors. Intermediates ramify into multiple states over time, and specialize their behaviors, ultimately producing a lineage tree that defines each mature cell type by a particular sequence of intermediates. This was first appreciated through classical lineage tracing and cell ablation studies. These studies showed that specifically labeled intermediate states generate stereotyped sets of downstream cell types, and that these downstream cell types fail to form if an intermediate that is upstream in their lineage is ablated (Takebayashi et al., 2002; Kretzschmar and Watt, 2012). Furthermore, most mature cell states are not produced by multiple differentiation paths in the embryo; the only known exceptions are rare and viewed as special cases, such as in the neural crest lineage.
In contrast to this rigid and hierarchical process, recent protocols that experimentally directly program cell fate suggest that the exact sequence of intermediates defining a lineage may be more flexible (Mazzoni et al., 2013; Li and Kirschner, 2014; Son et al., 2011; Vierbuchen et al., 2010; Ieda et al., 2010; Szabo et al., 2010; Cohen and Melton, 2011; Zhou and Melton, 2008; Treutlein et al., 2016). These studies reveal that mature cell states can be reached through paths that do not involve activation of the intermediate progenitor genes that are essential in embryos. Mouse embryonic stem cells (mESCs), for example, can be converted into motor neurons (MNs) by a process that involves overexpression of three transcription factors, Ngn2+Isl1+Lhx3 (Mazzoni et al., 2013; Velasco et al., 2017), and that never expresses the neural progenitor transcription factors Sox1 and Olig2 (Mazzoni et al., 2013). mESCs can also be driven rapidly into a terminal muscle phenotype without normal upregulation of intermediate genes such as Pax7 and Myf5, through a combination of cell-cycle inhibition and MyoD overexpression (Li and Kirschner, 2014). This plasticity of differentiation extends further, to the interconversion of mature cell states. Fibroblasts can be converted into mature neuron phenotypes (Vierbuchen et al., 2010; Treutlein et al., 2016), including MNs (Son et al., 2011), seemingly without completely dedifferentiating and retracing the embryonic lineage, as indicated by lack of expression of specific core pluripotency (Oct4, Sox2 and Nanog) and neural progenitor genes (Nestin) (Son et al., 2011).
Although these direct programming (DP) experiments imply the existence of differentiation paths that differ from those in embryos, much of what actually occurs in these new programs remains mysterious. Does DP bypass normal intermediates by short-circuiting the natural lineage, or does it transition through alternative intermediates (Figure 1A)? Does it diverge only briefly to bypass specific early or late states, or does it utilize an entirely distinct path (Figure 1B)? And can DP converge fully to the same final state that is produced in embryos despite taking an alternative path, or does it simply approximate that state (Figure 1C)? These questions have been challenging to answer in part due to the high degree of heterogeneity in direct programming experiments, where unbiased bulk measurements of global gene expression obscure changes, and also because marker genes allowing the isolation of new but potentially important DP-specific intermediates are not known in advance. Here we aimed to overcome these issues by applying single cell RNA sequencing to compare the gene expression trajectories of DP and growth factor guided differentiation of mESCs into MNs over time. Our core research questions are summarized in Figure 1.
We compared two in vitro differentiation protocols that convert mESCs into spinal MNs, with the goal of identifying and comparing alternate differentiation trajectories between the same start state and similar – if not the same – end states. The first, standard protocol (SP) is a widely used method that attempts to drive cells into a spinal MN state through sequential addition of developmental signals that normally guide MN differentiation in the embryo (Fgfs, Retinoic Acid, and Sonic hedgehog)(Figure 1A) (Wichterle et al., 2002; Wu et al., 2012). Many variations of this method exist, including protocols that use 2D or 3D culture formats, serum or serum-free basal culture medias, and even partially different cocktails of growth factors. Yet, each variant is believed to shuttle cells through the same sequences of intermediate states (Sances et al., 2016), namely those that exist in the embryo. As a contrast to SP, we used an alternative direct programming (DP) strategy as our second method. The DP protocol drives cells toward a spinal MN fate by forcing the expression of spinal MN transcription factors (Ngn2+Isl1+Lhx3) (Mazzoni et al., 2013; Velasco et al., 2017) in cells grown in a simple growth factor free medium (Li and Kirschner, 2014). Previous studies have shown that this DP protocol does not induce key marker genes of some embryonic MN intermediate states; for this reason DP was among the first differentiation protocols for which it was explicitly argued that a distinct differentiation process may be occurring (Mazzoni et al., 2013; Velasco et al., 2017). Our focus in this study was to understand how a non-embryonic differentiation trajectory differs from the SP, which drives cells through an embryo-like sequence of states (Figure 1).
We used single-cell RNA sequencing (InDrops) (Klein et al., 2015) to profile the differentiation in both protocols over time (Figure 1B and C). Single-cell transcriptomics was our method of choice because of its ability to parse cell states and differentiation trajectories within populations that are not pure or contain rare intermediates (Treutlein et al., 2016; Treutlein et al., 2014; Scialdone et al., 2016; Trapnell et al., 2014; Setty et al., 2016; Bendall et al., 2014), as expected for both DP and SP. We initially profiled 4590 single-cell transcriptomes sampled from early (day 4/5) and late (day 11/12) timepoints for each protocol, and also used our previously published data from 975 mES cells as a day 0 reference (Klein et al., 2015). To visualize the combined single cell time course data for each protocol, we followed the same procedure described in Klein et al (Klein et al., 2015).: we identified highly-variable genes in each data set, retained non-noise dimensions after a principal component analysis on the standardized gene expression scores, and then employed t-distributed stochastic neighbor embedding (tSNE) (van der Maaten and Hinton, 2008) to visualize the data (Figure 2B–C; further details in Supplementary Methods). The most immediate feature revealed by tSNE was a gene expression continuum for each protocol, which correlated with chronology, suggestive of a differentiation trajectory (Figure 2B–C, inset). Several disconnected cell clusters were also visible outside of the main trajectory.
By examining gene expression on the continua and disconnected cell clusters, we annotated the observed cell states produced in DP and SP. We applied unsupervised density gradient clustering (DBSCAN) to each tSNE map to partition cells into clusters, as described in Macosko et al (Macosko et al., 2015). This approach robustly partitions each of the disconnected cell clusters, but can also identify boundaries between states within continua, although the latter are less reliably defined. The cells states identified in this manner possessed unique transcriptional signatures, including expression of specific transcription factors, signaling molecules, and structural genes such as synaptic components and neurotransmitters. To annotate each cell state cluster, we identified marker genes for each subpopulation, and used prior knowledge about the characteristics of these genes to interpret the identity of each state. The criteria used for each state are provided in Figure 2—source datas 1–2, and the expression of specific marker genes is shown in Figure 2D–E and Figure 2—figure supplements 1–2.
These cell state annotations reinforced our interpretation of the main gene expression continua in DP and SP as representing differentiation trajectories. In addition to correlating with chronology, the ordering of cell states appeared sensible in both cases. Each trajectory begins with mESCs (Pou5f1+/Esrrb+/Nanog+), passes through familiar neural progenitor states (including for example neural progenitors (NPs) expressing Sox1, Pax6, and Lin28), progresses to early motor neurons (EMNs) (expressing Mnx1, Isl1, Isl2, Nefl, Nefm, Tubb3, and Map2), and terminates at a mature, or late, MN state (LMN) (downregulating MN progenitor genes Mnx1, Isl1, and Isl2, while upregulating genes indicative of neuronal maturation Robo3, Onecut2, Syn1, Vacht, and Gabr1) (Figure 2B–E). Cells along both trajectories overlap between consecutive time points, which would be expected due to asynchrony in differentiation. This was reassuring to us that our experiments sampled all intermediate states. Below we define and compare the specific cell-state transitions that occur during MN differentiation in DP and SP, after first examining the efficiency of MN generation in each protocol.
Our single cell data allow us to quantify the dynamics of MN generation in DP and the SP over time, by counting the percentage of cells in each state at each time point. Previous reports have claimed up to 50% MN differentiation efficiency for the SP (Wichterle et al., 2002; Wu et al., 2012), and up to 98% efficiency for DP (after excluding cells that failed to activate the transgene) (Mazzoni et al., 2013). However, both of these studies identified MNs as Mnx1+ cells, and for SP this was measured as GFP+ cells in an Mnx1:GFP reporter system. Since GFP is a stable molecule, and Mnx1 is first expressed in early committed MN progenitors, this definition includes cells from early MN progenitors (between PVNP and MNP in our classification) through to LMNs. If we apply these criteria to our day 5 SP data, we find a similar SP efficiency of 42.2%, and an efficiency of 66% for the DP protocol (not pre-filtering for transgene activating cells).
Using single cell data, it is possible to construct more refined criteria for cell states during MN differentiation, that resolve between the least mature MNPs, maturing EMNs, and mature LMNs, by focusing not on a single marker gene, but on the entire transcriptional state (Figure 2B–E; Figure 2—source datas 1–2). As EMNs mature into LMNs they downregulate progenitor markers including Mnx1, and upregulate markers of terminal differentiation. For DP, MN production was observed as early as day 4 (17.8% EMN), and increased over time to 3.4% EMN and 62.7% LMN by day 11 (Figure 2F). A minority of off-target neuron subtypes, glia, mesoderm, and endoderm cells were also identified, together accounting for <25% of the total population at day 11. We cannot rule out that these fates may emerge from cells that failed to activate the DP transgenes. In the SP, by contrast, we observed significantly lower efficiency in production of the most mature MN states. The SP population contained 18% EMNs and 1.1% LMNs at day 5, and just 6.5% EMNs and 2.6% LMNs at day 12 (9.1% total) (Figure 2G). These figures are lower than those obtained from Mnx1-GFP scoring, because they do not include the earliest committed MN progenitors, which are yet to downregulate progenitor genes such as Mnx1 and Isl1 or to upregulate neuronal genes such as Tubb3 and Map2, and because they distinguish between early and late MNs by their expression of terminal differentiation markers such as Syn1, Mapt, Celf4 or Gabr1. Our single cell classifications therefore capture varying degrees of maturity among differentiating MNs.
An unexpected result from this analysis was that, in the SP protocol, off-target populations increase in fractional abundance to overshadow the intended MN products. Off-target products included oligodendrocytes (7.2%), astrocytes (8.6%), muscle (30.6%), and stroma (8.9%), that together accounted for 55.3% of the population by day 12. The increase in these populations possibly reflects the ability of cell states such as oligodendrocytes and astrocytes to proliferate, which could progressively dilute the initial output of post-mitotic MNs, severely reducing MN generation efficiency. To test this idea and to gain additional EMN and LMN transcriptional data, we profiled an additional 1372 single cell transcriptomes from a replicate SP experiment at day 10 (Figure 3—figure supplement 1; Supplementary methods). At this time point, the SP population contained 34.4% EMNs and 14.9% LMNs, which was higher than at either day 5 or day 12. This is consistent with the hypothesis that SP has an initial output of MNs, which may then be diluted by proliferating off-target lineages.
What are the differentiation paths taken by each protocol? In the SP differentiation path, cells transit through seven states (Figure 2C and E). These state transitions parallel patterning events in the embryo (Wichterle et al., 2002; Davis-Dusenbery et al., 2014; Jessell, 2000): cells first commit to the neural lineage (NP; Sox1+/Sox2+), then are posteriorized (PNP; Sox3+/Hoxb8+/Hoxd4+), ventralized (PVNP; Nkx6-1+/Olig2+), enter the committed MN progenitor state (MNP; Mnx1+), and then mature through early (EMN; Mnx1+/Tubb3+/Map2+) and late (LMN; Mnx1-/Tubb3+/Map2+/Syn1+) MN states. This is not a surprise as the growth factor cocktail defining this method was designed to reflect the signaling events taking place in the embryo. By contrast, we found that the path produced by DP was condensed relative to the SP path (Figure 2B and D), consisting of only four apparent states as opposed to seven. After neural commitment (NP; Sox1+), cells immediately began expressing committed MN markers (EMN; Mnx1+/Tubb3+). This immediate transition was also registered by time course qPCR measurements for a panel of MN genes across bulk populations (Figure 2—figure supplement 3).
A major difference between the DP and SP cell states is the seemingly complete absence of typical spinal embryonic intermediates (PNP or PVNP) in DP. These states are normally recognized by expression of the marker genes Olig2, Nkx6-1, Hoxb8, and Hoxd4, but these genes were not detected by scRNA-Seq in DP intermediates. Olig2 is necessary for MN development in embryos (Takebayashi et al., 2002), indicating that DP drives differentiation from ESC into MNs through a fundamentally new route.
To more confidently demonstrate that DP does not transit through an Olig2+ state, even transiently, we conducted a high-sensitivity gene expression analysis. In SP, our scRNA-Seq data reveals that Olig2 expression peaks (in the PVNP cells) at a level 5.4-fold higher than Gapdh, a constitutively expressed ‘housekeeping’ gene. We performed a dense qPCR time course measuring the expression of Olig2 and a panel of MN genes every day during DP across populations of ~106 cells. These measurements showed that committed MN markers were upregulated immediately following early neural progenitor genes in DP, but Olig2 was not detected at any time point during DP (no qPCR signal at CT > 40), showing that the average expression levels are 106-fold or more lower than Gapdh (Figure 2—figure supplement 3). Accordingly, on any day, the fraction of cells expressing Olig2 at levels seen in SP was less than one in 106. Given that the total protocol duration was 11 days, this would require that an Olig2 transitional state would last less than 1 s (11 days/106), far shorter than the lifetime of mRNA molecules in a cell. Therefore, there is no significant Olig2+ transitional state in DP. This conclusion holds even after allowing for the possibility that Olig2 might act through just a single mRNA transcript. From the qPCR data, the fraction of cells expressing Olig2 at one molecule per cell cannot be more than 0.1% (see Supplementary Methods), and the corresponding lifetime of Olig2+ cells would correspondingly be 15 min (11 days/103), still inconsistent with the lifetime of a single mRNA molecule. These results are consistent with a previous study that reported Olig2 absence during DP (Mazzoni et al., 2013).
Since the DP omits spinal embryonic intermediates characteristic of the SP path, there must be one of two possible trajectories. Either DP must discontinuously transition from an early neural progenitor into a MN, or it must transit through alternate intermediate state(s). To determine which of these possibilities was the case, we employed a data visualization technique called SPRING (Weinreb et al., 2016) to directly compare the topology of both paths. While tSNE is a powerful method for identifying discrete cell states, SPRING provides a complementary description emphasizing continuum gene expression topologies. SPRING builds a k-nearest-neighbor graph over cells in high-dimensional gene expression space, and then renders an interactive 2D visualization of the cell graph using a force directed layout. This representation revealed that the DP and SP trajectories overlap during early neural commitment, but that they then bifurcate and take distinct paths that converge independently to the same MN state (Figure 3A). The dynamics of gene expression over these trajectories resembled the behavior inferred using tSNE, with DP omitting intermediate progenitor genes following its bifurcation from the SP path (Figure 3B).
The bifurcation and subsequent convergence of the two differentiation paths is also revealed by two other complementary analyses. Pairwise cosine similarities between the cell states from both trajectories (Figure 3C; Supplementary methods) indicate similarities between the early states (ESC and NP; cosine similarity >0.69) and late states (LMN; cosine similarity = 0.47), but not the intermediate states (PNP, PVNP, and MNP; cosine similarity <−0.15, –0.03, and 0.04 respectively). We also assigned every individual cell along the DP path to its most similar cluster in the SP path using a maximum likelihood method (Figure 3D; Supplementary methods). This showed that it was virtually impossible to find a single cell resembling the SP intermediate progenitors in the DP approach. Similarity was again observed only at the early and late states.
The bifurcation of the SP and DP trajectories leads to different intermediate cell states. A total of 26 transcription factors (TFs) are differentially expressed between the DP and SP intermediate states (Figure 4A). A majority of these (61%) were involved in an anterior-posterior positional gene expression axis. The SP intermediates were enriched more than 6-fold for nine posterior and spinal TFs including Olig2, Nkx6-1, Lhx3, and six posterior Hox genes with a corrected p-value<0.001. Each of these TFs is expressed in embryonic MNs. By contrast, the DP intermediates were enriched for seven forebrain TFs including Otx2, Otx1, Crx, Six1, Dmrta2, Zic1, and Zic3 at the same stringency, despite the absence of MNs in the forebrain of embryos. Anterior gene expression was previously observed through bulk measurements of DP (Mazzoni et al., 2013), and our results reveal that it occurs within a specific subpopulation of cells in the process of differentiating into MNs.
The positional gene expression signature that characterizes the DP intermediate state appears to be transient. Forebrain gene expression is upregulated along the DP differentiation path, as cells exit the early NP state into the EMN intermediate state (Figure 4B). This transition is accompanied by the downregulation of proliferation-associated genes (Figure 4B; Figure 4—figure supplement 1). By the time cells exit the EMN state and transit into the more differentiated LMN state, they downregulate forebrain genes and replace this abnormal positional signature with a spinal Hox expression signature characteristic of normal MNs (Figure 4B). Thus cells converge to the MN state in positional as well as neuronal identity gene expression in the final stages of DP.
We validated the existence of the abnormal transitional state in DP in two ways using an Mnx1:GFP reporter cell line. Mnx1 serves as a useful reporter in this context since its expression is localized precisely to the distinct intermediate populations of each differentiation path (Figure 3B). First, we performed bulk microarray comparisons of Mnx1+ DP and SP intermediates to compare their gene expression states by an independent method. The comparison confirmed the enrichment of forebrain TFs, and depletion of spinal progenitor and positional genes in the DP intermediates with just one exception – Zic1 was enriched in DP by our single cell comparison but in SP by microarray (Figure 4A). Second, we isolated Mnx1+ DP intermediate cells on day 3 of differentiation, the earliest time point at which Mnx1+ cells appear, and prior to maturation of any Mnx1- GFP+ cells. We cultured this population for an additional 10 days, and validated by immunostaining for Vacht and Map2 that it indeed generates motor neurons (Figure 4C). This experiment confirms that the abnormal transitional state in DP gives rise to MNs. The abnormal DP intermediate thus simultaneously has a forebrain positional expression state and possesses MN fate potential.
Given that the two protocols induce distinct – and in the case of DP, unnatural – differentiation paths, we were curious how their final products compared with primary MNs (pMNs). We harvested MNs from the embryo of a Mnx1:GFP reporter mouse and performed inDrops measurements on 874 Mnx1+ cells that were FACS purified from whole E13.5 spinal cords. Though the majority of Mnx1+ sorted cells were MNs (73.8%, n = 645), this population also contained glia (20.1%), fibroblast-like cells (1.8%), and immune-type cells (1.2%; Figure 5A; Figure 5—figure supplement 1). Using only the cells identified as bona-fide MNs, we compared the differentiating DP and SP cells to pMNs by both global transcriptome similarity of cell states centroids, and a nearest neighbor analysis of single cells. Global transcriptome comparisons confirmed that each state along the DP and SP differentiation paths becomes progressively more similar to pMNs (Figure 5B). The clusters most similar to pMNs were the LMN state from the DP protocol (cosine similarity = 0.60), and the LMN state from the SP (cosine similarity = 0.47). Since subsets of LMNs from DP and the SP might vary in similarity to pMNs, we analyzed the similarity of single cells from all three experiments using SPRING, by embedding all three data sets onto a single kNN graph. We performed this analysis including all cells (Figure 5C–i), and then including only EMNs, LMNs, and pMNs (Figure 5C–ii). Both approaches showed that pMNs closely associate with the LMNs of both DP and SP. It was also apparent that DP and SP LMNs are themselves heterogeneous, with particular subsets associating more closely with pMNs. Overall, a higher fraction of DP LMNs resembled primary MNs, as seen by calculating the fraction of cells in each state that had at least one pMN nearest neighbor out of its 50 most similar cells (Figure 5C–iii; 64% for DP, 6% for SP). DP LMNs therefore appear if anything more related to pMNs in gene expression than SP LMNs, despite their unusual developmental path.
The differences in similarity of DP and SP derived MNs to primary MNs can be understood in terms of differences in maturity and identity of the cells. As expected, both DP and SP EMNs were less mature than pMNs, as seen by expression of the MN progenitor marker Mnx1 (Figure 5C–iv). Among the SP LMNs we detected multiple subpopulations. One subset appeared more anterior than the spinal E13.5 pMNs, as indicated by expression of the hindbrain marker Gata3; the second appeared less mature than pMNs as indicated by residual expression of progenitor genes such as Mnx1 or Nkx2-2; and the final population was closely associated with pMNs (Figure 5C–iv). Most DP LMNs by contrast were closely related to pMNs. They lacked the residual progenitor gene expression seen in SP, and they shared with pMNs the expression of genes important to terminal MN differentiation (e.g. Ebf2 and Ebf3). Beyond individual marker genes, we also systematically compared DP LMNs, SP LMNs, and pMNs by pairwise differential gene expression analysis (Figure 5D). Overall, both DP and SP LMNs show limited gene expression differences to pMNs. In both cases there were 11 differentially expressed TFs. Of these, half were positional Hox genes that mark the most posterior MN pool and were enriched in pMNs, indicating a more anterior identity of the DP and SP LMNs. The remaining TFs included Gata3, which marks a hindbrain subpopulation of SP LMNs as noted above, and other genes related to neuronal maturation including Neurod6, Ebf2, and Mid1. These comparisons confirm the broad similarity of DP LMNs and SP LMNs, to pMNs, but reveal subtle differences between these terminal states that relate to positional identity and state of maturation.
Having established that MNs derived via both gene expression trajectories reach roughly the same MN transcriptional state, we wished to validate that their function and structural organization was also independent of their distinct developmental histories. The SP has been characterized extensively as giving rise to functional MNs (Wichterle et al., 2002), so here we examined structural and functional characteristics following DP. We confirmed that selected protein content matches the mRNA markers by immunostaining for Tubb3, Map2, VACht, Isl1, and Hb9 (Figure 6A). Tubb3 and Map2 were present, and VACht was seen at discrete puncta on the axons (suggesting localization to acetylcholine secretory vesicles). TFs Isl1 and Hb9 were localized in the nucleus. Finally, the GFP from the Mnx1:GFP reporter was activated and expressed in the cytoplasm. To test the functional properties of the DP MNs, we performed whole-cell patch clamp recordings. Depolarization induced single or multiple action potentials in current-clamp experiments (Figure 6B). Depolarizing voltage steps induced fast inward currents and slow outward currents characteristic of sodium and potassium channels, respectively (Figure 6C). Exposure to 500 nM Tetrodotoxin (TTX) blocked the inward current, indicating sodium channel involvement. We then tested whether our DP neurons would respond to neurotransmitters that act on MN (Figure 6D). Exposure of the neurons to AMPA, kainate, GABA, and glycine (100 µM each) induced in each case inward currents similar to that seen in primary embryonic MNs. To see if the DP neurons could also form neuromuscular junctions, we co-cultured the neurons with differentiated C2C12 skeletal muscle myotubes and incubated them for 7 days. We observed clustering of acetylcholine receptors on the C2C12 myotubes near contact points with the DP neurons, which can be seen with alpha-bungarotoxin (α-BGT), which binds to acetylcholine receptors (Figure 6E). We then observed regular contractions of some C2C12 myotubes that began after several days in co-culture (Figure 6E, Figure 6—video 1). These contractions could be stopped by the addition of 300 µM Tubocurarine (curare), an antagonist of acetylcholine receptors, indicating that the contractions were induced by the acetylcholine release from the MNs. Similarly, we noticed that the DP MNs could induce contractions in DP muscle myotubes that we previously generated with MyoD (Figure 6—video 2) (Li and Kirschner, 2014). These results confirm that DP MNs have the expected functional properties of bona-fide MNs.
We provide evidence at the single cell level that differentiation can proceed by multiple routes and yet converge onto similar transcriptional states. While cells differentiated by growth factor stimulation in the SP retrace the embryonic lineage, cells differentiated by terminal state transcription factors in DP take a dramatically different path. The DP path bypasses multiple intermediate progenitor-states that are produced in the embryo, and yet they still converge to the same discrete and recognizable MN phenotype. This convergence occurs via an abnormal intermediate state, and does not appear to involve a shared set of terminal cell state transitions. Moreover, as cells converge they manage to not only establish gene expression related to MN functions, but they also correct positional gene expression defects (exchanging forebrain for spinal gene expression) in the absence of external signals. We conclude that DP of mESCs into MNs occurs via a late bypass that involves alternative intermediate states not seen in the embryo, and that this new route converges near perfectly to the same final state (Figure 1). Convergence into a MN therefore does not depend rigidly on the precise history of intermediate states through which cells differentiate.
This ‘history independence’ of the final state is consistent with a dynamical view of gene regulation in which cell states correspond to ‘attractor basins’, i.e. stable states of gene expression that are robust to modest perturbations. If attractor basins do not exist, the precision of the observed overlap between DP and SP MNs would require a special coincidence. The concept of cell states behaving as attractors has been proposed previously to explain several properties of blood cell types (Chang et al., 2008; Huang et al., 2007; Huang et al., 2005). There are at least two important corollaries of this behavior applied to motor neuron development. From a practical perspective, it is a common concern that DP methods may generate cell types with subtle defects due to their unusual developmental histories (Cohen and Melton, 2011). Attractors would be robust to this vulnerability and indeed our results show that it is not necessary to recreate the precise sequences of steps taken in embryos to generate MNs with the highest similarity to primary cells. It could also hint at a mechanism that might help animal body plans evolve flexibly. Specifically, by decoupling the identity of mature cell state attractors from their developmental histories evolution would be able to act on each independently. In principle this could contribute to evolvability by allowing mature cell states to be transposed onto new lineages in new body locations, without massively rewiring the internal differentiation circuit.
The mechanisms that define the MN attractor basin and allow the artificial DP trajectory to converge onto the correct final state are unclear. The MN state is thought to be stabilized by a network of self-reinforcing TFs (Jessell, 2000), involving Mnr2, Mnx1, Lim3, Isl1, Isl2, and Lhx3, Ngn2, Myt1l, Nefl, and Nefm. DP aims to kick-start this network by activating a subset of important components. Yet, far from immediately activating this network, our data show that DP initially drives cells to differentiate into an early NP state through the same pathway as the SP trajectory, seemingly oblivious to the DP TFs, and then even activates non-MN genes in the transitional state. Understanding why the activation of the MN program lags behind TF induction may provide important clues into how the DP factors act. One possible source for a lag is that activating a complete neuronal program requires first activating additional core TFs (so-called ‘feed-forward’ circuitry). Indeed, recent studies have shown that Ebf and Onecut are activated by Ngn2 (one of the DP TFs), and that both are required to subsequently direct binding of Isl1 and Lhx3 (the other two DP TFs) to MN target genes across the genome during DP (Velasco et al., 2017; Rhee et al., 2016). A second possible source of lag is that extracellular signaling provides inputs that immediately affect cell state, but take time to sensitize cells to the DP factors. For example, signaling changes might activate DP TFs through post-translational modifications, by activating co-factors, or by inducing chromatin state changes. We have indeed observed that MNs are not generated if DP TFs are induced in cells cultured in pluripotency media, indicating a requirement for changes in signaling (not shown). Conversely, when mES cells are transferred to minimal media without inducing DP factors, they acquire a forebrain neural progenitor identity by default (Pankratz et al., 2007; Smukler et al., 2006; Okada et al., 2004). This suggests that the early dynamics and abnormal forebrain/MN expression of the DP transitional state might in fact be driven by the signaling environment and not the DP TFs. These alternatives suggest future experiments to better resolve the mechanisms driving the DP, by re-mapping the trajectory induced during DP after changing signaling conditions, or the choice of DP TFs.
As a methodology, DP is significantly more efficient than the SP without loss of quality in the MN populations produced (Figure 2F–G; Figure 5; Figure 6). The high-efficiency of DP most likely derives from both its more uniform experimental conditions as well as its more direct differentiation path. Experimentally, DP: relies on 2D rather than 3D tissue culture (as in the SP), minimizing uncontrolled cell-cell communication; forces every individual cell to express MN TFs from a genetically integrated construct, increasing uniformity; and employs a defined-media without growth factors that may minimize proliferation of competing progenitor states. The more direct differentiation path induced by DP might also itself increase MN conversion efficiency by minimizing error propagation through off-target fate choices. Long sequences of intermediate cell state transitions offer many sequential opportunities for off-target fate choices, which multiplicatively reduce efficiency, as compared to the more compact DP trajectory. Targeting terminal attractor basins through the shortest possible differentiation paths may prove to be a generally effective strategy for generating desired cell states, more quickly and more enriched in the desired product.
The mouse ES cell line containing doxycycline-inducible Ngn2+Isl1+Lhx3 (NIL) and the Hb9::GFP reporter was graciously provided by Esteban Mazzoni. ESCs were cultured in standard media (DMEM with LIF + 15% fetal bovine serum) on 0.2% gelatin-coated dishes.
Twenty-four hours before starting differentiation, ESCs were trypsinized and seeded onto plates pre-coated with a mix of poly-d-lysine (100 μg/ml) and laminin (50 μg/ml) instead of gelatin for adherence. ESCs were counted by a Beckman Coulter Counter and seeded at a density of approximately 200,000 cells per well of a 6-well dish. At day 0 (ESCs), the media was switched from standard ES media to N2B27 media (Invitrogen). Doxycycline was also added at 3 μg/ml starting to induce expression of the NIL transcription factors. Media was changed daily. For the standard programming protocol, the steps described in Wu et al. were followed strictly. The Wu et al. protocol involves a period of neural induction using embryoid body culture before subsequent replating on laminin. The culture media is serum containing base media supplemented with specific cocktails of growth factors that change during the early, middle, and late stages of differentiation (Figure 2A).
The B6.Cg-Tg(Hlxb9-GFP)1Tmj/J mice (JAX# 005029) were bred with C57BL/6J (JAX# 000664) for embryonic motor neurons dissection. All animal protocols were approved the Institutional Animal Care and Use Committee at Boston Children’s Hospital. On gestational day 13 (E13), the female mice were anesthetized and all embryos were collected by caesarian section. Only GFP embryos were used for further dissection. The spinal cords were isolated and their meninges were removed. Each isolated spinal cord was dissociated by trypsin and mechanical trituration. After filtering the cells with 100 μm strainers, the cells were spun down and re-suspended in PBS, and subjected to flow cytometry. Cells were run through a 100 μm nozzle at low pressure (20 p.s.i.) on a BD FACSaria II machine (Becton Dickinson, USA). A neural density filter (2.0 setting) was used to allow visualization of large cells.
We dissociated differentiating mESC cultures using a 0.25% Trypsin 2 mM EDTA solution (Gibco). Primary HB9+ sorted motor neurons were dissociated as above. Dissociated cell suspensions were verified to be monodisperse and of viability >95% using a coulter counter with trypan blue staining (BioRad). We then performed droplet-based barcoding reverse transcription (RT) reactions and prepared massively multiplexed sequencing libraries using InDrops as described in Klein et al. Briefly: cells, lysis and RT reagents, and barcoding primers attached to hydrogel beads are combined in nanoliter sized droplets suspended in an oil emulsion using the InDrops microfluidic device. A barcoding RT reaction is then carried out in the droplet emulsion, uniquely labeling the RNA contents of every cell using a cell barcode and unique molecular identifier (UMI). Following RT, barcoded emulsions are split into batches, and the emulsion is broken. Combined material is amplified into a nanomolar Illumina sequencing library through a series of bulk reactions: second strand synthesis, in vitro transcription (IVT), fragmentation, RT2, and a final low cycle number PCR. The majority of the amplification takes place during IVT ensuring uniform library coverage. Single cell libraries were then sequenced on either the Illumina HiSeq or NextSeq platforms. Reads were demultiplexed using an updated version of the custom bioinformatics pipeline described in Klein et al.. A python implementation of this pipeline is now publically available on GitHub. Briefly, it filters for reads with the expected barcode structure, splits reads according to their cell barcode, aligns them to a reference transcriptome (we used GRCm38 with some added mitochondrial genome transcripts), and then counts the number of different UMIs appearing for each gene in each cell. The final output is a counts matrix of cells vs. genes that we loaded into MATLAB for further analysis.
Before performing the analyses described below, three steps were taken to ensure that the data was of high quality. First, we required all cells to have at least 1000 UMIs detected. This removed any signal potentially coming from empty droplets. Second, data were total count normalized to ensure differences between cells were not due to technical variation in mRNA capture efficiency or cell size. Finally, cells that had a high stress gene signature were excluded from analysis. Stressed cells were initially recognized as a small percentage of cells (<10%) that clustered apart from everything else and specifically expressed very high levels of a mitochondrial gene set that is associated with cellular stress. Masks that convert raw counts into our filtered set are provided online.
To visualize high-dimensional single cell data, dimensionality reduction is essential. We chose to implement tSNE as described in Klein et al. The core steps are summarized as follows. Steps 1–2 preceed tSNE, and focus the algorithm on genes that best describe differences between cell populations.
Extract the top 1000 highly variable genes. We do this using a statistical test derived specifically for InDrops data (Klein et al.).
Extract principal variable genes. Principal variable genes are a subset of the highly variable genes from step two that we find best describes the cell population structure. The steps to find principal variable genes are:
a.Perform PCA using the top 1000 biologically variable genes.
b.Identify the number of non-trivial principal components. This is done by comparing the eigenvalue of each principal component from a. to the eigenvalue distribution for the same data after being randomized. Only principal components with eigenvalues higher than those observed on random data are retained.
c.Extract genes that contribute most highly to these principal components by imposing a threshold on the gene loadings for each non-trivial PC.
3.Perform tSNE. We used the MATLAB implementation of tSNE from Van de Marteen et al. As input, we supplied z-score normalized principal variable genes, and asked tSNE to perform an initial PCA to a number of dimensions equal to the number of non-trivial principal components in step 4. tSNE then takes cells embedded in this PCA space and nonlinearly projects them onto two dimensions for visualization.
The first goal of our single cell analysis was to describe the identity and proportions of cell states generated by each of two motor neuron differentiation strategies. For this purpose we found that good results were achieved using local-density gradient clustering on the 2D tSNE representation of the data. This approach provided a clear and natural cell state classification that was well aligned with prior knowledge about marker gene expression domains. The steps we took are summarized as follows:
1.Perform tSNE (as above) on cells pooled from all timepoints for each protocol (e.g. Figure 2B).
2.Apply local density gradient clustering to define cell states.
3.Identify genes specifically expressed by each cell state, and use prior knowledge on their expression domains to generate a cell type annotation (e.g. Figure 2D).
4.Quantify the fraction of cells in each state at each timepoint (e.g. Figure 2F).
We identified genes that were specifically expressed by each subpopulation through pairwise t-tests and visually inspected their expression over the tSNE embedding. In Figure 2D–E of the main text, we show the z-score normalized expression of a selection of marker genes that we used as a heatmap. Z-score normalization preserves differences between populations while putting the expression level of every gene on the same scale. In Figure 2—figure supplements 1–2, we show the un-normalized expression of each of these genes individually so that the reader may compare.
During our subpopulation analysis we observed, for both protocols, a continuous progression of cells that spanned the initial ES cell state through the early and late differentiation timepoints. This progression was punctuated by familiar progenitor states that were ordered in a way that was consistent with known events in motor neuron differentiation (Figure 2D–E), and was correlated with chronology (Figure 2B–C inset). We interpreted these progressions as differentiation trajectories. Each is reconstructed from three population snapshots (day 0, day 4/5, and day 11/12). Because differentiation in vitro is asynchronous, the single cell data overlapped from one timepoint to the next. We deduce from this that intermediate states were not missed due to the spacing of our timepoints. We also validated that important intermediate genes were not detected over a densely sampled timecourse using qPCR (see below).
We identified differentially expressed genes between cell states by using two-tailed t-tests with a multiple hypothesis testing correction. We defined differentially expressed genes conservatively at a FDR of 5% and a significance level of p<0.0001. We only considered genes where at least one of the states being compared had >= 10 cells with non-zero expression. To find marker genes of a population we asked for genes that were enriched in that population versus everything else. For several comparisons we restricted our analysis to Riken transcription factors. This list contains ~1500 genes with manually annotated transcription factor activity. We represented differential expression data using volcano plots, and colored the expression intensity of each gene using a colorbar; the mean of the higher expressing state was used.
Due to the low efficiency of EMN/LMN production in our initial SP experiment, we performed a replicate experiment with the goal of extracting additional SP terminal state cells. We repeated SP differentiation and performed InDrops on an additional 1372 cells sample at day 10 of differentiation. To identify EMN or LMN cells within this population we: projected replicate SP cells into the PC-space from the original experiment; calculated a Euclidean distance matrix between all cells; extracted replicate SP cells which had >= 3 of their five nearest neighbors belonging to EMN or LMN clusters; then finally assigned cluster identities to these cells such that cells were called EMN if a majority of their neighbors were EMN, or otherwise were called LMN (majority LMN nearest neighbors or equal EMN and LMN nearest neighbors). This final rule slightly favors LMN over EMN cluster assignment but this helps to counter the low LMN efficiency of the original experiment, which reduces the possible nearest neighbors available during label assignment. For all comparisons of DP and SP we pooled EMN and LMN cells across the two replicate experiments for SP, improving the statistics of the end state comparisons.
After our initial identification of the differentiation paths for the standard protocol and for direct programming we wished compare their routes. One of the most powerful ways to begin addressing such a problem is to simply visualize the data. Yet, we found tSNE gave unclear results for a direct comparison of the paths; visualizing both protocols together resulted in some mixed clusters, and some distinct clusters, but no overall coherence to the representation as we had found looking at one or the other trajectory separately. The limitations of tSNE for analyzing continuous processes are well known.
We therefore turned to a new method developed in parallel to this work in our lab called SPRING (Weinreb et al., 2016), that in our experience does better in analyzing continuous processes in single cell data. SPRING has four core steps: first, it filters for genes with average expression >0.02 UMIs/cell, and Fano factor >2; second, it reduces dimensionality to a 50 dim PCA space; third, it constructs a k-nearest-neighbor (kNN) graph in this space; finally, it renders an interactive 2D visualization of this kNN graph using a force directed layout. In this visualization edges of the kNN graph are literally springs that pull together similar cells, while every cell has a magnetic repulsion force that pushes it away from other cells and an intrinsic gravity. The balance of these forces illuminates the topology of how cells are positioned in high-dimension with respect to one another. In this visualization nodes can be moved around to rotate the projection and find the clearest representation; in each new position the graph re-relaxes according to the underlying physics of the force directed layout. The specific steps we took to generate Figure 3A are as follows. Note that a small manual correction was made to the position of the ESC and LMN populations to reduce white space, but did not distort relationships between populations.
1.Load the cells along each differentiation path into SPRING separately.
Remove doublets. Doublets were identified by three criteria. (1) they are rare, in line with our experimental expectation for cell doublets, (2) they formed long range connections between large cell groups in the same timepoint and sample, (3) they do not possess any unique marker genes; all genes they express are a linear combination of two other cell states. We identified approximately 80 doublets in the standard protocol (<3% of all cells), and approximately 40 in the direct programming approach (<2% of all cells).
Load the filtered cells from both protocols into SPRING together. To make plots the coordinates from the SPRING representation were exported into MATLAB; cells were either colored by cell state (Figure 3A), or gene expression (Figure 3B).
We also asked how the direct programming and standard motor neuron differentiation paths were related to each other by performing a pairwise comparison of cell state centroids. Centroids are the average or center of mass in gene expression space for a collection of cells. Centroids can provide a very accurate estimate of global gene expression that averages over the noise intrinsic to single-cell RNA sequencing at the level of individual cells. By computing the cosine similarity between two cell state centroids we are asking how similar these states are in average gene expression. If the paths do indeed split and then reconverge our expectation was to see similarity between early and late states in both progressions, but not between the intermediate states. We chose to perform these comparisons in a PV-gene space constructed from cells in both protocols to make our comparisons as sensitive as possible. We obtained similar but less sensitive results using all genes above a minimum expression value (not shown). A summary of the steps of this analysis is:
Combine all cells from both protocols, extract PV-genes (as described above), and z-score normalize their expression.
Calculate the centroid of every cluster from each trajectory in this PV-gene space.
Compute the pairwise cosine similarity for all direct programming clusters versus all standard protocol clusters.
Visualize: we chose to use a heatmap (Figure 3C).
Finally, we performed an independent comparison of the differentiation paths that did not depend on our definition of cluster boundaries in the direct programming trajectory. We asked whether any potentially rare individual cells or subpopulations in the direct programming trajectory resemble the intermediate states of the standard trajectory. Because we were now dealing with single cells, not averages of many measurements, care was necessary in this analysis to remain robust relative to the noise in single cell measurements. We therefore took a Bayesian approach and reasoned as follows. In our data, each cell is a vector of counts, with one element for every gene. These counts can be viewed as a multinomial sample from some underlying distribution of gene expression. Since each state of the standard differentiation trajectory is defined by a particular gene expression distribution, we can ask: what is the probability a given direct programming protocol cell was sampled from each standard differentiation trajectory cluster? In this usage, probability amounts to a measure of similarity, with high probability indicating high similarity. We obtained similar results working in either a PV-gene expression space, or considering all genes expressed above a minimum counts threshold. The results of this analysis using PV-genes are presented in Figure 3D. We also obtained similar (but more noisy) results using cosine similarity as an alternative distance metric (not shown). The specific steps of our computation are as follows:
Extract a set of genes with which to make comparisons (either PV-genes or all genes above a minimum expression threshold).
For each cluster of the standard trajectory, calculate the probability of observing a given gene (i.e. the fraction of counts in that cluster from the gene). For genes that are not detected add 1e-07 total counts.
For each cell in the direct programming trajectory, calculate the log-likelihood that it was drawn from each of these clusters. This log-likelihood is from the multinomial distribution function using the probabilities obtained in step 2.
4.Identify and tally the maximum likelihood assignments of all direct programming cells. Normalize raw assignments so that they sum to 100 (giving the percentage). Plot the percentage of direct programming cells assigned to each standard protocol state.
Cell-cycle activity can be estimated from a cell cycle associated gene expression signature; populations that express higher average levels of cell cycle genes are most likely cycling at higher frequency than a population with lower level expression of these genes. In Figure 4b and Figure 4—figure supplement 1 we performed an analysis in this spirit to determine which parts of the DP and SP differentiation trajectories appeared to be proliferative, and to estimate where cells are exiting the cell cycle. We computed a proliferation score that was the aggregate expression of a panel of 21 cell cycle related genes: Aurka, Top2a, Ccna2, Ccnd1, Ccnd2, Ccnd3, Ccne1, Ccne2, Ccnb1, Cdk4, Cdk6, Cdk2, Cdk1, Cdkn2b, Cdkn2a, Cdkn2c, Cdkn2d, Cdkn1a, Cdkn1b, Cdkn1c, Mcm6, Cdc20, Plk1, and Pcna. We also computed a cell cycle exit score on the basis of the aggregate expression of a panel of 4 tumor suppressor genes that inhibit the cell cycle: Cdkn1c, Cdkn1b, Cdkn1a, and Cdkn2d. In Figure 4—figure supplement 1 we show the expression of representative individual genes from this score; in general cell cycle genes were correlated with each other in their expression over cells, as were cell cycle exit genes.
We were interested to convert a bound on the population expression of Olig2, obtained by qPCR, into a bound on the lifetime of a hypothetical rare subpopulation that expressed appreciable levels of this transcript. To proceed, we noted that in SP, Olig2 is expressed 5-fold higher than Gapdh. By contrast, qPCR indicated that in DP, Olig2 was expressed 106-fold lower than Gapdh. Since we took qPCR measurements on every day during differentiation, we next reasoned that the maximum possible number of these intermediate cells would exist if they were spread uniformly across our timecourse (otherwise there would be a spike in their number, increasing the chances of their detection). This is unrealistic, but sets a conservative upper bound. Since the total differentiation protocol took 11 days, any Olig2+ intermediate that expressed the gene at the levels seen in SP (five fold higher than Gapdh) must thus exist for less than 11 days/106/5, or less than 0.2 s. To estimate what this lifetime would be for an Olig2+ population that expressed just one molecule per cell, we conservatively allowed for Gapdh expression levels of 1000 molecules per cell (smRNA-FISH studies have estimated ~200–500 Gapdh copies per cell). The maximum fraction of Olig2+ cells expressing one molecule per cell at any timepoint in DP is correspondingly 0.1% (1000/106*100%). Our lifetime calculation becomes 11 days/103, or less than 15 min. Since the timescales of mRNA production and degradation are typically on the order of hours, we therefore concluded that an Olig2+ subpopulation must not exist during DP.
How does the transcriptional state of motor neurons produced by both protocols compare to that of motor neurons in vivo? To answer this question we leveraged the ability of single-cell RNA sequencing to compare cell states even within populations that are not pure (also see below for functional comparisons of the phenotypes). We performed three analyses. First, we computed the cosine similarity between the centroids of each cell state in both protocols and primary motor neurons (Figure 5B). The specific steps of this analysis were as follows:
1.Combine all cell states from both protocols, and from HB9+ E13.5 primary tissue, extract PV-genes (as described above), and z-score normalize their expression.
2.Calculate the centroid of every cluster from each trajectory, and from primary motor neurons, in this PV-gene space.
Compute the cosine similarity for all in vitro populations versus primary motor neurons, and visualize as a bar graph.
Second, we compared single cells from the DP and SP terminal states to pMN reference single cells. The goal of this analysis was to detect whether heterogeneity within the DP and SP terminal states correlated with similarity to pMNs. To gain a broad understanding of the relationship between the DP and SP trajectories to our pMN reference, we initially projected pMN single cells into the PCA-space used in the visualization from Figure 3A (Figure 5C–i). The projection was performed by multiplying the pMN counts matrix (after quality filtering cells and z-score normalizing genes) by the PC-loadings matrix that was built on the DP and SP cells as described above. The combined coordinates of all cells in this space were then used as input to SPRING, which builds and visualizes a kNN graph. To generate a more refined visualization of how the terminal populations were related to each other we next constructed a new PC-space on EMN and LMN DP and SP cells, combined with pMNs, following the same steps of SPRING visualization described above (Figure 5C–ii). This differs from the original projection in that the second PC-space was constructed to reflect variation within only the mature populations of cells (EMN, LMN and pMN). It thus has higher sensitivity to detect subtle differences between MN subpopulations, because the PC-space does not contain dimensions that describe unrelated process that occur in the process of becoming a neuron from an mES transcriptional state, and thus only introduce noise into the comparison of terminal states. To quantify how similar individual cells were to the pMN reference we performed nearest neighbor analysis on the distance matrix underlying this second SPRING visualization. For each DP or SP terminal state, we asked of every cell whether it had a pMN cell among its 50 most similar cells, and found the fraction of cells with similarity to pMNs. The number of nearest neighbors chosen, in this case 50, sets the size of the region in which we required a pMN cell to fall in order for a cell to be classified as similar to pMN. 50 cells equated to roughly 2% of the total cells in this visualization, and so provided an appropriate and stringent threshold.
Third, we performed differential gene expression analysis, comparing the most mature motor neuron states from each protocol with primary HB9+ motor neurons (Figure 5D). This analysis was performed as described above.
RNA was isolated using Qiagen Rneasy Plus Kit. Purified RNA was then reverse transcribed using Bio-rad's iSCRIPT cDNA synthesis kit. Quantitative PCR was performed using Bio-rad's SYBR green Supermix on a CFX96 Real-Time PCR system.
Antibodies used for immunostaining were anti-Tubb3 (Cell Signaling D71G9; 1:100), anti-Map2 (Sigma M9942; 1:200), anti-VACht (Synaptic Systems 139 103; 1:200), anti-Isl1 (Abcam ab109517; 1:1000), and anti-Hb9 (DSHB 81.5C10; 1:100). Differentiated cells were fixed in 4% paraformaldehyde for 20 min and then permeabilized with 0.1% Triton-X for 15 min. After primary incubation for 1 hr, samples were labeled with a secondary antibody conjugated to AlexaFluor647. Samples were co-stained with DAPI before imaging on a Nikon Eclipse TE2000-E microscope.
All recordings were carried out at room temperature within 6 days of plating the neurons in 35 mm dish. Whole-cell voltage clamp recordings were made with a Multiclamp 700B amplifier (Molecular Devices) and patch pipettes with resistances of 2–3 MΩ. Pipette solution was 135 mM K-Gluconate, 10 mM KCl, 1 mM MgCl2, 5 mM EGTA, 10 mM HEPES, pH 7.2, adjusted with NaOH. The external solution was 140 mM NaCl, 5 mM KCl, 2 mM CaCl2, 2 mM MgCl2, 10 mM HEPES, and 10 mM D-glucose, pH 7.4, adjusted with NaOH. We used gravity perfusion system connected with Perfusion Pencil with Multi-Barrel Manifold Tip (AutoMate Scientific) to externally apply 0.5 µM tetrodotoxin, 100 µM AMPA, kainate, GABA, or glycine to the cells. Command protocols were generated and data was digitized with a Digidata 1440A A/D interface with pCLAMP10 software.
C2C12 myoblasts were grown in 10%FBS + DMEM media and then differentiated into myotubes by incubating in differentiation medium (2% horse serum + DMEM). After myotubes were formed, the neurons were dissociated by trypsinization and reseeded on top of the differentiated muscle to allow contractions to develop. Video of contractions were taken using Metamorph software and manually counted over 5 s intervals. For stopping assay, 300 µM Tubocurarine (Sigma) was added to the media as an acetylcholine competitor. For labeling of acetylcholine receptors, bungarotoxin (Invitrogen) was used after cells were fixed with paraformaldehyde. Similarly to C2C12 cells, ES cells over-expressing MyoD were also differentiated to myotubes using differentiation medium and subjected to co-culture with neurons.
The mRNA of undifferentiated iNIL ES cells and NIH 3T3 cells (grow in DMEM + 10% FBS) were collected and purified by RNA extraction using RNeasy Plus Extraction Kit (Qiagen). Neurons differentiated by either protocol were first sorted by flow cytometry on a BD FACSaria II machine (Beckton Dickinson, USA) to collect Hb9::GFP+ cells, and then were subjected to RNA extraction in a similar fashion. Collected RNA was then amplified and hybridized to Affymetrix GeneChip Mouse Transcriptome Arrays (MTA 1.0). Results were processed by the Children’s Hospital Microarray Core Facility, and were analyzed using Affrymetrix’s Transcriptome Analysis Console and Expression Console software.
Turning straw into gold: directing cell fate for regenerative medicineNature Reviews Genetics 12:243–252.https://doi.org/10.1038/nrg2938
Cell fates as high-dimensional attractor states of a complex gene regulatory networkPhysical Review Letters 94:128701.https://doi.org/10.1103/PhysRevLett.94.128701
Bifurcation dynamics in lineage-commitment in bipotent progenitor cellsDevelopmental Biology 305:695–713.https://doi.org/10.1016/j.ydbio.2007.02.036
Neuronal specification in the spinal cord: inductive signals and transcriptional codesNature Reviews Genetics 1:20–29.https://doi.org/10.1038/35049541
Modeling ALS with motor neurons derived from human induced pluripotent stem cellsNature Neuroscience 19:542–553.https://doi.org/10.1038/nn.4273
Wishbone identifies bifurcating developmental trajectories from single-cell dataNature Biotechnology 34:637–645.https://doi.org/10.1038/nbt.3569
Embryonic stem cells assume a primitive neural stem cell fate in the absence of extrinsic influencesThe Journal of Cell Biology 172:79–90.https://doi.org/10.1083/jcb.200508085
Visualizing Data using t-SNEJournal of Machine Learning Research : JMLR 9:2579–2605.
Efficient differentiation of mouse embryonic stem cells into motor neuronsJournal of Visualized Experiments e3813.https://doi.org/10.3791/3813
Extreme makeover: converting one cell into anotherCell Stem Cell 3:382–388.https://doi.org/10.1016/j.stem.2008.09.015
Martin PeraReviewing Editor; University of Melbourne, Australia
In the interests of transparency, eLife includes the editorial decision letter and accompanying author responses. A lightly edited version of the letter sent to the authors after peer review is shown, indicating the most substantive concerns; minor comments are not usually included.
Thank you for submitting your article "Mouse embryonic stem cells can differentiate via multiple paths to the same state" for consideration by eLife. Your article has been reviewed by three peer reviewers, one of whom is a member of our Board of Reviewing Editors and the evaluation has been overseen by Marianne Bronner as the Senior Editor. The reviewers have opted to remain anonymous.
The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission. Essential revisions are listed below and we have also included the detailed reviewer comments for your information.
This manuscript uses single cell transcriptomics to compare the cell state transitions that occur when mouse embryonic stem cells are directly reprogrammed into motor neurons using defined transcription factors (direct protocol) to those that occur during a directed differentiation designed to simulate normal embryological induction processes (standard protocol). The main conclusion is that the two pathways converge on a similar endpoint via distinct transition states. The authors assess this endpoint through functional assays and through comparison with bona fide neurons isolated from animals.
All three reviewers are enthusiastic about the general concept and approach of your study. However, two Reviewers note that the induction of forebrain markers in the direct protocol might be an artifact of culture conditions rather than an obligate intermediate in the pathway. The reviewers suggest that 1) knockdown of the patterning factors in the direct and standard protocol or 2) isolation of intermediates to follow their fate prospectively or 3) early induction of factors by dox would strengthen the conclusions of the study and provide more mechanistic insight. Two reviewers also note that the efficiency of induction of motor neurons in the standard protocol is very low. Repetition of some of the key experiments with a higher yield would overcome these concerns. Finally, the reviewers request more clarity around the definition of particular cell types. Specific points are listed below.
1) Though the overall concept behind this study is interesting, there are some significant issues with the experimental design. First, it is not clear to what extent the authors' observations are a function of the differences in the particular culture conditions employed-in particular transfer of pluripotent cells to serum-free medium in adherent culture in the direct protocol (which can induce forebrain specification as a default), versus three-dimensional culture or EB formation in the standard protocol (which can markedly influence heterogeneity and outcomes). Perhaps a greater concern with this study is the very low yields of the desired cell type from the standard protocol. The motor neuron yields peak at Day 5 with the standard protocol (about 25% early motor neurons) then a variety of other cell types are generated such that by Day 12, there are very few cells with properties of motor neurons. By contrast the direct protocol is highly efficient (62.7% versus 2.7% motor neurons). It is questionable whether any meaningful comparison can be drawn between these two protocols; when one yields only 3% of the desired endpoint using a different culture platform, it is hardly surprising to observe different intermediates. The protocol the authors follow for standard differentiation (Wu et al., 2012) reports ~50% motor neuron yield at 7 days; much closer to the direct protocol yields the authors report; it is not clear why they did not achieve a better yield in this study from the protocol of Wu et al., 2012.
2) The percentage of MNs in the SP protocol is surprisingly low. The culture conditions, crucial for this conclusion, should be clearly described in methods. Recapitulating development, SD protocol produces oligodendrocytes after MNs, thus efficiency as percentage of total cells it is expected to decrease.
3) MN are patterned along the rostro-caudal (R-C) axis of the spinal cord to innervate the different target by Hox genes. Thus, MN subtype identity varies greatly along the R-C axis. It is expected that neither of the in vitro differentiation protocol decapitates the diversity of embryonic derived MNs. Specifically, the SP was designed to make a very restricted set of MNs within the R-C axis. Without description of the method, it is left to the reader to assume MNs were collected from all R-C levels.
4) Single cell RNA seq data is biased towards detecting highly expressed genes. While the DP protocol appears to follow a very artificial trajectory it is possibly that key genes that are crucial for generation of the SD intermediates that the DP approach seems to skip are transiently expressed at very low levels. The authors could/should show that shRNA mediated knockdown of the key TF that are important for in vivo/SD generation of MNs but appear to get skipped in the DP protocol, reduce differentiation outcomes for the SD but not the DP protocol. This will be crucial to support that DP MNs are indeed derived by a process that skips/bypasses key developmental intermediates.
5) It will be important to investigate if the generation of Sox1/Sox3+ neural progenitors during the DP protocol is not simply an artefact of the authors experimental conditions. The authors cited (3) and stated in the introduction that MNs are generated by the DP approach by skipping the a Sox1 positive neural progenitor state. As neural commitment in N2B27 is a very rapid process and mESCs are exposed to the dox and N2B27 media at the same time it is possible that the cells start to differentiate before the NIL factors become active on a protein level. While the authors state that MNs cannot be generated by dox addition if cells are not removed from mESC maintenance media, pulsing the cells for 12, 24, or 48 hours with dox in mESC media before on set of differentiation very well might (the original DR publication shows that NIL factors are present on protein level by 48 hours). I understand that asking the authors to repeat the single cell transcriptome under this condition will be excessive, however I propose that the authors use qPCR (single cell) for neural progenitor markers and mature MN markers to test whether pre-exposure to dox abolishes the generation of neural progenitors during N2B27 culture. If it does this would indicate that SD/DR approaches are entirely different processes that do not even share the neural progenitor state and more importantly that mESCs do directly transition to MN-like state.
6) As commented above, it would be interesting to test whether knock-down of forebrain transcription factors that get transiently upregulated during DP mediated MN generation reduce differentiation outcomes for the DP protocol. This would help to demonstrate more clearly and in a mechanistic way, whether the forebrain signature is important for the process or merely an artefact imposed by the N2B27 media which generates under normal differentiation conditions forebrain fates.
7) Finally, to definitely proof that both the SD and DP protocol give rise to MN via different intermediate stages, it will be critical to isolate the presumed intermediate cells by FACS (either by using specific cell surface markers inferred from the single cell data or at least they should use their MNX1:GFP line) and show that these "presumed" intermediates are indeed the source of the mature MNs. That is, the authors should re-culture the isolated cells and demonstrate by IF characterisation that they give rise to the final end product.
This manuscript uses single cell transcriptomics to compare the cell state transitions that occur when mouse embryonic stem cells are directly reprogrammed into motor neurons using defined transcription factors (direct protocol) to those that occur during a directed differentiation designed to simulate normal embryological induction processes (standard protocol). The main conclusion is that the two pathways converge on a similar endpoint via distinct transition states. The authors assess this endpoint through functional assays and through comparison with bona fide neurons isolated from animals.
The manuscript addresses an important question in cellular differentiation and stem cell biology. It is fair to point out that the direct reprogramming of embryonic stem cells (which are pluripotent) is liable to be significantly different in important ways to direct reprogramming of differentiated adult somatic cells (the more likely starting point for many direct reprogramming protocols). The latter type of conversion most often requires crossing or upwards traverse of cell lineages whereas the former might be expected to be more closely related to normal developmental pathways. All the same this study is certainly a valid model for direct reprogramming, and having a common starting point for the standard and direct protocols simplifies the analysis.
Though the overall concept behind this study is interesting, there are some significant issues with the experimental design. First, it is not clear to what extent the authors' observations are a function of the differences in the particular culture conditions employed-in particular transfer of pluripotent cells to serum-free medium in adherent culture in the direct protocol (which can induce forebrain specification as a default), versus three-dimensional culture or EB formation in the standard protocol (which can markedly influence heterogeneity and outcomes). Perhaps a greater concern with this study is the very low yields of the desired cell type from the standard protocol. The motor neuron yields peak at Day 5 with the standard protocol (about 25% early motor neurons) then a variety of other cell types are generated such that by Day 12, there are very few cells with properties of motor neurons. By contrast the direct protocol is highly efficient (62.7% versus 2.7% motor neurons). It is questionable whether any meaningful comparison can be drawn between these two protocols; when one yields only 3% of the desired endpoint using a different culture platform, it is hardly surprising to observe different intermediates. The protocol the authors follow for standard differentiation (Wu et al., 2012) reports ~50% motor neuron yield at 7 days; much closer to the direct protocol yields the authors report; it is not clear why they did not achieve a better yield in this study from the protocol of Wu et al., 2012.
1) Results section–sometimes embryos do reach a common terminal differentiation endpoint through different pathways. Cartilage (derived from neural crest or mesoderm) is one example and there are others.
2) Results section–be more explicit about how the various cell states are defined.
3) Results section–how exactly were these efficiencies calculated? Do they take into account cell proliferation and loss at various stages?
4) Results section–is the intermediate DP state a necessary stage on the journey? Maybe it represents a default (early anterior stages) induced by growth factor depletion that then gets overtaken by the reprogramming factors.
5) Discussion section– because the MN endpoint is not characterised at a functional level by the authors in this study, the conclusions rely on comparison to historical controls, and therefore depend on the degree to which the authors' protocol mimics earlier efforts. Direct comparison of the two endpoints is only possible at the transcriptional level.
Briggs, Li and coworkers aim to investigate if the developmental path taken by differentiating cells determines the final cellular state or if two alternative trajectories are able to produce a similar terminal fate. To do so, the collect single cell RNA-seq data from differentiating motor neurons induced by extracellular signals or by forced expression of transcription factors. An important difference between these samples rests on the fact that direct differentiation by forced transcription factor expression bypasses the canonical motor neuron progenitor stages. The results show that direct programming and standard differentiating cells diverge at early progenitor like state to take different path and converge on a final similar cell state. Although this was postulated before, these results test the hypothesis with the appropriate and up to date technology.
There are some important points that in my opinion should be addressed;
1) A better description of the criteria used to assign cell types will be useful. For example, motor neuron progenitors (MNP) are characterized by Olig2 expression, low or no Isl1, Mnx1 and Lhx3. In Figure 1E, MNP have an average Olig2 and high Mnx1, Isl1 and Lhx3. Thus, some of the cells in this cluster might be early MNs than progenitors. Are the markers in Figure 1D-E some of the genes driving the clustering?
2) The percentage of MNs in the SP protocol is surprisingly low. The culture conditions, crucial for this conclusion, should be clearly described in methods. Recapitulating development, SD protocol produces oligodendrocytes after MNs, thus efficiency as percentage of total cells it is expected to decrease. Where there antimitotic added to both (DP and SP) conditions to restrict the generation of mitotic byproducts?
3) It is my opinion that Figure 5 should be divided into two figures. The first two panels are intrinsically linked to the goal of this manuscript, the characterization of DP MNs is less relevant and novel. The three-way comparison of SP, DP and PM is extremely valuable and on topic. Thus, I believe I deserves to be expanded and explained better. The supplementary figures associate with this figure should be in the main figure. On that note, it will be useful to color Figure 5—figure supplement 2A based on DP, SP and pMN to understand where cells "land" in these dimensions. Describe what is gene set that take DP closer to pMN than SP for example? Is it possible to add pMN and remake here Figure 3A?
4) MN are pattern along the rostro-caudal (R-C) axis of the spinal cord to innervate the different target by Hox genes. Thus, MN subtype identity varies greatly along the R-C axis. It is expected that neither of the in vitro differentiation protocol decapitates the diversity of embryonic derived MNs. Specifically, the SP was designed to make a very restricted set of MNs within the R-C axis. Without description of the method, it is left to the reader to assume MNs were collected from all R-C levels. Thus, it this is not a surprised or anything that should be left as a speculation and stated that way.
5) Although with some novel detail, the rest of Figure 5 is a recapitulation of previous work. Finally, the Tubb3+BTX staining is not informative. Also, clustering of acetylcholine receptors when cells are platted without neurons is required to conclude that MNs induce clustering. A Vacht or SV2 staining plus BTX and quantification of the overlap will be a better measure of active NMJs.
6) The finding that in the DP protocol cells don't transition through MNP stages is not novel. However, I do understand that for the purpose of comparing DP to SP it is a necessary description. Single cell RNA-seq and differentiation path is the appropriate and cutting edge tool to investigate the original question: does developmental path matters? The cellular system is appropriate too since three types of MNs can be compared: DP, SP, and pMN. The manuscript is clearly written following a logical argument. I believe it suffers from a very poor description of the experiment and methods that hurst the interpretation of the data. Also, it is superficial in important description as those single cell RNA-seq in Figure 5.
The study, "Mouse embryonic stem cells can differentiate via multiple paths to the same state" by Briggs et al. uses single cell sequencing to demonstrate that mouse ESCs can give rise to motor neurons via alternative transcriptional trajectories. In detail the authors compare a traditional cytokine driven differentiation protocol (SD) with a transcription factor based direct "programming" approach. While the authors were able to detect discreet developmental intermediates in the SD protocol that correspond to known in vivo counter parts, the DP approach followed a mostly alternative path. While both DP and SD cultures appeared to transition into neural progenitors initially, beyond that they seemed to directly transition into an early motor neuron state bypassing a number of intermediate stages. The authors postulate that the same cell fate can be derived from mESCs via alternative routes that diverge after the neural commitment stage.
This is an important and highly interesting study that tries to unveil how direct programming produces its target cell type and whether this is achieved by an accelerated form of differentiation (transcription factor guided) or whether this occurs via an alternative route/more direct cell state transition. While the authors use state of the art bioinformatical analyses to mine their data set is well done, and allowed them reach interesting results, I believe that a few key experiments are required for the authors to mechanistically support their key findings.
– Single cell RNA seq data is biased towards detecting highly expressed genes. While the DP protocol appears to follow a very artificial trajectory it is possibly that key genes that are crucial for generation of the SD intermediates that the DP approach seems to skip are transiently expressed at very low levels. The authors could/should show that shRNA mediated knockdown of the key TF that are important for in vivo/SD generation of MNs but appear to get skipped in the DP protocol, reduce differentiation outcomes for the SD but not the DP protocol. This will be crucial to support that DP MNs are indeed derived by a process that skips/bypasses key developmental intermediates.
– It will be important to investigate if the generation of Sox1/Sox3+ neural progenitors during the DP protocol is not simply an artefact of the authors experimental conditions. The authors cited (3) and stated in the introduction that MNs are generated by the DP approach by skipping the a Sox1 positive neural progenitor state. As neural commitment in N2B27 is a very rapid process and mESCs are exposed to the dox and N2B27 media at the same time it is possible that the cells start to differentiate before the NIL factors become active on a protein level. While the authors state that MNs can not be generated by dox addition if cells are not removed from mESC maintenance media, pulsing the cells for 12, 24, or 48 hours with dox in mESC media before on set of differentiation very well might (the original DR publication shows that NIL factors are present on protein level by 48 hours). I understand that asking the authors to repeat the single cell transcriptome under this condition will be excessive, however I propose that the authors use qPCR (single cell) for neural progenitor markers and mature MN markers to test whether pre-exposure to dox abolishes the generation of neural progenitors during N2B27 culture. If it does this would indicate that SD/DR approaches are entirely different processes that do not even share the neural progenitor state and more importantly that mESCs do directly transition to MN-like state.
– As commented above, it would be interesting to test whether knock-down of forebrain transcription factors that get transiently upregulated during DP mediated MN generation reduce differentiation outcomes for the DP protocol. This would help to demonstrate more clearly and in a mechanistic way, whether the forebrain signature is important for the process or merely an artefact imposed by the N2B27 media which generates under normal differentiation conditions forebrain fates.
– Finally, to definitely proof that both the SD and DP protocol give rise to MN via different intermediate stages, it will be critical to isolate the presumed intermediate cells by FACS (either by using specific cell surface markers inferred from the single cell data or at least they should use their MNX1:GFP line) and show that these "presumed" intermediates are indeed the source of the mature MNs. That is, the authors should re-culture the isolated cells and demonstrate by IF characterisation that they give rise to the final end product.
An experiment that I have not suggested because it maybe an excess, but it would really enhance this manuscript, would be to use lineage tracing for Olig2 or Nkx6-1 during the DP protocol, and in this way show that they are never upregulated during the generation of MN by the DP approach.https://doi.org/10.7554/eLife.26945.027
- Victor C Li
- Marc W Kirschner
- Allon Klein
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
We are grateful to Esteban Mazzoni for providing the transcription factor cassette containing ES cell line used in the DP experiments presented here.
Animal experimentation: This study was performed in strict accordance with the recommendations in the Guide for the Care and Use of Laboratory Animals of the National Institutes of Health. All of the animals were handled according to approved institutional animal care and use committee (IACUC) protocols (#16-01-3080R) of Boston Children's Hospital, and (#IS00000137, #1648, #1648a, #1648b) of Harvard Medical School. The Hb9-GFP embryos were collected from a pregnant mouse after euthanasia by inhalation of CO2, and every effort was made to minimize suffering.
- Martin Pera, Reviewing Editor, University of Melbourne, Australia
© 2017, Briggs et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.