Olig2 and Hes regulatory dynamics during motor neuron differentiation revealed by single cell transcriptomics

During tissue development, multipotent progenitors differentiate into specific cell types in characteristic spatial and temporal patterns. We addressed the mechanism linking progenitor identity and differentiation rate in the neural tube, where motor neuron (MN) progenitors differentiate more rapidly than other progenitors. Using single cell transcriptomics, we defined the transcriptional changes associated with the transition of neural progenitors into MNs. Reconstruction of gene expression dynamics from these data indicate a pivotal role for the MN determinant Olig2 just prior to MN differentiation. Olig2 represses expression of the Notch signaling pathway effectors Hes1 and Hes5. Olig2 repression of Hes5 appears to be direct, via a conserved regulatory element within the Hes5 locus that restricts expression from MN progenitors. These findings reveal a tight coupling between the regulatory networks that control patterning and neuronal differentiation and demonstrate how Olig2 acts as the developmental pacemaker coordinating the spatial and temporal pattern of MN generation.


Introduction
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][4][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][9][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 [11]. 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 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 [13], and the absence of Olig2 results in a characteristically slower tempo of neuronal differentiation [14]. 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][20][21][22][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][25][26].
Although these studies suggested that Olig2 promotes motor neurogenesis, subsequent studies ascribed antineurogenic and pro-proliferative functions to Olig2 (reviewed in [27]). These conclusions were based on the Olig2-mediated repression of Mnx1 (Hb9), a homeodomain transcription factor expressed by postmitotic MNs [28], and the cell cycle inhibitor p21 [29]. Olig2 also has the capacity to form heterodimers with Ngn2 that inhibit Ngn2's neurogenic activity [28], and oppose the functions of the tumor suppressor protein p53 [30]. Furthermore, addition of Olig2 to TF reprogramming cocktails inhibits conversion of fibroblasts to MNs [31], 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][33][34][35][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 [37]. 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 Oli-g2 LOW state, during which Hes1 expression decreases and Olig2 is coexpressed with Hes5, and a subsequent Olig2 HIGH 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 [37]. 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 [37]. 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][38][39][40].
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][42][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 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 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. 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 [37], 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][45][46][47].
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).  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. 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][42][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. 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.
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 [48]. 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  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 [51]. 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 [52]. 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 Oli-g2-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.
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 ( 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][58][59]. Single cell transcriptomics of motor neuron differentiation 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 Oli-g2 LOW expression to Olig2 HIGH . We therefore tested whether inhibiting Notch signalling upregulated 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 (

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) [60]. 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 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 mKate2 HIGH 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. https://doi.org/10.1371/journal.pbio.2003127.g004 Single cell transcriptomics of motor neuron differentiation  [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.
We next asked whether Olig2 might be responsible for the repression of Hes1 and Hes5 using Olig2 Cre knock-in mice [56,61]. In these mice, the Olig2 coding sequence has been replaced with bacteriophage P1 Cre recombinase (Cre) [56]. Thus, Cre protein expression demarcates the presumptive pMN domain in heterozygous control and in homozygous Olig2 Cre/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 0 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.
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 [13]. Whereas mammals have a single Hes5 gene, birds contain three Hes5 paralogs, termed HES5-1, HES5-2, and HES5-3 [62], 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; [13]). Based on these results, we conclude that Olig2 expression suffices to repress Hes gene expression in NPs.
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 [13]. 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   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][9][10][11][12][13][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. 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) ( [63], 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 [66]. 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 0 to the transcriptional start site for Hes5 in mouse, 10.5 kb 5 0 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 [28]. 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, 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 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).
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 Ebox 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).

Discussion
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).

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][34][35][36]. Here, we use pseudo-temporal ordering of 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.  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 [48]. 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 [28]. 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 [28]. 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][71][72][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 [74]. This phosphorylation appears to regulate the preference of Olig2 for open or closed chromatin and, thus, strongly influences its biological activity [71]. 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][76][77][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 [31]. 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 [82]. In cortical progenitors, fluctuations in Olig2 levels have also been documented, but these oscillations occur at a significantly slower frequency [65] 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 [84]. Similarly, sustained activation or inhibition of the Notch pathway causes, respectively, a ventral expansion or recession in p3 progenitors, located ventral to the pMN [61]. 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.

Animal welfare
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. Olig2 Cre and Ngn2 KIGFP 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 [85].

Differentiation of NPs from mouse ESCs
NPs were differentiated as described previously [37]. 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 10 6 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. qPCR analysis 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  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.

Flow cytometry
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 mKate2 HIGH cells in Fig 4L and S6G Fig, a threshold  Electrophoretic mobility shift assays pCS2+ plasmid expression vectors for Olig2, E12, Ngn2, and Id1 [13] 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. 32 P-dCTP-labeled probes were generated by Klenow end labeling of double-stranded oligonucleotides containing the native mouse Hes5(e1) sequence (forward 5 0 -ggccgCTCCCAAAAGACCATCTGGCTCCGTGTTATAA-3 0 ; reverse 5 0 -actag TTATAACACGGAGCCAGATGGTCTTTTGGGAG-3 0 ) or an E-box mutated version (forward 5 0 ggccgCTCCCAAAAGAggATCccGCTCCGTGTTATAA-3 0 ; reverse 5 0 -actagTTATAA CACGGAGCggGATccTCTTTTGGGAG-3 0 ). 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 0 -gaggcggccgcCGGTTCCCACACTTTGGT-3 0 ; reverse 5 0 -gagactagtCACAGTCCCAAGCTGCTTAAA-3 0 , and chick, forward 5 0 -gaggcggccgcTGCGTTT CCCATACTTTTCC-3 0 ; reverse 5 0 -gagactagtTCTGGCCTTGAAGCTAGGAG-3 0 . 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 [92]). 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.