Citation: Sagner A, Gaber ZB, Delile J, Kong JH, Rousso DL, Pearson CA, et al. (2018) Olig2 and Hes regulatory dynamics during motor neuron differentiation revealed by single cell transcriptomics. PLoS Biol 16(2):
Academic Editor: Marianne Bronner, California Institute of Technology, United States of America
Received: May 25, 2017; Accepted: January 5, 2018; Published: February 1, 2018
Copyright: © 2018 Sagner et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Data Availability: Single cell RNA sequencing data are available via ArrayExpress (http://www.ebi.ac.uk/arrayexpress/experiments/E-MTAB-5466). All other relevant data has been uploaded as Supporting information files.
Funding: Wellcome Trust (grant number FC001051; WT098326MA). The funder had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. EMBO Long-term fellowship (grant number 1438-2013). The funder had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. Cancer Research UK (grant number FC001051). The funder had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. Whitehall Foundation (grant number 2004-05-90-APL). The funder had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. March of Dimes Foundation (grant number 5-FY06-7). The funder had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. NINDS (grant number R01NS053976, R01NS072804, R01NS085227). The funder had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. UK Medical Research Council (grant number FC001051). The funder had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. Human Frontiers Science Program (grant number LT000401/2014-L). The funder had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. People Programme (Marie Curie Actions) of the European Union’s Seventh Framework Programme (grant number FP7-2013 624973). The funder had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Competing interests: The authors have declared that no competing interests exist.
basic helix-loop-helix; Chat,
choline acetyltransferase; Chip-Seq,
chromatin immunoprecipitation-sequence; CHIR,
Wnt pathway activator CHIR99021; Cre,
bacteriophage P1 Cre recombinase; DBZ,
Delta-like 1; e,
embryonic day; E-box,
bHLH transcription factor binding site; EGFP,
enhanced green fluorescent protein; EnR,
Engrailed transcriptional repression domain; EP,
embryonic stem cell; Fabp7,
fatty-acid binding protein 7; FGF,
fibroblast growth factor 2; FLAG,
FLAG epitope tag; FP,
floor plate; GFP,
green fluorescent protein; GSK3,
glycogen-synthase kinase 3; GO,
gene ontology; Hes,
Hairy/Enhancer of Split; HH,
Hes5(e1) genomic element; IN,
internal ribosomal entry site; Isl1,
K shortest paths; mKate2,
monomeric far-red fluorescent protein Katushka-2; MN,
motor neuron; nEGFP,
nuclear EGFP; Ngn2,
Neurogenin 2; NMP,
neuromesodermal progenitor; NLS,
nuclear localization signal; NMP,
neuromesodermal progenitor; NP,
neural progenitor; ns,
not significant; N2,
Ngn2 protein; N2B27,
N2 and B27 media supplements; Olig2,
oligodendrocyte transcription factor 2; O2,
Olig2 protein; PKA,
protein kinase A; pMN,
MN progenitor; p2,
V2 interneuron progenitor; p3,
V3 interneuron progenitor; RA,
retinoic acid; RT-qPCR,
real-time quantitative polymerase chain reaction; SAG,
Smoothened/Shh signaling agonist; scRNA-seq,
single cell RNA sequencing; sgRNA,
short guide RNA; Shh,
sonic hedgehog; TF,
transcription factor; Tubb3,
neuronal class III beta-tubulin; T2A,
Thosea asigna virus 2A peptide sequence; V2,
V2 interneuron; V3,
V3 interneuron; WT,
The orderly development of embryonic tissues relies on gene regulatory networks that control patterns of gene expression, tissue growth, and cell differentiation [1,2]. Genetic and molecular studies have identified many of the constituents of these networks and have begun to define the regulatory hierarchy between them. Nevertheless, how cell fate assignment is coordinated with proliferation and differentiation remains poorly understood.
An experimentally well-characterized tissue that exemplifies this problem is the vertebrate spinal cord. In ventral regions of the developing spinal cord, proliferating progenitors are exposed to a gradient of sonic hedgehog (Shh) signalling that controls the expression of a set of homeodomain and basic helix-loop-helix (bHLH) transcription factors (TFs) [3–5]. These TFs form a gene regulatory network that progressively allocates progenitor identity, dividing the spinal cord into molecularly discrete domains arrayed along the dorsal-ventral axis [6,7]. This combinatorial transcriptional code determines the subtype identity of the postmitotic neurons generated by progenitors in each domain, thereby controlling the position at which motor neurons (MNs) and interneurons emerge [3,8–10].
Among the first neurons to differentiate in the ventral spinal cord are MNs. In mouse and chick, these are formed over a 2–3-day period . During this time, most if not all MN progenitors exit the cell cycle and differentiate, whereas the adjacent progenitor domains that give rise to interneurons continue to divide and, consequently, differentiate at a much slower pace [11,12]. These differences in differentiation rate play an important role in the elaboration of spinal cord pattern and ensure that appropriate numbers of MNs are generated. This raises the question of how the regulatory mechanisms defining MN progenitors prime these cells to differentiate rapidly.
The induction and differentiation of MNs are characterized by a series of gene expression changes. Initially, Shh signaling induces the bHLH protein oligodendrocyte transcription factor 2 (Olig2), resulting in the repression of the homeodomain protein Irx3 and bHLH protein Bhlhb5 (also known as Bhlhe22) normally expressed in neural progenitors (NPs) dorsal to MNs [13–16]. Ectopic expression of Olig2 represses both Irx3 and Bhlhb5, resulting in ectopic MN production [13,16,17]. Conversely, in the absence of Olig2, MN generation fails, and instead, Irx3 and Bhlhb5 expression are maintained, and NPs differentiate into ventral interneurons [14–16,18].
The gene regulatory mechanisms that are responsible for the higher rate of neurogenesis of MN progenitors compared to other NPs in the spinal cord are not well understood. Whether Olig2 functions as an activator or inhibitor of neurogenesis is also unclear. Initial studies indicated that expression of Olig2 accelerates cell cycle exit , and the absence of Olig2 results in a characteristically slower tempo of neuronal differentiation . Olig2 promotes the expression of the proneural bHLH TF Neurogenin 2 (Ngn2), and ectopic expression of Ngn2 causes progenitor cells to exit the cell cycle and prematurely differentiate into neurons [13,17,19–23]. These studies also showed that Olig2 acts as a transcriptional repressor to promote Ngn2 expression [13,17], implying that Olig2 elevates Ngn2 expression by negatively regulating the expression of Ngn2 repressors. Candidate Ngn2 repressors include members of the Hairy/Enhancer of Split (Hes) family of TFs, which act downstream of the Notch signaling pathway to prevent neuronal differentiation and maintain progenitors in a dividing, undifferentiated state [24–26].
Although these studies suggested that Olig2 promotes motor neurogenesis, subsequent studies ascribed antineurogenic and pro-proliferative functions to Olig2 (reviewed in ). These conclusions were based on the Olig2-mediated repression of Mnx1 (Hb9), a homeodomain transcription factor expressed by postmitotic MNs , and the cell cycle inhibitor p21 . Olig2 also has the capacity to form heterodimers with Ngn2 that inhibit Ngn2’s neurogenic activity , and oppose the functions of the tumor suppressor protein p53 . Furthermore, addition of Olig2 to TF reprogramming cocktails inhibits conversion of fibroblasts to MNs , supporting the idea that Olig2 interferes with the differentiation of MNs. Thus, although the genetic evidence establishes Olig2 as a key determinant of MN identity, these apparently contradictory findings leave unexplained how Olig2 coordinates specification of neuronal identity while determining the rate of differentiation.
Single cell RNA sequencing (scRNA-seq) is emerging as a novel and powerful technology to identify distinct cell types in complex mixtures and to define developmental trajectories during differentiation [32–36]. Here, we took advantage of an in vitro model that allows the generation of ventral spinal cord cell types from embryonic stem cells (ESCs) to perform scRNA-seq analysis of developing NPs . We used these data to reconstruct and validate the differentiation trajectory of MN progenitors and to infer the gene regulatory mechanisms by which Olig2 promotes MN differentiation. Both in vivo and in vitro cells commit to MN differentiation asynchronously. This limits the temporal resolution of conventional gene expression assays, potentially obscuring details of the sequence of events during MN differentiation. Here, we developed a method to reconstruct the differentiation trajectory from scRNA-seq data that provides much greater temporal resolution of the transcriptional dynamics during MN differentiation than previously available. This approach identified a sequence of distinct phases in MN differentiation, including two distinct Olig2 expression states: an initial Olig2LOW state, during which Hes1 expression decreases and Olig2 is coexpressed with Hes5, and a subsequent Olig2HIGH state, in which high levels of Olig2 promote differentiation by repressing Hes5, thereby indirectly inducing Ngn2. We validated this two-phase model using quantitative image analysis of a fluorescent Olig2 reporter and provide in vitro and in vivo evidence that Olig2 acts directly on Hes genes to promote cell cycle exit and neurogenesis in the MN progenitor domain (pMN domain). Together, the data provide a comprehensive view of the regulatory network that controls the specification of MN progenitors and identify a molecular mechanism coordinating the specification of positional identity with differentiation.
In vitro generation of MN and V3 interneuron progenitors
To define the sequence of events that leads to the generation of somatic MNs, we took advantage of ESCs, which can be directed to differentiate into spinal NPs in vitro . This method relies on the exposure of ESCs, cultured as a monolayer, to a brief pulse of Wnt signalling prior to neural induction (Fig 1A). This induces the caudalizing TFs Cdx1,2,4 . Subsequently, removal of Wnt signalling and concomitant exposure to retinoic acid (RA) and the Smoothened/Shh signalling agonist (SAG) results in the generation of NPs expressing progenitor markers characteristic for the ventral spinal cord, such as Olig2 and Nkx2.2 (Fig 1B), and MNs expressing postmitotic markers, including Islet1 (Isl1), Mnx1, and neuronal class III beta-tubulin (Tubb3) (Fig 1B and 1C and S1A Fig). These NPs initially express Hoxb1 and, later, Hoxb9 (S1B Fig) and differentiate into Hoxc6-positive MNs, characteristic of forelimb level spinal cord MNs (S1B and S1C Fig) [37–40].
Fig 1. Characterization of MN differentiation from ESCs.
(A) Scheme outlining the differentiation protocol. ESCs are plated in N2B27 + FGF for 2 days before being exposed to N2B27 + FGF/CHIR, resulting in the production of NMPs at day 3. Cells are subsequently exposed to RA and SAG to promote differentiation into ventral NPs and MNs. (B, C) Expression of NP (Pax6, Olig2, Nkx2.2, Sox1) and MN (Isl1/2) markers between day 4 and day 7 in differentiating ESCs. (D) RT-qPCR analysis of Irx3, Pax6, Nkx6.1, Olig2, and Nkx2.2 expression from day 3 to day 7 reveals progressive ventralization in response to increasing duration of Shh signaling. Underlying data are provided in S1 Data. (E) MN induction after day 5, revealed by RT-qPCR analysis of Sox1, Ngn2, Isl1, and Tubb3. Underlying data are provided in S1 Data. Scale bars = 40 μm. CHIR, Wnt pathway activator CHIR99021; ESC, embryonic stem cell; FGF, fibroblast growth factor 2; MN, motor neuron; NMP, neuromesodermal progenitor; NP, neural progenitor; N2B27, N2 and B27 media supplements; RA, retinoic acid; RT-qPCR, real-time quantitative polymerase chain reaction; SAG, Smoothened/Shh signalling agonist.
In vivo NPs respond to both the levels and duration of Shh signalling by transitioning through a succession of progressively more ventral gene expression states (S2A Fig) [6,41–43]. To further characterize the behavior of ESC-derived NPs in vitro, we first asked whether treatment of NPs with increasing concentrations of SAG (0, 10, 50, 100, 500, and 1,000 nM SAG) leads to progressively more ventral cell fates. Generation of NPs in the absence of SAG resulted in the expression of Pax3, Pax7, and Dbx1, indicative of a dorsal and intermediate NP identity (S2B Fig). Treatment with 10 nM SAG resulted in the down-regulation of these genes and induction of the pan-ventral marker Nkx6.1 (S2B and S2C Fig). Between 50 and 500 nM SAG, expression of the MN progenitor marker Olig2 was observed, while treatment with 500 and 1,000 nM SAG resulted in further ventralization and induction of the V3 interneuron progenitor (p3) determinant Nkx2.2 (S2B and S2D Fig). Induction of ventral markers coincided with the successive down-regulation of Irx3 and Pax6, consistent with their in vivo expression patterns (S2B and S2D Fig). Thus, in vitro NPs respond to different levels of Shh pathway activity by induction of the same progenitor markers that demarcate NP domains in the embryonic ventral spinal cord.
We next tested whether in vitro NPs also displayed progressive ventralization in response to increasing exposure durations to a constant concentration of SAG. To this end, we treated cells with 500 nM SAG from day 3 and quantified gene expression by real-time quantitative polymerase chain reaction (RT-qPCR) over the course of the next few days. At day 3.5, 12 hours after the cessation of Wnt signalling and addition of RA and SAG, cells expressed Sox1, Pax6, and Irx3, consistent with the acquisition of NP identity (Fig 1B and 1D). The absence of ventral markers at this stage indicates that these NPs initially adopt a dorsal/intermediate positional identity [42,43]. By day 4, the expression of Pax6 and Irx3 were maintained and Nkx6.1, which is expressed broadly in the ventral third of the neural tube, was induced (Fig 1B and 1D). Within 12 hours of this time point, Olig2 expression commenced and both Pax6 and Irx3 declined (Fig 1B–1D). Over the next 48 hours, Pax6 and Irx3 were further repressed, Nkx2.2 increased, and Olig2 expression began to decline (Fig 1B–1D). The order in which these genes were activated and repressed closely resembles the temporal-spatial sequence of progenitor domains in the embryonic spinal cord [6,42,43] and suggests that, under these conditions, MN progenitors are generated in vitro between days 4.5 and 6.
Consistent with the generation of MN progenitors in vitro, Ngn2 was induced following Olig2 (Fig 1E), and with an approximately 12-hour delay, we observed markers characteristic of postmitotic MNs, including Isl1 and Tubb3 (Fig 1C and 1E). Concomitantly, the expression of the NP marker Sox1 declined (Fig 1C and 1E). Taken together, these data indicate that this method of directing ESC differentiation recapitulates in vivo dynamics of neural tube patterning between embryonic days (e)8.5 and 10.5, and results in the production of MN progenitors and MNs characteristic of those normally found at forelimb levels.
Single cell transcriptome analysis of in vitro NPs
We reasoned that analyzing the transcriptome of individual cells would provide insight into the transitions in gene expression associated with the differentiation of MNs and allow the construction of a detailed developmental timeline. We therefore performed scRNA-seq analysis using the Fluidigm-C1 platform on 236 cells isolated from day 4 to day 6 of the differentiation protocol. After applying quality filters (see S1 Text), transcriptomes of 202 cells were retained for subsequent analysis (25 cells from day 4, 68 cells from day 5, and 109 cells from day 6). To identify the cell states present in the dataset, we established a data-driven analysis pipeline based on hierarchical clustering and association of gene modules with specific gene ontology (GO) terms (see S1 Text). In brief, the data were first filtered by removing genes that did not exceed a Spearman correlation of r > 0.4 with at least two other genes (retaining 2,287 genes). A combination of hierarchical clustering and automated selection criteria identified 22 gene modules that represent distinct patterns of gene expression across the dataset (see S1 Text and S1 Table). Further functional characterization of these gene modules based on GO terms resulted in the identification of 10 gene modules that were sufficient to assign a cell type classification to each cell in the dataset using hierarchical clustering (S3A Fig and S2 Table). Cells in these clusters showed comparable read counts and number of expressed genes per cell, suggesting that these properties did not bias the clustering (S3B Fig).
Consistent with our previous finding that the spinal NPs generated by differentiation of ESCs share a developmental lineage with trunk mesoderm , we observed two mesodermal cell populations in our dataset: paraxial presomitic mesoderm characterized by the expression of Meox1 and Foxc1 and a vascular endothelial population expressing Dll4 and Cdh5 (S3A Fig). The remaining cell clusters corresponded to different stages of NPs and differentiating MNs (Fig 2A). Five gene modules were associated with these cells (comprising 306 genes; see S1 Table). Module 1 was enriched for genes up-regulated in early NPs, including the TF Irx3. Module 2 contained genes expressed in MN progenitors, including the ventral progenitor markers Olig2 and Nkx6.1 and the neural-specific POU domain TF Pou3f2 (aka Brn-2). Module 3 comprised a set of genes transiently expressed in MNs as they differentiate, such as the bHLH TFs Ngn2, Neurod1, Neurod4, and Hes6; the homeodomain TFs Isl1 and Lhx3; and the Notch ligand Delta-like 1 (Dll1). Modules 4 and 5 revealed two successive waves of neuronal gene induction. Module 4 contained genes induced early in differentiated MNs such as Tubb3, the RNA-binding protein Elavl3 (aka HuC) and the SoxC TF Sox4, while Module 5 consisted of genes characteristic of more mature MNs, represented by choline acetyltransferase (Chat) and the TFs Isl2 and Onecut1 [44–47].
Fig 2. Reconstruction of transcriptional changes during MN differentiation.
(A) Identification of NP cell states using hierarchical clustering of gene expression profiles of the individual cells. (B) Cell state graph constructed from minimum spanning trees, color coded for the cell populations identified in Fig 2A. Stars indicate start and end cells for the reconstruction of transcriptional changes along pseudotime. Shading of edges between cells indicates how often the edge was used in the reconstruction of gene expression along pseudotime (see S1 Text). (C) Cell state graph color coded for expression levels of Irx3, Olig2, Ngn2, Lhx3, Isl1, and Chat. (D) Inferred changes in gene expression over pseudotime from 9,000 shortest paths connecting start and end cells (stars in B). Each shortest path was resampled to a length of 41 pseudo–time points to enable statistical measurements of gene expression. Cell IDs are color coded according to cell states in (A). Quantification of the global rate of change in gene expression identifies three metastable states (light gray) separated by transition states, during which the rate of change in gene expression is increased (dark gray). Transition phases are defined as intervals along the pseudo-temporal timeline at which the second derivative of the global gene variation is negative, while metastable states are characterized by a positive second derivative. (E) Gene expression profiles along pseudo-time for NP TFs (Irx3, Pax6, Nkx6.1, and Olig2), genes associated with the transition to MNs (Ngn2, Lhx3, and Neurod4) and MN markers (Isl1/2, Tubb3, and Chat). (F) Levels of gene expression for Hes1/5, Olig2, and Ngn2 over pseudotime. Note that Olig2 expression appears biphasic, with up-regulation of Olig2 concommitant to Ngn2 induction and repression of Hes1/5 in the transition phase from NP to MN. Chat, choline acetyltransferase; KSP, K shortest paths; MN, motor neuron; NP, neural progenitor; pMN, MN progenitor; p3, V3 interneuron progenitor; Tubb3, neuronal class III beta-tubulin.
Whereas the five cell clusters defined by these modules represented a progressive shift of cell states from early progenitor cells to MNs, the remaining cell cluster exhibited a divergent gene expression signature. In this cluster, many genes contained in Modules 1 and 2 were down-regulated, but neuronal gene expression was not increased. This cluster exclusively consisted of day 6 cells (Fig 2A). Nkx2.2 could be detected in some cells of this cluster (S3C Fig), suggesting that it was composed of cells progressing from a pMN to a more ventral p3 identity. Further differential gene expression analysis on this population identified fatty-acid binding protein 7 (Fabp7) as enriched in these cells (S3C Fig). Fabp7 levels are markedly up-regulated in p3 progenitors in vitro and at cervical and brachial levels in embryonic spinal cords at e10.5 (S3D Fig). We therefore conclude that this cluster contains cells progressing from MN to p3 progenitors. Taken together, this suggests our scRNA-seq analysis identifies cells along the MN developmental timeline and partitions these into specific cell types from early NPs to postmitotic MNs and p3 progenitors.
We next asked whether it was possible to reconstruct the developmental timeline from the transcriptome data. For this, we used the 306 genes contained in the five neural gene modules to visualize the developmental trajectory as a pseudo-temporal ordering derived from a consensus of a large number of randomized minimum spanning trees (see S1 Text). The resulting cell graph represents the predicted developmental order of cells based on their transcriptome profile and, hence, differentiation state (Fig 2B). Strikingly, the five previously characterized cell clusters were ordered on the cell state graphs as expected from the characterization of their gene expression profile (Fig 2C). The graph revealed developmental trajectories originating from Irx3 expressing early NPs to MN progenitors characterized by Olig2 expression (Fig 2C). These progenitors then differentiated into MNs via the sequential expression of Ngn2, Lhx3, Isl1, and Chat (Fig 2C) or into p3 progenitors characterized by Nkx2.2 and Fabp7 expression (S3C Fig). To investigate these trajectories in more detail, we focused on the developmental trajectory leading from NPs to MNs. To represent changes in gene expression in an unbiased manner, we reconstructed the average gene expression program along pseudotime from the 9,000 shortest paths connecting Irx3-expressing progenitor cells to differentiated MNs on the cell state graph (starred cells in Fig 2B, see S1 Text). Each individual path was resampled to a constant length of 41 pseudotime points (Fig 2D), allowing statistical measurements along the developmental timelines. The outcome was predicted gene expression dynamics during MN differentiation.
Characterization of transcriptional changes during MN differentiation
As a first validation, we asked if the pseudo-temporal ordering reproduced the temporal sequence of well-characterized gene expression changes that lead to MN differentiation. The inferred trajectory correctly predicted the induction sequence of the homeodomain and bHLH TFs Irx3, Pax6, Nkx6.1, and Olig2 involved in ventral patterning of the spinal cord (Fig 2E) [41–43]. Next, we focused on the transition from progenitors to MNs. As expected, this transition was associated with the transient expression of Ngn2, Neurod4, and Lhx3, followed by the expression of MN markers, including Isl1/2, Tubb3, and Chat (Fig 2E). To assess the robustness of these gene expression dynamics, we utilized a bootstrapping approach to ask how dependent these are on individual cells with particular gene expression values (see S1 Text). A total of 1,000 bootstrapped datasets were constructed by randomly drawing cells, with replacement, while maintaining original sample size. Then, expression profiles were calculated for each gene in each replicate (S4 Fig). To statistically quantify their robustness, we asked how well these profiles were correlated between each pair of replicates (see S1 Text). This analysis revealed a mean Spearman correlation value greater than 0.85 for most genes (S4 Fig). This suggests that the observed gene expression dynamics do not depend on the levels of gene expression in specific cells along the pseudo-temporal trajectory and are a robust representation of the gene expression dynamics during MN differentiation.
The process of cell development has been characterized as a series of metastable states defined by a relatively homogenous gene expression program connected by stereotypic transitions . During these transitions, coordinated changes in gene expression occur, often induced in response to a change in signalling. We reasoned that metastable states and transition phases should be evident in the pseudo-temporal ordering. Quantifying the variation in gene expression by averaging the normalized derivative of the most dispersed genes’ expression profiles identified these phases (Fig 2D). The three metastable states in which gene expression changes were relatively modest corresponded to early NPs, MN progenitors, and MNs. Linking these states were transitions characterized by an increased change in the global gene expression profile. The first transition corresponded to the switch from Irx3-expressing intermediate progenitors to Olig2-expressing MN progenitors (Fig 2E), while the second captured the transition of progenitors to postmitotic neurons (Fig 2E).
We asked whether signatures of signalling pathways driving these transitions could be identified. To this end, we examined the induction and disappearance of canonical target genes for different signalling pathways. As expected, the transition from Irx3 to Olig2 coincided with the induction of well-known Shh target genes Ptch2, Hhip1, and Gli1, consistent with Shh signalling mediating this transition (S3E Fig). By contrast, the second transition was accompanied by a loss of Notch signalling, marked by the disappearance of Hes1/5 and induction of markers causing or characteristic of a loss of Notch signalling, including Numbl, Hes6, Dll1, Ngn2, and Neurod4 (Fig 2E and S3F Fig). Strikingly, the beginning of this stage coincided with peak expression levels of Olig2 (Fig 2E and 2F). This finding raised the possibility that high levels of Olig2 promote neurogenesis, potentially by directly regulating levels of Notch signalling. In summary, the characterization of changes in the transcriptional profile in pseudotime identified distinct metastable cell states and the signalling pathways associated with the transitions between these states.
In vitro and in vivo validation of the pseudo-temporal ordering
To extend this approach and validate the predicted timeline, we asked whether the data were sufficient to capture fine-grained temporal information that could be tested experimentally. Examination of the transition from Olig2-expressing progenitors to Isl1-expressing MNs predicted the transient expression of first Ngn2, then Lhx3, and finally Isl1 (Fig 2C and 2E). This is consistent with in vivo data indicating that Lhx3 precedes the expression of other MN markers in the spinal cord [47,49], and a similar sequence of gene expression has been described in an in vitro MN differentiation protocol based on embryoid bodies [45,50]. To confirm this sequence of events in vitro, we assayed Olig2, Ngn2, Lhx3, and Isl1 on day 6 of differentiation and quantified the levels of expression in individual nuclei (S5A Fig). Comparison of Olig2 and Isl1 levels in individual nuclei revealed a clear trajectory from Olig2-positive, Isl1-negative NPs to Isl1-positive, Olig2-negative MNs. Overlaying the levels of Ngn2 and Lhx3 in the same cells revealed that both proteins are only transiently expressed along the differentiation trajectory (S5B and S5C Fig). To confirm the absence of Lhx3 in more mature MNs, we assayed Lhx3, Isl1, and the pan-neuronal marker Tubb3 (S5D Fig). Consistent with the pseudo-temporal ordering, most Tubb3-expressing cells displayed high levels of Isl1 expression but only low levels of Lhx3, while cells with high levels of Lhx3 did not express high levels of Isl1 or Tubb3 (S5D Fig). In summary, these two observations confirm the predictions from the pseudo-temporal ordering and validate the approach for predicting fine-grained changes in the transcriptional program of cells along the differentiation trajectory to MNs.
To further test the reliability of the timeline and demonstrate the validity of the approach for understanding MN differentiation dynamics, we asked if we could predict novel genes involved in MN formation. To this end, we selected genes positively correlated with Olig2 and Ngn2 (S5E and S5F Fig). One gene with a particularly strong relationship was Zbtb18 (also known as RP58 or Zfp238). Zbtb18 is a zinc-finger TF with a BTB domain. In the brain, its loss causes microcephaly and decreased neuronal and increased glial differentiation . Less is known about its expression pattern and role in the spinal cord, although in situ hybridization analyses have suggested it is predominantly expressed in ventral progenitors . As expected, when we assayed Zbtb18 using immunohistochemistry, it was expressed in cells that also expressed Olig2 and Ngn2 (S5G–S5I Fig). Consistent with this, its expression was detected in vivo in the pMN domain at e9.5 in cells that also expressed high levels of Olig2 and Ngn2 (S5J Fig). At e10.5, it was still predominantly expressed ventrally, although no longer confined to the pMN domain (S5K Fig). In summary, this expression pattern further validates the computationally reconstructed MN differentiation timeline.
Olig2 expression increases as cells commit to MN differentiation
The MN differentiation timeline indicated that Olig2 expression was induced as Irx3 was repressed (Fig 2E), consistent with the cross-repressive interactions between these two genes [13,17,53]. This transition demarcated the transition from the first to the second phase identified in the MN timeline. It was noticeable that the expression of Olig2 appeared biphasic with a marked increase in levels of Olig2, which coincided with the transition from the second to the third phase. Moreover, this transition corresponded to the induction of Ngn2. This predicted that Olig2 levels peak at the onset of differentiation before being down-regulated as MN identity is elaborated (Fig 2E). To test this prediction, we first examined the levels of Olig2 and Ngn2 in the neural tube of e9.5 and e10.5 embryos during the period of MN production (Fig 3A–3D″). Consistent with previous studies, we found that at both stages, a proportion of Olig2-expressing cells also expressed Ngn2, while a much lower proportion of cells expressed Ngn2 outside the pMN domain [13,17,19]. To test if the levels of Olig2 expression varied in the way predicted by the scRNA-seq data, we quantified levels of Olig2, Ngn2, and the MN marker Mnx1 in nuclei of the pMN domain (Fig 3E–3H). This revealed a striking correlation between Olig2 and Ngn2 protein levels in individual cells throughout the pMN domain (Fig 3E and 3G). Moreover, cells expressing high levels of Olig2 and Ngn2 were differentiating into MNs, as measured by the induction of Mnx1 (Fig 3H). This quantification also indicated that Olig2 protein persisted longer than Ngn2 in MNs, as cells coexpressing high levels of Olig2 and Mnx1, but not Ngn2, were observed (Fig 3H). Taken together, these data suggest that high levels of Olig2 correspond to the induction of Ngn2 and the onset of neurogenesis within the pMN domain.
Fig 3. Olig2 expression is higher in Ngn2-expressing progenitors in the pMN domain.
(A–B″) Staining for Ngn2 (A, B) and Olig2 (A′, B′) merged with Mnx1 (green in A″, B″) in spinal cords at e9.5 (A–A″) and e10.5 (B–B″). (C–D′) Higher magnification images of the spinal cords shown in (A–B″). Red arrowheads indicate nuclei with elevated levels of Ngn2 and Olig2. (E, G) Positive correlation between Olig2 and Ngn2 protein levels in individual nuclei at e9.5 (E) (n = 464 nuclei) and e10.5 (G) (n = 1,078 nuclei). Underlying data are provided in S1 Data. (F, H) Levels of Olig2, Mnx1, and Ngn2 in individual nuclei throughout the pMN domain at e9.5 (F) and e10.5 (H). Plotting Olig2 versus Mnx1 protein levels reveals a clear differentiation trajectory from Olig2-positive pMN cells to Mnx1-positive MNs. Note that high levels of Ngn2 are only observed in cells with high levels of Olig2 expression. In addition, Olig2 protein perdures much longer in Mnx1-positive MNs than Ngn2. Underlying data are provided in S1 Data. Scale bars = 50 μm. e, embryonic day; MN, motor neuron; pMN, MN progenitor.
These data prompted us to test directly whether progenitors that expressed high levels of Olig2 were committed to MN differentiation. Since endogenous Olig2 protein disappears rapidly from differentiated MNs, we took advantage of an ESC line in which we fused monomeric far-red fluorescent protein Katushka-2 (mKate2) to the C-terminus of endogenous Olig2 via a self-cleaving peptide (Fig 4A) [54,55]. In these cells, the expression of mKate2 provides a readout of Olig2 levels, but the increased stability of the fluorescent protein offers a way to mark the progeny of Olig2-expressing cells and estimate Olig2 levels in the progenitor. Control ESC differentiations indicated that Olig2 expression dynamics, protein levels, and MN formation were similar in cells containing the engineered or wild-type Olig2 allele (Fig 4B and S6A–S6C Fig). Quantification of the mKate2 and Olig2 protein levels in individual nuclei revealed a positive correlation in most cells (Fig 4C–4C″ and 4G and S6D–S6F Fig). However, we noted a cohort of cells with much higher levels of mKate2 relative to Olig2. Assaying Isl1/2 expression revealed that these cells were MNs (Fig 4D–4D″). Consistent with this, high levels of Isl1/2 and Tubb3 expression were only detected in cells with high levels of mKate2 (Fig 4H and S6A–S6C Fig). Moreover, mKate2 levels negatively correlated with levels of the NP marker Sox1 (Fig 4E–4E″ and 4I). Thus, MNs indeed progress through a distinct Olig2HIGH state as they exit from the NP state.
Fig 4. Quantification of a fluorescent Olig2 reporter reveals a marked up-regulation of Olig2 prior to MN differentiation.
(A) Design of the Olig2-mKate2 reporter. A 3xNLS-FLAG-mKate2 reporter was fused to the C-terminus of endogenous Olig2 via a T2A self-cleaving peptide. (B) Western blot analysis reveals that the targeted allele shows the same expression dynamics and levels as endogenous Olig2. The targeted allele runs at slightly increased molecular weight due to addition of the T2A peptide (see A). Note that both alleles are targeted in this cell line and, consequently, no protein of wild-type size was detected. (C-F″) Immunofluorescence for mKate2 with Olig2 (C–C″), Isl1/2 (D–D″), Sox1 (E–E″), and Nkx2.2 (F–F″) at day 6 of the differentiations. (G–J) Quantification of protein levels of mKate2 and Olig2 (G, n = 2,851 nuclei), Isl1/2 (H, same dataset as G), Sox1 (I, n = 2,049 nuclei) and Nkx2.2 (J, n = 2,034 nuclei) in individual nuclei. Note the positive correlation between mKate2 and Olig2 and Isl1/2, and the negative correlation between mKate2 and Sox1 and Nkx2.2. Underlying data are provided in S1 Data. (K) Inhibition of Notch signaling using 10 ng/μl DBZ causes an increase of neurogenesis. Immunofluorescent staining for Olig2, mKate2, and Tubb3 in control or after 24 hours DBZ treatment at day 6 of the differentiation. (L) Frequency plots of mKate2 fluorescence intensity obtained by flow cytometry reveal a strong increase in the number of mKate2HIGH cells after 24 hours DBZ treatment. Scale bars = 25 μm. DBZ, dibenzazepine; FLAG, FLAG epitope tag; mKate2, monomeric far-red fluorescent protein Katushka-2; MN, motor neuron; Tubb3, neuronal class III beta-tubulin; T2A, Thosea asigna virus 2A peptide sequence; 3xNLS, 3 copies of a nuclear localization sequence peptide.
To address whether the transient up-regulation of Olig2 expression was specific for the transition from pMN cells to MNs, we quantified levels of mKate2 in Nkx2.2-expressing p3 progenitors (Fig 4F–4F″). During development, these progenitors transit through an Olig2-expressing pMN intermediate state before losing Olig2 expression and inducing Nkx2.2 [41,42,56]. In contrast to the positive correlation between Isl1/2 and mKate2 (Fig 4H), cells expressing high levels of Nkx2.2 had low or undetectable levels of mKate2 expression (Fig 4J and S6F Fig). Thus, distinct Olig2 expression dynamics underlie the progression of pMN cells to MNs and p3 progenitors.
Inhibiting Notch signaling increases Olig2 expression
These observations raise the question of what up-regulates Olig2 prior to MN formation. The Notch signalling pathway is implicated in controlling the rate of neurogenesis, and inhibition of Notch signalling in NPs is well known to trigger neuronal differentiation [25,57–59]. Furthermore, the inferred MN differentiation trajectory indicated that two canonical effectors of the Notch pathway, Hes1 and Hes5, decreased as cells switched from the early phase of Olig2LOW expression to Olig2HIGH. We therefore tested whether inhibiting Notch signalling up-regulated Olig2. As expected, inhibition of Notch signalling through the addition of the γ-secretase inhibitor dibenzazepine (DBZ) for 24 hours between day 5 and day 6 of differentiation caused a substantial increase in the number of neurons observed (Fig 4K, compare S6C and S6I Fig). Quantifying mKate2 levels using flow cytometry revealed a similarly substantial increase in the number of cells expressing high levels of mKate2 (Fig 4L and S6G Fig). Furthermore, co-staining these cells with the pan-neuronal marker Tubb3 revealed that most of the mKate2HIGH cells were neurons (S6I Fig). To test whether the increase in mKate2 fluorescence is due to increased Olig2 expression upon Notch inhibition, we quantified mRNA levels of Olig2 and other progenitor and neuronal markers using RT-qPCR after 0, 12, and 24 hours of Notch inhibition (S6J–S6L Fig). In contrast to other progenitor markers (Hes1/5, Sox2, Pax6), which decreased upon Notch inhibition (S6J Fig), Olig2 levels peaked at 12 hours before decreasing after 24 hours (S6K Fig). The observed Olig2 expression dynamics are strikingly similar to those of other genes previously implicated in MN formation, including Ngn2 and Pou3f2 (S6K Fig). Consistent with the increase in the expression of neurogenic markers after 12 hours, we also observed an increase in the expression of neuronal genes 12 hours and 24 hours after Notch inhibition (S6L Fig). Taken together, these data suggest that Notch signalling controls the transition between the distinct phases of Olig2 expression by restraining Olig2 expression in MN progenitors.
Olig2 represses the expression of Hes1/Hes5
To test whether the up-regulation of Olig2 coincided with the down-regulation of Hes1 and Hes5 in vivo, we examined the expression of these proteins in mouse embryos. Hes1 is broadly expressed by dorsal progenitors that express the homeodomain protein Pax3, as well as floor plate and p3 cells, marked by the expression of Foxa2 and Nkx2.2, respectively (Fig 5B, 5C, 5F and 5G) . By contrast, Hes5 is expressed by cells in the intermediate spinal cord, marked by the expression of Irx3 and Pax6 (Fig 5B, 5D and 5E). Olig2 expression was first detectable at e8.5, a time at which Hes5 was broadly expressed throughout the ventral neural tube (Fig 5M). Shortly thereafter, Olig2 and Hes5 showed a high degree of coexpression, which coincided with an increase in the number of Olig2-expressing MN progenitors (Fig 5N–5O). However, coexpression of Olig2 and Hes5 appeared to be transient, as by e9.5, Hes5 was down-regulated in many Olig2+ cells, particularly at rostral levels (Fig 5P), and few coexpressing cells could be found by e10.5 (Fig 5Q). During this time, Hes1 expression was low or absent in most MN progenitors (Fig 5H–5L). The progressive decrease in Hes5 expression from Olig2+ cells was mirrored by a reciprocal increase in Ngn2 expression and, subsequently, the exit of these cells from the cell cycle and the onset of MN differentiation marker expression (Fig 5M′–5Q′, [13,17,28]). Thus, the transient coexpression of Olig2 and Hes5 in vivo marks the pMN state, while the clearance of Hes5 from Olig2-positive cells coincided with the onset of Ngn2 expression and MN differentiation.
Fig 5. Olig2 and Hes are dynamically expressed in the mouse neural tube.
(A–D) Expression patterns of Ngn2 (green in A), Olig2 (red in A, C, D), Hes1 (red in B, green in C), and Hes5 (green in B, D) in the neural tube at e10.5. Note the low expression levels of Hes1/5 and high expression levels of Ngn2 in the pMN domain (compare A, B). (E) Hes5 (green) expression coincides with the expression of high levels of Pax6 (red) in the intermediate neural tube. (F, G) Hes1 expression (green) is readily detected in both Nkx2.2+ p3 progenitors (red in F) and floor plate cells labelled by Foxa2 expression (red in G). (H–Q′) Time course of Olig2 (blue), Hes1 (red), Hes5 (red), and Ngn2 (green) expression in neural tubes between e8.5 and e10.5. Multiple panels shown for e9.5 reflect developmental progression from caudal to rostral positions along the neuraxis. Hes1 expression appears to recede from the ventral neural tube upon the onset of Olig2 expression at e8.5 (H) and is thereafter absent from most Olig2+ cells (I–L). Olig2 and Hes5 are initially coexpressed (M, N). Over time, Hes5 expression progressively disappears from the pMN domain (N–Q), and Ngn2 concomitantly increases (N′–Q′). Insets show single channel images of the outlined area for the respective markers. Scale bars = 50 μm. e, embryonic day; pMN, MN progenitor; p3, V3 interneuron progenitor.
We next asked whether Olig2 might be responsible for the repression of Hes1 and Hes5 using Olig2Cre knock-in mice [56,61]. In these mice, the Olig2 coding sequence has been replaced with bacteriophage P1 Cre recombinase (Cre) . Thus, Cre protein expression demarcates the presumptive pMN domain in heterozygous control and in homozygous Olig2Cre/Cre mutant embryos, which entirely lack Olig2 activity (Fig 6A and 6E). In controls, the pMN domain was flanked dorsally and ventrally by Hes5 and Hes1 expression, respectively, with little overlap of Cre with either protein (Fig 6C–6D′ and 6I). By contrast, Olig2 mutant spinal cords displayed a marked dorsal expansion of Hes1 and ventral expansion of Hes5 into the pMN domain, such that their expression domains appeared to contact one another (Fig 6G–6I). This juxtaposition was associated with a substantial decrease in the number of cells expressing Ngn2 in the pMN domain (Fig 6B, 6F and 6I). Thus, Olig2 is required to maintain the boundaries of Hes1 and Hes5 and allow Ngn2 to accumulate within MN progenitors.
Fig 6. Repression of Hes1/5 in the pMN domain depends on Olig2 activity.
(A–D) Expression of Cre (green in A–D), Olig2 (red in A), Ngn2 (red in B), Hes1 (red in C, grey in C′), and Hes5 (red in D, grey in D′) in e10.5 Olig2Cre heterozygous embryos. (E–H) In Olig2Cre/Cre homozygous mutants, Hes1 expands dorsally (G, G′) and Hes5 ventrally (H, H′) into the pMN domain, marked by Cre expressed from the Olig2 locus. The expansion of Hes1/5 coincides with a loss of the high levels of Ngn2 normally seen in the pMN domain. (I) Quantification of Hes1, Hes5, and Ngn2 expression in Olig2Cre heterozygous and homozygous embryos. The overlap between Cre and Hes1/5 significantly increases in Olig2Cre homozygotes, while overlap between Ngn2 and Cre is strongly reduced. Plot shows the mean ±SEM from multiple sections collected from 3–5 embryos for each group. Each section is represented by a single dot, with n = 8–11 for each group. Underlying data are provided in S1 Data. **** p < 0.0001, unpaired t test. (J–S) Electroporation of myc-tagged OLIG2 and an OLIG2-bHLH-Engrailed repressor domain fusion protein in chick neural tubes represses expression of the Hes5 homologues HES5-1–HES5-3 (K–M; Q–R) and the Hes1 homologue HAIRY1 (N, S). “+” indicates transfected side of the spinal cords. Results are representative of >5 successfully transfected embryos collected from two or more experiments. Scale bars = 50 μm. bHLH, basic helix-loop-helix DNA binding domain; Cre, bacteriophage P1 Cre recombinase; e, embryonic day; EnR, Engrailed transcriptional repression domain; EP, electroporation; pMN, MN progenitor.
To address whether Olig2 expression was sufficient to repress Hes1 and Hes5, we used in ovo electroporation to deliver retroviral expression constructs driving the expression of a myc epitope-tagged form of OLIG2 into the developing spinal cord of Hamburger-Hamilton (HH) stage 11–13 chick embryos. These conditions have been previously shown to increase NGN2 expression . Whereas mammals have a single Hes5 gene, birds contain three Hes5 paralogs, termed HES5-1, HES5-2, and HES5-3 , clustered at a common genomic locus (Fig 7A). When Olig2 was misexpressed, all three chick HES5 genes were substantially reduced, as was the chick Hes1-related gene HAIRY1 (Fig 6J–6N). Similar results were achieved with the misexpression of a dominant repressor form of OLIG2 containing its bHLH DNA-binding domain fused to a heterologous Engrailed transcriptional repression domain (Fig 6O–6S; ). Based on these results, we conclude that Olig2 expression suffices to repress Hes gene expression in NPs.
Fig 7. Olig2 binds to an evolutionarily conserved element near Hes5.
(A) Identification of an evolutionarily conserved element containing an E-box in the vicinity of the Hes5 genomic locus in chick, mouse, and human (Hes5(e1)). (B) Analysis of Olig2 Chip-Seq data from  reveals Olig2 binding sites in the vicinity of the Hes1 and Hes5 genes. The peak corresponding to the Hes5(e1) element is highlighted in red. (C) Electrophoretic mobility shift assays show that both Olig2 and E12 homodimers can individually bind to the Hes5(e1) E-box and do not form any heterodimeric complexes (lanes 1–4). Positions of the different protein complexes are indicated by colored arrows. Binding depends on the E-box, as both proteins fail to bind probes containing an E-box mutation (Hes5(e1ΔE)) (lanes 5–7). Olig2 binding to Hes5(e1) can be abolished by the addition of unlabelled Hes5(e1) probes, but not those containing the E-box mutation (lanes 8–14). (D) Id1 inhibits binding of E12, but not of Olig2 or Ngn2, to the Hes5(e1) element. Olig2, E12, and Ngn2 alone or Ngn2/E12 heterodimers can bind the Hes5(e1) element. Mixing Olig2 or Ngn2 with Id1 does not inhibit their homodimeric binding activities (lanes 2, 5, 8, and 10). In contrast, Id1 strongly inhibits binding of both E12/E12 and Ngn2/E12 complexes (lanes 6 and 10). The addition of E12 without and with Id1 does not affect Olig2 binding efficiency (lanes 2, 4, and 7). ATG, translational initation codon; Chip-Seq, chromatin immunoprecipitation-sequence; E-box, bHLH transcription factor binding site; N2, Ngn2 protein; O2, Olig2 protein.
Hes proteins and proneural bHLH TFs such as Ngn2 act antagonistically in multiple developmental contexts [25,26], and ectopic Olig2 expression has been shown to promote Ngn2 expression . Two sequences of events could explain the repression of Hes genes and induction of Ngn2 in MN progenitors. Either Olig2 represses Hes1/5 and thereby indirectly induces Ngn2 or, alternatively, Olig2 induces Ngn2, which then antagonizes the expression of the Notch effectors. To distinguish between these possibilities, we investigated if Olig2 represses Hes5 in Ngn2 null mutants (S7A–S7F″ Fig). To this end, we utilized a Ngn2 knock-in GFP (Ngn2KIGFP) mouse line, in which an internal ribosome entry sequence-Green Fluorescent Protein (IRES-GFP) construct has been inserted into the coding sequence of Ngn2 . Assays at e9.5 and e10.25 revealed that Olig2 expression was maintained (S7A‘–S7F’ Fig) and Hes5 repressed in MN progenitors lacking Ngn2 (S7A”–S7F” Fig). Consistent with the observed repression of Hes5, GFP expression from the endogenous Ngn2 locus was strongly elevated in MN progenitors (S7C–S7F Fig). We therefore conclude that Olig2, not Ngn2, is the main repressor of Hes genes in MN progenitors.
These findings raise the question of how important the Olig2-mediated repression of Hes genes is for the pattern of neurogenesis in the ventral spinal cord. To address this, we investigated the consequences of preventing Olig2-mediated repression by ectopically expressing chick HES5-2, an ortholog of murine Hes5, in the ventral spinal cord of chick embryos. Ectopic HES5-2 expression did not have a noticeable effect on the levels of the progenitor markers SOX2, NKX6.1, PAX6, and OLIG2 (S7G–S7J Fig). By contrast, and consistent with the well-known antineurogenic role of Hes proteins [24,60,65], ectopic HES5-2 resulted in the down-regulation of the proneuronal TFs NGN2 and NEUROD4 and of the pan-neuronal marker NEUN (S7K–S7Q Fig). Of note, cells that maintained NGN2 and NEUROD4 in these experiments usually contained little if any GFP, marking transfected cells, suggesting they were not electroporated (S7N and S7O Fig). Consistent with this, only a minor fraction of GFP-electroporated cells left the ventricular zone and activated NEUN expression (compare S7M and S7P Fig). These results suggest that the repression of Hes genes is the key mechanism by which Olig2 promotes neurogenesis and that the antineurogenic function of Hes proteins needs to be overcome before Ngn2 can be induced and neurogenesis initiated. Together, these data indicate that Olig2 plays a critical role in repressing the expression of Hes genes within MN progenitors to promote the expression of proneurogenic bHLH TFs, such as Ngn2 and Neurod4, and thereby increases the rate of neuronal differentiation in MN progenitors.
Olig2 acts directly on a Hes5 regulatory element
The striking effects of Olig2 on Hes1 and Hes5 expression prompted us to ask whether Olig2 might directly regulate these genes. Examination of chromatin immunoprecipitation data from mouse NPs revealed several prominent binding sites of Olig2 in the vicinity of the two loci (Fig 7B) (, http://www.ebi.ac.uk/ena/data/view/ERX628418). Furthermore, some of these binding sites are in close proximity to previously mapped binding sites for the Notch signalling cofactor RBPJ . Bound regions included sites close to the transcription start sites of the genes and in putative distal regulatory elements (Fig 7B). Aligning genomic sequences of the Hes5 locus from chick, mouse, and human indicated that one of the binding sites for Olig2 and RBPJ coincided with a highly conserved, approximately 200–base pair element, hereafter termed Hes5(e1), that is 80% identical between mouse and human and 53% identical between chick and mouse (Fig 7A and 7B). This element is 7.9 kilobases (kb) 5′ to the transcriptional start site for Hes5 in mouse, 10.5 kb 5′ to the transcriptional start site in human, and in the middle of the HES5 gene cluster in chick.
Like many bHLH proteins, Olig2 binds to canonical bHLH transcription factor binding site (E-box) DNA response elements with the palindromic sequence CANNTG . This motif was found within the most conserved central region of the Hes5(e1) element (87% identity between chick and mouse over 46 bp; 98% identity between mouse and human) (Fig 7A). To confirm that Olig2 could bind to the Hes5(e1) element, we performed in vitro binding experiments using a probe comprising the conserved central region. In vitro translated Olig2 readily bound to the Hes5(e1) E-box, as did other bHLH proteins such as E12 and Ngn2 (Fig 7C and 7D). These binding activities were abolished when the conserved E-box sequence was mutated (Fig 7C). To test if Olig2 binding activity is enhanced by the presence of E proteins, we mixed Olig2 protein with E12 but found no evidence of either an Olig2:E12:DNA complex or enhanced binding affinity to the Hes5(e1) E-box (Fig 7C). In addition, mixing Olig2 with Id1, a potent competitor for E protein binding, did not diminish Olig2 binding, although mixing E12 and Ngn2 with Id1 completely abolished both E12/E12 and Ngn2/E12 binding activities (Fig 7D). The binding of both Olig2 and Ngn2 to Hes5(e1) in vivo was further confirmed through chromatin immunoprecipitation experiments (S8A Fig). Taken together, these data indicate that Olig2 homodimers bind directly to a highly conserved Hes5(e1) regulatory element through a single E-box site that may be targeted by other bHLH proteins.
The Hes5(e1) element restricts gene expression from the pMN
The observation that Olig2 could bind to the conserved element within the Hes5 locus prompted us to test whether this element restricted gene expression selectively from the pMN domain. To test this, we generated reporter constructs consisting of Hes5(e1), with or without an intact E-box, upstream of a β-globin minimal promoter driving expression of a nuclear enhanced GFP gene (Hes5(e1)-βG::nEGFP) (Fig 8C). We co-electroporated these constructs into the chick spinal cord, together with a plasmid encoding nuclear β-galactosidase (βgal) driven by the CMV/β-actin promoter (Fig 8A, 8B, 8D and 8E). βgal expression appeared to be uniform throughout the dorsal-ventral axis of the neural tube (Fig 8A, 8D and 8F). By contrast, Hes5(e1)-βG::nEGFP activity was spatially restricted, with high levels of expression in the intermediate portions of the neural tube but little, if any, expression in both ventral and dorsal regions (Fig 8B and 8F). Strikingly, the ventral limit of Hes5(E1)-βG::nEGFP expression coincided with the dorsal border of the Olig2 expression domain (Fig 8B and 8F).
Fig 8. The Hes5(e1) element is required for repression of reporter genes in the pMN domain.
(A, B) Co-electroporation of CMV/β-actin::nLacZ and Hes5(e1) reporter plasmids into chick spinal cord. Although electroporation (revealed by β-Gal antibody staining, magenta in [A]) is uniform along the dorsal-ventral axis, expression of the EGFP reporter is confined to intermediate parts of the neural tube (A, B), and little coexpression of Olig2 and EGFP was detected (B). (C) Design of Hes5(e1) and Hes5(e1ΔE) reporters. The Hes5(e1) element was cloned in front of β-globin minimal promoter to drive EGFP reporter gene expression. To test the importance of the E-box in the Hes5(e1) element, critical base pairs for Olig2 binding were mutated (red). (D, E) Co-electroporation of CMV/β-actin::-nLacZ and Hes5(e1ΔE) reporter plasmids into chick spinal cord. In contrast to the Hes5(e1) reporter plasmid, significant coexpression of Olig2 and GFP in the pMN domain is detected (E). Note that E-box mutation reduced the basal activity of the reporter such that longer exposure times were needed to achieve the signal levels seen in the intermediate spinal cord with the nonmutated Hes5(e1) reporter (S8B and S8C Fig). (F) Scatter dot plots display the dorsal-ventral positions (distance from the roof plate) of individual cells expressing the Hes5(e1) and Hes5(e1ΔE) reporters, relative to CMV/β-actin::-nLacZ and Olig2. Results are aggregated from five representative sections taken from five well-electroporated and stage-matched spinal cords. The Hes5(e1ΔE) reporter exhibits a significant ventral shift in its activity and considerable overlap with Olig2 expression (blue dotted box). Box plots include the median and whiskers represent 5th and 95th percentiles. Data points that lay outside the DV scale used to assess these experiments were excluded from this analysis. ** p = 0.0005, Mann-Whitney test; p = 0.6649. Underlying data are provided in S1 Data. (G) EGFP expression in Hes5(e1)-nEGFP whole mount embryos at e10.5. (H–H″) Cryosections of Hes5(e1)-nEGFP embryos at e10.5 assayed for GFP, Olig2, and Hes5. EGFP expression colocalizes with Hes5 expression (H″) but not with Olig2 (H). (I–N) Hes5(e1)-nEGFP expression in Olig2 heterozygous (I, K, L) and homozygous mutants (J, M, N). In Olig2 heterozygotes, little nEGFP expression can be detected in the Olig2 expression domain, resulting in a pronounced gap between the expression domains of EGFP, Nkx2.2, and Hes1 (K, L). By contrast, the EGFP, Nkx2.2, and Hes1 expression domains directly abut each other in Olig2 homozygous mutants (M, N). β-Gal, beta-galactosidase; βGlob, beta-globin; CMV/β-actin::nLacZ, cytomegalovirus/chick beta-actin promoter driving nuclear LacZ gene expression; ΔE, E-box deletion; E-box, bHLH protein binding site; EGFP, enhanced green fluorescent protein; FP, floor plate; GFP, green fluorescent protein; H5(e1), Hes5(e1) genomic element; ns, not significant; pMN, motor neuron progenitor; WT, wild-type.
To determine whether the E-box within Hes5(e1) was essential for this spatially restricted expression pattern, we compared the activity of a Hes5(e1)-βG::nEGFP reporter construct in which the E-box had been mutated (Hes5(e1ΔE)-βG::nEGFP) (Fig 8D–8F). Loss of the E-box substantially reduced the overall activity of the nEGFP reporter, compared to the original construct (S8B and S8C Fig). In addition, nEGFP expression now showed abundant overlap with Olig2 in the ventral spinal cord (Fig 8E and 8F). Together, these results indicate that the Hes5(e1) element integrates both positive and negative regulatory information through its E-box.
Finally, to test whether Olig2 is responsible for the restriction of Hes5(e1)- βG::nEGFP from the pMN domain, we generated transgenic mice containing this construct that displayed activity throughout the neuraxis (Fig 8G). In agreement with the chick electroporation data, Hes5(e1)-βG::nEGFP activity was spatially restricted, with high levels of expression seen only in intermediate regions of the spinal cord, where high levels of Hes5 were expressed (Fig 8H–8J). The ventral extent of Hes5(e1)-βG::nEGFP activity coincided with the dorsal border of the pMN domain with little overlap between Olig2 and GFP (Fig 8H). By contrast, in Olig2 mutant embryos, the expression of Hes5(e1)-βG::nEGFP extended ventrally to reach the dorsal boundary of Nkx2.2, a result that was not seen in control embryos (Fig 8K–8P). Together, these data provide evidence that Olig2 represses expression of Hes5 in the pMN domain, at least in part through direct interactions with the E-box site within Hes5(e1).
Here, we provide a detailed molecular description of somatic MN differentiation. Single cell transcriptomics defines distinct phases of differentiation and reveals the regulatory relationships that drive progression from NPs to postmitotic MNs. Experimental validation confirmed these predictions and demonstrated that Olig2 plays a pivotal role coordinating growth and patterning by integrating differentiation and fate determination signals (Fig 9).
Fig 9. Olig2 coordinates patterning and neuronal differentiation.
(A) Proposed model of the Olig2-controlled gene regulatory network. Olig2 not only acts as central organizer for dorsal-ventral patterning in the spinal cord but also controls the rate of MN differentiation through direct repression of Hes TFs. This leads to a higher levels of Ngn2 expression and, consequently, a higher rate of neuronal differentiation in the pMN domain, compared to adjacent progenitor domains. (B) Olig2 is a core component of the Shh-controlled gene regulatory network that patterns the ventral spinal cord [6,67]. (C) Olig2-mediated down-regulation of the Notch effectors Hes1/5 relieves repression of Ngn2 in the pMN domain. (D) Consolidated activities of Ngn2 and Olig2 cause differentiation of NPs to MNs. Olig2 promotes differentiation of MNs through repression of alternative IN cell fates. bHLH, basic helix-loop-helix; IN, interneuron; MN, motor neuron; NP, neural progenitor; pMN, MN progenitor; p2, V2 interneuron progenitor; p3, V3 interneuron progenitor; RA, retinoic acid; Shh, sonic hedgehog; TF, transcription factor; V2, V2 interneuron; V3, V3 interneuron.
The trajectory of NP to MN differentiation
Single cell mRNA sequencing is emerging as a powerful tool to reconstruct transcriptional changes in cells during tissue development [33–36]. Here, we use pseudo-temporal ordering of cells based on their expression profile to obtain a high-resolution map of the developmental trajectory of MN differentiation (Fig 2). Examining gene expression along this timeline highlighted the dynamics of signalling pathways and transcriptional networks as cells transit from proliferative progenitors to postmitotic neurons. This computationally reconstructed trajectory accurately recapitulated the known changes in gene expression associated with MN generation in vivo and identified features of the dynamics not previously evident. This provides evidence that, by exploiting the inherent heterogeneity and asynchrony of differentiating cells that confound population-based assays, scRNA-seq allows the inference of transcriptional dynamics during developmental cell state transitions with high resolution. Moreover, the data illustrate how scRNA-seq analysis of a defined developmental process in vitro accurately predicts gene regulatory interactions and transcriptional dynamics in vivo.
Examination of the timeline revealed periods of relatively stable gene expression. Punctuating these were transition phases with marked differences in gene expression profiles, which coincided with changes in the signalling status within the cells. This is consistent with the saltatory view of cell fate specification, in which differentiation proceeds through a series of metastable states separated by coordinated signal-driven changes in gene expression . Based on this observation, we used the global rate of change in gene expression in the pseudo-temporal orderings to develop a principled approach that objectively defines these phases. These distinct phases identified in MN differentiation corresponded to known cell types. Expression of Irx3 marked early, uncommitted NPs, normally located in the intermediate spinal cord. These progenitors transition to pMN cells in response to Shh and retinoid signaling, and this was identifiable by the up-regulation of Olig2 and down-regulation of Irx3. In addition, distinct phases in the acquisition of postmitotic MN identity could be identified with cells expressing markers such as Lhx3, Isl1/2, and Chat correctly positioned in the pseudo-temporal ordering.
The temporal ordering provided much greater resolution of the sequence of events leading to MN commitment than previously available. In particular, the transition from MN progenitor to MN was associated with a series of distinct and transient expression changes. This included the induction of well-known proneurogenic factors such as Ngn2, Neurod1, Neurod4, and Hes6 (Fig 2E and S3F Fig). Increased expression of Olig2 was also associated with this stage (Fig 2F). Consequently, the level of Olig2 expression distinguished two sequential stages in MN progenitors during their differentiation. In the earlier phase, initiated as Irx3 is down-regulated, pMN cells express low or moderate levels of Olig2. This is followed by the second phase, in which the levels of Olig2 substantially increase and Ngn2 becomes expressed at high levels (Fig 2C). In vivo analysis, together with the short-term lineage tracing afforded by the Olig2-mKate2 reporter, confirmed that Olig2 up-regulation coincided with the commitment to differentiate into postmitotic MNs. By contrast, in the earlier phase of pMN development, the lower levels of Olig2 appeared compatible with the transition of cells to Nkx2.2 and Fabp7 expressing p3 progenitors (S3C and S3D Fig). Together, these data provide new insight into the process of MN specification, identifying a series of distinct phases in NP differentiation to fate commitment and highlighting the changes in gene expression that characterise phase transitions.
Olig2 as a coordinator of neurogenesis
Previous studies have shown that both Olig2 and Ngn2 are required for the elaboration of MN identity and that Olig2 activity induces Ngn2 expression [13,17,18,21]. Consistent with these observations, progenitors in the Olig2 expression domain differentiate at a much higher rate than cells in other progenitor domains of the neural tube [11,13]. Our results suggest a mechanism for this enhanced rate of neuronal differentiation. The canonical Notch effectors Hes1 and Hes5 act to suppress neurogenesis by inhibiting the expression of neurogenic bHLH proteins, thus maintaining NPs in an undifferentiated state [24,25]. Olig2 activity represses Hes1 and Hes5, thereby allowing expression of the proneural gene Ngn2 and downstream effectors such as Neurod4 (Fig 6 and S7 Fig). The ability of Olig2 to repress Hes5 appears to be direct, as Olig2 binds to a conserved regulatory element within the Hes5 locus that restricts gene expression from MN progenitors (Fig 7 and Fig 8). Similarly, Olig2 binding sites are found in putative regulatory elements associated with the Hes1 gene, raising the possibility that this regulatory interaction is also direct. Consistent with a role for Olig2 in promoting neuronal differentiation, the levels of Olig2 transcript and protein peak at the onset of neurogenesis, concomitant with the induction of Ngn2 in vivo and in vitro (Fig 2 and Fig 3). These findings therefore establish a mechanism by which neural patterning and neurogenesis intersect. In this view, by modulating the Notch pathway, the Shh and retinoid-dependent induction of Olig2 not only specifies MN identity but also determines the rate at which these progenitors differentiate, thus imposing the distinctive kinetics of MN production (Fig 9).
This model is surprising, as previous studies suggested antagonistic activities for Olig2 and Ngn2 during the induction of neuronal target genes . Both Olig2 and Ngn2 have been shown to heterodimerize with E47 and bind to E-box elements, but with opposing activities [28,68]. In addition, similar to Id proteins, Olig2 proteins could potentially sequester E proteins (E12 and E47) from forming heterodimeric Ngn2/E-protein complexes that activate transcription [28,69]. The sequential expression of Olig2 and Ngn2 has been proposed as a potential mechanism to reconcile the inhibitory activity of Olig2 on neurogenesis with the high rate of neurogenesis in the pMN domain . However, our results suggest that this mechanism is unlikely to apply to the differentiation of MNs in the spinal cord for several reasons. The higher temporal resolution provided by the pseudo-temporal ordering indicated that primary Ngn2 target genes such as Dll1 and Neurod4 are induced when the rate of Olig2 expression is maximal in cells (Fig 2E and S3F Fig). Furthermore, Olig2 protein perdures longer in differentiating MNs than Ngn2, resulting in significant coexpression of Olig2 and early MN markers such as Lhx3 and Mnx1 in Ngn2-negative cells (Fig 3H and S5B and S5C Fig). Hence, instead of sequential expression of these proteins, these results suggest that Ngn2 is capable of mediating neurogenesis despite the presence of high Olig2 levels.
A potential solution to this puzzle may be that the activities of both Olig2 and Ngn2 are regulated by phosphorylation. Olig2 phosphorylation at specific Ser/Thr residues regulates its choice of dimerization partner, intracellular localization, and DNA binding preference for open and closed chromatin [70–73]. Indeed, homodimeric complexes of Olig2 appear to mediate Hes5 repression (Fig 7C and 7D). Furthermore, the cyclin-dependent kinases CDK1/2 have been proposed to phosphorylate Olig2 at Ser14, priming Olig2 for further phosphorylation at multiple Ser residues . This phosphorylation appears to regulate the preference of Olig2 for open or closed chromatin and, thus, strongly influences its biological activity . The declining levels of CDK1/2 during neurogenesis may similarly affect Olig2 activity to overcome its inhibitory role on neuronal gene expression. Similarly, phosphorylation of Ngn2 affects its stability and interaction with Lim-homeodomain TF complexes and E-proteins [75–78]. Thus, additional posttranslational events extend the regulatory interactions between both proteins beyond stochiometric interactions through protein–protein binding and competition for DNA binding sites. Notably, some of the relevant phosphorylations are performed by protein kinase A (PKA) and glycogen-synthase kinase 3 (GSK3), kinases linked to the activity of the Shh pathway, which appears to peak at the initiation of MN differentiation [6,11,72,75]. Connecting the activity of these neurogenic TFs to the activity of the Shh pathway would allow a tight coupling between MN generation and overall developmental dynamics dictated by the dynamics of morphogen signalling.
The pseudo-temporal ordering indicated that although Olig2 levels peaked at the onset of MN differentiation, expression then decreased rapidly, prior to MNs reaching the next metastable state along the differentiation trajectory, which is characterized by the induction of mature MN markers such as Isl2 and Chat (Fig 2). This suggests that Olig2 may need to be down-regulated to allow progression of MN differentiation. Consistent with this, overexpression of Olig2 has been shown to inhibit the generation of MNs and to directly repress genes associated with MN identity, such as Mnx1 [28,68]. Furthermore, the addition of Olig2 to canonical reprogramming factors decreases the efficiency of conversion from fibroblasts to spinal MNs . Thus, Olig2 up-regulation can promote MN generation by initiating differentiation, but its down-regulation is needed to complete the switch from progenitor to postmitotic neuron. These dynamics of Olig2 expression may help impose directionality to differentiation and ensure that the correct temporal sequence of gene expression occurs as MNs mature.
Oscillation of bHLH TFs in the spinal cord
The maintenance of NPs in the brain has been ascribed to the oscillatory expression of Hes and proneural bHLH TFs [65,79]. The Hes proteins are proposed to generate the oscillations by negatively regulating their own expression as well as Ngn2 and Ascl1 [25,80,81]. This phenomenon results in Hes1 and proneural bHLH TFs exhibiting reciprocal expression phases at an equivalent frequency [65,79]. Oscillations in the levels of the Notch ligand Dll1 have been reported in spinal cord progenitors . In cortical progenitors, fluctuations in Olig2 levels have also been documented, but these oscillations occur at a significantly slower frequency  and may therefore be regulated by a different mechanism.
Although we did not specifically investigate the occurrence of bHLH oscillations in the spinal cord in vitro or in vivo, our results may shed light on this question. The Hes5(e1) element can be bound by both Olig2-Olig2 repression dimers and Ngn2-E12 activation complexes (Fig 7B and 7C and S8A Fig). It is notable that mutation of the E-box in this element reduced the overall level of reporter activity in the spinal cord (S8B and S8C Fig) at the same time as it disrupted its spatial restriction from the pMN (Fig 8). These data are consistent with a model in which positive activators, such as Ngn2 or other E-box binding factors, could interact with Hes5(e1) to directly elevate Hes5 expression, which would in turn serve to repress Ngn2 expression, thereby contributing to alternating phases of Hes5 and Ngn2 expression. In this regard, Olig2 binding and repressing Hes5 through this element would interrupt the oscillator, allowing Ngn2 expression to reach its maximal levels and neuronal differentiation to commence. Thus, by inserting itself into the Notch-regulated neural differentiation program, Olig2 shuts down Notch activity, ensuring MN development proceeds in a spatially and temporally controlled manner. This reconciles stochastic and oscillatory models of neuronal differentiation with the spatially predetermined pattern of neuron production observed in the spinal cord.
To examine Olig2 expression, we used the relatively long half-life fluorescent protein mKate2, introduced into the Olig2 genomic locus. Quantification of the levels of mKate2 and Olig2 revealed a striking correlation between both proteins in NPs (S6D–S6F Fig). This argues against Olig2 oscillations in these cells, as oscillatory behavior would be expected to decrease the correlation between both proteins. Although further investigation is necessary, the data are consistent with out-of-phase oscillations between Ngn2 and Hes5, while Olig2 levels steadily increase in MN progenitors over time. Understanding these relationships will provide insight into the transition from MN progenitor to differentiation.
The Notch pathway regulates Olig2 expression and Shh signaling
Besides promoting neurogenesis, inhibition of Notch effectors also appears to be important for dorsal-ventral patterning of the neural tube and the consolidation of pMN identity. Patterning of the ventral neural tube is mediated by a gene regulatory network that interprets both levels and duration of Shh signalling [6,42,67,83]. Previous studies have suggested that Notch signalling influences patterning of the ventral spinal cord by promoting the activity of the Shh pathway [61,84]. Consistent with this, overexpression of HAIRY2, the chick homologue of murine Hes1, causes a down-regulation of Olig2 and induction of Nkx2.2 in the pMN domain . Similarly, sustained activation or inhibition of the Notch pathway causes, respectively, a ventral expansion or recession in p3 progenitors, located ventral to the pMN . Here, we show that besides modulating Shh activity, Notch signalling can also regulate expression levels of Olig2. Conversely, Olig2 represses the canonical Notch effector Hes5 and could thereby negatively regulate the levels of Shh signalling in the pMN domain. Thus, Olig2 may consolidate pMN identity not only by direct repression of other progenitor markers but also indirectly, by modulating levels of Shh signalling through its effect on Notch pathway.
Taken together, our data reveal a tight coupling between the gene regulatory networks that control patterning and differentiation in the ventral spinal cord. This highlights the pivotal role of Olig2 in this process, which not only acts as central organizer of dorsal-ventral patterning in the spinal cord but also as developmental pacemaker for MN formation. The Olig2-mediated repression of Notch pathway targets provides a molecular mechanism for the much higher rate of neurogenesis observed in the pMN domain, compared to the rest of the spinal cord, and thereby explains the spatial and temporal patterns of neurogenesis observed in the neural tube. These findings raise the question of whether similar mechanisms also apply in other progenitor domains in the neural tube.
Materials and methods
Animal experiments in the Briscoe lab were performed under UK Home Office project licenses (PPL80/2528 and PD415DD17) within the conditions of the Animal (Scientific Procedures) Act 1986. Animals were only handled by personal license holders. Olig2Cre and Ngn2KIGFP knock-in/knock-out mice were used as previously described [56,64] and interbred to create Olig2 or Ngn2 mutant embryos. All mice in the Novitch lab were maintained and tissue collected in accordance with guidelines set forth by the UCLA Institutional Animal Care and Use Committee. Fertilized chicken eggs were acquired from AA McIntyre Poultry and Fertile Eggs, incubated, and electroporated as previously described .
Differentiation of NPs from mouse ESCs
NPs were differentiated as described previously . In brief, HM1 (Thermo Scientific), DVI2, and Olig2::T2A-mKate2 ESCs were maintained in ES cell medium with 1,000 U/ml LIF on mitotically inactivated mouse embryonic fibroblasts (feeder cells). DVI2 cells were generated by integrating an 8xGBS-H2B::Venus Shh pathway reporter into the HPRT locus of HM1 cells and used for all ESC experiments except 4-color stainings in S5A and S5B Fig, which rely on HM1 cells, and experiments in Fig 4 and S6 Fig, which were conducted using the Olig2::T2A-mKate2 reporter cell line.
For differentiation, cells were dissociated in 0.05% Trypsin (Gibco) and replated onto tissue culture plates for 25 minutes to remove feeder cells. Cells staying in the supernatant were spun down and resuspended in N2B27 medium at a concentration of 106 cells/ml. Forty-five thousand cells were plated onto 35-mm CellBind dishes (Corning) precoated with 0.1% Gelatine solution in 1.5 ml N2B27 + 10 ng/ml bFGF. At D2, the medium was replaced with N2B27 + 10 ng/ml bFGF + 5 μM CHIR99021 (Axon). At D3 and every 24 hours afterwards, the medium was replaced with N2B27 + 100 nM RA (Sigma) + 500 nM SAG (Calbiochem). For Notch inhibition, differentiations were treated at day 5 with N2B27 + 100 nM RA + 500 nM SAG + 10 ng/μl DBZ (Tocris Biosciences) for 24 hours. Cells were washed with N2B27 medium at later medium changes, when many dead cells were detected in the dish.
mRNA was extracted using RNeasy Mini Kit (Qiagen) according to the manufacturer’s instructions. 1.5–2 μg of RNA was used for reverse transcription using Super-Script III First-Strand Synthesis kit (Invitrogen) with random hexamers. Platinum SYBR Green qPCR mix (Invitrogen) was used for amplification on a 7900HT Fast Real Time PCR machine (Applied Biosystems). Expression values were normalized against β-actin. Three independent repeats of each RT-qPCR time course were performed and three independent samples at each time point of each repeat were analyzed. For a complete list of used primers, see S3 Table. qPCR data presented in Fig 1, S1 Fig, and S6 Fig show one representative repeat and show mean ± standard deviation. The heat map in S2B Fig was plotted using Graphpad Prism 7.
Cells were washed using N2B27 medium and PBS (Gibco) and then fixed in 4% paraformaldehyde in PBS at 4°C for 20 minutes. After fixation, cells were washed twice with PBS and stored in a refrigerator until stainings were perfomed. For staining, cells were washed three times in PBS containing 0.1% Triton X-100 (PBS-T). Primary and secondary antibodies were diluted in PBS-T + 1% BSA. Cells were incubated with primary antibodies overnight at 4°C, then washed three times for 5–10 minutes in PBS-T, incubated with secondary antibodies for 1 hour at room temperature, and washed again three times in PBS-T. Stainings were mounted using ProLong Gold Antifade reagent (Life Technologies). Mouse and chicken spinal cord tissues were fixed with 4% paraformaldehyde, cryoprotected in 30% sucrose, sectioned, and processed for immunohistochemistry or in situ hybridization, as previously described [85,86].
Antibodies against a peptide in the C-terminal portion of mouse Hes5 (APAKEPPAPGAAPQPARSSAK, aa 127–147) were raised in rabbits and guinea pigs (Covance). The rabbit serum was affinity purified and used at 1:8,000 and the crude guinea pig serum at 1:16,000. Additional primary antibodies were used as follows: goat anti-β-galactosidase (Biogenesis 4600–1409 1:2,000), mouse anti-Cre (Covance MMS-106P, 1:2,000), rabbit anti-Dbx1 (kind gift of Susan Morten and Thomas Jessell, 1:8,000), rabbit anti-Fabp7 (Abcam ab32423, 1:2,000 or Chemical AB9558, 1:2,000), rat anti-FLAG (Stratagene 200474, 1:1,500), chicken anti-GFP (Abcam ab13970, 1:20,000), sheep anti-GFP (AbD Serotec 4745–1051, 1:800), rabbit anti-Hes1 (, 1:1,000), mouse anti-Hoxc6 (Santa Cruz Biotechnology sc-376330, 1:250), mouse anti-Hb9/Mnx1 (DSHB, 1:40), mouse anti-Isl1/2 (DSHB, 1:100), goat anti-Isl1 (R&D AF1837, 1:1,000), rabbit anti-Lhx3 (Abcam ab14555, 1:500), mouse anti-NeuN (Rbfox3, Chemicon/Millipore MAB377, 1:1,000), rat anti-chick Neurod4 (NeuroM ), goat anti-Ngn2 (Santa Cruz Biotechnology sc-19233, 1:500), mouse anti-Ngn2 (5C6, , 1:50), guinea pig anti-chick Ngn2 (, 1:32,000), mouse anti-Nkx2.2 (DSHB, 1:25), mouse anti-Nkx6.1 (DSHB, 1:100), rabbit anti-Olig2 (Millipore AB9610, 1:1,000), guinea pig anti-mouse Olig2 ( 1:20,000), guinea pig anti-chick Olig2 (, 1:8,000), rabbit anti-Pax6 (Millipore AB2237, 1:1,000), mouse anti-Pax6 (DSHB, 1:25), goat anti-Sox1 (R&D AF3369, 1:500), goat anti-Sox2 (Santa Cruz Biotechnology sc-17320, 1:2,000), rabbit anti-Sox2 (, 1:2500), rabbit anti-TagRFP (Evrogen AB233, 1:1,000), rabbit anti-Tubb3 (Covance PRB-435P, 1:2,000), mouse anti-Tubb3 (Covance MMS-435P, 1:1,000), rabbit anti-Zbtb18 (Proteintech 12714-1-AP, 1:1,000).
Secondary antibodies used throughout this study were raised in donkey. Alexa488, Alexa568, Cy3, and Dylight 647-conjugated antibodies (Life Technologies or Jackson Immunoresearch) were diluted 1:1,000 and Alexa647 conjugated antibodies (Life Technologies) 1:500. Cy5-conjugated antibodies (Jackson Immunoresearch) were diluted 1:700 and CF405M donkey anti-guinea pig secondary antibody (Sigma) 1:250.
Image acquisition and analysis
Immunofluorescent images of ESC-derived NPs were acquired using a Zeiss Imager.Z2 microscope equipped with an Apotome.2 structured illumination module and a 20× magnification lens (NA = 0.75). Five phase images were acquired for structured illumination. For each image, z-stacks composed of 12 sections separated by 1 μm were acquired. Maximum intensity projection was performed in Fiji.
Cryosections were documented using a Leica SP5 confocal microscope equipped with a 40× oil objective, or Zeiss LSM5, LSM700, or LSM800 confocal microscopes and Zeiss Apotome imaging systems equipped with 10×, 20×, and 40× oil objectives. For nuclear staining intensity measurements, 3–4 individual sections separated by 1 μm were analyzed. Nuclei segmentation and intensity measurement were performed in CellProfiler. Data were normalized and plotted using R. Other images were processed and manually quantified using Fiji and Adobe Photoshop imaging software.
Single cell sequencing
NPs were dissociated using 0.05% Trypsin (Gibco) spun down in ES-medium, resuspended, washed, and spun down in 10 ml PBS (Gibco). Afterwards, cells were resuspended in 1 ml N2B27 and filtered into a FACS tube (Falcon). The Fluidigm C1 platform was used to capture individual cells using 96 small or medium IFC chips. Cells were diluted in the range of 250,000–400,000 cells per ml for chip loading. Capturing efficiency was evaluated by manually inspecting each capture site on the chip using the automated NanoEntek JuLi cell imager. Only capture sites containing single cells were processed for library preparation and sequencing. Single cell full-length cDNA was generated using the Clontech SMARTer Ultra Low RNA kit on the C1 chip using manufacturer-provided protocol. ArrayControl RNA Spikes (AM1780) were added to the cell lysis mix, as recommended in the Fluidigm protocol. Libraries were prepared using the Illumina Nextera XT DNA Sample Preparation kit, according to a protocol supplied by Fluidigm, and sequenced on Illumina Hiseq 2500 or 4000 using 50- or 75-bp paired-end runs.
Generation of Olig2::T2A-mKate2 ESC line by CRISPR
pNTKV-T2A-3xNLS-FLAG-mKate2 was generated by cloning a T2A-3XNLS-FLAG-mKate2 cassette into pNTKV using HpaI and HindIII restriction sites. To integrate the T2A-3xNLS-FLAG-mKate2 cassette at the 3′ end of the Olig2 open reading frame, a donor vector comprising app. 2.8 kb upstream and 5 kb downstream of the stop codon was constructed. For CRISPR/Cas9-mediated homologous recombination, a short guide RNA (sgRNA) sequence (CGGCCAGCGGGGGTGCGTCC) was cloned into pX459 (Addgene), according to .
For electroporation, DVI2 ESCs were grown in 2i medium + LIF. 4 μg of both plasmids were electroporated into 4×106 cells using Nucleofector II (Amaxa) and mouse ESC Nucleofector kit (Lonza). Afterwards, cells were replated onto 10-cm CellBind plates (Corning) and maintained in 2i medium + LIF. For selection, cells were first treated with 1.5 μg/ml Puromycin (Sigma) for two days and afterwards with 50 μg/ml Geniticin (Gibco) until colonies were clearly visible. Individual colonies were picked using a 2-μl pipette, dissociated in 0.25% Trypsin (Gibco), and replated onto feeder cells in ES-medium + 1,000 U/ml LIF in a 96-well plate. Correct integration of the T2A-3xNLS-FLAG-mKate2 transgene was verified using long-range PCRs.
Cells were lysed in RIPA buffer supplemented with protease inhibitors. 10 μg of total protein was loaded per lane. Antibodies used were rabbit anti-Olig2 (Millipore AB9610, 1:3,000) and mouse anti-β-tubulin (Sigma T4026, 1:2,000). Secondary antibodies were donkey anti-mouse IRDye 800CW and donkey anti-rabbit IRDye 680RD (both Licor). Blots were scanned using an Odyssey Scanner (Licor).
To quantify the number of Olig2::T2A-mKate2 positive cells, NPs were trypsinized in 0.05% Trypsin (Gibco) and spun down in ES medium. Cells were resuspended in PBS (Gibco) supplemented with the live-cell staining dye Calcein Violet (Life Technologies), according to the manufacturer’s instructions. Flow analysis was performed using a Becton Dickinson LSRII flow cytometer. HM1 or DVI2 cells were differentiated in parallel and analyzed using the same settings to estimate the number of mKate2 positive cells in Fig 4L and S6H Fig. The threshold for counting cells as mKate2 positive was set to the intensity for which 0.5% of cells without the mKate2 transgene were counted as positive. Thirty thousand events were recorded for each replicate. To count cells as mKate2HIGH cells in Fig 4L and S6G Fig, a threshold was set to the shoulder visible in the histograms of the control conditions depicted in Fig 4L. This shoulder typically appears at day 6 and is not present in day 5 cells. Based on the flow cytometry data in S6C Fig and the correlation between mKate2 fluorescence and Isl1/2 expression shown in Fig 4H, we infer that most mKate2HIGH cells are MNs. Cells were counted as mKate2HIGH if their mKate2 intensity exceeded this threshold.
For analysis of Tubb3 stainings by flow cytometry (S6 Fig), cells were fixed for 20 minutes in 2% PFA on ice, washed three times in PBS-T, and incubated with Alexa647-conjugated anti-Tubb3 antibody (BD Pharmingen, 1:10) in PBS-T + 1% BSA for one hour at room temperature. Cells were afterwards washed 3× in PBS-T and endogenous mKate2 fluorescence, and Tubb3 staining in 10,000 cells was quantified by flow cytometry. Data were analyzed using FlowJo.
Electrophoretic mobility shift assays
pCS2+ plasmid expression vectors for Olig2, E12, Ngn2, and Id1  were transcribed and translated in vitro using the Promega TNT Coupled Wheat Germ Extract System. Programmed extracts were mixed as indicated in a buffer containing 100 mM Hepes pH 7.6, 25 mM KCl, 1.5 mM MgCl2, 0.2 mM EDTA, 2.5% glycerol, and 100 ng poly dIdC and incubated for 15 minutes at room temperature. 32P-dCTP-labeled probes were generated by Klenow end labeling of double-stranded oligonucleotides containing the native mouse Hes5(e1) sequence (forward 5′- ggccgCTCCCAAAAGACCATCTGGCTCCGTGTTATAA-3′; reverse 5′- actagTTATAACACGGAGCCAGATGGTCTTTTGGGAG-3′) or an E-box mutated version (forward 5′ ggccgCTCCCAAAAGAggATCccGCTCCGTGTTATAA-3′; reverse 5′-actagTTATAACACGGAGCggGATccTCTTTTGGGAG-3′). The E-box sequence is underlined. Lowercase indicates substitutions and added flanking sequences. Samples were incubated with labeled probes for 15 minutes before resolving on a 4.5% polyacrylamide gel and subsequent autoradiography. Probe competition was achieved by incorporating unlabeled oligonucleotide probes in the binding reaction mix.
Hes5(e1) transgenic assays
Mouse and chick Hes5(e1) genomic DNA fragments were amplified by PCR using the following primers: mouse, forward 5′-gaggcggccgcCGGTTCCCACACTTTGGT-3′; reverse 5′-gagactagtCACAGTCCCAAGCTGCTTAAA-3′, and chick, forward 5′-gaggcggccgcTGCGTTTCCCATACTTTTCC-3′; reverse 5′-gagactagtTCTGGCCTTGAAGCTAGGAG-3′. Lowercase indicates added flanking sequences with restriction enzyme sites (NotI and SpeI) underlined. PCR products were digested and cloned into a reporter construct containing the β-globin basal promoter, a nuclear EGFP coding sequence, and a bovine growth hormone polyadenylation sequence (βG::nGFP ). Mutations of the Hes5(e1) E-box were generated through splicing by overlap extension PCR. Chick embyos were co-electroporated with chick Hes5(e1) constructs along with a plasmid vector, producing nuclear-tagged β-galactosidase under the control of the cytomegalovirus enhancer and β-actin promoter. Embryos were collected, fixed, and cryosectioned. Sections were stained using antibodies to β-galactosidase and Olig2 and fluorescent secondary antibodies. GFP levels were measured based on its native fluorescence. Images were collected and positions of cells expressing each marker rendered using the spots function of the Bitplane Imaris 8.4 imaging suite. Calculated positions were exported and processed using Microsoft Excel and Graphpad Prism 7 software.
Transgenic mice expressing the mouse Hes5(e1)-βG::nGFP reporter were generated with the assistance of the University of Michigan Transgenic Animal Model Core by microinjection of purified plasmid DNA into fertilized eggs obtained by mating (C57BL/6 X SJL)F1 female mice with (C57BL/6 X SJL)F1 male mice and subsequent transfer to pseudo-pregnant recipients. Analysis was conducted on both embryos collected from the primary reporter injections and offspring collected from matings of a transgenic line that passed through the germ line.
Electroporation and in situ hybridization
Chick embryos were electroporated at HH stages 11–13 with RCAS-myc-tagged chick Olig2 and Olig2-bHLH-Engrailed plasmid constructs  and collected at HH stages 21–23. Spinal cord sections were hybridized with digoxigenin-UTP labeled RNA probes generated from plasmid DNA or templates generated by PCR. 3′ UTR sequences for chick HES5-1 and HES5-3 were amplified from spinal cord cDNA with the following primers: HES5-1, forward 5′-GCGGAATTCAGGGAAGCTCTCACTTAGTGAAC-3′ and reverse 5′-GCGCTCGAGATACCCTCCTGCTGAAGACATTTGC-3′; HES5-3, forward 5′-GCGGAATTCGCCAAGAGCACGCTCACCATCACCT-3′ and reverse 5′-GCGCTCGAGCTACACAGCTTGAGTTATGGTTTAG-3′ and directionally cloned into the pBluescript. Underlined sequences indicate restriction enzyme sites incorporated into the primers. Chick HAIRY1 and HES5-2 3′ UTR riboprobes were generated as previously described .
ESC-derived MN progenitors were dounce homogenized and sonicated in 3 ml lysis buffer (1% SDS, 50 mM Tris pH 8.0, 20 mM EDTA, 1 mM PMSF, and 1X Complete protease inhibitor cocktail [Roche]). 150 μg of lysate DNA were used per immunoprecipitation reaction and mixed with 3–5 μl of either rabbit anti-Olig2 antibodies (Millipore AB9610), rabbit anti-Ngn2 serum (generous gift of Dr. Soo-Kyung Lee, Oregon Health Sciences University), normal rabbit sera, or purified rabbit IgG. Antibody-chromatin complexes were collected using Dynabeads protein A (Invitrogen), washed and eluted DNA used for RT-qPCR in triplicate using the following primer pairs: Hes5(e1) forward 5′-CTGCTTCTGAATGAATGAGGGCGG-3′ and reverse 5′-AGCAGACGAGCCCTTTATTGCTCT-3′; Hes5(e2), a nonconserved element 3′ to the Hes5 coding exons that contains an E-box, forward 5′-AGATGGCTCAGCGGTTAAGAG-3′ and reverse 5′- CCATGTGGTTGCTGGGATTTG-3′. Fold enrichment for each region was calculated as compared to normal rabbit serum or purified IgG.
S1 Fig. Characterization of Hox gene expression in NPs and MNs.
(A) Expression of the somatic MN marker Mnx1 in MNs at day 6. (B) RT-qPCR analysis of Hox genes expression levels from day 3 to day 7. Underlying data are provided in S1 Data. (C) Hoxc6 expression in MNs characterized by Isl1 and Tubb3 expression from day 6 to day 8. Scale bars = 40 μm. MN, motor neuron; NP, neural progenitor; RT-qPCR, real time-quantitative polymerase chain reaction; Tubb3, neuronal class III beta-tubulin.
S2 Fig. Establishment of different NP identities by different levels of Shh pathway activation.
(A) Schematic of the embryonic spinal cord. Expression domains of TFs defining NP domains are indicated. (B) RT-qPCR analysis of day 6 differentiations treated with 0–1,000 nM SAG after day 3. (C) Expression of Dbx1, Nkx6.1, and Isl1 in day 6 differentiations treated with the indicated concentrations of SAG. (D) Expression of Pax6, Olig2, and Nkx2.2 in day 6 differentiations treated with the indicated concentrations of SAG. Scale bars = 40 μm. NP, neural progenitor; Shh, sonic hedgehog; TF, transcription factor; RT-qPCR, real time-quantitative polymerase chain reaction; SAG, Smoothened/Shh signalling agonist.
S3 Fig. Identification of gene modules and cell states by hierarchical clustering of single cell sequencing data.
(A) Identification of cell states by hierarchical clustering from 202 cells based on 10 identified gene modules. Genes characteristic for the individual modules are indicated. The boxed region corresponds to the heat map in Fig 2A. (B) Quantifications of read counts per cell (top) and number of expressed genes per cell (bottom) for neural cell states identified by hierarchical clustering. Colors of the graphs match cell states in Fig 2A and S3A Fig. (C) Cell state graphs color coded for the expression levels of Olig2, Nkx2.2, and Fabp7. (D) Analysis of Olig2, Nkx2.2, and Fabp7 expression in differentiations at day 6 (top row) and e10.5 embryonic spinal cords (bottom row) confirms higher Fabp7 expression levels in p3 progenitors. (E) The transition phase from early Irx3 NPs to Olig2 NPs correlates with the induction of Shh target genes Ptch2, Gli1, Hhip. (F) Inhibition of Notch signalling, revealed by decreasing expression levels of Hes1/5 and expression of the Notch ligand Dll1 and the pathway inhibitors Hes6 and Numbl, identifies the cell state transition from NPs to MNs. Scale bars = 25 μm (D, top row) and 50 μm (D, bottom row). e, embryonic day; NP, neural progenitor; Shh, sonic hedgehog.
S4 Fig. Robustness analysis of gene expression dynamics by bootstrapping.
Solid lines indicate original gene levels shown in Figs 2E, 2F and S3E and S3F. At each pseudo–time point, box plots indicate median, first, and third quartiles of the gene level distribution obtained from 1,000 bootstrapped datasets. Whiskers indicate the largest and smallest values no further than 1.5 times the interquartile range, taken from the hinge. The associated correlation coefficient, r, is the average of the Spearman correlation coefficients over all pairs of boostrap replicates.
S5 Fig. Validation of predictions from the pseudo-temporal ordering.
(A) Sequential expression of Olig2, Ngn2, Lhx3, and Isl1 during MN differentiation, revealed by immunofluorescent staining for these markers at day 6 of differentiation. (B,C) Quantification of levels of Olig2, Isl1, Ngn2 (color code in B), and Lhx3 (color code in C) reveals a clear differentiation path from Olig2HIGH cells to MNs and sequential induction of Ngn2 and Lhx3 during this process (n = 2,236 nuclei). Underlying data are provided in S1 Data. (D) Staining for Isl1, Lhx3, and Tubb3 reveals high levels of Tubb3 expression in Isl1-positive but not Lhx3-positive MNs at day 6 of differentiation. This is consistent with the earlier MN stage of Lhx3 MNs. (E, F) Most positive and negative Spearman-correlated transcription factors for Olig2 (E) and Ngn2 (F) reveal Zbtb18 (green in E–G) as a novel gene involved in MN formation. Underlying data are provided in S1 Data. (G–G″) Immunofluorescent staining for Olig2 (G), Zbtb18 (G′), and Ngn2 (G″) at day 6 of differentiation. (H, I) Quantification of levels of Olig2, Ngn2 (H), and Zbtb18 (I) in individual nuclei reveals a good correlation between these markers (n = 1,431 nuclei). Underlying data are provided in S1 Data. (J–L) Analysis of Olig2, Ngn2, and Zbtb18 expression in neural tubes at e9.5 (J), e10.5 (K), and e11.5 (L). Note that Ngn2 and Zbtb18 are expressed in cells with high levels of Olig2 at e9.5 and e10.5, but not at e11.5 (left insets in K–M). In addition, Zbtb18 and Ngn2 are coexpressed in nuclei at the edge of the progenitor domain in dorsal areas of the neural tube at e10.5 (L) and e11.5 (M) (right insets). Scale bars = 25 μm in (A,D), 10 μm in (G) and insets in J–L, 50 μm in J–L. e, embryonic day; MN, motor neuron; Tubb3, neuronal class III beta-tubulin.
Characterization of the Olig2-mKate2 reporter cell line by flow cytometry and upon Notch inhibition (A–C) Quantification of mKate2 and Tubb3 fluorescence intensity by flow cytometry at day 4 to day 6 of differentiation. Note that high levels of Tubb3 are predominantly detected in mKate2HIGH cells at day 6 (C). (D–F) Correlation between Olig2 and mKate2 levels in individual nuclei quantified from images in Fig 4C–4F. Plots are color coded for levels of Isl1/2 (D), Sox1 (E), and Nkx2.2 (F). (G) Quantification of the fold change in mKate2HIGH cells (see Fig 4L) upon 24 hours Notch inhibition for five experimental repeats by flow cytometry. Each repeat consists of the measurement of three independent dishes for control and Notch inhibition from the same differentiation. Underlying data are provided in S1 Data. *** p < 0.001, unpaired t test. (H) Fold change of mKate2-positive cells (see Fig 4L) upon Notch inhibition (grey) relative to untreated control differentiations (black). Notch inhibition does not cause an overall change in the number of mKate2 positive cells. Underlying data are provided in S1 Data. * p < 0.05; ** p < 0.01, unpaired t test. (I) Quantification of mKate2 and Tubb3 fluorescence intensity by flow cytometry upon 24 hours Notch inhibition. Note that most mKate2HIGH cells differentiated into MNs (compare to S6C Fig). (J–L) RT-qPCR quantification of expression levels of progenitor markers Sox1, Hes5, Hes1, and Pax6 (J); neurogenesis markers Olig2, Ngn2, Zbtb18, and Pou3f2 (K); and MN markers Tubb3 and Isl1 (L) after 0, 12, and 24 hours of Notch inhibition (grey) and in untreated controls (black). Note that Olig2 expression increases in contrast to other progenitor markers after 12 hours of Notch inhibition. Underlying data are provided in S1 Data. MN, motor neuron; RT-qPCR, real time-quantitative polymerase chain reaction; Tubb3, neuronal class III beta-tubulin.
S7 Fig. Olig2, not Ngn2, is required for the repression of Hes5, which is necessary for the formation of MNs.
(A–D″) Staining of wild-type (A–B″) and Ngn2KIGFP mutant (C–D″) embryonic spinal cords for Olig2, Ngn2, GFP, and Hes5 at e9.5 (A,C) and e10.25 (B,D). (E–F″) GFP expression (E,F) is still increased and Hes5 expression (E″,F″) is still reduced in the pMN domain in Ngn2KIGFP mutant spinal cords (same sections as C–D″). The yellow dotted line indicates the dorsal boundary of the pMN domain. (G–J) Ectopic expression of cHES5-2 does not affect levels of the progenitor markers SOX2 (G), NKX6.1 (H), PAX6 (I), and OLIG2 (J). (K-P) Ectopic expression of cHES5-2 (N–P) leads to a reduction of NGN2 (K,N), NEUROD4 (L,O), and NEUN (M,P). (K–M) These show control electroporations with a nuclear EGFP (nEGFP) expression construct. (Q) Quantification of the effect of ectopic cHES5-2 expression on expression levels of NGN2, NEUROD4, and NEUN relative to nEGFP controls. Underlying data are provided in S1 Data. **** p < 0.0001; *** p < 0.0005, Mann-Whitney test; Scale bars = 20 μm (A–A″,C–C″,G–P), 50 μm (B–B″,D–D″). cHES5-2, chick HES5-2 gene; e, embryonic day; EGFP, enhanced green fluorescent protein; GFP, green fluorescent protein; MN, motor neuron; nEGFP, nuclear EGFP; pMN, MN progenitor.
S8 Fig. Olig2 and Ngn2 bind to Hes5(e1) in vivo.
(A) Both Olig2 and Ngn2 antibodies precipitate the Hes5(e1) genomic element from ESC-derived MN progenitors, but not Hes5(e2), an unrelated genomic element 3′ to the Hes5 coding exons that also contains an E-box. Fold enrichment relative to normal rabbit sera or purified IgG is displayed. **** p < 0.0001; ** p < 0.01, Mann-Whitney test. Underlying data may be found in S1 Data. (B,C) Comparison between Hes5(e1)-βG::nEGFP (top row) and Hes5(e1ΔE)-βG::nEGFP (bottom row) reporter activities. Sections were imaged using identical settings for each pair. The overall activity of the Hes5(e1ΔE)-βG::nEGFP reporter is lower than that of the Hes5(e1)-βG::nEGFP reporter. Scale bars = 50 μm. E-box, bHLH protein binding site; ESC, embryonic stem cell.