# A cell fitness selection model for neuronal survival during development

Yiqiao Wang, Haohao Wu, Paula Fontanet, Simone Codeluppi, Natalia Akkuratova, Charles Petitpré, Yongtao Xue-Franzén, Karen Niederreither, Anil Sharma, Fabio da Silva, et al.

► To cite this version:

Yiqiao Wang, Haohao Wu, Paula Fontanet, Simone Codeluppi, Natalia Akkuratova, et al. A cell fitness selection model for neuronal survival during development. Nature Communications, 2019, 10 (1), pp.4137. 10.1038/s41467-019-12119-3. hal-02299154

HAL Id: hal-02299154

https://hal.science/hal-02299154

Submitted on 9 Jan 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

ARTICLE

https://doi.org/10.1038/s41467-019-12119-3

OPEN

# A cell fitness selection model for neuronal survival during development

Yiqiao Wang<sup>1,11</sup>, Haohao Wu<sup>1,11</sup>, Paula Fontanet<sup>1</sup>, Simone Codeluppi<sup>2</sup>, Natalia Akkuratova<sup>3</sup>, Charles Petitpré<sup>1</sup>, Yongtao Xue-Franzén<sup>1</sup>, Karen Niederreither<sup>4</sup>, Anil Sharma<sup>1</sup>, Fabio Da Silva<sup>5</sup>, Glenda Comai<sup>6</sup>, Gulistan Agirman<sup>1</sup>, Domenico Palumberi<sup>1</sup>, Sten Linnarsson<sup>2</sup>, Igor Adameyko<sup>3,7</sup>, Aziz Moqrich<sup>8</sup>, Andreas Schedl<sup>5</sup>, Gioele La Manno<sup>9</sup>, Saida Hadjab<sup>1,12</sup> & François Lallemand<sup>1,10,12</sup>

Developmental cell death plays an important role in the construction of functional neural circuits. In vertebrates, the canonical view proposes a selection of the surviving neurons through stochastic competition for target-derived neurotrophic signals, implying an equal potential for neurons to compete. Here we show an alternative cell fitness selection of neurons that is defined by a specific neuronal heterogeneity code. Proprioceptive sensory neurons that will undergo cell death and those that will survive exhibit different molecular signatures that are regulated by retinoic acid and transcription factors, and are independent of the target and neurotrophins. These molecular features are genetically encoded, representing two distinct subgroups of neurons with contrasted functional maturation states and survival outcome. Thus, in this model, a heterogeneous code of intrinsic cell fitness in neighboring neurons provides differential competitive advantage resulting in the selection of cells with higher capacity to survive and functionally integrate into neural networks.

<footnote n="1">¹Department of Neuroscience, Karolinska Institutet, 17177 Stockholm, Sweden. ²Unit of Molecular Neurobiology, Department of Medical Biochemistry and Biophysics, Karolinska Institutet, 17177 Stockholm, Sweden. ³Department of Physiology and Pharmacology, Karolinska Institutet, Stockholm 17165, Sweden.</footnote>

<footnote n="4">⁴Institut de Génétique et de Biologie Moléculaire et Cellulaire (IGBMC), CNRS UMR7104, Inserm U964, Université de Strasbourg, Illkirch, France.</footnote>

<footnote n="5">⁵Université Côte d'Azur, Inserm, CNRS, iBV, 06108 Nice, France. ⁶Stem Cells & Development - Institut Pasteur - CNRS UMR3738, 75015 Paris, France.</footnote>

<footnote n="7">⁷Center for Brain Research, Medical University Vienna, Vienna, Austria. ⁸Aix-Marseille-Université, CNRS, Institut de Biologie du Développement de Marseille (IBDM), UMR 7288, 13288 Marseille, France. ⁹Brain Mind Institute, School of Life Sciences, Swiss Federal Institute of Technology (EPFL), Lausanne, Switzerland. ¹⁰Ming Wai Lau Centre for Reparative Medicine, Stockholm node, Karolinska Institutet, Stockholm, Sweden. ¹¹These authors contributed equally: Yiqiao Wang, Haohao Wu. ¹²These authors jointly supervised this work: Saida Hadjab, François Lallemand. Correspondence and requests for materials should be addressed to S.H. (email: saida.hadjab@ki.se) or to F.L. (email: francois.lallemand@ki.se)</footnote>
Programmed cell death plays a fundamental role in animal development and tissue homeostasis. In the mammalian nervous system, massive cell death and axonal degeneration represent an essential process in the formation of functional neural circuits during development¹,². The accepted mechanism driving the death of neurons during this critical period involves a stochastic competition between neurons for limiting access to target-derived neurotrophic factors³,⁴. Yet, how or which given cells are chosen to survive or die is still poorly understood.

In the vertebrate nervous system, the majority of neurons form during embryonic development when, after a phase of neurogenesis, a vast number of them are quickly removed by apoptosis in order to optimize quantitatively the connectivity between neurons and their targets (systems matching) and control the final architecture of the neuronal structures¹,³,⁴. In the peripheral nervous system, survival of selected neurons depends on target-derived neurotrophins (NTs)⁵,⁶, which bind to specific tropomyosin receptor kinase (TRK) receptors on the neuronal cell surface⁷. Typically, the binding of TRK receptors by their cognate NTs at the nerve endings induces a retrograde transport signaling that supports the survival, growth and differentiation of the neurons⁸,⁹. The neurotrophic theory is the prevailing model for the regulation of neuronal cell death during development³,⁴. This model proposes that neurons produced in excess compete for a limited supply of NTs synthesized in their target fields⁴,¹⁰⁻¹². Built into this model is also the assumption that the probability of individual neurons to survive the competition for NTs cannot be predicted¹⁰,¹³,¹⁴, implying that all neurons are endowed with an equal potential to compete. However, very early on, neurons exhibit intrinsic differences that control their differentiation and axon targeting¹⁵⁻¹⁷, emphasizing the need to re-examine the stochastic nature of the neurotrophic theory and to determine whether cell-to-cell molecular variability and survival predictability represent strong selection mechanisms during neuronal competition for NTs-induced survival (Fig. 1a). In this context, a more deterministic life-or-death regulation, coupled with these intrinsic differentiation mechanisms, would provide a way for selecting the fittest cells during development.

Here, we tested this hypothesis using tropomyosin receptor kinase C (TRKC)⁺ proprioceptive sensory neurons (PSNs) of the dorsal root ganglia (DRG) and demonstrate that the survival probability is not equal amongst neurons and hence some neurons are more likely to die. We show that before cell death period starts, the presumptive PSNs display different molecular blueprints initiated by the morphogen retinoic acid (RA) and which lead through transcriptional regulation to distinct subgroups of neurons with different functional fitness and survival outcome. Thus, our data demonstrate a predisposition to cell death during development, where molecular fitness heterogeneity in sensory neurons predicts their survival probability and the selection of the fittest cells for the establishment of functional neural networks.

## Results

### Heterogeneity of TRKC expression in PSNs before cell death.
To investigate whether neurons are equivalent prior to the cell death period, we used the DRG that have been extensively studied for neuronal cell death regulation during development. We focused on the development of the PSNs, which later innervate skeletal muscles and represent a major neuronal population of the DRG of the trunk¹⁸. During embryogenesis, PSNs are characterized by the expression of the transcription factors RUNX3 and of TRKC, the specific receptor for neurotrophin-3 (NT3) (Fig. 1b)¹⁷⁻¹⁹. At embryonic day (E) 11.5, most TRKC⁺ DRG neurons at brachial levels are RUNX3⁺ (and all RUNX3⁺ neurons are TRKC⁺), the vast majority of them being actual PSNs, as

assessed by genetic tracing of TRKC-expressing neurons from E11.5 to E17.5 using TrkCCreER;R26tdTOM mice (Fig. 1c and Supplementary Fig. 1)¹⁸,²⁰⁻²². On the other hand, NT3 is synthesized in the mesenchyme and the early muscle mass surrounding the growing axons of TRKC⁺ PSNs during development²³. Upon binding to TRKC, NT3 retrogradely activates pro-survival signaling events that are required for PSNs to survive the cell death period, which occurs from E12.0 to E12.5 and before target innervation (Fig. 1d)²⁴⁻²⁶.

The classic neurotrophic hypothesis considers that neurons have equal capacity to compete and thus, to survive the developmental cell death period, implying similarity in their molecular features before they reach their target¹³,¹⁴. However, at E11.0–11.5, analysis of the molecular identity of PSNs showed a high variability in their TRKC expression, with no obvious spatial pattern or a correlation with the cell size (Fig. 1e–i and Supplementary Fig. 2). This heterogeneity was maintained in vitro, with TRKC levels correlating between the soma and nerve endings where NT3 engages its receptor in vivo (Supplementary Fig. 3). The functional, full length (FL) tyrosine kinase isoform of *Ntrk3*⁷ was the predominant isoform (64.34 ± 1.03%) in E11.5 brachial DRG when analyzed by real-time quantitative polymerase chain reaction (PCR), and results from immunostainings and single-molecule RNA fluorescence in situ hybridization (smFISH) confirmed the high cell-to-cell variable expression of the *Ntrk3* FL isoform (Fig. 1j, k and Supplementary Fig. 4). Thus, our results clearly indicate that prior to cell death period, TRKC FL expression in individual PSNs is highly heterogeneous.

To examine whether the different levels of TRKC in PSNs arise from protein accumulation during neurogenesis (from E9.5 to E10.5 in mice)²⁷, with the early born neurons having accumulated more protein over time than later born counterparts, we fate-mapped early born TRKC neurons. For this, we induced recombination in TrkCCreER;R26tdTOM mice at ~E9.75 with a single injection of 4-hydroxytamoxifen (4-OHT, 0.06 g/kg). Convincingly, recombination in PSNs was not correlated with their TRKC levels observed at E11.5 as tdTomato (TOM)⁺ cells analyzed at this stage were similarly distributed among the TRKCHigh and TRKCLow categories of PSNs (Fig. 1l–n), indicating that TRKC heterogeneity is independent of birthdate.

### TRKC levels are associated with competitive advantage.
The high variability of TRKC levels among PSNs at early developmental stages suggests a target-independent control of their molecular heterogeneity. Here, we explored the possible influence of the environment on TRKC levels, focusing on NT3, which is expressed around the projecting axons in the limb²³. Using real time quantitative PCR and immunohistochemistry, we showed that TRKC expression is completely independent of NT3 both in vitro and in vivo (Fig. 2a–c); *Etv1* expression, known to depend on NT3 in vivo²⁸, showed, however, a 4.4-fold increase in NT3 condition (*P* < 0.001, unpaired Student's *t*-test). In line with that hypothesis, a previous work in embryonic chicken DRG showed that ablation of the limb does not alter TRKC expression in the corresponding DRG neurons prior to cell death, indicating that TRKC levels in PSNs are established independently of target-derived signals¹⁷. To study the transcriptional control of TRKC levels, we focused on the transcription factor RUNX3, which is known to transactivate the *Ntrk3* locus²⁹ (*Ntrk3* codes for TRKC), and showed a tight correlation with TRKC levels in E11.5 PSNs (Fig. 2d). Analysis of *Runx3*⁻¹⁻¹⁻ mice at E11.5 showed significantly lower levels of TRKC in PSNs (now similar to the TRKCLow category in Ctr animals), with however similar number of neurons (Fig. 2e, f). This confirmed earlier studies on the role of RUNX3 on TRKC levels but not on its induction¹⁷,¹⁹. In
<figure><img src="image_3.png" /><figcaption>**Fig. 1** Differential expression of TRKC in PSNs prior to the cell death period. **a** Scheme of our working hypothesis. **b, c** Temporal fate mapping of TRKC PSNs by 4-OHT induction. *TrkC<sup>CreER</sup>* mice allow temporary activation of CreER in the TRKC<sup>+</sup> cells 2 h after 4-OHT injection<sup>21,22</sup>. Immunostaining for PV, RFP and RUNX3 on E17.5 DRG sections (c) and graph showing distribution of PV<sup>+</sup>/RUNX3<sup>+</sup> PSNs among the TOM<sup>+</sup> cells (n = 4). Scale bar: 20 µm. **d** Quantification of PSNs at C5 and C7. ***P < 0.001, one-way analysis of variance (ANOVA) with Sidak's multiple comparisons test (n = 2-3). The window of PSNs cell death is shown. **e** TRKC expression in E11.5 ISL1<sup>+</sup> (and RUNX3<sup>+</sup>, whose staining is not shown for more visibility) DRG neurons. Scale bar: 50 µm. **f** TRKC levels in PSNs of **e** illustrated by color coding; dark blue indicates the lower and red the higher TRKC levels. From here, all observations are done at brachial levels (C5-8). **g** Distribution of TRKC levels in PSNs from **e. h** Distribution of TRKC levels in PSNs in E11.5 DRG neurons (from **g**). The data exhibit a Poisson-like distribution (one representative animal), with the mean used to define the two different categories of TRKC intensity (TRKC<sup>High</sup> and TRKC<sup>Low</sup>). **i** Projection of seven images of RUNX3<sup>+</sup>/TRKC<sup>+</sup> PSNs from one brachial DRG; dots indicate TRKC-labeled neurons and color codes reveal TRKC intensity as shown in **h. j** Projection image of smFISH for pan Ntrk3 and Ntrk3 full length (FL) transcripts in E11.5 DRG, visualized at high magnification in (1) and (2) (images show full projection); right panel shows color coding of Ntrk3 FL levels in red; the brighter, the higher levels. **k** Distribution of the number of Ntrk3 FL molecules in E11.5 DRG neurons by smFISH, normalized to pan Ntrk3 (Ntrk3 FL represent 68% of all Ntrk3 transcripts). **l** *TrkC<sup>CreER</sup>; R26<sup>tdTOM</sup>* mice were injected at E9.75 with 4-OHT and analyzed at E11.5 (n = 3). **m, n** Frequency distribution (**m**) and pie chart (**n**) of TOM<sup>+</sup>/TRKC<sup>+</sup> neurons from **l** according to their level of TRKC intensity. Source data are available as a Source Data file</figcaption></figure>

contrast, RUNX3 levels remained unchanged up to E12.5 in the absence of NT3-TRKC signaling (Supplementary Fig. 5), confirming that NT3 does not induce or regulate RUNX3 expression¹⁹. Collectively, our data support a model in which the variable levels of TRKC expression in PSNs are intrinsically regulated prior to cell death period in a target- and NT-independent manner.

Given the essential role of TRKC signaling in the survival of developing PSNs, we sought to analyze the possible relationship between TRKC levels in PSNs and their survival capacity. We found a positive correlation between the levels of TrkC and the expression of the phosphorylated form of the serine/threonine-specific protein kinase AKT, which acts within the TRKC signaling survival pathway during the cell death period⁷ (Fig. 2g),

suggesting more survival advantage for neurons with high relative TRKC levels. For our hypothesis to hold, lower TRKC-expressing PSNs would be more vulnerable to limited amount of NT3. We thus tested whether the levels of TRKC in PSNs impact their responsiveness to NT3 and hence their survival capacity. We took advantage of *Runx3<sup>-/-</sup>* mice, which have very low levels of TRKC in PSNs and in which all PSNs die at E12.5 (Fig. 2e, f)¹⁷,³⁰. We challenged PSNs from *Runx3<sup>-/-</sup>* DRG in vitro with different concentrations of NT3. Compared to wild-type neurons, PSNs from *Runx3<sup>-/-</sup>* mice did not respond to NT3-dependent survival at low concentration (1 ng/ml, EC50, Supplementary Fig. 6b) while a saturating concentration of NT3 (50 ng/ml) could promote survival, although to a lesser extent than in wild-type (Fig. 2h). This indicates that in the *Runx3<sup>-/-</sup>* animals (with low
<figure><img src="image_4.png" /><figcaption>**Fig. 2** Regulation of TRKC levels in PSNs and their correlation with differential survival. **a** Quantification of *Ntrk3* in E11.5 DRG cultured for 6 h with or without NT3 (50 ng/ml) (*n* = 4). **b** TRKC/RUNX3 immunostaining on E12.5 DRG sections from *Bax*<sup>-/-</sup> and *Bax*<sup>-/-</sup>;*NT3*<sup>-/-</sup> animals<sup>34</sup>. Scale bar: 50 µm. **c** Quantification of TRKC levels in individual PSNs from E12.5 *Bax*<sup>-/-</sup> and *Bax*<sup>-/-</sup>;*NT3*<sup>-/-</sup> animals in **b** (*n* = 2), with similar number of PSNs in both *Bax*<sup>-/-</sup> and *Bax*<sup>-/-</sup>;*NT3*<sup>-/-</sup> DRG (*P* = 0.83, unpaired Student's *t*-test). **d** Correlation between RUNX3 and TRKC levels in individual E11.5 PSNs. **e** TRKC expression in E11.5 DRG sections from *Runx3*<sup>+/+</sup> and *Runx3*<sup>-/-</sup> animals. Scale bar: 20 µm. **f** Quantification of TRKC intensity in individual PSNs in E11.5 *Runx3*<sup>+/+</sup> and *Runx3*<sup>-/-</sup> animals. **\*\*P* < 0.01. Bar graph: number of TRKC<sup>+</sup> neurons in E11.5 *Runx3*<sup>+/+</sup> and *Runx3*<sup>-/-</sup> animals; ns, not significant, unpaired Student's *t*-test (*n* = 3). **g** Immunostaining for TRKC, RUNX3 and p-AKT on E11.5 DRG sections and quantified correlation between TRKC and p-AKT levels in individual PSNs (right panel). **h** Survival of E11.5 DRG neurons from *Runx3*<sup>+/+</sup> and *Runx3*<sup>-/-</sup> animals treated with NT3 for 48 h (*\*\*\*P* < 0.05, *n* = 4–6, one-way analysis of variance (ANOVA) with Sidak's multiple comparisons test). **i** Scheme representing the in ovo electroporation (HHst 13/14, or E2). **j** Cross section of electroporated HHst26/27 embryos (E5) showing transfected cells in half spinal cord (SC) and in the ipsilateral DRG. **k** Quantification at E5 (before cell death period) of TRKC intensity among transfected PSNs (*RUNX3+/TRKC+*). \*\*\*\*\* *P* < 0.0001, unpaired Student's *t*-test. **l** Comparison before (E5) and after cell death period (E7) of the % of PSNs in transfected *GFP+* neurons. \*\*\* *P* < 0.01, unpaired Student's *t*-test (*n* = 4). **m** TRKC and RUNX3 expression in PSNs of E13.5 *Bax*<sup>+/+</sup> and *Bax*<sup>-/-</sup> DRG sections. (Right) Boxed regions (1) and (2) are magnified. Arrowheads show PSNs with lower size (see Supplementary Fig. 6c). Scale bar: 50 µm. **n** Quantification of TRKC levels in PSNs from E13.5 *Bax*<sup>+/+</sup> and *Bax*<sup>-/-</sup> DRG neurons (*n* = 4). Dashed line distinguishes TRKC<sup>High</sup> and TRKC<sup>Low</sup> populations. Source data are available as a Source Data file</figcaption></figure>

TRKC levels in PSNs), the NT3-TRKC survival signaling is still responsive but only to high concentrations of NT3. To further support a role for TRKC levels in predisposition of neurons to apoptosis during competition for NTs, we set out to reduce but not abolish the level of TRKC in PSNs by electroporating chick embryos at Hamburger and Hamilton stage 13 (HHst13; before neural crest migration, around E2) with siRNA either directed against *Ntrk3* and a GFP plasmid (*Ntrk3* siRNA) or with a negative control siRNA and GFP plasmid (*Ctr* siRNA) (Fig. 2i). This strategy allows fate tracing of neural crest-derived *Ntrk3* siRNA-expressing PSNs, which exhibit low levels of TRKC (Fig. 2j). It also enables the targeting of only a small proportion of the neurons (~7% of PSNs), hence avoiding possible indirect role on neurogenesis as previously shown in mice and chicken embryos in the absence of NT3-TRKC signaling<sup>23,31,32</sup>. Analysis

of DRG sections at HHst26/27 (E5), i.e. before the cell death period of PSNs<sup>33</sup>, revealed no noticeable loss of *GFP+* PSNs in *Ntrk3* siRNA condition (Fig. 2k). In contrast, at HHst30/31, after the cell death period of PSNs, there was an 80% decrease in the number of *GFP+* PSNs electroporated with *Ntrk3* siRNA, when compared with DRG sections from embryos transfected with *Ctr* siRNA (Fig. 2l). This suggests that low TRKC-expressing PSNs are less fitted to compete for NTs during the cell death period and predict a greater survival of neurons with higher TRKC levels. Inhibition of survival competition by means of cell death blockade using *Bax*<sup>-/-</sup> mice (in which all sensory neurons survive)<sup>34</sup> should then prevent this selection and result in a proportional increase of neurons with low TRKC. In *Bax*<sup>-/-</sup> DRG, the number of PSNs at E13.5 is about 2.2-fold higher (TRKC+/RUNX3+ neurons from C5–C8; wild-type (WT):
100 ± 13.7, Bax⁻¹⁻: 228 ± 12.6% of WT), reflecting a rate of cell survival that approximates 50% in wild-type PSNs. Consistent with our prediction, this increase was accompanied by an overrepresentation of neurons with low TRKC expression (Fig. 2m, n), a large proportion of them being atrophied (Supplementary Fig. 6c), suggesting that prevention of natural PSNs cell death results in the accumulation of neurons with relative lower TRKC content, which would have been eliminated in a wild-type context. Together, these data indicate a strong, positive relationship between TRKC levels, trophic signaling state and the capacity to compete and survive.

**Differential survival of PSN subpopulations in vivo.** To directly assess the higher survival probability of PSNs with high TRKC levels in vivo, we developed a Cre recombinase-dependent methodology to genetically trace PSNs as a function of TRKC expression. Indeed, by using the CreER driver and a low dose of 4-OHT, the fraction of cells undergoing recombination should depend on the amount of Cre recombinase they express³⁵. To test this in vivo, we induced limited recombination in *TrkCCreER*; R26tdTOM mice at E11.0 with a single injection of 4-OHT (0.018 g/kg) and analyzed the traced cells before (E11.5) and after the cell death period (E13.5) (Fig. 3a). At E11.5, recombination was observed in 11% of PSNs, predominantly labeling neurons

with high TRKC levels (with 75% of TOM⁺ PSNs expressing levels of TRKC above the average) (Fig. 3b-d), supporting our strategy to increase the probability of tracing high TRKC⁺/high Cre⁺ neurons by limiting the Cre activity. Importantly, this preferential tracing of the high TRKC neurons also demonstrate that the difference in TRKC levels in PSNs is relatively stable during this time period. Strikingly, when assessing the survival of the labeled PSNs after the cell death period (E13.5), we found an 1.88-fold increase in the proportion of labeled PSNs (from 11.6% at E11.5 to 20.9% at E13.5; Fig. 3e, f), confirming our prediction and supporting our hypothesis. In sharp contrast, injection of *TrkCCreER*; R26tdTOM mice at E14.0 (after the cell death period of PSNs) with a similar dose of 4-OHT led to an identical proportion of recombined PSNs at both E14.5 (13.2%) and E16.5 (13.5%) (Fig. 3g). Moreover, similar experiments using an R26CreERT2 driver line, which labels PSNs independently of their TRKC levels, did not show any enrichment in the proportion of recombined PSNs (TOM⁺ amongst TRKC⁺/RUNX3⁺ neurons) between E11.5 and E13.5 (Fig. 3g and Supplementary Fig. 7).

Altogether, our data demonstrate that neuronal competition selects for neurons with higher TRKC levels during development. Since NT3 expression in the limb does not seem limited to particular sub-regions of PSNs targets at E11.5²³, it is highly unlikely that this selection mechanism could reflect the possibility that PSNs with differential levels of TRKC would project

<figure><img src="image_5.png" /><figcaption>**Fig. 3** PSNs with high TRKC levels preferentially survive the cell death period. **a** Temporal fate mapping of TRKCHigh PSNs by 4-OHT (low dose, 0.02 g/kg). **b–d** Injection of *TrkCCreER*; R26tdTOM mice with low dose of 4-OHT at E11.0; DRG analyzed at E11.5 with recombination in few (b), preferentially high TRKC PSNs (c, d) (*P* < 0.001). Frequency distribution of TOM⁺ PSNs according to TRKC intensity (c) and pie charts (d) illustrating the large proportion of TOM⁺ cells among TRKCHigh PSNs. Scale bar: 50 µm. **e, f** Percentage of recombined PSNs at E11.5 and E13.5 in DRGs from *TrkCCreER*; R26tdTOM animals after 4-OHT injection at E11.0 (*P* < 0.05, Student's t-test; *n* = 2 litters with 6 embryos, E11.5; 2 litters with 5 embryos, E13.5). **g** The percentage of labeled PSNs does not change between E11.5 and E13.5 in R26CreERT2; R26tdTOM embryos injected at E11.0 with 0.032 g/kg 4-OHT (*n* = 4). Similarly, the recombination rate in TOM⁺ PSNs does not change between E14.5 and E16.5 in *TrkCCreER*; R26tdTOM embryos after injection at E14.0 with 0.02 g/kg 4-OHT (*n* = 2). Unpaired Student's t-test. **h** Whole-mount immunostaining for TRKC, NF160 and RFP of E13.5 forelimb from *TrkCCreER*; R26tdTOM embryos injected with low dose 4-OHT at E11.0. Insert shows restricted number of TOM⁺ fibers dispersed amongst TRKC⁺ axons. Scale bar: 200 µm. **i** Pattern and color-coded depth (in micrometers) of innervation of TRKC⁺, NF160⁺ and RFP⁺ nerve fibers (processed from **h**) of E13.5 forelimb from *TrkCCreER*; R26tdTOM embryos injected with a low dose of 4-OHT (0.02 g/kg) at E11. The pattern and depth color code reveal similar territories (in all dimensions, xyz) of innervation of the TOM⁺ PSNs compared to all axons (NF160). Scale bar: 200 µm. **j** Scheme illustrating the preferential selection of TRKCHigh PSNs during the cell death period. Source data are available as a Source Data file</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>b</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>c</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>d</figcaption></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /><figcaption>e</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>f</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>g</figcaption></figure>

<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /><figcaption>h</figcaption></figure>

<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /><figcaption>i</figcaption></figure>

<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /><figcaption>j</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>k</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>l</figcaption></figure>

<figure><img src="image_21.png" /></figure>

peripherally to distinct target fields, with distinct availability to
NTs. In support of this, defasciculation of sensory projections
within the limb in mutant mice for axon guidance molecules has
no impact on the neuronal cell death process³⁶,³⁷, and peripheral
projections of PSNs with high TRKC levels did not show any
preference of spatial distribution within peripheral nerves (Fig. 3h,
i). Altogether, this strongly suggests a selection model in which
neighboring neurons with differential fitness compete for NTs in

the target tissue, resulting in the elimination of cells with lower
advantage to compete and survive (Fig. 3j).

**Cell fitness characters of PSN subpopulations.** Our findings of the survival selection of one PSNs group (expressing high levels of TRKC and of RUNX3) over the other group suggest that these two groups might represent two genetically distinct PSNs
**Fig. 4** Single-cell RNA-seq identifies two molecularly distinct populations of PSNs prior to cell death period. **a** Workflow for single-cell RNA-seq (scRNA-seq) processing of DRG neurons. **b** Heat map showing single-cell expression of the top 20 differentially expressed genes in the different clusters of DRG neurons, identifying Clusters 1 and 2 as PSNs. **c** t-SNE plot showing two distinct populations of *Ntrk3+/Runx3+* PSNs at E11.5, the Clusters 1 and 2 in **b**. eTRKA: early TRKA neurons; MR: mechanoreceptive neurons¹⁸; both populations transiently express TRKC at an earlier time point³⁸. **d** List of Cluster 1-(top panel) and Cluster 2-enriched (bottom) genes. **e** The two clusters express similar levels of classic neuronal markers but differ in their level of *Runx3* and *Ntrk3* transcripts and are named *Runx3High* and *Runx3Low*. Box plots visualize the summary of the dataset (minimum, lower quartile, median, upper quartile and maximum). Wilcoxon rank sum test († indicates significant difference, see Methods section for details). **f, g** In vivo validation of the correlation between TRKC protein and *Runx3* mRNA (**f**) and between TRKC protein and *Ntrk3* mRNA (**g**). **\*\*P** < 0.01, unpaired Student's t-test (n = 3). **h-j** Marker genes of the *Runx3Low* and *Runx3High* subgroups of PSNs; t-SNE plots (**h**) and in vivo validation by immunostaining or in situ hybridization (**i**). Graphs in **i** show quantification of the expression of marker genes for *Runx3High* (High) and *Runx3Low* (Low) PSNs (n = 3 animals for each experiment), expressed in average number of dots for transcripts and average immunostaining intensity for proteins (and normalized to the maximum average value). **\*\*\*P** < 0.001, unpaired Student's t-test. Scale bar: 5 µm. **j** In situ hybridization on E11.5 DRG sections showing the (almost) mutually exclusive expression of *Fxyd7* and *Pou2f2*. Scale bar: 15 µm. **k** Quantitative illustration of **j**, where for each *TrkC⁺* neurons, their level of *Fxyd7* and *Pou2f2* is represented. **l** Schematic representation of neuronal types with their key markers and their average soma size (in square micrometers). Source data are available as a Source Data file

subgroups. To test this hypothesis, a total of 350 TOM⁺ cells from *TrkC⁺Cre;R26TOM* E11.5 DRG, in which the majority of the neurons generated between E9.5 and E10.5 are traced³⁸, were processed for single-cell transcriptome analysis (Fig. 4a). The cell expression data from *Ntrk3+/Runx3⁺* PSNs were clustered using R package SEURAT and visualized using bi-dimensional *t*-distributed stochastic neighbor embedding (*t*-SNE) (Fig. 4b, c). Interestingly, two clusters of PSNs were identified unbiasedly, and were associated with distinct levels of *Runx3*: a *Runx3High* group and a *Runx3Low* group (Fig. 4b–e and Supplementary Data 1). It is important to note that both subgroups are PSNs since their genetic tracing at E11.5 lead at E17.5 to 95% of labeled cells with a PSN phenotype (Fig. 1c and Supplementary Fig. 1) and also because they all can develop into PSNs when NT3 is overexpressed in developing muscles³⁹. These two PSNs subgroups were further distinguished based on their differential expression of marker genes (Fig. 4d). *Runx3High* group was characterized by the expression of *Hmx1*, *Fxyd7*, *Fam19a4* and an enrichment in *Ntrk3* while the *Runx3Low* group expressed *Sst*, *Pou2f2* and *Tbx3* (Fig. 4d, h). The correlation between *Runx3* mRNA level and TRKC protein, and between TRKC protein and *Ntrk3* mRNA, was then validated in vivo (Fig. 4f, g), as well as the existence of the two PSNs subgroups (Fig. 4h–k). Our results thus confirmed the existence of two PSNs subgroups based on their gene expression and that may correspond to the two PSNs groups with distinct survival outcome (Fig. 4l). We explored further our dataset and performed a gene ontology analysis to reveal benchmark biological categories behind the genetic blueprint of the *Runx3High* and *Runx3Low*. The top significant categories shared by both subgroups were biological categories associated with neuron development, axogenesis and neuron projection (Fig. 5a–c and Supplementary Fig. 8a). The two subgroups also showed similar levels of expression of genes involved in energy/metabolism (Supplementary Fig. 8b). Though, the two subgroups were embodied by distinct categories that underlined their transcriptome differences (Fig. 5d). Thus, the *Runx3High* group showed enrichment in neurotransmission and synaptic signaling categories (Fig. 5d, e), which could underlie a higher maturation of the functional properties of this subgroup. Notably, the levels of expression of the voltage-gated sodium channels *Scn3b* (*Na⁺β3*), *Scn7a* (*Na⁺*) and *Scn9a* (*Na⁺1.7*), which contribute to the generation of a mature action potential⁴⁰, and of *Piezo2*, the main mechanotransduction channels in non-noxious mechanosensitive sensory neurons⁴¹,⁴², were much higher in the *Runx3High* subgroup (Fig. 5f). To visualize the direction of developmental change within our population of cells, we applied RNA velocity to our single-cell measurements (PSNs from E11.5 and E12.5) to predict transcriptional dynamics and the cell state trajectory of

our neurons⁴³. We found similar clustering at E11.5, with the *Runx3Low* to the *Runx3High* subgroups, but only one homogeneous cluster of PSNs at E12.5. The RNA velocity analysis showed a strong directional flow from the *Runx3Low* to the *Runx3High* subgroups and then to the E12.5 stage (Fig. 5h–k), confirming that the subgroups of PSNs observed at E11.5 represent neurons of various cell states along their maturation trajectory. In support, marker genes of *Runx3Low* are not found in E12.5 PSNs, confirming also our genetic tracing data in Fig. 3 indicating that competition selects for TRKCHigh/RUNX3High PSNs, the other population being eliminated.

Altogether, those results reveal an unexpected genetic variability amongst developing PSNs that distinguishes two main subgroups prior to the cell death period (Fig. 5g). They also suggest that these transcriptomic differences could represent the genetic substrates defining both the survival potential of individual neurons during competition with neighboring cells and their functional, maturation states before the neurons acquire transduction competence⁴⁰.

**RA controls RUNX3 expression in PSNs.** The postmitotic induction of RUNX3 in PSNs prompted us to examine the potential extrinsic signaling pathways regulating RUNX3 expression in presumptive PSNs before the cell death period. Interestingly, the RA-synthetizing enzyme RALDH2 is found expressed around the nascent DRG during the early phase of DRG neuron generation (E9.5–E10.0) and precedes the appearance of RUNX3/TRKC PSNs observed at E10.5 (ref. ⁴⁴). To define the contribution of RA signaling to PSNs differentiation, we first used chicken DRG, which can clearly be distinguished and dissected at the early stages of DRG formation. Similar to mouse, RALDH2 enzyme was expressed in the mesenchyme surrounding the nascent DRG following a dorso-ventral gradient and in the motor neurons (MNs) area of HHst23 chicken embryos (Fig. 6a), before RUNX3 induction¹⁷,⁴⁵,⁴⁶. HHst23 DRG explants cultured for 12 h did not express RUNX3 (Fig. 6b, Ctr), even in the presence of NT3 (data similar to Ctr), while in contrast, all-trans RA induced RUNX3 expression in ISL1⁺ neurons (Fig. 6b, RA). Interestingly, increasing concentration of RA increased the levels of RUNX3 expression per individual cells (Fig. 6c). Thus, these results show that RA can induce a dose-dependent expression of RUNX3.

To study the role of RA on Runx3 expression in PSNs in vivo, we used *Raldh2⁻/⁻* embryos treated by a short maternal RA administration (E7–E9) to allow fetal development⁴⁷, creating thus an RA deficiency before RUNX3 onset. In line with the in vitro results, a significant proportion of TRKC⁺ neurons failed
<figure><img src="image_22.png" /></figure>

j

<table><caption>Cell survival model of sensory neurons during development.</caption><thead><tr><th>1-postmitotic PSNs fate</th><th>2-before cell death</th><th>3-cell death period</th><th>4-after cell death</th></tr></thead><tbody><tr><td>Extrinsic and Intrinsic dependent</td><td>Intrinsic differential survival and<br/>functional properties</td><td>Extrinsic and Intrinsic dependent:<br/>selection of the high expressing receptor cells (TrKC<sup>High</sup>)</td><td></td></tr></tbody></table>

**Fig. 6** Retinoic acid signaling regulates induction of RUNX3 expression in PSNs. **a** Section of the brachial region of an HHst23 chicken embryo showing expression of the RA-synthesizing enzyme RALDH2 at the time of physiological induction of RUNX3. Scale bar: 50 µm. **b** RUNX3 expression is induced in the presence of RA (500 ng/ml) in a subset of ISL1⁺ DRG neurons from HHst23 embryos within 12 h of treatment. Scale bar: 50 µm. **c** Concentration-dependent response of RA treatment on RUNX3 levels (n = 3). **d** Immunohistochemistry for ISL1, TRKC and RUNX3 on brachial DRG sections from E11.5 WT and *Raldh2<sup>-/-</sup>* mice. Scale bar: 20 µm. **e** Quantification of *d*, showing no difference in the number of ISL1 or TRKC positive neurons in *Raldh2<sup>-/-</sup>* mice (n = 3, unpaired Student's t-test). **f** Quantification of RUNX3 expression in TRKC⁺ neurons of *d* shows an overall decrease in intensity (*P* < 0.05, unpaired Student's t-test) and a significant *F* test (*F* = 0.0009) indicating a larger spreading of the variance between *Raldh2<sup>+/+</sup>* and *Raldh2<sup>-/-</sup>* mice (n = 3). **g** Quantification of TRKC expression in *d* shows significant decrease in its level of intensity per cell in the *Raldh2<sup>-/-</sup>* compared to *Raldh2<sup>+/+</sup>* mice (*** *P* < 0.001, unpaired Student's t-test). **h** Quantification of TRKC⁺/RUNX3⁺ DRG neurons on sections from E12.5 *Raldh2<sup>+/+</sup>* and *Raldh2<sup>-/-</sup>* mice shows a twofold decrease in the number of PSNs in the mutants (*** *P* < 0.01, unpaired Student's t-test, n = 4). **i** Number of PSNs is also largely decreased at E13.5 in the *Raldh1/2/3cKO* (*RCKO*) mutant mice compared to the control mice *Raldh1/2/3* (*Cre<sup>-/-</sup>*) (n = 3, *P* < 0.05, unpaired Student's t-test). **j** Proposed model for the selection of the TRKC DRG neurons during development, in which neurons, soon after being generated, feature different functional molecular signatures that are predictive of their survival or death during the cell death period (E12-E12.5). Source data are available as a Source Data file

to induce RUNX3 expression in *Raldh2<sup>-/-</sup>* mice (29%, n = 4, *P* = 0.005, t-test, Fig. 6d, f), while the total number of neurons and TRKC⁺ cells remained unchanged (Fig. 6e). Consistent with the role of RUNX3 on TRKC levels, presumptive PSNs that failed to acquire RUNX3 expression in *Raldh2<sup>-/-</sup>* mice showed a significantly lower expression of TRKC per cell than PSNs in *Raldh2<sup>+/+</sup>* mice (Fig. 6d, g). Notably, analysis of *Raldh2<sup>-/-</sup>* mice after the cell death period (E12.5) showed an almost twofold decrease in the number of PSNs compared to *Raldh2<sup>+/+</sup>* mice (Fig. 6h). To prevent possible compensatory mechanisms for the loss of RALDH2 by ectopic expression of RALDH1 and

RALDH3, we used a mutant mouse model in which the three members of the *Raldh* family are deleted conditionally: *CAGG<sup>CreER</sup>*; *Raldh1fl/fl*; *Raldh2fl/fl*; *Raldh3fl/fl*, hereafter named *Raldh1/2/3cKO*. When induced by tamoxifen at E10, *Raldh1/2/3cKO* mice do not need RA supplementation to survive⁴⁸. However, this timing of induction does not permit the analysis at the brachial level, where RUNX3 expression is observed from E10.5. Since RUNX3 induction follows the rostral-to-caudal gradient of normal DRG development, we characterized the phenotype at lumbar levels, where RUNX3 is first observed at E11.5. Similarly to the rescued *Raldh2<sup>-/-</sup>* mice, E13.5 *Raldh1/2/*
$^3$cKO mutant mice showed a large decrease in the number of TRKC$^+$/RUNX3$^+$ neurons (47%, Fig. 6i). Altogether, these results provide evidence that RA signaling in vivo is sufficient to induce RUNX3 in presumptive PSNs. The observation of some RUNX3 induction remaining in both mutant models could be due either to compensatory mechanisms, through activation of parallel signaling, or most likely, to an incomplete reduction of RA signaling during the critical period of RUNX3 induction using our strategies (incomplete efficiency of the recombination or of the clear-out of the RA supplement). Importantly, our results suggest a model in which the instructive role of RA on PSNs differentiation might eventually impact their differential capacities to interact with and integrate environmental cues for their survival (Fig. 6j). It is worth noting that although distinct levels of trophic factor receptors between neurons might participate in the fitness characteristic of the cell, other differentially expressed genes could be contributing to preferential cell survival.

Discussion

The prevailing view explaining cell death mechanisms in the
developing mammalian nervous system is mostly built on the
principle of a one-sided control of cell death by the environment
and a stochastic selection of the surviving neurons. This view has
been supported indirectly by observations that experimental
manipulations of neurons' environment impact their survival<sup>2,3</sup>.
Our findings support an alternative model in which sensory
neurons are endowed with distinct molecular features that are
predictive of their probability to respond to environment-derived
survival signals and thus, to survive the cell death period. We find
that prior to cell death, sensory neurons exhibit distinct func-
tional and survival states that are genetically encoded, involve RA
signaling and RUNX3 transcriptional activity. These states are
independent of signals derived from the target tissue and of NT
signaling and eventually lead to the survival selection of a defined
neuronal population. Thus, an early integration of cellular sig-
naling and transcriptional programs driving cell fate decisions are
coupled to a genetically distinct functional and survival fitness
and thus life-death fate outcomes of neurons during competition
for NTs (Fig. 6j).

Past studies have on the developmental cell death of neurons in
vertebrates have assumed that all TRK-positive sensory neurons
are initially equivalent and that their individual probability to
survive the competition period is unpredictable$^{4,13,14}$. The selec-
tive nature of the survival-or-death choice of PSNs during
development that we now reveal contrasts with this principle. Our
results show that two subgroups of PSNs display distinct intrinsic
capacities to compete and thus to survive prior to the cell death
period. Strikingly, our data show that the molecular fitness
characteristics of sensory neurons emerge very early, i.e. soon
after cell cycle exit, and before the selection of nerve trajectories.
Moreover, cell death of sensory neurons is a relatively fast process
(~12 h) and occurs days before the nerve endings contact their
final target$^{26}$. Thus, contrary to the classic competition model, in
which the size of the end organ target field would dictate that of
the neuronal population by competition for target-derived
trophic factors, the final set of sensory neurons, which is
defined by the cell death program, would determine the number
of peripheral sensory organs being innervated. In this context, a
more selective mechanism of cell death might provide an
advantage to guaranteeing a rapid clearance process of neurons
during early development. The systems matching principle has
been proposed to endow the nervous system with the capacity to
adapt quantitatively to environment changes. It has remained
unclear, however, how neurons are selected during this devel-
opmental period, especially if neurons had equal competitiveness

(stochastic selection) and if the neural network has not fully
developed yet. In our study, the distinct genetic features of the
two neuronal subgroups reveal molecular differences that are
associated with their competitive advantage, and hence defining a
fitness attribute to survive the cell death period. The advantage of
this survival selection process seems to lie within the difference in
the functional state of these two populations of neurons, which
co-varies with the survival advantage status of the cells. Indeed,
the fittest neurons displayed a more mature functional phenotype,
as illustrated by the large representation of genes involved in
neurotransmission in that group. Such profile is consistent with
their acquisition of mature biophysical properties and mechan-
osensitivity during this period, while concomitantly number of
neurons with immature action potentials decreases dramati-
cally⁴⁰. Interestingly, voltage-gated sodium channels, whose high
density is essential for the generation and propagation of action
potentials, were highly enriched in the fittest neuronal subgroup,
as was Piezo2 which encodes for the major mechanotransduction
channel PIEZO2 in mechanosensitive sensory neurons⁴¹,⁴².
Hence, our data reveal a selection model during nervous system
development where the competitive advantage of a neuron and its
functional state are consubstantial elements of a general fitness
profile that is established by specific transcriptional programs.
The result of this is the emergence of two distinct subgroups of
neurons, a more fitted and a less fitted, the former out-competing
the latter, which is selectively removed during competition for
survival signals. The mechanisms responsible for the differential
acquisition of these molecular properties by early sensory neurons
remains, however, unclear. Interestingly, those are linked to the
maturation trajectory of PSNs, yet our data show that the timing
of neurogenesis would not contribute significantly to these dif-
ferences. Provided that RA is presumably symmetrically dis-
tributed around the ganglion, stochastic process might play a role
in exposing early born neurons to various levels of RA or some
other molecular signaling pathways that would regulate the levels
of RUNX3 and other genes that comprise the two molecularly
distinct clusters of neurons.

Thus, in addition to scaling the final number of neurons to match the size of the targets, cell death in developing sensory neurons might also participate in scaling their physiological fitness. The genetic distinction of populations of neurons with various intrinsic molecular settings and fitness observed in our study could serve as a more general mechanism preceding the selection step (through cell death), which would eventually favor the most suitable cells for the establishment of functional neural networks.

Our study also suggests similarities with the developmental cell
death mechanisms in *Caenorhabditis elegans* where neuron pre-
cursors feature complex death codes that predict their survival or
death⁴⁹. A particularity of the development of *C. elegans* is,
however, its invariant character, with 105 cells from the neuronal
lineage being programmed to die at specific places and times⁴⁹,
underlying the importance of lineage programming in this
organism. In this context, it remains to be determined in verte-
brates whether the predisposition of some neurons to survive or
die during development is also lineage dependent or acquired
postmitotically.

Sensory neurons of the trunk derive from neural crest stem
cells which delaminate from the dorsal neural tube, with some
migrating ventrally and settling between the neural tube and the
dermomyotome to generate the DRG⁵⁰. Later during develop-
ment, these neurons will connect to peripheral organs or tissue
and to specific neurons in the spinal cord¹⁸. In the case of PSNs,
their peripheral projections reach their muscle targets by E14.0,
while their central projections contact MNs in the ventral spinal
cord at about E17.5⁵¹,⁵². The establishment of these networks
occurs several days after the cell death period in the DRG. How then is the development of each of these spatially disconnected structures coordinated for quantitative systems matching? These three structures (MNs, sensory neurons and skeletal muscles) develop simultaneously from adjacent precursor regions⁵⁰,⁵³,⁵⁴. Hence, the use of a common source of signaling molecules acting locally could participate in their synchronous development. RA, which is produced by the paraxial mesoderm⁴⁴–⁴⁶,⁵⁵, has been shown to control the cellular specification and number of limb muscle cells and MNs⁵⁴⁻⁵⁷. Our results provide evidence that RA regulates the specification of PSNs cell fate, which might impact their final number. Therefore, the paraxial mesoderm, through RA production, could serve as a signaling center for the coordinated development of sensory neurons, MNs and skeletal muscles and play an important part in the establishment of a future functional sensory-motor circuit.

Finally, cell-based therapies for neuronal replacement after cell loss have great potential for treating neurological disorders⁵⁸,⁵⁹. For this, in vitro culture systems are used to derive specific neuron types from stem cells, which are then transplanted locally to disperse in the targeted neural tissue. A limitation of this therapeutic strategy is, however, to optimize the survival and the integration of the transplanted (newly specified) neurons into functional networks⁵⁹,⁶⁰. Recent studies and ours suggest a heterogeneity in newly developed neurons, which affects their fitness and would eventually impact their survival and functional integration⁶⁰,⁶¹. Therefore, a better understanding of the distinct molecular features of individual developing neurons or stem-cell-derived neurons and their correlation with in vivo cell survival and function might provide a critical step in ensuring optimal and reproducible cell preparations for successful transplantation.

## Methods

**Animals.** Wild-type C57BL6 mice were used unless specified otherwise. *Runx3*⁻¹/⁻¹, *TrkC*CreER, *Raldh2*⁻¹/⁻¹, *Raldh1*fl/fl, *Raldh2*fl/fl, *Raldh3*fl/fl and the *CAGG*CreERTM mouse strains have been described elsewhere³⁰,⁴⁸,⁶²–⁶⁵. *Bax*⁻¹/⁻¹, *R26*CreERT2, *TrkC*⁻¹/⁻¹ and *Ai14* mice have been purchased from Jackson Laboratories, and *NT3*⁻¹/⁻¹ and *TrkC*Cre from MMRRC. 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.

Fertile white Leghorn eggs were incubated at 38 °C and embryos were staged according to Hamburger–Hamilton (HH) tables.

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ämd*.

**Treatments.** Rescued RA animals. For the RA-rescue experiments, pregnant females were treated with RA (all-*trans*-RA, Sigma) in powdered food between E7 to E9 at a dosage of 100 µg/g of food.

Conditional KO with tamoxifen treatment. Pregnant mothers of *CAGG*CreEM, *Radh1*fl/lox/lox; *Raldh2*fl/lox/lox; *Raldh3*fl/lox/lox embryos (referred as *Raldh1/2/3cKO*) were administered by gavage a dose of tamoxifen (0.1 mg per g of b.w.; Sigma) dissolved in corn oil (25 mg/ml).

Cell fate tracing experiments with 4-OHT. *R26*tdTom mice were crossed to the *TrkC*CreER and to *R26*CreERT2 mice for cell fate tracing of the PSNs and for timing of early and late born neurons from the first wave of neurogenesis in the DRGs. 4-OHT induction of the *TrkC*CreER mouse line was performed by single intraperitoneal injection of either low doses (0.018 g or 0.032 g per kg bw) or higher doses (0.06 g or 0.08 per kg bw) depending on the hypothesis tested (4-OHT, Sigma, H6278) into pregnant mothers at the indicated stages. 4-OHT was dissolved in 99.5% ethanol and further diluted 1:1 in Cremophor® EL (Sigma C5135) to obtain a concentration of 20 mg/ml for storage at −20 °C. 4-OHT stock solution was further diluted 1:3 in DPBS to a working concentration of 5 mg/ml upon usage. Analysis was performed on 2–8 animals from 2–3 separate experiments.

**In ovo electroporation.** A control plasmid pCAGeGFP and siRNAs were injected in ovo into the neural tube of HHst13/14 chick embryos (plasmids concentration: 1 µg/µl, siRNA concentration: 0.07 nM, Ambion). Five pulses of 50 V/cm was performed using a square wave electroporator (BTX). Embryos were collected 3 or 5 days post electroporation and processed for immunostaining.

**In vitro cultures of DRG neurons/whole DRGs.** Dissociated DRG neurons were obtained after enzymatic digestion of brachial DRG from E11.5 embryos using trypsin (0.05% Trypsin-EDTA; Gibco). Cells were plated in 24-well plates precoated with poly-p-lysine (0.01%; Sigma)/laminin (10 µg/ml; R&D Systems). Cells or whole DRG were cultured in N2 medium (DMEM-F12/glutamax medium with N2 supplement; Gibco) supplemented with pen/strep, gentamicin and with NT3 (Peprotech) when specified. For in vitro RUNX3 induction in whole chicken DRGs, NT3 (10 ng/ml, R&D systems), RA (500 ng/ml, Sigma) and the pan-caspase inhibitor Q-VD-OPh (10 µM, Sigma) were used. For the survival assay using mouse E11.5 DRG neurons, the majority of neurons die during the first 24 h in the absence of NTs⁶⁶, and the addition of NT3 promotes selectively the survival of PSNs⁶⁷. For transcripts (qPCR) measurements, whole DRGs were cultured in 2-ml open tubes in closed Petri dishes for the indicated time.

**Immunostainings and RNAscope® in situ hybridization.** Animals were collected, decapitated and fixed for 1–4 h at +4 °C (4% PFA in PBS) depending on the stage, washed in PBS, cryopreserved in 30% sucrose in PBS, embedded in OCT (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-*Fyxd7* and Mm-*Pou2f2* (ACDBio). Dissociated cultures were fixed for 5–20 min in cold PFA (4% in PBS). Sections or dissociated cultures were incubated for 24 h at +4 °C with primary antibodies diluted in blocking solution (2% donkey serum, 0.0125% NaN₃, 0.5% Triton X-100 in PBS). Primary antibodies used were: rabbit anti-RUNX3 (gift from Jessell TM), 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), rabbit anti-pAKT (1:100, Cell Signaling 4060S), chicken anti-RFP (1:250, Rockland 600-901-379S). After washing with PBS, Alexa Fluor secondary antibodies (Live Technology; 1:500 in blocking solution) were applied overnight (at +4 °C). Phalloidin conjugated with Alexa 488 was purchased from Life Technology. 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 neuronal counting, soma area and intensity quantifications ImageJ software was used. Only neurons with a visible nucleus were used. Quantification of molecular markers in the dora root ganglia was carried out on a minimum of five DRGs sections or more randomly selected sections per animal, per condition (see figure legends for n's and genotypes). In Fig. 1a, for accurate quantification and comparison of the number of PSNs across stages, we quantified the total number of PSNs in every third sections of C5 and C7 DRGs. Also, note that for the βIII-tubulin quantification, the two outer lateral rows of cells in the different DRGs sections were not included in the analysis due to saturation of signal caused by axonal track passage that are heavily stained.

**Cell intensity quantification.** Quantification of the intensity of protein of interest from DRG sections from mice was carried out by manually drawing the cell area of interest (defining the regions of interest (ROIs)); soma, nuclei and growth cone depending on the experiment. To define RUNX3 intensity in the RUNX3⁺/TRKС⁺ neuron, for a given RUNX3⁺/TRKС⁺ cell (exp_cell), the raw integrated density of a RUNX3⁺/TRKС⁻ cell (control_cell) (representing the background intensity), normalized to its area was subtracted from the raw integrated density of the RUNX3-positive cell (exp_cell), normalized to its area, as follows: RUNX3 intensity (exp_cell) = [RawIntDensity (exp_cell) / Area (exp_cell)] – average [RawIntDensity (control_cell) / Area (control_cell)]. Reads were averaged by section and plotted by genotype.

**Statistics.** Data were analyzed using GraphPad Prism 5 and expressed as mean ± S.E.M., except for dot plots showing S.D. 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. Data distribution was assumed to be normal. No animals or data points were excluded from the analyses. No statistical methods were used to pre-determine sample size but our sample sizes are similar to those generally employed in the field. Note also that each measurement was taken from distinct samples.

**qPCR.** Brachial level DRGs were dissected from six E11.5 embryos. DRGs were arranged so that each biological replicate contained five DRG (C4–C8) from one embryo, and the same from another embryo. Experimental conditions were Zero (DRG placed directly into lysis buffer after dissection), and +NT3/−NT3 (DRG cultured in for 6 h in with or without recombinant human NT3 (Peprotech, cat. no. 450-03). DRG were cultured in DMEM:F12 (Gibco) + 10% FBS (Gibco) supplemented with 5 µM of caspase inhibitor Q-VD-OPh (Calbiochem, cat. no. 551476) in 2-ml tubes. After either dissection (Zero) or culture (+NT3/−NT3), RNA was extracted using the Qiagen RNeasy Mini Kit according to the manufacturer's instructions, including DNase I to degrade potential genomic DNA contamination. RNA quality was quantified using the
Biorad Experion Standard Sensitivity RNA Kit, which showed the RNA had RINs between 9.7 and 10. RNA was quantified using the Qubit RNA BR Assay Kit, and 50 ng of each biological replicate used for reverse transcription using Biorad iScript in a 20-µl reaction according to manufacturer's instructions. No reverse transcriptase controls were carried out at the same time by substituting water for the reverse transcriptase. Resulting cDNAs were used to analyze transcript levels using real-time PCR in a Biorad CFX96 Real-time System and using Biorad iTaq Universal SYBR® Green Supermix in 20-µl reactions for 40 cycles (95 °C 10 s, 60 °C 20 s, 72 °C 30 s followed by fluorescence measurement at 75 °C to reduce the influence of primer dimer on quantification). Primers used (5' and 3'): *B2m*, GGT CGC TTC AGT CGT CAG and TC AGT ATG TTC GGC TTC CC; *Etv1*, TAC AAG AAA CAT GGC TTG CTG AAG and CAT GAA AAG CCA AAC TTT CAG CCT; *Gapdh*, AAC TCC CAC TCT TCC ACC TTC and GAT AGG GCC TCT CTT GCT CAG; *Hprt*, GAA TCT GCA AAT ACG AGG AGT CCT and CTT TAC TAG GCA GAT GGC CAC A; *Mapk6*, TGT TGT GCA CGA GGA GGA TC and TGT GTG GCT CAG TCA TGT CG, *Ntrk3* PAN, AGT AAC CGG CTC ACC ACA CTC and AGC GGA TGT CAC AGC TGC AGT; *Ntrk3* FL, TGA TCC TCG TGG ATG GAC AG and CTT CAC TAG TAG ATT GGC TCC. *Mapk6* and *Gapdh* were included as reference genes. No template and no reverse transcriptase controls did not amplify, or were detected at least five cycles after the least strong biological sample. In addition, amplicons were sequenced to confirm specificity. Data were analyzed using qBase+ (Biogazelle, Belgium).

The truncated isoforms of *Ntrk3* cannot be uniquely measured by qPCR method. Hence, the quantity of *Ntrk3* T transcripts is deduced by subtracting the value of *Ntrk3* FL from the value of pan *Ntrk3*.

**smFISH.** Ten-micrometer-thick sections were mounted on cover glasses. Fresh sections selected for the experiment were post-fixed with 4% PFA (in PBS) for 10 min at room temperature, rinsed with PBS and permeabilized with methanol. After permeabilization, the sections were incubated for 10 min at 70 °C in Tris-EDTA pH 8.0. After washing twice with SSC 2x, the sections were incubated with hybridization buffer containing 250 nM fluorescent label probes (LGC Biosearch Technologies; see Supplementary Table 1) for 4 h at 38.5 °C. After four 20% formamide-SSC 2x washes, the slides were counterstained with Hoechst, washed with SSC 2x and mounted with Prolong Gold mounting medium (ThermoFisher Scientific). Image stacks (0.3 µm distance) were acquired using a Nikon Ti-E with motorized stage (Nikon).

The images were analyzed with a custom python script using the numpy, scipy, ndimage (http://www.scipy.org/) and scikit-image libraries⁶⁸. Briefly, after background removal using a large kernel Gaussian filter, a Laplacian-of-Gaussian was used to enhance the RNA dots in the maximum projected image. Background objects significantly larger than the smFISH dots were removed after image thresholding and the remaining RNA dots counted in user-defined ROIs. The images were stitched, aligned and pseudocolored in Fiji. Pseudocolored images were generated by mapping the number of RNA molecules counted in each ROI to the value of a HSV color scale. For each image, a color value of zero corresponds to zero RNA molecules and the maximum color value is used to represent the highest molecular count.

Arlecchino image: after manual segmentation of the cells, the signal intensity was measured using imageJ. The signal intensity was binned in six expression classes and each class was identified with a different hue. The green color represents the lowest expression level and cyan the highest.

The truncated isoforms of *Ntrk3* cannot be uniquely labeled using the smFISH method. Hence, the number of *Ntrk3* T molecules in Supplementary Fig. 4a, b is deduced by subtracting, in individual cells, the number of *Ntrk3* FL from the number of pan *Ntrk3*.

**Single-cell isolation for single-cell analyses.** Brachial DRGs were dissected and collected in Leibovit'z L-15 medium (Life technologies) on ice. Then the DRGs were incubated in 0.05% trypsin-EDTA (1x) (Life technologies) for 5, 7, 10 and 15 min for DRGs 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). DRGs were 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.

**Fluorescence-activated cell sorting (FACS).** Single Tomato-positive cells were sorted by FACS into individual wells of 384-well plates containing lysis buffer. The plates containing single cells were frozen immediately on dry ice and stored at -80 °C. We used six animals (ten DRGs per animal) per experiment.

**Single-cell RNA sequencing.**

* Smart-Seq2 protocol was performed on single isolated cells by Eukaryotic Single Cell Genomics Facility at SciLifeLab, Stockholm⁶⁹. Raw data processing:

the raw data of single-cell RNA sequencing was processed with standard pipeline by Eukaryotic Single Cell Genomics Facility at SciLifeLab, Stockholm. Most of the downstream analysis was done using R software package Seurat v2.0.

* Pre-processing (Initial filtration, normalization):

In our initial analysis of the data, we noticed contamination of non-sensory neurons in each dataset. But sensory neurons could be easily distinguished by their high expression of *Pou4f1*. Thus, we only selected cells with *Pou4f1* transcript count over 5 for further analysis.

Considering the distribution of genes identified in every cell, we selected cells with at least 4000 unique genes detected for downstream analysis. This step removed empty wells and low-quality cells. Genes expressed in less than three cells were also removed from analysis.

We used a global normalization strategy incorporated in SEURAT package to make the data comparable between different cells and between different SMART-seq2 runs. The gene expression measurements within each cell were scaled by a constant factor 10,000, then natural-log transformed. This generated a new gene expression matrix $y_{i,t} = \log (x_{i,j} / Z_j \times 10,000)$, where $x_{i,j}$ is the count of gene *i* in cell *j* and $Z_j$ is the total counts of all genes in cell *j*.

We calculated the dispersion (ratio of variation to mean) of genes, excluding those lowly expressed on average (normalized value <0.0125) or less variable (variation <0.5). The remaining 4008 variable genes were used for dimensional reduction. We also regressed out the cell-cell variation in gene expression driven by percent.mito (the percentage of mitochondrial gene content) using ScaleData function.

* Linear dimensional reduction:

We used principal component analysis (PCA) to reduce the dimensionality of the data. Twenty PCs were computed and first 15 PCs were determined as statistically significant through the JackStraw method. We then applied *t*-SNE to visualize the data in a two-dimensional space.

* Identify differentially expressed genes:

To identify differentially expressed genes in each cluster, we used likelihood-ratio test implemented in function FindMarkers. The parameters min.pct = 0.25 and thresh.use = 0.1 are chosen to select genes expressed in more than 25% of cells and with average expression level over 0.1 (normalized data) in the given cluster. The genes were ranked by *P*-value from low to high, with lower *P*-value indicating more differentially expressed genes.

Boxplots: For boxplots shown in Supplementary Fig. 8, the lower and upper hinges represent the first and third quartiles, respectively, so the box spans the inter-quartile range (IQR). The horizontal line inside the box corresponds to median. The upper whiskers extend to the largest value no further than 1.5 × IQR from the upper hinge, while the lower whiskers extend to the smallest value at most 1.5 × IQR from the lower hinge. Dots represent the outliers beyond the range of whiskers.

**Gene Set Enrichment Analysis visualized by network.** All differentially expressed genes in *Runx3*<sup>flight</sup> and *Runx3*<sup>Low</sup> population were used as input for Gene Set Enrichment Analysis (GSEA) using Cytoscape and its plugin Bingo. The most up-to-date gene ontology file go-basicobo and gene association file gene-association.mgi were downloaded directly from Gene Ontology Consortium website instead of using Bingo's default ontologies or associations, which is not updated frequently. The following evidence codes were discarded from analysis: IEA (inferred from electronic annotation), ISS (inferred from sequence similarity).

Significant GO terms (*P*-value < 0.05) were used as input for Cytoscape's plugin EnrichmentMap to generate a network where mutually overlapping gene sets cluster together. Following parameters were used: *P*-value cutoff 0.001, FDR *Q*-value cutoff 0.05, similarity cutoff with Jaccard coefficient 0.25.

**RNA velocity.** We performed a pre-analysis that included the entirety of cells obtained by our sorting strategy; the aim was to isolate the *Runx3+/Ntrk3+* population in an unbiased fashion. This choice was motivated by our intention to avoid explicit filtering based on the expression of those genes, a procedure that could have excluded different part of the cells of interest (because of sequencing drop-out events) or included unrelated cells (because of cross-contamination events).

As a first preprocessing step, we excluded cells with <150,000 total reads counts and red blood cells, these were easily recognized by their low number of detected genes or the extremely high hemoglobin read counts. Furthermore, we removed genes that were not reliably detected, our condition for retention was: a gene needs to have at least 38 read counts across the whole dataset and at 18 counts in at least one cell. To select the most informative genes for downstream analysis, we applied a feature selection procedure⁷⁰. Briefly, a noise model was fit on the data and genes are ranked by the residuals to the fit, the top 1400 genes were selected. We log-transformed the data matrix and we generated a 20 nearest-neighbors graph using Euclidean distance in the space spanned by the first six principal components. This graph was, then, imputed to the Louvain community detection algorithm to cluster the cells.

For the core analysis on the *Runx3+/Ntrk3+* populations (*n* = 274), we used a QC and filtering procedures that identical to the one described above with the following
differences: that we selected more genes (1900 genes), more PCA components (65) and 25 nearest neighbors to construct the graph and perform Louvain clustering. These choices were aimed at extracting the most accurate representation of the dataset. To visualize the data, we used UMAP with parameters n_neighbors = 40, learning_rate = 0.6, min_dist = 0.5 and spread = 0.55. RNA velocity analysis was then performed using our tool *velocyto.py*<sup>43</sup>. We used the default parameters that we used in the original publication<sup>43</sup> with three differences:

* We pre-processed the data matrix only minimally (since the data is deeply sequenced and high quality), by calling *velocyto.knn_imputation* with parameters k = 5, b_sight = 50 and b_maxl = 20

* Only the most mature cell population (Cluster 3 in Fig. 5h) consisting entirely of cells sampled at the time point E12.5 was used to fit the degradation rate parameter (i.e. gamma). This procedure corresponds to the assumption that these cells are less rapidly changing and closer to a steady state.

* We adapted the parameter analysis to the number of cells in this dataset, in particular the *estimate_transitionProb* function was passed the arguments *n_neighbors* = 90 and *calculate_grid_arrows* the arguments smooth = 0.4, steps = (20, 20), *n_neighbors* = 40.

**Reporting summary.** Further information on research design is available in the Nature Research Reporting Summary linked to this article.

## Data availability

The accession number for the data used in this article is GEO: GSE135181. Data that support the findings of this study are also available from the corresponding authors upon reasonable request. Source data underlying Figs. 1-6 and Supplementary Figs. 2-7 are available as a Source Data file.

Received: 12 October 2018 Accepted: 16 August 2019

Published online: 12 September 2019

## References

1. Yamaguchi, Y. & Miura, M. Programmed cell death in neurodevelopment. *Dev. Cell* **32**, 478–490 (2015).

2. Oppenheim, R. W. Cell death during development of the nervous system. *Annu. Rev. Neurosci.* **14**, 453–501 (1991).

3. Buss, R. R., Sun, W. & Oppenheim, R. W. Adaptive roles of programmed cell death during nervous system development. *Annu. Rev. Neurosci.* **29**, 1–35 (2006).

4. Davies, A. M. Regulation of neuronal survival and death by extracellular signals during development. *EMBO J.* **22**, 2537–2545 (2003).

5. Cowan, W. M. Viktor Hamburger and Rita Levi-Montalcini: the path to the discovery of nerve growth factor. *Annu. Rev. Neurosci.* **24**, 551–600 (2001).

6. Levi-Montalcini, R. Developmental neurobiology and the natural history of nerve growth factor. *Annu. Rev. Neurosci.* **5**, 341–362 (1982).

7. Segal, R. A. Selectivity in neurotrophin signaling: theme and variations. *Annu. Rev. Neurosci.* **26**, 299–330 (2003).

8. Huang, E. J. & Reichardt, L. F. Trk receptors: roles in neuronal signal transduction. *Annu. Rev. Biochem* **72**, 609–642 (2003).

9. da Silva, S. & Wang, F. Retrograde neural circuit specification by target-derived neurotrophins and growth factors. *Curr. Opin. Neurobiol.* **21**, 61–67 (2011).

10. Levi-Montalcini, R. & Angeletti, P. U. Nerve growth factor. *Physiol. Rev.* **48**, 534–569 (1968).

11. Thoenen, H. & Barde, Y. A. Physiology of nerve growth factor. *Physiol. Rev.* **60**, 1284–1335 (1980).

12. Hamburger, V. & Yip, J. W. Reduction of experimentally induced neuronal death in spinal ganglia of the chick embryo by nerve growth factor. *J. Neurosci.* **4**, 767–774 (1984).

13. Deppmann, C. D. et al. A model for neuronal competition during development. *Science* **320**, 369–373 (2008).

14. Miguel-Aliaga, I. & Thor, S. Programmed cell death in the nervous system-a programmed cell fate? *Curr. Opin. Neurobiol.* **19**, 127–133 (2009).

15. Shirasaki, R. & Pfaff, S. L. Transcriptional codes and the control of neuronal identity. *Annu. Rev. Neurosci.* **25**, 251–281 (2002).

16. De Marco Garcia, N. V. & Jessell, T. M. Early motor neuron pool identity and muscle nerve trajectory defined by postmitotic restrictions in Nkx6.1 activity. *Neuron* **57**, 217–231 (2008).

17. Lallemand, F. et al. Positional differences of axon growth rates between sensory neurons encoded by runx3. *EMBO J.* **31**, 3718–3729 (2012).

18. Lallemand, F. & Ernfors, P. Molecular interactions underlying the specification of sensory neurons. *Trends Neurosci.* **35**, 373–381 (2012).

19. Kramer, I. et al. A role for Runx transcription factor signaling in dorsal root ganglion sensory neuron diversification. *Neuron* **49**, 379–393 (2006).

20. de Nooij, J. C., Doobar, S. & Jessell, T. M. Etv1 inactivation reveals proprioceptor subclasses that reflect the level of NT3 expression in muscle targets. *Neuron* **77**, 1055–1068 (2013).

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

22. Guenther, C. J., Miyamichi, K., Yang, H. H., Heller, H. C. & Luo, L. Permanent genetic access to transiently active neurons via TRAP: targeted recombination in active populations. *Neuron* **78**, 773–784 (2013).

23. Farinas, I., Yoshida, C. K., Backus, C. & Reichardt, L. F. Lack of neurotrophin-3 results in death of spinal sensory neurons and premature differentiation of their precursors. *Neuron* **17**, 1065–1078 (1996).

24. Ernfors, P., Lee, K. F., Kucera, J. & Jaenisch, R. Lack of neurotrophin-3 leads to deficiencies in the peripheral nervous system and loss of limb proprioceptive afferents. *Cell* **77**, 503–512 (1994).

25. Klein, R. et al. Disruption of the neurotrophin-3 receptor gene trkC eliminates the muscle afferents and results in abnormal movements. *Nature* **368**, 249–251 (1994).

26. Liebl, D. J., Tessarollo, L., Palko, M. E. & Parada, L. F. Absence of sensory neurons before target innervation in brain-derived neurotrophic factor-, neurotrophin 3-, and TrkC-deficient embryonic mice. *J. Neurosci.* **17**, 9113–9121 (1997).

27. Lawson, S. N. & Biscoe, T. J. Development of mouse dorsal root ganglia: an autoradiographic and quantitative study. *J. Neurocytol.* **8**, 265–274 (1979).

28. Patel, T. D. et al. Peripheral NT3 signaling is required for ETS protein expression and central patterning of proprioceptive sensory afferents. *Neuron* **38**, 403–416 (2003).

29. Inoue, K. et al. The transcription factor Runx3 represses the neurotrophin receptor TrkB during lineage commitment of dorsal root ganglion neurons. *J. Biol. Chem.* **282**, 24175–24184 (2007).

30. Levanon, D. et al. The Runx3 transcription factor regulates development and survival of TrkC dorsal root ganglia neurons. *EMBO J.* **21**, 3454–3463 (2002).

31. ElShamy, W. M., Fridvall, L. K. & Ernfors, P. Growth arrest failure, G1 restriction point override, and S phase death of sensory precursor cells in the absence of neurotrophin-3. *Neuron* **21**, 1003–1015 (1998).

32. Lefort, F., Clary, D. O., Rusoff, A. C. & Reichardt, L. F. Inhibition of the NT-3 receptor TrkC, early in chick embryogenesis, results in severe reductions in multiple neuronal subpopulations in the dorsal root ganglia. *J. Neurosci.* **16**, 3704–3713 (1996).

33. Carr, V. M. & Simpson, S. B. Jr Proliferative and degenerative events in the early development of chick dorsal root ganglia. I. Normal development. *J. Comp. Neurol.* **182**, 727–739 (1978).

34. White, F. A., Keller-Peck, C. R., Knudson, C. M., Korsmeyer, S. J. & Snider, W. D. Widespread elimination of naturally occurring neuronal death in Bax-deficient mice. *J. Neurosci.* **18**, 1428–1439 (1998).

35. Klinger, M., Chmura, S. A. & Killeen, N. Reporter alleles that inform on differences in Cre recombinase expression. *J. Immunol.* **184**, 6170–6176 (2010).

36. Haupt, C., Kloos, K., Faus-Kessler, T. & Huber, A. B. Semaphorin 3A-Neuropilin-1 signaling regulates peripheral axon fescapes and pathfinding but not developmental cell death patterns. *Eur. J. Neurosci.* **31**, 1164–1172 (2010).

37. Huettl, R. E., Soellner, H., Bianchi, E., Novitch, B. G. & Huber, A. B. Npn-1 contributes to axon-axon interactions that differentially control sensory and motor innervation of the limb. *PLoS Biol.* **9**, e1001020 (2011).

38. Hadjab, S. et al. A local source of FGF initiates development of the unmyelinated lineage of sensory neurons. *J. Neurosci.* **33**, 17656–17666 (2013).

39. Wright, D. E., Zhou, L., Kucera, J. & Snider, W. D. Introduction of a neurotrophin-3 transgene into muscle selectively rescues proprioceptive neurons in mice lacking endogenous neurotrophin-3. *Neuron* **19**, 503–517 (1997).

40. Lechner, S. G., Frenzel, H., Wang, R. & Lewin, G. R. Developmental waves of mechanosensitivity acquisition in sensory neuron subtypes during embryonic development. *EMBO J.* **28**, 1479–1491 (2009).

41. Woo, S. H. et al. Piezo2 is the principal mechanotransduction channel for proprioception. *Nat. Neurosci.* **18**, 1756–1762 (2015).

42. Ranade, S. S. et al. Piezo2 is the major transducer of mechanical forces for touch sensation in mice. *Nature* **516**, 121–125 (2014).

43. La Manno, G. et al. RNA velocity of single cells. *Nature* **560**, 494–498 (2018).

44. Niederreither, K., McCaffery, P., Drager, U. C., Chambon, P. & Dolle, P. Restricted expression and retinoic acid-induced downregulation of the retinaldehyde dehydrogenase type 2 (RALDH-2) gene during mouse development. *Mech. Dev.* **62**, 67–78 (1997).

45. Swindell, E. C. et al. Complementary domains of retinoic acid production and degradation in the early chick embryo. *Dev. Biol.* **216**, 282–296 (1999).
46. Berggren, K., McCaffery, P., Drager, U. & Forehand, C. J. Differential distribution of retinoic acid synthesis in the chicken embryo as determined by immunolocalization of the retinoic acid synthetic enzyme, RALDH-2. *Dev. Biol.* **210**, 288–304 (1999).

47. Niederreither, K., Vermot, J., Fraulob, V., Chambon, P. & Dolle, P. Retinaldehyde dehydrogenase 2 (RALDH2)- independent patterns of retinoic acid synthesis in the mouse embryo. *Proc. Natl Acad. Sci. USA* **99**, 16111–16116 (2002).

48. Kumar, S., Dolle, P., Ghyselinck, N. B. & Duester, G. Endogenous retinoic acid signaling is required for maintenance and regeneration of cornea. *Exp. Eye Res.* **154**, 190–195 (2017).

49. Ellis, H. M. & Horvitz, H. R. Genetic control of programmed cell death in the nematode C. elegans. *Cell* **44**, 817–829 (1986).

50. Marmigere, F. & Ernfors, P. Specification and connectivity of neuronal subtypes in the sensory lineage. *Nat. Rev. Neurosci.* **8**, 114–127 (2007).

51. Tourtellotte, W. G., Keller-Peck, C., Milbrandt, J. & Kucera, J. The transcription factor Egr3 modulates sensory axon-myotube interactions during muscle spindle morphogenesis. *Dev. Biol.* **232**, 388–399 (2001).

52. Ozaki, S. & Snider, W. D. Initial trajectories of sensory axons toward laminar targets in the developing mouse spinal cord. *J. Comp. Neurol.* **380**, 215–229 (1997).

53. Bentzinger, C. F., Wang, Y. X. & Rudnicki, M. A. Building muscle: molecular regulation of myogenesis. *Cold Spring Harb. Perspect. Biol.* **4**, a008342 (2012).

54. Appel, B. & Eisen, J. S. Retinoids run rampant: multiple roles during spinal cord and motor neuron development. *Neuron* **40**, 461–464 (2003).

55. Novitch, B. G., Wichterle, H., Jessell, T. M. & Sockanathan, S. A requirement for retinoic acid-mediated transcriptional activation in ventral neural patterning and motor neuron specification. *Neuron* **40**, 81–95 (2003).

56. Reijntjes, S., Francis-West, P. & Mankoo, B. S. Retinoic acid is both necessary for and inhibits myogenic commitment and differentiation in the chick limb. *Int. J. Dev. Biol.* **54**, 125–134 (2010).

57. Ji, S. J. et al. Mesodermal and neuronal retinoids regulate the induction and maintenance of limb innervating spinal motor neurons. *Dev. Biol.* **297**, 249–261 (2006).

58. Kikuchi, T. et al. Human iPS cell-derived dopaminergic neurons function in a primate Parkinson's disease model. *Nature* **548**, 592–596 (2017).

59. Steinbeck, J. A. & Studer, L. Moving stem cells to the clinic: potential and limitations for brain repair. *Neuron* **86**, 187–206 (2015).

60. Kirkeby, A. et al. Predictive markers guide differentiation to improve graft outcome in clinical translation of hESC-based therapy for Parkinson's disease. *Cell Stem Cell* **20**, 135–148 (2017).

61. Southwell, D. G. et al. Intrinsically determined cell death of developing cortical interneurons. *Nature* **491**, 109–113 (2012).

62. Bai, L. et al. Genetic identification of an expansive mechanoreceptor sensitive to skin stroking. *Cell* **163**, 1783–1795 (2015).

63. Hayashi, S. & McMahon, A. P. Efficient recombination in diverse tissues by a tamoxifen-inducible form of Cre: a tool for temporally regulated gene activation/inactivation in the mouse. *Dev. Biol.* **244**, 305–318 (2002).

64. Niederreither, K., Subbarayan, V., Dolle, P. & Chambon, P. Embryonic retinoic acid synthesis is essential for early mouse post-implantation development. *Nat. Genet.* **21**, 444–448 (1999).

65. Vermot, J. et al. Conditional (loxP-flanked) allele for the gene encoding the retinoic acid-synthesizing enzyme retinaldehyde dehydrogenase 2 (RALDH2). *Genesis* **44**, 155–158 (2006).

66. Lentz, S. I., Knudson, C. M., Korsmeyer, S. J. & Snider, W. D. Neurotrophins support the development of diverse sensory axon morphologies. *J. Neurosci.* **19**, 1038–1048 (1999).

67. Friedel, R. H., Schnurch, H., Stubbusch, J. & Barde, Y. A. Identification of genes differentially expressed by nerve growth factor- and neurotrophin-3-dependent sensory neurons. *Proc. Natl Acad. Sci. USA* **94**, 12670–12675 (1997).

68. van der Walt, S. et al. scikit-image: image processing in Python. *PeerJ* **2**, e453 (2014).

69. Picelli, S. et al. Full-length RNA-seq from single cells using Smart-seq2. *Nat. Protoc.* **9**, 171–181 (2014).

70. La Manno, G. et al. Molecular diversity of midbrain development in mouse, human, and stem cells. *Cell* **167**, 566–580.e519 (2016).

## Acknowledgements

We thank D. Ginty for the *TrkCCreER* mice; Y. Groner and D. Levanon for the *Runx3* mice; C. Ibanez for the *R26tdTOM* mice; P. Ernfors, H. Abdo and L. Calvo Henrique for the *R26CreERT2;R26tdTOM* mice; T. Jessell for the RUNX3 antibody; N. Ghyselinck and P. Dollé (IGBMC, Strasbourg) for the *Raldh* floxed mice; the CLICK imaging Facility supported by the Knut and Alice Wallenberg Foundation. We are also grateful to Prof. A. El Manira and P. Ernfors for critical reading of the manuscript. This work 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), the Ming Wai Lau Foundation and Åke Wiberg Foundation (F.L.); by RARENET V No. 1.7 EU Regional Development Fund, University of Strasbourg Institute of Advanced Studies (USIAS) (K.N.); and by ERC (BRAINCELL grant 261063), VR, the Wellcome Trust (grant 108726/Z/15/Z), and the EU (FP7/DDPDGENES) (S.L.). F.L. is a Ragnar Söderberg fellow in Medicine, a Wallenberg Academy Fellow in Medicine and a MWLC investigator. Open access funding provided by Karolinska Institute.

## Author contributions

S.H. and F.L. designed and supervised the experiments and either has the right to list her/himself last in bibliographic documents. Y.W., H.W., G.A. and S.H. carried out the immunostaining experiments and microscopy analysis. Y.W and S.H. performed the in vitro experiments. S.H. and H.W. performed the genetic tracing experiments. H.W. and G.L.M. performed the bioinformatics analyses, helped with related figures and wrote associated methods section. S.H. and G.L.M. conceived and supervised the RNA velocity and fate bias analyses. Y.W., A. Sh., C.P., Y.X.-F. and D.P. performed the gene expression experiments and provided help with genetic lines breeding. Y.W., H.W., G.L.M., S.H. and F.L. performed data analysis. S.C. and S.L. carried out and analyzed smFISH experiments. N.A. and I.A. performed chicken electroporation, analyzed by P.F. and Y.W. K.N., F.D.S., G.C., A. Sc. and A.M. processed and provided mouse tissues for analysis. S.H. prepared the figures. S.H. and F.L. wrote the paper with input from all authors.

## Additional information

Supplementary Information accompanies this paper at https://doi.org/10.1038/s41467-019-12119-3.

**Competing interests:** The authors declare no competing interests.

**Reprints and permission** information is available online at http://npg.nature.com/reprintsandpermissions/

**Peer review information** *Nature Communications* thanks Bruce Carter and the other, 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.

<figure><img src="image_23.png" /></figure>

**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/.

© The Author(s) 2019