Mouse RGCs are born between approximately embryonic days (E) 11 and 17 with newborn RGCs exiting the mitotic cycle near the apical margin, then migrating basally to form the ganglion cell layer (Dräger, 1985; Marcucci et al., 2019; Voinescu et al., 2009; Figure 1b). Reported birthdates differ among publications and are complicated by naturally occurring cell death and the central-peripheral developmental gradient, but a detailed analysis concludes that >95% of RGCs in adult mouse retina are born after E12.8 and >85% before E16 (Farah and Easter, 2005). Shortly after they are born, RGCs extend axons through the optic nerve, with some reaching retinorecipient areas by E15 (Godement et al., 1984; Osterhout et al., 2011) and forming diverse projection patterns (Martersteck et al., 2017). During early postnatal life, they extend dendrites apically into the inner plexiform layer of the retina, receiving synapses from amacrine cells by postnatal day (P)4 and bipolar cells a few days later (Kim et al., 2010; Lefebvre et al., 2015). Light responses are detected in RGCs by P10 but image-forming vision does not begin until eye-opening, around P14 (Hooks and Chen, 2020).

To determine when and how RGCs diversify, we used droplet-based scRNA-seq (Macosko et al., 2015; Zheng et al., 2017) to profile them at five stages: E13 and E14 (during the period of peak RGC genesis), E16 (by which time RGCs axons are reaching target retinorecipient areas), P0 (as dendrite elaboration begins), and P5 (shortly after RGCs begin to receive synapses). As RGCs comprise ≤1% of retinal cells (Jeon et al., 1998), we enriched them with antibodies to two RGC-selective cell surface markers, Thy1/CD90 (Barres et al., 1988) and L1cam (Demyanenko and Maness, 2003; Figure 1—figure supplement 1a).

We obtained a total of 98,452 single-cell transcriptomes with acceptable quality metrics (Materials and methods). Of these, we identified 75,115 (76%) as RGCs based on their expression of canonical RGC markers, including Rbpms (an RNA-binding protein) and Slc17a6 (the vesicular glutamate transporter VGLUT2) (Figure 1c–e, Figure 1—figure supplement 1b and c). Non-RGCs included amacrine cells (Tfap2a+ Tfap2b+), cone photoreceptors (Otx2+ Crx+ ), microglia (P2ry12+ C1qa+), anterior segment cells (Mgp+ Bgn+ ), and RPCs. Anterior segment cells were found only in E13 and E14 samples because whole eyes were dissociated at these stages. RPCs formed a continuum, containing both ‘primary’ RPCs expressing cell cycle-related genes (e.g., Mki67, Ccnd5, and Birc5) and previously described RPC regulators (e.g., Sfrp2, Vsx2, and Fgf15), and ‘neurogenic’ RPCs expressing proneural transcription factors (TFs) (e.g., Hes6, Ascl1, and Neurog2) (Clark et al., 2019). Importantly, these markers were not expressed in cells annotated as RGCs (Figure 1—figure supplement 1d). These stringent criteria ensured that our dataset comprised postmitotic committed RGCs, allowing us to focus on their diversification and maturation.

Overall, we recovered ~5900–18,500 RGCs at each of the five time points. Of the two surface markers used for enriching RGCs, Thy1 was effective at later stages as shown previously (Kay et al., 2011; Rheaume et al., 2018; Tran et al., 2019), whereas L1cam expression was more selective at E13 and E14 (Figure 1—figure supplement 1b and c). However, identical clusters were observed with both methods at E13, E14, and E16, albeit with different frequencies. This concordance supports the idea that neither marker fails to capture particular RGC types. To further evaluate the effectiveness of our enrichment strategy at early stages, we compared our data with two recent studies in which developing retinal cells were profiled using scRNA-seq without any enrichment (Clark et al., 2019; Lo Giudice et al., 2019). A joint analysis of these datasets at embryonic time points showed consistency in the transcriptional signatures of major cell groups without discernible biases (Figure 1—figure supplement 1e–g). However, our enrichment protocols increased the fractional yield of RGCs by >3× at E14 and E16 and by >100× at P0 (Figure 1—figure supplement 1h), which enabled us to resolve heterogeneity within this class at immature stages. We also compared our P5 data with those from an earlier study in which P5 RGCs were profiled (Rheaume et al., 2018) and found a good correspondence (Figure 1—figure supplement 1i). For the analysis that follows, we combined precursor RGCs (E13–P5) with a previously described dataset of 35,699 mature RGCs at P56 (Tran et al., 2019).

One can envision two extreme models of RGC diversification. In one, RGC type would be specified at or before mitotic exit, with each type arising from a distinct set of committed precursors that are transcriptomically defined. At the other extreme, all precursor RGCs would be identical when they exit mitosis and gradually acquire distinct identities as they mature (Figure 2a). Intermediate models could involve multiple groups of precursor RGCs, each biased towards a distinct set of terminal types.

Figure 2 with 1 supplement see all

Download asset Open asset

To distinguish among these alternatives, we analyzed the transcriptomic diversity of RGCs at each developmental stage using the same dimensionality reduction and graph clustering approaches devised for analysis of adult RGCs (Tran et al., 2019; see Materials and methods). This analysis led to three main results.

First, RGCs were already heterogeneous soon after mitotic exit. There were 10 transcriptionally defined precursor clusters at E13 (Figure 2b) before or at the peak time of RGC birth. The number of discrete clusters increased only slightly by E14 (from 10 to 12; Figure 2c), arguing against a model in which the number of precursor types extrapolated back to 1. No single cluster dominated the frequency distribution at either time as would be expected if a totipotent precursor RGC were to exist shortly after terminal mitosis.

Second, the number of transcriptionally defined clusters increased gradually, between E13 and adulthood, reaching 45 only after P5 (Figure 2b–g). Several arguments indicate that this increase is biologically significant rather than being an artifact of the data or computational analysis. (1) We used the same clustering procedure at all ages. (2) The qualitative trends were robust against variations in clustering parameters. (3) All embryonic clusters contained cells isolated with both cell markers, L1cam and Thy1 (Figure 2—figure supplement 1a–c), indicating that lower cluster numbers at early stages did not result from biased collection methods. (4) The increase in the number of effective molecular types was robust as demonstrated by three diversity indices – Rao, Simpson, and Shannon – all of which buffer against artificial inflation of diversity due to small clusters (Figure 2h, Figure 2—figure supplement 1d; see Materials and methods). (5) There was no systematic dependence of the number of clusters on the number of cells. For example, we identified 12 clusters from 17,100 cells at E14 and 38 clusters from 17,386 cells at P5.

Third, the transcriptomic variation became increasingly discrete with age. We quantified this increase in inter-cluster separation by calculating (1) the average cross-validation error of a multi-class classifier, and (2) the ratio of mean cluster diameter to mean inter-cluster distance in the low dimensional embedding (Materials and methods). Both metrics decrease in numerical value as the clusters are more well defined. From these trends, we conclude that the boundaries between RGC clusters become sharper as development proceeds (Figure 2i and j).

Taken together, our results show that transcriptomic clusters of RGCs increase in number and distinctiveness over time, making it unlikely that RGC-type identity is fully specified at the progenitor stage.

We next investigated the temporal relationships among precursor RGC clusters identified at different ages. We again consider two extreme models. In a ‘specified’ model, each terminal type arises from a single cluster at every preceding developmental stage (Figure 3a, left). In this model, distinct transcriptomic states among precursor RGCs correspond to distinct groups of fates. At the other extreme, distinct clusters would share similar sets of fates (Figure 3a, right). In an intermediate model, fates of precursor clusters would exhibit partial overlap.

Figure 3 with 1 supplement see all

Download asset Open asset

As a first step in discriminating among these scenarios, we used transcriptome-wide correspondence among clusters as a proxy for fate association. We identified mappings between clusters across each pair of consecutive developmental stages (E13–E14, E14–E16, E16–P0, P0–P5, and P5–P56) using gradient boosted trees (Chen and Guestrin, 2016), a supervised classification approach (Materials and methods). In each case, a classifier trained on transcriptional clusters at the older stage was used to assign older cluster labels to cells at the younger stage (e.g., E16 labels assigned to E14 RGCs). Patterns expected for the extreme models are schematized as ‘confusion matrices’(Stehman, 1997) in the lower panels of Figure 3a.

Correspondence fell between the two extremes (Figure 3d–h, Figure 3—figure supplement 1a–d). We quantified the extent of correspondence using two metrics: normalized conditional entropy (NCE) and the adjusted Rand Index (ARI) (Materials and methods). Both NCE and ARI are restricted to the range (0,1), with lower values of NCE and higher values of ARI consistent with a specified mode of diversification. Both metrics exhibited an increased degree of specificity with age (Figure 3b and c). Since NCE and ARI provide a single measure of specificity for the entire datasets being compared, we also computed a ‘local metric,’ the occupancy fraction (OF), which quantifies mapping specificity for each cluster (Materials and methods). Results based on this metric were consistent with increased specificity of correspondence with age (Figure 3—figure supplement 1e). Overall, this analysis of transcriptomic correspondence suggests that poorly specified relationships among transcriptomic clusters at early stages are gradually refined to yield increasingly specific associations at later stages.

The classification analysis presented so far relied on comparing clusters between ages and was therefore unable to link individual precursors to specific terminal fates. At one extreme, individual precursor clusters might contain several groups of cells, each committed to a distinct, small number of fates. Alternatively, individual cells might be as multipotential as the clusters in which they reside (Figure 4a).

Figure 4 with 2 supplements see all

Download asset Open asset

Unfortunately, this approach does not afford a straightforward way to explore variations in patterns of fate associations within clusters. We therefore turned to Waddington optimal transport (WOT), a computational method rooted in optimal transport theory (Kantorovich, 2006; Monge, 1781) that utilizes scRNA-seq measurements at multiple stages to infer developmental relationships (Schiebinger et al., 2019). Briefly, WOT computes a ‘transport matrix’ $\mathrm{\Pi }$ between each pair of consecutive ages with elements ${\mathrm{\Pi }}_{ij}$ encoding fate associations between a single RGC i at the younger age and RGC j at the older age (see Materials and methods). WOT directly computes fate associations at the level of individual cells without requiring clustering as a prior step. We conducted computational tests to assess the numerical stability of associations reported by WOT (Figure 4—figure supplement 1). We also determined that when collapsed to the level of clusters the WOT-inferred transport maps strikingly mirrored the confusion matrices obtained from multi-class classification (Figure 4—figure supplement 2).

Based on the success of these tests, we applied WOT to compute the ‘terminal fate’ for each precursor RGC. We leveraged the fact that in WOT fate associations between RGCs at nonconsecutive ages (e.g., E16 and P56) can be estimated in a principled way by multiplying the intermediate transport matrices. This yielded a fate vector $\overrightarrow{f}$ for each of the 75,115 immature RGCs, whose kth element $f_k$ represents the predicted probability of commitment to adult type $k\in (1,2,\dots ,45)$ (Materials and methods). A fully committed precursor would have all but one element of $\overrightarrow{f}$ equal to zero, whereas a partially committed precursor would have multiple nonzero elements in $\overrightarrow{f}$ . Since the elements of $\overrightarrow{f}$ are interpreted as probabilities, they are normalized such that ${\sum }_{\mathrm{k}}{f}_k=1$.

We quantified the commitment of each precursor by computing its ‘potential’ $P=\frac{1}{{\sum }_{\mathrm{k}}{f}_k^2}$ , which is defined analogously to the ‘inverse participation ratio’ in physics (Fyodorov and Mirlin, 1992). In our case, the value of P for a given RGC ranges continuously between 1 and 45, with lower values implying a commitment to specific fates and higher values reflecting indeterminacy. Importantly, this measure of commitment does not rely on arbitrary thresholding of the $f_k$ values to assign precursors to types.

Five results emerged from this analysis.

• Nearly all prenatal RGCs (i.e., on or before P0) were multipotential rather than committed to a single terminal fate, with individual potentials distributed across a range of values (Figure 4b).

• Multipotentiality was a general feature of immature RGCs, being present in cells of all clusters at E13, E14, and E16 (Figure 4c–f).

• At early stages, the average value of P varied among transcriptomic clusters, reflecting asynchronous specification (Figure 4c). The tempo of commitment is further explored in the next section.

• Although they were multipotential, no precursor RGC was totipotential (i.e., completely unspecified, corresponding to P = 45). At E13, the average value of P was 11.6 ± 4.9, which was fourfold lower than the maximum possible value of 45, and no precursor had P > 30.

• Finally, inferred multipotentiality decreased gradually during development, and some persisted postnatally (average P = 3.4 ± 2.1 at P0, and 1.6 ± 0.9 at P5; Figure 4g and h).

From these results, we conclude that early postmitotic RGCs are multipotential but not totipotential, and that type identity is specified gradually via progressive restriction.

As a first step in understanding the progressive restriction of RGC fate, we analyzed the extent to which pairs of mature types were likely to have arisen from the same set of immature precursors. To this end, we computed a ‘fate coupling’ value $C(l,m;age)$ for each pair of terminal RGC types (l and m), defined as the Pearson correlation coefficient between the values of $f_l$ and $f_m$ across all precursors at a given age (Materials and methods). $f_l$ and $f_m$ are fate probabilities corresponding to types l and m as defined in the previous section. Values of $C(l,m;age)$ in our data ranged from –0.11 to 0.95. Higher values of $C(l,m;age)$ indicate strong fate coupling between types l and m, implying the existence of common postmitotic precursors, whereas low $C\left(l,m;age\right)$ values suggest that types l and m arose from largely nonoverlapping sets of precursors. We visualized the pattern of fate couplings as network graphs, where the nodes represent types and the edge weights represent values of $C(l,m;age)$. The arrangement of nodes was determined at E13 using a force-directed layout algorithm (Fruchterman and Reingold, 1991), with pairwise distances being inversely proportional to the values of $C(l,m;E13)$, the fate coupling values at E13 (Figure 5a). To visualize the temporal evolution of these fate couplings, we retained the same layout of nodes while updating edge weights according to $C(l,m;age)$ (Figure 5b–e).

Figure 5 with 1 supplement see all

Download asset Open asset

Types that were coupled in fate at the earliest time point gradually decoupled as development proceed. For example, at E13, 118/990 pairs (12%) were strongly coupled (threshold of $C(l,m;age)$>0.2 as determined by randomization tests; see Materials and methods), while at P5, only 8/990 (<1%) passed this criterion (Figure 5a and e). Lowering this threshold for coupling to 0.05 increased the number of strongly coupled pairs at P5 to only 2% (20/990).

Different pairs of types decoupled at different rates (Figure 5f). As they decoupled, RGC precursors became increasingly restricted to a single type (i.e., $f_k\gg f_{l\ne k}$ for a precursor favoring type k). This corresponded to a ‘localization’ of precursors in transcriptomic space and is a proxy for specification (see Materials and methods). We modeled the extent of localization vs. age via a logistic function (Figure 5g, Figure 5—figure supplement 1d) and used this to calculate a specification time for each type (${\tau }_{sp}$) (see Materials and methods for details). Based on this analysis, 7/45 types are specified postnatally. The average ${\tau }_{sp}$ for RGCs was E17.8, but individual RGC types exhibited a wide range from E13.9 to P5.2 (Figure 5h). The inferred specification time was not correlated with adult frequency (Figure 5i).

We illustrate this range by considering three pairs of RGC types in Figure 5—figure supplement 1. C12 and C22 (numbered as in Tran et al., 2019; see Figure 2g) exhibit low fate coupling at all ages profiled (Figure 5—figure supplement 1a), indicative of separate precursor populations. In contrast, C19 and C20 decouple only at P0, implying the existence of a common precursor throughout embryogenesis (Figure 5—figure supplement 1b). C21 and C34 display an intermediate pattern, decoupling around E16 (Figure 5—figure supplement 1c). Taken together, these results suggest that RGC types emerge by asynchronous fate decoupling of multipotential precursors.

Because fate coupling is a metric of inferred overlap of developmental history, it is likely that tightly coupled types share common precursors. This relationship implies that tightly coupled types might also be specified by common transcriptional programs. As a step towards identifying candidate fate determinants, we identified eight TFs that are expressed by distinct groups of mature RGC types (Figure 6a, Figure 6—figure supplement 1a). Three of these are well-characterized RGC-selective TFs: Foxp2, expressed by five F-RGC types (Rousso et al., 2016); Tbr1, expressed by five T-RGC types (Liu et al., 2018); and Eomes (also known as Tbr2), expressed by seven types (Mao et al., 2020; Tran et al., 2019). The seven Eomes/Tbr2 types include the melanopsin-expressing intrinsically photosensitive (ip) RGC types (Berson et al., 2002). The remaining five were Neurod2, Irx3, Mafb, Tfap2d, and Bnc2, which label 8, 5, 4, 6, and 3 types, respectively. Eomes types also co-expressed Tbx20 and Dmrbt1 while Neurod2 types also co-express Satb2. Together, 40/45 mature types expressed at least one of these TFs in a manner that was, with a few exceptions, mutually exclusive. In many cases, the fate proximity of types that shared TF expression was obvious (Figure 6a).

Figure 6 with 1 supplement see all

Download asset Open asset

We refer to these TF-based groups as fate-restricted RGC subclasses – an intermediate taxonomic level between class and type based on inferred fate relationships. Consistent with their definition, the pairwise fate coupling among types from different subclasses was significantly lower than among types from the same subclass (Figure 6b). Thus, precursor RGC states associated with any two subclasses are more distinct than those associated with any two types. This is evident by the significant separation at E13 and negligible overlap at P5 for precursors favoring the Eomes, Mafb, and Neurod2 subclasses, respectively, as shown in Figure 6c–e.

We also asked whether the TF-based subclasses differed in inferred transcriptomic specification time ${\tau }_{sp}$ , as defined in Figure 5g. As shown in Figure 6f–h and Figure 6—figure supplement 1b–e, four subclasses were specified within a narrow interval (E16.8–E17.2), but three others differed substantially. The average specification time for the Eomes group was E14.6 (p<0.0001, Student’s t-test, compared to the mean for all types), while that for the Mafb and Neurod2 groups were E16.9 (p<0.001) and E18.5 (p<0.0001), respectively. The early specification of the Eomes group is consistent with birthdating studies showing the average earlier birthdate of ipRGCs compared to all RGCs (McNeill et al., 2011).

In summary, our results suggest the existence of fate-restricted RGC subclasses that arise from distinct sets of precursors and diversify into individual types. This method of defining RGC groups, which relies on inferred proximity of precursors in transcriptomic space, is distinct from previous definitions of RGC subclass based on shared patterns of adult morphology, physiology, or gene expression (see Discussion). Accordingly, the fate couplings at E13 were only weakly correlated with transcriptomic proximity in the adult retina (Figure 6—figure supplement 1f). Further, while TF-based groups align with some previously defined subclasses (e.g., ipRGCs or Tbr1+ RGCs), they do not map to other subclasses such as alpha-RGCs (four types) or T5-RGCs (nine types) (Figure 6—figure supplement 1g).

Finally, we considered the origin of two RGC groups defined by their projections: those with axons that remain ipsilateral at the optic chiasm (I-RGCs) and those that cross the midline to innervate contralateral brain structures (C-RGCs). The proportion of I-RGCs varies among vertebrates, in rough correspondence to the extent of binocular vision, ranging from none in most lower vertebrates to ~50% in primates. In mice, 3–5% of RGC axons remain ipsilateral, with most I-RGCs residing in the ventrotemporal (VT) retinal crescent (Mason and Slavi, 2020). While some I-RGCs have been observed to project from the dorsocentral retina during embryonic stages, these are rapidly eliminated so-called ‘transient’ I-RGCs (Soares and Mason, 2015). Thus, in adulthood, C-RGCs are present throughout the retina while ‘permanent’ I-RGCs are confined to the VT crescent.

The zinc-finger TF Zic2 is expressed in a subset of postmitotic RGCs in VT retina and is both necessary and sufficient for establishing their ipsilateral identity (Herrera et al., 2003); transient dorsolateral I-RGCs do not express Zic2 (Pak et al., 2004). Isl2 marks a subset of C-RGCs throughout the retina and appears to specify a contralateral identity in part by repressing Zic2 (Pak et al., 2004). These two TFs were expressed in a mutually exclusive fashion in RGC precursors between E13 and E16 (Figure 7a); Zic2 was downregulated at later ages (Figure 7—figure supplement 1a). Furthermore, Zic2 expression at E13 correlated with Igfbp5 and Zic1, and anticorrelated with Igf1 and Fgf12, consistent with recent reports (Wang et al., 2016; Figure 7b, Figure 7—figure supplement 1b and c). We scored each cell at E13 based on its expression of ipsilateral genes (Materials and methods), confirming that the expression of ipsilateral and contralateral gene signatures was anticorrelated (Figure 7c). Together, these observations support the idea that at E13 Zic2+ cells represent I-RGCs and Isl2+ cells represent some but not all C-RGCs.

Figure 7 with 1 supplement see all

Download asset Open asset

Using WOT, we then identified the descendants of presumptive I-RGCs at later ages. We found that I-RGCs comprised 4.3% of adult RGCs, consistent with the range of 3–5% noted above. We queried these cells to identify the genes that distinguished putative I-RGCs and C-RGCs throughout the developmental time course. At a fold change of ≥1.5, we found 59 differentially expressed (DE) genes at E13 and 89 at E14 (Figure 7e and f). In addition to Zic2, Igfbp5, Isl2, and Igf1, which had been used to define I-RGCs and C-RGCs at E13, they included Igfbpl1, Pou3f1, and Cntn2 enriched in I-RGCs, and Lmo2, Pcsk1n, and Syt4 enriched in C-RGCs. The number of genes differentially expressed between I- and C-RGCs decreased after E14, with 20, 9, and 0 significant genes at E16, P0, and P5, respectively (Figure 7—figure supplement 1d and e), presumably reflecting the downregulation of axon guidance programs once retinorecipient targets have been reached (see Discussion).

We also asked which RGC types included I-RGCs. At E13, putative I-RGCs were highly enriched in 2 of 10 clusters, comprising 38–40% of clusters 2 and 9, 9–14% of clusters 3 and 5, and <2% of the other six clusters (Figure 7d). In adults, RGCs expressing Tbr1, Mafb, Foxp2, and Neurod2 contained 3–4× more I-RGCs than RGCs expressing Eomes, Irx3, or Tfap2d. These results are consistent with previous observations that I-RGCs are morphologically and physiologically heterogenous but not uniformly distributed across types (Hong et al., 2011; Johnson et al., 2021). Lastly, the WOT-predicted relationship between E13 precursor RGC clusters and I-RGC-rich or -poor terminal types was consistent with these patterns. The top six I-RGC-rich types (C4, C15, C19, C20, C38, and C45) derived 14 and 4% of their relative fate association from E13 clusters 2 and 9, while the top six I-RGC-poor types (C8, C14, C18, C22, C31, and C41) derived only 3.8 and 0.2% of their relative fate association from E13 clusters 2 and 9. Thus, E13 clusters 2 and 9 are preferred precursors of adult types that are relatively rich in I-RGCs.

Definitions of neuronal class, subclass, and type have been contentious (Yuste et al., 2020). In general, classes share general features of structure, function, molecular architecture, and location, whereas types comprise the smallest groups within classes that can be qualitatively distinguished from other groups based on these and other criteria. Subclasses lie in-between. For RGCs, class identity has been clear for a century, but inventories of subclasses and types have emerged only over the last few decades as high-throughput methods have been implemented for quantifying structural (primarily dendritic morphology), functional (responses to an array of visual stimuli), and molecular properties (gene and transgene expression) of large numbers of RGCs. Fortunately, to the extent that they have been compared, there is excellent concordance among types defined by molecular, structural, and physiological criteria (Bae et al., 2018; Goetz et al., 2021; Tran et al., 2019; see http://www.rgc.types.org). Moreover, RGCs of a single type exhibit a regular spacing, called a mosaic arrangement, in that they tend to avoid other members of the same cell type, whereas their association with members of other types is random (Kay et al., 2012; Keeley et al., 2020; Rockhill et al., 2000). The molecular basis of this property is poorly understood, but it provides an additional criterion for defining a type. Thus, while no two RGCs are identical, and variation may be continuous in some other structures (Cembrowski and Spruston, 2019), there is strong reason to believe that RGC types are discrete.

Developmental trajectories of cell types cannot be better than the adult types at which they are aimed. We have two reasons to believe that our adult RGC atlas (Tran et al., 2019) is accurate and complete.

First, the atlas is based on a detailed analysis of 35,699 cells and is therefore powered to detect types occurring at ~0.1% frequency (>40 cells per type; https://satijalab.org/howmanycells/). Results were stable over a variety of parameters (Tran et al., 2019). Moreover, in the course of studies on responses of RGCs to injury, we recently profiled an additional ~120,000 cells (A. Jacobi, N. Tran, W Yan, and J.R.S, in preparation), without identifying additional types.

Second, RGCs have now been classified by functional and structural properties, based on physiological responses to visual stimuli (Baden et al., 2016; Goetz et al., 2021) and serial section electron microscopy (Bae et al., 2018). The numbers of types defined in these ways (47 in Bae et al., 2018, 42 in Goetz et al., 2021, and >32 in Baden et al., 2016) match well to the 45–46 defined molecularly (Tran et al., 2019).

The multipotentiality of dividing progenitor cells can be demonstrated by indelibly labeling a progenitor and then examining its progeny at a later stage. For mammals, this was initially done by infecting single cells with a recombinant retrovirus encoding a reporter gene that could be detected following multiple cell divisions (Price et al., 1987; Sanes et al., 1986; Turner and Cepko, 1987). More recently, it has become possible to greatly increase throughput by tracking scars or barcodes introduced by CRISPR/Cas9 (Baron and van Oudenaarden, 2019; Espinosa-Medina et al., 2019; McKenna et al., 2016). In sharp contrast, conclusively demonstrating that a single postmitotic cell is multipotential would require following a cell from an unspecified to a specified state, then turning back time, watching it again, and asking if it acquired the same mature identity. Since this is impossible, we used computational methods to draw tentative conclusions about the extent to which newly postmitotic RGCs are committed to mature into a particular type based on their transcriptomic profiles. Consequently, our results are based on inferred rather than experimentally determined lineages.

Our analysis proceeded in three steps. First, to ask whether RGCs were committed to a particular fate before or shortly after they were born, we assessed transcriptomic heterogeneity at a time when a large fraction was newly postmitotic (E13 and E14). We found that heterogeneity was present but limited: 10 transcriptomic clusters were distinguishable at E13 and 12 at E14. Thus, some heterogeneity is present in precursor RGCs, but far less than would be required to specify type identity before or immediately after their birth. A relevant issue is whether at these early stages the transcriptomic variation among cells could reflect variation in their stages of differentiation, perhaps as a consequence of different intervals between birthdate and sampling. Although we cannot exclude this possibility, inter-cluster variability of key RGC-specific genes (e.g., Rbpms) was no greater at early times than in adults (data not shown).

Second, we used a supervised classification approach to ask whether precursor RGC clusters mature into mutually exclusive sets of adult types. This model would imply an orderly, step-wise restriction of cell fates. However, our results indicate substantial overlap in the types derived from cells in different immature clusters. This result argues against a deterministic model of diversification and suggests that precursor RGCs are incompletely committed to a specific type for a substantial period after they are generated.

Third, we used optimal transport inference (WOT) to ask whether the multipotentiality observed at the level of groups was also a property of individual cells. This approach circumvents the limitation of the supervised classification approach, which compares similarity only at the level of clusters. WOT utilizes time-course scRNA-seq snapshots to infer fate associations between individual cells sampled at different time points, without reference to the clusters in which they reside (Schiebinger et al., 2019). While being consistent with supervised classification results at the cluster level, WOT indicated that the majority of individual RGCs were multipotential at E13 and E14. Of equal importance, immature RGCs were not totipotential: the average predicted potential (P) was 11.6 at E13, or ~25% of the maximum possible value of 45, and no RGCs had P > 30. We conclude that single multipotential immature RGCs are biased in favor of particular groups of adult RGC types.

Further analysis provided insight into the structure of multipotentiality among RGCs. The adult RGC types associated with a precursor RGC were not a randomly chosen subset; rather some were more likely to arise from a common precursor state (‘fate coupled’) than others. This suggests a model in which RGC types arise via a progressive decoupling of fates within multipotential precursors. Decoupling is asynchronously, with different types emerging at different times. By modeling the temporal kinetics of fate decoupling, we were able to estimate a tentative specification time for each type – that is, the time at which precursors acquire a distinct transcriptomic identity. Analysis of transcriptomic changes that occur during this process, and the effects of visual experience on maturation, will be presented elsewhere (K.S., I.E.W., S.B. and J.R.S., in preparation).

Our conclusions about fate restriction are based on analyzing cell states defined by expression patterns of highly variable genes (HVGs) identified in the data (Figure 2h). An alternative and common view is that if a small set of genes is sufficient to define cell state, they should be the focus of analysis. We believe this is an incorrect argument based on confusion between genes that determine a cell type or state and those that define it. The only way a small number of genes, whatever their functional role, can exclusively define cell state is if they are expressed at very high levels. If this were the case, they would be detected in our data and drive the clustering.

For RGCs, class identity has been clear for a century, and type identity has been solidified during over the last few decades, but criteria for defining subclasses remain unclear. Tentative classifications have used molecular, physiological, and morphological criteria (Sanes and Masland, 2015; Tran et al., 2019). In general, these criteria correlate imperfectly with each other, a main exception being that ON and OFF RGCs (responding preferentially to increases and decreases in illumination, respectively) have dendrites that arborize in the inner and outer portions of the inner plexiform layer (Famiglietti and Kolb, 1976).

The pattern of fate couplings between RGC types at E13–14 provides an alternative way to define RGC subclasses – groups of RGC types that arise from the restriction of a common transcriptionally defined precursor state. We identified TFs selectively expressed within these subclasses. Our rationale was that among them would be fate determinants, an idea that could be tested by conventional genetic manipulations. Support for this idea is that there is already strong evidence that one such factor is a fate determinant in mouse: Eomes is selectively expressed by ipRGCs (and a few other types), and Eomes mutants fail to form ipRGCs although their retinas are normal in many respects (Mao et al., 2014). This encourages the hope that some of the other TFs in this set are also fate determinants. It will also be interesting to determine whether members of fate-restricted subclasses share structural or functional properties.

The TFs Isl2 and Zic2 are selective markers of embryonic RGCs that project contralaterally or ipsilaterally, respectively, and are critical determinants of this choice (Herrera et al., 2003; Pak et al., 2004). We found that their expression was largely nonoverlapping in RGCs at E13, and that they were co-expressed with previously reported markers of contralaterally and ipsilaterally projecting RGCs, respectively. Because few RGC axons reach the optic chiasm before E14, our results are consistent with genetic evidence that this differential expression is a cause rather than a consequence of the divergent choices the axons make at the chiasm. Among many genes co-expressed with Isl2 or Zic2 may be others that play roles in this choice.

Zic2 is downregulated later in embryogenesis, so we selected some RGCs as putative I-RGCs using genes known to be expressed in them (Wang et al., 2016), then used WOT to infer the RGC types to which they gave rise. Our analysis suggests that I-RGCs comprise many differentiated types, consistent with previous results (Hong et al., 2011; Johnson et al., 2021). Surprisingly, however, there were few if any genes differentially expressed between the putative mature I- and C-RGCs. Assuming that WOT results are valid – an assertion that can be tested directly in the future – this result suggests a model in which newborn RGCs are doubly specified – by laterality and type – but that once axonal choice has been made the laterality program is shut down.

Generation of neuronal classes has been analyzed in many parts of the vertebrate nervous system, but we are aware of few reports on how classes diversify into types. A recent study addressed this issue for primary sensory neurons and reached the conclusion that newborn neurons in dorsal root ganglia are transcriptionally unspecialized and become type-restricted as development proceeds (Sharma et al., 2020). Similarly, both excitatory neuronal subclasses appear to diversify postmitotically in the mouse cerebral cortex (Di Bella et al., 2021; Lodato and Arlotta, 2015), and there is suggestive evidence that the same is true for interneuronal subclasses (Wamsley and Fishell, 2017). In all of these cases, it is attractive to speculate that diversification may occur by a process of fate decoupling in subpopulations of distinct multipotential precursors, akin to that documented here for RGCs. Our study provides a computational framework for investigating this issue.

All animal experiments were approved by the Institutional Animal Care and Use Committees (IACUC) at Harvard University. Mice were maintained in pathogen-free facilities under standard housing conditions with continuous access to food and water. Animals used in this study include both males and females. A meta-analysis (not shown) did not show any systematic sex-related effects in either DE genes or cell-type proportions. For scRNA-seq and histology, we used C57Bl/6J (JAX #000664). Embryonic and early postnatal C57Bl/6J mice were acquired either from Jackson Laboratories (JAX) from time-mated female mice or time-mated in-house. For timed-matings, a male was placed with a female overnight and removed the following morning (with the corresponding time recorded as E0.5).

RGCs were enriched from dissociated retinal cells as previously described with minor modifications (Tran et al., 2019). All solutions were prepared using Ames' Medium with L-glutamine and sodium bicarbonate (equilibriated with 95% O₂/5% CO₂), and all spin steps were done at 450 × g for 8 min. Retinas were dissected out in their entirety immediately after enucleation and digested in ~80 U of papain at 37°C, with the exception of some E13 and E14 eyes that were digested whole, followed by manual trituration in ovomucoid solution. Clumps were removed using a 40 μm cell strainer and the cell suspension was spun down and resuspended in Ames + 4% BSA at a concentration of 10 million cells per 100 μl. Cells from E13, E14, E16, and P0 were incubated for 15 min at room temperature with antibodies to Thy1 (also known as CD90) and L1CAM pre-conjugated to the fluorophores APC (Thermo Fisher Scientific #17-0902-82) and PE (Miltenyi Biotec #130-102-243), respectively. Cells were washed with 6 ml of Ames + 4% BSA, spun down and resuspended at a concentration of ~7 million cells/ml, and calcein blue was added to label metabolically active cells.

Viable Thy1- or L1CAM-positive cells were sorted using a MoFlo Astrios sorter into ~100 μl of AMES + 4% BSA. Sorted cells were spun down a final time and resuspended in PBS + 0.1% BSA at a concentration of 500–2000 cells/μl. P5 RGCs were enriched using only CD90, with either magnetic-activated cell sorting (MACS) using large cell columns and CD90 pre-conjugated to microbeads (#130-042-202 and #130-049-101, Miltenyi Biotec) or fluorescence-activated cell sorting with CD90 pre-conjugated to PE/Cy7 (Thermo Fisher Scientific #25-0902-81), or both.

All single-cell libraries were aligned to the UCSC mm10 transcriptomic reference (Mus musculus), and gene expression matrices were quantified using standard protocols described previously. For the single-cell libraries generated using the 10X platform (E13, E14, E16, P0, and P5), these steps were performed using cellranger v2.1.0 (10X Genomics). For the single-cell libraries generated using Drop-seq (P5), we used Drop-seq tools (v1.12; Macosko et al., 2015), following the procedures described earlier (Shekhar et al., 2016). Alignment and quantification was done for each sample library separately to generate a genes × cells expression matrix of transcript counts. These matrices were column-concatenated for further analysis.

We retained cells that expressed at least 700 genes, resulting in 98,452 cells. We also removed genes expressed in fewer than 10 cells. The resulting M genes × N cells matrix of UMI counts $C_{mn}$ was normalized along each column (cell) to sum to 8340, the median of the column sums resulting in a normalized matrix $X_{mn}$ . This was followed by the transformation $X_{mn}\leftarrow \mathrm{l}\mathrm{o}\mathrm{g}(X_{mn}+1)$.

The following procedure was adopted to perform batch correction, dimensionality reduction, and clustering throughout the article. The procedure was first applied on the entire dataset to separate RGCs from other cell classes, and then to RGCs at each age to identify transcriptomically distinct groups.

1. Identification of HVGs: We used the Gamma-Poisson framework described previously to identify HVGs (Pandey et al., 2018). Briefly, we compute for each gene the mean (${\mu }_m$) and the coefficient of variation ($C{V}_m$) for the UMI counts $C_{mn}$ ,

   • ${\mu }_m=\frac{1}{N}{\sum }_n{C}_{mn}$

   • ${\sigma }_m^2=\frac{1}{N}{\sum }_n{\left(C_{mn}-{\mu }_m\right)}^2$

   • $C{V}_m=\frac{{\sigma }_m}{{\mu }_m}$

   • For a given ${\mu }_m$ , the Gamma-Poisson model predicts a ‘null’ coefficient of variation ($C{V}_m^{null}$) arising from a combination of Poisson ‘shot’ noise and large variations in library size, assumed to be due to technical reasons,

   • $C{V}_m^{null}=\frac{1}{{\mu }_m}+\frac{1}{\alpha }$

   • Here, $\alpha$ is the shape parameter of a Gamma-distribution fit to the distribution of normalized library sizes $T_n$ (using the R package MASS),

   • $T_n=\frac{{\sum }_m{C}_{mn}}{{\sum }_{m,n}{C}_{mn}}$

   • In practice, $C{V}_m^{null}$ serves as a tight lower bound for empirically observed values of $C{V}_m$ across the full range of ${\mu }_m$ . This enables us to compute for each gene m, a deviation score $d_m=\log \frac{C{V}_m}{C{V}_m^{null}}$ , quantifying the extent to which its observed coefficient of variation exceeds the predicted null model. HVGs are selected if they satisfy ${\displaystyle d_m>Mean\left(d_m\right)+0.8Std\left(d_m\right)}$.

2. Batch correction and dimensionality reduction: We subsetted the rows of the expression matrix $X_{mn}$ to the HVGs identified in Step 1. As our data comprised cells sampled at different developmental ages as well as multiple biological replicates within each age, we used Liger, a non-negative matrix factorization technique, to embed the data in a reduced dimensional latent space of shared factors (Welch et al., 2019). Liger computes a factorized representation for each matrix that separates ‘shared’ and ‘dataset-specific’ gene expression modules (factors). We use Liger’s normalized H factor loadings for cells to build a nearest-neighbor graph and define clusters.

   • As in any matrix factorization technique, Liger requires the user to choose k, the dimensionality of the latent space. To find k, we use a Random Matrix Theory approach (outlined in Peng et al., 2019). Briefly, k is estimated as the number of eigenvalues of the sample gene–gene correlation matrix that exceed the 99th percentile of the distribution of the largest eigenvalue of a random Hermitian matrix of the same dimensions. This is given by the Tracy–Widom distribution. For these calculations, we used the R package RMTstat.

3. Clustering and 2D visualization: To cluster cells based on transcriptomic similarity, we first built a nearest-neighbor graph on the cells based on their normalized H factor coordinates computed using Liger. The number of nearest neighbors was chosen to be 30. The edges were weighted based on the Jaccard overlap metric, and graph clustering was performed using the Louvain algorithm, as described previously (Shekhar et al., 2016). The normalized H factor coordinates were also used as input to project cells on to a nonlinear 2D space using the Uniform Manifold Approximation and Projection algorithm (UMAP; Becht et al., 2018). Graph construction, clustering, and the UMAP projection were performed using the R packages FNN, igraph, and umap, respectively.

We began by clustering the full dataset combining all ages using the procedure outlined above. We identified groups of clusters corresponding to cell classes, which included RGCs (Rbpms, Slc17a6, Sncg, Nefl), microglia (P2ry12, C1qa-c, Tmem119), photoreceptors (Otx2, Gngt2, Gnb3), amacrine cells (Tfap2a, Tfap2b, Onecut2), anterior segment cells (Mgp, Col3a1,Igfbp7), cycling progenitors (Ccnd1, Fgf15, Hes5), and neurogenic progenitors (Hes6, Ascl1, Neurog2). Deeper annotation (e.g., of RGC type) was not done at this stage. No other cellular classes were identified. Three clusters comprising fewer than 1.2% of the cells expressed markers of more than one class. These were flagged as doublets and removed from further analyses.

RGC precursors at each age were separately analyzed following the clustering pipeline outlined previously. When implementing Liger, each biological replicate was regarded as a separate batch. The nominal clusters identified by the Louvain algorithm were refined as follows:

1. Removing contaminants: Clusters were flagged for further examination if they did not exclusively express RGC-specific markers (e.g., Rbpms, Slc17a6, Sncg, Nefl). These clusters were small (typically <1–2% of cells) and in all cases expressed non-RGC markers (e.g., P2ry12 or Tfap2b). These cells, which likely reflect trace contaminants, were discarded from further analysis.

2. Merging proximal clusters: Transcriptomic relationships between nominated clusters were visualized on a dendrogram computed using hierarchical clustering, as noted above. Neighboring clusters on the dendrogram, which were leaves in a terminal branch, were assessed for differential expression using the MAST DE test (Finak et al., 2015). A gene g was regarded as significantly DE between clusters ${\displaystyle {\mathcal{C}}_a}$ and ${\displaystyle {\mathcal{C}}_b}$ if it satisfied ${\displaystyle \left|\log F{C}_g\left({\mathcal{C}}_a,{\mathcal{C}}_b\right)\right|>0.5}$ and MAST p-value was less than ${10}^{-5}$ (false discovery rate [FDR] corrected), where

   • ${\displaystyle \begin{array}{c}\log F{C}_g\left({\mathcal{C}}_a,{\mathcal{C}}_b\right)=\ln \left(\frac{\left|{\mathcal{C}}_b\right|\sum _{n\in {\mathcal{C}}_a}{X}_{gn}}{\left|{\mathcal{C}}_a\right|\sum _{n\in {\mathcal{C}}_b}{X}_{gn}}\right)\end{array}}$

is defined to be the log-fold change in expression. Clusters that showed fewer than 10 significant DE genes were merged. In this manner, we identified 10 RGC clusters at E13, 12 at E14, 19 at E16, 27 at P0, and 38 at P5. Using MAST, we identified DE genes that distinguished each cluster against the rest at any given age.

We quantified the molecular diversity of RGCs based on clusters at each stage using three measures of population diversity – the Rao index (Figure 2), the Shannon index, and the Simpson index (Figure 2—figure supplement 1). For N clusters with relative frequencies $p_1,p_2,\dots ,p_N$ , these indices are defined as follows:

1. Let $d_{ij}$ be a distance measure between clusters i and j ($0\le d_{ij}\le 1$). The Rao index is defined as

   • We used varying number of genes (≈1200–3000) to calculate $d_{ij}$ . The computed Rao index was insensitive to these variations.

2. The Shannon index is defined as

   • ${\displaystyle H=-\sum _i{p}_i\log p_i}$

3. The Simpson index is defined as

   • ${\displaystyle S=\sum _i{p}_i^2}$

4. While the Rao and Shannon indices increase with diversity, the Simpson index decreases with diversity.

We quantified the mutual separation of clusters at each age using two approaches:

1. Multi-class classification: We trained a multi-class classifier (R package xgboost) at each age on 50% of the cells using their cluster IDs. The remaining 50% of the cells were used to test the learned classifier and estimate a classification error per cluster, which were averaged at each age. As clusters become better separated, the average classification error decreases.

2. Relative positions in PC space: At each age, the top 20 principal component analysis (PCA) coordinates were first standardized by z-scoring. For each cluster C at a given age, we computed two quantities:

   a. $r_C$ , the median of Euclidean distances of each cell from the cluster centroid in the standardized PCA coordinates.

   b. $d_C$ , the median of Euclidean distances of each cell from the centroid of the nearest external cluster.

For a cluster C, a low of value $r_C/d_C$ indicates a higher degree of separatedness. Averaging this metric across all the clusters at a given age quantifies the degree to which clusters are separated in the UMAP representation.

Let $N_{ij}$ denote the number of cells in ${\displaystyle {\mathcal{A}}^T}$ that are part of transcriptomic cluster ${\displaystyle {\mathcal{C}}_j^T}$ and are assigned by $Clas{s}^R$ to reference cluster ${\displaystyle {\mathcal{C}}_i^R}$ . Thus,

$N_{ij}$ defines a contingency table, whose marginal sums are defined as

Let ${\displaystyle N=\sum _{i,j}{N}_{ij}=\left|{\mathcal{A}}^T\right|,}$ the number of cells in ${\displaystyle {\mathcal{A}}^T}$ . Then, the ARI corresponding to the assignments can be evaluated using the following equation:

The ARI ranges from 0 and 1, with extremes corresponding to random association and 1:1 correspondences between ${\displaystyle {\mathcal{A}}^R}$ and ${\displaystyle {\mathcal{A}}^T}$ , respectively. The ARI can technically also take on negative values for certain scenarios, but these are not observed in our data.

As an alternative, we also used the NCE, an information-theoretic measure. The NCE quantifies the extent to which knowledge of the value of $cluste{r}^{`}\left(\mathit{u}\right)$ reduces the uncertainty (measured in information bits) about the value of $cluster\left(\mathit{u}\right)$ for ${\displaystyle \mathit{u}\in {\mathcal{A}}^T}$ .

We introduce probability weights $q_{ij}$ and the corresponding marginals $q_{i,.}$ and $q_{.,j}$ as follows:

The conditional entropy (CE) is then given by the expression:

Note that CE is asymmetric, that is, ${\displaystyle H\left(cluster\left(\mathit{u}\right)|cluste{r}^{\prime }\left(\mathit{u}\right)\right)\ne H\left(cluste{r}^{\prime }\left(\mathit{u}\right)cluster\left(\mathit{u}\right)\right)}$. One notes that H = 0 if for each cluster $i\in \{1,\dots ,r\}$, $q_{ij}={\delta }_{i,k_i}$ for a single cluster $k_i$ , where ${\displaystyle {\delta }_{ij}}$ is the Kronecker delta defined as

Finally, NCE is defined as

where $H\left(cluster\left(\mathit{u}\right)\right)=-{\sum }_j{q}_{.,j}\log q_{.,j}$ is the Shannon entropy. Due to this normalization, NCE values range from 0 to 1, with extremes corresponding to fully specific mapping or random association, respectively, between ${\displaystyle {\mathcal{A}}^R}$ and ${\displaystyle {\mathcal{A}}^T}$ . ARI and NCE are inversely related. Unlike ARI, however, NCE is able to detect specificity in both many:1 and 1:1 mappings. ARI returns a value lower than 1 for specific mappings if the number of clusters in ${\displaystyle {\mathcal{A}}^R}$ and ${\displaystyle {\mathcal{A}}^T}$ is not equal.

ARI and NCE quantify global correspondences between ${\displaystyle {\mathcal{A}}^R}$ and ${\displaystyle {\mathcal{A}}^T}$ . We also computed a local metric, the OF that quantified whether individual reference labels ${\displaystyle {\mathcal{C}}_i^R}$ were distributed in a ‘localized’ or ‘diffuse’ manner between test clusters${\displaystyle \in \left\{{\mathcal{C}}_1^T,{\mathcal{C}}_2^T,\dots {\mathcal{C}}_t^T\right\}.}$

Note that the term ${\displaystyle \frac{q_{ij}}{q_{i,.}}}$ is simply the fraction of the total test cells belonging to test cluster ${\displaystyle {\mathcal{C}}_j^T}$ that are assigned to reference cluster ${\displaystyle {\mathcal{C}}_i^R}$ by the classifier. Defined this way, the term in the square brackets computes an occupation number that ranges from 1 to t and can be interpreted as the number of test clusters that are specifically associated with ${\displaystyle {\mathcal{C}}_i^R}$ . Division by t the number of test clusters, therefore, converts this number into a fraction.

To identify fate relationships among maturing RGCs, we used WOT (Schiebinger et al., 2019), a recently developed framework that is rooted in optimal transport theory (Villani, 2009). WOT does not rely on clustering, and therefore is able to identify ancestor–descendant relationships between any pair of temporally separated RGCs in our data.

At its heart, WOT models cellular transcriptomes $\mathit{u}$ measured at a given age t as a probability distribution in gene expression space ${\displaystyle P_t(\mathit{u})}$ (note that $\mathit{u}$ may represent the original gene expression space or a reduced dimensional embedding estimated via PCA or diffusion maps). This probability distribution evolves with time as cells differentiate and mature. Different temporal measurements collected at times $\dots ,t_{i-1},t_i,t_{i+1},\dots$ represent temporal snapshots of the corresponding cell distributions $\dots ,{\mathbb{P}}_{t_{i-1}},{\mathbb{P}}_{t_i},{\mathbb{P}}_{t_{i+1}},\dots$ . Unfortunately, as each cell can only be measured once, the measurement at different times is from different cells. Therefore, for a particular cell $\mathit{u}$ at time $t_i$ , it is not clear which cell(s) at time $t_{i-1}$ are likely to be its ancestor(s) and which cell(s) at time $t_{i+1}$ are likely to be descendant(s). It is this problem that WOT addresses.

Briefly, for a given pair of consecutive transcriptomic snapshots ${\displaystyle P_{t_i}(\mathit{u})}$ and ${\displaystyle P_{t_{i+1}}(\mathit{v})}$, we wish to estimate the joint distribution ${\mathbf{\Pi }}_{t_i,t_{i+1}}(\mathit{u},\mathit{v})$, representing the probability that a cell having an expression vector $\mathit{u}$ at time $t_i$ transitions to a cell with an expression vector $\mathit{v}$ at time $t_{i+1}$ . ${\mathbf{\Pi }}_{t_i,t_{i+1}}(\mathit{u},\mathit{v})$ is also called the temporal coupling, which, owing to the destructive nature of scRNA-seq assays, is not directly observable. Under the assumption that cells move short distances in transcriptomic space when $\mathrm{\Delta }{t}_i=t_{i+1}-t_i$ is ‘reasonably close,’ WOT estimates ${\mathbf{\Pi }}_{t_i,t_{i+1}}(\mathit{u},\mathit{v})$ as the solution to the following convex optimization problem:

In the above equation,

• ${\displaystyle {\hat{P}}_{t_{i+1}}\left(\mathit{v}\right)}$ is an empirical distribution constructed from ${\displaystyle {\mathcal{A}}^{t_{i+1}}}$ , which denotes the scRNA-seq atlas at $t_{i+1}$ ,

  • ${\displaystyle {\hat{P}}_{t_{i+1}}\left(\mathit{v}\right)=\frac{1}{|{\mathcal{A}}^{t_{i+1}}|}\sum _{{\mathit{x}}_i\in {\mathcal{A}}^{t_{i+1}}}\delta \left(\mathit{v}-{\mathit{x}}_{\mathit{i}}\right)}$

• where ${\displaystyle \delta (\mathit{v}-\mathit{x})}$ denotes the Dirac delta function, a probability distribution placing all its mass at the location $\mathit{x}$.

• ${\displaystyle {\hat{Q}}_{t_i}\left(\mathit{u}\right)}$ is the cell distribution at $t_i$ rescaled by the relative growth rate to account for cell division/death,

  • ${\displaystyle {\hat{Q}}_{t_i}\left(\mathit{u}\right)={\hat{P}}_{t_i}\left(\mathit{u}\right)\frac{g{\left(\mathit{u}\right)}^{t_{i+1}-t_i}}{\int g{\left(\mathit{u}\right)}^{t_{i+1}-t_i}d{\hat{P}}_{t_i}}}$

• Here, ${\displaystyle g\left(\mathit{u}\right)}$ represents the relative growth rate of cell $\mathit{u}$ in the time interval $(t_i,t_{i+1})$ and is estimated within the framework of unbalanced optimal transport (Chizat et al., 2018). For more details, we refer the reader to the supplementary information of Schiebinger et al., 2019.

• $c\left(\mathit{u},\mathit{v}\right)$ is a cost function defined as the Euclidean distance ${\left|\mathit{u}-\mathit{v}\right|}^2$ . The first term of the objective function minimizes the cost function weighted by the temporal couplings, which may be interpreted as the transport distance between the distributions ${\displaystyle {\hat{P}}_{t_i}}$ and ${\displaystyle {\hat{P}}_{t_{i+1}}}$ (also known as the Wasserstein distance).

• The second term on the RHS represents entropic regularization, and $ϵ$ is the corresponding strength. Classic OT identifies ‘deterministic’ couplings in that one cell at $t_i$ is transported to a single cell at $t_{i+1}$ . Introduction of the entropic regularization term makes this problem nondeterministic, capturing the notion that there may exist immature cells whose fate is not completely determined. Our inferences of multipotentiality is directly a consequence of adding this entropic regularization term. Additionally, entropic regularization also makes the problem strongly convex, which is computationally beneficial.

• The third and the fourth terms are features of unbalanced optimal transport, where equality constraints on the marginals (a consequence of mass conservation) are relaxed. ${\lambda }_1$ and ${\lambda }_2$ are the corresponding Lagrange multipliers.

We note that the values of the hyperparameters $ϵ,{\lambda }_1$, and ${\lambda }_2$ are held fixed for all pairwise transport map calculations (E13, E14), (E14, E16), ….

We apply WOT to each pair of consecutive ages $t_i$ and $t_{i+1}$ to estimate the transport map ${\displaystyle {\hat{\mathbf{\Pi }}}_{t_i,t_{i+1}}}$ . Transport maps connecting nonconsecutive time points $t_i$ and $t_{i+k}$ are estimated through a simple matrix multiplication of intermediate transport maps

The transport matrices ${\displaystyle {\hat{\mathbf{\Pi }}}_{t_i,t_j}}$ encode fate relationships between cells at $t_i$ and cells at a later time ${\displaystyle t_j>t_i}$ . These relationships can be analyzed at the level of clusters at $t_j$ to associate each cell ${\displaystyle \mathit{u}\in {\mathcal{A}}^{t_i}}$ with transcriptomically defined cluster. This is particularly useful in estimating the terminal identity of immature RGCs.

Operationally we compute for each cell ${\displaystyle \mathit{u}\in {\mathcal{A}}^{t_i}}$ a ‘cell fate vector’ $f_{t_j}\left(\beta ;\mathit{u},t_i\right),(\beta =1,2,\dots )$ encoding the probabilities that $\mathit{u}$ is associated with cluster ${\displaystyle {\mathcal{C}}_{\beta }^{t_j}}$ at time $t_j$ ,

It is easy to verify that

The cell fate vector $f_{t_j}\left(\beta ;\mathit{u},t_i\right)$ encodes probabilistic associations between the cell $\mathit{u}$ and terminal clusters at ${\displaystyle t_j>t_i}$ indexed by $\beta$. The ‘cluster ancestry vector’ at an earlier time $t_i$ of a cluster ${\displaystyle {\mathcal{C}}_{\beta }^{t_j}}$ at time ${\displaystyle t_j>t_i}$ , denoted ${\mathrm{\Gamma }}_{t_i}(\mathit{u};C_{\beta }^{t_j})$ , is defined as follows:

In a similar vein, the ‘cluster descendant vector’ at a later time $t_o$ of a cluster ${\displaystyle {\mathcal{C}}_{\beta }^{t_j}}$ at a time ${\displaystyle t_j<t_o}$ , denoted ${\mathrm{\Gamma }}_{t_o}\left(\mathit{u};C_{\beta }^{t_j}\right)$ , is defined as

These equations can be used to compute the putative ancestral or descendent cells associated with a cluster $C_{\beta }^{t_j}$ at time $t_j$ .

RGC vectors from all ages were combined, median-normalized, and log-transformed. 2854 HVGs were identified using the Gamma-Poisson model, and WOT was run on this matrix as follows:

wot optimal_transport --matrix RGC_mat.mtx --cell_days cell_day.txt --growth_iters 3 --epsilon 0.005 --out tmaps/RGC

Cell days were specified in cell_day.txt as 0, 1, 3, 6, 11, and 20 for E13, E14, E16, P0, P5, and P56, respectively. We computed trajectories and fates for each age using the following command illustrated for P0:

wot trajectory --tmap tmaps/RGC --cell_set cell_sets.gmt --day 6 --out tmaps/traj_RGC_P0.txt

Fates were computed as

wot fates --tmap tmaps/RGC --cell_set cell_sets.gmt --day 6 --out tmaps/fate_RGC_P0.txt

The above process was repeated for each age.

We define

as the fate coupling between RGC types $\alpha$ and $\beta$ at age $t$. Clearly, $C\left(\alpha ,\beta ;t\right)$ is simply the Pearson correlation coefficient between $f_{P56}\left(\alpha ;\mathit{u},t\right)$ and $f_{P56}\left(\beta ;\mathit{u},t\right)$ , the probabilities that a cell ${\displaystyle \mathit{u}\in {\mathcal{A}}^t}$ is a precursor of $\alpha$ and $\beta$ precursor. Here,

is the mean probability that a cell at age t is a precursor of type $\alpha$. We computed $C\left(\alpha ,\beta ;t\right)$ across all 990 pairs of RGC types at each immature age $t\in \left(E13,E14,\dots ,P5\right).$ The values $C\left(\alpha ,\beta ;E13\right)$ were used as edge weights to visualize the fate coupling network of RGC types using the force-directed layout method (Fruchterman and Reingold, 1991) as implemented in the R package igraph. The node layout was computed using $C\left(\alpha ,\beta ;E13\right)$ values. For other ages, the node layout at E13 was retained but the edges were replotted based on $C\left(\alpha ,\beta ;t\right)$ values at the corresponding age.

We computed a null distribution of $C\left(\alpha ,\beta ;t\right)$ by randomizing the values of $f_{P56}\left(\alpha ;\mathit{u},t\right)$ within each cell $\mathit{u}$ across types. The null values of $C\left(\alpha ,\beta ;t\right)$ rarely exceeded 0.1 and never exceeded 0.2, so only the edges with larger weights were visualized in Figure 5.

For each pair of RGC types $\alpha$ and $\beta$, we fitted a logistic equation to model the decay of pairwise couplings as

The values of t corresponding to E13, E14, E16, P0, and P5 were t = 0, 1, 3, 6, and 11, with $C\left(\alpha ,\beta ;t\right)$ computed as above. We also assumed that $C\left(\alpha ,\beta ;t\right)=0$ at t = 36, corresponding to P30. Thus, six data points were used to estimate two parameters for each of the 180 pairs of RGC types that had nonzero values of $C\left(\alpha ,\beta ;t\right)$ . The nls function from the R package stats was used to estimate ${\beta }_0$ and ${\beta }_1$ . The results are plotted in Figure 5f.

We hypothesized that the specification of a type $\beta$ corresponds to the localization of its precursors in transcriptomic space. The extent of localization for a RGC of type $\beta$ across the time course was calculated as follows. At each age t, we identified the set of precursor RGCs $Prec(\beta ;t)$ showing the highest fate probability corresponding to type $\beta$:

Next, we calculated how the precursors of $\beta$ were distributed across clusters at time t. We computed the OFs (see above) of precursor cells for type $\beta$ across all clusters ${\displaystyle {\mathcal{C}}_k,k=1,2,\dots ,N\left(t\right)}$ at a particular time t ($N\left(t\right)$ is the number of transcriptomically defined clusters at time t):

The localization score for each type $\beta$ at a given time t was defined as

where the index k ranges over the number of clusters at time t. As defined, $Localization(\beta ;t)$ is restricted to be between 0 and 1, with higher values representing a greater specification. We used a logistic model to approximate the localization of each type as

As in the previous section, the nls function was used to estimate the logistic parameters ${\gamma }_0$ and ${\gamma }_1$ . We consider a type $\beta$ as specified if its specification crosses the line $y\left(t\right)=0.95\left(1-1/N\left(t\right)\right)$ . Thus, the specification time for a type $\beta$ is defined as

Note that, as defined, ${\tau }_{sp}$ can be any time point in the interval (E13, P30) corresponding to $t\in \left(\mathrm{0,36}\right)$.

To identify putative ipsilateral- and contralateral-specified RGC precursors at E13, we scored each precursor RGC based on their expression of bona fide ipsilateral genes (Zic2, Zic1, and Igfbp5) and bona fide contralateral genes (Isl2, Fgf12, Igf1) as in Wang et al., 2016. We refer to these as I-RGC and C-RGC scores. Putative I-RGCs were those cells that expressed the I-RGC score at 1.5 standard deviations higher than the mean across all cells, and those that express the C-RGC score at 1.5 standard deviations lower than the mean across all cells. C-RGCs were defined analogously. Many cells did not express either of these marker sets as shown in Figure 7c. These are likely to be RGCs that have not declared their laterality or C-RGCs that are not defined by the expression of Isl2, Fgf12, and Igf1.

WOT was then used to compute the descendants of E13 I-RGCs at all subsequent ages through P56 using the wot fates command introduced above. These descendants were used for two purposes. First, we assessed the proportion of putative I-RGCs across types as in Figure 7d. We also performed a differential gene expression test between putative I-RGCs and the remaining RGCs at all ages, as shown in Figure 7e and f and Figure 7—figure supplement 1d and e.

We compared our data with three scRNA-seq studies that profiled the whole retina during development:

1. (Clark et al., 2019): Count matrices and cell-/gene-level annotations were downloaded from the author’s public repository https://github.com/gofflab/developing_mouse_retina_scRNASeq (Goff Lab, 2021). This dataset contains whole retinal cells sampled at 10 time points (E11, E12, E14, E16, E18, P0, P2, P5, P8, P14) with four of these (E14, E16, P0, P5) common with our study. We excluded P5 from our analysis as only N = 11 RGCs were identified by the authors at this time point.

2. (Lo Giudice et al., 2019): Count matrix corresponding to E15.5 retinal cells was kindly provided by the authors.

3. (Rheaume et al., 2018): Count matrix corresponding to P5 RGCs was downloaded from the online submission.

For consistency with our filtering parameters, we extracted cells based on a cutoff of 700 genes/cell from each of the above datasets. For the Clark et al. dataset, this selected 17,827 cells at E14, 1674 cells at E16, and 8343 cells at P0, respectively (N = 27,844 cells). In these samples, RGCs comprised 19, 28, and 0.45%. For the Giudice et al. dataset, this selected 5218 cells, of which 23% were RGCs.

The Rheaume et al. dataset was directly compared with P5 RGCs collected in this study using supervised classification (Figure 1—figure supplement 1i). Clark et al. and Lo Giudice et al. data were combined with the retinal cells profiled in this study at corresponding time points (25,685 cells at E14; 21,274 cells at E16; and 23,251 cells at P0). Together, this resulted in a 14,350 genes × 103,272 cells expression matrix that was analyzed following the steps outlined previously. In the alignment step, cells from each combination of age and study were considered as a separate ‘batch.’ We visualized the transcriptional heterogeneity of the full dataset using UMAP and used the expression of canonical markers to confirm the co-clustering of cell classes in Figure 1—figure supplement 1 (Rbpms for RGCs, Tfap2b for ACs, Fgf15 for RPCs, and Gngt2 for RPCs).

All scRNA-seq data collected in this study were submitted to the Gene Expression Omnibus (GEO) under GSE185671. The data can be visualized on the Broad Institute’s Single Cell Portal under the identifier SCP1706.

The scripts (written in R) generated for this study are shared at https://github.com/shekharlab/mouseRGCdev, (copy archived at swh:1:rev:ca6a97adabb7bc4ffb2fb1187c78cb277513665c; Shekhar, 2022).