Single-cell transcriptional logic of cell-fate specification and axon guidance in early-born retinal neurons

Understanding how diverse neuronal cell types emerge in the mammalian central nervous system (CNS) is essential to delineate the logic of neuronal network assembly. The timing of production of neurons is crucial to the establishment of functionally efficient networks (Rossi et al., 2017). Single-cell RNA-seq has now allowed the identification of an increasing number of cell types (Poulin et al., 2016; Tasic et al., 2016; Yuan et al., 2018; Zeisel et al., 2015), but untangling the transcriptional features involved in the chronological generation of diverse classes of neurons remains a conundrum. The mouse retina was one of the first CNS tissues analyzed using single-cell transcriptomics, establishing its status as a model of sequencing analysis (Macosko et al., 2015; Shekhar et al., 2016; Trimarchi et al., 2008). The retina also represents an excellent model for deciphering how neuronal types can emerge from a single pool of progenitor cells, as it is one of the simpler parts of the CNS with only six classes of neurons, all having been extensively characterized both molecularly and morphologically (Cherry et al., 2009; Livesey and Cepko, 2001; Sanes and Masland, 2015).

The retinal classes of neurons are produced in two waves (Rapaport et al., 2004). The first wave gives rise to the early-born cell types from embryonic days (E) 10 to 17, which comprise the retinal ganglion cells (RGCs), the horizontal cells (HCs), the amacrine cells (ACs) and the cones. In the second wave, ACs are still produced, together with bipolar cells (BCs) and rod photoreceptors, from E14 to postnatal day 5. Recently, BCs were shown to consist of distinct transcriptionally defined subclasses in the adult mouse, which were associated with morphological features (Shekhar et al., 2016). With more than 30 subtypes based on their dendritic morphologies, early-born RGCs are the sole output from the retina to the brain, reaching up to 46 targets in the mouse brain (Martersteck et al., 2017; Morin and Studholme, 2014; Rivlin-Etzion et al., 2011; Seabrook et al., 2017). Recent work in adult and juvenile RGCs revealed a high degree of transcriptional heterogeneity (Macosko et al., 2015; Rheaume et al., 2018). However, the manner in which retinal progenitors give rise to this extreme heterogeneity is not fully understood.

Seminal work using vector lineage tracing identified fundamental aspects of the logic that allows various retinal cell fates (Turner and Cepko, 1987; Turner et al., 1990; Wetts and Fraser, 1988; Wetts et al., 1989). This process involves intrinsic components, with sequential expression of transcription factors playing a key role in driving competence in progenitors to generate distinct cell types (Cayouette et al., 2003; Cepko, 2014). Recent advances in single-cell genomics have proven the capacity to delineate neuronal lineages using either transcriptomic profiles (single-cell RNA-seq) or enhancer signatures (ATAC-seq) (Kester and van Oudenaarden, 2018; La Manno et al., 2018).

Here, we used single-cell transcriptomic reconstructions to unveil the programs at work in the early specification of mouse retinal neurons. Taking advantage of their temporally organized production, couples of ligand-receptors were identified as putative axon guidance pairs that can guide RGC growing axons on their paths to the brain. Exploiting the spatial distribution of RGCs, dorso-ventral and temporo-nasal scores were attributed to each RGC, a refinement that allowed us to identify RGCs from the ventro-temporal crescent of the retina, in which was found a significant subset of cells expressing ipsilateral-projecting RGC genes. Finally, this analysis was validated by isolating and deep-sequencing ipsilateral-projecting RGCs, giving new insight into the molecular identity of this RGC subset that expresses specific axon guidance transcriptional programs.

The goal of this study was threefold: to trace the origin of early-born retinal cell fates, to infer the spatial relationships across retinal neurons, and, finally, to identify RGC groups with transcriptional signatures linked to different projections in the developing brain.

We produced transcriptional profiles of early-born retinal progenitors obtained from 5348 single cells of E15.5 retina using 10x Genomics (Fig. 1A,B). Cells were distributed in 14 clusters in a t-distributed stochastic neighbor embedding (t-SNE) analysis obtained by merging two replicates (Fig. 1C, Fig. S1). Each cluster was characterized using known marker genes (Fig. 1D-I, Fig. S2) (Bassett and Wallace, 2012). The central clusters (0-3) were composed of retinal cycling cells that were referred to as the retinal progenitor cells (RPCs) based on the expression of Sox2, Fos and Hes1 (Fig. 1D). Emanating from this core unit, we identified a narrower group (cluster 4) expressing cell cycle exit genes (Top2a, Prc1) as well as the neuronal-specific genes Sstr2, Penk and Btg2 (Fig. 1E, Fig. S2A,B). The latter gene was previously shown to induce neuronal differentiation (el-Ghissassi et al., 2002) and was identified in retinal progenitors (Trimarchi et al., 2008). These cycling cells progressed to a cluster aggregating early neuroblast transcription factors including Neurod4 and Pax6 (cluster 5) as well as genes known to promote axonal growth (Pcdh17) (Fig. S2C-E) (Hayashi et al., 2014). Following this neuroblast bottleneck furrow were three distinct branches. The main one (clusters 9-12) was characterized by RGC markers (Isl1, Pou4f2, Pou6f2 and Elavl4) (Fig. 1F, Fig. S2M-O). The second neuronal group (clusters 7-8) was composed of both ACs in the root part (Onecut2⁺, Prox1⁺) and HC at the extremity (cluster 8, Onecut1⁺, Lhx1⁺) (Fig. 1G, Fig. S2G,H). The third (cluster 6) was positive for Otx2, Crx and the early marker Thrb, a signature for cones (Fig. 1H, Fig. S2I,J). In addition to Thrb and Crx, Rbp4 transcripts were also highly enriched (Fig. S2J,K). The specific expression of Rbp4 was confirmed using mice expressing Cre under the control of the Rbp4 promoter (Fig. S3). Both morphologies and positions of Rbp4⁺ cells followed cone hallmarks, as at this stage they have a distinctive morphology aligned along the apical side of the retina (Decembrini et al., 2017). The smallest cluster (2.9% of the cells) was positive for mitochondrial genes and lacking Rps/Rpl genes (cluster 13) (Fig. 1I, Table S1); this cluster was designated as unknown/RGC-like (U/RGC) as it is positive for RGC markers (Pou6f2, Pou4f1, Isl1, Islr2, Syt4, Ebf1/3 and L1cam) but these cells may be RGCs with poor viability outcome, or originating from an alternative source such as the ciliary margin zone (Bélanger et al., 2017; Marcucci et al., 2016). These main neuronal clusters were validated with in situ hybridization (ISH) and immunohistochemistry (Fig. S3).

Finally, all clusters were represented on a heatmap with their top 10 differentially expressed genes (Fig. 1J, top 15 in Table S1). This clustering respected neuronal type similitudes and highlighted the major transcriptional modules defining the differentiation process into early retinal cell types (Fig. 1K). Together, these results provide a comprehensive characterization of all the early neurons with their transcriptional signatures, thus allowing further exploration of each cell type.

The global organization of the retinal clusters, as represented in Fig. 1C, was colored with the clusters identified by cell type, and the RGC groups split into young, mid1, mid2 and old RGC based on their expression profile, resembling a lineage tree where most neuronal fates emerge from a common neuroblastic state (Fig. 2A). In the retina, neural fates follow a gradual shift of competences (Boije et al., 2014), a continuum of transcriptional diversities that was represented using a uniform manifold approximation and projection (UMAP), which facilitates the reconnection of divergent clusters (Fig. 2B, Movie 1) (McInnes et al., 2018). We observed that the progenitor pool is organized by cell cycle phases (Fig. S4). Because neurons originated from the progenitor branch, we hypothesized that the extensions of these clusters were associated with their birthdate. As it was shown that the progenitor-to-neuron differentiation follows pseudo-timed incremental steps in the days following cell cycle exit (Telley et al., 2016), we asked whether the extension of the clusters observed in the t-SNE and the UMAP would match their order along a pseudo-time axis (Fig. 2C). A progression from the pool of RPCs toward each of the clusters was observed, with the first branch (purple) giving rise to photoreceptors (PRs), the second branch to AC/HCs and finally the most mature group, the RGCs. Although this is an E15.5 snapshot visualization of retinal development, the transcriptional continuum reflects similarities in the mode of production of the different neuronal types.

Because early-born retinal neurons are produced from multipotent progenitors going through asymmetric divisions (see Fig. 2D illustrating known modes of production) (Belliveau and Cepko, 1999; Livesey and Cepko, 2001; Turner et al., 1990), we explored the hierarchical mode of production to unveil the establishment of generic programs that may drive the main differentiation paths. To achieve this, a reconstruction of linear and branched trajectories based on Monocle 2 was performed (Fig. S5A). This tree representation only reflects common transcriptional shifts in the generation of the four cell types, with five nodes found with drastic changes in transcriptional programs oriented toward the maturation axis of the main clusters (Fig. S5A,B). Whereas specific programs were oriented toward the PRs (node 3), the ACs and HCs were found to follow similar trajectories (node 4), indicating shared transcriptional dynamics during their maturation (Fig. 2C, Fig. S5). Further down the tree, RGCs segregated into two groups (node 5): one from a pool of neuroblasts expressing Stmn1 and Sncg (right branch), and one containing cells that were denoted as U/RGC cluster (many mitochondrial genes; left branch) (Figs S3 and S5). In order to link cells from a progenitor pool to their corresponding fate, branched expression analysis modeling (BEAM) was used to trace back the genes correlated with the temporal transitions along differentiation paths (Hanchate et al., 2015; Qiu et al., 2017b). Dynamic expressions of the genes identified by BEAM were the most likely to contribute to the balance toward one fate. For node 3, Casz1, Thrb and Meis2 were among the top cone-oriented transcripts, whereas Stmn2 and Klf7 were scored toward the RGC/AC/HC fates (Fig. 2E, Fig. S5). Node 4 showed a misbalance with Onecut2, Casz1 and Ldhb shifted to AC/HC group whereas Stmn3 and Elavl4 were scored as specific to RGC fates (Fig. S5). Interestingly, we observed unbalanced expression of Prdm1 transcripts for cones, and two miRNAs, miR124a-1hg and miR124-2hg, for which preferential expression pairs with AC/HC fates (Fig. S6). Finally, mitochondrial genes were scored in the U/RGC-like branch versus Ptma, H3f3a, H3f3b and Sncg in the RGC-enriched branch. Expression patterns of some of the influential transcripts were validated by ISH (Fig. 2F). Importantly, node 3 represents the junction from the progenitor pool to one of the four neuronal cell types and appears to emerge from a conserved transcriptional root. Overall, this analysis reveals transcriptional sequences being shared across the early-born retinal cell types.

Because the convergence in the neuroblast pool is accompanied by a reduction of heterogeneity in the t-SNE and the UMAP (Fig. 2A,B), we sought to explore the transcriptional dynamic by extrapolating a time depth by means of mRNA maturity levels with the RNA-velocity method (La Manno et al., 2018). Based on the ratio between nascent and mature mRNAs, RNA-velocity analysis allowed us to tag each cell by a velocity vector corresponding to its putative near future transcriptional state, further validating the progenitor-to-neuron organization (Fig. 3, Fig. S7). The transcriptional dynamics were particularly decreased at the intersection prior to neuronal cell-type specification (Fig. 3A). Low velocities were also observed at the root of the progenitor cluster, followed by robust directional flow towards each of the neuronal branches, finally slowing down again at their extremities (ends) (Fig. 3A, right). A striking feature of the progenitors was a robust wheel with well-defined progression of all cells transitioning from one cluster to another. Across the circular pattern we detected a step-wise expression of cell cycle genes (Fig. 3B). Analyzing this process using one-dimensional trajectory modeling revealed four transcriptional waves corresponding to cell cycle phases (G1, G1/S, S/G2, G2/M, M) (Fig. 3C). This dominant organization matched with a core set of cell cycle genes (Fig. 3D). We reasoned that gene dynamics segregating apart from this main cell cycle force may orient progenitors toward one of the main cell-type paths. To extract this information, we unweighted the cell cycle components and re-clustered the progenitor pool into six clusters (Fig. 3E). In this framework, we selected genes that over this trajectory expressed high variance (Fig. 3F). Among those were found early genes such as Odf2 (PR specification) and Vax2os (ventral retina transcript also enriched in rods PR; Corbo et al., 2007). ISH validated the spatially restricted expression of Crym, Slc1a3, Pak1 and Vax2os (Fig. S7C). Importantly, we found Vsx2, a transcript associated with retinal progenitor proliferation and the production of a late-born cell type: the BC (Fig. 3C,I, Fig. S7B).

In order to see whether these genes would depict a fate-oriented dynamic, RNA-velocity analysis was performed on post-mitotic neurons (Fig. 3G). An increased number of complexities emerged in which we observed preferential path for the transcriptional trajectories. Prdm1, Otx2 and Neurod1 were found enriched in their unspliced form in the G2/M group of cells with an orientation leading to PR fates (Fig. 3H). Conversely, Pou factors (Pou2f2, Pou3f1, Pou3f2, Pou4f1, Pou4f2, Pou4f3 and Pou6f2), Ccer2 and Cntn2 followed either RGC or AC/HC paths. To reconcile gene dynamics from both progenitors and neurons, we showed the profiles of velocity variances with their global expression levels for the genes linking the two groups together, with most of the genes bursting in the bottleneck prior to cell-type branches (Fig. 3I). In this zone, we took three neighboring cells and showed that we could assign them predictive trajectories as transition probabilities of joining one of the three main cell types (RGC, PR or AC/HC) (Fig. 3J).

Next, our analysis focused on RGCs to study in-depth the transcriptional programs driving their differentiation. In order to sub-classify this large group of 1312 cells, we first performed a principal component analysis (PCA) to identify the main genes responsible for their transcriptional heterogeneity (Fig. 4A). We found that the principal gene ontology (GO) categories segregating the RGCs were ‘cell differentiation’ (PC1), ‘nervous system development and transcription regulation’ (PC2) and ‘cell differentiation and survival’ (PC3). PC1 was responsible for >7% of the variance with an apparent temporal progression allowing compartmentalization of the RGCs based on their differentiation level, i.e. reflecting their birthdate within the retina. This parameter allowed the identification of gene expression patterns along the temporal progression with the demarcation of six transcriptional waves from 383 genes (Fig. 4B). These genes are important candidates for the coordination of the differentiation of RGCs, including six well-characterized neuronal markers that we plotted on an RGC UMAP: Dlx1, Pou4f1, Tshz2, Pou6f2, Dnm3 and Eomes (also known as Tbr2) (Fig. 4C,D, Table S2). Fifty-eight of these genes encode transcription factors that are highly interconnected in gene networks and segregated into three groups (Fig. 4E). The first group are cell cycle exit genes, with Fos, Jun and Notch1/2 in its core (Fig. 4E, top cluster, green), mostly connected to a second group of developmental transcription factors (such as Neurod1 and Isl1, previously shown to specify RGCs; middle cluster, pink), themselves slightly connected through Pou4f2 to a more distant cluster of genes that pattern synaptic terminals (e.g. Snap25), representing an advanced level of maturation (Fig. 4E). Taking Pou4f2 and Eomes, two genes from the starting and final waves, respectively, we observed a peripheral and neuroblast layer-specific pattern (Fig. 4F). We produced the same analytical framework for the PRs in which we found a similar organization with four transcriptional waves (Fig. S8). Overall, this analysis showed that, at this stage (E15.5), we can recover a large spectrum of differentiation stages of early retinal cell types.

As retinal neurogenesis follows a central-to-peripheral gradient, we reasoned that molecules that are surface-bound or secreted in higher amounts by mature RGCs located in the central retina would be the best candidates to attract axons of the more recently produced RGCs located at the periphery (Fig. 5A). Conversely, cues produced in higher amounts by young RGCs may play a repulsive role in an autocrine manner for their axons fleeing toward the optic nerve entry in the central part of the retina. In order to classify the central-to-peripheral position of RGC cells, we took advantage of our differentiation analysis to attribute a ‘birthday’ status to each of them.

Eleven ligand-receptor pairs were identified (Fig. 5B). Two of them (Slit2/Robo2 and Dcc/netrin) are known to play a role in these guidance steps (Fig. 5B,C, Table S4) (Deiner et al., 1997; Thompson et al., 2009). Among the ‘old’ genes, Igf1 was found as a positive control, a secreted ligand shown to be strongly expressed in the most central part of the retina (Wang et al., 2016). ‘Young’ genes included Igfbp5 and Igfbpl1, two essential components of IGF1 signaling (Fig. 5C, Tables S3, S4) (Guo et al., 2018). Among other secreted ligands, reelin mRNA was expressed in a gradual manner from periphery to center, as did Kitl transcripts (Fig. 5B-E). The expression of Kit, its receptor, was validated by immunohistochemistry to be expressed by peripheral RGCs with a strong enrichment on RGC axons (arrowheads) growing across more mature RGCs (EOMES⁺) (Fig. 5D). Of note was the strong detection in RGCs of Tenm3 and EphB1, two genes involved in the segregation of ipsilateral RGCs at the optic chiasm (Leamey et al., 2007; Williams et al., 2003). Beyond the molecules known to play a direct role in axon guidance were also identified genes encoding for putative downstream machinery, including the central-high gene Pde4d and the periphery-enriched gene Btg2 (Tables S3, S4). Complementary patterns for Slit2/Robo2 and the netrin 1 receptors (DCC, Unc5) were plotted in the same t-SNE (Fig. 5C). Together, these results validate previously known genes and identify new ones with expression patterns that suggest a role in retinal axon guidance.

The most striking axon guidance choice for RGCs occurs at the optic chiasm. Axons that do not cross the chiasm originate from ipsilaterally projecting RGCs (I-RGCs), which have cell bodies in the ventro-temporal (VT) part of the retina. As E15.5 is the peak stage for axon segregation (Colello and Guillery, 1990; Drager, 1985; Guillery et al., 1995), RGCs projecting contralaterally (C-RGCs) and ipsilaterally may have distinct molecular profiles. To capture the molecular profiles of I-RGCs, we proceeded in two steps. First, we exploited their prominent spatial feature, i.e. their position along the dorso-ventral (DV) and temporo-nasal (TN) axis (McLaughlin and O'Leary, 2005), to identify genes enriched in the VT RGCs. To do so, RGCs were positioned in a pseudo-spatial orientation (see Materials and Methods). From these micro-clusters, transcriptional signatures of each quadrant were extracted (Fig. 6). Second, as the VT is a mixture of both I- and C-RGCs, an ipsilateral-specific gene, Slc6a4, was used to recover the distinct subpopulations without any spatial considerations (Fig. 7).

To visualize whether spatially distributed genes can explain the genetic diversity, cumulative expression of known DV and TN genes were plotted in tSNEs from five dimensions. In a configuration in which spatial patterns aggregate was found (Fig. 6A,B; see Materials and Methods), we attributed VT, dorso-nasal (DN), dorso-temporal (DT) and ventro-nasal (VN) scores to every cell that was mapped (Fig. 6C). Because the aggregation was not completely obvious, we ordered the cells in two ordinal axes with color-coded maturation index (Fig. 6D). This classification allowed a comparison between groups of transcriptomes between the quadrants. We identified 92 VT single cells, with 19 genes enriched in these cells, including Crabp1, Nefl, Gal, Eno1, Cbx1, Rassf4 and Btg2 (Fig. 6E). When comparing temporal versus nasal cells, we detected >1.5 enrichment for five temporal genes (including Crabp1, EphB1 and Gal) and five nasal genes (including Cntn2 and Pou4f1) (Table S5). When comparing dorsal and ventral RGCs, we identified 423 genes that were differentially expressed (Table S5) from which six were validated by ISH: three dorsal genes (Fstl4, Cnr1 and Nr2f2) and three ventral genes (Cntn2, Irx2 and Vax2os) (Fig. S9). Interestingly, the ventrally enriched gene Unc5c was reported earlier this year to be specific to the ventral retina and play a role in RGC axon guidance (Murcia-Belmonte et al., 2019).

Next, to understand how the VT genes would be linked to previously known ipsilaterally associated genes, we used an ipsilateral-specific gene, Slc6a4, to define subpopulations without any spatial considerations (Fig. 7). A Cre mouse line was used to induce the expression of tdTomato signals in Slc6a4-expressing cells, which is specific to I-RGC at E15.5 (Slc6a4-Cre x tdTomato Ai14) (Peng et al., 2018; Zhuang et al., 2005). tdTomato-positive and -negative cells (tdT⁺ and tdT⁻) were sorted by fluorescence-activated cell sorting, captured in a C1-chip and deep-sequenced at about 500,000 reads/cell to reach optimal sensitivity to detect specific genes expressed by this rare subpopulation (Fig. 7A).

We scored 719 cells with specificity for both Slc6a4 and tdTomato transcripts in the tdT⁺ fraction (Fig. 7B), and found a large fraction of RGCs but also progenitors and non-RGC neurons (Fig. 7C,D). Among the RGCs, 148 were positive for Slc6a4, which we compared with RGCs belonging to the same maturation window those expressing Slc6a4 (Fig. 7E). Differential expression analysis between Slc6a4⁺/tdT⁺ and Slc6a4⁻/tdT⁻ cells showed an enrichment of 15 I-RGC-specific genes and nine C-RGC genes (logFC>0.4=FC>1.5; P<0.05) (Fig. 7F,G, Fig. S10, Table S6). Among the I-RGC-specific genes we found Slc6a4, but also EphB1, a gene known for its role in the segregation of I-RGC axons (Petros et al., 2009, 1; Williams et al., 2003). Interestingly, we found three genes coding for neuropeptides: Tac1, Sst and Gal. ISH demonstrated SST expression specifically in the ventral part of the retina toward the periphery (Fig. 7H, Fig. S11). At the single-cell level, we observed that neuropeptide genes were only co-expressed in a subset of the Slc6a4⁺ cells (Fig. 7H,I). Finally, analysis of other known I-RGC genes in a supervised clustering showed their covariance in a Spearman correlation matrix (Fig. 7J). Among these genes, together with Slc6a4, were Zic2, Zic1, FoxD1, EphB1, Tenm2 and Tenm3 (Herrera et al., 2003, 2; Herrera et al., 2004, 1; Leamey et al., 2007; Wang et al., 2016; Williams et al., 2003, 1; Young et al., 2013). A central cluster with higher correlations for ten genes, Slc6a4, Zic5, Gal, Tac1, Lgi2, Rdh10, EphB1, Pcsk2, Crabp1 and Sst, was identified further validating the association of neuropeptides with the I-RGCs. Surprisingly, Tenm2 and Tenm3 were segregated with Isl2, a transcription factor known for its role in patterning C-RGCs located within the VT retina (Pak et al., 2004).

As a control, a co-expression analysis was performed to show the extent to which they are co-expressed in wild-type single cells using the drop-seq procedure. Among the I-RGC genes, eight were sufficiently expressed (> 0.5% RGCs) to be included in the analysis: Slc6a4, EphB1, Tenm2, Tenm3, Zic1, Zic2, FoxD1 and Ptch2 (Clark et al., 2019) (Fig. S12). Among them, very few overlaps were found in single cells, as shown with the six most highly expressed genes (Fig. S13). Among nine genes found co-expressed in I-RGCs (including Crabp1, Zic5, Ephb1, Gal and Tac1) Crabp1 was previously shown to be strongly enriched in the VT (Díaz et al., 2003), strengthening its potential role in I-RGC maturation. Altogether, this in-depth analysis provides a new set of genes with specificity of expression in a well-defined RGC subpopulation.

Here, we provide a resource of 6067 retinal single-cell transcriptomes that allows the stepwise progression of gene expression during the generation of early-born retinal cells to be deciphered. We used RGCs, the most populated and robust group, to delineate the logic of this progressive maturation. We discovered that transcriptionally inferred spatial information can be used to define subgroups of RGCs, leading to the identification of novel candidate genes that may play a role in their identity and connectivity.

Our cell fate-oriented analysis consolidated the existence of a transcriptional logic associated with the hierarchical production of early-born retinal cell types. It is important to note that on its own it does not directly imply causal lineage relationships. However, searching for candidate genes that are biased along binary fate decisions, we observed a segregation of progenitors giving rise to the four classes of early-born neurons versus those giving rise to other cell types (Fig. 2C, Fig. S5A). This aligns well with previous findings showing a shift in competence between the two phases of production, a phenomenon shown to be intrinsically driven (Belliveau and Cepko, 1999; Reh and Kljavin, 1989). In the PR branch, which appears homogeneous, we further emphasize a role for Otx2, Neurod4, Casz1 and Rpgrip1 in the generation of PRs (Emerson et al., 2013, 2; Hameed et al., 2003; Mattar et al., 2015). Primary activation of Prdm1, Otx2, Cplx2 and Neurod4 in progenitors was found to be a potential predictor of early fate bias towards the PR (Fig. 3, Fig. S5). The progenitors upstream of node 3 seem to have the competence to generate both cones and other early-born cell types such as HCs and RGCs. This latter observation supports previous lineage analysis from labeled clones in zebrafish, in which a common progenitor can give rise to both RGCs and cones (Poggi et al., 2005), or HCs and cones, whereas ACs and HCs seem to originate from a common transcriptional program (Boije et al., 2014). In our study, Olig2⁺ postmitotic neurons had trajectories coinciding with both Ascl1 and Prdm1, two genes associated with ACs and cones, thus supporting evidence in mouse where such a path was found for Olig2 expression history (Hafler et al., 2012) (Fig. S14). Finally, post-mitotic neuroblasts giving rise to RGCs tend to segregate in two groups (Fig. 2C, Fig. S5A), possibly due to distinct origins, since lineage tracing and time-lapse microscopy have shown in transgenic mice that a subset of RGCs originates from the ciliary margin zone (Marcucci et al., 2016; Wetts et al., 1989).

Relaxing chromatin prior to transcriptional activation is a balance between transcription factors and chromatin remodelers. Our study reveals upstream actors dynamically expressed during the initial neuronal maturation, including chromatin remodelers such as Hdacs, Kdm6b, Ezh1/2, Hmga1/2 (ACs), Hmgb2, and several Prdm transcripts known to regulate transcription (Fig. 3, Figs S6, Tables S2, S3). Their expression pattern precedes the remodeling of chromatin and may serve as a permissive signal to restrict multipotency and cell-fate commitment (Chen and Cepko, 2007; Fabre et al., 2015, 2018; Iida et al., 2015, 2; Kohwi et al., 2013; Pereira et al., 2010, 2). Among the Prdm genes (with catalytic SET methyltransferase domain), whereas Prdm1 and Prdm13 are found specific to cones and AC cluster, respectively, Prdm2, Prdm10 and Prdm16 define RGC subpopulations (Groman-Lupa et al., 2017). Interestingly, we saw enrichment of metabolic genes (Ldhb in AC/HCs and RGCs, Ldha in PRs). This may highlight involvement in aerobic glycolysis in cones, allowing them to meet their high anabolic needs during their maturation (Chinchore et al., 2017; Zheng et al., 2016) (Table S1, Fig. S6). In the U/RGC cluster we observed a misbalance of mitochondrial transcripts that could indicate unhealthiness (Fig. S6). However, this could reflect increased numbers of mitochondria linked to a differential mode of metabolism or asymmetric distribution of aged mitochondria that may influence the stemness of progenitors (Katajisto et al., 2015). How this differential expression influences RGC fate or origin remains to be established.

To follow the establishment of cell-fate specification, it is essential to delineate the chronological order in which transcription factors are produced, i.e. to understand their hierarchical organization in order to reconstruct gene regulatory networks. In this study, we took advantage of a recently developed procedure to unveil forms of pre-mRNA produced by a cell giving a prediction that can be modeled by UMAP (Fig. 3, Fig. S7) (La Manno et al., 2018). Neuroblasts access a transcriptional bottleneck with velocity vectors that decrease near the main branching area, implying a collective loss of expression dynamics followed by secondary multi-directional divergence, as shown for oligodendrocyte precursor cells in the neocortex (Zeisel et al., 2018). Using this framework, we established a chronological sequence between transcription factors, with associated transition probabilities that, although only predictive, will be key for further studies to validate post-mitotic trajectories prior to their definite cell-fate acquisition.

Our results led us to postulate that older neurons that project in an untracked environment are equipped with a particular set of guidance receptors in order to reach their target. One well-characterized neuronal cell type that is produced early, with diverse, long-range projections, is the RGC (Petros et al., 2008; Sanes and Masland, 2015; Seabrook et al., 2017; Trimarchi et al., 2008). We hypothesized that E15.5 RGCs are already connected with distant brain territories, where they are known to segregate based on their location in the retina (Seabrook et al., 2017). Taking advantage of our DV/TN classification, we compared their transcriptional signatures (Fig. 6D,E). Previous studies showed that retinotopy of the retina is associated with differential expression of ephrin genes (Cang et al., 2008; McLaughlin and O'Leary, 2005). Validation of these groups of cells was exemplified by Crabp1, a gene shown to be highly enriched in the ventral segment of the retina (Díaz et al., 2003). Beside EphB1 and Crabp1, we also identified Gal and Nefl as strongly enriched in VT RGCs. These four genes were also found to be upregulated in our second approach using the selection of Slc6a4⁺ RGCs, reinforcing their putative role in the development of I-RGCs and supporting the validity of our approach based on spatial inferences.

Our study led to the identification of a dozen of ligand/receptors, including a few for which functional roles were validated in mutant retina (Fig. 5B). For instance, the peripheral-Slit2/central-Robo2 complementary pattern, already detected by ISH (Erskine et al., 2000), was implicated in the fine-wiring of RGC axons on their path toward the optic nerve, exclusively affecting axon guidance within the peripheral retina (Erskine et al., 2000; Thompson et al., 2006). This phenotype now seems compatible with a lack of repulsion from the periphery-to-middle part of the retina, where attractive cues may take the relay to stir RGC axons.

Two of the ligand-receptor systems we identified are the Kitl/Kit pair and netrin 1/Dcc/Dscam/Unc5c signaling (Fig. 5B-D). The former pair encodes for SCF and its receptor and has never been shown to be involved in guiding retinal neurons, thus further functional characterization will be required to confirm this effect (Williams et al., 1990). However, their function as axon outgrowth-promoting cues has been established in the peripheral nervous system (Gore et al., 2008; Hirata et al., 1993). In the retina, as expression of the ligand (Kitl) increases along the RGC maturation axis (Fig. 5C), expression of its receptor (Kit) is mostly detected in younger RGCs. Interestingly, the combined and gradual expression of Adam10, which encodes a protease that processes SCF in its soluble form, from the early-to-late RGCs, may participate in the gradient formation of SCF along the central-to-periphery domains. Unc5/Dcc/Dscam receive signals from Ntn1-expressing cells from the optic nerve entry (Deiner et al., 1997), with Dcc strongly expressed in ‘young’ RGCs, at the time that they send axons toward the center of the retina (Fig. 5C). Although the outgrowth-promoting activity of netrin 1 was shown to be abolished by a DCC-blocking antibody (Deiner et al., 1997), it was unclear then how RGC axons would escape this netrin attraction to pursue their route in the optic nerve. Axons have opposite turning responses to netrin 1 depending on the status of cytosolic cAMP-dependent activity (Ming et al., 1997), and the switch in netrin 1 activity was then shown to rely on the age of RGCs (Shewan et al., 2002). In our study, we identified two other netrin receptors that may influence this switch. Whereas Dcc is expressed in young RGCs, Unc5c and Dscam, encoding for two receptors that have been shown to heterodimerize to mediate growth cone collapse (Purohit et al., 2012), are switched on while RGCs mature, concomitant with the extension of their axons in the optic nerve (Fig. 5C). How the switch from Dcc to Unc5c/Dscam is orchestrated remains to be established. Whether the other pairs identified in this study play an actual role also remains to be validated by a combination of in vitro and in vivo functional studies.

Our analysis of the transcriptional signature of ipsilateral RGCs led to the identification of ten genes using complementary methodologies (Fig. 7F-J, Fig. S12A,B), including a set of neuropeptide-encoding genes, potentially implicating their role in the maturation of I-RGCs. Importantly, all of the genes identified in our study except EphB1 differ from those found in a recent study that used retrograde tracing to distinguish crossed from uncrossed RGCs (Wang et al., 2016). However, as I- and C-RGC production timing differ significantly, a caveat of that study was the difficultyof comparing RGCs from the same maturation level, leading to the identification of many factors reflecting the immature state of I-RGCs (Marcucci et al., 2019). To overcome this issue, we extracted transcriptomes from single cells for which we could classify the maturation level. Furthermore, we used Slc6a4 as an early marker of ipsilateral RGCs to identify the transcriptional set of mRNA that accumulated at the time of (un)crossing. Although we failed to detect Zic2, a transcription factor first identified for its patterning of ipsilateral RGCs (Herrera et al., 2003), we did detect Zic5, a gene flanking Zic2 and thus sharing many regulatory elements that may be responsible for its expression in I-RGCs. Of note, other genes involved in I-RGC guidance, such as Boc (Fabre et al., 2010; Peng et al., 2018) and FoxD1 (Carreres et al., 2011, 1; Herrera et al., 2004), were not detected in a sufficient number of cells. Along with Zic2, these genes might be either expressed at levels below detection, or subjected to high drop-out events, a caveat of the single-cell procedure. Of note, Zic5 and Zic1 (also found enriched in I-RGCs) were detected with both approaches (C1 and 10x) in more cells than Zic2, suggesting either a higher rate of drop-out events or a more transient expression. The latter possibility may relate to the fact that Zic2 is expressed in bursts to control the transient expression of Slc6a4 in I-RGCs (García-Frigola and Herrera, 2010). In agreement with our spatial inference in the RGCs from the VT, our analysis showed enrichment in Crabp1, the most specific factor that we found in Slc6a4⁺ cells in both our C1 and 10x experiment. Altogether, these new genes represent potential candidates that play a role in the establishment of I-RGC connectivity. However, functional validation of the dozens of candidate genes found in this study will be crucial to establish their exact role during RGC maturation and axon guidance.

This study provides a fundamental resource of single-cell transcriptomes in the developing retina, showing that retinal progenitors exhibit a high level of transcriptional heterogeneity, and unveils its meaning with the identification of all early-born cell types. Moreover, fine-scale diversity can be resolved among these cell types. For RGCs, the identification of differentiation waves, spatial position and meaningful genes for functions such as patterning and circuit formation can be recovered. This study will facilitate our understanding of how the retina develops, especially in terms of cell-fate specification for early-born neurons. It paves the way for functional perturbations related to sequential generation of cones and RGCs, and the transcriptional dynamics leading to cell differentiation. Finally, the refinement of gene networks and co-expression analysis has revealed deep relationships between retinal genes, and thus represents an important basis for a better understanding of retinal development at the single-cell level.

Experiments were performed using C57Bl/6J (Charles River), and Ai14 Cre reporter (Jackson Laboratory #007914), B6.129(Cg)-Slc6a4tm1(cre)Xz/J (Stock No: 01455) and Rbp4-Cre (Tg(Rbp4-Cre)KL100Gsat/Mmcd; GENSAT RP24-285K21) mice bred on a C57Bl/6J background. E0.5 was established as the day of vaginal plug. All experimental procedures were performed at E15.5 and were approved by the Geneva Cantonal Veterinary Authority.

Coordinated pregnant mice of 12-30 weeks were ethically sacrificed to extract E15.5 embryo eyes. Thirty retinas were extracted in ice-cold L-15, micro-dissected under a stereomicroscope and incubated in 200 µl single-cell dissociation solution consisting of papain (1 mg/ml)-enriched HBSS at 37°C for 12 min with trituration every 2 min. At the end of the incubation time, cells were further dissociated via gentle up-and-down pipetting. The reaction was stopped with the addition of 400 µl of 2 mg/ml ovalbumin-enriched cold HBSS and the cell suspension was then passed through a 40-µm cell strainer to remove cellular aggregates. Cells were then centrifuged for 5 min at 500 g at 4°C. After removal of the liquid, the pellet was suspended in 250 µl of cold HBSS and the resulting solution was finally sorted on a MoFloAstrios device (Beckman) to reach a concentration of 410 cells/µl.

The libraries of single cells were prepared using the Chromium 3′ v2 platform following the manufacturer's protocol (10x Genomics). Briefly, single cells were partitioned into gel beads in EMulsion (GEMs) in the GemCode instrument followed by cell lysis and barcoded reverse transcription of RNA, amplification, shearing and 5′ adaptor and sample index attachment. On average, approximately 5000 single cells were loaded on each channel with 2675 cells recovered in the first replicate (index F2), and 2673 cells recovered in the second (index E2). Libraries were sequenced on a HiSeq 4000 (paired end reads: Read 1, 26 bp, Read 2, 98 bp).

Slc6a4-Cre;tdTomato cells were dissociated from E15.5 retina and were sorted on MoFloAstrios (8 μl) to obtain both positive (tdT⁺) and negative (tdT⁻) populations. Each condition was then mixed with the C1 suspension reagent (2 μl; Fluidigm) yielding a total of 10 μl of cell suspension mix with ∼500 cells/μl. The cell suspension mix was then loaded on a C1 Single-Cell AutoPrep integrated fluidic circuit (IFC) designed for 10- to 17-μm cells (HT-800, Fluidigm, 100-57-80). cDNA synthesis and preamplifications were processed following the manufacturer's instructions (C1 system, Fluidigm) and captured cells were imaged using the ImageXpress Micro Widefield High Content Screening System (Molecular Devices). Single-cell RNA-sequencing libraries of the cDNA were prepared using the Nextera XT DNA library prep kit (Illumina). Libraries were multiplexed and sequenced according to the manufacturer's recommendations with paired-end reads using HiSeq4000 platform (Illumina) with an expected depth of 0.5 M reads per single cell, and a final mapping read length of 98 bp. Captures and sequencing experiments were performed within the Genomics Core Facility of the University of Geneva. The sequenced reads were aligned to the mouse genome (GRCm38) using the read-mapping algorithm STAR version 2.6.0c. The number of raw counts obtained were reads per million mapped reads normalized (RPM).

Both replicates (2675 and 2673 cells) were preliminarily analyzed separately using the Chromium v2 analysis software Cell Ranger (Version 2.1.1) giving similar t-SNE-based clustering (Fig. S1D,E). Seurat package version 2.3 (Butler et al., 2018) was used to import both datasets in R version 3.4.4 (R Core Team). Cells considered during the creation of the Seurat objects were expressing at least 200 genes, and genes kept are expressed in a minimum of three cells. Mitochondrial gene effect was regressed out for the whole dataset. Merging of replicates was performed by aggregation with the in-built Seurat function MergeSeurat and batch-effect correction was made possible by linear regression of the transcriptomic expression of the two batches during the scaling and centering of the dataset by the ScaleData function. 1648 variable genes were defined on a variability plot by the FindVariableGenes Seurat function as an RPM mean expression above 0.1 and dispersion above 1.

Seurat package version 2.3 (Butler et al., 2018) was used to import in R version 3.4.4 (R Core Team) 800 cells. 719 cells considered during the creation of the Seurat objects expressed at least 200 genes, and genes kept are expressed in a minimum of three cells. Mitochondrial gene effect was regressed out for the whole data. 2716 variable genes were selected by an RPM mean expression above 0.01 and a dispersion above 0.5.

PCA was performed on variable genes to reduce dimensionality of the dataset. Spectral t-SNE was based on the reduced dimensional space of the five most significant dimensions of the PCA using the Rtsne package Barnes-Hut version of t-SNE forked in Seurat with a perplexity set at 30. Dimensions used 1:5 for 10x and 1:8 for C1, default parameters (van der Maaten, 2014). A t-SNE-based clustering analysis was then performed by the shared nearest neighbor (SNN) modularity optimization algorithm (van der Maaten, 2014). Differentially expressed genes between clusters were obtained by Seurat-implemented Wilcoxon rank sum tests with default parameters. The cluster identities in this t-SNE-space were uncovered by feature plots of typical cell type marker genes. Another approach of the dimensionality reduction problem was explored by UMAP (McInnes et al., 2018), generated with the help of Seurat and the UMAP-learn python package on the ten most significant dimensions of the PCA. Three embedded dimensions of the UMAP were outputted for further use to represent the pattern of various features during differentiation. The minimal distance parameter of the UMAP was set to 0.3 and the number of neighboring points used in local approximations was set to 30. GO term analyses were performed using DAVID bioinformatics resources 6.8. (Huang et al., 2008, 2009).

Pseudo-time analysis was performed with Monocle 2 using genes that have passed the quality control of the Seurat object creation (Qiu et al., 2017a). Genes considered as defining the progression of the pseudo-time were those that were detected as having an expression above 0.5 by Monocle 2. Negative binomial was considered for the model encoding the distribution that describes all genes. During the pseudo-time processing, the dimensionality of the dataset was reduced by the Discriminative Dimensionality Reduction with Trees (DDRTree) algorithm on the log-normalized dataset with ten dimensions considered. Branched expression analysis modeling was performed on major branching points with default parameters.

Analysis of RGC and PR clusters (merging, normalization, batch correction, dimensionality reduction techniques and differential expression) was carried out as previously described using Seurat with the exception of variable genes defined from the variability plot as genes with mean expression above 0 and dispersion above 0.8. Transcriptional wave analysis was performed with default parameters as described previously (Telley et al., 2016). Genes presenting interesting variations were regrouped along a pseudo-time axis, forming clusters composed of similar time-dependent gene expressions. Theses clusters of patterns were labeled as waves one to six (RGCs, Fig. 4) or one to four (PRs, Fig. S8, Table S7).

Cell cycle analysis was performed by isolating the cells of the proliferating progenitor clusters and analyzing them independently. Briefly, we considered for analysis the 2000 highly variable genes selected using a model to the CV-mean relation feature selection procedure previously described (La Manno et al., 2016); dimensionality of the dataset was then reduced to the first 32 principal components using PCA. We used the principal components as input to compute a two-dimensional embedding using the UMAP algorithm (McInnes et al., 2018) and clustered the cells using the Louvain community detection algorithm (Blondel et al., 2008). The velocyto pipeline (La Manno et al., 2018) was used to compute the RNA velocity field on the manifold obtained by UMAP. The parameters used were k=120 for knn imputation and n_neighbours=300 to estimate the transition probability and the scaling of the vectors adjusted using a randomized control.

We fitted a pseudo-time model of the cell cycle using the non-parametric principal curve algorithm described by Ozertem and Erdogmus (2011). The pseudo-time coordinate correspondent to each cell was computed by summing the distance between consecutive cells starting from an arbitrary point and progressing through a full period. For every gene expression, we fitted a smooth function by support vector regression using the pseudo-time coordinate as the only predictor (scikit-learn implementation was used and parameters were set to kernel=”rbf” gamma=0.2, C=0.5 and periodic boundary conditions imposed).

We devised an analysis strategy to identify initial cell-type commitment within the cells that were proliferating. This strategy is based on two ideas: (1) that the genes involved in commitment of progenitor subsets will show a periodical pattern and (2) that a stronger difference between these progenitors should be detectable by focusing on the unspliced fraction of the RNA counts. Therefore, we based our analysis on unspliced RNA level expression, using the counts obtained by velocyto. We considered genes that were among the top 1400 variable genes of either the spliced or the unspliced count matrixes (2311 genes in total). Then, we excluded genes for which the pseudo-time fit had a coefficient of determination higher than 0.7. The resulting set of genes was used to compute a UMAP embedding and to cluster the cell using the Louvain community detection algorithm. Clusters were then used to calculate a gene enrichment score (μ^_cluster+0.1)/(μ^_all+0.1)* (f^_cluster+0.1)/(f^_all+0.1) where μ^ is the sample average for the gene and f ^ is the fraction of cells expressing the gene. Genes with high enrichment score were inspected and a selection was presented in Fig. 3 and Fig. S7.

Computation of RNA velocities uses the abundances of spliced and unspliced RNA to fit a model of gene expression and to estimate the rate of change of expressions in time across the whole transcriptome. This model is then used to extrapolate the short future gene expressions of each cell of the dataset. We used velocyto CLI to obtain separated counts for spliced and unspliced molecules and we used velocyto analysis module to calculate and visualize velocities following the same filtering and pre-processing procedure as done in the original study (La Manno et al., 2018). We processed the reads of the single-cell experiment with the run10x function on the loom file outputted from the CellRanger Pipeline, with reference genome of the mouse, mm10 (Ensembl 84), from 10x Genomics, and masked the corresponding expressed repetitive elements that could constitute a cofounding factor in the downstream analysis with the UCSC genome browser mm10_rmsk.gtf. Downstream analysis and representation of velocities were performed by the velocyto python package for main figures (La Manno et al., 2018). We used the Jupyter notebook DentateGyrus.ipynb as a guideline for the latter with default parameters for fitting gene models, normalizing and representing the data with following exceptions: the 15 first principal components of the PCA were used in analysis and the matrix were smoothed by k-nearest neighbors (knn) using 180 neighbors. RNA velocity transcriptional dynamic was embedded on a UMAP space computed with the UMAP python package, parameters set for a number of neighbors of 120, a learning rate of 0.5 and a min distance of 0.4.

To establish spatial positional information on single cells, we plotted the expression levels of known marker genes of the four different orientations on RGCs (Table S4) (Behesti et al., 2006, 4; McLaughlin and O'Leary, 2005; Takahashi, 2003). The selection of dimensions 1 and 3 was based on the visualization of the first five t-SNE dimensions (Fig. S15, DVTN t-SNE). For each RGC cell, a DV and a TN score was calculated as the log2 ratio of the quantile 80 of the dorsal genes (or temporal genes), with the quantile 80 of the ventral genes (or nasal).

A pseudo-expression of 1 was used to avoid infinite values. Those scores were used to classify the RGC cells into four groups: ‘DT’ (DV>0 and TN>0), ‘VT’ (DV<0 and TN>0), ‘VN’ (DV<0 and TN<0) and ‘DN’ (DV>0 and TN<0). Only cells with a non-null DV or TN score were classified (606 cells). Differentially expressed genes between the cell groups were obtained by Seurat-implemented Wilcoxon rank sum tests with default parameters. These genes were validated by ISH with quantification using Fiji (Fig. S9). ISH images were converted to 8-bit and thresholded using the RenyiEntropy. A mask on a selection of the retinal area to exclude the retinal pigmented epithelium was generated using the ‘Analyse Particles’ function. Finally, the gradient distribution was shown as a ‘Plot Profile’.

Co-expression analysis was performed as previously described (Fabre et al., 2018). For C1 cells (Fig. 7J), we analyzed only the RGC cluster without the old-RGC (Fig. 7B, red). Thirty-seven genes showed a minimum absolute Spearman's correlation of 0.25 with the 358 cells that expressed at least one of the nine ipsilateral-RGC known genes. Genes expressed in less than 5% of those 358 cells were excluded from the analysis. The seven clusters shown on the left side of the heatmap were identified by complete-linkage clustering based on Euclidian distance. Similarly, for RGCs captured by 10x (Fig. S13), 176 genes were found with a minimum absolute Spearman's correlation of 0.2 with the 781 cells expressing at least one of the nine ipsilateral-RGC known genes. On this heatmap, eight clusters identified previously are shown.

Embryonic heads from E15.5 or E17.5 mouse embryos were fixed overnight at 4°C with 4% paraformaldehyde, rinsed with PBS and then cryopreserved through 30% sucrose and frozen in optimal cutting temperature (OCT; Sakura). Eyes were cryosectioned as 16μm-thick sections that were dried 1 h prior to immunostaining. Immunostaining was performed using rabbit anti-dsRed (1:1000; Clontech, 632496), rabbit anti-Sox2 (1:500; Abcam, ab97959), rabbit anti-BRN3B (Pou4f2) (1:200; Santa Cruz Biotechnology, N15; sc-31987), rabbit anti-cKIT (1:200; Santa Cruz Biotechnology, C-19; sc-168), rat anti-SST (1:400; Millipore, MAB354) and rat anti-TBR2 (1:500; Invitrogen, 14-4875-82) and images were acquired using an Eclipse 90i epifluorescence microscope (Nikon, full retina view) or confocal Zeiss LSM800 (high magnifications in airyscan mode). All ISH data was retrieved either from the Allen Developing Mouse Brain Atlas [Figs 2F, 4F, 5E (Reln), 7H (Sst), S3A-F, S7C, S8F (Rorb and Lhx4) and S9 (Cnr1, Cntn2 and Irx2); www.brain-map.org] or from the digital atlas of gene expression patterns in the mouse [Figs 5E (Kitl), 7H (Gal and Tac1), S7C (Slc1a3, Pak1 and Vax2os), S8F (Tgfb2) and S9 (Fstl4, Nr2f2 and Vax2os); https://gp3.mpg.de/].

We thank Nicolas Lonfat and members of the Jabaudon lab for critical reading of the manuscript, as well as Alexandre Dayer and Denis Jabaudon for their advice and for sharing mice and reagents. We also thank Audrey Benoit, Christelle Borel and Wafae Adouan for assistance in the 10x procedure. We thank Polina Oberst and Philipp Abe for their help with the CAG-RFP and Rbp4-Cre mouse line, and N. Baumann for the generation of the 3D UMAP animation. We thank the Geneva Genomics Platform (M. Docquier and D. Chollet, University of Geneva) and the Bioimaging Core Facility (F. Prodon and O. Brun, University of Geneva).