Evolution of neuronal cell classes and types in the vertebrate retina

Abstract

The basic plan of the retina is conserved across vertebrates, yet species differ profoundly in their visual needs¹. Retinal cell types may have evolved to accommodate these varied needs, but this has not been systematically studied. Here we generated and integrated single-cell transcriptomic atlases of the retina from 17 species: humans, two non-human primates, four rodents, three ungulates, opossum, ferret, tree shrew, a bird, a reptile, a teleost fish and a lamprey. We found high molecular conservation of the six retinal cell classes (photoreceptors, horizontal cells, bipolar cells, amacrine cells, retinal ganglion cells (RGCs) and Müller glia), with transcriptomic variation across species related to evolutionary distance. Major subclasses were also conserved, whereas variation among cell types within classes or subclasses was more pronounced. However, an integrative analysis revealed that numerous cell types are shared across species, based on conserved gene expression programmes that are likely to trace back to an early ancestral vertebrate. The degree of variation among cell types increased from the outer retina (photoreceptors) to the inner retina (RGCs), suggesting that evolution acts preferentially to shape the retinal output. Finally, we identified rodent orthologues of midget RGCs, which comprise more than 80% of RGCs in the human retina, subserve high-acuity vision, and were previously believed to be restricted to primates². By contrast, the mouse orthologues have large receptive fields and comprise around 2% of mouse RGCs. Projections of both primate and mouse orthologous types are overrepresented in the thalamus, which supplies the primary visual cortex. We suggest that midget RGCs are not primate innovations, but are descendants of evolutionarily ancient types that decreased in size and increased in number as primates evolved, thereby facilitating high visual acuity and increased cortical processing of visual information.

Fig. 1. Conserved retinal structure across vertebrates.

a, Cartoon of a section through a vertebrate retina showing the arrangement of its six major cell classes: photoreceptors (including rods (r) and cones (c)), horizontal cells (HC), bipolar cells (BC), amacrine cells (AC), retinal ganglion cells (RGC) and Müller glia (MG). The outer segments of rods and cones (OS), outer nuclear layer (ONL), inner nuclear layer (INL) and ganglion cell layer (GCL)—which contain cell somata—are indicated, as are the outer (synaptic) layer (OPL) and inner plexiform layer (IPL). b, Phylogeny of the 17 vertebrate species analysed in this work. The scale bar on the right indicates estimated divergence time. c, Sections from retinas of eight species immunostained for RBPMS (a pan-RGC marker), CHX10 (also known as VSX2) (a pan-bipolar cell marker) and AP2A (also known as TFAP2A) (a pan-amacrine cell marker) and stained with the nuclear stain DAPI. Scale bars, 25 µm. Figures are representative of images from three retinas.

Fig. 2. Class- and subclass-specific transcriptomic signatures.

a, Heat map showing average expression of marker genes (columns) within each major cell class in 17 species (rows). Rows are grouped by cell class (left). Within each class, species are ordered as in Fig. 1b. Grey tiles indicate data that are missing owing to the absence of the corresponding orthologue in the species annotation. Colours indicating cell class are uniform in a–e. PR, photoreceptor. b, Cross-correlation matrix (Spearman) of pseudobulk transcriptomic profiles for the 16 jawed vertebrates. Rows and columns are grouped by class, and then ordered by phylogeny within a class. A total of 4,560 1:1 gene orthologues were used to calculate the correlation values. c, As in b, with rows and columns grouped by species instead of class. Matrices including lamprey are shown in Extended Data Fig. 7c,d. d, Left, uniform manifold approximation and projection (UMAP) embedding of integrated cross-species data, with points indicating class identity (left) or expression levels of subclass-specific markers (right). GAD1, a marker for GABAergic amacrine cells, is also expressed by some horizontal cells, and ISL1, a marker for ON bipolar cells, is also expressed by some RGCs, horizontal cells, and amacrine cells. Details of gene expression by species are shown in Extended Data Fig. 8d. e, Pairwise mean squared divergence of class-specific pseudobulk gene expression profiles between species (y axis) increases with evolutionary distance, as estimated by substitutions per 100 bp (x axis). Data from mammals, chicken and lizard are included. Data including zebrafish are presented in Extended Data Fig. 7e. Solid lines represent power law (y = ax^b) regression fits. Across the graphs, $a\in \left[0.33,0.47\right]$ and $b\in [0.23,0.35]$. The coefficient of determination (R²) values range from 0.75 to 0.92. Source Data

Fig. 3. Multispecies integration of bipolar cells.

a, UMAP of mammalian bipolar cells computed with the raw (left) and integrated (right) gene expression matrices. Cells are coloured by species of origin. b, Feature plots showing expression within the integrated space of the rod bipolar cell marker PRKCA, the ON bipolar cell marker ISL1, and the OFF bipolar cell marker GRIK1. c, As in a, but with cells coloured by orthotype identity. d, Left, confusion matrix showing the percentage of cells from each mouse bipolar cell-type mapping to each mammalian bipolar cell orthotype. Each column sums to 100%. See Extended Data Fig. 9a for a higher magnification view. Centre, confusion matrix showing specific mapping between mammalian bipolar cell orthotypes and bipolar cell clusters within each mammalian species (Extended Data Figs. 1 and 3–6). Right, confusion matrices showing the mapping of bipolar cell clusters in lizard, chick and zebrafish to the mammalian bipolar cell orthotype. Mapping that includes non-mammalian orthotype is shown in Extended Data Fig. 9e. e, Dot plot showing differentially expressed genes within each bipolar cell orthotype. The size of the dot represents the number of mammalian species (out of 13 mammalian species in total) that express the gene in at least 30% of cells in the corresponding orthotype, and the colour represents normalized expression level. f, Confusion matrix showing the species bipolar cell clusters (columns) that map specifically to the orthotype oRBC and oBC1B. Bipolar cell types are named on the basis of their species of origin and within-species bipolar cell cluster ID (Extended Data Figs. 1 and 3–6); for example, Peromyscus bipolar cell cluster 1 is called ‘Per_1’. g,h, Confusion matrix showing mapping of mammalian photoreceptor (g) and horizontal cell (h) types to orthotype. Format as in d, centre. Source Data

Fig. 4. Multispecies integration of retinal ganglion cells.

a, Integrated UMAP of RGCs from 12 mammals (cow was excluded owing to the paucity of RGC data.). Cells are labelled by species of origin. For primates, cells from fovea and periphery are plotted separately. b, As in a, with RGCs labelled by orthotype. c, Dot plot showing differentially expressed genes within each RGC orthotype. Representation as in Fig. 3e. d, Left, confusion matrices showing that species-specific RGC clusters (Extended Data Figs. 1 and 3–6) map to mammalian RGC orthotypes in a specific fashion. Representation as in Fig. 3d, centre, except that clusters from fovea (F) and periphery (P) are mapped separately for primates. Right, confusion matrices showing the mapping of RGC clusters (columns) in lizard, chick and zebrafish to the 21 mammalian RGC orthotypes. Mapping to the single non-mammalian RGC orthotype is shown in Extended Data Fig. 10d. e, Left, confusion matrix showing that mouse RGC types (rows; naming as in ref. ) belonging to transcription factor-based subsets (colours) map to the same orthotypes (columns). Right, dot plot showing specific expression patterns of subclass-specific transcription factor-encoding genes in orthotypes. Representation as in Fig. 3e. Source Data

Fig. 5. Mammalian orthologues of midget and parasol RGCs.

a, Confusion matrix showing RGC clusters from different species that map specifically to oRGC1, oRGC4, oRGC5 and oRGC2, which contain OFF and ON midget RGCs (MGCs) and OFF and ON parasol RGCs (PGCs). Representation as in Fig. 3f. Column names corresponding to primate midget and parasol types are shown in red, and mouse α-RGC types are shown in blue. b, Schematic delineating morphological and physiological similarities between primate and midget RGCs and their α-RGC orthologues. Orthotypes (OTs) of each pair as well as the orthology among bipolar cell types that innervate them are also shown. Morphologies of neuronal types were created on the basis of published data (Supplementary Note 1). Within each pair, the left column corresponds to primate types and the right column corresponds to mouse types. c, FLDA projection of the scRNA-seq data for primate midget and parasol types and mouse α-RGC types onto the corresponding 3D space, with axes representing species, polarity and kinetics (see Supplementary Note 2). d, Matching MGCs and PGCs to mouse types by GAGE. Inset, given sets of mouse and primate RGC types, the model fits the arrangement of their cluster centroids in gene expression space by assuming a shape that is simply shifted to the other species via a linear translation. Symbols mark the four response types: circle, sustained; square, transient; open, ON; filled, OFF. The graph is a histogram of the fraction of explained variance showing for each proposed combination of four mouse cell types how well the resulting shape fits the macaque RGC geometry. The red bar shows the set of four α-RGC types. Green bars show combinations containing three α-RGC types. Grey bars, remaining sets of four mouse cell types as shown in Supplementary Table 4. e, Relative proportion of OFF and ON midget RGC orthologues in mammalian species based on frequencies of cells in oRGC1 and oRGC4. Source Data

Extended Data Fig. 1. snRNA-seq data from the fovea/macula and peripheral retina of healthy human donors (n = 18).

a. UMAP embedding of nuclei (n = 184,808) from the central and peripheral retina of healthy human donors, with individual points colored by cell class. PRs have been divided into rod and cone subclasses, and ACs have been divided into GABAergic and glycinergic subclasses. b. Same as a, with points colored by sample identity. c. UMAP embedding of RGC nuclei (n = 80,032) from the foveal and peripheral retina of healthy human donors, with individual points colored by type identity. Only ON and OFF midget ganglion RGCs are labeled. d. UMAP embedding of non-midget RGC nuclei (n = 6615) from c, with individual points colored by type identity. ON and OFF parasol ganglion cells are labeled. e. UMAP embedding of BC nuclei (n = 9126) from the fovea and peripheral retina of healthy human donors, with individual points colored by type identity. f. Dotplot showing expression of cell class-specific markers (columns) in the human clusters (rows). The size of each dot represents the fraction of cells in the group with non-zero expression, and the color represents expression level. The six classes are MG, HC, PR (subdivided into Rod and Cone), AC (subdivided into Gabaergic ACs (GabaAC) and glycinergic ACs (Gly AC)), BC and RGC. Only BCs and RGCs have been subclustered. Rows corresponding to BC and RGC clusters are ordered based on hierarchical clustering (dendrograms, left). Barplot on the right of the dotplot depicts the relative frequency of each cluster within a class (colors). The rightmost heatmap depicts the distribution of each cluster across biological replicates (columns).

Extended Data Fig. 2. Nuclear enrichment strategies for retinal ganglion cells (RGCs) and bipolar cells (BCs).

a. Examples of gating strategy in fluorescent activated cell sorting (FACS) experiments for collecting single nuclei labeled with either PE-conjugated NEUN, which enriches RGCs, or APC-conjugated CHX10 (also known as VSX2), which enriches BCs. Data shown are representative from experiments in the pig retina. NEUN and CHX10-based enrichment resulted in ~90% yield for RGCs and ~95% yield for BCs. b. Same as panel a, for human macular retina samples. NEUN-based enrichment resulted in ~90% yield for RGCs; BCs were not analyzed in this experiment. c. Brightfield image showing the morphology and integrity of FACS-purified nuclei. d. Confocal image of DAPI stained FACS-purified nuclei. e. Retinal sections from six species show that PE-conjugated NEUN (red) and APC-conjugated CHX10/VSX2 labels RGCs and BCs, respectively. Retinal sections were co-stained for DAPI (blue) to visualize nuclei. Scale bar, 50 μm. Images in panels a–e representative of n $\ge$ 3 experiments.

Extended Data Fig. 3. Summary of cell type atlases for tree shrew, sheep, cow, and pig.

a. Dotplot showing expression of cell class-specific markers (columns) in the tree shrew (n = 3 animals; 71,571 nuclei) clusters (rows). The size of each dot represents the fraction of nuclei in the group with non-zero expression, and the color represents expression level. The six classes are MG, HC, PR (subdivided into Rod and Cone), AC (subdivided into GABAergic AC (GabaAC) and glycinergic AC (Gly AC)), BCs and RGCs. Only BCs and RGCs have been subclassified through a within-species integration and clustering analysis (Methods). Rows corresponding to BC and RGC clusters are ordered based on a hierarchical clustering analysis (dendrograms, left). Barplot on the right of the dotplot depicts the relative frequency of each cluster within a class (colors). The rightmost heatmap depicts the distribution of each cluster across samples (columns). Panels b-d depict the same information as panel a for sheep (n = 6 animals; 65,490 nuclei) (b), cow (n = 6 animals; 75,794 nuclei) (c), and pig (n = 4 animals; 49,955 nuclei) (d). Note that in this figure, as well as Extended Data Figs. 1 and 4–6, the proportions shown accurately report our data but do not necessarily represent the true endogenous proportions. This is because in many cases we depleted photoreceptors or enriched BCs or RGCs to obtain sufficient numbers of rare cell types (see Methods).

Extended Data Fig. 4. Summary of cell type atlases for Peromyscus, ferret, opossum, and brown anole lizard.

Panels a-d depict the atlases (as in Extended Data Fig. 3) for peromyscus (n = 3 animals; 44,223 cells) (a), ferret (n = 2 animals; 49,972 cells) (b), opossum (n = 5 animals; 76,763 nuclei) (c), and brown anole lizard (n = 3 animals; 42,848 nuclei) (d).

Extended Data Fig. 5. Summary of cell type atlases for Rhabdomys, squirrel, marmoset and sea-lamprey.

Panels a-d depict atlases (as in Extended Data Fig. 3) for Rhabdomys (n = 2 animals; 65,338 nuclei) (a), squirrel (n = 1 animal; 22,821 cells) (b), marmoset (n = 2 animals; 52,559 cells) (c), and Sea-lamprey (n = 2 animals; 18,928 cells) (d).

Extended Data Fig. 6. Summary of cell type atlases for macaque, mouse, chick and zebrafish.

Panels a-d depict atlases (as in Extended Data Fig. 3) for macaque (n = 4 animals; 146,054 cells) (a), mouse (n = 10 animals; 51,162 cells) (b), chick (n = 4 animals; 34,788 cells) (c), and zebrafish (n = 15 biological replicates; 60657 cells) (d). Cluster labels are consistent with published annotations^,,,,. Each biological replicate in zebrafish involved a pooling of eyes from multiple (5-8) fish.

Extended Data Fig. 7. Evolutionary conservation of retinal classes.

a. Dendrogram showing transcriptional relationships among pseudobulk expression vectors following integration. Each node is a cell class within a particular species. Dendrograms were computed via hierarchical clustering analysis (correlation distance, average linkage). b. Same as Fig. 2d, with cells colored by species of origin. Inset shows a magnified region containing samples from all species. c. Cross-correlation matrix (spearman) of class- and species-specific cell-averaged profiles for all 17 vertebrates (compare with Fig. 2b). Rows and columns are grouped by class, and then ordered by phylogeny within a class. d. Same as panel c, but rows and columns grouped based on species instead of class (compare with Fig. 2c). e. Pairwise mean-squared distance of class-specific cell-averaged gene expression profiles between all 16 jawed vertebrate species (y-axis) increases with evolutionary divergence, as estimated by substitutions per 100 bp (x-axis) (compare with Fig. 2e). Gray shaded regions demarcate species pairs involving zebrafish. Solid lines represent power law (y = ax^b) regression fits. Across the panels, $a\in [0.34,0.47]$ and $b\in [0.29,0.45]$. The coefficient of determination (R²) values range from 0.79-0.93.

Extended Data Fig. 8. Evolutionary conservation of retinal subclasses.

a. UMAP embedding of integrated cross-species data (as in Fig. 2d), highlighting PR subclasses cones and rods. Insets show feature plots of cone-specific (top) and rod-specific (bottom) transcription factors (TFs). b. Same as panel a, for AC subclasses GABAergic ACs (GabaAC) and glycinergic ACs (GlyAC). Insets show feature plots of a GABAergic TF MEIS2 and a glycinergic TF TCF4. c. Same as panel a, for BC subclasses ON BCs and OFF BCs. Insets show feature plots of OFF BC-specific (top) and ON BC-specific (bottom) transcription factors (TFs). d. Heatmap showing average expression of subclass-specific genes (columns) within the six subclasses across 17 species (rows). Rows are grouped by subclass (annotation bar, left). Within each subclass, species are ordered as in Fig. 1b, with top and bottom nodes in each dendrogram corresponding to lamprey and human, respectively (corresponding to right and left in Fig. 1a). Gray tiles correspond to missing orthology information. e. Cross-correlation matrix (spearman) of subclass- and species-specific pseudobulk transcriptomic profiles for all 16 jawed vertebrates. Rows and columns are grouped by subclass, and then ordered by phylogeny within a class. Lamprey was excluded due to paucity of shared orthologs. f. Same as panel d, but rows and columns grouped based on species instead of subclass.

Extended Data Fig. 9. Bipolar Cell OrthoType analysis including non-mammals.

a. Confusion matrix showing the rationale behind naming mammalian BC OTs (rows) based on the mapping patterns of mouse BC types (columns). Representation as in Fig. 3d, with each column summing to 100%. OT BC8/9 contains mappings from both mouse BC8 and BC9, which are transcriptionally proximal. b. Barplot showing within-species relative frequencies (y-axis) of the 13 cone BC OTs within each mammalian species (x-axis). The foveal and peripheral data from primates are plotted separately. c. Integrated UMAP of BCs from all 16 jawed vertebrates. Cells are colored by species of origin. Lamprey, a jawless vertebrate, was excluded from the analysis due to the paucity of shared orthologous genes. d. Same as c, with cells colored by OT identity. The integration of all jawed vertebrates recovers all the mammalian BC OTs listed in Fig. 3c, but additionally identifies two OTs enriched for non-mammalian BCs from chick, lizard and zebrafish. The two OTs, named NM_OFF and NM_ON, are enriched for OFF and ON BCs from non-mammals (also see panel e). e. Confusion matrices showing the mapping of species-specific BC clusters (columns) to BC OTs (rows) identified by integrating BCs from all jawed vertebrates (panel c). Representation as in Fig. 3d’. Mammalian BC clusters predominantly map to the mammalian OTs (rows 1-14), and the pattern of mapping is similar to Fig. 3d. Chick, Lizard and Zebrafish BCs largely map to the non-mammalian OTs NM_OFF and NM_ON (rows 15-16). f. Dotplot showing species-specific genes (columns) expressed in RBC orthologs in mammals (rows). The size and color of each dot represent the percentage of cells within the species cluster expressing the gene and the average expression level, respectively.

Extended Data Fig. 10. Retinal Ganglion Cell OrthoType analysis including non-mammals.

a. Barplot showing within-species relative frequencies (y-axis) of the 21 RGC OTs within mammalian species (x-axis) (Fig. 4b). The foveal and peripheral data from primates are shown separately. Cow is excluded due to the paucity of data. b. Integrated UMAP of RGCs from all 15 jawed vertebrates (excluding cow). Cells are colored by species of origin. For primates, fovea and periphery are plotted separately. c. Same as b, with cells colored by RGC OT. OTs 1-21 map 1:1 to the mammalian OTs in Fig. 4b, but we recover an additional OT (NM) predominantly containing non-mammalian RGCs from chick, lizard and zebrafish (also see panel d). d. Confusion matrices showing the mapping of species RGC clusters (columns) to RGC OTs (rows) identified by integrating RGCs from all jawed vertebrates (panel c). Representation as in Fig. 4d. Mammalian RGC clusters predominantly map to the mammalian OTs (rows 1-21), and the pattern of mapping is similar to Fig. 4d. Except for ipRGCs, chick, lizard and zebrafish RGCs largely map to oRGC_NM (row 22). e. Confusion matrix showing the species-specific RGC clusters (columns) that map to the oRGC8 and 9, corresponding to ipRGCs. Representation as in Fig. 3f. Annotation bar (bottom) highlights species-specific RGC clusters that express OPN4 and EOMES, a transcription factor expressed selectively by ipRGCs^,.

Extended Data Fig. 11. Midget and Parasol OTs.

a. Dotplot showing examples of DE genes across OT1-4 and their expression across orthologous species-specific clusters. The size and color of each dot represent the percentage of cells within the species cluster expressing the gene and the average expression level, respectively. Column order as in Fig. 5a. b. Relative proportion of parasol RGC orthologs in mammalian species based on the frequencies of cells in oRGC2 and oRGC5.

Extended Data Fig. 12. Factorized Linear Discriminant Analysis (FLDA) and Geometric Analysis of Gene Expression (GAGE).

a. FLDA workflow and eigenvalue analysis. The gene expression matrices of primate and mouse RGCs were combined by their shared orthologous genes. Highly variable genes were selected, and PCA was applied to remove multicollinearity. FLDA was performed on different combinations of mouse RGC candidates with known polarity and kinetics listed Supplementary Table 4. The combinations were ranked based on their FLDA eigenvalues, which measures the variance along each attribute captured in the projection. b. Visualization of the FLDA projection (Fig. 5c) along the 2D subspace corresponding to polarity (x-axis) and kinetics (y-axis). c. Scatter plot of the FLDA eigenvalues for the kinetics (y-axis) vs. polarity (x-axis), measuring the magnitude of the variance corresponding to these attributes captured in the projection. Inset highlights the top four matches (numbered 1-4) from the 432 combinations of 4 mouse types shown in Supplementary Table 4. d. Mouse RGC types present within the top four combinations out of the 432 combinations in panel c. The top matched set contains all four α-RGC types; the next three include 3 α-RGC types. e. Geometric analysis of gene expression (GAGE) in which primate MGCs and PGCs are compared to all combinations of 4 mouse RGC types (45 choose 4 * 4! = 3,575,880) rather than only the 432 curated combinations used to generate Fig. 5d. Grey bars: histogram of scores for all sets of 4 mouse types. Red bar highlights the set of 4 α-RGC types with the correct matching of polarity and kinetics with the primate types, also marked by the red arrow located at a score of x = 0.657. The bulk of the distribution is approximated as a Gaussian with mean 0.50 and standard deviation 0.0374 (blue line). The 4 α-RGC fit has the second highest score among ~3.6 million candidates. The null hypothesis that this arises by chance has a p-value of p < 10⁻⁶ based on a one-sided Student’s t-test. The top scoring combination with a score of 0.658 involves mouse RGC types C18, C7, C39 and C8 corresponding to the ON PGC, ON MGC, OFF PGC and OFF MGC respectively. Of the four mouse types, two – C18 and C8 - have been physiologically characterized to exhibit sustained ON responses, which violates their expected phenotypic correspondence to ON PGC (ON transient) and OFF MGC (OFF sustained).