Single cell rna sequencing identifies early diversity of sensory neurons forming via bi-potential intermediates

ABSTRACT Somatic sensation is defined by the existence of a diversity of primary sensory neurons with unique biological features and response profiles to external and internal stimuli.

However, there is no coherent picture about how this diversity of cell states is transcriptionally generated. Here, we use deep single cell analysis to resolve fate splits and molecular

biasing processes during sensory neurogenesis in mice. Our results identify a complex series of successive and specific transcriptional changes in post-mitotic neurons that delineate

hierarchical regulatory states leading to the generation of the main sensory neuron classes. In addition, our analysis identifies previously undetected early gene modules expressed long

before fate determination although being clearly associated with defined sensory subtypes. Overall, the early diversity of sensory neurons is generated through successive bi-potential

intermediates in which synchronization of relevant gene modules and concurrent repression of competing fate programs precede cell fate stabilization and final commitment. SIMILAR CONTENT

BEING VIEWED BY OTHERS A DEVELOPMENTAL ATLAS OF SOMATOSENSORY DIVERSIFICATION AND MATURATION IN THE DORSAL ROOT GANGLIA BY SINGLE-CELL MASS CYTOMETRY Article 27 October 2022 SINGLE-CELL

RNA-SEQUENCING ANALYSIS OF THE DEVELOPING MOUSE INNER EAR IDENTIFIES MOLECULAR LOGIC OF AUDITORY NEURON DIVERSIFICATION Article Open access 05 July 2022 SINGLE-CELL GENOMICS OF THE MOUSE

OLFACTORY CORTEX REVEALS CONTRASTS WITH NEOCORTEX AND ANCESTRAL SIGNATURES OF CELL TYPE EVOLUTION Article 08 April 2025 INTRODUCTION In the developing nervous system, various neuronal and

glial cell types are generated from multipotent progenitor cells in a specific chronological sequence. The continuous change in the cellular potential of the progenitors is a key aspect to

this developmental process1. As neuronal progenitors exit cell cycle, they differentiate into distinct neuron types that already provide the underpinning for the identity of neuronal circuit

in which they will operate1,2. Hence, the formation of specific fates critically depends on the timing and diversity of expression of molecular factors in discrete parts of the nervous

system3,4. However, following this process, how the expression of molecular determinants is progressively acquired and restricted to particular neurons as they further diversify into

distinct neuronal classes remains unclear. A holistic understanding of the stepwise execution of elaborate transcriptional programs for neuronal diversification in heterogeneous system is

still elusive. It is therefore essential to dissect and understand the precise temporal progression of the gene-regulatory networks that produce and assemble neuronal complexity. Indeed,

such complexity is found in all parts of the developing nervous system, including the highly heterogeneous population of sensory neurons—cells that live in dorsal root ganglia (DRG) and

provide us with sensations of touch, pain, itch, temperature and position in space5. The diversity of sensory modalities emerges during early embryonic development as subtypes of sensory

neurons are generated from a common progenitor population, the neural crest cells (NCCs). Two main waves of sensory neurogenesis in DRG occur between E9.5 and E13. The first wave of

neurogenesis (E9.5−10.5)6 gives rise to myelinated neurons, the cutaneous low-threshold mechanoreceptors (here after named mechanoreceptors) and the muscle proprioceptive afferents (here

after named proprioceptors) and a small population of Aδ-fibers that can be distinguished by E11.5 based on their specific expression of transmembrane receptors and transcription factors6,7.

Mechanoreceptors express RET and MAFA and/or TRKB, proprioceptors express tropomyosin receptor kinase C (TRKC) and RUNX3, and the small Aδ-nociceptor population expresses TRKA. Later, these

populations will diversify even further, giving rise to additional subtypes representing the medium-to-large diameter DRG neurons responsible for muscle proprioceptive feedback and skin

mechanosensation modalities5. A second wave of neurogenesis (that peaks at E12 in a mouse embryo) generates the small diameter unmyelinated nociceptive neurons5,6,7, which co-express TRKA

and RUNX1 and later will diversify into subtypes necessary for pain, itch and thermal sensation8 and innervate the skin and deep tissues. Although the heterogeneity of DRG neurons in adult

mice are increasingly characterized8,9, the logic and molecular mechanisms that allow progenitors and sensory neurons to fate split and unfold specific sensory neuronal subprograms from

sensory-biased NCCs remain to be understood. Here we use deep single-cell RNA sequencing (scRNAseq) together with other experimental approaches to investigate fate choices, lineage

relationships and molecular determinants during early stages of sensory neurogenesis in the mouse. Overall, our data provide insights into the structure of the Waddingtonian landscape of the

sensory lineage and suggest that counteracting gene-regulatory networks operate within immature postmitotic neurons, resulting in the emergence of alternative-lineage transcriptional

programs defining the neuron fate choice. RESULTS DIFFERENTIATION TRAJECTORIES AND MAJOR FATE SPLITS To obtain a map of lineage splits from progenitors towards specific sensory subtypes

during mouse embryonic development, we sequenced the transcriptome of individual single cells (progenitors and neurons) of the somatosensory neuro-glial progeny of the trunk. For this, we

FACS isolated tdTomato positive (TOM+) cells isolated from mouse lines which selectively represent all NCCs (_Wnt1__Cre_;_R26__tdTOM_ and _Plp1__CreERT2_;_R26__tdTOM_) and from

neuron-specific Cre mouse lines (_Isl1__Cre_;_R26__tdTOM_ and _Ntrk3__Cre_;_R26__tdTOM_) to obtain an enrichment of somatosensory neuronal populations (Supplementary Fig. 1a, b) from E9.5,

E10.5, E11.5 and E12.5 and sequenced mRNAs of single cells with high coverage using Smart-seq2 protocol (median of 8070 genes detected per cell) (Fig. 1a and Supplementary Fig. 1c, d; see

“Methods”). The selected stages correspond to key events of sensory lineage development starting from migrating NCCs (E9.5) and finishing with the end of the second wave of sensory

neurogenesis and the early specification of neuronal subtypes (E12.5). Our dataset contains cells from _Isl1__Cre_;_R26__tdTOM_ at E9.5 and E10.5 and _Ntrk3__Cre_;_R26__tdTOM_ at E11.5 and

E12.5, both lines label sensory neurons, _Wnt1__Cre_;_R26__tdTOM_ was used to label NCCs, glial and neuronal progenitors at E9.5 and E10.5 and _Plp1__CreERT2_;_R26__tdTOM_ at E12.5. The

combination of cells collected after tracing with these Cre lines should provide a complete compendium of neuronal progenitors and early neuronal subtypes in mouse DRG (for FACS gating

strategies, see Supplementary Fig. 2). Co-embedding of the cells from all tracing strategies was made without tracing-based batch correction and lead to overlapping of cells of different

tracings following both a developmental time-wise trajectory and a progression from progenitors to specified neurons (Fig. 1a–e and Supplementary Fig. 1a, b; see below for trajectories)

which together convincingly validated the use of the dataset for further computational analysis. We used RNA velocity10, an unbiased approach that leverages the distinction of spliced and

unspliced RNA transcripts from the aligned sequences, allowing us to obtain an additional timepoint for each cells (_t_: expression of old spliced RNA and _t_ + 1: expression of new

unspliced RNA). Using the proportion between these two values for a given gene, and under some assumptions, it is possible to infer whether the expression of a gene is being initiated or

downregulated. By combining this information for all genes in each cell, as well as comparing it to its neighbor cells, it is possible to infer a direction vector indicating the putative

future transcriptional state of the cell. In our dataset, RNA velocity revealed a clear differentiation directionality from NCCs to neuronal progenies. In order to recapitulate such

transitions, the dataset was first projected onto diffusion space to denoise the underlying geometry, and a principal tree has been fitted using ElPiGraph11 in a semi-supervised way with the

help of clustering results (Fig. 1b, c). ElPiGraph is a manifold learning method which aims at inferring a principal graph (such as a tree) “passing through the middle of data” in high

dimension. Cells were ordered along the principal tree, and for each cell the distance on the graph to a chosen root is considered as a pseudotime value. Pseudotime analysis of gene

expression showed significant and robust changes along the reconstructed tree (8432 genes at a false discovery rate (FDR) of 0.05), leading to the identification of major clusters (Source

Data file). It also captured tree root populations and endpoints as well as intermediate states as sub-clusters, leading to the identification of two main neurogenic trajectories named

branch A and branch B (Fig. 1c). Those branches evolved from a clear transition from progenitors to post-mitotic newborn neurons. This was determined using Pagoda2 (Fig. 1d, e)12, which

separated the lineage tree into two major states, cycling (_Sox10__+_ and _Sox2__+_) versus non-cycling (_Isl1__+_), corresponding to non-neuronal versus neuronal populations, the latter

being defined by specific expression of axon-related cytoskeleton and function genes (including _Tubb3_, coding for βIII-tubulin) (Fig. 1d, e). Consistently, SOX10+ cells incorporated EdU

(cycling cells) and did not co-localize with the post-mitotic neuronal marker ISL1 (Supplementary Fig. 3a). Similarly, βIII-tubullin staining labeled sensory neurons and did not co-localize

with SOX2+ cells, which were also EdU+ (Supplementary Fig. 3b). PC scores for positive (GO:0030335) and negative (GO:0030336) regulation of cell migration GO term gene sets showed a

transition from migratory NCCs to neuronal progenitors (Fig. 1f). Interestingly, as NCCs migrate and coalesce into sensory ganglia in vivo, they showed enriched expression of genes

associated with interactions with the local microenvironment, including integrin subunits, the principal receptors for binding and interacting dynamically with the extracellular matrix,

which represents a necessary property of motile cells (Supplementary Fig. 3c). In contrast, neuronal progenitors and neurons were defined by a switch in expression of cadherin genes that

promote cell-to-cell adhesion (Supplementary Fig. 3d), which is important for cell positioning and border delimitation at a tissue level. Thus, these findings reveal a general principle and

molecular pathways governing the change in cell adhesiveness associated with migrating precursors and post migratory neuron progenitors. Within the cycling root populations, some cells

expressed the pro-neural basic helix loop helix transcription factors _Neurog1_ and _Neurog2_ (coding for neurogenins 1 and 2, NGN1 and 2), identifying the neurogenic niche prior to their

differentiation into post-mitotic sensory neurons. From these progenitors, branches A and B differentiated to endpoint neuronal clusters that express by E12.5 the molecular markers

characteristic of the main early DRG neuronal populations: the proprioceptive, mechanoreceptive and nociceptive populations. The proprioceptive lineage is characterized at E12.5 by the

co-expression of _Ntrk3_ (_Ntrk3_ coding for TRKC) with _Runx3_13,14,15 while the mechanoreceptive neurons express _Ret_ and _Ntrk2_ (_Ntrk2_ coding for TRKB)5,16 (Fig. 1g–i and Source Data

file). The mechanoreceptors further branches into _Ntrk2_+ only cells and _Ret_+/_Ntrk2_+ cells, which constitute subpopulations of this lineage with distinct connectivity and function in

the adult5,17,18. Altogether, those populations represent branch A. In contrast, branch B differentiated into nociceptive _Ntrk1_+ (TRKA) cells and is characterized by the expression of

_Ntrk1_ (and the beginning of induction of _Runx1 _in the TRKA+ cells) (Fig. 1g–i), a specific marker of the nociceptive lineage5. Combinatorial use of the above cited population markers

covers all (ISL1+)19 sensory neurons at this stage (Fig. 1i). Hence, by representing the main early sensory neuron categories, our dataset and analysis generated a refined order of the

somatosensory neurons diversification tree whereby progenitors undergo rapid transcriptional changes that progress through stepwise branching points and discrete sensory neuron states

consistent with a developmental ordering of the emergence of the different neuron types (Fig. 1j). This refined map helps in understanding the temporal sequence and lineage relationship

along the differentiation trajectories of the main somatic sensory neurons observed at this stage and will allow further explorations on the molecular principles underlying their generation.

STEM CELL POPULATIONS AND WAVES OF SENSORY NEUROGENESIS During neural crest migration, sequential pools of progenitors depend on either NGN2 or NGN1 to generate large-size neurons (mostly

TRKC+, TRKB+ and RET+ neurons) of the first wave, then small-size TRKA+ neurons of the second wave of neurogenesis, respectively20. A third wave of neurogenesis arises from boundary cap

cells (BCCs, which derive from NCCs too) and generates mostly small-size TRKA+ neurons21. One major question in the field is how cells progress from one state to another during neurogenesis

and how neurogenic programs specify the two major waves of neurogenesis leading to specific neuronal subtypes. To address those questions, we first identified the different clusters of cells

and trajectories. In our dataset, the clusters 1–3 (cycling cells), which define the origin of the differentiation tree, showed progressive downregulation of the remaining neural plate

border genes (e.g. _Zic_ and _Msx_, found in only few cells of cluster 1), upregulation of the NCCs specifiers (_Sox9_, _Foxd3_ and _Ets1_) and expression of genes involved in NCCs migration

(_Cdh11_ and _Itga4_, the latter being highly expressed in the bridge between clusters 2 and 3) (Fig. 2a, b). Also, following the RNA velocity stream, clusters 2 and 3 started to express

pro-neurogenic markers (_Neurog1_ and _Neurog2_) (Fig. 2c), defining those populations as neuronal progenitors. While transcriptionally close to cluster 2, the cluster 3 was marked by the

expression of genes specific of the BCCs (_Egr2_, _Ntn5_, _Prss56_, _Hey2_ and _Wif1_)22,23. We therefore labeled the three clusters as early NCCs (cluster 1, top differentially expressed

genes (TDEGs): _Prtg_, _Crabp1_, _Snai1_ and _Dlx2_), late NCCs (lNCCs, cluster 2, TDEGs: _Hoxd10_, _Hoxd11_, _Hoxa11_ and _Hoxc12_) and lNCCs/BCCs (cluster 3, TDEGs: _Fabp7_, _Sparc_,

_Serpine2_ and _Ednrb_) (Fig. 2a, b, Supplementary Fig. 4a and information on the BCCs in cluster 3 found in Source Data file). The alignment of the sequencing data revealed two trajectories

(or branches) emerging from the highly heterogeneous clusters NCCs and lNCCs/BCCs and named hereafter branch A and branch B leading to the onset of neurogenesis (Fig. 2c). Our analysis also

indicates that while neuronal progenitors from lNCCs/BCCs expressed only _Neurog1_, the two main neurogenesis trajectories (branches A and B, Fig. 2c) could not be distinguished based on

their specific NGN expression. Indeed, in both branches, progenitors evolved from a _Neurog2__+_ to a _Neurog1__+_ state (Fig. 2d) and numerous single cells were captured expressing both

transcripts simultaneously (Fig. 2e–g). These results contrast with previous data suggesting specific expression of NGN2 and of NGN1 within the first and second waves of neurogenesis,

respectively20. To confirm our results, we performed genetic cell-lineage tracing analysis of the NGN1 and NGN2 expressing progenitor cell lineages using _Neurog2__CreERT2_;_R26__tdTOM_ and

_Neurog1__Cre_;_R26__tdTOM_ mice. In both mice, recombination (TOM+) was observed in all classes of DRG neurons, regardless of their neurogenesis wave (Fig. 2h–l)24,25, confirming results

from a recent study26. Hence, we conclude that most sensory neuron progenitors sequentially express _Neurog2_ then _Neurog1_ with co-expression in an intermediate state where both genes are

co-expressed and suggest that the neuronal precursors generating the first or the second waves of neurogenesis are timely poised to preferentially require either NGN2 or NGN1 for neuronal

differentiation. NEUROGENESIS OF MYELINATED AND UNMYELINATED NEURONS We next investigated the genesis and key molecular determinants of myelinated versus unmyelinated neuronal

differentiation. To this end, we used diffusion maps to computationally compare the two neurogenesis trajectories, i.e. branches A and B (Fig. 2m). A common point representing the position

in diffusion space where the two branches were at their closest was inferred (Fig. 2n). Comparison was performed between the two branches starting from this inferred point. This analysis

revealed the modular pattern of gene expression along the trajectories. We could identify specific clusters of co-expressed genes uniquely distributed along or shared between each trajectory

and which constituted divergent and convergent genes, respectively (Fig. 2o, p and Source Data file). Comparing our results with the known regulatory interactions between pro-neural genes

from previous studies25, we identified a cascade of TFs expressed in each branch and which are likely to act in the acquisition of a neuronal fate. Hence, after _Sox10_ expression, _Neurog2_

was followed by _Neurog1_, then _Neurod4_, _Pou4f1_ (BRN3A), _Neurod1_ and _Isl1_. Convergent genes common to the two main trajectories defined a generic sensory differentiation program

(Fig. 2n). The common genes included _Neurod1_ and _Neurod4,_ which are conserved pro-neurogenic transcription factors and _Sox2_ and _Sox5_, which represent previously undescribed findings

and might play a specific role in the development of the somatosensory system (Fig. 2p and Supplementary Fig. 2). At the same time, divergent genes specific to either trajectory

distinguished the myelinated (branch A) from the unmyelinated lineage (branch B) (Fig. 2o). Specifically, we found the involvement of _Nfia_, _Gli2_ and _Neurod2_ in branch B (Fig. 2o and

Supplementary Fig. 2), suggesting critical role of distinct TFs during sensory neurogenesis of branches A and B. Overall, differentiation of myelinated versus unmyelinated neurons proceed

via different intermediate states demarcated by the expression of TFs, chromatin modifiers (_Prdm12_ as previously shown25 and _Prdm8_) and signaling molecules. Presumably, the dynamic local

signaling landscape might influence and bias the neuronal progenitors towards either branch A or B. Finally, separate GO term analysis on each of the gene sets (common and branch specific)

for cellular components showed concomitant enrichment for intracellular part (GO:0044424, GO:0005622). This, in line with the previous conclusion, implies that these waves of neurogenesis

are characterized by drastic reconfiguration of the internal components of the progenitor cells. FROM GENE MODULES TO SPECIFIC CELL FATE CHOICE Next, we questioned how the emergence of the

main neuronal subtypes is dictated at the transcriptomic level around the two last, most downstream, hierarchical fate split points in our dataset. For this, we examined the course of

transcriptional changes and the branching points representing neuronal subtype diversification events, focusing our analysis on the three downstream bifurcations along branch A (Fig. 3a).

Branch A gives rise to three mechanosensitive cell types, namely the RUNX3/TRKC proprioceptors and the two mechanoreceptors subpopulations TRKB/RET and TRKB only. Although the identity of

the sensory neuron clusters observed at E12.5 is consistent with previous knowledge in the field5, our aim was to explore the early genes possibly driving fate choice and fate biasing before

actual fate commitment. Each bifurcation could be broken down into three intermediate stages: pre-bifurcation (root), bifurcation point and post-bifurcation part of the trajectory (Fig. 

3b). The pre-bifurcation and post-bifurcation stages reveal the dynamic of transcriptional changes leading to fate choice and to its consolidation, respectively. We identified gene modules

(groups of genes that change in the same direction and tend to synchronize along the pseudotime) that we qualified as early and late and which characterized the pre- and post-bifurcation

stages of the first bifurcation respectively and therefore corresponded to fate choice and fate commitment. Within branch A, the first bifurcation corresponded to the binary cell fate choice

between mechanoreceptors and proprioceptors (Fig. 3b, c). We identified two main early gene modules associated with the fate split towards either proprioceptor (45 genes) or mechanoreceptor

(60 genes) fates at the pre-bifurcation stage (Fig. 3d) and late genes modules that define cell commitment (Fig. 3e, Supplementary Fig. 4 and Source Data file). Interestingly, before the

bifurcation node, the pattern of expression of those gene modules changed dynamically along the differentiation trajectory. Indeed, as cells progressed on the pseudotime axis, the two gene

modules were found originally co-expressed in single cells and then gradually mutually exclusive towards the bifurcation point correlating to preference in fate choice (Fig. 3d). GO term

analysis for cellular component showed enrichment in postsynaptic membrane (GO:0045211) and synapse (GO:0045202) for mechanoreceptors fate early gene module. Intrinsic components of membrane

(GO:0031224) and membrane parts (GO:0044425) are the top hits for the early gene module of the proprioceptors. Compared to transcriptomic hallmarks of the neurogenesis waves, this suggests

that remodeling in post-mitotic immature neurons likely affects the protein composition of the external parts of the cells, presumably for differential interactions with the environment.

Other genes are related to cytoskeleton remodeling, including _Cdh4_, _Pcdh9_, _Dock5_, _Dock3_ for the proprioceptors (mutation of _Dock3_ is known to results in ataxia) and _Stom_ and

_Pdlim1_ for mechanoreceptors, as well as ion channels (_Asic1_ and _Kcnip2_ for proprioceptive fate). We next analyzed the inference of cell fate decision governed by TF activity and found

TFs differentially expressed for both fates, where _Runx3_ and _Nfia_ are pro-proprioceptors and _Pou6f2, Nr5a2, Hoxb5, Egr1 and Junb_ are pro-mechanoreceptors (Fig. 3d and Source Data 

file). Interestingly, _Nfia_ and _Nr5a2_ were previously described in other systems to be linked to cell fate decision during neuronal development27 and _Runx3_, a main regulator of

proprioceptive neurons development13,28, was one of the two TFs that belong to the early gene module of the pro-proprioceptor trajectory. Hence, our data support a role for RUNX3 in cell

fate choice decision well before the bifurcation point between proprioceptor and mechanoreceptor neurons and suggest co-activation of competing programs prior to fate commitment. The second

bifurcation in the A-fibers low-thresold mechanoreceptors differentiation trajectory represented the fate choice between _Ntrk2_ only and _Ntrk2_/_Ret_ (Fig. 3f–h). In adults, those two

populations are further diversified and contribute to touch sensation and are distinguished by discrete innervation patterns of cutaneous end organs and different electrophysiological

properties29. Our data could identify early and late gene modules defining their fate split biasing fate towards either _Ntrk2_ only or _Ntrk2_/_Ret_ populations (Fig. 3i, j). GO term

enrichment analysis revealed enrichment in axon guidance receptor activity (GO:0008046) for _Ntrk2_ only fate. Similar to the previous bifurcation, _Ntrk2/Ret_ early gene module was enriched

in plasma-membrane-related genes (GO:0005886). Altogether, along the trajectory following the pseudotime, we identified bursting of modules of highly specific and coherent genes. Although

limited at this point to a transcriptomic analysis, these data suggest a competition between these early modules within the early cells of the myelinated lineage that would likely result in

biasing the future fate split towards either mechanoreceptors or proprioceptors (bifurcation 1, Fig. 4a) (or towards either two different subpopulations of tactile mechanoreceptors,

bifurcation 2, Fig. 4b) before the actual split occurs (Fig. 4a, b). Those early biasing modules appear initially positively correlated locally at the single-cell level before the

bifurcation (Fig. 4a, b, intra-module, see arrows), and present an increase of the negative correlation inter-modules at the pre-bifurcation stage (Fig. 4a, b, inter-module, see arrows, Fig.

 4c), hence suggesting co-activation of competing biasing programs prior to fate commitment (Fig. 4c–f). As a result, the retained program would participate then in the unfolding of the fate

towards one cell type versus the other. CELL FATE-BIASING PROGRAMS IN VIVO To investigate how TFs among early gene modules can affect fate choice in vivo, we focused on the role of RUNX3,

which was found expressed in part of the sensory neurons within the pre-bifurcation stage of the proprioceptor/mechanoreceptor trajectory and in all proprioceptors after bifurcation (Figs. 

3c, d, 5a, b). Using _Runx3__−/_−;_Bax__−/−_ animals, in which the deletion of _Bax_ allows the study of RUNX3 function in the absence of cell death28, we observed a threefold increase in

the proportion of TRKB+ cells in E12.5 DRG (Fig. 5c), similar to previous results13. This suggests a fate change toward mechanoreceptor fate in the absence of RUNX3. To investigate further a

potential cell fate change, we mapped the genes composing the early modules that we previously identified within the proprioceptor trajectory onto a published differential gene expression

analysis of (FACS isolated and RNA sequenced) _Runx3_ null and _Runx3__+/−_ GFP+ presumptive proprioceptors from E11.5 _Runx3_-_P2__GFP/GFP_ knock-in mice30. In those _Runx3_

promoter-specific knock-in mice, GFP expression pattern faithfully reproduces RUNX3 expression during development. This mouse model can thus be used to trace and specifically analyze neurons

of the RUNX3 lineage, expressing (in heterozygous state) or lacking (in homozygous state) _Runx3_ expression. In GFP+ cells lacking _Runx3_, a large number of early and late genes

identified in the pro-proprioceptor decision were significantly downregulated (i.e. _Runx3_, _Fam19a4_, _Cntnap2_, _Shisa6_, _Anxa2_, Fig. 5d, blue dots) while some identified as

pro-mechanoreceptor were significantly upregulated (i.e. _Ntrk2_, _Gfra2_, _Pou6f2_, _Gpm6a_, _Plxna2_, Fig. 5d, red dots), suggesting the existence of competing genetic programs prior to

cell fate choice during sensory neuron diversification. Interestingly, these differentiation programs not only activate genes specific of one cell fate but also repress gene programs of the

alternative fate. Despite the observation that the loss of RUNX3 in the early myelinated lineage results in a shift of the cell identity which underlines the weight of a single factor in

maintaining the integrity of a differentiation program, our computational analysis suggests a combinatorial coding of the cell identity, as previously suggested27,31. Together, those data

strongly suggest an essential role of gene expression dynamic prior to and at the bifurcation point in unfolding the transcriptional programs that defines fate choice towards proprioceptor

and mechanoreceptor fates and that competing programs between opposing fate determinants is an effective means to simultaneously promote one fate and exclude the other. CONVERGENCE OF

TRANSCRIPTIONAL PROGRAMS DURING NEURONAL DIFFERENTIATION Our last investigation was to define the specific signature of the most differentiated clusters, the endpoints of the trajectories

(Fig. 6a). Also, based on intersectional gene module of all endpoints (i.e. common genes), we ought to reveal the possible existence of a sensory neuron module that would define the sensory

lineage and which would be passed-on by the mother cell (the cycling progenitors) to all the daughters cells (postmitotic sensory neurons). Our data show that although the differentiation

trajectories exhibited divergent intermediate paths, they all seem to converge again at a transcriptional level as they mature to distinct neuron types, as judged by the relatively low

number of transcriptional heterogeneity between neuron types at the endpoints of our analysis compared to the progenitors population (Fig. 6b) and the high overlapping of GO term categories

(Fig. 6c). Our analysis could identify around 120−180 differentially expressed genes between clusters constituting the endpoints of the trajectories in branches A and B (Source Data file and

the final recapitulating scheme (Fig. 7)). Interestingly, some of these genes are markers specific of subpopulations (Fig. 6a) identified later during development. For example, within the

nociceptor lineage, _Trpv1_, _Th_ and _Trpm8_ were found to be expressed already at E12.5, yet they will each participate in the classification and function of discrete nociceptive neuron

subtypes in adults8, confirming the results of a recent study26. Altogether our findings provide a general principle for diversification of sensory neurons into subtypes which involves

repression of specific genes within intermediate cell states, leading to emergence of neuronal types with defined functional properties. The convergence of the population at the

transcriptomic level led us to examine whether common module genes were detectable at E12.5 and if they could be found at earlier developmental stages and preserved throughout the life. To

identify such a program that defines the entire sensory lineage, we back-traced modules of genes related to neuronal function expressed at E12.5 (Fig. 6d, e) and already present at the root

of neurogenesis. We identified 44 genes as being transmitted from mitotic neuronal progenitors to all post-mitotic sensory neurons at E12.5, and which code for generic neuronal features that

determine the neurotransmission phenotype and morphological features specific to neurons (most represented GO terms: axon neuron projection, myelin sheath and cell projection) (Fig. 6f and

Source Data file). Importantly, those genes were maintained at all levels of the diversification tree and were still found in adult sensory neurons using dataset from Mousebrain.org8

composed of three sub-taxon: peripheral sensory neurofilament neurons, peripheral sensory non-peptidergic neurons and peripheral sensory peptidergic neurons (Fig. 6g, h). Our data thus

provide insight into a set of genes defining the entire sensory lineage from sensory progenitors to adult sensory neurons (Fig. 6i). DISCUSSION The molecular mechanisms regulating the fate

of the somatosensory neurons are not entirely understood mainly because of the specific challenges including a relatively fast neurogenesis process and substantial mixing of progenitors and

cell types at the location where the sensory ganglia coalesce. As a result, the current paradigm established only a fraction of molecular determinants essential for the development of

particular sensory fates5. Moreover, these previously established molecular programs are likely associated with the processes of implementation of a fate rather than they are related to a

process of a fate choice, which might include complex accumulation of cell history, activation of biasing factors and the integrative role of dynamic signaling landscape. In this study, we

provide a comprehensive scRNAseq-based analysis of the neurogenesis in the mammalian somatosensory system and identified unique and mixed-lineage transcriptomic states that evolve to

culminate in specific neuronal subtypes (Fig. 7). Moreover, our study provides an analysis of the developmental dynamic throughout early stages revealing the emergence of early somatosensory

neurons. Our findings identify a combinatorial process through which cell type-specific neuronal identities emerge and reveal the preceding fate-related changes in heterogeneous population

of early dividing progenitors and neurons prior to fate branching points in the diversification tree leading to concrete sensory subtypes. Our branching model helps in understanding the

hierarchy of the sensory neuron lineage and reveals sequential binary decisions yielding to the main sensory subtypes. Such binary fate decision models32 are recognized for increasing the

robustness and reproducibility of cell type distribution (for review, see ref. 33). Binary decision among cycling cells leading to neuronal diversity has been largely described in the

context of asymmetric division in progenitors or lineage pre-patterning of precursors earlier during development1,33,34,35,36. Our data indicate that newly born sensory neurons within a

trajectory derive from a relatively homogeneous population of cells and that molecular specificity in neurons is progressively acquired along the pseudo-time axis of differentiation by a

progression of consecutive binary decision. The notion of  branche segregation  by mutual repressive activity had already been suggested for cell types diversification during late embryonic

somatosensory development, for instance during the differentiation of TRKA+/RUNX1+ nociceptors into a peptidergic and a non-peptidergic population of sensory neurons5,37, yet the timing and

gene-regulatory strategies that regulate this differentiation process remain poorly understood. Our data suggest a mechanism of distributing the fates of early neuron types in the developing

nervous system whereby fate-biasing heterogeneity of the post-mitotic neuronal population is gradually increasing until the commitment point in binary decision-making38. This is generally

similar and conserved during the formation of other fates in the multipotent NCC lineage tree32. The molecular underpinning of the bifurcation events along the hierarchical decision tree of

somatosensory development revealed initial co-activation of bi-potential transcriptional states in cells prior to bifurcation, followed by gradual shifts toward cell fate commitment and

binary decisions. It suggests that determination of a specific sensory neuron fate during early development is therefore achieved by the increased synchronization of relevant programs (as

key lineage priming factors) together with concurrent repression of competing fate programs. This complements the view of an early segregation of cell state and identification of cell fate

based on differential but unique expression of genes or gene modules that are necessary for their neuronal differentiation4,39,40. While restricted neuronal expression (in an on−off pattern)

based on progenitors’ fate-limiting TFs have been largely demonstrated in the control of a binary fate choice39,40, the sequence of transcriptional events and in general, the molecular

mechanisms operating in neurons as they diversify have been difficult to appreciate. However, complex neuronal diversification programs might be considered as the collection of multiple

binary fate decisions integrated over time. Therefore, the high-resolution analysis of the dynamic of transcriptional profiles in developing neurons presented here may thus provide an

analytical basis for studying further how neuron diversity is generated in other parts of the nervous system. Our analysis also complements a study recently published26, which presented the

transcriptional profile of the somatosensory neurons across various stages covering the embryonic and postnatal development. In contrast to our findings, Sharma et al. defined the early

postmitotic sensory neurons as unspecialized, whereas in our study they are already differentiated into specific sensory sub-branches. Our data confirm previous results demonstrating

neuronal heterogeneity in DRG as early as E11.5, with clearly defined populations of proprioceptor and tactile  mechanoreceptor at this stage13,14,15,16,17,28,41. Overall, our studies mostly

differ in their focus: Sharma et al. opted for analysis of large number of cells for offering a detailed developmental transcriptional atlas of sensory cell types while our study, focusing

on fewer early time points, provides a mechanistic insight of the molecular events leading to the early emergence of sensory neuron diversity. The mosaic differentiation pattern with a

relatively constant proportion of generated sensory neurons and cell types within the somatosensory system makes it unlikely to be constrained by purely deterministic mechanisms. Instead, it

raises the possibility of a stochastic fate decision within immature neurons in which one of their bivalent molecular programs would get stabilized and would thus progressively drive one of

the two possible fates. Such mechanism has been demonstrated during the differentiation of photoreceptor neurons in _Drosophila_, where the stochastic induction of a TF, Spineless, controls

cell fate decision cell-autonomously and simultaneously coordinates the alternate fate identity in neighboring photoreceptors42. In our system, we show that during decision between

proprioceptor and mechanoreceptor lineages, high levels of RUNX3 expression are required to specify a proprioceptive fate. Hence, in _Runx3_ mutants, cells of the RUNX3 lineage express

molecular determinants characteristic of the early and late gene modules of the alternate tactile mechanoreceptive cell lineage. The mechanism by which _Runx3_ becomes heterogeneously

expressed in a subset of immature neurons remains however unknown and is most likely driven by differences in microenvironment including possible crosstalk with other differentiating

subtypes or access to diffusing molecules. Upstream transcriptional regulators are likely functioning in coordinating the expression state of these key neuronal TFs. It will thus be

interesting to see whether and how upstream mechanisms tightly control the activity of the various promoters and regulatory elements that promote _Runx3_ expression specifically in the

presumptive proprioceptor lineage30. Importantly, by integrating multiple steps of cell identity regulation, such mechanism operating along the hierarchical differentiation tree could help

ensuring a correct representation of all neuronal subtypes developing from a common progenitor pool, especially via influencing the fate choice dynamics. In conclusion, our work has provided

a detailed transcriptional roadmap of neurogenesis and sensory neurons development that captured the dynamic of molecular events participating in cell fate choice, with progression through

the bifurcation and lineage specification of the sensory neurons. We identified dynamic burst of modules of genes, which led to branching point within lineages (Fig. 7). The identification

of genes modules that prime sensory neurons to specific fate might be of benefit for engineering-induced pluripotent stem cell and embryonic stem cell technologies, and help advancing our

ability to engineer specific neuronal populations for basic research, tissue regeneration, or for screening drugs against neurodegenerative or pain disorders. Also, this might help

deciphering the combinatorial code of TFs responsible for cell state and to manipulate it in a context of a loss of cellular homeostasis with the hope to reset the cells to a healthy state.

METHODS ANIMALS Wild-type C57BL6 mice were used unless specified otherwise. _Plp1__CreERT2_, _Wnt1__cre_, _Isl1__Cre_, _Ntrk3_Cre, _Runx3__−/−_, _Fam19a4__YFP_, _Neurog2__CreERT2_,

_Neurog1__Cre_, _Bax__−/−_ and _Ai14_ (_Rosa26__tdTOM_) mouse strains have been described elsewhere14,25,28,43,44. Animals of either sex were included in this study. Animals were

group-housed, with food and water ad libitum, under 12 h light–dark cycle conditions. All animal work was performed in accordance with the national guidelines and approved by the local

ethics committee of Stockholm, Stockholms Norra djurförsöksetiska nämnd. TREATMENT For cell fate tracing experiments with inducible mouse lines, tamoxifen (Sigma, T5648) was dissolved in

corn oil (Sigma, 8267), and delivered via intraperitoneal (i.p.) injection to pregnant females at E9.5 (100 µg/g of bodyweight). For cell cycle study, intraperitoneal injection of pregnant

females with EdU (100 mg/kg, Invitrogen) was performed at E10 and E12 days of gestation. Injected females were killed 2 h after injection for analysis. EdU incorporation was subsequently

resolved using Alexa Fluor 488 azide according to the manufacturer’s instructions (Invitrogen). IMMUNOSTAININGS AND RNASCOPE® IN SITU HYBRIDIZATION Animals were collected, decapitated and

fixed for 1–4 h at +4 °C with 4% paraformaldehyde in phosphate-buffered saline (PBS) depending on the stage, washed in PBS, cryopreserved in 30% sucrose in PBS, embedded in optimal cutting

temperature (Tissue-Tek) and cryosectioned at 14 μm. In situ hybridization was performed using standard RNAscope protocol (ACDBio). The RNAscope probes used in this study are Mm-Fxyd7,

Mm-Neurog1, Mm-Neurog2, Mm-Spp1, Mm-Ahr, Mm-Girk1, Mm-Cbln4, Mm-Chodl and Mm-Anxa2 (ACDBio). Sections were incubated for 24 h at +4 °C with primary antibodies diluted in blocking solution

(2% donkey serum, 0.0125% NaN3, 0.5% Triton X-100 in PBS). Primary antibodies used were: rabbit anti-RUNX3 (gift from T. M. Jessell), rabbit anti-RUNX144, goat anti-TRKA (1:400, R&D

Systems AF1056), goat anti-TRKB (1:500, R&D Systems AF1494), rabbit anti-TRKC (1:1000, Cell Signaling 3376), mouse anti-ISL1 (1:250, Developmental Studies Hybridoma Bank 39.4D5), goat

anti-RET (1:100, R&D Systems AF482), goat anti-TRKC (1:500, R&D Systems AF1404), mouse anti-βIII-tubulin (1:1000, Promega G712A), chicken anti-RFP (1:250, Rockland 600-901-379S),

mouse anti-NF200 (Sigma, N0142; 1:500), mouse anti-SOX10 (1:1000, Santa-Cruz, #sc-36569). After washing with PBS, Alexa Fluor secondary antibodies (Live Technology; 1:500 in blocking

solution) were applied overnight (at +4 °C). Samples were then washed in PBS and mounted in DAKO fluorescent mounting medium. Staining was documented by confocal microscopy (Zeiss LSM700)

using identical settings between control and experimental images. Optical sections were 2 μm in ×20 overview pictures unless specified. QUANTIFICATIONS For cell type quantifications, ImageJ

software was used. Only neurons with a visible nucleus were considered for analysis. Quantification of molecular markers in the DRG was carried out on five DRG sections/animal, selected from

the most equatorial region of each DRG and covering the segments C5-T1 (see figure legends for _n_ values and genotypes). STATISTICAL ANALYSIS Data were analyzed using GraphPad Prism 6 and

expressed as mean ± s.e.m. The statistical test performed is reported in the figure legend. _t_ tests were two-sided. Legend for significance: *_P_ ≤ 0.05, **_P_ ≤ 0.01, ***_P_ ≤ 0.001. No

animals or data points were excluded from the analysis. Our sample sizes are similar to those generally employed in the field. SINGLE-CELL ISOLATION FOR SINGLE-CELL ANALYSIS Brachial DRGs

were dissected (from E11 embryos) or whole embryonic trunks collected (E10.5 and younger) in Leibovit’z L-15 medium (Life Technologies) on ice. The tissue was incubated in 0.05% trypsin-EDTA

(1×) (Life Technologies) for 5, 7, 10 and 15 min for tissue obtained from E9.5, E10.5, E11.5 and E12.5, respectively, at 37 °C, in a thermomixer comfort (Eppendorf) at 700 RPM. After

spinning down the samples at 100 RCF for 5 min, the supernatant was removed and replaced by Leibovit’z L-15 medium (Life Technologies). The tissue following enzymatic digestion was

physically triturated using two different sizes of pipettes previously coated with 0.2% bovine serum albumin until the solution homogenized. The cell suspension was then filtered through a

70-μm cell strainer (BD Biosciences) to remove the clusters of cells. FACS For the scope of this work, we sorted TOMATO+ cells from different transgenic lines (both inducible and

non-inducible) using either instruments from BD (manufacturer) (BD FACSDiva 8.0.2 for _Wnt1__cre_;_Ai14_, _Plp1__creERT2_;_Ai14_ and _Ai14_ and BD Influx for _Isl1__cre_;_Ai14_ and

_Ntrk3__cre_;_Ai14_. The software used is FlowJo v10 (BD) and BD FACS Software 1.2.0/142/Utopex 1.2.0.108 respectively. The gating strategy is shown in Supplementary Fig. 2. Debris and

erythrocytes were gated out from the total events using FSC-A versus SSC-A plotting. Doublets were further gated out using FSC-A versus FSC-H plotting. Finally, we plotted FSC-A versus PE-A

(corresponding to TOMATO+ signal) to select cells that were TOMATO+ (Supplementary Fig. 2a, b, e) or we plotted the FSC versus the 585/29 [561]-tdtomato (corresponding to TOMATO+ signal,

Supplmentary Fig. 2c, d). A cell preparation negative for TOMATO was used to define the negative population (Supplementary Fig. 2e). Note, the BD Influx instrument was manually set up and

calibrated on daily basis using BD CST and BD FACS Accudrop beads. Single TOMATO+ cells were sorted by fluorescence-activated cell sorting (FACS) into individual wells containing lysis

buffer in 384-well plates provided by the facility (Karolinska Insitutet). All plates were immediately placed on dry ice and stored at −80 °C before processed for Smart-seq2 protocol at the

single-cell facility (Karolinska Institutet). SINGLE-CELL SEQUENCING AND GENERATION OF COUNT MATRICES, QC AND FILTERING The single-cell transcriptome data were generated, using Smart-seq2

protocol at the Eukaryotic Single-cell Genomics facility at Science for Life Laboratory in Stockholm, Sweden. This sequencing approach does not use UMI and hence does not correct for PCR

duplicates. The samples were analyzed by first demultiplexing the fastq files using deindexer (https://github.com/ws6/deindexer) using the nextera index adapters and the 384-well layout.

Individual fastq files were then mapped to mm10_ERCC genome using the STAR aligner using 2-pass alignment45. Reads were filtered for only uniquely mapped and were saved in BAM file format,

count matrices were subsequently produced. Estimated count matrices were gathered and converted to a SingleCellExperiment object, QC metrics were computed using calculateQCMetrics function.

Cells having less that 5 × 104 transcripts, less than 2500 genes and more than 25% of proportion of ERCC reads were filtered out. For the selection for the developing sensory cells, the main

analysis of the count matrix was mainly performed using pagoda2 R package12. The filtered overview count matrix gene expression variance was adjusted (pagoda2, _k_ = 10) and 5801

overdispersed genes were detected. PCA was performed on the overdispersed genes (pagoda2, nPcs = 100, maxit = 1000) and 13 PC were retained using elbow curve selection. KNN graph (pagoda2,

_k_ = 40, centered, cosine distance) and UMAP visualization (umap python, _n__neighbors = 30, min_dist = 0.5) were generated in PCA space. Clusters were then identified on the KNN graph

using leiden algorithm46, and pathway overdispersion analysis was performed (pagoda2, correlation.distance.threshold = 0.95) to identify relevant biological aspects. Differential gene

expression was performed for each of the detected clusters using Wilcoxon rank test. Cells in clusters with _Sox10_ gene positively differentially expressed were kept, as well as linked

neurogenic clusters (using aspect1). To cluster and visualize the developing sensory data, the subsetted count matrix was further processed using pagoda2 with similar parameters as for the

overview (4887 detected overdispersed genes), clusters were detected using infomap algorithm in addition to leiden for a finer clustering. A biological aspect linked to mitochondrial

respiratory chain has been identified to be different among batches from the same condition. This aspect was regressed on the raw count matrix using ScaleData from Seurat package47, the

regressed count matrix was further processed through the Seurat 2 pipeline, variable genes identification (default parameters) and PCA (pcs.compute = 100), eight PCs were retained using

elbow curve selection. UMAP was performed on the retained PCs (_n__neighbors = 100, min_dist = 0.5). For the heterogeneity analysis, a bootstrapping analysis was performed for each timepoint

separately to validate that cell type heterogeneity decreases upon sensory neuron differentiation. One hundred cells were randomly one after the other sampled with replacement, top 100 most

expressed genes per selected cell were identified and the number of newly identified genes were added to the number of unique genes seen before. This process was repeated 100 times for each

timepoint. Count matrix correction was performed prior to pseudotime analysis, the raw count matrix was corrected using scde R package. scde fits error models for each cell in order to

estimate the drop-out and amplification biases on gene expression magnitude48. In order to infer trajectories on the resulting data, diffusion maps were first computed using Palantir49

python package (run_diffusion_maps, default parameters) on PCA space. Principal graphs were fitted on the first five diffusion components using ElPiGraph50. To avoid the E12.5 clusters

describing the proprioceptors and the mechanoreceptors to be wrongly associated with the E12.5 nociceptive lineage, the principal tree was performed in several steps. First, a principal tree

of 30 nodes was generated on the E10.5/11.5 subset of the data, by excluding leiden clusters 5, 6, 7, 9, 10, and infomap cluster 23. Proprioceptor cluster (leiden 10) was added by linking

its geometrical center in diffusion space to the tree. Mechanoreceptors cluster (leiden 9) was added by inferring a 4-node principal tree in diffusion space on this subset, and by linking it

to the E10.5/E11.5 tree using closest node. Trajectory for nociceptive lineage was recovered by inferring a 15-node principal curve in diffusion space on leiden clusters 5 and 6, the

resulting principal curve was then linked to the whole tree by using closest nodes. In order to have branches that have at least more than one node, additional nodes were added by separating

each edge of the tree in two equidistant ones. Downstream analysis was performed using a slightly modified crestree R package (https://github.com/LouisFaure/crestree) provided by the

previous neural crest study32, in which precise description of the underlying statistics is mentioned. The difference in its usage is first the inference of a principal tree with ElPiGraph,

which contains two roots and is mapped only once, and second the usage of the updated JASPAR2018 database50 for TF activity inference. Associated genes with the tree were detected using

test.associated.genes (default parameters) on the log transformed gene expression. Associated genes were fitted using fit.associated.genes (gamma = 5). The fitted profiles were clustered

into 30 major patterns with hierarchical clustering using Ward method and Euclidean distance. Inference of TF activity was performed on 101 selected TF using activity.lasso (default

parameters). Bifurcation analysis was performed in three steps, genes differentially upregulated after bifurcation point were detected using test.fork.genes (default parameters), and

differentially expressed genes were then assigned between two post-bifurcation branches using branch.specific.genes. Optimum expression and time of activation were estimated for each of the

detected gene using activation.fork (default parameters), allowing the separation between early and late modules by setting a pseudotime threshold before the bifurcation. To analyze

molecular mechanisms of cell fate selection, cell composition was approximated by a sliding window of cells along the pseudotime axis, cells were manually selected in order to represent the

different steps of differentiation. The local gene−gene correlation reflecting the coordination of genes around a given pseudotime _t_ was defined as a gene−gene Pearson correlation within

each window of cells. The local correlation of a gene _g_ with a module was assessed as a mean local correlation of that gene with the other genes comprising the module. Similarly,

intra-module and inter-module correlations were taken to be the mean local gene−gene correlations of all possible gene pairs inside one module, or between the two modules, respectively51.

For RNA velocity analysis, BAM files from each plate were processed using python command-line velocity tool10, using run-smartseq2 command with GENCODE M21 genome and repeat masker

annotation files, leading to a loom file for each plate containing spliced and unspliced transcript counts, which are then combined in one loom dataset, and cells are filtered out according

to the cells kept in the final developing sensory dataset. Using scvelo python tool52, genes having <20 spliced counts or genes having <10 unspliced counts were excluded and the 4000

top highly variable genes were kept. In addition, cell cycle genes are also filtered out from the analysis. PCA was performed on the spliced matrix, keeping 30 principal components and kNN

neighbor graph was produced with _k_ = 30. Moments of spliced/unspliced abundances, velocity vectors and velocity graph were computed using default parameters. Extrapolated states are then

projected on the UMAP embedding produced during the initial analysis step. In order to compare the two neurogenic waves, the early E10.5 and late E11.5/12.5 neurogenic waves, associated

genes with each waves were identified using test.associated.genes (default parameters) on their respective subtrees, having as a starting point the first crossing point of the tree, and the

endpoint the next bifurcation for the early wave, or the endpoint of the late wave. Common genes were identified by the intersection of the two set of detected genes, the other genes were

considered as wave specific. For common and differing genes respectively, early and late genes were identified separating the fitted profiles in two clusters via hierarchical clustering

using Ward method and Euclidean distance. Endpoints trajectories identity were defined by performing differential gene expression using Wilcoxon rank test on infomap clusters 8, 10, 15, 20

and 22, genes with log2fc higher than 1 and expressed in more than 90% of the cells of a given cluster were considered as a marker for identity. GO term enrichment analysis was performed

with topGo R package on these markers for ontologies “biological process” as well as “cellular compartment”. Test was performed using “elim” algorithm and “fisher” statistical testing. For

the comparison with Runx3 KO bulk dataset, differential gene expression results of bulk data from Apple et al.30 were obtained from GEO database (GSE81140), early genes for each bifurcation

were overlaid on the volcano plot to identify their distribution of the genes between wild type and Runx3 mutants. To identify the common early neuronal identity, the tree was separated in

two subtrajectories, with one having root in the E10.5 NCCs population, and endpoints E10.5/E11.5 and E12.5 mechanoreceptors and proprioceptors, and the other having the E12.5 NCCs

population as root and the E12.5 nociceptive neurons as endpoint. A gene defining neuronal identity is defined as a gene being activated in the NCCs pool (pseudotime < 6) and present in

more than 98% of the cells which at pseudotime > 6 in all each single linear trajectory (from root to end), activation.statistics was used to identify early genes being activated in the

NCCs pool. Expression of each of the detected gene markers was checked in adult mice, by looking at their levels in aggregated cluster data of peripheral sensory neurons from mouse brain

atlas8 (http://mousebrain.org/loomfiles_level_L6.html). REPORTING SUMMARY Further information on research design is available in the Nature Research Reporting Summary linked to this article.

DATA AVAILABILITY The authors declare that all data supporting the findings of this study are available within the article and its supplementary information files or from the corresponding

author upon reasonable request. Raw sequencing data have been deposited in the GEO database under accession code: GSE150150. This includes a pagoda2 web file (p2w_sensory.bin) allowing

exploration dataset on a web browser. The file can be opened on the following link: http://pklab.med.harvard.edu/nikolas/pagoda2/frontend/current/pagodaLocal/index.html. For the Comparison

with Runx3 KO bulk dataset, differential gene expression results of bulk data30 were obtained from the GEO database (GSE81140). Source data are provided with this paper. CODE AVAILABILITY

All codes for scRNAseq data preprocessing and pseudotime analysis are deposited under the form of notebooks on the following github repository:

https://github.com/LouisFaure/sensoryfates_paper. Source data are provided with this paper. REFERENCES * Telley, L. et al. Temporal patterning of apical progenitors and their daughter

neurons in the developing neocortex. _Science_ 364, eaav2522 (2019). Article  CAS  PubMed  Google Scholar  * Mi, D. et al. Early emergence of cortical interneuron diversity in the mouse

embryo. _Science_ 360, 81–85 (2018). Article  ADS  CAS  PubMed  PubMed Central  Google Scholar  * Wamsley, B. & Fishell, G. Genetic and activity-dependent mechanisms underlying

interneuron diversity. _Nat. Rev. Neurosci._ 18, 299–309 (2017). Article  CAS  PubMed  Google Scholar  * Greig, L. C., Woodworth, M. B., Galazo, M. J., Padmanabhan, H. & Macklis, J. D.

Molecular logic of neocortical projection neuron specification, development and diversity. _Nat. Rev. Neurosci._ 14, 755–769 (2013). Article  CAS  PubMed  Google Scholar  * Lallemend, F.

& Ernfors, P. Molecular interactions underlying the specification of sensory neurons. _Trends Neurosci._ 35, 373–381 (2012). Article  CAS  PubMed  Google Scholar  * Marmigere, F. &

Ernfors, P. Specification and connectivity of neuronal subtypes in the sensory lineage. _Nat. Rev. Neurosci._ 8, 114–127 (2007). Article  CAS  PubMed  Google Scholar  * Lawson, S. N. &

Biscoe, T. J. Development of mouse dorsal root ganglia: an autoradiographic and quantitative study. _J. Neurocytol._ 8, 265–274 (1979). Article  CAS  PubMed  Google Scholar  * Zeisel, A. et

al. Molecular architecture of the mouse nervous system. _Cell_ 174, 999–1014.e1022 (2018). Article  CAS  PubMed  PubMed Central  Google Scholar  * Usoskin, D. et al. Unbiased classification

of sensory neuron types by large-scale single-cell RNA sequencing. _Nat. Neurosci._ 18, 145–153 (2015). Article  CAS  PubMed  Google Scholar  * La Manno, G. et al. RNA velocity of single

cells. _Nature_ 560, 494–498 (2018). Article  ADS  PubMed  PubMed Central  CAS  Google Scholar  * Albergante, L. et al. Robust and scalable learning of data manifolds with complex topologies

via ElPiGraph. _Entropy_ 22, 1–27 (2020). Article  MathSciNet  Google Scholar  * Fan, J. et al. Characterizing transcriptional heterogeneity through pathway and gene set overdispersion

analysis. _Nat. Methods_ 13, 241–244 (2016). Article  CAS  PubMed  PubMed Central  Google Scholar  * Kramer, I. et al. A role for Runx transcription factor signaling in dorsal root ganglion

sensory neuron diversification. _Neuron_ 49, 379–393 (2006). Article  CAS  PubMed  Google Scholar  * Wang, Y. et al. A cell fitness selection model for neuronal survival during development.

_Nat. Commun._ 10, 4137 (2019). Article  ADS  CAS  PubMed  PubMed Central  Google Scholar  * Levanon, D. et al. The Runx3 transcription factor regulates development and survival of TrkC

dorsal root ganglia neurons. _EMBO J._ 21, 3454–3463 (2002). Article  CAS  PubMed  PubMed Central  Google Scholar  * Bourane, S. et al. Low-threshold mechanoreceptor subtypes selectively

express MafA and are specified by Ret signaling. _Neuron_ 64, 857–870 (2009). Article  CAS  PubMed  Google Scholar  * Luo, W., Enomoto, H., Rice, F. L., Milbrandt, J. & Ginty, D. D.

Molecular identification of rapidly adapting mechanoreceptors and their developmental dependence on ret signaling. _Neuron_ 64, 841–856 (2009). Article  CAS  PubMed  PubMed Central  Google

Scholar  * Li, L. et al. The functional organization of cutaneous low-threshold mechanosensory neurons. _Cell_ 147, 1615–1627 (2011). Article  CAS  PubMed  PubMed Central  Google Scholar  *

Sun, Y. et al. A central role for Islet1 in sensory neuron development linking sensory and spinal gene regulatory programs. _Nat. Neurosci._ 11, 1283–1293 (2008). Article  CAS  PubMed 

PubMed Central  Google Scholar  * Ma, Q., Fode, C., Guillemot, F. & Anderson, D. J. Neurogenin1 and neurogenin2 control two distinct waves of neurogenesis in developing dorsal root

ganglia. _Genes Dev._ 13, 1717–1728 (1999). Article  CAS  PubMed  PubMed Central  Google Scholar  * Maro, G. S. et al. Neural crest boundary cap cells constitute a source of neuronal and

glial cells of the PNS. _Nat. Neurosci._ 7, 930–938 (2004). Article  CAS  PubMed  Google Scholar  * Coulpier, F. et al. Novel features of boundary cap cells revealed by the analysis of newly

identified molecular markers. _Glia_ 57, 1450–1457 (2009). Article  PubMed  Google Scholar  * Garrett, A. M. et al. Analysis of expression pattern and genetic deletion of Netrin5 in the

developing mouse. _Front. Mol. Neurosci._ 9, 3 (2016). Article  PubMed  PubMed Central  CAS  Google Scholar  * Zirlinger, M., Lo, L., McMahon, J., McMahon, A. P. & Anderson, D. J.

Transient expression of the bHLH factor neurogenin-2 marks a subpopulation of neural crest cells biased for a sensory but not a neuronal fate. _Proc. Natl Acad. Sci. USA_ 99, 8084–8089

(2002). Article  ADS  CAS  PubMed  PubMed Central  Google Scholar  * Bartesaghi, L. et al. PRDM12 is required for initiation of the nociceptive neuron lineage during neurogenesis. _Cell

Rep._ 26, 3484–3492.e3484 (2019). Article  CAS  PubMed  PubMed Central  Google Scholar  * Sharma, N. et al. The emergence of transcriptional identity in somatosensory neurons. _Nature_ 577,

392–398 (2020). Article  ADS  CAS  PubMed  PubMed Central  Google Scholar  * Bikoff, J. B. et al. Spinal inhibitory interneuron diversity delineates variant motor microcircuits. _Cell_ 165,

207–219 (2016). Article  CAS  PubMed  PubMed Central  Google Scholar  * Lallemend, F. et al. Positional differences of axon growth rates between sensory neurons encoded by runx3. _EMBO J._

31, 3718–3729 (2012). Article  CAS  PubMed  PubMed Central  Google Scholar  * Abraira, V. E. & Ginty, D. D. The sensory neurons of touch. _Neuron_ 79, 618–639 (2013). Article  CAS 

PubMed  Google Scholar  * Appel, E. et al. An ensemble of regulatory elements controls Runx3 spatiotemporal expression in subsets of dorsal root ganglia proprioceptive neurons. _Genes Dev._

30, 2607–2622 (2016). Article  CAS  PubMed  PubMed Central  Google Scholar  * Hobert, O., Carrera, I. & Stefanakis, N. The molecular and gene regulatory signature of a neuron. _Trends

Neurosci._ 33, 435–445 (2010). Article  CAS  PubMed  PubMed Central  Google Scholar  * Soldatov, R. et al. Spatiotemporal structure of cell fate decisions in murine neural crest. _Science_

364, eaas9536 (2019). Article  CAS  PubMed  Google Scholar  * Jukam, D. & Desplan, C. Binary fate decisions in differentiating neurons. _Curr. Opin. Neurobiol._ 20, 6–13 (2010). Article

  CAS  PubMed  PubMed Central  Google Scholar  * Bertrand, V. & Hobert, O. Linking asymmetric cell division to the terminal differentiation program of postmitotic neurons in C. elegans.

_Dev. Cell_ 16, 563–575 (2009). Article  CAS  PubMed  PubMed Central  Google Scholar  * Poole, R. J. & Hobert, O. Early embryonic programming of neuronal left/right asymmetry in C.

elegans. _Curr. Biol._ 16, 2279–2292 (2006). Article  CAS  PubMed  Google Scholar  * Li, Q. et al. A functionally conserved gene regulatory network module governing olfactory neuron

diversity. _PLoS Genet._ 12, e1005780 (2016). Article  PubMed  PubMed Central  CAS  Google Scholar  * Gascon, E. et al. Hepatocyte growth factor-Met signaling is required for Runx1

extinction and peptidergic differentiation in primary nociceptive neurons. _J. Neurosci._ 30, 12414–12423 (2010). Article  CAS  PubMed  PubMed Central  Google Scholar  * Taverna, E., Gotz,

M. & Huttner, W. B. The cell biology of neurogenesis: toward an understanding of the development and evolution of the neocortex. _Annu. Rev. Cell Dev. Biol._ 30, 465–502 (2014). Article

  CAS  PubMed  Google Scholar  * Peng, Y. R. et al. Binary fate choice between closely related interneuronal types is determined by a Fezf1-dependent postmitotic transcriptional switch.

_Neuron_ 105, 464–474.e6 (2019). Article  PubMed  CAS  PubMed Central  Google Scholar  * Delile, J. et al. Single cell transcriptomics reveals spatial and temporal dynamics of gene

expression in the developing mouse spinal cord. _Development_ 146, dev173807 (2019). Article  CAS  PubMed  PubMed Central  Google Scholar  * Wende, H. et al. The transcription factor c-Maf

controls touch receptor development and function. _Science_ 335, 1373–1376 (2012). Article  ADS  CAS  PubMed  Google Scholar  * Wernet, M. F. et al. Stochastic spineless expression creates

the retinal mosaic for colour vision. _Nature_ 440, 174–180 (2006). Article  ADS  CAS  PubMed  Google Scholar  * Adameyko, I. et al. Schwann cell precursors from nerve innervation are a

cellular origin of melanocytes in skin. _Cell_ 139, 366–379 (2009). Article  CAS  PubMed  Google Scholar  * Hadjab, S. et al. A local source of FGF initiates development of the unmyelinated

lineage of sensory neurons. _J. Neurosci._ 33, 17656–17666 (2013). Article  CAS  PubMed  PubMed Central  Google Scholar  * Dobin, A. et al. STAR: ultrafast universal RNA-seq aligner.

_Bioinformatics_ 29, 15–21 (2013). CAS  PubMed  Google Scholar  * Traag, V. A., Waltman, L. & van Eck, N. J. From Louvain to Leiden: guaranteeing well-connected communities. _Sci. Rep._

9, 5233 (2019). Article  ADS  CAS  PubMed  PubMed Central  Google Scholar  * Butler, A., Hoffman, P., Smibert, P., Papalexi, E. & Satija, R. Integrating single-cell transcriptomic data

across different conditions, technologies, and species. _Nat. Biotech._ 36, 411–420 (2018). Article  CAS  Google Scholar  * Kharchenko, P. V., Silberstein, L. & Scadden, D. T. Bayesian

approach to single-cell differential expression analysis. _Nat. Methods_ 11, 740–742 (2014). Article  CAS  PubMed  PubMed Central  Google Scholar  * Setty, M. et al. Characterization of cell

fate probabilities in single-cell data with Palantir. _Nat. Biotechnol._ 37, 451–460 (2019). Article  CAS  PubMed  PubMed Central  Google Scholar  * Chen, H. et al. Single-cell trajectories

reconstruction, exploration and mapping of omics data with STREAM. _Nat. Commun._ 10, 1903 (2019). Article  ADS  PubMed  PubMed Central  CAS  Google Scholar  * Khan, A. et al. JASPAR 2018:

update of the open-access database of transcription factor binding profiles and its web framework. _Nucleic Acids Res._ 46, D260–D266 (2018). Article  CAS  PubMed  Google Scholar  * Bergen,

V., Lange, M., Peidli, S., Wolf, F. A. & Theis, F. J. Generalizing RNA velocity to transient cell states through dynamical modeling. _Nat. Biotechnol._

https://doi.org/10.1038/s41587-020-0591-3 (2020). Download references ACKNOWLEDGEMENTS We thank Y. Groner and D. Levanon for the _Runx3_ mice; C. Ibanez for the _R26__tdTOM_ mice; A. Moqrich

for the _Fam19a4__YFP_; T. Jessell for the RUNX3 antibody; the CLICK imaging Facility supported by the Knut and Alice Wallenberg Foundation; Single Cell Facility at SciLife Laboratory. We

are thankful to Jaromir Mikes for sorting the cells. We are grateful to Prof. Patrik Ernfors and Prof. Abdel El Manira for critical reading of the manuscript. We are also grateful to Dr.

Ruslan A. Soldatov (Kharchenko Lab) for sharing codes. I.A. was supported by ERC Consolidator grant “STEMMING-FROM-NERVE”, Bertil Hallsten Foundation, Swedish Research Council,

Paradifference Foundation. F.L. was supported by the Swedish Research Council (VR), the Ragnar Söderberg Foundation, Knut and Alice Wallenberg Foundation, Swedish Brain Foundation,

Karolinska Institutet, the Karolinska Institutet Strategic Research program in Neuroscience (StratNeuro) and the Ming Wai Lau Foundation. L.F. is supported by the Austrian Science Fund DOC

33-B27. S.H. is supported by the Swedish Research Council (VR), the Swedish Brain Foundation and Karolinska Institutet. Open access funding provided by Karolinska Institute. AUTHOR

INFORMATION AUTHORS AND AFFILIATIONS * Department of Molecular Neurosciences, Center for Brain Research, Medical University Vienna, 1090, Vienna, Austria Louis Faure, Maria Eleni Kastriti 

& Igor Adameyko * Department of Neuroscience, Karolinska Institutet, Stockholm, Sweden Yiqiao Wang, Paula Fontanet, Kylie K. Y. Cheung, Charles Petitpré, Haohao Wu, Lynn Linyu Sun, 

François Lallemend & Saida Hadjab * INMED INSERM U1249, Aix-Marseille University, Marseille, France Karen Runge & Antoine de Chevigny * Università Vita-Salute San Raffaele, 20132,

Milan, Italy Laura Croci & Gian Giacomo Consalez * Department of Neuroscience, UT Southwestern Medical Center, 5323 Harry Hines Boulevard, Dallas, TX, 75390, USA Mark A. Landy & 

Helen C. Lai * Ming-Wai Lau Centre for Reparative Medicine, Stockholm node, Karolinska Institutet, Stockholm, Sweden François Lallemend * Department of Physiology and Pharmacology,

Karolinska Institutet, 17177, Stockholm, Sweden Igor Adameyko Authors * Louis Faure View author publications You can also search for this author inPubMed Google Scholar * Yiqiao Wang View

author publications You can also search for this author inPubMed Google Scholar * Maria Eleni Kastriti View author publications You can also search for this author inPubMed Google Scholar *

Paula Fontanet View author publications You can also search for this author inPubMed Google Scholar * Kylie K. Y. Cheung View author publications You can also search for this author inPubMed

 Google Scholar * Charles Petitpré View author publications You can also search for this author inPubMed Google Scholar * Haohao Wu View author publications You can also search for this

author inPubMed Google Scholar * Lynn Linyu Sun View author publications You can also search for this author inPubMed Google Scholar * Karen Runge View author publications You can also

search for this author inPubMed Google Scholar * Laura Croci View author publications You can also search for this author inPubMed Google Scholar * Mark A. Landy View author publications You

can also search for this author inPubMed Google Scholar * Helen C. Lai View author publications You can also search for this author inPubMed Google Scholar * Gian Giacomo Consalez View

author publications You can also search for this author inPubMed Google Scholar * Antoine de Chevigny View author publications You can also search for this author inPubMed Google Scholar *

François Lallemend View author publications You can also search for this author inPubMed Google Scholar * Igor Adameyko View author publications You can also search for this author inPubMed 

Google Scholar * Saida Hadjab View author publications You can also search for this author inPubMed Google Scholar CONTRIBUTIONS Study design and supervision: S.H. Tissue collection and

acquisition of data: Y.W., M.E.K., P.F., K.K.Y.C, C.P., H.W., L.L.S., K.R., L.C., M.A.L., H.C.L., A.d.C. and S.H. Analysis and interpretation of data: L.F. and S.H., with inputs from M.E.K.,

F.L. and I.A. Figures: S.H. Drafting of manuscript: L.F., F.L., I.A. and S.H. with inputs from all co-authors. CORRESPONDING AUTHOR Correspondence to Saida Hadjab. ETHICS DECLARATIONS

COMPETING INTERESTS The authors declare no competing interests. ADDITIONAL INFORMATION PEER REVIEW INFORMATION _Nature Communications_ thanks the anonymous reviewer(s) for their contribution

to the peer review of this work. Peer reviewer reports are available. PUBLISHER’S NOTE Springer Nature remains neutral with regard to jurisdictional claims in published maps and

institutional affiliations. SUPPLEMENTARY INFORMATION SUPPLEMENTARY INFORMATION REPORTING SUMMARY PEER REVIEW FILE SOURCE DATA SOURCE DATA FILE RIGHTS AND PERMISSIONS OPEN ACCESS This

article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as

you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party

material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s

Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/. Reprints and permissions ABOUT THIS ARTICLE CITE THIS ARTICLE Faure, L., Wang, Y., Kastriti, M.E. _et al._

Single cell RNA sequencing identifies early diversity of sensory neurons forming via bi-potential intermediates. _Nat Commun_ 11, 4175 (2020). https://doi.org/10.1038/s41467-020-17929-4

Download citation * Received: 03 January 2020 * Accepted: 23 July 2020 * Published: 21 August 2020 * DOI: https://doi.org/10.1038/s41467-020-17929-4 SHARE THIS ARTICLE Anyone you share the

following link with will be able to read this content: Get shareable link Sorry, a shareable link is not currently available for this article. Copy to clipboard Provided by the Springer

Nature SharedIt content-sharing initiative