A Shh/Gli-driven three-node timer motif controls temporal identity and fate of neural stem cells

A Shh/Gli-driven incoherent feed-forward loop circuit encodes time and counterbalances noise in temporal neural patterning.


INTRODUCTION
Time is a central axis of information during embryogenesis, but few mechanisms explaining the timing of developmental events have been resolved (1)(2)(3)(4). In the forming central nervous system (CNS), defined pools of multipotent neural stem cells (NSCs) produce distinct cell types in a specific sequential order and over defined time frames. In this process, aging NSCs become progressively restricted in their developmental potential by losing competence to generate earlyborn cell types (5), and genome-wide analyses have revealed that NSCs undergo dynamic transcriptional changes over time (6). However, the composition and functional properties of time-encoding circuitries determining time frames and point of transitions have not been resolved in any model system (7,8). Temporal neural patterning in vertebrates is a slow process progressing over days or even weeks depending on species (9,10). Yet, in several lineages, progenitors undergo coordinated and fast temporal transitions (9)(10)(11)(12)(13). Biological timers regulating temporal neurogenesis are therefore likely to exhibit properties that counterbalance noise in regulatory networks, but how this is achieved at the molecular level remains unknown.
Temporal patterning contributes to the generation of neural cell diversity at all axial levels of the CNS, but only few transcription factors (TFs) and/or signaling molecules regulating temporal fate and potency have been defined in various temporal lineages of the vertebrate CNS (11,13). To approach the question of time, we focused on a relatively well-defined lineage in the ventral brainstem that sequentially produces motor neurons (MNs), serotonergic neurons (5HTNs), and oligodendrocyte precursors (OPCs) (12,14). The lineage is induced by Sonic hedgehog (Shh) and defined by the expression of the TF Nkx2.2 (12), and the temporal progression of differentiation is easy to monitor as the NSC pool remains at a fixed position and does not comprise the specification of proliferative intermediate progenitors as in the developing neocortex (13). Nkx2.2 + NSCs show, otherwise, many common features to cortical progenitors, genetic lineage-tracing studies support that young NSCs are multipotent and competent to generate early-and late-born neurons, and aging NSCs become progressively restricted in their potential over time (12,13,15,16). In addition, the transition from early to late phases of neurogenesis is governed by late-acting extrinsic signals whose activation is intrinsically programmed within the lineage (15,17,18). Young Nkx2.2 + NSCs coexpress early-and late-acting fate determinants (12,15,(19)(20)(21), but the activity of the TF Phox2b predominates by specifying MN fate (15,21). Once Phox2b is downregulated or genetically ablated, MN production is terminated and 5HTNs are generated by default (12), suggesting that Phox2b functions as temporal effector output. Activators of Phox2b have not been defined, but a self-sustained and temporally delayed activation of Tgf operates as an extrinsic signal that triggers MN-to-5HTN fate switch by repressing Phox2b (15). Thus, Phox2b and Tgf are important regulatory components of a timer circuitry, but to understand how time is set by the network, it is necessary to define activators of Phox2b and resolve how the temporally gated activation of Tgf is mechanistically implemented. Here, we define a Shh/Gli-driven three-node circuitry, which explains how time is encoded in the Nkx2.2 + lineage and provide evidence that intrinsically programmed activation of extracellular switch signals fulfill an important role to counterbalance noise through spatial averaging, which is unattainable with temporal networks based exclusively on intrinsic transcriptional regulators.
To define genome-wide transcriptional changes over time, we determined the transcriptome of Nkx2.2 + NSCs isolated at different time points by RNA sequencing (RNA-seq). In this differentiation paradigm, >90% of Sox3 + NSCs express Nkx2.2 at 3.5 days in differentiation conditions (DDC) (15). Since the temporal patterning process is slow and progresses over days (15), cells isolated at a given time were examined at the population level to even out changes of gene expression associated with cell cycle progression and neuronal differentiation (22). We hypothesized that activators driving Phox2b are progressively lost over time since Phox2b becomes down-regulated and cells undergo a MN-to-5HTN fate switch in the absence of Tgf signaling, but on a delayed temporal schedule (15). Therefore, we determined genes down-regulated between 3.5DDC, when Phox2b expression and MN production has been initiated, and 5.5DDC when Phox2b becomes down-regulated in NSCs (Fig. 1A) (15). Approxi-mately 2200 genes showed a significant decline over this time frame [P ≤ 0.05; log 2 fold change (FC) ≥ 0.22]. The number of genes was markedly reduced with increased FC stringency, and only 27 genes remained at a FC of ~5.6 [log 2 (FC) ≥ 2.5] (fig. S1A). Among these were Gli1, Gli2, and Gli3 (fig. S1A), which encode zinc-finger TFs that transduce Shh signaling within the nuclei (23).
Gli2 and Gli3 are bifunctional TFs that are processed into repressors (GliR) in the absence of Shh but stabilized as full-length activators (GliA) within cilia by Smo (Smoothened) in response to binding of Shh to Ptch1 (Fig. 1B) (23). Gli1, in turn, is an obligate activator and a direct target of the Shh pathway (Fig. 1B) and is together with Ptch1 commonly used as an indicator of ongoing Shh signaling (24). Gli1 and Ptch1 exhibited a biphasic temporal expression: they were up-regulated during ventralization of NSCs and induction of the Nkx2.2 + lineage (~1DDC to 3DDC) followed by a progressive decline over the Scatter plot of nuclear protein levels of Gli1 with Nkx2.9 or Phox2b in Nkx2.2 + nuclei at 3.5DDC. Trend line in red; r, correlation coefficient. (F, I, and M) Error bars, means ± SEM; asterisks, Student's t test, ***P ≤ 0.001. MNw, motor-neuron window; 5HTNw, serotonergic-neuron window; A.U., arbitrary units.
Biochemical analyses showed that the down-regulation of Gli transcription was translated into a progressive decay of full-length activator forms of Gli2 (Gli2A, 185 kDa) and Gli3 (Gli3A, 190 kDa) between 2.5DDC and 5.5DDC (Fig. 1, F and G). Gli1 showed a similar profile but reached peak values at 3.5DDC, consistent with the notion that Gli1 is a target of Gli2A/Gli3A (Fig. 1, F and G). Gli1-3 proteins were detected in nuclei and Arl13b + cilia at 3DDC but not at 5.5DDC ( Fig. 1H and fig. S1, C and D). Smo translocates to cilia in response to binding of Shh to Ptch1 (23) and localized to Arl13b + cilia between 3DDC and 5.5DDC (Fig. 1H). Thus, down-regulation of Gli transcription mediates a progressive desensitization of the Shh pathway through loss of GliA, despite that Shh/Smo signaling is active over this period. The progressive loss of Gli1-3 expression correlated with the down-regulation of Phox2b and Nkx2.9 (fig. S1B). Inhibition of Shh/Smo signaling by CyC subsequent to the establishment of the Nkx2.2 + lineage resulted in premature down-regulation of Phox2b and Nkx2.9 in Nkx2.2 + NSCs and increased generation of Gata3 + 5HTNs (Fig. 1I). Expression of Nkx2.2 was not affected by CyC treatment initiated at 2.5DDC (Fig. 1I) consistent with data suggesting that induction but not maintenance of Nkx2.2 expression relies on Shh signaling (28). These data show that Shh/Gli signaling is required for sustained expression of Phox2b and Nkx2.9 and suggest that inhibition of Shh signaling is sufficient to trigger a MN-to-5HTN fate switch.
Despite the overall decline of Gli proteins over time, there was notable fluctuations of Gli1 and Gli2 expression levels between individual Nkx2.2 + nuclei at a given time examined (Fig. 1, J and K). Gli1 is a readout of Shh signaling activity (29), and there was a direct correlation between Nkx2.9 and Gli1 expression levels at 3.5DDC (Fig. 1N). Moreover, the expression of Nkx2.9 declined over time in a manner similar to Gli proteins (Fig. 1, J to L, and fig. S1B), suggesting a GliA dose-dependent regulation of Nkx2.9. In contrast, there was no correlation between Phox2b and Gli1 expression levels ( Fig. 1O), implying that regulators, in addition to GliA, are likely to influence Phox2b expression. Nevertheless, the fraction of Nkx2.2 + NSCs expressing Phox2b remained largely constant between 3DDC and 4.5DDC and only dropped markedly at 5DCC to 6DDC when Gli expression approached undetectable levels ( Fig. 1, J, K, and M, and fig. S1B) implying that the GliA level sufficient to sustain Phox2b transcription is very low.

Ptch1-independent establishment of parallel temporal GliA and GliR gradients
Immunoprecipitation assays revealed the presence of processed repressor forms of Gli2 (75 kDa; Gli2R) and Gli3 (83 kDa; Gli3R), which declined similarly to their respective full-length activator forms between 2.5DDC and 5.5DDC ( Fig. 2A). Gli3A/Gli3R ratios remained largely constant, while Gli2A/Gli2R ratios decreased somewhat over time (Fig. 2B). Ptch1 is up-regulated in response to Shh signaling, in which feedback can inhibit the Shh pathway, as free Ptch1 promotes GliR formation through Smo inhibition (23). Analysis of differentiating Ptch1 −/− ESCs showed that the temporal expression profiles of Gli genes, Gli1-3 activator proteins, Nkx2.2, Nkx2.9, and Phox2b, and induction of Tgf2 at 6.5DDC were similar in Ptch1 −/− and wild-type (WT) ESC cultures (Fig. 2, C and D, and figs. S1B and S2, A and D). These data establish that Ptch1 function is largely dispensable in the temporal differentiation process. Gli2R and Gli3R were detected in Ptch1 −/− cultures and at GliA/GliR ratios similar to controls (Fig. 2, E and F). Ptch2 (30) was expressed at extremely low levels in control and Ptch1 −/− cells ( fig. S2, B and C). Gli2R and Gli3R were present also in Ptch1 −/− ESC cultures treated with Shh Ag 1.3 (SAG) (Fig. 2E), which activates Smo downstream of Ptch1 and Ptch2 (30). These data reveal that a fraction of bifunctional Gli proteins are processed into GliR forms in conditions of fully activated Smo and independent of Ptch1-mediated feedback. The amount of processed Gli2R and Gli3R as well as full-length Gli2 and Gli3 produced at a given time is instead determined by the level of Gli2, Gli3 transcription, respectively ( fig. S3). Thus, downregulation of Gli genes produces parallel declining GliA and GliR gradients in the lineage (Fig. 2G).

The temporal network conforms a three-node incoherent feed-forward loop circuitry
To define the effect of constant GliA input on temporal output, we generated mice in which Gli1 was constitutively expressed in the Nkx2.2 + lineage by crossing a ROSA26-Gli1 FLAG transgene (31) with a Nkx6.2-Cre mouse line (15) (hereafter termed Gli1 ON mice) ( fig. S4A). Ptch1 is down-regulated in Nkx2.2 + NSCs by embryonic day 11.5 (E11.5) but was sustained at high levels in Gli1 ON mice over this period (Fig. 3B), consistent with continuous GliA expression at high levels. Forced Gli1 expression did not result in any overt spatial patterning phenotype ( Fig. 3A and fig. S4B), and the number of Isl1 + trigeminal MNs was similar to controls at E10.5 ( fig. S4A), a time when most MNs have been specified (12). In controls at E11.5, MN production is terminated, and the generation of Pet1 + 5HTNs is initiated ( Fig. 3A) (12). Most sl1 + MNs have migrated laterally to form the trigeminal nuclei but some late-born MNs are still present in the migratory stream (MS) (Fig. 3A), and Phox2b, Ptch1, and Gli1 have been down-regulated and Tgf2 up-regulated in Nkx2.2 + NSCs (Fig. 3B). There was a surplus of Isl1 + MNs in the MS and reduction of Pet1 + 5HTNs in Gli1 ON mice at E11.5 (Fig. 3A), but the expression of Phox2b in Nkx2.2 + NSCs was repressed (Fig. 3B). This suggests a mild temporal extension of MN production in Gli1 ON mice but underpins that Phox2b is suppressed and cells undergo a MN-to-5HTN fate switch on an almost normal temporal schedule in conditions of constant GliA input (Fig. 3C).
The mild temporal phenotype in Gli1 ON mice was associated with a marked up-regulation of Tgf2 at E11.5 (Fig. 3B), revealing a feedforward propagation of Tgf2 transcription by GliA. This implies a three-node circuitry forming an incoherent feed-forward loop (IFFL), whereby GliA activates not only Phox2b but also the suppressive Tgf node negatively regulating Phox2b (Fig. 3H). This predicts that the MN window should be extended if GliA is maintained and Tgf concurrently inactivated, and to test this, we crossed Gli1 ON mice onto a Tgfbr1 mutant background. Removal of one or both copies of Tgfbr1 resulted in progressively more pronounced temporal phenotypes, as determined by quantification of Phox2b + NSCs and Isl1 + MNs and by 5-Ethynyl-2'-deoxyuridine (EdU)-birth dating experiments (Fig. 3, D to F, and fig. S5A). In Gli1 ON :Tgfbr1 −/− at E11.5, there was a massive accumulation of premigratory and migrating MNs, a complete lack of Pet1 + 5HTNs and maintained expression Phox2b in NSCs (Fig. 3, A and B). The extension of MN production was more pronounced in Gli1 ON :Tgfbr1 −/− as compared to Tgfbr1 −/− mutants (Fig. 3, E and G, and fig. S5B), supporting that the delayed fate switch in Tgfbr1 mutants (15) occurs due to depletion of GliA. In addition, Tgf2 expression was reduced in Gli1 ON :Tgfbr1 −/− relative to Gli1 ON mice ( Fig. 3B), suggesting that positive feedback signaling by Tgf is necessary for robust induction of Tgf2 downstream of GliA (Fig. 3H). Collectively, these data support an IFFL regulatory topology for the GliA-Phox2b-Tgf motif (Fig. 3H) and establish that Tgf predominates over the Shh pathway by suppressing Phox2b even if cells express GliA at levels sufficient to sustain Phox2b transcription.

Different GliR sensitivities account for sequential activation of Phox2b and Tgf2
GliA promotes both Phox2b and Tgf2, raising the key question how Tgf2 induction can be circumvented at early stages when GliA activity peaks? As the down-regulation of Gli genes produces parallel declining GliA and GliR temporal gradients (Fig. 2G), we considered that the delayed activation of Tgf2 could be explained by an inhibitortitration regulation, a regulatory motif known to convey nonlinear switch-like responses (32,33). To explore this possibility, we characterized the spatiotemporal dynamics of gene expression in Nkx2.2 + NSCs in vivo by semiquantitative in situ hybridization and quantitative polymerase chain reaction (qPCR). Induction of the Nkx2.2 + domain and subsequent MN-to-5HTN fate switch progress in a ventral-to-dorsal (V → D) manner (12,34,35), and ventrally located progenitors are thereby older than more dorsally located siblings at a given time. Consistent with this, down-regulation of Gli1-3, Ptch1, Phox2b, and Nkx2.9 progressed in an overall V → D fashion but at different kinetics (Fig. 4A). Gli2 and Gli3 were most rapidly down-regulated and had, by E10.5, become constrained to the dorsal third of the Nkx2.2 + domain, both in WT and Gli1 ON mice (Fig. 4, A and B, and fig. S6A). Consistent with the ESC differentiation system, immunoblot analysis revealed that Gli2 and Gli3 proteins were also down-regulated over time in ventral hindbrain tissue ( fig. S6C). There was a mutually exclusive relationship between Tgf2 and Gli2,Gli3 both in control and Gli1 ON mice (Fig. 4B). Up-regulation of Tgf2 in Gli1 ON mice was first detected at E10.5 in ventral cells that ceased to express Gli2 and Gli3 and which had initiated Tgf2 expression also in controls but at lower levels ( Fig. 4, A and B, and fig. S6B). This shows that a young Gli2 + /Gli3 + context is nonpermissive for GliA-mediated Tgf2 induction, while Tgf2 responds to GliA in a dose-dependent man-ner in older Gli2 − /Gli3 − cells. Accordingly, Gli2 + /Gli3 + cells must express an inhibitor that acts dominant negative over GliA and that is removed over time, as is the case for GliR (Fig. 2, A and G). In functional support for such a repressor function of GliR, we found that forced expression of GliR was sufficient to suppress GliA-mediated activation of Tgf2 in epistasis experiments (Fig. 4G). Removal of GliR is thus necessary for activation of Tgf2, and the kinetics of its elimination therefore determine the kinetics of Tgf2 induction. Because of the bifunctional nature of Gli2 and Gli3 proteins, we have not been successful in selectively eliminating GliR in the Nkx2.2 + lineage and cannot therefore rule out that early suppression of Tgf2 may involve repressors in addition to GliR. However, this would not change the overall topology of the timer circuitry, since these hypothetical repressors would then need to be removed either with the same kinetics or faster than GliR to permit activation of Tgf2 when the GliR concentration declines.
Unexpectedly, in Gli1 ON mice, Phox2b and Nkx2.9 expression were augmented in ventral Gli2 − /Gli3 − cells but reduced in dorsal Gli2 + / Gli3 + cells at E10.5 (Fig. 4B). These transcriptional outputs can only be explained if Phox2b and Nkx2.9 display certain GliR sensitivity, and that the failure to down-regulate Ptch1 in Gli1 ON mice results in reinforced Ptch1-mediated feedback inhibition of the Shh pathway. We observed a notable displacement of Smo out of cilia in Gli1 ON conditions over time (Fig. 4, C to E) and increased formation of Gli2R and Gli3R in Gli1 ON ESC cultures (Fig. 4F). Thus, increased GliR formation due to Ptch1-mediated feedback begins to suppress Phox2b and Nkx2.9 expression in young Gli2 + /Gli3 + cells in Gli1 ON conditions. This cannot occur once cells ceased to express bifunctional Gli proteins, thereby resulting in the anticipated Gli1-mediated up-regulation of Phox2b and Nkx2.9 in older Gli2 − /Gli3 − cells. Collectively, these data support a model whereby a high GliR sensitivity for Tgf2 prohibits GliA-mediated activation until GliR has been titrated out, thereby establishing a delayed activation of the Tgf node (Fig. 4H). Phox2b and other Shh-regulated genes display lower GliR sensitivities (Tgf2 > Phox2b, Nkx2.9 > Ptch1, Gli1), allowing early activation in the lineage (fig. S7). These data further establish that the decay of GliA followed by downregulation of Gli genes is functionally important, as it is necessary for evading Ptch1-mediated feedback inhibition, which otherwise begins to interfere with Phox2b expression during the MN window.

Gli1 is required for Tgf2 induction and prompt termination of MN production
A Gli inhibitor-titration regulation of Tgf2 requires that GliA levels remain high enough to activate Tgf2 once GliR has been titrated away, raising the possibility that the feed-forward activation of Gli1 by bifunctional Gli proteins functions to boost GliA activity late in the differentiation process. In direct support for this, we found that Tgf2 failed to be induced on time in Gli1 −/− mice (Fig. 5A). This was accompanied by prolonged expression of Phox2b in progenitors (Fig. 5, B and C, and fig. S8A), a moderate overproduction of Isl1 + MNs and concurrent reduction of Lmx1b + 5HTNs at E11.5 (Fig. 5, D and E) and Pet1 + 5HTNs at E12.5 ( fig. S8B). Tgf2 expression had been initiated at low levels by E12.5 in Gli1 mutants (Fig. 5A), presumably reflecting input by remaining Gli2/3A activity followed by propagation of Tgf2 expression by positive feedback signaling (Fig. 3B) (15). This shows that Gli1 is required for timed activation of Tgf2 and prompt termination of Phox2b and MN production (Fig. 5F), suggesting that Gli2A and Gli3A are the primary activators of Phox2b, while timed induction of Tgf2 also depends on Gli1. In addition, considering that Gli1 binds to DNA with lower affinity than Gli2 (36) and Tgf2 is a direct target of Gli proteins (37), this provides a mechanistic rationale for why Tgf2 is not accessible for Gli1-mediated induction until Gli2R/Gli3R are depleted.

Robust temporal switch by hysteresis and spatial averaging
The timer motif suggested by experimental data relies on the parallel temporal decay of a repressor and an activator, and the delayed activation of a diffusible repressor (Fig. 4H). For all progenitors to undergo a coordinated MN-to-5HTN switch in a short time window, cells have to down-regulate Phox2b synchronously, despite that Gli protein levels vary between cells (Fig. 1, J and K). If we model declining GliA and GliR levels with similar experimental noise levels and a noisy readout threshold (Fig. 6, A and B), then we predict a noisy decline of Phox2b (Fig. 6B) and a transition period from MNto-5HTN production that lasts for almost 2 days (Fig. 6C), both if we consider a direct regulation of Phox2b by GliA (Fig. 6, B and C, red line) or a relay via GliR and Tgf (Fig. 6, B and C, cyan line); we ignore the direct negative impact of GliR on Phox2b as it must be weak, given the strong expression of Phox2b while GliR levels are high. Tgf is selfactivating, and the introduction of this positive feedback steepens the response (Fig. 6B, blue line), especially when the positive feedback results in a hysteretic switch (Fig. 6B, black line), but does not shorten the transition period t, which we define as the time that it takes for the likelihood of MN differentiation to decrease from 0.95 to 0.05 (Fig. 6C). In case of a hysteretic switch, the system is bistable, and as the GliR levels decline, the system reaches a critical point (Fig. 6D, red dot) where it jumps from a low to a high Tgf steady state (Fig. 6, D and E, black line). The critical GliR level at which the system returns to the low Tgf steady state is at a much higher GliR level, making the switch effectively one way (38). Such a hysteretic switch is much more robust to noise in the readout process (Fig. 6F), but this advantage disappears when also considering the impact of kinetic noise (Fig. 6G and fig. S9A), as hysteretic switches are particularly sensitive to kinetic noise (Fig. 6, B and E). Therefore, how can the system achieve a fast, reliable transition? Tgf is a short range diffusible protein (39) that permits spatial averaging, and this leads to simultaneous fate switching on the population level despite noise ( Fig. 6H and fig. S9B). The diffusion range has little impact, and even next-neighbor interactions result in robust synchronized switching at the population level. Note that spatial averaging improves the robustness to noise for all network architectures, but the hysteretic switch achieves a substantially shorter transition period (Fig. 6H and fig. S9B).
Phox2b expression is controlled by both GliA and Tgf. Addition of the GliA link makes the switch less robust to noise compared to the GliR link alone (via Tgf relay), given the lack of spatial averaging. Accordingly, the highest robustness is achieved for a low GliA threshold, as supported by experimental data. Given that generation of Gli2A/Gli3A and Gli2R/Gli3R are coupled (Fig. 2), their levels may be somewhat correlated. However, even perfect correlation (Fig. 6I and  fig. S9C, green line) does not improve robustness compared to the uncorrelated case (orange) but does lengthen the time to the switch (Fig. 6, J and K). Thus, a low GliA threshold for activation of Phox2b is important to ensure that the time of Phox2b down-regulation is set by Tgf signaling, whose noisiness is reduced by spatial averaging, and reduces the impact of variations of GliA and GliR levels on the duration of the switching interval.
While the timer network is robust to noise, it remains sensitive to structural changes, which allows us to computationally recapitulate the various mutants discussed above (Fig. 6K), thereby confirming the above verbal reasoning that led to the modeled network architecture (Fig. 4H). The combination of a decaying activator and repressor introduces a further save guard as becomes visible when either component is removed (Fig. 6K). Changes in the Gli decay rate (Fig. 6A) alter the time point at which the fate switch occurs (Fig. 6L). In case of an exponential loss of the signal (Fig. 6A), the time point changes proportionally to the inverse of the decay rate (Fig. 6L).

DISCUSSION
How time is measured by NSCs during temporal patterning of neurons has remained unresolved. In this study, we define a Shh/Gli-driven three-node timer circuitry underlying the sequential specification of MNs and 5HTNs by NSCs. The network is founded on a parallel temporal decay of GliA and GliR established through a progressive down-regulation of Gli1-3 transcription. Regulatory interactions conform an IFFL circuitry in which GliA not only promotes Phox2b expression and MN fate but also accounts for a delayed activation of a suppressive Tgf node that triggers a MN-to-5HTN fate switch by repressing Phox2b. Data suggest a model in which activation of the Tgf node is temporally gated by a Gli inhibitor-titration mechanism, whereby induction of Tgf2 by Gli1 is prohibited until GliR has been titrated out. We show that GliR is generated also in the absence of Ptch1-mediated feedback inhibition and when the Shh pathway is fully activated (Fig. 2, E and F). Down-regulation of Gli genes is consequently necessary for elimination of GliR and initiation of a fate switch by Tgf (Figs. 2 and 4H). Changes in the rate of Gli decay alter temporal output by the circuitry and conceptually explain how time is encoded in the lineage. In relation to this, Nkx2.2 has been shown to feedback inhibit Gli genes (35,40), and the MN window is extended in Nkx2.2 −/− mice (12), suggesting that negative feedback by Nkx2.2 itself impacts on the pace of Gli decay and time output. Such feedback modulation of Gli genes illustrate how time can be tuned in a lineage-specific manner, which is needed for the generic deployment of the core-timer circuitry in multiple Shh-induced lineages, producing distinct neural progenies over different time frames (15,41,42).
In the adult brain, the transition of dormant NSCs into transient amplifying cells is accompanied by down-regulation of Gli genes and is inhibited by GliR (43), implying that removal of GliR through adaptation of Gli transcription accounts for the transition between these cellular states. Together with our study, this suggests that Gli inhibitor titration and the idea that down-regulation of Gli genes facilitates induction through removal of GliR may provide a common mechanism to attain switch-like responses in Shh-regulated differentiation processes.
Biological timers based on accumulation or titration of activators or repressors have previously been reported (1,2,4), and Gli decay is formally sufficient to mediate timer function without Tgf, raising the question why timing of the MN-to-5HTN fate switch involves a more complex network architecture. Computational modeling suggest that spatial averaging enabled by the diffusible and selfactivating properties of Tgf, in combination with hysteresis, produces prompt suppression of Phox2b and a coordinated switch at the population level (Fig. 6, H and K). Integration of the Tgf node thereby counterbalances noise and generates a more precise timer mechanism as compared to a timer based only on Gli decay. Spatial averaging further provides potential to reset temporal synchrony at the population level, which should benefit coordination of subsequent fate transitions, as sequentially occurring switches has been shown to be temporally coupled (15). These community features are not attainable with temporal networks based exclusively on intrinsic transcriptional regulators. Our data consequently provide a functional basis for the intrinsically programmed activation of extrinsic signals in temporal neural patterning processes, and it seems likely that late-acting extrinsic cues implicated in temporal neurogenesis in the cortex (18) function in a manner analogous to Tgf.

Mouse models
The following published mouse strains were used in this study: Nkx6.2::Cre and Tgfbr1 fl/fl (15), ROSA26-Gli1 FLAG (31), Gli1 −/− (44), and Nkx2.2 −/− (12). WT and experimental mice were maintained in a C57BL/6 background. Male and female mice were both used depending on availability and were housed in breeding pairs or grouphoused with littermates of the same sex after weaning (two to five mice per cage) on a 12-hour light/dark cycle, with food and water provided ad libitum. All experimental procedures were conducted in accordance to the Swedish Animal Agency guidelines for animal experimentation and approved by the regional animal ethics committee of Northern Stockholm.

p-MACS of NPCs
For RNA-seq, qPCR, immunoblotting, and immunoprecipitation experiments neural progenitor cells (NPCs) from 3.5DDC and onward were isolated from differentiating ESCs cultures by Prominin1-based magnetic sorting (p-MACS) according to the manufacturer's instructions (Miltenyi Biotec). Briefly, cells were dissociated using TrypLE Express, washed twice with ice-cold PBS supplemented with 0.5% bovine serum albumin and 2 mM EDTA (Sigma-Aldrich), and passed through a 40-m cell strainer. Dissociated cells were incubated with magnetic beads conjugated with anti-prominin-1 antibody. Prominin1 + cells were separated on magnetic columns and processed for biochemical analysis. To assess purity of sorted fraction, cells were plated, incubated for 2 hours, and fixed for immunocytochemistry.

Chick electroporation
Fertilized chick eggs were stored up to 1 week at 16°C, and at the start of experiment, eggs were placed horizontally in a 38°C humidified chamber for ~40 hours until the developing embryo reached a Hamilton-Hamburger stages 11 and 12. Approximately, 5 ml of white was removed to lower the embryo, and an opening in the shell on the top of the egg was made to expose the embryo, which was visualized under a magnifying microscope. A mixture of DNA plasmid(s) and Fasta green (used for coloring) was injected into the central canal of the chick embryo using a pulled capillary needle attached to a mouth aspirator tube. Several drops of PBS were placed on top of the embryos, the electrodes were positioned manually one at each side of the embryos at the hindbrain level, and two 25 ms 6-V electric pulses 1 s apart were applied to the embryo using an Electro square Porator ECM 830. The eggs were sealed and incubated at 38°C for 40 hours before embryo collection and processing for immunofluorescence or in situ hybridization.

Immunohistochemistry
Embryonic stage of collected mouse embryos was determined by somite number (E9.5 to E11.5 embryos) and by limb anatomy (E10 to E12.5 embryos). Mouse embryos from E9.5 to E11.5 5 and chick embryos collected 40 hours after electroporation were fixed in 4% paraformaldehyde in 0.1 M PBS for 2 hour on ice with shaking, and E12.5 mouse embryos were fixed overnight at 4°C. Fixed tissue was washed with PBS, cryoprotected by equilibration in 30% sucrose in PBS, embedded in O.C.T. TM , frozen on dry ice, and cryostat-sectioned in the transverse plane at 12 m. Tissue sections were collected at hindbrain level, and rhombomeric level was determined by anatomical landmarks (i.e., otic vesicle) and Phox2b/Isl1 immunohistochemistry to visualize MNs. Immuhistochemistry was carried out in water-humidified chamber. The tissue on slides was incubated with blocking solution (3% FCS/0.1% Triton X-100 in PBS) at room temperature (RT) for 1 hour followed by incubation with primary antibodies overnight at 4°C and fluorophore-conjugated secondary antibodies for 1 hour at RT. Both primary and fluorophore-conjugated secondary antibodies were diluted in blocking solution. Sections were mounted using Vectashield and coverslipped for imaging.

In situ hybridization
Mounted tissue sections were postfixed in 4% paraformaldehyde (PFA) in phosphate buffer for 10 min, washed three times for 5 min with PBS, incubated with proteinase K solution [1 g/ml proteinase K in 50 mM tris-Cl (pH7.5) and 5 mM EDTA] for 5 min and PFAfixed and washed with PBS as before. Slides were subsequently treated with an acetylating solution [2 mM HCl water solution containing triethanolamide (14 l/ml) and acetic anhydride (2.5 l/ml)] for 10 min, washed three times with PBS, and incubated with hybridization solution [50% formamide, 5× saline-sodium citrate (SSC), 5× denharts, yeast RNA (250 g/ml), herring sperm DNA (500 g/ml), and blocking reagent from Roche (20 mg/ml)] for 1 hour in a chamber humidified with a 50% formamide/5× SSC solution. All steps were performed at RT. Slides were subsequently incubated overnight at 70°C with labeled probe diluted in hybridization solution, followed by washes with a 0.2× SSC solution, first at 70°C for 1 hour and then by another for 10 min at RT. Slides were incubated in B1 solution [0.1 M tris-Cl (pH7.5), 0.15 M NaCl, and 10% heat inactivated fetal bovine serum (FBS)] for 1 hour at RT, followed by an overnight incubation with anti-digoxigenin (DIG)-alkaline phosphatase Fab fragments antibody diluted in B1 solution at 4°C in a water humidified chamber. Excess antibody was washed away with three washes at RT with B1 solution and equilibrated in B3 solution [0.1 M tris-Cl (pH 9.5), 0.1 M NaCl, and 50 mM MgCl 2 ] before incubation of tissue with developing solution [10% polyvinyl alcohol, 100 mM tris-Cl (pH 9.5), 100 mM NaCl, 5 mM MgCl 2 , and levamisol (0.24 mg/ml)] containing nitro blue tetrazolium/bromochloroindolyl phosphate (Roche) at RT. The developing reaction was stopped by placing the slides in water and after several washes, and slides were mounted using Aquatex mounting medium (Merck). Labeled probes were produced from in vitro transcription of linearized plasmid containing a sequence specific for the desired gene using a DIG RNA labeling kit (Roche) according to the manufacturer's protocol. Probes were purified using Microspin G-50 columns (GE Healthcare) and stored at −20°C.

EdU Labeling
EdU (Life Technologies) was injected intraperitoneally in pregnant female mice at 0.04 mg/g of body weight and detected using the Click-iT EdU Alexa Fluor 555 Imaging Kit (Life Technologies) according to the manufacturer's protocol.

Isolation of mouse neural tissue for qPCR and Western blot
Neural tube tissue was isolated as described before (45). In brief, mouse embryonic tissue was dissected in L-15 medium and staged based on somite number and limb morphology, and hindbrain region was roughly dissected out with the help of tungsten needles. To completely remove the mesenchymal tissue surrounding the neural tissue, the dissected hindbrain piece was treated sequentially with dispase solution [dispase (1 mg/ml) in L-15 medium] for 5 min, 10% FBS solution (in L-15 medium) for 5 min and washed twice in L-15 medium. The remaining mesenchymal tissue was then removed with the help of tungsten needles. For qPCR analysis, the ventral region of rhombomeres 2 and 3 of the hindbrain was isolated, while for Western blot analysis, the ventral and intermediate domains of the hindbrain were collected. At the stages analyzed, the rhombomeric units of the hindbrain can be easily defined by grooves in the neural tissue. RNA for each individual piece was isolated using the RNeasy Micro Kit with a deoxyribonuclease digestion step (Qiagen) according to the manufacturer's protocol. Whole-cell protein extracts and Western blot were performed as described in the "Immunoblotting and immunoprecipitation" section.

Immunoblotting and immunoprecipitation
For whole-cell extracts, p-MACS isolated NPCs were lysed in radioimmunoprecipitation assay buffer (Sigma-Aldrich) complemented with protease and phosphatase inhibitor cocktail (Thermo Fisher Scientific) and incubated on ice with shaking for 30 min. Lysate was cleared by centrifugation (20,000g for 20 min at 4°C) and protein concentration determined by bicinchoninic acid assay. Protein lysate was resuspended in NuPAGE TM LSD sample buffer (Thermo Fisher Scientific) with 2.5% 2-mercaptoethanol and denaturated at 95°C for 5 min. Fifteen to 30 g of protein were loaded per lane of a 10% SDS polyacrylamide gel (Bio-Rad) and transferred onto nitrocellulose membranes (Bio-Rad) using a Trans-Blot Turbo System (Bio-Rad). Membranes were incubated 1 hour in blocking solution [tris-buffered saline with 0.1% Tween 20 (TBST) and 5% nonfat dry milk], followed by overnight incubation at 4°C with primary antibodies. After three washes with TBST, membranes were incubated with horseradish peroxidase (HRP)-conjugated secondary antibodies for 1 hour at RT. Detection of HRP was performed by chemiluminescent substrate SuperSignal West Dura substrate, and signal was detected on a ChemiDoc Imaging System (Bio-Rad).
For immunoprecipitation assays, 300 g of protein lysate was diluted with Pierce IP lysis buffer (Thermo Fisher Scientific) up to 1 ml and incubated overnight at 4°C with 2 g of primary antibody. Proteinantibody complex was isolated with protein A/G magnetic beads according to the manufacturer's instructions (Thermo Fisher Scientific). Precipitated material was used for immunoblotting as described above. Primary antibodies used are as follows: mouse anti-Gli1 (Cell Signaling Technology), goat anti-Gli2, Gli3 (R&D Systems), mouse antiglyceraldehyde-3-phosphate dehydrogenase (GAPDH) (Thermo Fisher Scientific), and mouse anti-actin (Seven Hills Bioreagents).

Image analysis
For quantification of protein level from immunofluorescence, confocal images were analyzed in ImageJ 1.48v. Nuclear area was defined by 4′,6-diamidino-2-phenylindole, Nkx2.2, or Sox1 immunostaining, and integrated density of fluorescent signal in individual nuclei was determined.
Images from in situ hybridization were acquired in a Zeiss Axio Imager.M2 microscope coupled to an Axio camera 503 mono and processed in Photoshop CS6. Images in Fig. 5 (A and B) were processed in ImageJ 1.48v using LUT-Fire function.
To determine relative gene expression levels along the Nkx2.2 domain (Fig. 4), images were analyzed in ImageJ 1.48v using "Plot Profile" function. A region of tissue that did not express the gene was used to determine background levels.
Western blots were quantified in ImageJ 1.48v using Gel Analysis Tool.
RNA-seq, gene expression quantification, and differential expression analysis RNA concentration and integrity were determined on an Agilent RNA 6000 Pico chip using Agilent 2100 BioAnalyzer (Agilent Technologies). RNA-seq of condition-specific samples was performed by the National Genomics Infrastructure at the Science for Life Laboratory in Stockholm on Illumina HiSeq2500 in RapidHighOutput mode with single-end setup (1× 50-base pair read length). Reads were mapped to the mouse genome assembly build GRCm38 using TopHat (v2.0.4). Next, the number of reads mapped to each gene was calculated using htseq-count (v0.6.1). The gene level abundances were estimated as fragments per kilobase of transcript per million mapped reads (FPKMs) using Cufflinks (v2.1.1). Plots and histograms were made using FPKM values. Furthermore, we processed the read count data with RNA-seq-specific functions of R package limma. For adjusting the mean-variance dependence, we used the function voom (46), which rendered the expression data into normally distributed log 2 counts per million values, accompanied with observationlevel gene-specific adjusting weights. The differential expression was estimated with functions lmFit, eBayes, and topTable. The variance estimates were obtained by treating all samples as replicates (design = NULL) and obtaining library sizes from counts (lib.size = NULL) without further normalization (normalize.method = "none"). The output FC values of differential expression were accompanied with P values. The latter were adjusted for multiple testing by calculating Benjamini-Hochberg's false discovery rate (47).

Identification of genes changing over time
To determine genes down-regulated in NSCs between 3.5DDC and 5.5DDC, FPKM values at each time point, P values, and FC values calculated as described above were used, and the following criteria were applied: FPKM (3.5DDC) ≥ 2, P ≤ 0.05, and downregulated (log 2 FC ≥ 0.22). In addition, only protein-coding genes were considered.
In addition, only protein-coding genes were considered.

Theoretical framework
Our theoretical framework considers the dynamics of GliA (A), GliR (R), Tgf (T), and Phox2b (P). We approximate the measured dynamics of GliA and GliR (Fig. 1N) by an exponential decay from initial values The equations for Tgf and Phox2b differ for the different regulatory networks. A generalized form is given by  (Table 1). Thus, in the case of  = 0, H w+ = 1, and H w− = 1, and the link represented by this function is not being considered. In the case of  = 1, H w+ = H + , and H w− = H − represent a positive and negative Hill function, respectively. We use v T as production and decay rates, so that the maximal concentrations Tgf and Phox2b is one and their values are restricted to the interval [0, 1].
To represent the effects of molecular noise, we added scalar white noise to the equations, i.e., a noise term that is proportional to the values of the variables. Therefore, the ordinary differential equations of the form dx = f(x, y, z ) dt become stochastic differential equations of the form dx = f(x, y, z ) dt +  x dW where W denotes a Wiener process and  denotes the kinetic noise level.
In Fig. 6 (I to K), we have considered both correlated and uncorrelated kinetic noise in GliA and GliR. In the correlated case (green lines), we used the same values for GliA and GliR, while in the uncorrelated case (yellow lines), the values are obtained by independently solving Eqs. 1 and 2, with a noise term.
To evaluate the impact of noise in the Phox2b concentration on the likelihood of a cell to differentiate into either a MN or 5HTN, we defined the likelihood of MN differentiation as the cumulative distribution function of a gamma distribution with different levels of variance. More specifically, we defined the shape (k) and the scale () parameters of the gamma distribution as k = 10/ and  = 0.04 , where  represents the threshold noise level. The average of this gamma distribution is given by k = 0.4 for all , while the variance is proportional to  and given by 1.6 10 −2 .
The effect of Tgf diffusion (Fig. 6, H to K) was introduced by replacing the variable T in the term H w+ (T, h TT ,  TT ) in Eq. 3 and the term H w− (T, h TP ,  TP ) in Eq. 4 with 0.5T + 0.5〈T〉 nn , where 〈T〉 nn represents the average values of Tgf in the nearest neighboring cells. Averaging over a wider range of neighbors yields similar results. In our simulations, we have considered an array of 100 × 100 cells with periodic boundary conditions.

Mutants
In the table, the full network represents the WT case. For the Tgfbr1 −/− mutant, we removed the Tgf modulation on Phox2b by setting  TP = 0.0. On the basis of data in Fig. 2, we used a 30% lower starting value for GliA in the Gli1 −/− mutant compared to WT. Last, in the Gli ON Tgfbr1 −/− mutant, the GliA concentration remains always well above the threshold for Phox2b production. Given that Tgf is absent, the Phox2b concentration always remains above the threshold for MN formation.

Parameter values
In the simulations, we have used  A =  R = 0.86, v T = 50, and v P = 90. The other parameter values vary according to the model and are presented in table below.

Statistical analysis
Quantitative data are reported as means ± SD or as means ± SEM or as box plots (with whiskers entending to 5th and 95th percentile) as specified in figure legends. The Student's t test or pairwise Wilcoxon test was used for evaluation of significance as specified in figure legends. The sample number (n) indicated in figure legends refers to the number of biological replicates. n.s. = not significant (P > 0.05), *P ≤ 0.05, **P ≤ 0.01, and ***P ≤ 0.001.

SUPPLEMENTARY MATERIALS
Supplementary material for this article is available at http://advances.sciencemag.org/cgi/ content/full/6/38/eaba8196/DC1 View/request a protocol for this paper from Bio-protocol.