Single-cell RNA-sequencing analysis of the developing mouse inner ear identifies molecular logic of auditory neuron diversification

Different types of spiral ganglion neurons (SGNs) are essential for auditory perception by transmitting complex auditory information from hair cells (HCs) to the brain. Here, we use deep, single cell transcriptomics to study the molecular mechanisms that govern their identity and organization in mice. We identify a core set of temporally patterned genes and gene regulatory networks that may contribute to the diversification of SGNs through sequential binary decisions and demonstrate a role for NEUROD1 in driving specification of a Ic-SGN phenotype. We also find that each trajectory of the decision tree is defined by initial co-expression of alternative subtype molecular controls followed by gradual shifts toward cell fate resolution. Finally, analysis of both developing SGN and HC types reveals cell-cell signaling potentially playing a role in the differentiation of SGNs. Our results indicate that SGN identities are drafted prior to birth and reveal molecular principles that shape their differentiation and will facilitate studies of their development, physiology, and dysfunction.

T he ability to detect and discriminate auditory stimuli depends on sensory coding of sound components realized in the cochlea by the spiral ganglion neurons (SGNs). These primary auditory neurons receive synaptic input from the auditory hair cells (HCs, the sensory receptors located in the organ of Corti) at their distal dendrites and convey all auditory information to the cochlear nuclei in the brainstem through their central projections that form the auditory nerve. Recent molecular classifications have identified at least four different types of SGNs in mice (I a -, I b -and I c -SGNs and II-SGNs) 1-3 , consistent with the physiological and anatomical diversity of primary auditory afferents 4,5 : I-SGN fibers exhibit large differences in their spontaneous discharge rates and sensitivities, which correlate with the location and structural features of their synaptic contact at the base of the IHC. Such diversity of I-SGN fibers makes them capable of collectively span all levels of sound intensity. Therefore, abnormal differentiation and connectivity of the four SGN types is expected to disrupt their functional organization and thereby critically impact auditory function. However, when and how SGN subtype identities are established and the transcriptional mechanisms by which they emerge are currently unknown.
SGNs are produced in a base to apex progression from neuroblasts that have delaminated from the otocyst 6 and migrated to form a dense ganglion along the medial side of the inner ear sensory epithelium between E10 and E12 in mice 7 . As the cochlear duct elongates, these early SGNs extend peripheral projections through the expanding otic mesenchyme (OM) population to reach the developing sensory epithelium around E15-17 8,9 . Thereafter, SGN projections continue to grow within the cochlear epithelium to form synaptic connections with the HCs 10 . Simultaneously, the central projections of SGNs reach the cochlear nuclei 11 , making contact with neurons of the central auditory pathway. Slightly before their innervation by SGNs, HCs differentiate into inner HCs (IHCs) and outer HCs (OHCs) 12,13 and complex nerve endings below OHCs suggesting differentiation of type II SGNs (II-SGNs)-like dendrites begin to be observed during late embryogenesis 11 . At birth, the characteristic morphologies of type I SGN (I-SGN) and II-SGN afferents below the HCs are clearly apparent. Accordingly, multiple I-SGNs innervate radially one IHC-which in the adult transduces the physical energy of sound into electrochemical signals-while, by contrast, single II-SGNs form en passant synaptic connections with a dozen of OHCs that in the mature cochlea modulate the sensitivity and therefore indirectly the output of the IHCs. Moreover, already by this stage, SGNs exhibit a large diversity of molecular profiles, varying connectivity patterns and distinct intrinsic physiological properties which to some degree prefigure the functional organization evident in the adult 1,14 . An unresolved question is therefore to understand the timing and diversity of the molecular events that produce and assemble sensorineuronal specificity in the developing cochlea and eventually enable animals to detect and respond to a plethora of auditory stimuli.
In this study, we used single-cell RNA-sequencing (scRNAseq) to define the molecular logic by which SGN subtypes diversify in the mouse. Our detailed analysis of transcriptional dynamics of SGNs from several embryonic and perinatal stages indicates that neuronal subtypes successively emerge during HC afferent innervation and provides a comprehensive collection of molecular states and candidate gene regulatory networks associated with each lineage and underlying fate divergence. Importantly, we uncover the functional consequence of Neurod1 deletion on the first binary decision between the presumptive I c -SGNs and the rest of the neuronal lineage. Finally, we catalogue chemotropic signaling that would define the complex cochlear wiring and identify deafness genes associated with distinct cell states of developing SGNs and HCs.

Results
Embryonic emergence of SGN diversity. To analyse the molecular changes associated with the diversification of SGN and HC types during development, we used flow cytometry to sort tdTomato + (TOM + ) cells isolated from E14.5, E15.5, E16.5, E17.5 and E18.5 (Fig. 1a) and from P3 (previously published) 1 cochlea of Ntrk3 Cre ;R26 tdTOM or PV Cre ;R26 tdTOM mice and sequenced their mRNA with high coverage using the Smart-seq2 protocol 15 . The selected stages cover the earliest time point (P3) SGN subtypes have been defined transcriptionally 1 as well as a key period in the embryonic development of the neurosensory elements of the cochlea characterized by the innervation of HCs by SGNs 11,16 .
A total of 2308 cells were pre-processed and clustered with the pagoda2 pipeline. We applied Harmony 17 , an algorithm to bridge time points, and combined it with Palantir 18 , to generate a diffusion space which was then used as a basis for Force Atlas 2 (FA) embedding and for subsequent trajectory analysis (see "Methods"). The dynamic of the gene expression on our multidimensional integrated dataset identified 19 clusters (Cl., Fig. 1a), including intermediate states with varying degrees of cell fate biases/lineage restrictions of the SGNs, HCs and OM compartments (Fig. 1b). The neuronal compartment was characterized by the expression of general neuronal markers such as Tubb3 (βIII-tubulin) and Elavl3 (HuC) (Fig. 1b), as well as Actl6b, Map2, Gap43 or Uchl1 (PGP9.5) (Supplementary  Fig. 1) 19,20 . Finally, Cl.19 expressed Pou3f4 and Tbx18, which are known to be enriched in OM cells (Supplementary Data 2) 9, 21 . Tbx1 and Car3 have also been shown for the OM and are expressed in our dataset 22,23 , together with additional gene markers for the OM population, amongst which Prrx1 and Twist1 (Fig. 1b, 1145 genes found differentially expressed in OM cluster compared to the E16.5 SGN clusters, Supplementary Fig. 1). This suggests that Cl.19 is likely an OM contaminant population from the FACS purification. This cluster was removed from the dataset for downstream analysis.
Focusing on the neuronal compartment (Fig. 1c, d), alignment of the time points revealed a neuronal differentiation progressing from Cl.1 to Cl.10, 13, 15 and 17 which at P3 represent specified SGN types as identified by the expression of previously published cell-type specific genes for II-, I a -, I b -and I c -SGNs 1 (Fig. 1f). The top differentially expressed genes for each cluster is provided in Fig. 1e. Of note, the I b and I c trajectories almost joined at late time-points due to their high similarity in gene expression. However, overlaying them onto FA embedding revealed two paths leading to these two branches, one (the I c ) directly from the unspecialized population (see also Fig. 2a) and the other (I b ) going through several intermediate states and also giving rise to I a and II populations (see "Methods" section). Moreover, along the differentiation trajectories, consecutive developmental time-points partially overlapped (Fig. 1d). This is likely explained by the basal-to-apical gradient of neuronal differentiation in the developing cochlea and the fact that whole spiral ganglia (from base to apex) were sampled. Therefore, a continuum of neuronal differentiation states and branching can be observed in a single developmental time point.   The cell state/type specific marker expression was validated in situ using RNAscope and immunostaining (Fig. 1g), confirming at E16.5 the emergence of a I c identity (I c :Mgat4c + and Lypd1 + ; I a / I b /II: Islr2 + ) and at E18.5 or P0 (Fig. 1h) the existence of the three main subtypes of I-SGNs [I a : CR + (CR for calretinin); I b : Lypd1 + / CR + and Calb1 + ; I c : Lypd1 + /CR − ; at this stage, Runx1 marks both I b and I c and Igfbpl1, the I a , I b and II] and the II-SGNs [Plk5 + and PRPH + (peripherin)]. Overall, this progression recapitulates with the selected experimental time points the dynamic of transcriptional changes occurring in SGNs from E14.5 to P3.
Molecular codes of neuronal diversification. To study the flux of molecular changes that occur within the neuronal differentiation continuum, we fitted a principal tree in diffusion space, excluding the HC and OM clusters. The tip of the branch that was enriched with E14.5 cells was selected as a root and pseudotime was subsequently calculated as the distance on the tree from that root. The tree was then represented into a dendrogram recapitulating a branched trajectory based on the transcriptional similarity of pseudotime-ordered cells (Fig. 2a). The resulting tree accurately reflected developmental stages and the branching features of the FA representation (Fig. 2b). Plotting gene expression of known neuronal maturation genes, including Slc17a7 (VGLUT1), Nsf, Syp (synaptophysin), Stxbp1 (Munc18-1), Snap25, Grin1 (NMDAR1) and Nefh (neurofilament heavy chain) on the tree supported the developmental progression from an immature, unspecialized state to differentiated neuron subtypes (Supplementary Fig. 2d). Other neuronal maturation genes showed instead a cell-type restricted expression pattern (Supplementary Data 1). Interestingly, some genes such as Cplx2 (complexin 2) and Stx1a (syntaxin 1a), which are members of the SNARE family and associated with the synaptic release of neurotransmitters, exhibited a reversed expression pattern with enrichment in immature neurons, during extension of SGN afferents to their targets 24 (Supplementary Fig. 2d), which suggests a role in the expansion of the plasma membrane during axonal growth 25 towards the sensory epithelium. Therefore, the neuronal differentiation tree could be divided into 2 major states, unspecialized versus differentiated neurons (the branching tree leading to the known terminal cell types). Notably, the tree showed unspecialized neurons diverging into I c -SGNs and intermediate I a /I b /II-SGNs around E15-16, which marks the period of afferent innervation of HCs 26 . Knowing the first SGNs are generated at E10 in mice, this indicates that neuronal diversification in the cochlea is a relatively late event that might require cell-cell communication with their peripheral target field, but also, that it is initiated prior to E18.5 in mice, i.e., before functional synapses first emerge and correlated activation of SGNs is observed in the cochlea [27][28][29] .
To investigate the dynamic transcriptional changes along the individual segments of the trajectories (Fig. 2a, c-g, dendrograms) and discover transcription and signalling genes with likely functions in SGN differentiation, we used differential gene expression approaches that characterize pseudotime-ordered molecular trajectories ( Fig. 2c-i, Supplementary Data 3-6). Focusing on transcriptional regulators (Fig. 2i, j), we observed that the shared, unspecialized neurons' trajectory was characterized by rapid downregulation of genes commonly associated with early neuronal differentiation processes (Neurod1, Nhlh1 or Ebf2), and upregulation of genes that are involved in specification events in cochlea and other systems (Zfhx3, Runx1, Mafb, Meis2, Id1 and Myt1l) (Fig. 2i, j). Those genes were later found downregulated, restricted to or over-represented in select neuron type trajectories. Runx1 and Meis2 for instance were found further increased in the I c -SGNs as they differentiate from the unspecialized state. In contrast to Runx1 however, Meis2 was not maintained in I c -SGNs, and Runx1 was later found increased also in I b -SGNs (see below) ( Fig. 2i, j). Pou4f1 (BRN3A) showed a similar trend, i.e., maintained or slightly increased in the I c -SGN trajectory, while dropping in the intermediate I a /I b /II population trajectory, which is in line with a recent study 30 . At the opposite, Gata3, which increased in the unspecialized neurons, as previously shown 31 , was downregulated in all type I SGNs when each population diverged from their shared trajectory (which happens earlier for I c -SGNs) (Fig. 2i, j). This confirms previous observations 31 and suggests that its master regulator function in specifying a generic identity of SGNs 32 might be incompatible with the final differentiation program of type I subtypes. Gata3 was however maintained in II-SGNs, together with Mafb (which in parallel progressively decreases in I-SGNs) (Fig. 2i, j), confirming earlier studies on their expression 1-3,33,34 and known molecular interactions-with MafB acting downstream of GATA3 to regulate auditory synaptogenesis 35 .
Following the first branching point (or split), the intermediate I a /I b /II trajectory was primarily characterized by a general maintenance of specific transcription factors (TFs) expression found before bifurcation and the upregulation of Id1 and of Gfra1 (GDNF family receptor alpha 1) ( Fig. 2h-j). This suggests that this intermediate state is marked by a progression of neuronal specification events that had already started during the last period of the unspecialized neuron trajectory. In contrast, the further distinction between type II and the transient I a /I b trajectories was marked in the nascent type II neurons by the absence/decrease of expression of TFs specific to the transient I a /I b population such as Runx1 or Shox2 and the upregulation of Sox9 and Tshz3 (Fig. 2i,  j). Also, and similar to the I c population, Rora (retinoid-related orphan receptor alpha) and Prox1, a generic marker of adult type I SGNs 1 , were both upregulated in transient I a /I b SGNs (Fig. 2i, j). This suggests that Prox1 could be part of a terminal differentiation program 36 that would be continuously required to maintain a type I specific differentiated state. Finally, the last branching event identified I a -and I b -SGN differentiation coinciding in the I b trajectory with an upregulation of Runx1 and Pou4f1 (although Pou4f1 showed lower expression relative to I c neurons), and in I a -SGNs, with a downregulation of Runx1 and the maintenance, although at lower levels, of Id1 (Fig. 2i, j). Overall, the data provide a compendium of genes associated with identity divergence, and describe a series of transient, up-and downregulation of TFs whose expression dynamics may play a significant role in neuronal diversification and cell identity maintenance (see also Fig. 4j).
Regulon analysis identifies NEUROD1 as a I c fate regulator in vivo. Cell type trajectories are shaped by underlying gene regulatory networks (GRNs) that are centred on a limited number of TFs (or master regulators) and co-factors that interact with cisregulatory genomic regions (target genes) to mediate a specialized transcriptional programme that governs individual cell type/state definition. To comprehensively reconstruct GRNs along the neuronal differentiation tree, and reveal the master regulators and co-factors (i.e., regulons) that might govern individual cell type identities, we used SCENIC (Single-Cell rEgulatory Network Inference and Clustering) 37 , a computational workflow that enables inference of GRN (regulon) activities. We identified multiple regulons, each representing a TF, along with a set of coexpressed and motif-enriched target genes, and the regulon activity scores for each neuronal trajectory (Fig. 3a). While some of these regulons were shared among multiple trajectories, others were highly specific and non-overlapping and defined the developmental progression of select cell identities. Moreover, many regulons were found either up-or downregulated, or transiently expressed (e.g., Pou2f1(+), Onecut2(+) or Rora(+)), along the developmental trajectories, highlighting the need of a decrease in specific GRNs activity for driving proper SGN differentiation. Interestingly, the upregulation of Ppargc1a(+) (PGC-1α gene network), that is associated with mitochondrial biogenesis 38 , in all emerging type I neurons-with a large energy demand later in life 39 -suggests their possible priming with an increased metabolism already before birth. This co-existence and temporal segregation of various GRNs in each trajectory produce thus combinatorial codes for SGN fate decisions during development.  Our analysis also revealed that the highest representation of regulons defined the unspecialized group of neurons. Because the enrichment and number of effector genes, together with the expression level of the master TF, are key parameters in predicting the activity of a regulon, a potential limitation of the above analysis is the possibility that a regulon might significantly change activity, and thus visibility, depending on the cellular context (co-factors and gene modules involved) such as between two temporally and molecularly separated events, e.g., during neurogenesis (immature state) or neuronal diversification. By reanalysing the activity of neurogenesis-related regulons along the pseudotime axis, we observed that Neurod1(+), active in all cells during early specification of the unspecialized SGNs, was also and specifically active at the beginning of the I c trajectory (Fig. 3b, c). This temporal analysis of Neurod1(+) identified two different regimes of Neurod1-asssociated GRNs. A first one was associated with the early stage of neuronal differentiation, included targets such as Nhlh1, Neurod6 and Ebf2, and decreased along the first, unspecialized neurons, trajectory (Fig. 3d). A second one progressively increased along this trajectory, diverged to the I c path, and was marked by the expression of targets including Runx1, Pou4f1 and Lypd1, which all characterize the emergence of a I c identity. To determine if NEUROD1 could influence I c -SGN differentiation, we crossed Neurod1 loxP/loxP mice 40 with Isl1 Cre mice 41 (ISL1 is expressed from E8.5 in neurons of the developing inner ear 42 ) to delete Neurod1 in SGNs (Neurod1 cKO ), therefore avoiding the nearly complete loss of SGNs observed in the full Neurod1 knockout mice 43,44 (Fig. 3e). Neurod1 cKO mice have previously been shown to exhibit normal organization of the organ of Corti and to enable development of SGNs at basal and mid-basal regions, albeit in reduced number 45 . Analysis of cochlea in P3 Neurod1 cKO mice revealed a loss of I c -SGN marker staining (Fig. 3f, g). This phenotype was similar at P0 and already visible at E16.5 ( Fig. 3h-j). To determine if NEUROD1 is required for I c -SGN specification, or if it regulates the survival of this emerging population, or both, we examined cleaved caspase-3 a progression of regulons over the tree  (c-casp-3) immunostaining at E16.5 and P0 in Ctrl and Neurod1 cKO mice cochlea. Analysis of the basal and mid-basal regions of the cochlea revealed in average about 1 neuron positive for c-casp-3 labelling every two sections in Neurod1 cKO mice, with no labelling in Ctrl animals ( Supplementary Fig. 3). Therefore, while we cannot rule out a potential role of NEUROD1 in early differentiating I c -SGN cell survival, we believe the small presence of SGN cell death observed during the period of neuronal diversification in the absence of NEUROD1 cannot account for the great loss of I c -SGNs seen at E16.5 and around birth. Taken together, these results suggest that NEUROD1 is required for the early specification of I c -SGNs from the unspecialized pool of SGNs.
Molecular regulation of trajectories and cell states. Emergence of cell types during development is a continuous process of differentiation where more immature cells progressively become fate restricted in a series of stepwise bifurcation events. These branch points however only indicate an overall change or switch of transcriptional identity along the trajectory, and do not explain the critical events that characterize the pre-and post-transition states and are responsible for the emergence and consolidation of cell fate 46 . To analyse the dynamical behaviour of the hierarchical fate split points that represent diversification (decision) events, we analysed and identified for each branching points gene modules (groups of genes that change in the same direction and tend to synchronize along the pseudotime) possibly driving fate choice (early modules, pre-bifurcation) and fate biasing (late modules, post-bifurcation) before actual fate commitment and during consolidation, respectively. Each bifurcation event was preceded by a period of increasing heterogeneity in cell type-specific module expression, suggesting higher transcriptional differences between the alternative cell fates in each single cell while approaching the bifurcation point which correlated with preference towards a specific fate choice (Fig. 4a- Fig. 4j. Ligand-receptor interactions between HC and SGN types. To reveal potential communication processes between developing SGNs and their neighbouring cells, we first manually curated genes coding for receptors of well-known pathways that showed changing expression along the pseudotime in our dataset. We noted that, with a few exceptions, most changes in genes associated with morphogens (WNT, BMP, RA and SHH)-, growth factors-and hormones-related signalling were mostly specific to either all type I or type II SGNs (Fig. 5a and Supplementary  Figs. 4 and 5). Also, the specific increase of Nbl1, Smurf2 and Smad6/7/9 in II-SGNs suggests a selective inhibition of BMP signalling in this cell type, as previously suggested 1 . Interestingly, not only receptors and modulators or intracellular signalling components were differentially expressed in SGNs, but also ligands, with for instance Wnt5a and Shh, which were detected in the type II lineage and the unspecialized population (confirming an earlier study 52 ), respectively, suggesting instructive function from the neurons themselves, as previously demonstrated for SGN-derived SHH in HC differentiation 53 . Genes linked to different types of cell-cell adhesion and axon guidance molecules showed, in contrast, a strong cell type trajectory and/or temporal heterogeneity ( Fig. 5a and Supplementary Fig. 5). For instance, in the cadherins and protocadherins superfamily, Cdh23 was found differentially expressed among differentiating type I SGNs at late stages, with higher levels in I b neurons, Pcdh17 was initially expressed in immature SGNs and became restricted to the type II trajectory, and Cdh9 was gradually increased in the transient I a /I b and maintained in postnatal I a and I b neurons, though with higher levels in I a -SGNs. In axon guidance molecule families, while some genes showed high specificity to specific cell types (e.g., Epha3 in I b /I c /II trajectories or Epha6 and Plxna2 in II-SGNs), most genes were broadly expressed but with various levels of expression in distinct trajectories ( Fig. 5 and Supplementary Fig. 5). While transcripts at low levels might eventually not be translated or could result in protein levels too low to be functional, these results may also indicate that the construction of the peripheral auditory circuits is at least partly defined by a combinatorial and temporal code of levels of genes, as previously suggested 26 . In this study, the authors showed that, before birth in mice, Semaphorin-3F (SEMA3F, expressed in the OHC region) acts as a repulsive signal for I-SGN processes via the activation of neuropilin-2 (NRP2), and suggested that possible differential expression of NRP2 and/or co-factors (e.g., plexins) might be accounted for the lack of responsiveness of II-SGNs to SEMA3F 26 . In our data, Nrp2 is progressively upregulated in all I-SGNs from E16.5, but less in II-SGNs ( Fig. 5 and Supplementary Fig. 5). Moreover, Nrp1 levels are at the same time increased in II-SGN trajectory, as well as that of Plxna2 and Plxnc1, suggesting other possible semaphorin signalling, together with select ephrin signalling, acting specifically on II-SGN afferents.
While the distinct innervation patterns of I a -, I b -and I c -SGNs might be temporally regulated, with a gradient of development from the modiolus to the pillar side of the IHCs, a gradient of signalling molecules was also observed, for instance in the axon guidance-related genes, e.g., Epha3, Ntng1/2, Plxna4, Sema3a/4d, and Slit1/2, but also in the cell adhesion families of molecules, e.g., including cadherins and protocadherins ( Fig. 5a and Supplementary Figs. 4 and 5). Overall, our data thus suggest that during their innervation of the sensory epithelium and their central target neurons, the different classes of SGNs might assemble various types of receptors and adhesion molecules in time and space to mediate the selective contributions of diverse attractants and repellents from cells of their surrounding environment (e.g., epithelial cells, HCs or neurons).
We next explored the potential receptor-ligand pairing possibilities between the different HC and SGN types during the early stages of HC innervation. To differentiate between IHC and OHC in our dataset, we re-clustered the 39 cochlea HCs of E18.5 from the HC cluster (Cl.18 in Fig. 1a) and obtained two clusters whose gene expression was consistent with the known molecular identity of immature IHCs (Fgf8, Rprm, and Trh) and OHCs (Bcl11b, Scn11a, and Insm1) 12,13 (Fig. 5b-d and Supplementary Fig. 6). At this early stage of differentiation, a high number of differentially expressed genes (359) distinguished OHCs from IHCs ( Fig. 5e and Supplementary Data 8, 9). The relatively low representation of OHCs in our dataset, despite they represent 75% of cochlea HCs in vivo, might be explained by the delayed differentiation of OHCs compared to IHCs 54 coupled to the fact that PV (used to drive Cre expression in our tracing strategy) is a late HC differentiation marker 55 Fig. 4 Gene modules defining cell fate choice and commitment along the differentiation trajectories. a-i Analysis of the bifurcations of the differentiation tree representing cell fate selection between distinct neuronal cell lineages: I c versus I a /I b /II (a-c), II versus I a /I b (d-f) and I a versus I b (g-i). a, d, g Scatter plots show average expression of lineage-specific modules in each cell along the trajectory. Early competing modules show gradual co-activation, followed by selective upregulation of one fate-specific module and downregulation of the alternative fate-specific module. Late modules show almost mutually exclusive expression within the two lineages after bifurcation reflecting commitment to either fate. Colors encode tree branches as in Fig. 2a. b, e, h Representative genes in each module; asterisks highlight TFs (see also Supplementary Data 7). c, f, i Average local correlations of early gene modules with branch-specific correlations, in cells with similar developmental pseudotime (in black on the trajectories); the difference between intra-and inter-module correlations is shown in the upper right corner of the correlation plots and would reflect the repulsion between modules. j Summary scheme of the cellular diversification of developing SGNs via distinct trajectories. The transcriptional regulators predicted to be involved in the unfolding of the different lineages are represented on the different trajectories and are derived from the analysis of Figs. 2-4. Genes in colored boxes indicate TFs of early modules from the bifurcation analysis and potentially driving specific cell fate choice (arrows). The arrow symbols after the genes indicate an upregulation or downregulation within a trajectory, while they indicate a high or low expression, respectively, in the final states (I a , I b , I c and II states at the far-right end of the differentiation tree). Int., intermediate.

interrogated potential cell-cell interactions between HCs and
SGNs using the permutation-based tool CellPhoneDB 56 that computes a communication score and evaluates the significance of each known ligand-receptor pair (from public resources) found within the scRNAseq data. Importantly, this database considers the multisubunit nature of the protein complexes to identify likely functional ligand-receptor interactions. This method revealed potentially active ligand-receptor pairs from HCs to SGNs and vice versa ( Fig. 5f and Supplementary Fig. 7). However, although in other systems instructive signalling from neurons to their peripheral target cells are known to be crucial for the development of the targeted cells 57 , both IHC and OHC have been shown to develop in the absence of innervation 58,59 , we therefore here only focused on the outgoing HC to SGN signalling as it is more likely to represent an actual signalling network. Moreover, while OHCs are known to interact with I-SGNs, mostly through repulsive signals 26,60,61 , IHC to II-SGN negative interactions are more unlikely to occur at E18.5 60 and were then discarded from the analysis. Also, we filtered ligandreceptor modules based on their expression level and on the number of cells expressing the specific interactors as to provide higher confidence of their biological relevance. Eventually, only cell-type specific ligand-receptor pairs with the highest score are shown (Fig. 5f). More than half of the ligand-receptor pairs involved the ephrin family confirming their important role in shaping cochlear wiring [61][62][63] . Also, a population of transient I a /I b was still observed at E18.5 and showed almost systematically a distinct cell-signalling score with IHCs compared to either I aand/or I b -SGNs, illustrating a clear parallel between cell type differentiation and their changing pattern of cell-cell communication with HCs. When focusing on OHCs, an interesting observation is the higher score observed for the repellent  Fig. 4 and 5). The color bars at the top indicate cell states as in Fig. 2a. b tSNE plot of the IHC and OHC clusters at E18.5. c Representation of selected genes specific to IHC (Fgf8 and Trh) and OHC (Bcl11b and Scn11a) on the tSNE plot. d In vivo confirmation of Scn11a and Trh expression in HCs at E16.5 and E18.5 using immunohistochemistry and in situ hybridization (RNAscope). Note that Scn11a is a general marker for both HC types at E16.5. e Dot-plot of the top 25 DEGs between IHC and OHC at E18.5. f CellPhoneDB analysis of the top selected potential outgoing cell-to-cell signalling from HCs to SGNs (see also Supplementary Fig. 7). Data in a, c and f represent max normalised fitted gene expression, Knn smoothed log10(fpm) expression and log2 mean of (interacting molecule 1, interacting molecule 2), respectively. Genes referenced through the text are highlighted in red in a and f. Scale bar: 20 µm. OC organ of Corti.
SEMA3F-NRP2 signalling from OHCs to I-SGNs, confirming previous observations 26,64 , while OHCs to II-SGNs signalling was marked by the SEMA3A-NRP1 pair, which is known to be chemoattractant for neurites 65 .
While providing a series of potentially active signalling pathways between HCs and SGNs during a period of important morphological events that characterize the innervation of HCs, these only represent a snapshot of the many possible communications during SGN development and do not consider the possible signalling between SGNs and other cell types. However, because of the depth of our sequencing, these data will be highly valuable to functionally access the molecular basis of the development of the exquisite organization of HC innervation by SGNs.
Cell-type-specific deafness-associated genes. To confirm our dataset but also provide insights into the spatio-temporal expression profile of genes associated with hereditary hearing loss (HHL), we extracted HHL-related genes from the Deafness Variant Database (http://deafnessvariationdatabase.org/references) and performed systematic analysis of their expression profile in both developing HCs and SGNs from our databases (Fig. 6a, b and Supplementary Data 10,11). This analysis identified ninety-two known genes, amongst which about half were either uniquely expressed or enriched in HCs. This indicates that many genes previously associated with HCs function are also potentially expressed, at least transiently, in SGNs. For instance, in our dataset, Espn (DFNB36), necessary for late HC stereociliogenesis and postnatal stability/maintenance 66,67 , is specifically and progressively enriched in I c -SGNs as they differentiate. Similarly, Lhfpl5 (also known as TMHS, DFNB66/67), which is critical for HC mechanotransduction 68 , is progressively increased in all SGNs during their differentiation, and Myo6 (DFNA22/ DFNB37), necessary for HC structural integrity 69 , is also expressed in all SGNs at early stages.
Our dataset allowed us to also refine knowledge of expression of genes within SGN and HC types. Pdzd7 (USH2C), Pnpt1 (DFNB70) 70 , Prps1 (DFNX1) 71 , Tmprss3 (DFNB8/10) 72 and Wfs1 (DFNA6/14/38) 73 , which were known to be expressed in SGNs, showed enrichment in particular subclasses of neurons and/or stages of differentiation. On the other hand, genes known for their expression in other structures of the cochlea were also found in HCs and SGNs. For example, Tecta (DFNA8/12/ DFNB21) 74 , a major non-collagenous component of the tectorial membrane, was found expressed in IHCs at a surprisingly high level. Moreover, the connexin genes Gjb2 (Cx26, DFNA3/3A/B1/ 1A) and Gjb6 (Cx30, DFNA3/3B/1B) which are crucial for the functional differentiation of HCs 75 and yet known to be expressed in non-sensory cells of the organ of Corti postnatally, were also expressed in IHCs themselves at E18.5 (Fig. 6c), arguing for a potential earlier function of connexins in this cell type.
Together, this analysis uncovers spatio-temporal distributions of deafness genes within developing HCs and SGNs, with some having a transient and cell-type specific expression during embryogenesis. Although deletion or mutations of these genes may not lead to observed dysfunction in either SGNs or HCs in vivo, this updated pattern of expression might help in future studies on cellular mechanisms linked to HHL.

Discussion
The recent demonstration of diverse molecular types of SGNs has provided a solid basis for the distinct anatomical and physiological properties of auditory afferents [1][2][3]76 , yet how these neuronal identities are generated has remained unexplored. Using scRNAseq-based analysis of developing sensorineural cells of the cochlea, we demonstrate a continuous representation of changes in gene expression that define the course of SGN diversification. We further identify the molecular framework associated with their neuronal fate choices or implementation and cell map expression of genes associated with deafness in developing SGN and HC types. Moreover, our study identifies neuron type specific and temporal differences in expression of chemotropic signaling, notably from HCs, that could potentially act as strong determinants of neuronal differentiation. This work likely contributes very importantly with mechanistic understanding of SGN diversification and will allow researchers to uncover potentially important biological pathways that shape the functional architecture of the primary auditory afferents.
Our study reveals that the fates of SGN subtypes are defined before the first coordinated bursts of action potentials are observed in the auditory nerve, which in rodent occurs perinatally [27][28][29] . This suggests that SGN diversity does not emerge in response to activity-dependent mechanisms, but instead by specific transcriptional programs that unfold over the course of few days during late embryogenesis. Early spontaneous network activity seen in the postnatal cochlea undoubtedly participates in the maturation and plasticity of the ascending auditory pathways 3,77,78 , however most aspects that are directly linked to the functional diversity of SGNs seem to be intrinsically defined before birth. Moreover, unlike the main somatosensory neuron types which differentiate from the neural crest stem cells almost immediately after they exit cell cycle 46 , the differentiation of SGNs into subtypes resembles the second phase of neuronal diversification in the somatosensory system in which distinct extrinsic cues act on axons of sensory neurons for generating   further cellular diversity 47,[79][80][81][82] . It appears that in the inner ear, this transition through subsequent maturation steps, during which neuronal cells acquire different transcription factor networks, is a relatively long process. It is in this context that NEUROD1, which is essential with other co-factors for the development of a spiral ganglion neuron fate during neurogenesis 43,83 , was later found transiently expressed and necessary for the emergence of a I c phenotype through certainly the recruitment of a distinct gene regulatory network. In the developing cochlea, this timing of diversification coincides with the innervation of the sensory epithelium by the unspecialized SGNs 7 . IHCs and OHCs, which can already be transcriptionally distinguished at E14.5 in mice 13 , are important candidate cell types for releasing extrinsic signals essential for the differentiation of SGNs into subclasses. While this is well illustrated by the many potentially active cell-cell signaling cassettes between HCs and SGNs in this study, cell-cell communication with OM 9 , surrounding glial cells and non-sensory cells of the developing organ of Corti will need to be further studied by cross-comparison computational analysis with recent 13 and future single cell data. Indeed, although a possible spatial (lateral polarity, pillar versus modiolar sides), subcellular localization of signaling in the developing IHCs might play a crucial role in shaping neuronal identities and connectivity within the type I neurons, direct or indirect interactions with pillar cells (facing the I a -type) or cells of the Kölliker's organ (facing the I c -type) for instance might also contribute to SGN differentiation.
Another interesting aspect of this study is the developmental history of SGN type lineage. In addition to the comprehensive insights it provides into the molecular programs of SGN differentiation, it also reveals a basic temporal outline of neuronal diversification in the cochlea wherein a I c identity differentiates first from a common unspecialized pool of SGNs, followed by differentiation of a type II identity. The I b identity emerges later in the transcriptional hierarchical tree from a common I a/b lineage, which also leads in parallel to the I a fate. Hence, we propose that diversification into two main types of neurons that differ in their threshold of sensitivity (highversus low-threshold cochlear neurons) and innervation pattern in birds 84 have preceded the emergence of the type II neurons (and presumably the I b type) that appeared in mammals. Moreover, we suggest that the I a identity could represent a default path since no specific regulons have been identified to define their differentiation and the profile of expression of specific TFs such as GATA3 (linked to a general program of auditory neuron development) 32 and of peripherin, commonly associated with all immature neurons of the peripheral sensory system during embryogenesis 50,85 , were still found in I a -SGNs around birth. Therefore, the I a trajectory might represent a relatively plastic differentiation path that could have been amenable throughout evolution to the implementation of genetic programs leading to cell type diversification and innovation, including the II-and I b -SGNs.

Methods
Ethics and experimental animals. All animal care and procedures were performed in accordance with the national guidelines published by the Swedish Board of Agriculture and approved by the local ethics committee of Stockholm, Stockholms Norra djurförsöksetiska nämnd. Mice were housed in groups, with standardized pellet food and water ad libitum and under 12 h light-dark cycle conditions. PV Cre ;R26 tdTOM and Ntrk3 Cre ;R26 tdTOM mice were crossed from PV Cre (from The Jackson Laboratory, stock No: 017320, C57Bl/6J background) or Ntrk3 Cre (from MMRRC, stock No: 000364-UCD, C57Bl/6J background) and R26 tdTOM (Ai14, from The Jackson Laboratory, stock No: 007914, C57Bl/6J background) and used to genetically label SGNs and HCs for scRNAseq experiments. Neurod1 loxP/loxP (C57Bl/6 background) was published elsewhere 40 and Isl1 Cre was obtained from The Jackson Laboratory (stock No: 024242; C57Bl/6J background). Wild-type C57Bl/6J mice were obtained from The Jackson Laboratory (stock #000664) and used for most experiments unless otherwise specified.
Tissue collection and preparation. To obtain embryos, time-mating of the respective mouse line was performed. Pregnancy was verified by performing a plug-check on the following day. Plugs were assumed to occur at midnight, wherefore noon of the plug date was designated as embryonic day 0.5 (E0.5) and the date of birth was considered as postnatal day 0 (P0). For embryos, pregnant females were euthanized on E14.5 (Ntrk3 Cre ;R26 tdTOM ), E15.5, E16.5, E17.5 and E18.5 (all from PV Cre ;R26 tdTOM ) by CO 2 . After decapitation, the embryonic and neonatal tissues were processed depending on future use.
RNA in-situ hybridization and immunohistochemistry. For E16.5 embryos, the whole heads were processed. For E18.5 embryos to P3 pups, the cochleae were surgically dissected from the temporal bone under a stereomicroscope. Briefly, the skull was exposed and split into two halves by cutting along the ventral and dorsal axis from caudal to rostral. The brain and connective tissue were removed, exposing the inner ear, which was then dissected out. Following dissection, the whole head and cochlea were immediately fixed in fresh 4% paraformaldehyde (PFA, Sigma Aldrich) in PBS for 2 h or overnight (O/N) rolling at 4°C. After fixation, the tissue was washed three times in PBS for 10-15 min each and incubated in sucrose, rolling O/N at 4°C. Tissues were cryoprotected in 30% sucrose O/ N before embedding in OCT. Frozen tissues were kept at -80°C until sectioning. Blocs were sectioned at 14-16 µm with a Leica cryostat and the slides, kept at −20°C until further use.
RNA in situ hybridization experiments were performed using RNAscope ® . Paired double-Z oligonucleotide probes were designed by the manufacturer against target RNA and are available from Advanced Cell Diagnostics (Newark, CA). The RNAscope® Reagent Kit (Advanced Cell Diagnostics) was used according to the manufacturer's instructions (kit version 2). Frozen fixed tissue sections were prepared according to the manufacturer's recommendations. Each sample was quality controlled for RNA integrity with a probe specific to the housekeeping gene Ppib. Negative control background staining was evaluated using a probe specific to the bacterial DapB gene. The following probes were used in this study: Mm-Islr2- Image acquisition and analysis. Images were acquired using Zeiss confocal microscope LSM700, LSM800, LSM880 and LSM800 airy equipped with 5×, 10×, 20× and 40× objectives.
Single cell isolation. The same dissociation protocol was used for all stages. First, the presence of tdTomato signal was verified under a fluorescence stereomicroscope. From tdTomato positive animals, the cochlea was carefully dissected from the temporal bone as described above. Upon removal of the cochlear capsule, the spiral ganglia (SG) were dissected and collected on ice in Leibovit'z L-15 medium (Life technologies). For E14.5 and E15.5 tissue samples, the cochlea was cut open, and the modiolus region was dissected out. Tissue samples were then digested in a papain-DNAse solution (1.5 ml of papain at 1 mg/ml, 0.5 ml of DNAse at 0.1%) for 20 min, shaking at 700 RPM. After centrifugation of the partially digested tissue at 400 RCF for 10 min, the dissociation mix-solution was removed, and the cell pellet, gently resuspended in Dulbecco's modified Eagle's medium (DMEM) F-12 (Life technologies). Subsequently, cells were mechanically triturated using fire polished Pasteur pipettes coated with 0.2% bovine serum albumin until a homogenized solution was achieved. To remove residual cell aggregates and to generate a single cell suspension, the cell homogenate was then passed through a 70 μm nylon cell strainer (BD Biosciences). Single RFP + cells were sorted by fluorescence-activated cell sorting (FACS) into individual wells containing lysis buffer in a 384-well plate. The plates were immediately placed on dry ice and stored at −80°C before being processed for Smart-seq2 protocol. We Generation of count matrices, QC and filtering. 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 alignment. Reads where filtered for only uniquely mapped and were saved in BAM file format, count matrices were subsequently produced. Estimated count matrices from all plates were combined into one data object, QC metrics were computed using scanpy function 86 . Cells having more than 1000 detected genes and less than 5% of proportion of ERCC reads were kept. A median of 8469 genes were detected per cell in an initial analysis, and 7851 genes after glial code cleanup (see below) ( Supplementary Fig. 9).
Removal of glial contamination. First, developmental timepoints were analyzed separately. Genes being expressed in less than 3 cells were removed, count matrices were normalized per cell to a target sum of 1000 reads and then log1p was applied. High variable genes detection was performed using pagoda2 approach via scFates package (pp.find_overdispersed, default parameters). PCA was performed on the scaled matrix of over-dispersed genes (scanpy, default parameters). KNN graph (scanpy, pp.neighbors, n_neighbors = 30,n_pcs = 30, metric = "cosine") was generated from the PCA space, and was used as basis for cluster identification using leiden algorithm (scanpy, tl.leiden, default parameters), as well as UMAP embedding (scanpy, tl.umap, default parameters). First inspections of the expression profiles revealed glial contamination, as some clusters were mirrors of other neuronal clusters, while being positive for both neuronal and glial markers. E18.5 timepoint displayed the greatest amount of glial contamination and was used to calculate the glial code. A gene is considered part of the glial code if its expression has a correlation with Sox10 expression of more than 0.3. This threshold represents the best trade-off between removing the mirrored clusters while keeping most of the possibly informative genes (1371 genes were removed, see Supplementary Data 12 for gene names and their expression level).
Alignment of the timepoints and main analysis. The cleaned datasets were combined into one, and the same pipeline was employed as initial analysis, with different parameters for the KNN graph generation (scanpy, pp.neighbors, n_neighbors = 15, n_pcs = 15, metric = "cosine"). The first 15 PCs and developmental time annotation were then used for aligning the data using Harmony python package 17 . E17.5 and E18.5 were merged into one timepoint as there were too few cells from E17.5, which was affecting the results. Augmented affinity matrix was generated (core.augmented_affinity_matrix, n_neighbors=20). Diffusion maps was then generated from this affinity matrix (Palantir, run_diffusion_maps, default parameters) and multiscale space was determined (Palantir, determine_multisca-le_space, n_eigs=10) ( Supplementary Fig. 8). To generate the Force Atlas embedding, a t-SNE embedding was first generated from the multiscale diffusion space (scanpy, tl.tsne, perplexityt = 100, learning_rate = # of cells/12). Second, from the same multiscale diffusion space, a KNN graph was generated (scanpy, pp.neighbors, n_neighbors=30). Force atlas was then generated from the neighbors graph using the tSNE coordinates as initialization (scanpy, tl.draw_graph, init_pos = "X_tsne",maxiter = 500). Finally, for a selection of leiden clusters, cells being far away from their members on the FA embedding were considered doublets. To detect them, pairwise distances were calculated and standardized, cells hazing a z score distance of more than 3 were considered doublets. Differential gene expression was performed using scanpy-python package on this corrected count matrix, using Wilcoxon rank-sum test. We separately analyzed hair cells by reclustering the HC subset via pagoda2 and generated an UMAP embedding of the cells; vestibular hair cells were discarded from the analysis. SCENIC analysis. SCENIC pipeline was performed using the python package pySCENIC. First, the log-normalized count matrix was used as input, combined with a list of known TFs, to generate regulons based on correlation with putative target genes. Second, using the generated adjacency matrix combined with cis-Traget databases (mm10 500bpUp100Dw and TSS ± 10 kbp), the regulons were refined by pruning targets that do not present an enrichment for a corresponding motif of the TF. Third, cells were scored for each regulon with a measure of recovery of target genes from a given regulon.
Pseudotime tree inference. These steps were performed using scFates v0.4.0, a python package built in continuity of the crestree R package 87 . Trajectory inference was performed on SGN cells only.
First, 100 principal graphs composed of 400 nodes were generated with a different random initialization at each run, with SimplePPT approach on the multiscale diffusion space (scFates, Nodes = 400, method = "ppt", ppt_lambda = 1000, ppt_sigma = 0.2). While all trees merged the I b /I c clusters due to their high similarity, overlaying them onto FA embedding revealed two paths leading to this merged branch, one from immature and one going through other biasing (II and I a ). This led us to separately construct two trees to capture both I b and I c fates. A first principal graph composed of 600 nodes was fitted with SimplePPT approach on the multiscale diffusion space of a subset containing all clusters except 14 and 15 (scFates, tl.tree, Nodes = 600, method = "ppt", seed = 42, ppt_lambda=50, ppt_sigma = 0.15). A second 100-node principal graph, composed of only Cl. 14 and 15 was then fitted with the same parameters. The two graphs were then manually attached by the tips linking I a trajectory to I b trajectory. A root was automatically selected on the resulting merged tree, by selecting the tip with the lowest mean aggregated developmental time value (scFates, tl.root, tips_only = True, min_val=True). From the root was then calculated the pseudotime (scFates, tl.pseudotime, default parameters), from which was generated the dendrogram (scFates, tl.dendrogram, crowdedness = 0.2). Note also that a different Seurat based SNN algorithm was used previously to separate with success the two populations I b and I c at P3 1 .
Testing for features associated with the tree. Feature expression was modeled as a function of pseudotime in a branch-specific manner, using cubic spline regression exp i $ t i for each branch independently. This tree-dependent model is then compared with the unconstrained model exp i $ 1 using F-test. P-values were then corrected for multiple testing, features were considered significant if FDR < 0.0001. log10(fpm) count matrix was used to test which genes are significantly changing along the whole tree (scFates, tl.test_association, default parameters), with significant genes being fitted using GAM to obtain smoothed trends (scFates, tl.fit, default parameters). Whole tree was also used in combination with SCENIC derived AUC score, to detect significantly changing regulon activities. To do so, AUC scores were tested (scFates, test_association, A_cut = 0.025) and fitted (scFates, tl.fit, default parameters) via GAM.
Validating I b and I c separate trajectories. To validate that I b and I c biased cells can be separately fitted into two different branches (see Pseudotime tree inference section), both populations from E18 and P3 were fitted together with a single curved trajectory in diffusion space using ElPiGraph 88 (scFates, tl.tree, Nodes=30, epg_mu=200). Significantly changing genes along this trajectory were determined by testing for association separately for I b and I c cells, then by taking the union of significant genes (scFates, tl.test_association_covariate, A_cut = 0.5, fdr_cut = 0.01). This list of genes was then used for covariate testing, inspired by a recent preprint 89 . Genes were first tested for amplitude using the following GAM model: where s(.) denotes the penalized regression spline function and s(pseudotime):Covariate denotes interaction between the smoothed pseudotime and covariate terms. From this interaction term, p-values were extracted and then corrected for multiple testing (scFates, tl.test_covariate, fdr_cut=0.1). Genes were then tested for trend differences, comparing the model described in (1) to the following reduced one: Comparison was performed with ANOVA and p-values were corrected for multiple testing (scFates, tl.test_covariate, fdr_cut = 0.1).
Per trajectory analysis. For the per trajectory analysis, transitions from parts of the tree to endpoints were subsetted from the whole tree (scFates, tl.subset_tree) and a test for significance was reapplied (scFates, tl.test_association, A_cut = 0.3) to obtain genes changing on that part of the trajectory. Genes were considered specific to a trajectory if they were significantly changing exclusively in one of the analyzed trajectories.
Bifurcation analysis. Branch specific genes were first detected via amplitude testing using the following GAM model: From s(pseudotime):Branch interaction term, p values were extracted and then corrected for multiple testing (scFates, tl.test_fork, fdr_cut = 0.1). Then, each significant gene was tested for its upregulation along the path from progenitor to terminal state, using the linear model g i $ pseudotime. Differentially expressed genes were then assigned between two post-bifurcation branches with fdr < 0.05 and defined differences in expression cutoffs (scFates, tl.branch_specific, cutoffs were specifically set for each bifurcations). Finally, pseudotime of activation was estimated by separating the trajectory into 10 bins, and by calculating the relative expression rate at a specific bin: where f(b) is the mean fitted expression at a specific bin, if the rate was higher than a defined threshold, the gene was considered to activate at the pseudotime value of the related bin.
To analyze molecular mechanisms of cell fate biasing, 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 (window sizes were 90 cells for bifurcation A, and 50 cells for the others). 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 genegene correlations of all possible gene pairs inside one module, or between the two modules, respectively (scFates, tl.slide_cors, default parameters).
Cell to cell communication. First, cell to cell communication was performed using CellphoneDB python package 56 using leiden clustering and corrected log10(fpm) matrix as inputs. As this package uses human databases of interactions, genes were converted to human format using biomaRt database. Pipeline was performed on main E18.5 SGn and HC leiden clusters with default parameters, with means and p values for each pair of interaction between two clusters being the output. For the main figure, manual curation of the results was performed using biological knowledge.
Statistics and reproducibility. Statistical analyses were performed using Graph-Pad Prism 8. All validations of sequencing results by in situ hybridization and immunostaining were replicated across at least 5 sections from multiple animals. All micrographs are representative images.
Reporting summary. Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Data availability
Raw sequencing data are available on GEO database under accession code GSE165502 (GSE165502), pagoda2 web application can be explored on a browser via the following link: https://adameykolab.srv.meduniwien.ac.at/SGN/. The data file is also available on the Lallemend laboratory website (lallemend-laboratory) or via the link lfaure/SGN. Moreover, all the analyzed data are available by browsing and analysis via the gene Expression Analysis Resource (umgear). Source data are provided with this paper.

Code availability
Tree fitting, pseudotime and bifurcation analysis was performed using scFates v0.4.0 python package, available via pypi: scFates. All codes and data for downstream analysis are deposited on the following github repo: sgnfates.