Single-cell transcriptomic atlas of zebrafish retina

To establish a benchmark from which to understand changes associated with retinal degeneration, we first developed a cell atlas of the WT adult zebrafish retina using single‐cell 3′ RNA sequencing (scRNA‐seq). The single-cell library was derived from pooled cells from retinas of three WT fish, two females and one male, at the age of 7 months. The library was developed using the 10X Genomics Chromium platform. We have repeated the experiments twice with different single-cell chemistries V2 and V3. Both runs provided similar cluster patterns and cell types. The data presented here are from the V3 chemistry as this has improved Median Genes per Cell and Median UMI Counts (unique transcripts) per Cell. Following quality control and filtering using the Seurat (V3) package [28], the final dataset contained 13,552 cells. The scRNA‐seq data were initially analyzed using an unsupervised graph clustering approach implemented in Seurat to classify individual cells into cell populations according to similarities in their transcriptome profiles. Figure 1A shows a tSNE plot visualization of the clusters in two-dimensional space. Overall, the cells were classified into 23 transcriptionally distinct clusters.

Fig. 1

Single-cell RNA seq analysis of WT zebrafish retina. A tSNE dimensionality reduction and visualization of transcriptional profiles of 13,552 retinal cells from WT zebrafish. B Heatmap showing the top 5 genes enriched for each cluster

The clusters were classified into different retinal cell types using markers identified in previous studies (Table S1). We were not able to assign any cellular identity to one cluster as it expressed markers from multiple retinal cell types. This cluster is not shown in the figure and was not used for further analysis. All the major retinal cell types are represented, including five clusters of amacrine cells, three clusters of microglia, two clusters each of cones, bipolar cells, retinal pigmented epithelium (RPE) and Müller glia, and one cluster each of rods, horizontal cells, retinal ganglion cells, neurogenic progenitor cells (NPCs), oligodendrocytes and retinal progenitor cells (RPCs). Each cell type cluster expresses a unique set of genes along with some genes shared between clusters that define it as a specific cell type. The top 5 genes that are specifically enriched in each cell type are displayed in the heatmap shown in Fig. 1B; a list of the top 10 genes in each cluster is listed in Table S2.

Cellular changes in a zebrafish model of retinitis pigmentosa

To identify transcriptional changes that occur during chronic retinal degeneration and regeneration, we performed the same type of analysis on a single-cell transcriptome library prepared from the adult retina of a transgenic zebrafish model of retinitis pigmentosa expressing P23H mutant rhodopsin [26]. Following quality control and filtering, our final dataset contained 15,511 cells. We compared the WT and P23H datasets using an integration algorithm, which aims to identify shared cell states that are present across different datasets (Fig. 2A). Using the integrated analysis, we detected all the major cell types that are present in the WT dataset, but there were differences in the number of cells in different cell types (Fig. 2A, B). There was a reduction in the number of mature rod photoreceptors (cluster 3) in the P23H retina, reflecting the ongoing degeneration occurring in this model system. In addition, a class of newly differentiated rods (cluster 13) that shows both progenitor and rod markers that did not appear in the clustering of the WT dataset was present. While some cells from the WT retina fell into this cluster, this class was threefold more abundant in the P23H retina, reflecting the continuous regeneration underway in the P23H retina (Fig. 2B, C). In keeping with this, we also detected a fourfold higher number of retinal progenitor cells (RPCs) in the P23H retina compared to the WT (Fig. 2C, cluster 14). Interestingly we saw an eightfold decrease in retinal pigmented epithelial cells (RPE, clusters 10 and 15) compared to WT (Fig. 2C).

Fig. 2

Single-cell RNA seq analysis of integrated WT-P23H zebrafish retina datasets. A Overlap of tSNE dimensionality reduction and visualization of single-cell transcriptomes from WT and P23H transgenic fish retina. B Separate display of P23H and WT tSNE maps from integrated data set with cell types labeled. C Relative proportions of each cell type found in the P23H and WT retinal samples

Numerical comparisons of cell numbers in single-cell transcriptome data present challenges for interpretation. Cells of different types are not captured with equal efficacy and tissue dissociation conditions can favor the preservation of some cell types and the loss of others (e.g., the larger number of bipolar cells than photoreceptors in our datasets). To reinforce the estimates of relative changes in the numbers of cells of certain types, we also examined the numbers of cells in an earlier dataset using 10 × V2 chemistry prepared using different animals of comparable age. Table S3 shows that trends in the relative ratios of cell numbers are conserved between datasets. In this dataset, we find that rods are more than twofold reduced in the P23H, whereas the newly formed rods are ninefold higher in P23H. Our dataset from V2 chemistry also showed a reduction of RPE cells in P23H compared to the WT (2.7-fold decrease in P23H) (Table S3).

Rod photoreceptors display enhanced oxidative metabolism, oxidative stress, and synaptic remodeling

Our previous study showed that rod photoreceptors undergo continuous cell death and that the number of rods is almost fourfold less in the P23H retina than in the WT [26]. Our current single-cell study also reflects the decrease in the number of rods in the P23H dataset compared to the WT (Fig. 3A). Violin plots of the expression of individual genes in a cluster allow us to compare gene expression levels in P23H and WT conditions. Figure 3B shows violin plots of the rod dominant genes rhodopsin (rho) and rod arrestins (saga and sagb). The expression of all of these genes was conserved between WT and P23H. We performed functional pathway analysis using Cystoscope Cluego [29] and Metascape [30] analysis software (Fig. S1). The P23H rod dataset revealed an increased level of markers for oxidative stress (sirt2, atf4a, and oxr1b), response to misfolded proteins (hsp70.3, hspa5, dnajb1b, and xbp1), actin de-polymerization (eps8 and sptbn1) (Fig. 3C), and lipid phosphorylation (dgkzb, pik3ca, and pik3r3a) (Fig. 3D), reflecting physiologic responses to the P23H mutant rhodopsin that ultimately lead to photoreceptor degeneration. Interestingly we have also seen the increased transcriptome of certain ciliary proteins (cep290, bbs5, rpgrip1, and rp1) and calcium channels (cacna1c and tmco1) (Fig. 3D).

Fig. 3

Transcriptomic analysis of rods and newly formed rods in the WT-P23H integrated dataset reveal an increased response to oxidative stress and misfolded proteins. A Overlap of UMAP dimensionality reduction and visualization of rod transcriptomes from WT and P23H transgenic fish. B Violin plots showing levels of expression of canonical marker genes in WT and P23H rods. C–E Violin plots showing expression levels of differentially expressed genes (DEGs) involved in various functional pathways including oxidative stress, response to misfolded proteins, ciliary transport, synaptic remodeling, and metabolism in WT and P23H rods. F Dot plot showing the expression level and the percentage of cells expressing DEGs involved in different functional pathways in WT and P23H. G Violin plots showing expression levels of DEGs involved in proteasome and ubiquitin-mediated degradation pathway and stress response markers in WT and P23H newly formed rods

Metabolism-related differentially expressed genes (DEGs) were upregulated in the P23H rods including oxidative phosphorylation (sdha, sdhb, and ndufv1), glycolysis (pgam1b, eno3, pfkfb4b, and slc3a2b) and TCA cycle (cs and ndufv1) (Fig. 3E). It is possible that the stress due to misfolded protein response leads to increased energy demand that is met by increased oxidative metabolism. Notably, nxnl1, the zebrafish homolog of Rod-derived Cone Viability factor, is strongly decreased in P23H rods (Fig. 3E). This factor promotes glucose uptake in cones [32] and potentially rods [33]; its loss could result in reduced glucose uptake forcing a switch from oxidative glycolysis to oxidative phosphorylation to supply cellular energy demand. Our pathway analysis depicts an increase in pathways involved in the synthesis of precursor metabolites, glycolysis, pyruvate metabolism, and TCA cycle to cope with increased stress and energy demand in the P23H rods (Fig S1). It has been shown recently by a multi-omic study that early metabolic imbalance and mitochondrial stress in neonatal photoreceptors lead to cell death in the Pde6b rd1 mouse model of retinal degeneration [34]. Our results correlate with the metabolic changes reported in that study.

Importantly, we noticed an increase in gpx4b and ppifb, which are involved in the cellular response to oxidative damage, suggesting that rods are experiencing oxidative stress (Fig. 3F). DEGs involved in ciliary transport including tulp1b and arl13a are also increased in P23H. Mutations in tulp1 are reported to contribute to  ~ 5% of total RP cases [35] and mutations in arl13 lead to the failure of rod outer segment formation [36]. Interestingly, we noted a significant increase in the transcript level of rod arrestins (saga, sagb) in P23H compared to the WT (Fig. 3B, F). The arrestins may be upregulated to compensate for the constitutive activation of P23H mutant rhodopsin, as it has been shown previously that in the mammalian retina stable complexes are formed between mutant rhodopsin and rod arrestin leading to rod cell death [37]. It is implicated that tulp1b is responsible for the ciliary transport of arrestins [38] and the increased arrestins may accumulate in cilia leading to cell death. DEGs involved in cyclic nucleotide-gated channels including cnga1 and cngb1a are both reduced in P23H. Both are involved in the photo-transduction cascade, and their decrease suggests a downregulation of photo-transduction. The DEGs involved in the regulation of synaptic plasticity (atf4a and rims2b) are upregulated in the P23H, suggesting ongoing remodeling of synapses in P23H rods (Fig. 3F). Previous studies in zebrafish models of rod damage do not reveal such comprehensive transcriptional changes [39,40,41] and it is important to understand these changes at the molecular level to target therapeutic interventions appropriately.

A newly forming rod cell cluster is evident in P23H

We identified cluster 13 as newly forming rods, as this cluster expresses both progenitor cell markers and rod photoreceptor markers. The progenitor markers including stmn1a, ccnd1, and hes6 are present at a low level, whereas the rod differentiation markers including fabp7a, cxxc5a, arl13a, and neurod1 are present at high levels in the transcriptome, representing the shift of cells from progenitor state to a differentiated state. This cluster of newly forming rods was not evident in the WT dataset when clustered alone (Fig. 1A), and can be only seen in WT after integrating with P23H for the cluster analysis (Fig. 2A, B). In keeping with the persistent loss and regeneration of rods in the P23H animals, there was a larger number of newly forming rods in P23H (3.8% of cells in the dataset) compared to the WT (1.1% of cells). Functional analysis with ClueGo and Metascape showed that pathways involved in the proteasome and ubiquitin-mediated degradation (UPR), eukaryotic translation, rod cell development, and eye morphogenesis are enhanced in the P23H transcriptome (Fig S1C). The transcripts of genes involved in proteasome and ubiquitin-mediated degradation including psma1, psma8, psmc2, psmc5, psmd3, psmd7, and stress response markers including ciarta, cry1aa and hsf2 are more highly expressed in the P23H dataset compared to WT (Fig. 3G), suggesting that even newly forming rods are experiencing stress associated with expression of the P23H mutant rhodopsin.

Cone photoreceptors show markers of stress and changes in metabolic and circadian genes

We initially identified two cone clusters, cluster 2 and cluster 6. The short-wavelength opsin opn1sw1 is seen in cluster 6 and the long-wavelength opsin opn1lw1 is enriched in cluster 2 (Fig. 4A). While expression of the P23H mutant rhodopsin is limited to the rods [26], cones in the P23H retina show an increase in the transcriptome of hsf2-associated stress response pathway (Fig. 4B), suggesting that cones in this model experience elevated physiological stress. There may be increased oxidative stress in cones due to increased exposure to oxygen in the absence of fully functioning rods, which may lead to increased hsf2 and subsequent activation of the MAPK signaling pathway [42]. Our data further confirmed the increased expression of mapk6 in P23H cones compared to WT (Fig. 4B). Cones in the P23H retina also show increased expression of genes for glycolysis and the TCA cycle (pfkfb3, pfkfb4b, eno3, cs, aco2, and got1) (Fig. 4B), but like rods, expression of nxnl1 is reduced in cones. These results imply that cones rely on aerobic metabolism, but that glucose supply could be limiting. It has been shown previously that aerobic glycolysis in photoreceptors is essential for normal rod function and is a metabolic choice that augments cone function and survival during nutrient-stress conditions [43]. We also noticed increased expression in DEGs involved in ribose phosphate biosynthesis (socs3s and gcdha) and RNA polymerase II transcription and elongation (cdk7, polr2a, and gatad2ab) (Fig. 4B). Sufficient supplies of NADPH, ATP, and the metabolic intermediates ensure rapid macromolecular synthesis underlying the continuous self-renewal of cone OS and the demand may be increased to cope with the increased stress.

Fig. 4

Comparative transcriptomic analysis of cones in WT-P23H integrated dataset reveals changes in metabolism and circadian regulation. A Violin plots showing levels of expression of short- and long-wavelength opsins in cone clusters 2 and 6. B Dot plot showing the expression level and the percentage of cells expressing DEGs involved in different metabolic functional pathways between WT and P23H. C Dot plot showing the expression level and the percentage of cells expressing DEGs involved in circadian regulation between WT and P23H. D Hiplex fluorescent in situ hybridization for tefa, per1, and per2 in sections of adult retina. The fluorescent product is shown in the left third of each image, while the remainder displays a binary mask in the fluorescent in situ hybridization product channel. All three genes show increased expression throughout the retina in P23H compared to WT, reflecting the changes captured by single-cell transcriptome analysis. Scale bar applies to all panels

The hsf2-associated stress response pathway includes increased levels of cry1a and nr1d2b, which are shown to be involved in circadian function [44]. It has been shown that in zebrafish larvae, the expression of hsf2 was concomitant with the expression of genes involved in the response to oxidative stress and chaperone genes, and it occurs in phase with cry1a and other genes that belong to the negative arm of the transcriptional regulation network of the circadian rhythm [45]. We examined the transcriptome profile of other molecular circadian regulators and found that the circadian regulators cry1aa, cry1bb, per1b, per2, and tefa were increased in the P23H cones compared to the WT (Fig. 4C). To validate the single-cell transcriptome data, we examined the expression of per1, per2, and tefa using multiplex in situ hybridization. The results demonstrate the increased expression of per1, per2, and tefa in the P23H retina compared to the WT (Fig. 4D). Increased levels of expression of all three genes were detected in both the outer nuclear layer (ONL), harboring photoreceptors, and in the inner nuclear layer (INL), occupied by numerous other cell types. The overall increased expression of these circadian cycle genes in P23H tissue samples compared to WT tissue samples is consistent with our single-cell data. Disruption of behavioral and retinal circadian rhythms has recently been reported in a mouse model of RP due to an RNA splicing factor mutation [46], suggesting commonality in pathology. Overall, these changes suggest that the cones undergo comprehensive modifications including metabolism and circadian activity during the degeneration/regeneration scenario.

RPE tight junction genes and stress-protective mechanisms are lowered in the RP model

The zebrafish retinal transcriptome showed two different RPE clusters, clusters 10 and 15 (Fig. 5A). Both RPE cell clusters express the canonical RPE genes rpe65a, dct and lrp1aa (Fig. 5B) but are segregated due to the presence of transcriptome that enriches for different functional pathways (Fig S2). Differential expression analysis revealed that cluster 15 has many unique transcripts that were not expressed in cluster 10 (Fig. 5C). The 5 most up-regulated genes in RPE cluster 10 were RNA-binding proteins including rpl3, rpl4, rpl5a, and rpl6, which enable protein translation; these were also present in cluster 15. Functional pathways involved in phagocytosis, endocytosis, lysosome, oxidative stress and drug metabolic process are all enriched in cluster 15, whereas cluster 10 showed pathways involved in translation and ribosome assembly (Fig S2). The receptor tyrosine kinase mertka involved in RPE phagocytosis of photoreceptor outer segments (POS) [47] is highly expressed in cluster 15, suggesting active phagocytosis in these RPE cells. RPE cells have been categorized into multiple subpopulations based on their location and morphology [48,49,50,51]. While our data do not have the spatial resolution to differentiate RPE subpopulations by location and morphology, the 2 subpopulations of RPE cells in our data are functionally consistent with subpopulations found in previous studies based on transcriptomic functional activities [52, 53].

Fig. 5

RPE transcriptome is altered in the RP model. A tSNE map showing the proportions of RPE cells found in the P23H and WT retinal samples. B Violin plots showing levels of expression of canonical RPE marker genes in WT and P23H RPE cells. C Violin plots showing expression levels of DEGs differentiating RPE clusters 10 and 15. D Violin plots showing expression levels of genes involved in tight junctions in WT and P23H RPE. E Immuno-fluorescent labeling of tight junction protein ZO1 in retina of adult P23H and WT zebrafish. The layer of retinal pigmented epithelium (RPE) is at the top. Additional strong ZO1 labeling below the photoreceptor inner segments (PR) is the outer limiting membrane formed by Müller cell end-feet tight junctions. Right panels show magnified view of ZO1 labeling in the RPE. Scale bars in WT apply to both WT and P23H. F Violin plots showing levels of expression of DEGs including those involved in circadian regulation, POS phagocytosis and retinoid signaling in WT and P23H RPE cells

We have seen a drastic decrease in the number of RPE cells in the P23H RP model compared to the WT (Fig. 5A, E), which is consistent with several previously reported studies that when the rods are lost oxygen continues to flow into the outer retina from the choriocapillaris, creating a hyperoxic environment that is presumably hostile for the remaining cells [54, 55]. In some retinal degeneration models, disruption of RPE tight junctions is observed, affecting the integrity of RPE cells [17, 56]. Analysis of the DEGs in RPE cells involved in tight junction formation, including tjp1a, tjp1b, ctnna1, and ctnnb1, revealed a decrease in the P23H compared to the WT (Fig. 5D). Immunostaining for tight junction protein ZO-1 confirmed that tight junctions in the P23H RPE were disrupted when compared to WT (Fig. 5E). In some areas of the P23H retina, the tight junctions were almost completely absent (Fig. 5E). This suggests that RPE tight junction disruption is a pathological effect in RP that is conserved across species. We also note that the tight junction pathology and reduced photoreceptor outer segment lengths may have influenced the number of RPE cells captured during preparation of the single-cell libraries, potentially contributing to the reduced number of RPE cells found in P23H fish.

Ongoing degeneration and regeneration in the P23H RP model had a significant impact on several other aspects of the RPE transcriptome. We have noticed an increase in circadian genes, such as per2, per1b, and cry1aa, in the P23H dataset (Fig. 5F), similar to that seen in cones, suggesting the change is coordinated between RPE and photoreceptors. It has been shown in previous studies that per2 regulates POS phagocytosis in RPE cells [57], suggesting a change in the rhythm of POS phagocytosis in the RP model. The phosphoinositide 3-kinase regulators pik3r2, pik3r3a, and socs3a are reduced in the P23H retina compared to the WT (Fig. 5F). It has been shown in previous studies that the activation of the PI3K-Akt pathway protects RPE cells against the deleterious effects of oxidative stress that occur due to POS phagocytosis [58, 59]. These protective mechanisms seem to be reduced in the RPE of P23H zebrafish.

The retinoid X receptor rxrga and nuclear receptors nr4a1 and nr4a3 are lowered in the P23H dataset compared to the WT (Fig. 5F). NR4A receptors heterodimerize with RXR receptors and activate transcription in a 9-cis retinoic acid-dependent manner. They can also repress inflammatory gene promoters by recruiting corepressor complexes [60]. This complex suite of transcriptome changes suggests that the RPE cells lose several stress-protective mechanisms during the degeneration–regeneration scenario.

Rod bipolar cells show evidence of stress and neuronal remodeling

We have identified three different bipolar cell clusters, 0, 20 and 23 in our dataset; each cluster expressed a specific set of markers (Fig. 6A). The rod bipolar cell marker prkca [22] is highly enriched in cluster 20 (avg. expression 2.7) and to a minor extent in cluster 0 (avg. expression 0.7) only in very few cells (Fig. 6B). Because of the high expression of prkca, we considered cluster 20 likely to be the rod-dominated mixed bipolar cell [61] that forms most synaptic contacts with rods. All of the bipolar cell clusters displayed differences in transcriptome profile between the WT and P23H (Fig. 6B), including downregulation of neuronal specification genes neurod4 and nrxn3b and upregulation of stress-related genes including ubl3a. A prominent change in all bipolar cell clusters is the dramatic upregulation of rgs16 (Fig. 6B), which accelerates the offset of G protein signaling [62, 63], potentially modifying signaling by metabotropic glutamate receptors, dopamine receptors, or other GPCRs.

Fig. 6

Comparative transcriptomic analysis of bipolar cell clusters in the WT-P23H integrated dataset shows evidence of stress and neuronal remodeling in the RP model. A Violin plots showing levels of expression of specific marker genes in different clusters of bipolar cells (0, 20, and 23). B Violin plots showing changes in expression of genes involved in differentiation and stress response in the three bipolar cell clusters. C Dot plot showing the expression level and the percentage of cells expressing DEGs involved in axon guidance, stress response and circadian regulation between WT and P23H in cluster 20

We focused on cluster 20 to examine the influence of rod degeneration and regeneration on gene expression. Genes involved in synaptic vesicle trafficking including snap25a, snap25b and nsfa are increased in the P23H dataset compared to the WT (Fig. 6C). This suggests that there may be compensation at the bipolar cell output synapses for the reduced rod input in P23H retina. We also detected increases in the transcriptome of some genes involved in the stress response including dnajb1b, hspd1, sod1, and gsto1 in the P23H dataset compared to the WT (Fig. 6C). As we found in cones and RPE, levels of per1b, per2, and tefa, which are involved in circadian regulation, were elevated (Fig. 6C). This suggests that rod photoreceptor loss leads to circadian changes that are reflected throughout the retina, consistent with the increased expression detected throughout the retina in the in situ hybridization experiments (Fig. 4D).

The overall comparison between WT rod bipolar (RBP) transcriptome and P23H RBP transcriptome with STRING, which infers protein–protein interaction networks, suggests that the interaction network changes in the P23H RBP cells when compared to WT. P23H RBP display far fewer annotated interactions, and the top enriched pathways are different between WT and P23H (Table S5). Pathways for ATP synthesis and reactive oxygen species (ROS), and reactive nitrogen species (RNS) production were enriched in the P23H RBP transcriptome, suggesting increased energy demand and stress response in the P23H RBP cells (Table S5). Thus, physiological remodeling in the P23H RBP cells involves a complex interplay of oxidative stress, synaptic remodeling and synaptic transmission.

Horizontal and amacrine cells show changes in genes regulating axonal remodeling in the P23H RP model

We identified the horizontal cells using the markers rprmb, aqp9a and mdka (Fig. 7A). The gap junction proteins cx52.6, cx52.9, and cx55.5 are all uniquely enriched in the horizontal cells, serving as additional marker genes (Fig. 7B), and showed no significant differences between WT and P23H retinas. We studied the differences between the horizontal cells in WT and P23H by functional analysis (Fig S3) and identified increased expression of genes involved in axon remodeling and microtubule stability including apc2, camsap1b, sptbn1, sqstm1, ndel1a, and zgc:92,606 (Fig. 7A). Previous studies have shown that apc2 and camsap1 play an essential role in axonal projections through the regulation of microtubule stability [64,