Research Article

Spatiotemporal structure of cell fate decisions in murine neural crest

See allHide authors and affiliations

Science  07 Jun 2019:
Vol. 364, Issue 6444, eaas9536
DOI: 10.1126/science.aas9536

Binary decisions refine fate decisions

Neural crest cells develop into tissues ranging from craniofacial bones to peripheral neurons. Combining single-cell RNA sequencing with spatial transcriptomics, Soldatov et al. analyzed how neural crest cells in mouse embryos decide among the various fates available to them (see the Perspective by Mayor). These multipotent cells become biased toward a given fate early on and step through a progression of binary decisions as their fate is refined. Competing fate programs coexist until increased synchronization favors one and repression disfavors the other.

Science, this issue p. eaas9536; see also p. 937

Structured Abstract

INTRODUCTION

Multipotent progenitors must choose among multiple downstream fates. In developing embryos, progenitor cells exhibit transcriptional or epigenetic heterogeneity that is related to early biases in cell fate choices, and can be externally induced or stochastic in nature. Molecular assessment of the transient states assumed by cells during these developmental progressions has the potential to illuminate how such fate-specific biases emerge and unfold to ensure fate commitment. With this aim, we examine multipotent neural crest cells—transient embryonic progenitors unique to vertebrates that build the head, teeth, neuroendocrine tissue, and autonomic and sensory nervous systems. Cranial neural crest preferentially gives rise to a multitude of mesenchymal types of facial cartilage and bones, in addition to neuronal, glial, and pigment cell–type progeny. By contrast, trunk neural crest does not form bone or cartilage derivatives in vivo. The logic and molecular mechanisms that allow neural crest to resolve multiple potential cell fates at each axial level remain poorly understood.

RATIONALE

Here we used single-cell and spatial transcriptomics with statistical analysis of branching trajectories to investigate lineage relationships in mouse neural crest. Combined with lineage tracing and functional perturbations, we addressed spatiotemporal dynamics associated with early cell fate decisions in mouse trunk and cranial neural crest cells with different fate potential.

RESULTS

We find that up to early migration, neural crest cells progress through a sequence of common transcriptional states, followed by fate bifurcations during migration that can be formalized as a series of sequential binary decisions. The first decision separates sensory neuro-glial fate from all other fates, whereas the second decision occurs between autonomic and mesenchymal lineages and reveals a bipotent Phox2b+/Prrx1+ subpopulation. Decision points uncover distinct roles of neural crest regulators: Neurog2 is involved in early repression of melanocytes and activation of sensory fate at later steps. Each decision consists of initial coactivation, gradual biasing, and commitment phases. Early genes of competing cell fate programs coactivate in the same cells, starting from premigratory stage. As cells approach cell fate bifurcation points, increased synchronization of fate-specific programs and repulsion of competing fate programs lead to gradual appearance of cell fate bias, which becomes pronounced upon neural crest migration. Cell fate commitment culminates with activation of mutually exclusive, fate-specific gene expression programs. Early transcriptional patterns reveal that fate biasing of neural crest is already detectable when neural crest cells delaminate from the neural tube. In particular, the neuronal bias of trunk and mesenchymal bias of cranial neural crest emerge during delamination, indicating that this might be the time when the mesenchymal potential, distinct between cranial and trunk neural crest, is installed. In support to this hypothesis, we find that sustained overexpression of a single gene, Twist1, normally activated upon delamination only in the cranial compartment, is sufficient to reverse the trunk crest developmental program to a mesenchymal route.

CONCLUSION

Our analysis resolved a branching transcriptional trajectory of the differentiating neural crest, illustrating transcriptional implementation of major cell fate decisions and pinpointing the key differences defining cranial versus trunk neural crest potential. Our results show that neural crest cells differentiate through a series of stereotypical lineage-restriction events that involve coexpression and competition of genes driving alternative fate programs.

Cell fate acquisition in neural crest.

(A) Stepwise model of neural crest developmental dynamics. EMT, epithelial-to-mesenchymal transition. (B) Key phases of fate decision in neural crest cells. (C) Different fate potential of cranial versus trunk neural crest is encoded by the mesenchymogenic genes Twist1 and Prrx2. A single factor, Twist1, confers mesenchymal potential to trunk neural crest.

Abstract

Neural crest cells are embryonic progenitors that generate numerous cell types in vertebrates. With single-cell analysis, we show that mouse trunk neural crest cells become biased toward neuronal lineages when they delaminate from the neural tube, whereas cranial neural crest cells acquire ectomesenchyme potential dependent on activation of the transcription factor Twist1. The choices that neural crest cells make to become sensory, glial, autonomic, or mesenchymal cells can be formalized as a series of sequential binary decisions. Each branch of the decision tree involves initial coactivation of bipotential properties followed by gradual shifts toward commitment. Competing fate programs are coactivated before cells acquire fate-specific phenotypic traits. Determination of a specific fate is achieved by increased synchronization of relevant programs and concurrent repression of competing fate programs.

Multipotent progenitors acquire one of the multiple downstream fates. In developmental systems, progenitor cells exhibit transcriptional or epigenetic heterogeneity related to early biases in cell fate choices and can be oscillatory or stochastic (13). Understanding how fate-specific biases emerge and unfold to ensure fate commitment requires deeper understanding of cell profiles during progenitor, transient, and derivative states. Toward that goal, we examined a continuum of states in the neural crest (NC) lineage at different points along the murine rostrocaudal body axis during embryogenesis. As neural crest cell fate decisions are traceable, irreversible, and produce well-known differentiated cell types, we were able to investigate the interplay of multiple fate-specific genetic programs.

Multipotent neural crest cells are a migratory embryonic cell population found only in vertebrates that confers traits such as complex head, teeth, elaborate endocrine regulation, and autonomic and sensory nervous systems. Neural crest cells have distinct developmental potential along the anterior-posterior axis. Cranial neural crest cells give rise to mesenchymal cell types of the head, including facial cartilage and bone, glia, and some neurons of cranial ganglia and pigment cells (4). By contrast, trunk neural crest cells do not form bone or cartilage derivatives in vivo, even after being grafted to the head (5). However, despite knowledge about genes and signals that regulate neural crest development (6), the mechanisms that enable neural crest cells to commit to a multitude of possible cell fates at each axial level remain unclear.

Here we combine single-cell RNA sequencing (scRNA-seq) with spatial transcriptomics and lineage tracing to examine the spatiotemporal transcriptional landscape and cell fate decisions involved in cranial and trunk neural crest differentiation in the mouse embryo at the embryonic day 8.5 (E8.5) to E10.5 stages.

scRNA-seq reveals the transcriptional changes in neural crest during delamination

Wnt1Cre recombines in the dorsal neural tube where the premigratory neural crest resides. Focusing first on trunk neural crest, we dissected, then dissociated, the cervical region and trunk areas posterior to the otic vesicle from E9.5 Wnt1Cre/R26RTomato mouse embryos (Fig. 1A), and measured mRNAs of single TOMATO+ cells with high coverage using Smart-seq2 protocol (median of 7025 genes detected per cell) (fig. S1A). Our design takes advantage of the developmental asynchrony of neural crest formation along the anteroposterior axis of the embryo to sample cellular states along neural crest maturation trajectories.

Fig. 1 Major aspects of trunk neural crest heterogeneity.

(A) Immunohistochemistry of E9.5 Wnt1Cre;R26RTOM+ embryo, showing SOX10+/Wnt1TOM+ neural crest cells (NCCs) migrating in the head and trunk. Note post-otic vesicle incision (dashed lines), showing separation between cranial and trunk portions. (B and C) t-SNE embedding shows 1107 cells (points) from mouse E9.5 Wnt1Cre/R26RTomato trunk. The embedding reflects spatiotemporal aspects of neural crest (NC) development. Twelve major clusters of transcriptionally similar cells (colors) correspond to different stages of NC and ventral neural tube (vNT) development. (D) Analysis of RNA velocity shows major directions of cell progression in transcriptional space. The arrow start and end points indicate observed-current and predicted-future cell states, respectively. (E) Gene markers (cluster-specific) link NC clusters to developmental states. (F) Anatomical localization of NC clusters inferred from in situ sequencing of 32 cluster-specific genes in E9.5 sections (representative trunk section). Cell identity of vNT, which was not captured in the single-cell dataset, is imputed using the observed NC and NT clusters. (G) Characteristic gene expression patterns as detected by scRNA-seq (upper) and in situ sequencing (bottom) of a representative trunk section (see data S11 for others). (H) Transcriptional dynamics around delamination illustrate relationships between previously unknown and known markers. Cells of pre-EMT, delaminatory, and early migratory clusters were ordered by their delamination pseudotime, estimated using a principal curve. (I) scRNA-seq data predicts specific Dlx5 expression around delamination (left), validated through immunochemistry (right).

Analysis of transcriptional heterogeneity (7) separates neural tube and neural crest populations into two compartments (6) (Fig. 1, B, C, and E, and fig. S1, F and G). These are connected by two “bridges,” corresponding to neural crest delamination [marked by activation of epithelial-to-mesenchymal transition (EMT) drivers such as Snai1] (6) (Fig. 1, G to I, and fig. S1C) and neurogenesis processes [marked by proneurogenic transcription factors (TFs)] (fig. S2, B to E). Estimates of RNA velocity, which predict the direction in which cells are moving in transcriptional space (8), show a convergent velocity pattern in the neurogenesis bridge, reflecting convergence of neural crest and neural tube neurogenesis programs (Fig. 1D). By contrast, the delaminating bridge shows pronounced movement from the neural tube toward differentiating neural crest cells. Although delamination of neural crest has been viewed as an abrupt transition of pre-EMT to migrating neural crest (9), our data reveal an extensive sequence of transcriptional events that unfold during delamination and early migration (Fig. 1, B, H, and I; fig. S1C; and table S2).

This enabled us to separate the premigratory neural crest into two distinct subpopulations. The earliest, pre-EMT population, is composed of cells that have not yet started delaminating from the neural tube. It is marked by peak expression of neural plate border specifiers previously tied to NC, such as Zic1/3/5, Msx1 (9, 10), Mafb (11), and Gdf7 (12) (Fig. 1, E and G, and fig. S1C). The pre-EMT population, however, still expresses neural tube markers such as Olig3 and FoxB1 and is localized to the neural tube (Fig. 1H and fig. S1, C and F). The second, delaminating subpopulation is marked by activation of the key EMT gene Snai1 (6) and absence of Atoh1, and is accompanied by sequential transient up-regulation of a battery of genes, including Dlx5, Pak3, Pdgfra, and Hapln (Fig. 1H,I). Up-regulation of classical neural crest specifiers (Sox9, Foxd3, and Ets1) and down-regulation of many neural plate border specifiers (Zic3, Mafb, Gdf7) show variable timing across this progression, with some neuroectodermal border specifiers, such as Msx1, Msx2, Wnt3a, and Lmx1a (9, 13, 14), retaining expression in the delamination cluster. On the side opposite the neural tube, the delamination bridge connects with a subpopulation of Sox9+ cells that express Cdh11 or Itga4, genes involved in neural crest migration (15, 16), but do not yet express known markers of downstream neural crest fates (Fig. 1E). We will refer to these populations as migrating progenitors. Overall, single-cell data identify substages within the pre- and postmigratory neural crest stages (6, 9) and yield sets of marker genes for various stages (table S1).

In situ sequencing shows spatial segregation of distinct neural crest states

Transcriptional profiles reveal a cascade of states from neural crest formation to fate commitment (Fig. 1, A to E). The early pre-EMT and delaminating crest gives rise to a pool of migrating progenitors, whose heterogeneity is linked to transcriptional properties of the downstream fates (fig. S1D), as they begin to express fate-specific genes. RNA velocity (Fig. 1D) shows a general flow of neural crest cells along the sequence of these developmental stages, followed by divergent flows of migrating progenitors toward more mature neural crest derivatives, including autonomic and sensory neuronal lineage and mesenchyme (Fig. 1, H and I, and data S1). All subpopulations demonstrated robust and uniform cell cycle signatures, except for a small number of cells committed to neuronal lineages (fig. S4, A to D). Immunohistochemistry-validated markers of major subpopulations (fig. S1, F to I) and RNAscope confirm migratory spatial patterns of classical lineage markers and new genes predicted to be implicated in fate specification (fig. S2 and data S12) (17, 18).

To establish the relative anatomical distribution of the identified subpopulations, we used in situ sequencing (19) to simultaneously detect 32 subpopulation-specific transcripts in a spatially resolved manner in 15 serial sections of the E9.5 embryo (Fig. 1G and data S1 and S11). These two-dimensional measurements confirmed the expected anatomical separation of the neural crest and neural tube populations, and a dorsal neural tube localization of pre-EMT cluster, as well as spatial segregation of mature fates along the dorsoventral axis (Fig. 1, F and G, and fig. S1E). The subpopulations of migratory progenitors showed dorsoventral spatial segregation, with progenitors transcriptionally adjacent to the sensory fate (yellow, Fig. 1B) biased toward dorsal regions, and the progenitors transcriptionally adjacent to autonomic and mesenchymal fates (violet, Fig. 1B) biased toward ventral regions (Fig. 1F and fig. S1E).

Differentiating neural crest undergoes sequential binary fate restrictions

To disentangle the transcriptional logic of neural crest differentiation in the trunk, we used a branching process to computationally model the progression of emigrating neural crest cells through the heterogeneous pool of progenitors toward committed fates. Specifically, we adapted a principal graph approach (20) to derive a statistical ensemble of contiguous tree trajectories that best explain the observed cell distribution in a high-dimensional transcriptional space (Fig. 2A and fig. S3, A to D; see methods). The resulting trees capture transcriptional changes associated with transitions between pre-EMT, delaminating, and migratory states, followed by multiple cell fate decision branches (figs. S3, A to C, and S4, H to J, and table S3). The developmental potential of neural crest cells has been the subject of previous investigations, showing that pre-EMT and early migrating neural crest cells are multipotent (21, 22). However, it remains unclear how neural crest cells lose multipotency and whether the choice of multiple fates occurs in a stochastic manner or follows a structured pattern of lineage restrictions (23, 24). Our results demonstrate a well-defined transcriptional structure of cell-fate splits, which appear as a sequence of binary bifurcations separated by additional robust transcriptional changes (Fig. 2A and figs. S3, A to C, and S4, H to J). The first stable bifurcation separates sensory lineage from common progenitors of autonomic and mesenchymal branches. The second stable fate split, separating autonomic neuronal fate from mesenchymal differentiation, captures a spatially restricted process that takes place within the cervical region. Additional branches can be attributed to glial differentiation as they show expression of early glial markers (Mpz, Fabp7, Zfp488, Plp1, Sox10) and transcriptional signatures of E10.5 Schwann cell precursors (SCPs) (fig. S4, E to G). Furthermore, removal of the top 50 cells associated with SCP development leads to disappearance of these glia-specific branches (fig. S4, H and I). Although the topology relating sensory, autonomic, and mesenchymal branches is stable, the anchoring of the glial branch is uncertain, with some of the trees in the ensemble placing it within the sensory branch, and others positioning it within the autonomic-mesenchyme branch (fig. S4J). This likely reflects both technical uncertainties arising from limited number of such glial precursors, as well as the fact that glial precursors arise in both sensory and autonomic lineages.

Fig. 2 Developmental portrait of trunk neural crest cells.

(A) Developmental progression of trunk NCCs modeled using principal trees. Non-NCCs (empty circles) were excluded from the analysis. Cells with high average expression of fate-specific markers are shown by distinct vibrant colors; a heterogeneous pool of migrating progenitors is shown in gray. (B) Projections of cells onto the principal tree yield smoothed pseudotime-associated gene expression profiles. Top: Color-coded segments of the principal tree; segment colors match the colors of corresponding unbiased clusters in Fig. 1A, where possible. (C) Clustering of genes based on similarity of NCC expression profiles (see data S1 for extended set of clusters). (D) Transcription factor (TF) regulatory activity along the NCC developmental progression. Lower: A pseudotime activity pattern is shown for a set of TFs predicted to be active during NCC differentiation (FDR < 0.25). (E) Patterns of expression (left) and activity (right) shown for representative TFs. (F and G) Immunohistochemistry of control Neurog2CreERT2/+;R26RYFP/+ (carrying one copy of CreERT2) and mutant Neurog2CreERT2/CreERT2;R26RYFP/+ (carrying two copies of CreERT2, see methods for details) tamoxifen-injected at E9.5 and analyzed at E15.5 shows fate partitioning of traced cells between sensory neurons (ISL1+) (white solid arrowheads), glia (SOX10+) (empty arrowheads), and melanocytes (MITF+) (white solid arrows). Note that no melanocytes are traced in control embryos. (H) Quantification of the percentage of traced sensory neurons, glia, and melanocytes between control and mutant embryos per organ shows doubling of tracing in neurons and glia, but significantly higher proportion in melanocytes (*p < 0.05, **p < 0.01, ***p < 0.001). DRG, dorsal root ganglion; SG, sympathetic ganglion.

Because it has been previously noted that Wnt1Cre line can activate the Wnt1 signaling pathway and induce a developmental phenotype (25), we generated another single-cell snapshot of the E9.5 trunk neural crest, using a different neural crest-specific Cre line (Sox10CreERT2/R26RTomato). Transcriptional dynamics, structure of bifurcations, and patterns of markers expression were similar between Wnt1Cre and Sox10CreERT2 snapshots, indicating that the effect of Wnt1 activation in early neural crest is relatively minor (fig. S5, A to I), consistent with observations by others (26). To evaluate the extent of spatiotemporal conservation of transcriptional program in posterior neural crest, we also generated a single-cell snapshot of hindlimb and tail crest at E10.5 using Wnt1Cre/R26RTomato line. The results recapitulate transcriptional dynamics at E9.5 in trunk, revealing more advanced stages of sensory neurogenesis (fig. S5, H and I). These combined results show that upon maturation, trunk NC proceeds through a stereotypical sequence of binary lineage-restriction decisions.

Distinct functions of early and late Neurog2 expression

Many genes exhibit statistically significant and robust changes along the reconstructed trunk neural crest developmental tree [1048 genes at a false discovery rate (FDR) of 0.05; Fig. 2, B and C, and table S4], which we grouped into 21 major clusters based on the similarity of their expression patterns (Fig. 2C and data S2). Although expression of many clusters recapitulated the canonical stages of neural crest (6), others captured more complex regulatory patterns, such as down-regulation after delamination followed by reactivation upon mesenchymal specification (cluster 20, data S2). Among the tree-associated genes, we identified 137 TFs, many of which were not previously described in the context of neural crest development (e.g., Maf, Ikzf2, Rfx4, Ldb2, Plagl1, Nhlh2) (data S2). Differential expression of a TF, however, does not establish its regulatory role, as TF activity may be modulated by cofactors or other regulatory machinery. We therefore looked for coordinated up-regulation or down-regulation of the predicted TF targets as an indicator of TF regulatory activity (Fig. 2D). We used regularized linear models (27) to analyze 50 of out 137 trajectory-associated TFs that had known nonredundant sequence binding sequence motifs, and identified TFs whose transcriptional changes show statistically significant regulatory impact on their corresponding target genes (FDR < 0.25) (Fig. 2D and fig. S6). A pattern of TF activity matching TF expression was observed for transcriptional activators, and a complementary pattern for repressors (e.g., Foxd3). A repressive effect of Foxd3 in non-ectomesenchymal developmental stages is in agreement with the results of FoxD3 loss-of-function experiments (28), including a biasing role of FoxD3 against the mesenchymal program (29). Several TFs showed expression in parts of the neural crest development tree without exhibiting a detectable regulatory impact on their target genes, suggesting modulation by other processes. For example, Insm1 is expressed in both autonomic and sensory branches, but shows detectable regulatory impact only in the autonomic branch. Indeed, activity of Insm1 is modulated by autonomic-specific factor Ascl1 in specification of autonomic lineage (30).

Proneurogenic Neurog2 exhibits two peaks of expression (Fig. 2B): early after the delamination stage that appears to lack direct regulatory impact (Fig. 2E), and late at the onset of sensory neurogenesis that can be linked to corresponding regulatory activity (31). Expression of Neurog2 around the dorsal neural tube has been previously assumed to mark sensory progenitors (32), and in that regard, the early expression of Neurog2 that precedes any bifurcation points is surprising. It suggests that instead of being limited to sensory lineage, all neural crest derivatives may exhibit transient expression of Neurog2 during differentiation (Fig. 2E). Indeed, lineage tracing in Neurog2CreERT2/+;R26RYFP mouse line starting from E9.5 confirms that cells expressing Neurog2 at this early stage end up being widely distributed across all types of the neural crest progeny in the trunk and are not limited to the canonical sensory fate (Fig. 2F). Earlier studies have shown that Neurog2 knockout results in transient delay of sensory lineage differentiation, which is subsequently compensated by its homolog Neurog1 (31) that also has detectable regulatory impact in the sensory branch (fig. S6A). To assess broader functional impact of early Neurog2 expression on fates distribution, we analyzed knockin Neurog2CreERT2/CreERT2;R26RYFP embryos (33). In the absence of Neurog2, the efficiency of genetic tracing appeared proportionally higher in all neural crest derivatives except melanocytes, indicating that early Neurog2 expression is not primarily associated with committed or biased sensory progenitors. Rather, it suggests a nonneurogenic role of Neurog2 in active repression of the melanocyte fate during early neural crest migration (Fig. 2, G and H).

Coactivation of alternative programs and biasing precede fate commitment

To investigate the process of cell fate commitment, we examined how transcriptional modules associated with alternative cell fates assemble and compete around decision points, focusing initially on the first bifurcation separating sensory and autonomic branches (Fig. 3A). We approximated fate-specific modules by sets of genes selectively up-regulated in one branch compared to another (Fig. 3, B and C). Whereas some fate-specific genes are up-regulated after the bifurcation point (late genes, 45 sensory, 36 autonomic), others show earlier up-regulation in the progenitor branch, before the actual bifurcation point (early genes, 53 sensory, 86 autonomic) (Fig. 3B and data files S3 and S4). For example, Neurod1, thought to be a master regulator of sensory fate, is induced late in the sensory branch, whereas another sensory TF, Pou4f1, is up-regulated earlier during delamination, in a manner similar to that of Neurog2 (Fig. 3B). This is expected, as bifurcation points capture trajectory positions where the extent of the transcriptional differences between the alternative cell fates becomes sufficiently large. The expression onset differences were also confirmed by the RNAscope measurements of select genes in both sensory and autonomic branches (fig. S2). We observed that early modules of competing cellular fates are coexpressed in progenitor cells and exhibit gradual coactivation as the cells progress toward the bifurcation point, followed by selective up-regulation of one module and down-regulation of another after the bifurcation point (Fig. 3C). By contrast, late-competing modules show mutually exclusive activation in their corresponding branches, without coexpression in individual cells, consistent with commitment to a particular fate (Fig. 3D). Many known critical master regulators do not show detectable expression around the bifurcation point (including autonomic Phox2b, Ascl1, and sensory Neurog1, Neurod4, and Neurod1), indicating contribution of other genes or mechanisms to the early steps of the decision process. These likely include exposure to extracellular signals such as Wnt signaling, BMP signaling, or Delta-Notch interactions (34).

Fig. 3 Logic of cell fate selection for sensory-autonomic bifurcation.

(A) Schematic overview of the sensory-autonomic split in NC differentiation. (B) Branch-specific (sensory versus autonomic) genes classified as activated early (before split) and late (after split). Right: Example of an early (Pou4f1) and a late (Neurod1) sensory-specific gene. Vertical line marks predicted time of up-regulation. Time of up-regulation of all genes can be found in data S3 and S4. Colors encode branches [as in (A)]. (C) Average expression pattern of early and late branch-specific modules. Early and late modules are set of genes activated before or after the bifurcation. (D) Scatter plots show average expression of sensory and autonomic modules in each cell. Left: Early competing modules show gradual coactivation, followed by selective up-regulation of one fate-specific module and down-regulation of another. Right: Late modules show almost mutually exclusive up-regulation. Colors encode tree branches. (E) Bifurcation point is characterized by correlated expression of genes within modules, and repulsion of competing modules. Left: Probability of a cell proximity to the sensory-autonomic bifurcation point shown (darker colors corresponding to higher probability). Right: Average correlation of each early sensory- or autonomic-specific gene (points) with the genes from the sensory- (x axis) and autonomic-specific (y axis) modules among cells around the bifurcation point. Genes that have low average correlation with their own module (<0.07) are shown with faded colors and are excluded from analysis in (F) for clarity. (F) Early modules show gradual intramodule coordination and intermodule repulsion. Upper: The plots show average local correlations of module genes with autonomic and sensory modules in a set of cells with similar developmental time (marked by black points). Bottom: Average local correlations of genes within and between branch-specific modules show gradual increase in coordination within each module and increasing antagonistic expression between the branch-specific modules. The difference between intra- and intermodule correlations, which characterize the extent of the antagonistic expression between the alternate modules, is shown in the upper right corner of tSNEs. (G) Time of gene inclusion in coordinated branch-specific module shows early signatures of fate-biased expression programs. Bottom: Schematic illustrating of formation of a branch-specific module, with connecting lines indicating onset of expression correlation between genes. Heatmap: Probabilities of gene inclusion in the branch-specific module during early stages of differentiation (see data S5 and S6 for gene names). Right: Examples of inclusion probability estimates for individual genes in the two modules. (H) Cell fate decisions formalized as a three-phase process: Initially, cells coactivate competing cell fate–specific programs, gradually switch to preferential expression of one module, and culminate by activating mutually exclusive fate-specific expression.

We next examined whether despite general coactivation of competing-fate modules, the progenitor cells exhibit an early transcriptional bias toward a particular fate beyond the effects of expression noise. In doing so, we controlled for apparent biases due to developmental asynchrony of the subpopulations being analyzed, which can limit interpretation of similar analysis based on targeted or bulk measurements (fig. S7I) (2). For the cells positioned around the actual bifurcation point, the transcriptional bias toward one of the possible fates manifests itself as a pronounced negative correlation between the competing gene modules (Fig. 3E). We find that as cells move toward the first bifurcation point following the initial coactivation stage, the degree of transcriptional coordination within each module increases, and antagonistic expression of the competing modules becomes more apparent (Fig. 3, E and F). Simulated controls, with expression profiles randomized across cells with similar pseudotime to preserve the overall gene activation profiles, fail to show such local coordination of gene modules (fig. S7, A and B). This indicates that as cells progress toward the bifurcation point, they undergo gradual transcriptional biasing toward one of the competing fates. The biasing phase does not appear to be driven by up-regulation of a single gene. Rather, we observe broad expression shifts between genes of the competing fate-specific modules (Fig. 3F). These results indicate the presence of lineage priming early in neural crest migration, prior to the predicted sensory-autonomic bifurcation point. The transcriptional correlations within the competing modules are reduced following the bifurcation point. This reduction illustrates that most of the observed intramodule correlation can be attributed to compositional heterogeneity of cells around the bifurcation point (35). When considering more homogeneous cell populations within the branches following the bifurcation point, correlation of genes within the module becomes relatively low, indicating that regulatory interactions between genes, such as correlations induced by stochastic covariation of a TF and its target genes, are not apparent in such data (36) (Fig. 3F).

To determine the initial transcriptional events in the assembly of the fate-specific modules, we monitored coordinated expression of gene subsets, tracking backward along the progenitor branch. In doing so, we controlled for the overall pattern of module activation, as such general trends can overshadow correlations within a subset (fig. S8). The results show that early correlated subsets of modules specific to either sensory or autonomic fate can be statistically distinguished at the delamination stage, but not in the pre-EMT neural crest (Fig. 3G and data files S5 and S6). We detect earlier formation of the sensory module, compared to the autonomic module. Overall, we observe three primary phases of cell fate decision, formalized as a sequence of initial coactivation, gradual biasing, and commitment phases (Fig. 3H). The initial coactivation stage is characterized by coexpression of competing modules within individual cells that gradually shifts in favor of one of the modules during the biasing phase, culminating in up-regulation of mutually exclusive fate-specific gene expression patterns during commitment. The three-phase scheme of cell fate decision also holds for the autonomic-mesenchyme bifurcation discussed below (fig. S7, C to H). This model elaborates on earlier observations of coexpression of key regulatory factors (3740) or programs (4143) from alternative fates prior to cell fate commitment noted in other tissues.

Autonomic-mesenchymal bifurcation reveals bipotent progenitors

Studies of neuroblastoma, a tumor of sympathetic nervous system that originates from neural crest cells, revealed that tumor cells can transit between adrenergic and mesenchymal phenotypes, creating inter- and intra-tumoral heterogeneity that poses a therapeutic challenge (44, 45). Our data uncovered the developmental bifurcation that separates autonomic and mesenchymal fates (Fig. 4A). After the initial stage of co-activation of both programs (Fig. S7C-E), commitment to fates is marked by expression of two fate-specific late genes: Phox2b, a key driver of autonomic neuronal lineage (46); and Prrx1, a marker of specified mesenchyme (47). However, some transient cells co-express both genes at high levels (Fig. 4A). Consistent with this, corresponding Prrx1 and Phox2b mRNA molecules colocalize and proteins are coexpressed in the same few cells in the ventral neural crest pathway of migration (Fig. 4, C and D). Nearly all late genes of the autonomic module are coexpressed with late genes of the mesenchymal module in at least a few cells, illustrating ubiquitous stochastic coactivation of the alternative fate modules (Fig. 4B).

Fig. 4 Phox2b and Prrx1 coexpression reveals bipotent progenitors.

(A) Infrequent coexpression of Phox2b and Prrx1 is observed despite mutually exclusive expression of modules. Bottom left: Late modules of the autonomic-mesenchyme bifurcation show mutually exclusive expression. Bottom right: Prrx1 and Phox2b show largely branch-specific expression pattern, but are coexpressed in some cells. (B) The majority of antagonistic pairs of branch-specifying genes are coexpressed at high levels in a few cells. Upper: Heatmap with color indicating numbers of cells with coexpression of two genes from competing late modules (rows, autonomic; columns, mesenchymal). Lower: Heatmap showing expression of branch-specific TFs in cells assigned to either branch. (C) Pattern of Prrx1 and Phox2b expression (by in situ sequencing) at E9.5. Coexpression of the two genes is observed at the prospective cranial ganglia. (D) Immunofluorescence for PHOX2B, PRX1, and SOX10 in E9.5 trunk (PHOX2B+/PRX1+ cells shown by white arrowheads, PHOX2B+/PRX1 by yellow arrows). (E) Immunohistochemistry on E17.5 embryo for Prrx1GFP, S100, and ISL1 shows absence of overlap between the mesenchyme (Prrx1GFP+), glia, and sensory neurons. Note that Prrx1GFP+ shows colocalization with COLIV, indicating that the vascular mesenchyme was traced. (F) Immunohistochemistry on E13.5 embryo for SOX10, Phox2bTOM+, neurofilaments, and PHOX2B (in heart-related panels). Note the presence of numerous Phox2bTOM+/ PHOX2B mesenchymal cells in the cervical region and surrounding Phox2bTOM+/ PHOX2B+ autonomic neurons in the heart. (G) Immunohistochemistry on E13.5 embryo for SOX10 or PHOX2B, Phox2bTOM+, and ISL1. Note the presence of numerous Phox2bTOM+/ mesenchymal cells in aortic areas, negative for ISL1, PHOX2B, or SOX10. (H) Analysis of the percentage contribution of Phox2bTOM+ and Prrx1GFP+ cells to mesenchyme, autonomic neurons, and glia over the total number of traced cells in the area of the heart (N = 3 embryos per strain).

To test whether the cells expressing Phox2b can end up in the mesenchymal domain, we performed lineage tracing with a Phox2bCre/R26RTOMATO mouse line (Fig. 4, F and G). As expected, the autonomic neurons and numerous glial cells proximal to autonomic components were traced. However, we also observed tracing of multiple cells of the mesenchymal phenotype in the cervical region, as well as in the heart, in the proximity of autonomic components (parasympathetic neurons of the heart) (Fig. 4H). Detailed analysis at the region of the heart and the proximal autonomic ganglia showed that 21.6 ± 3.1% of traced cells contributed to glia, 9.9 ± 1.6% to autonomic neurons, and 68.4 ± 4.7% to mesenchyme (Fig. 4H). These results confirm that some of the cells expressing late-autonomic signatures subsequently deviate toward mesenchymal specification in the cervical region. The opposite experiment with lineage tracing in Prrx1-Cre;R26RmTmG showed no traced autonomic neurons or glial cells in the entire embryo, including the cervical region where the large superior cervical ganglion is positioned (Fig. 4, E and H). Thus, cells expressing Prrx1 always differentiate into mesenchymal fates. As a result, when both Prrx1 and Phox2b master regulators are coactivated in the same cells, cell fate is resolved in the direction of mesenchymal derivates in vivo (Fig. 4B). No sensory neurons were traced from progenitors expressing Phox2b, demonstrating that Phox2b+ cells at the second bifurcation point are bipotent. The possibility of transdifferentiation of more committed Ascl1+/Phox2b+ autonomic progenitors into mesenchymal population is not supported by RNA velocity analysis (Fig. 1D) and lineage tracing from E10.5 with Ascl1CreERT2;R26RTOM (fig. S9).

Transcriptional specificity of cranial neural crest emerges during delamination

Migration and differentiation of cranial neural crest are spatiotemporally separated from events occurring in the trunk. The cranial crest primarily gives rise to mesenchymal populations responsible for building the head, which are not produced by the trunk crest (48). We analyzed single-cell snapshots of anterior cranial, pre-otic, and post-otic Wnt1-traced populations at E8.5 (Fig. 5A). Early neural crest was composed of three spatially distinct populations: Hox population that corresponds to the anterior cranial neural crest, Hoxb2+ that contributes primarily to the mandibular level, and HoxD3+ population that includes post-otic (including cardiac and vagal) streams of the neural crest (Fig. 5A, fig. S10A, table S5, and data S7) (49). The developmental portrait of Hox anterior cranial neural crest at E8.5 represented a single differentiation trajectory that encompasses diverse transcriptional events: down-regulation of a neural tube–associated program, including Zic3 and Pax8; transient expression of a battery of genes such as Pak3 and Bmper; coherent, but not simultaneous, up-regulation of neural crest markers, including Foxd3, Sox9 and Sox10, and mesenchymal fate–specifying genes, including the TFs Twist1 and Prrx2 (Fig. 5, B and C, table S6, and data S8). Consistent with this, analysis of Hox- neural crest variability shows trajectory-associated and cell cycle processes (fig. S10G). Although we did not observe fate bifurcations in the cranial crest at the E8.5 stage, an anticorrelation of the mesenchymal and neuronal genetic programs emerged, coinciding with the down-regulation of the early neural crest genes Foxd3 and Sox9 and activation of the mesenchymal factors (Fig. 5D and fig. S10B). Thus, cranial neural crest, similar to trunk cells, progresses through a biasing stage before commitment to a specific fate. This indicates that signs of ongoing cell fate bifurcations often can be detected in advance, similar to early detection of transitions noted in other fields (50). As in the case of the anterior cranial neural crest, the snapshot of Hoxb2+ neural crest cells underwent a differentiation progression through an intermediate step toward migrating cells expressing mesenchymal markers (fig. S10, C and D; Table S7; and data S9). Hoxb2/Hoxd3+ neural crest population was represented as an isolated cluster.

Fig. 5 Spatiotemporal analysis of cranial neural crest.

(A) t-SNE embedding shows subpopulations observed in 1345 cells collected from the cranial part of E8.5 Wnt1-Cre;R26TOM+ embryos (left: representative embryo shown). (B) Analysis of the developmental progression identifies timing of transcriptional events of early Hox cranial NCCs. (C) Pseudotime expression trends of representative genes [from (B)] are shown for early Hox cranial NCCs. (D) Correlation of fate-specific gene modules reveals emergence of mesenchymal versus neuronal biases in migration. Similar to Fig. 3F, the cells were arranged by pseudotime and binned in groups of 50 cells (upper panels) to estimate intra- and intermodule correlations between the mesenchymal and neuronal (sensory, autonomic) gene modules from E9.5 trunk NCC analysis (lower panels). (E) Transcriptional differences of migratory Hoxb2+ and Hox cranial NCCs emerge during delamination. Left: NCCs from Hox- and Hoxb2+ populations show more similar profiles in the migratory stage, compared to NT, whereas transcriptional differences between Hox- and Hoxb2+ in NT (ΔNT) and migratory stage (ΔNC) are not correlated. Right: The difference in transcriptional shifts accompanying delamination in Hoxb2+ (dNC+) and Hox- (dNC) NCCs correlates with the resulting expression difference between Hox and Hoxb2+ populations at the migratory stage (ΔNC), indicating that most migratory differences between levels emerge during delamination, instead of reflecting initial differences at the NT stage. (F) A joint embedding of NCCs from different stages (E8.5, E9.5, E10.5) and locations (cranial, trunk, cardiac, hindlimb, and tail) shows stereotypic landscape of NC development. Inset shows schematic tree of differentiation. (G) Distribution of transcriptional states of spatiotemporally distinct datasets (colored according to developmental time). (H) Three early ectomesenchyme developmental states of progressive maturation in cardiac NC population.

Cranial neural crest at all spatial levels expresses a shared core program (including classical markers Sox10, Sox9, Foxd3, and Ets1). As a result, Hox and Hoxb2+ neural crest subpopulations are more similar to each other than to the corresponding neural tube regions from which they originate (Fig. 5E). Nevertheless, Hox and Hoxb2+ neural crest shows expression differences (fig. S10, E and F), which may be a consequence of divergent differentiation trajectories, or simply reflect residual differences of their distinct starting states prior to delamination from the neural tube. We find evidence for the former, with the difference between Hox and Hoxb2+ neural crest correlating well with the differences in the neural crest developmental programs induced at the corresponding spatial levels (r = 0.49), but not correlated with the difference of their starting Hox and Hoxb2+ neural tube states (Fig. 5E). This indicates that the transcriptional divergence of the neural crest from different spatial levels arises primarily during downstream delamination and migration stages through the activation of new gene expression programs.

Cranial and trunk crest follow different paths through a shared differentiation landscape

To compare more-differentiated populations arising from different neural crest levels, we further sampled crest-derived cells of cranial region at E9.5, post-otic (vagal and cardiac) at E10.5, and lumbar and tail regions at E10.5. Joint alignment of all crest datasets (see methods) revealed that neural crest cells at various positional levels and time points navigate a stereotypic space of transcriptional states, with some notable biases (Fig. 5F; fig. S10, H and I; and tables S9 and S10). The distribution of mature cells of the cranial compartment is biased toward mesenchymal fate, whereas cells of the trunk compartments gravitate toward sensory and autonomic fates (Fig. 5G and fig. S10M). In line with this, the fraction of mesenchymal committed cells was high in the anterior population, whereas the fraction of sensory committed cells was high in the posterior population as inferred from Hox gene expression (fig. S10J). The cranial neural crest progressed to a more mature mesenchymal state at E9.5 as compared to the earlier E8.5 time point, as evident from activation of Prrx1 and Twist2. By contrast, trunk neural crest showed active maturation of sensory progenitors from E9.5 to E10.5, including activation of Dcx, Stmn2, and down-regulation of Neurog2 (fig. S10I). Post-otic neural crest occupied an intermediate position between anterior cranial and trunk neural crest, and shared some features with each of those populations (51). Snapshot of post-otic neural crest populations at E10.5 represented mostly Hoxd3+ cells, indicating a dominant contribution from the vagal neural crest (Fig. 5G). The population contains cells of autonomic and cardiac-mesenchyme fates with a small admixture of sensory-committed cells. Cardiac mesenchyme is necessary for development of cardiac valves and the outflow tract (52). Heterogeneity of this subpopulation revealed three initial steps of ectomesenchyme differentiation that are marked by sequential acquisition and down-regulation of Six2 followed by transient up-regulation of Dlx6 and activation of Msx2 and heart-specific Hand2, and then activation of downstream cardiac markers Hand1, Dkk1, and Gata6 (Fig. 5H, table S8, and data S10). Screening of the neural crest snapshot at E8.5, E9.5, and E10.5 stages did not reveal coexpression of the melanoblast markers (Mitf, Pmel, Dct) or a broader transcriptional signature of prospective melanoblasts at any of the spatiotemporal levels, indicating late consolidation of melanocytic fate (fig. S10, K and L).

During delamination and early migration stages, both cranial and trunk neural crest progress through a similar sequence of transcriptional steps, as illustrated by intermixing in joint expression space (Fig. 5G and fig. S11A). Nevertheless, in addition to such common maturation signatures, cranial and trunk neural crest also activate their own distinct modules (Fig. 6A and fig. S11, A and D). For cranial cells, this includes activation of a broad mesenchymal program, including regulatory factors such as Twist1, Prrx2, and Dlx2 (Fig. 6B), whereas activation of the sensory program that includes Neurog2, Pou4f1, and Nkx2.1 is specific to trunk cells. Therefore, we hypothesized that early activation of fate-specific programs after delamination can predispose cells toward specific fates and account for the distinct downstream fate biases of trunk and cranial neural crest cells.

Fig. 6 Activation of early mesenchymal and sensory transcriptional programs.

(A) Scatter plot contrasts expression changes accompanying NC transition from pre-EMT to migrating stages in cranial and trunk NC shows overall correlation (r = 0.5). However, more up-regulated genes upon transition in trunk or in cranial compartments end up preferentially expressed in sensory and mesenchymal-autonomic fates, respectively. (B) Mesenchymal regulators are activated early in cranial cells (top, only cranial expression shown), whereas sensory regulators are activated early in trunk cells (bottom, only trunk expression shown). (C) Immunofluorescence for SOX2, ISL1, and RFP (red fluorescent protein) after electroporation of control RFP-containing plasmid (left) versus after Twist1 electroporation (right) in the developing chicken embryo trunk shows ectopic mesenchymal induction upon Twist1 overexpression. NC, neural crest; NT, neural tube; NCC, neural crest cell; DRG, dorsal root ganglion.

Twist1 activation is sufficient to reroute trunk neural crest into a mesenchymal fate

To test if activation of mesenchymal regulators can bias trunk neural crest toward mesenchymal differentiation, we electroporated Twist1-expression constructs into the developing chick neural crest. Such targeted and sustained overexpression of Twist1 in the trunk neural crest starting from pre-EMT stage resulted in down-regulation of traditional trunk fates such as neuronal sensory, autonomic, and glial and alteration of the neural crest migration patterns (Fig. 6C). Instead of the expected migration toward the dorsal aorta, Twist1+ cells in the trunk predominantly followed a path to dermal regions, where they activated expression of a dermal mesenchymal marker Prrx1, not observed in wild-type trunk neural crest (Fig. 6C). Some Twist1+ cells migrated to sensory ganglia and the dorsal aorta, where they did not take neuronal or glial fates (fig. S12, A and B). Conversely, CRISPR-Cas9–based knockdown of Twist1 in chick cranial neural crest resulted in a reduction of mesenchymal populations and an increase in glial and neuronal derivatives (fig. S12, C to F), consistent with previous studies (53, 54).

Overexpression of the downstream mesenchymal factor Prrx2 in the trunk resulted in massive apoptosis of neural crest cells (fig. S11B). The rare surviving cells were not able to pick traditional trunk fates, as compared to controls [green fluorescent protein (GFP)–only DNA or neuro-glia–related genes, such as Crabp1; fig. S11B]. By contrast, regulator of sensory fate Pou4f1 is strongly activated upon delamination in the trunk, but is poorly expressed in cranial crest: Overexpression of the dominant-negative version of Brn3c, a homolog of Pou4f1 in chick, in cranial and trunk chick neural crest resulted in no obvious phenotype in the trunk, as expected, whereas it caused a decrease of Sox9+/Pax7+ cranial neural crest cells (fig. S11C). This suggests that Pou4f1 might be another factor differentially regulating fates in different NC populations along the anteroposterior axis.

Discussion

We used single-cell profiling to characterize spatiotemporal transcriptional dynamics of neural crest, showing progressive restrictions of cell fates, each realized through sequential phases of coactivation, biasing, and resolution of competing fates (Fig. 7A). Computational reconstruction of such differentiation trajectories from transcriptional similarity merely predicts likely lineage relationships, requiring additional validation (55). Here we relied on targeted lineage tracings of Neurog2, Phox2b, and Prrx1 to validate gradual restriction of the developmental potential predicted by the derived transcriptional tree (Fig. 7B). However, clonal tracing techniques that are being developed to progressively record clonal histories of individual cells (56) would be ideally suited for validation of a complex differentiation process such as the one that takes place in neural crest.

Fig. 7 Schematic representation of neural crest fate decisions.

(A) Compared to a classical model (left), the model of NC development proposed here resolves a sequence of stages around delamination, heterogeneity of migratory progenitors, and sequential structure of lineage restrictions. (B) Upper: Individual binary decisions are resolved via a three-phase mechanism involving coactivation, biasing, and commitment. Lower: Lineage tracing shows that early activation of Neurog2 is linked to repression of melanocytes in multipotent progenitors, whereas Prrx1 and Phox2b, involved in coactivation and biasing phases, mark downstream mesenchymal and autonomic fates with sporadic coexpression in bipotent population around fate bifurcation. (C) Distinct gene modules, activated upon delamination along the cranial to post-otic axis of the embryo, establish bias toward either mesenchymal (e.g., Twist1) or sensory neuronal fate.

The process of delamination largely erases the transcriptional signatures that distinguish different anteroposterior neural crest levels within the neural tube and activates fate-specifying gene programs. Such dissipation of pre-EMT signatures is consistent with the transient phenotype observed by Simoes-Costa et al., who identified a cranial-specific pre-EMT gene regulatory network capable of transient activation of the mesenchymal properties in trunk, but unable to overcome the environmental signals and generate trunk mesenchyme (57). It indicates that acquisition of a fate likely requires appropriate cell-intrinsic state and extrinsic signaling during delamination. Consistently, we find that mesenchymal transcriptional bias emerges during delamination, indicating that this might be the point where mesenchymal potential is imbued to the neural crest. In support of this hypothesis, we noted that sustained overexpression of Twist1, normally activated upon delamination only in the cranial compartment, is sufficient to reverse trunk crest developmental program to a mesenchymal route. This parallels previous observations in chordates, where misexpression of a tunicate Twist resulted in conversion of a cephalic melanocyte lineage into migrating ectomesenchyme (58). Twist1 might therefore be the ancient and sufficient factor defining the mesenchymal potential of migrating neural crest—a key step in the evolutionary development of the head in vertebrates.

The spatiotemporal data presented here may provide a resource for studies of neural crest biology and regulation, as well as investigations of neural crest–derived cancers and neurocristopathy-associated diseases. Although cell fate decisions in other systems may be based on different organizational principles, the molecular logic of cell fate decision characterized here for murine neural crest may also be shared by other organisms or tissues. Additional data, arising from analogous single-cell transcriptional studies or the new generalized lineage-tracing techniques, might soon allow such interorganismal comparisons to be performed.

Materials and methods summary

Neural crest cells were traced using the Wnt1TOMATO and Sox10CreERT2;R26RTOMATO strains from E8.5 to E10.5 mouse embryos. Embryos were collected and the TOMATO+ parts of the tissue were selected after stereoscopical observation using ultraviolet illumination. The tissue was dissociated and TOMATO+ cells were sorted by fluorescence-activated cell sorting on a 384-well plate in single-cell mode. The sorted cells were lysed and their transcriptomes sequenced using the Smart-seq2 protocol. Individual single-cell datasets were preprocessed using PAGODA1 and PAGODA2 packages, which includes normalization, clustering of cells, and t-distributed stochastic neighbor embedding (t-SNE) for visualization (7). Joint analysis of multiple datasets was performed using the CONOS package (59).

Individual populations of clustered cells were validated using a variety of approaches. Major subpopulations (pre-EMT, delaminating and migrating neural crest, sensory neurons, autonomic neurons, and mesenchyme) were validated on E9.5 embryos using immunofluorescence and comparison with public in situ databases (Allen Brain Atlas and Eurexpress) (17, 18). Anatomical localization of the subpopulations was established using multiplexed in situ sequencing of 32 genes (19), followed by computational alignment of the single-cell transcriptional profiles with the spatially resolved in situ measurements. Expression patterns of specific markers of the sensory lineages were validated on E9.5 and E11.5 stages using RNAscope.

RNA velocity was used to assess transcriptional dynamics of cells (8). To quantify segregation of lineages, we adapted a computational method to reconstruct transcriptional trajectory tree from single-cell profiles using principal graph approach (20). Expression of transcription factors and their potential targets, predicted on the basis of DNA binding motifs, was analyzed to assess regulatory activity along the tree. We developed a computational method to analyze patterns of cell fate decisions, which tracks formation of fate-specific gene modules prior to fate bifurcations. Dynamics around major branches of neural crest fate acquisition (mesenchyme, sensory and autonomic neurogenesis) were further analyzed using lineage tracing with the Prrx1-Cre;R26RmTmG, Phox2b-Cre;R26RTOMATO and Ascl1CreERT2;R26RTOMATO mouse strains.

Additional validation included the analysis of knockin Neurog2CreERT2;R26RYFP mice (where Neurog2 is replaced by CreERT2, constituting homozygous embryos Neurog2 knockouts). Overexpression of Twist1, as well as CRISPR-Cas9–mediated knockdown, performed as described in (60), in developing chicken embryos using in ovo electroporation, was used to validate Twist1 involvement in mesenchymal fate acquisition. Similarly, in ovo chicken electroporation was used to overexpress Prrx2 (mesenchymal driver) and Crabp1 (driver of neuron-glia fate acquisition) and a dominant-negative version of Brn3c (avian homolog of Pou4f1).

Supplementary Materials

science.sciencemag.org/content/364/6444/eaas9536/suppl/DC1

Materials and Methods

Supplementary Text

Figs. S1 to S12

Data S1 to S12

Data Tables S1 to S11

References (6175)

References and Notes

Acknowledgments: We thank the Eukaryotic Single Cell Genomics Facility for the single-cell RNA-seq and the in situ sequencing pilot facility at the Science for Life Laboratory, Sweden (for technical assistance with the in situ sequencing experiments). We also thank O. Kharchenko for help with illustrations and V. Petukhov for assistance with computational analysis of in situ sequencing data. Funding: I.A. was funded by a Swedish Research Council (Vetenskapsradet, VR) grant, ERC Consolidator grant “STEMMING-FROM-NERVE” N647844, the Paradifference Foundation, and the Bertil Hallsten Research Foundation. P.V.K. was supported by NSF-14-532 CAREER award and NIH R01HL131768. M.N. was funded by a Swedish Research Council (Vetenskapsradet) grant 2016-03645, the Knut and Alice Wallenberg Foundation, and Familjen Erling Perssons stiftelse. V.D. was supported by a Russian Science Foundation grant (18-75-10005, immunochemistry) and the Swedish Research Council (2015-03387). M.E.K. was supported by Stiftelsen Riksbankens Jubileumsfond (Erik Rönnbergs fond stipend). N.A. was supported by RSF grant 16-15-10237. Author contributions: M.K., M.E.K., J.P., L.E., N.A., Y.Y., M.H., V.D., M.F., M.L.P., F.B., C.Y., X.Q., W.-Y.H., and L.C. acquired all biological data and performed the relevant analysis. R.S., P.V.K., and T.C. performed computational analysis of single-cell data. R.S. and P.V.K. developed computational methods. R.S. and M.M.H. performed computational analysis of in situ sequencing data. C.B., M.N., M.E.B., D.A.G., J.-F.B., G.G.C., P.E., K.F., P.V.K., and I.A. gave feedback on experimental aspects, supervised experimental approaches, and implemented the data interpretation. R.S., M.K., M.E.K., and J.P. made all figures containing data and resulting analysis (except Fig. 7). R.S., M.K., M.E.K., P.V.K., and I.A. designed the study, organized experimental work, and wrote the manuscript. All authors provided feedback on figures, manuscript composition, and structure. Competing interests: M.N. holds shares in Cartana AB, a company commercializing in situ sequencing reagents. Other authors declare no conflicts of interest. Data and materials availability: All single-cell RNA-seq datasets have been deposited in the GEO under accession code GSE129114. Processed data, code, supplementary materials, and interactive views of datasets can be accessed on the authors’ website: http://pklab.med.harvard.edu/ruslan/neural.crest.html. These URLs will be maintained exactly as they are for at least 5 years with no changes. The Sox10CreERT2 strain is available from the laboratory of Vassilis Pachnis (The Francis Crick Institute, UK) under a material transfer agreement with the institution. All other data needed to evaluate the conclusions in the paper are present in the paper or the supplementary materials.
View Abstract

Stay Connected to Science

Navigate This Article