The retina scRNA-seq datasets of the three model organisms, namely, Homo sapiens (human), Mus musculus (mouse), and Danio rerio (zebrafish), were downloaded from Gene Expression Omnibus (GEO) and ArrayExpress. These datasets include developing human retina (GSE138002), adult human retina (E-MTAB-7316), developing mouse retina (GSE118614), adult mouse retina (GSE132229), and developing zebrafish and adult zebrafish retina (GSE122680 and GSE160140, respectively) (Clark et al., 2019; Heng et al., 2019; Lukowski et al., 2019; Lu et al., 2020; Xu et al., 2020; Liu et al., 2022). Specifically, adult human retinal samples, including the whole retina, macular region, and peripheral region, were separated for further analysis.

Homologous gene lists were downloaded from Ensembl BioMart. An aggregated table of “meta-genes” was created to account for gene paralogs and gene duplication events. Each meta-gene may include all gene symbols homologous to a specific human gene (Supplementary Table S1).

The raw unique molecular identifier (UMI) count matrices were converted into a Seurat object using Seurat (v5.0.1) of the R package (Hao et al., 2024). Batch effects among the samples were alleviated using Harmony (Korsunsky et al., 2019). Cells with UMI numbers <1,000 or with feature counts <200 or with mitochondrial percentage >10 were considered low-quality cells and removed. The resulting datasets were normalized using the ScaleData function, and principal component analysis (PCA) was applied. The main cell clusters were identified using the FindClusters function, and 30 PCs were used in subsequent cell cluster analysis with a resolution of 0.8. They were then visualized using 3D t-distributed stochastic neighbor embedding (tSNE) plots. The cell barcode information of datasets was added to categorize every cell into a known biological cell type. Marker genes and differentially expressed genes (DEGs) were identified using the FindConservedMarkers and Findallmarker functions. |avg_log2FC| > 1 and Bonferroni-adjusted p-value (p_val_adj) < 0.05 were the screening threshold of significantly DEGs. We performed Gene Set Enrichment Analysis (GSEA) using clusterProfiler of the R package (Subramanian et al., 2005; Yu et al., 2012).

A total of 277 IRD genes selected from RetNet (https://sph.uth.edu/retnet/) through manual text mining (Supplementary Table S2) included the top 2,000 highly variable genes as input in our gene set. Weighted gene co-expression network analysis (WGCNA) was conducted by scWGCNA to identify functional modules in the co-expression network (Langfelder and Horvath, 2008; Morabito et al., 2021). An adjacency matrix was generated using a soft threshold (β = 4) to reach a scale-free topology. Average gene significance (GS) was calculated to identify the correlation between module eigengenes and a certain cell type. For each gene, Pearson’s correlation between the gene and module eigengene was calculated as module membership (MM). We screened for the hub genes of each module using the threshold |GS| > 0.2 and |MM| > 0.8. The interactive network was visualized using Cytoscape v3.9.0 (Shannon et al., 2003).

Single-Cell rEgulatory Network Inference and Clustering (SCENIC) enabled us to reconstruct a GRN based on co-expression and DNA motif analysis (Aibar et al., 2017). Simple regulons comprise TFs and the set of genes they regulate. SCENIC was performed on 2,000 cells randomly selected from each of the human and mouse whole retina to identify cell-specified regulons. In addition, hg19-500bp-upstream-10species and mm9-500bp-upstream-7species databases were served as RcisTarget (Herrmann et al., 2012; Imrichova et al., 2015). Regulon activity and regulon OFF/ON status were determined according to the area under the curve (AUC) score.

Cell–cell communication analysis was performed using CellPhoneDB (Efremova et al., 2020). Ligand–receptor complexes play crucial roles in development, differentiation, inflammation, and other processes by orchestrating various biological actions. The CellPhoneDB results were visualized using the ggplot2 package.

Trajectory inference of photoreceptors (PRs) and BCs was performed using Monocle 2.24.0 and Monocle 3 1.3.1 (Trapnell et al., 2014; Qiu et al., 2017). After creating a Monocle object, we performed the analysis on the top 2,000 DEGs as ordering genes. The “DDTree” method was used for dimensionality reduction. Cluster and period annotations were projected on the inferred trajectories. Furthermore, Monocle developed BEAM to assess branch-dependent gene expression by formulating the problem as a contrast between the two negative binomial GLMs.

The list of data and method information was generated (Supplementary Methods).

By comparing the single-cell transcriptome, we analyzed a total of 259,087 cells and found that the correlation among most retinal cell types across species was significantly higher than that among cell types from different categories. Following quality control filtering, 107,037 cells from developing retinas were included, with 86,962 cells originating from the whole retina (Figure 1B), 9,070 cells (8.5%) from the macula, and 11,005 cells (10.3%) from the periphery (Figures 1C,D). Additionally, 83,994 cells from developing mouse retinas were categorized into 9 major retinal cell types, while 32,866 cells derived from 8 developing zebrafish retinas were classified into 10 cell types (Supplementary Figure S1A). Adult retinal samples from 3 species were divided into 8 major cell types, comprising 6,143 (human), 22,186 (mouse), and 4,861 (zebrafish) cells (Supplementary Figure S1B).

The top five GO terms enriched by each cell cluster were consistent with our current knowledge about the retina (Figure 1E). The known marker genes for retinal cell types showed a cluster-specific expression pattern (Figure 1F). Cross-species analysis revealed a highly conserved proportion of retinal cell types (Figure 1G). BCs and MGs demonstrated the highest degree of interspecific conservation between humans and mice, while there are obvious differences in the comparative analysis of GO terms of the top 30 markers among the three species (Figure 1H, Supplementary Figure S1C). A pairwise comparison between humans and other species showed that differences in the expression level of markers exist among species. In zebrafish MGs, atoh7 is the most significantly highly expressed gene, while GPX3 and DKK3 are highly expressed genes in the human and mouse retina, respectively. The expression levels of IRDs (SAG, RGR, ITM2B, etc.) varied across species (Supplementary Figure S1A, Supplementary Figures S2A,B). These expression differences may contribute to phenotypic variations between different species. For instance, the mouse model with an APOE mutation recapitulated the clinical feature of AMD patients rather than the zebrafish model (Fletcher et al., 2014; Vessey et al., 2022). Our results provide a basis for constructing IRD models. Notably, the stathmin family genes, stmn1a/b (STMN1_homo1/2) and stmn2a/b (STMN2_homo1/2), were highly expressed in zebrafish MGs. GPX3, APOE, and GADD45B were highly expressed in human MGs, while Dkk3 (DKK3_homo1/2), Vim (VIM_homo1/2), and Junb (JUNB_homo1/2) were highly expressed in mouse MGs.

Significant GO terms of IRD genes included cell development, differentiation, cell homeostasis, metabolic processes, and ion transport (Supplementary Figure S2C and Supplementary Table S3). Among these genes, most RP genes shared high homology across species and displayed similar expression patterns. For instance, the ceramide kinase-like (CERKL) gene showed high expression levels in the rods of humans, mice, and zebrafish. Additionally, we observed that certain genes exhibited species-specific expression, such as Fam161a and NEUROD1 (Figures 2A–C), and differences in cell-specific expression patterns among species, such as CFB and Elovl1 (Figures 2D–F, Supplementary Figures S2D,E).

Most CRD genes (such as GNAT2 and PDE6H) are mainly expressed in cones and rods, and the expression level of cones is higher than that of rods, which is conservative in humans, mice, and zebrafish. However, there were also a few genes that showed species-specific expression. For instance, in humans and mice, the AIPL1 gene was highly expressed in both cones and rods, while it was highly expressed in the RPCs and rods in zebrafish (Supplementary Figures S2D,E).

In addition to the rods and cones, MGs were also identified as significant sites for the distribution of IRD genes in humans and zebrafish, as opposed to mice. PR degeneration in geographic atrophy or choroidal neovascularization can be triggered during the onset of AMD (Deng et al., 2022). The complement pathway, which plays a crucial role in recognizing and mediating the removal of pathogens, debris, and dead cells, has been reported to be associated with AMD (Liu et al., 2010). Among the complement pathway genes, C2, C3, and CFB were highly expressed in MGs, highlighting the previously underestimated effect of MGs on retinal complement homeostasis (Figures 2D,F). However, we also identified that CFH, a complement pathway gene, was mainly expressed in BCs. Moreover, we observed that CSNB genes were highly expressed in either rods (e.g., GRK1 and PDE6B) or BCs (e.g., GPR179, NYX, and TRPM1), whereas CACNA1F and GNB3 were highly expressed in both BCs and PRs in humans (Figure 2A). These results showed that IRD gene expression exhibits cell specificity in humans, suggesting that different mutant genes leading to the same IRDs may be associated with their expression location in the retina. This insight provides clues for the subsequent exploration of the pathological mechanism of IRDs.

Cells were categorized into eight types based on their co-expression of the genes distinguished among these samples (Figure 3A). Subsequently, , 13 modules were identified according to WGCNA, with sizes ranging from 14 genes (salmon) to 348 genes (turquoise) (Figure 3B, Supplementary Figure 2F). Eight of these modules were highly correlated with retinal cell types: HCs (purple), ACs (pink), MGs (yellow), BCs (red), cones (black), rods (brown), RPCs (blue), and RGCs (turquoise) (Figure 3C, Supplementary Figure 3A).

The modules were identified to have significant associations with the recognized function of retinal cells (Figure 3D). Furthermore, many hub genes were identified among IRD genes, indicating a strong interconnection between different IRD genes. As shown n Figure 3E, GNAT1 was recognized as a hub gene in the brown module. Similarly, GNAT2, PDE6H, OPN1LW, RBP4, PDE6C, and CNGB3 were identified in the black module. It is worth noting that GNAT2, PDE6H, PDE6C, and CNGB3 have been associated with ACHM, an IRD characterized by impaired cone PR function. Additionally, MTTP, USH1C, KCNJ13, and RGR were found in the yellow module (Figures 3E–G).

TFs play a crucial role in determining cell fate and directly govern the transcription patterns. We discovered 32 cell-specific regulons (Figures 4A,B), many of which encompass IRD genes. Among the 44 IRD genes identified in the cell-specific regulons, 55% (22) were RP genes distributed in the rod-specific regulons (including CRX and RAX2 regulons). The functions of RP genes in the CRX regulon include visual system development and PR outer-segment formation (Supplementary Figure 3C). Patients with CSNB are characterized by the dysfunction of rods and impaired signaling from PRs to BCs (Varin et al., 2021). Our findings indicate that CSNB genes are predominantly present in rod- or BC-specific regulons, such as SLC24A1 in the CRX regulon, GNB3 in the RAX2 regulon, and GRM6 in the OTX2 regulon (Figure 4C, Supplementary Figure 2E), whereas GRMP was specifically highly expressed only in BCs (Figure 2A).

The MYB regulon was identified as the RPC-specific regulon with target genes NR2E1, ZHX2, LAMA1, and EPHB2 associated in the eye development term. The CEBPB regulon was identified as MG-specific with target genes CLU, TIMP3, VEGFA, etc. (Figure 4C, Supplementary Figure 3C, Supplementary Table S4). We also found the specificity of TF regulons such as that ONECUT1 and LHX specific in HCs, RAX2 and CRX in rods, MEIS2 in ACs, ZNF331 and CEBPB in MGs, and POU2F2 in RGCs (Figures 4C,D). Consistent with human results, the IRD gene Crx and Otx2 regulons were identified as PR-specific in mice (Supplementary Figures 4A–C). The Lhx9 regulon was identified as AC-specific, and the Pou4f2 regulons were identified as HC- and BC-specific, while Sox4 showed high AUC scores specific for both ACs and HCs (Supplementary Figures 4A, D–F). Remarkably, many TFs from the AP-1 family were identified as members of cell-specific regulons. The AP-1 complex is known to play an important role in neuron protection (Feldker et al., 2020). Regulons involving JUNB, MAFF, and ATF3 were specific to MGs, while the JUN regulon was specific to RPCs. The AP-1 family-involved regulons (Jun, Junb, and Fos) were also identified in the mouse retina (Supplementary Figures 4G–I).

The process of vision relies on the cooperation of various cells in the retina, highlighting the importance of understanding the communication between retinal cells. In humans, MGs exhibit the highest frequency of interaction with other cell types (Figure 5A). Conversely, RPCs showed an even higher frequency in mice and zebrafish (Figures 5B,C). Interestingly, we observed that high expression levels of G protein-coupled receptors (GPR37 and GPR37L1) interact with PSAP, a highly conserved glycoprotein that induces differentiation and prevents the death of various neuronal cells through an active domain of saposin C (Figures 5D,E) in MGs. This suggests the potential application of MGs in the development of retinal regeneration therapy. Furthermore, MGs may primarily receive signals from RGCs and HCs through neurotrophic tyrosine kinase receptor type 1 (NTRK1) (Figure 5E). The NTRK1 signal, along with the ligand–receptor pair, FAM3C and CXADR, exhibits striking activity between RGCs and other retinal cells (Figure 5F), indicating that it is involved in the epithelial-to-mesenchymal transition and retinal laminar formation processes (Katahira et al., 2010). Extensive signaling exchange is observed within the Ephs family (Figures 5D,E,G). The MDK signal in human MGs was far lower than that in RPCs, whereas the mdka signal was active in zebrafish MGs (Supplementary Figure 3B).

A total of 46 genes showed significantly different expression in the macular and peripheral regions, with 10 genes enriched in the macula and 36 genes enriched in the periphery (Figure 6A, Supplementary Table S5). Among these genes, HES1, which confers controlled competence of differentiation (Kageyama et al., 2008), was identified as the most upregulated gene in the periphery. We also observed variation in the abundance of cones and MGs from the macula to the periphery (Figures 6B,C). According to the results, IRD genes CUCA1B and CNGA1 were highly expressed in macular rods, while many mitochondrial genes, including MT-ND4, were highly expressed in macular cones, confirming the heightened oxygen consumption and metabolic activity in the macula (Figures 6B,D). CRABP2, the highest expression gene in the periphery, functions in transporting retinoic acid to the receptor in the nucleus and regulating cell proliferation and apoptosis (Feng et al., 2019). Furthermore, many ribosomal genes, such as RPS27, were highly expressed in peripheral cones, indicating potential active transcription and translation activities (Figure 6B) (Gnanasundram and Fahraeus, 2018). Variations in the gene expression of cones and MGs were also observed from the macula to the periphery (Figures 6B,C,E). Positive regulators of the cell population proliferation term were enriched in the periphery, whereas respirasome and ion transport terms were enriched in the macula (Figure 6F). Consistent with a large number of mitochondrial genes enriched in macular cones, the cellular respiration term was enriched in macular cones from the macula, whereas the translation and cytosolic ribosome term was enriched in the periphery (Figure 6G). The term vascular endothelial cell migration, related to the AMD phenotype, was enriched in macular MGs, while metal ion balance-related terms were enriched in the periphery (Figure 6H).

The cell trajectory was determined using cones as the initiating cells, with the sample period serving as a reference in humans (Figure 7A). An increase in the expression of genes associated with cellular respiration and oxidative phosphorylation was observed in the PR of cell fate 1 compared to cell fate 2. Both branches showed enrichment in visual system development and the metabolic process (Figure 7B). The same trajectory trend was identified in developing mice (Supplementary Figure 5A). In adult zebrafish and human retinas, PRs are characterized by unbiased trajectory trends (Supplementary Figures 5B,C).

We observed variations in the expression of many IRD genes along the pseudotime in PRs. For example, LCA genes OTX2 and AIPL1 exhibited high expression at the beginning, and the LCA gene LCA5, Usher genes ABHD12 and CIB2, and CRD genes RIMS1 and CNGA3 showed the earliest expression patterns. RP genes CRX and TULP1 exhibited a gradual increase over time, whereas many other RP genes (RHO, RBP3, PEX7, etc.) were identified as terminal expression patterns (Figure 7C, Supplementary Figures 5D,E). Our results further support the genetic overlap between LCA and late-onset PR degeneration due to the shared pathogenic genes (Wright et al., 2010). The IRDs (GNAT2, PDEFG, and CSPG4) displayed highly consistent expression in humans and mice (Supplementary Figure 5F). Some IRD genes, including GNAT1, MT-ND4, and PDE6A, exhibited highly dynamic expression during the differentiation of PRs in humans, whereas these patterns were less dynamic in mice.

We re-clustered BCs of human developing retinal samples, generated eight clusters, and annotated them based on the known maker genes. Clusters 4 and 6 consist of a mix of rod and cone BCs. Along the BC trajectory, the clusters containing rod BCs showed a terminal differentiation pattern (Figure 7D). It has been reported that a high level of OTX2 is associated with the development and differentiation of BCs (Yamamoto et al., 2020). Consistently, OTX2 (LCA gene) showed high expression at the beginning of the trajectory, as well as in the PR trajectory. We also observed similar expression patterns for NEUROD4 and RORB, suggesting their potential role in the development and differentiation of BCs (Figure 7D). In contrast, CSNB gene TRPM1 showed terminal expression patterns, which are essential for the development of rod BCs and their synaptic connections with subsequent neurons (Kozuka et al., 2017). Similarly, we identified that PCP2, ISL1, SEPINI1, and CPLX3 exhibited terminal mode characteristics (Figure 7E). PCP2 was highly expressed in rod BCs at the end of trajectory. The retina-specific splice variant of PCP2, Ret-PCP2, accelerates the light response of rod BCs by modulating the mGluR6 transduction cascade (Sulaiman et al., 2010). In mice, Isl1 orchestrates the early differentiation and maintenance of various cell types in the retina across different vertebrates (Bejarano-Escobar et al., 2015).

The tissue characteristics of the retina make it extremely difficult to obtain both healthy adult samples and diseased adult samples from humans, and the lack of patient samples also makes it impossible for us to verify the observation results with patient samples at the single-cell levels. With the advancement of iPSC and retinal organoid technology, it is expected that the difficulty of obtaining adult retinal samples can be greatly improved in the future. In addition, while most cell types are balanced among samples, it is difficult to distinguish physiologically significant processes or technical effects. The intrinsic sparsity of scRNA-seq data makes it difficult to annotate cell types with low RNA levels. In the future, enrichment of single-cell references can help refine labeling for less common cell types.