Single-cell RNA sequencing of the retina in a model of retinitis pigmentosa reveals early responses to degeneration in rods and cones
BMC Biology volume 20, Article number: 86 (2022)
In inherited retinal disorders such as retinitis pigmentosa (RP), rod photoreceptor-specific mutations cause primary rod degeneration that is followed by secondary cone death and loss of high-acuity vision. Mechanistic studies of retinal degeneration are challenging because of retinal heterogeneity. Moreover, the detection of early cone responses to rod death is especially difficult due to the paucity of cones in the retina. To resolve heterogeneity in the degenerating retina and investigate events in both types of photoreceptors during primary rod degeneration, we utilized droplet-based single-cell RNA sequencing in an RP mouse model, rd10.
Using trajectory analysis, we defined two consecutive phases of rod degeneration at P21, characterized by the early transient upregulation of Egr1 and the later induction of Cebpd. EGR1 was the transcription factor most significantly associated with the promoters of differentially regulated genes in Egr1-positive rods in silico. Silencing Egr1 affected the expression levels of two of these genes in vitro. Degenerating rods exhibited changes associated with metabolism, neuroprotection, and modifications to synapses and microtubules. Egr1 was also the most strongly upregulated transcript in cones. Its upregulation in cones accompanied potential early respiratory dysfunction and changes in signaling pathways. The expression pattern of EGR1 in the retina was dynamic during degeneration, with a transient increase of EGR1 immunoreactivity in both rods and cones during the early stages of their degenerative processes.
Our results identify early and late changes in degenerating rd10 rod photoreceptors and reveal early responses to rod degeneration in cones not expressing the disease-causing mutation, pointing to mechanisms relevant for secondary cone degeneration. In addition, our data implicate EGR1 as a potential key regulator of early degenerative events in rods and cones, providing a potential broad target for modulating photoreceptor degeneration.
The degeneration of light-sensitive rod and cone photoreceptor cells underlies various forms of inherited retinal disorders (IRDs), including retinitis pigmentosa (RP), the primary cause of blindness among children and working-age adults [1, 2]. Over 100 genomic loci are implicated in RP , a progressive and heterogeneous group of disorders often characterized by night blindness due to primary rod photoreceptor death, followed by the loss of high-acuity color vision upon a secondary wave of cone photoreceptor death, leading to complete blindness .
Due to the innate heterogeneity of the retina, animal models that recapitulate aspects of the disease while preserving the differences among retinal cells are essential for studying IRDs. The rd10 mouse is an RP model frequently used for therapeutic studies since it exhibits a later onset and slower progression of retinal degeneration than other similar models (e.g., rd1), resembling the human autosomal recessive form of the disorder more closely .
The rd10 mouse harbors a spontaneous, autosomal missense mutation in the rod-specific beta-subunit of the phosphodiesterase gene (Pde6b). PDE6 hydrolyses cyclic GMP (cGMP) upon absorption of light by rhodopsin, decreasing the open state probability of cGMP-gated Na+ and Ca2+ channels and hyperpolarizing the photoreceptor cell. The rd10 mutation causes the destabilization and mislocalization of PDE6B in rods, causing an increase in intracellular cGMP and Ca2+ levels . This leads to a primary wave of rod photoreceptor degeneration peaking around postnatal day (P) 21 [6, 7]. Even though cones do not express Pde6b, they are indirectly affected by rod degeneration and die in a secondary wave of cone death starting at around P40 [7, 8]. Various consequences of rod death may contribute to the secondary degeneration of cones. The loss of trophic factors secreted by healthy rods, such as rod-derived cone viability factor , may lead to an imbalance of glucose metabolism in cones , resulting in a slow mode of cell death . In addition, the loss of oxygen-consuming rods in the degenerating retina may impose hyperoxia on cones, causing increased oxidative stress . Despite the importance of cones to vision, studies of cone-specific responses to rod degeneration are hindered by the relatively low numbers of cones in the retina [13, 14]. Therefore, it remains unclear which early effects are induced in cones at the peak of rod degeneration, and how these may initiate or contribute to the events leading to their eventual degeneration.
Transcriptomic analyses of RP mouse models via microarrays [15, 16] and whole-retina RNA sequencing [17, 18] have been instrumental in furthering our understanding of RP, especially regarding the pathophysiological mechanisms involved in retinal degeneration. However, such techniques cannot differentiate between individual cells, potentially missing critical transcriptional changes restricted to specific cell types or subpopulations of cells in different stages of degeneration. Therefore, techniques with higher resolution are needed to overcome cellular heterogeneity. Recent advances in droplet-based single-cell RNA sequencing (scRNAseq) have enabled the detailed characterization of retinal cell types in mouse  and human retinas and organoids [19,20,21], as well as the elucidation of signaling events and critical factors controlling injury responses in microglia  and Müller glia , in addition to those regulating distinct phases of retinal development  at a single-cell resolution.
Here, we analyzed the degenerative process in the rd10 mouse retina at the single-cell level. Using the 10X Genomics platform, we investigated transcriptomic changes in rd10 rods with progressing degeneration using trajectory analysis. We also identified cone responses to rod degeneration at P21, which involves the long-term activation of EGR1 (early growth response 1), an immediate-early gene.
Droplet sequencing of rd10 and wild type retinas
To obtain transcriptomes of individual cells during retinal degeneration, we collected retinas from rd10 and C57BL/6 wild-type control mice at P21, a time point that falls within the primary wave of rod photoreceptor degeneration in rd10 mice . A male and a female littermate were sequenced for each genotype to ensure equal representation of sexes in the analysis. To avoid potential confounding effects of optical coherence tomography (OCT) on retinal transcriptomes, we confirmed the increased reflectivity of the degenerating outer retina in rd10 animals with OCT in littermates not used for sequencing (Additional file 1: Fig. S1).
After alignment and demultiplexing with the Cell Ranger pipeline , quality control with the Seurat package , and a marker gene co-expression-based doublet filtering algorithm , we retained 27,374 cells for downstream analyses. Using shared-nearest neighbor embedding and modularity optimization, followed by additional quality control measures (see Methods and Additional file 2: Fig. S2), we filtered our dataset to 25’484 single cells grouped into 31 clusters (Fig. 1A), with 1169 median genes and 2001.5 median unique molecular identifiers (UMIs) per cell. Clusters were annotated using known marker genes [23, 28] (Fig. 1B, Additional file 3). We identified microglia, Müller glia, and all types of retinal neurons except ganglion cells (Fig. 1A, B), which might not have tolerated the dissociation and capture protocols used in our experiments.
While clusters generally included cells from both wild-type and rd10 retinas, three rod clusters (C01-03) contained almost exclusively either wild-type (C01) or rd10 (C02-03) cells (Fig. 1C, Table 1). The specificity of the effect to rod clusters excluded a batch effect and instead suggested that wild type and rd10 rods had distinct transcriptomic profiles (see below and Methods).
Transcriptomic analysis of rod photoreceptors points to Egr1 as a critical regulator of early degenerative events
Preceding note: to increase readability, we do not define all gene abbreviations and refer to Ensembl  for detailed information on genes of interest.
Rod degeneration is the primary event in most IRDs and leads to the secondary loss of cones. For an in-depth analysis of rod transcriptomes during their degenerative process, we analyzed the 18,268 rod photoreceptors in clusters C01-C03 (Fig. 2A). Rod degeneration in the rd10 model exhibits a slow progression [4, 6, 8] with a subtle center-to-periphery gradient , suggesting heterogeneity among degenerating rd10 rods in the same retina, with cells simultaneously at different stages of degeneration. To resolve this heterogeneity in our rod dataset, we used the Slingshot package, which constructs one-dimensional trajectories that connect cells at different stages of continuous processes (“pseudotime”) . We retrieved a single lineage that started from the predominantly wild-type cluster C01, followed sequentially by the predominantly rd10 clusters C02 and C03 (Fig. 2B, C), with significantly different pseudotime values between the clusters (Fig. 2C). This trajectory indicated that rd10 clusters C02 and C03 comprised rods at different stages of degeneration, with C03 cells further along the process (“late degeneration”) than C02 cells (“early degeneration”).
To study potential transcriptomic changes during rod degeneration, we restricted our gene set by performing pairwise differential expression (DE) analysis between clusters adjacent in pseudotime: C01-C02 (“Early DE analysis”: Additional file 4) and C02-C03 (“Late DE analysis”: Additional file 5). We then investigated the pseudotemporal behavior of all significant DE transcripts from the two analyses in binned pseudotime (Fig. 2D—bar graph) and defined four major DE patterns with k-means clustering : (1) downregulation, primarily at the onset of C01-C02 transition (potential early degeneration); (2) early permanent upregulation, in which expression was induced early and remained elevated; (3) early transient upregulation, in which expression levels peaked during early pseudotime and returned to near baseline later; and (4) late upregulation, in which expression peaked late in pseudotime (potential late degeneration) (Fig. 2D—heatmap, Additional file 6).
At the potential onset of degeneration, 155 transcripts were significantly upregulated and 72 were downregulated in C02 rd10 rods compared to wild-type C01 rods (early DE analysis, |logFC| > 0.25, adjusted p-value < 0.05, Additional file 4). Various transcripts linked to cellular stress were induced in this stage (Fig. 2D, pseudotime bins 8-12), which was further characterized by the activation of transcription factors that belong to the class of immediate-early genes (Egr1, Nr4a1, Ier2). Egr1 was the most highly (but transiently) upregulated transcript in early DE analysis (Fig. 2D, Additional file 4). We noted the transient upregulation of several factors potentially involved in the regulation of water (Aqp1), ion (Kcnb1, Slc17a7, Slc12a5 [permanently upregulated]) and, notably, lipid/retinoid metabolism (Abca4, Sdc4) and homeostasis (Fig. 2D, Additional file 4). Importantly, both the Ca2+ extruder PMCA1 (Atp2b1) expressed in mouse photoreceptors  as well as the calcium signaling mediators Camk2g and Pcp4 were upregulated early and remained high, indicating an attempt to counteract the increase in intracellular Ca2+ levels due to the rd10 mutation  (Fig. 2D, Additional file 4). Furthermore, signaling mediators like Akap11 and Mapk8 and several transcripts associated with microtubules (Tubb5b, Arhgef2, Agtpbp1), synaptic structure and transmission (Cpe, Snap25, Slc1a2, Slc17a7, Slc12a5) as well as anti-apoptotic roles (Gadd45b , Krt18 ) were upregulated in early degeneration (Fig. 2D, Additional file 4).
Meanwhile, early degenerative rd10 rods exhibited downregulation of the linker histone H1 variant Hist1h1c, followed by Hspa1a (HSP70-1) and Hspa1b (HSP70-2) (Additional file 4). HIST1H1C has been associated with increased stress and autophagy in diabetic mouse retinas , while members of the heat shock protein 70 family are involved in retinal neuroprotection . Similar to models of induced retinal degeneration , multiple transcripts involved in phototransduction (Pde6g, Cnga1-b1, Guca1a-1b) decreased after the onset of degeneration in rd10 rods (Fig. 2D, Additional file 4). In addition, we noted the early downregulation of rod-derived cone viability factor Nxnl1  and transcripts encoding glycolytic enzymes, particularly hexokinase 2 (Hk2), which mediates aerobic glycolysis in photoreceptors and is involved in survival upon stress and aging [39, 40], and enolase-1 (Eno1). This suggests an early functional and metabolic disruption in rd10 rods.
Further along the degenerative process, C03 rods (Fig. 2C) had 78 transcripts upregulated and 39 downregulated compared to C02 rods (late DE analysis, |logFC| > 0.25, adjusted p-value < 0.05, Additional file 5). In these cells, the most highly upregulated transcript was Cebpd (Additional file 5), a stress-response transcription factor that was shown to be induced in the rd10 retina  but whose cell-type localization was not described. Another strongly regulated transcription factor was Junb, an AP-1 component found at increased levels in the light-damaged retina  (Additional file 5). These changes were paralleled by a decrease in Egr1 (Fig. 2D, bins 13-16). Many transcripts associated with survival, specifically with leukemia inhibitory factor-mediated endogenous retinal neuroprotection in models of inherited  and induced  degeneration (Fgf2, Edn2, Stat3, Socs3), were strongly upregulated in these cells (Additional file 5). Concomitantly, levels of transcripts with potential anti-apoptotic functions (Gadd45b , Clu , Cerkl ) (Fig. 2D), were increased (Additional file 5). Furthermore, we observed the continued increase in tubulin deglutamylase Agtpbp1 and the upregulation of Arl6, which is critical in protein trafficking to the connecting cilium . Together with the early upregulation of microtubule-associated transcripts, this may indicate progressive changes to the photoreceptor outer segments (Additional file 4). Potential changes also to synaptic structures were suggested by the upregulation of Rock2 [47, 48] and Agrn, as well as a decrease in VGLUT1 (Slc17a7) (Additional file 5).
In late degeneration, the levels of phototransduction mediators (e.g., Rho, Gnb1, Gnat1) decreased further in rd10 rods (Fig. 2D, Additional file 5) and the potential metabolic disruption observed already in early degeneration was followed by the significant downregulation of genes involved in oxidative phosphorylation in late degenerative rods (Additional file 5). Especially affected were the mitochondrial respiratory complex I NADH dehydrogenase (mt-Nd1-5), a crucial subunit of complex III coenzyme Q-cytochrome c reductase (mt-Cytb) , and complex IV cytochrome c oxidase (mt-Co1-3), pointing to mitochondrial respiratory dysfunction at the late degenerative stage (Fig. 2D, Additional file 5).
Several genes connected to retinal degeneration in humans were differentially regulated in degenerating mouse rods, including the downregulation of Fam161a (RP ) and the upregulation of Cerkl (RP ), Arl6 (Bardet-Biedl syndrome ), and Samd7 (RP ), a CRX-regulated transcriptional repressor  which is essential for rod photoreceptor identity  and a paralog of the RP disease gene SAMD11 . Differential regulation of these genes in mouse rods may additionally validate the rd10 mouse as a valuable RP model.
In addition to pairwise DE analyses, we utilized the tradeSeq package  to test the association of gene expression changes to pseudotime. All but 7 significant DE transcripts as determined by pairwise DE analyses had a statistically significant correlation to pseudotime (Additional file 7), further supporting the relationship of these genes to potential damage progression.
Since the most highly upregulated transcripts in the two distinct stages of rod degeneration encode for transcription factors EGR1 (early) and CEBPD (late), we asked whether they could be the main drivers of transcriptomic changes within their associated stages. To explore this possibility, we performed an in silico promoter scan using Enrichr (TRANSFAC and JASPAR PWMs [position weight matrices])  and searched for exact matches to known transcription factor motifs in a − 2000:+ 500 base range of the transcriptional start site of significantly differentially expressed genes (see above). EGR1 was the transcription factor most strongly associated with the promoters of early degeneration DE genes, based on the Enrichr combined score (Additional file 8). No significant association could be made for late degeneration DE genes (Additional file 9). To validate the in silico analysis, we directly tested whether EGR1 is involved in regulating the expression of Gadd45b and Agtpbp1, two of the genes predicted to be EGR1 targets (Additional file 8). Using small interfering RNA (siRNA), we downregulated the expression of Egr1 in photoreceptor-like 661W cells . While the expression of Apoe, a gene not predicted to be a target for EGR1, was not affected, the levels of Gadd45b and Agtpbp1 increased significantly upon siRNA-mediated knockdown of Egr1 (Additional file 10: Fig. S3), suggesting that EGR1 might function as a repressor of these two genes, an activity that was ascribed to EGR1 already by others [59, 60]. This potential regulation is particularly interesting, as we found the levels of both transcripts, especially Agtpbp1, to increase further in rods with progressing degeneration, when Egr1 expression was lost (Fig. 2D).
The cone response to rod degeneration at P21
Even though the mutant Pde6b gene that leads to rod degeneration in rd10 mice is not expressed in cones, cones respond to the loss of rods and degenerate in a second wave of cell death. To address early mechanisms in cones that may lead to this secondary cone cell death, we analyzed the transcriptomes of the 1190 cone cells in cluster C04 (Fig. 1A, B). Differential gene expression analysis between rd10 and wild type cones in this cluster (Fig. 3A) returned 245 upregulated and 87 downregulated transcripts (|logFC| > 0.25, adjusted p-value < 0.05; Additional file 11) and revealed that Egr1 was the most highly upregulated gene in rd10 cones (Fig. 3B, C), similar to rd10 rods early in degeneration (Fig. 2D).
Among upregulated transcripts in rd10 cones, we noted other immediate-early genes (Nr4a1, Junb, Ier2) and various signaling mediators (e.g., Inhbb, Bambi [TGF-beta signaling], Tnfrsf12a [TNF signaling], Ctnnbip1 [beta-catenin signaling], Map3k8) (Fig. 3C, Additional file 11). Since glucocorticoid receptor activation is neuroprotective against light damage , it is interesting that glucocorticoid signaling regulators (Ncoa6 , Calr [63, 64], Arid1a ) were upregulated in cones (Additional file 11). In addition, ER stress and unfolded protein response (Foxred2, Xbp1, Ddit3, Eif2ak3, Pdia6, Hsp90b1) as well as autophagy-related transcripts (Gabarap, Gabarapl1, and a component of the BRAF35-HDAC complex [Gse1]) were expressed at higher levels in rd10 cones. Similar to rods, the levels of several tubulin cytoskeleton-associated transcripts (including Tubb5, Tubb4b, and Agtpbp1) were also increased in rd10 cones. We also observed a notable increase in Capns1, which encodes the regulatory calpain subunit, in rd10 cones but not rods (Additional file 11). Calpain activity may be involved in photoreceptor degeneration in multiple RP models .
Several genes involved in phototransduction (Gnat2, Opn1sw, Opn1mw, Guca1b, Pde6g, Pdc) were downregulated in rd10 cones prior to the onset of their degeneration (Fig. 3C; Additional file 11). This observation is in line with our previous report showing a transient downregulation of cone-specific genes early after light exposure, which primarily damages rods . Interestingly, transcripts of the ‘rod-specific’ Pde6g gene were detected in cones, an observation already made by others in various mouse models [66,67,68]. Aside from genes required for phototransduction, rod degeneration significantly affected mitochondria in neighboring cones. As in degenerating rods, various mitochondrially encoded genes involved in oxidative phosphorylation were downregulated in rd10 cones (Fig. 3C, Additional file 11). Among those were multiple subunits of the respiratory complex I (mt-Nd1-5, Ndufaf5, Ndufa10), mt-Cytb, all mitochondrially encoded subunits of the respiratory complex IV (mt-Co1-3), subunits 6 and 8 of the ATP synthase Fo region (mt-Atp6, mt-Atp8) as well as the nuclear-encoded hexokinase 2 (Hk2) (Additional file 11).
In general, we found an extensive overlap of differentially regulated genes between rods and cones. Among 332 statistically significant DE genes in cones, 118 were also DE among rod clusters, including Egr1, Krt18, Nr4a1, Junb, Agtpbp1, Fgf2, Hspa1a, Hspa1b, Hk2, and various mitochondrial genes (Additional file 11). Correlation analysis between rd10 cones, wild-type cones, and rod clusters confirmed that rd10 cones were transcriptomically more similar to degenerating rods than wild-type cones to wild-type rods (Fig. 3D).
Expression levels for transcripts of interest in rd10 retinas with increasing age
For additional validation of our single-cell transcriptomic findings, we determined the levels of selected transcripts with semi-quantitative real-time PCR (qPCR) in whole retinas of rd10 and control mice at eight time points, from 2 weeks up to 6 months of age (Additional file 12: Fig. S4).
Egr1 and Agtpbp1 were selected as representative transcripts activated in damaged rods and cones. Egr1 reached its peak expression in P21 and remained slightly above controls up to P49. This effect was statistically significant only in P21, in line with a transient activation of Egr1 in rods, the most abundant cell type of the retina (Fig. 2D). Whereas Agtpbp1 levels were similar between rd10 and control retinas at P21, it was significantly upregulated at P28, a time point when most rods are in late stages of degeneration . This correlated well with the strongly increased expression in rd10 rods of cluster C03 that contained cells further along the degenerative process (Fig. 2D). Interestingly, Agtpbp1 levels were significantly lower in rd10 retinas compared to controls at P14, prior to the onset of rod degeneration. This may be linked to pre-degenerative outer segment defects in these retinas .
To test our findings regarding the decrease in mitochondrial transcripts, we analyzed mt-Atp6, mt-Co2, and mt-Nd4 levels in the retina. While they were not strongly regulated when analyzed in whole retinas, we noted significantly lower values for all three genes in rd10 retinas during the phase of strong rod degeneration (P21 for mt-Atp6 and mt-Co2, P28 for mt-Nd4). No significant differences were detected in the expression of Tfam, a key mitochondrial transcriptional factor.
While droplet sequencing analyzes gene expression changes in distinct subpopulations of cells, qPCR averages expression levels over all cells in a tissue, which may explain differences between the two sets of data. Regardless, qPCR recapitulated the sequential upregulation pattern in Egr1 and Agtpbp1, supporting the pseudotime analysis where the peak of Egr1 expression was an early but transient event that preceded the peak of Agtpbp1 expression (Fig. 2D). Since Agtpbp1 was among the predicted target genes of EGR1 (Additional files 8 and 9) and EGR1 may act as a repressor of Agtpbp1 expression (Additional file 10: Fig. S3), it is tempting to speculate that the late upregulation of Agtpbp1 in rd10 rods is directly connected to the downregulation of Egr1 at this point in pseudotime (Fig. 2D).
EGR1 constitutes an early response to degeneration in both rods and cones
Based on our findings that Egr1 mRNA was the most highly upregulated transcript in both early degenerating rods and in cones responding to primary rod degeneration, we hypothesized that EGR1 might be involved in an early stress response shared by both types of photoreceptors. To investigate the localization of EGR1 protein during degeneration, we combined immunostainings against EGR1 with short-wavelength opsin (OPN1SW) or rhodopsin (RHO) to discriminate between EGR1-expressing cones and rods, respectively (Fig. 4), and with glutamine synthetase (GS) to label Müller glia (Additional file 13: Fig. S5A-C).
In wild-type retinas, EGR1 immunoreactivity was confined mainly to the inner nuclear layer (INL) and ganglion cell layer (GCL), with few positive cells in the outer nuclear layer (ONL) (Fig. 4A). Co-immunostaining with OPN1SW and the location of the positive nuclei in the distal ONL identified these cells as cones (Fig. 4B, arrow). Müller glia, identified by GS staining, did not co-localize with the EGR1 signal in the INL (Additional file 13: Fig. S5A). Findings for wild-type neurons were in agreement with scRNAseq, which detected Egr1 mRNA in few cones, subsets of bipolar and amacrine cells, and microglia (Additional file 14: Fig. S6). In contrast, the relatively high levels of Egr1 mRNA detected in wild-type Müller glia did not correlate to immunoreactivity. The successful detection of EGR1 in some Müller glia in rd10 retinas (Additional file 13: Fig. S5B, C) suggested a potential cell-type-specific regulation of EGR1 at the post-transcriptional level [70, 71]. EGR1 localization in wild-type retinas was comparable at all time points tested (Fig. 4A and Additional file 13: Fig. S5A for P21, S5D for P35).
At P21, rd10 retinas displayed a major increase in the number of EGR1-immunopositive cells in the ONL compared to wild-type retinas (Fig. 4C). Since the EGR1 signal was present both in OPN1SW (Fig. 4D, arrows) and RHO (Fig. 4E, triangles) positive cells, this indicates that both rods and cones activated EGR1 in degenerating rd10 retinas. However, some rods with an intense RHO signal in the perikaryon were only faintly positive for EGR1 immunoreactivity (Fig. 4E, asterisks). We additionally observed some Müller glia in rd10 retinas at P21 that were EGR1-positive, including few with their nuclei shifted towards the ONL (Additional file 13: Fig. S5B, C; arrows). A similar pattern of EGR1 localization was observed in rd10 retinas at P25 (Fig. 4F–H).
The pattern of EGR1 immunoreactivity in the ONL changed with progressing degeneration from P28 onward (Fig. 4I, L). Cones remained EGR1-positive (Fig. 4J, M; arrows), whereas many rods now lacked EGR1 signal (Fig. 4K, N; asterisks), in support of our transcriptomic data showing that Egr1 expression in rods decreased with increasing degeneration. At P59, rd10 retinas lost almost all rods , and the ONL mostly contained a single row of cone nuclei. While the number of cells with EGR1 immunoreactivity was generally low at late time points, several cones and few cells of the inner retina remained EGR1-positive (Fig. 4O, P). Some cone cells with no outer segments and OPN1SW staining only in the perikaryon had no EGR1 immunoreactivity (Fig. 4P, arrow). These cells were likely at advanced stages of secondary cone degeneration.
Overall, the number of EGR1-positive cells in the inner retina seemed to decrease with the progression of degeneration (Fig. 4C, F, I, L, O). The dependence of EGR1 localization on the degenerative state of the rd10 retina was especially apparent at later time points such as P35, where we occasionally observed mildly affected peripheral areas with higher numbers of photoreceptors and better conserved outer segments, adjacent to regions with severe degeneration (Fig. 4Q). In mildly affected areas, the number of EGR1-positive cells was high in the inner retina and low in the ONL, whereas the pattern was largely reversed in the strongly affected regions (Fig. 4Q). This pattern might further bolster previous a hypothesis proposing that EGR1 expression is regulated by a complex communication between photoreceptors and inner retinal cells through their synaptic network .
Due to the heterogeneity of the retina and its cell-type-specific degenerative processes, high-resolution approaches such as scRNAseq are invaluable for mechanistic studies of inherited retinal degenerations. Using droplet sequencing, we identified two transcriptomically distinct stages of rod degeneration, characterized by the sequential activation of Egr1 and Cebpd. Furthermore, we identified cone responses to primary rod cell death and showed that Egr1 is also a part of early transcriptomic changes in cones, which include potential respiratory dysfunction.
Rod transcriptomes at the peak of degeneration
In order to resolve the heterogeneity among degenerating rods, we performed trajectory analysis and defined a single pseudotemporal trajectory that originated from wild-type C01 cells, followed by rd10 C02 cells, and ended at rd10 C03 cells. Therefore, we hypothesized that C02 rods represented an early stage of degeneration, whereas C03 rods were further along in the same process.
Early events in rod degeneration
Hallmarks of early degenerative events in rd10 rods included the activation of IEGs, changes in genes associated with water, ion and retinoid homeostasis as well as aerobic glycolysis, a potential induction of signaling events for cellular protection, and the downregulation of genes involved in phototransduction, potentially to save energy and/or prevent oxidative stress.
IEG activation constitutes the first response to a variety of intrinsic and extrinsic stimuli, including stressors. Many are transcription factors that orchestrate cellular responses by regulating downstream targets before they are rapidly inactivated . Egr1 was the most highly upregulated IEG early in rod degeneration (Additional file 5). Egr1 is a downstream target of many signaling pathways , including those relevant for retinal neuroprotection [73, 74] and photoreceptor synaptic stability [48, 75, 76]. As such, Egr1 levels were found to be differentially regulated by diverse forms of retinal stress and injury including photoreceptor degeneration in rd1  and rds mice , in retinoschisis  and after optic nerve crush in rats . In rds mice, Egr1 was detected specifically in rods further supporting a general role for Egr1 in the photoreceptor stress response. While Egr1 expression was upregulated in all of the above models, it was reduced in mouse retinas after hypoxia  and in the Akita mouse model of diabetes . Thus, regulation of Egr1 may depend on the stimulus and, potentially, the cell types affected by the stress.
Our in silico analysis pointed to EGR1 as a potential regulator of genes involved in rod degeneration (Additional file 8), including Gadd45b and Agtpbp1, genes with steadily increasing expression levels with the progression of degeneration (Fig. 2D, Additional file 6). Downregulating Egr1 in photoreceptor-like 661W cells in vitro led to the upregulation of these two transcripts. While counterintuitive at first sight, this result is in line with EGR1’s dual function as both activator and repressor of gene expression, sometimes acting in opposite ways on the same target gene depending on the cellular context . While speculative, EGR1 might assume either activity in retinal cells according to their degenerative state. The late upregulation of Agtpbp1 and Gadd45b (Fig. 2D, Additional file 6) in degenerating rods might thus be directly connected to the reduced expression of Egr1 and thus attenuated repression of some target genes at this point in pseudotime. Together, our results highlight a potential role of EGR1 in the reaction of rods to stress. Importantly, Egr1 was also the highest upregulated transcript in cones. However, while Egr1 expression declined in rods along the trajectory and in whole retinas as degeneration progressed, it remained high in cones that persisted in the rd10 retina for some time after primary rod death (see below).
The activation of IEGs accompanied the induction of transcripts linked to calcium (Camk2g, Pcp4) and apoptosis/neuroprotective signaling (Gadd45b, Krt18, Fgf2, Stat3) as well as to ion, water, and retinoid homeostasis (Atp2b1, Aqp1, Abca4, Kcnb1, Atp2b1) (Additional file 5). Changes in transcripts related to calcium signaling validated the degenerative state of the early rd10 rods, as the rd10 mutation induces a detrimental Ca2+ influx in affected cells . Together, these transcripts point to the activation of defense mechanisms in early rod stress, potentially including the initiation of an anti-apoptotic response, suggested by the early upregulation of Gadd45b and Krt18, which have anti-apoptotic functions in other contexts [34, 35, 81].
Another hallmark of early rod degeneration was the downregulation of transcripts related to phototransduction (Fig. 2D, Additional file 5), potentially to reduce energy consumption and channel available resources towards survival. Such energy-saving strategies may be particularly important since glycolytic enzymes Hk2 and Eno1 were also downregulated in early degenerative rods (Additional file 5), indicating early changes in energy metabolism. In rod-specific Hk2 knockout mice, the resulting metabolic deficiency led to increased oxidative phosphorylation and a decrease in photoreceptor-specific proteins . However, in rd10 rods, oxidative phosphorylation is likely affected with further degeneration (see below), making energy conservation crucial. Consistent with downregulated phototransduction, transcripts that regulate rod photoreceptor fate and expression of phototransduction genes (Nrl , Nr2e3 ; Additional file 5) were also downregulated in early degenerative rods.
Late events in rod degeneration
The late phase of degeneration in rods was characterized by a shift of transcription factors from Egr1 to Cebpd, the remodeling or maintenance of rod cell structures, and the modulation of respiratory functions.
While Egr1 was the most prominently upregulated transcript in early degeneration, Cebpd may be involved in late degenerative events. In lipopolysaccharide-activated microglia in vitro, EGR1 interacts with the Cebpd promoter, suggesting a link between the two transcription factors . The induction of Cebpd was previously detected in the rd10 model  and light damage , and its upregulation is generally associated with neuroprotection [85,86,87] (reviewed in ). Since Cebpd mRNA was shown to remain high in rd10 retinas even after photoreceptor degeneration was complete , CEBPD may act in other cells in advanced degeneration. In this respect, it is notable that CEBPD was shown to be crucial for axon regeneration in the peripheral nervous system . However, the roles CEBPD may play in degenerating photoreceptors are yet to be characterized.
Neural remodeling in the inner retina is a hallmark of both human RP [90,91,92] and its mouse models [92, 93] including rd10 [6, 94], where subtle changes in synaptic transmission are detectable as early as P18 . Synapse-related transcripts (e.g., Rock2, Agrn, Slc17a7 ) were differentially regulated in late degenerative rd10 rods (Additional file 5). RhoA-ROCK signaling mediates photoreceptor axon retraction in retinal detachment [48, 75] and ROCK plays a pro-apoptotic role in the Royal College of Surgeons rat model of RP . Agrin is a heparan sulfate proteoglycan present in developing chick  and mouse retinas  with functions in neurogenesis and synapses [97,98,99]. As EGR1 was shown to repress Agrn expression at the neuromuscular junction , the downregulation of EGR1 in the late phase of rod degeneration may contribute to synaptic changes through releasing Agrn from a repressed state.
Modulation of photoreceptor cellular structures in late degeneration may happen not only in synapses but also in the photoreceptor segments. The photoreceptor outer segment is responsible for phototransduction and is connected to the cell body through a connecting cilium formed by doublet microtubules. Upregulation of transcripts encoding for tubulin proteins and the tubulin deglutamylase NNA1 (Agtpbp1) during late degeneration may indicate structural alterations to the rod outer segment. Defects in rod outer segments were observed already prior to photoreceptor degeneration at around P15 in the rd10 model . The upregulation of the tubulin deglutamylase suggests that post-translational tubulin modifications may be affected in degenerating photoreceptors at P21; genetic defects in these modifications lead to photoreceptor degeneration [101, 102]. Whether tubulin dysregulation also plays a role in photoreceptor degeneration in RP, and whether restoring tubulin modifications may help to ameliorate the degenerative process, remains to be explored. Interestingly, transcripts for tubulins and AGTPBP1 were also upregulated in cones, suggesting that these cells may also remodel their structures or attempt to support the structural integrity of inner and outer segments in response to rod death.
Although metabolic changes started in the early phase (see above), significant downregulation of mitochondrial gene expression (mt-Atp6, mt-Co1-3, mt-Cytb, mt-Nd1-5; Fig. 2D, Additional file 5) in the late phase indicates that respiratory dysfunction develops mainly late during degeneration in rd10 rods. Since this effect primarily involved mitochondrially encoded genes, impairments in mitochondrial DNA stability might be involved. This is supported by the protective effect of metformin in rd10  and rd1  mice. Metformin is an antihyperglycemic guanidine derivative  that activates AMP-activated protein kinase . Applied to the rd10 retina, metformin increased the copy number of mitochondrial DNA and induced expression of mitochondrial genes . Thus, a decrease in mitochondrial activity and respiration may be one of the most significant hallmarks of the late stage of rod death.
Cone responses to rod degeneration
Little is known about how primary rod degeneration affects cones, which do not express the causative mutant gene but nevertheless die in a secondary wave of degeneration. To rescue cone function and vision in patients, it is crucial to identify the early responses of cones to rod degeneration. Three major events may be of significance for cones at this stage. These include increased expression of EGR1, changes in signaling, and reduced mitochondrial gene expression.
Like in the early phase of degeneration in rods, rd10 cones upregulated several IEGs including Egr1, Nr4a1, Junb and Jund, with Egr1 as the highest upregulated transcript. Intriguingly, pseudotime analysis suggested that rods lose Egr1 expression as their degeneration progresses in older rd10 mice. Comparably, rd10 cones exhibited persisting expression of EGR1 until secondary cone degeneration. Together, these findings indicate that EGR1 levels are elevated very early in stressed rods and cones, potentially to protect cellular integrity and function prior to the full activation of cell death mechanisms. Considering the diversity of roles played by EGR1 depending on cellular context , it will be of importance to define the function of EGR1 by artificially modulating its activity specifically in rods or cones to evaluate EGR1 as a treatment target to preserve photoreceptor function and vision.
Furthermore, we observed a striking dynamic regulation of EGR1 protein localization in degenerative retinas. While EGR1 was prominent in the inner retina of wild types, its distribution changed with retinal degeneration in rd10 mice. While photoreceptors became EGR1-positive early in degeneration, cells in the inner retina concomitantly lost EGR1 immunoreactivity. This change of EGR1 localization was especially obvious in sections of rd10 retinas where strongly degenerated regions neighbored less affected areas. The reduced expression of EGR1 in the inner retina of older rd10 mice with advanced photoreceptor degeneration and loss of visual function may indicate that EGR1 expression in the inner retina may depend on signals from photoreceptors. Simon et al. previously showed that EGR1 protein levels increased in the inner retina shortly after the onset of light. Since the increase was independent of circadian gene expression , this suggests a communication between light-activated photoreceptors and cells of the INL. Our results corroborate and extend this finding, by indicating that healthy and functional photoreceptors need to be present to regulate EGR1 expression in the inner retina.
The second emerging pattern in rd10 cones was the differential regulation of transcripts involved in signaling pathways, especially transforming growth factor-beta (TGF-beta). Inhbb and Bambi were two strongly upregulated transcripts that encode proteins associated with TGF-beta signaling. While BAMBI is a silencer of TGF-beta signaling , Inhbb is a member of the TGF-beta superfamily and encodes inhibin beta, which may have diverse, even opposing effects on downstream targets . Inhbb mRNA is also found at higher levels in the P10 Nrl-/- retina, which is devoid of rods , and thus may have a role in cone development. However, the functions of either gene may have in degenerating photoreceptors is unclear.
Perhaps the most important finding related to rd10 cones in our data was the downregulation of multiple transcripts involved in mitochondrial respiration. This is similar to the potential respiratory dysfunction we observed in rd10 rods, and in line with previous work that describe metabolic dysfunction in cones at the onset of cone cell death, when rods have cleared from the retina . Here, we show the reduction in the levels of mitochondria-related transcripts before cone death commences. Interestingly, Nxnl1 was downregulated in early degenerating rd10 rods, suggesting that cones may already experience impaired glucose metabolism  at P21. Metabolic changes in cones are further suggested by the decreased expression of Hk2, a key enzyme for aerobic glycolysis and required for survival of photoreceptors during aging [39, 40]. Together with the downregulation of mitochondrially encoded transcripts, our data indicate that cone metabolism is likely affected already before the completion of rod degeneration.
In conclusion, our study provides insight into transcriptomic changes in rods and cones during the primary wave of rod degeneration in the rd10 mouse, a model of human RP. Our results highlight EGR1 as a key protein during the early phase of the cell death process in rods, and of the response of cones to primary rod degeneration. We also provide evidence for mitochondrial dysregulation in rods and especially in cones as a reaction to rod loss. Mitochondrial and metabolic changes might have crucial implications for cone photoreceptor survival in retinas that degenerate due to mutations in rod-specific genes.
All animal experiments adhered to the standards of the ARVO Statement for the Use of Animals in Ophthalmic and Vision Research and the regulations of the Veterinary Office of Canton Zurich, Switzerland, under the animal experimentation license ZH141/2016 (single-cell sequencing), and ZH091/2019 (tissue collection). rd10 mice on a C57BL/6 background (IMSR Cat# JAX:004297, RRID:IMSR_JAX:004297) from Jackson Laboratories (Bar Harbor, USA) and wild type C57BL/6 controls (IMSR Cat# CRL:027, RRID:IMSR_CRL:027), from Charles River (Wilmington, MA, USA) were maintained in the Laboratory Animal Services Center facilities of the University of Zürich, with a 14 h:10 h light to dark cycle and access to food and water ad libitum. Euthanasia was performed with CO2 exposure followed by cervical dislocation for ZH141/2016 animals, and with CO2 exposure followed by decapitation for ZH091/2019 animals. All samples were collected between 10 am and 12 am to prevent a potential influence by the circadian rhythm or time since light onset .
Fundus imaging and optical coherence tomography
Pupils of the mice were dilated using drops of 1% Cyclogyl (Alcon Pharmaceuticals, Fribourg, Switzerland) and 5% Neosynephrine (Ursapharm Schweiz GmbH, Roggwil, Switzerland) 20 min before anesthesia with a subcutaneous injection of 51 μl ketamine (50 mg/kg, Parke-Davis) and 6 μl Rompun (2%, Bayer AG, Leverkusen, Germany) per 30 g body weight. Eyes were kept moist throughout the experiment with the topical application of 2% Methocel (OmniVision AG, Neuhausen, Switzerland). OCT and fundus imaging were carried out with the Micron IV system (Phoenix Research Labs, Pleasanton, CA, USA). Mice were kept on a heating pad throughout the procedure.
Preparation of single-cell suspensions from mouse retinas
A modified version of the Worthington Papain Dissociation System (Worthington, Lakewood, NJ, USA) was used for tissue dissociation. All components of the dissociation system were equilibrated with 95% O2: 5% CO2 before use. Retinas were isolated from euthanized P21 mice through a slit in the cornea and kept in equilibrated 1X EBSS (Cat. No. 24010043, ThermoFisher Scientific, Waltham, MA, USA) on ice. After removing the vitreous with surgical forceps, both retinas from each mouse were combined and transferred to an equilibrated solution containing 4.76 units of papain (Worthington) and 23.8 units of DNase (Worthington) with 0.95 mM L-cysteine and 0.47 mM EDTA in 1X EBSS for 15 (rd10) to 20 minutes (wild type) of incubation in a 37 °C water bath. The incubation time was shortened for rd10 retinas due to the increased fragility of the degenerative retina. After papain digestion, 30 units of DNase and 1.034 mg/ml ovomucoid + albumin (Worthington) were added to the retinal suspensions, which were then layered on top of a 10 mg/ml ovomucoid + albumin solution for a discontinuous density gradient and centrifuged at 70 × g for 6 min to pellet intact cells. Pelleted cells were resuspended in 1X PBS (-Ca2+, -Mg2+; ThermoFisher Scientific) + 1% BSA (ThermoFisher Scientific) for calcein-propidium iodide staining.
Calcein - propidium iodide staining
Calcein-AM (Sigma-Aldrich, St. Louis, MO, USA) at a final concentration of 2 μM and propidium iodide (Sigma-Aldrich) at a final concentration of 3 μM were added to a 1:20 diluted aliquot of each retinal suspension for cell count and viability measurements using a Neubauer cell counting chamber (Brand GmbH, Wertheim, Germany). Five million cells were used for the subsequent methanol fixation.
Methanol fixation was performed to prevent cell death and limit changes to single-cell transcriptomes during handling and transport to the sequencing facility. All steps of the methanol fixation were performed on ice. Five million cells were pelleted at 300 × g for 3 min and resuspended via dropwise addition of 1 ml 1X PBS (-Ca2+, -Mg2+; ThermoFisher Scientific) followed by gentle pipetting using a wide-bore P1000 pipette tip (ThermoFisher Scientific), centrifuged at 300 × g for 5 min at 4 °C, and resuspended again in 1 ml 1X PBS as described above. Centrifugation was repeated at 300 × g for 5 min at 4 °C. The pellet was resuspended in 100 ul 1X PBS (-Ca2+, -Mg2+), resuspended as described above with a wide-bore P200 tip. While gently vortexing the tube, 900 μL of chilled 100% methanol was added dropwise to cells. The cells were fixed on ice for 15 min, and fixation was confirmed using trypan blue staining (NanoEntek Inc, Seoul, Korea).
Methanol-fixed cells were pelleted at 3000 × g for 10 min at 4 °C and resuspended in 1X PBS with 1% BSA (ThermoFisher Scientific) and 0.5 U/μL RNasin (Promega, Madison, WI, USA) and filtered through a 40-μm Flowmi Cell Strainer (Sigma-Aldrich). From each sample generated from an individual mouse, 10,000 cells (maximum recovery) were loaded onto an individual lane of a 10X Genomics Chromium Single Cell microfluidics chip using the v3 chemistry (10X Genomics, Pleasanton, CA, USA). Libraries were pooled and sequenced on an Illumina NovaSeq 6000 platform (Illumina, San Diego, CA, USA) at the Functional Genomics Center of the University of Zurich/ETH Zurich.
Bioinformatic analysis of single-cell RNA sequencing data
Sequencing raw data was converted to FASTQ files, mapped to the GRCm38.p5 genome assembly, and used to build a cell-gene matrix in the Cell Ranger pipeline. The combined gene expression matrix generated in CellRanger, which contained data from all samples, was analyzed using the R package Seurat v3.1. Quality control was performed in multiple steps. Initially, genes expressed in less than 3 cells, and cells with less than 200 transcripts were removed from the analysis. Based on plots of the number of unique genes vs. the number of transcripts for each cell, only cells with 200–6000 unique genes and 1000–20,000 transcripts were retained. Data were normalized with log normalization and a scale factor of 10.000. Five thousand most highly variable genes with 0.1–8 average expression and > 1 dispersion were selected for downstream analyses. Prior to dimensionality reduction, the data were scaled to ensure each gene had a mean expression of 0 and variance of 1 across all cells and differences based on total UMI count were regressed out. The first 50 principal components were used for the main analyses. Cells were embedded on a k-nearest neighbor graph, and clustering was performed with modularity optimization on a shared nearest neighbor graph with a resolution parameter of 0.8. tSNE was used for visualizations. Cell cluster identity was assigned by the expression of known marker genes.
After clustering, additional filtering and cluster merging steps were applied to two specific clusters (Additional file 2: Fig. S2—cluster names containing “omitted”). We noted the presence of a small rod cell cluster composed primarily of wild-type rods. These cells had lower levels of total features per cell (Additional file 2: Fig. S2B) and transcripts per cell (Additional file 2: Fig. S2C), as well as a higher percentage of mitochondrial transcripts (Additional file 2: Fig. S2D). Differential gene expression analysis between these rods and the remaining rod clusters (Additional file 2: Fig. S2E) revealed that these cells expressed higher levels of ribosomal transcripts. Together, these features indicated the low-quality of cells in this cluster as cytoplasmic RNA leaks out of the damaged cells, while mitochondrial and ribosomal transcripts remain part of the respective organelles . A cluster co-expressing Müller glia and rod photoreceptor markers was also excluded (Additional file 2: Fig. S2F), as this likely represented a technical artifact from Müller glia in close contact with rod photoreceptors bringing transcripts from both cell types into a single droplet. Differential gene expression analysis was performed with Seurat using a Wilcoxon rank sum test and Bonferroni correction for multiple testing. Hemoglobin transcripts were removed from gene sets prior to plotting. Batch correction was not applied, as batch-like effects in the data (Fig. 1C) correlated well with known biological variability between degenerating and wild-type rods. For the rod cell subset, variable feature detection (first 3000) and PCA (first 12 dimensions) were repeated prior to trajectory analysis. For the cone cell subset, data was renormalized and rescaled prior to plotting.
Removal of doublets based on marker gene co-expression
Doublet removal based on the number of UMIs was performed as described above. For filtering doublets based on marker gene co-expression, the code for a previously published algorithm  was kindly provided by Jie Wang (Johns Hopkins University). Briefly, this approach identifies marker genes with high statistical power to define specific clusters and removes cells that express such markers for multiple cluster types. For further details, see Hoang et al. .
Pseudotime analysis for rod photoreceptors
The R package Slingshot  was used to determine the pseudotime trajectory among rods, with the wild-type C01 cluster set as the starting cluster. A transcript list for further analysis was constructed by performing differential gene expression analysis as described above on clusters adjacent in pseudotime (C01 vs C02, and C02 vs C03) and including all genes with an adjusted p-value < 0.05 (Bonferroni correction). Then, rods were separated into 20 initial bins along pseudotime, and any adjacent bins with < 100 cells were merged, resulting in 16 bins. The average expression of each significant DE transcript in each pseudotime was obtained with Seurat . The binned pseudotime matrix was then scaled for each gene, and the scaled expression values were plotted using the ComplexHeatmap package + , with k = 4 for k-means clustering. Heatmap colors were generated with the viridis palette with the R package circlize .
The expression profiles of all detected transcripts among rod photoreceptors were also fit to pseudotime values for each cell using the R package tradeSeq , which fits a generalized additive model for each gene. The associationTest function in the tradeSeq package was used to perform a Wald test to evaluate the null hypothesis that gene expression does not change along pseudotime. P-values were adjusted with Bonferroni correction.
In silico position weight matrix (PWM) analysis
Lists of differentially expressed genes with adjusted p-value < 0.05 were provided as input to Enrichr (https://maayanlab.cloud/Enrichr/#)  to analyze the presence of transcription factor binding sites in the promoter regions (− 2000 to + 500 base pairs of the transcription start site) of these genes using the PWMs found in the TRANSFAC and JASPAR databases of transcription factors. The combined score, which is corrected for the statistical bias to the differences among numbers of potential targets of transcription factors [57, 112], was used to assess the results.
Correlation analysis for rd10 and wild-type cone photoreceptors
Rod clusters C00, C01, C02 and the cone cluster C04 were isolated from the main dataset, and the cone cluster C04 was separated into rd10 and wild-type cones. The resulting five groups were used for the correlation analysis as follows: the average expression of 5000 most variable genes determined during the analysis of the main dataset was calculated for each group, and Spearman’s correlation coefficient was calculated using the average expression values with the R function cor.
Cell culture and RNA interference
661W cells (RRID:CVCL_6240) were a kind gift from Dr. Muayyad al-Ubaidi . Cells were cultured in high-glucose Dulbecco’s modified Eagle’s medium (DMEM) (Sigma-Aldrich) with 10% fetal bovine serum (FBS) and 100 U/ml penicillin and 100 μg/ml streptomycin (ThermoFisher) and maintained in a 37 °C humid incubator with 5% CO2. For siRNA transfection, 20,000 cells were seeded into each well of a 24-well plate (TPP, Trasadingen, Switzerland) 1 day prior. Transfections of siEgr1 (SMARTPool siGENOME Mouse Egr1 siRNA M-040286-01, Horizon Discovery Ltd., Waterbeach, UK) and a control siRNA (AllStars Negative Control siRNA, Qiagen, Hilden, Germany) were performed with Lipofectamine RNAiMAX reagent, according to kit instructions (ThermoFisher). Cells were collected 24 h after transfection for RNA extraction according to RNA extraction kit instructions (Macherey-Nagel NucleoSpin RNA kit, Oensingen, Switzerland).
RNA extraction, cDNA synthesis, and semi-quantitative real-time PCR
Retinas were isolated from euthanized mice through a slit in the cornea, and snap-frozen in liquid nitrogen. Total RNA was isolated with a commercial, column-based RNA isolation kit with on-column DNaseI treatment (Macherey-Nagel NucleoSpin RNA kit). cDNA synthesis was carried out with oligo-(d)T selection and M-MLV reverse transcriptase (Promega). For semi-quantitative real-time PCR, 10 ng cDNA was used alongside primers targeting the gene of interest and the PowerUp SYBR Green Master Mix (ThermoFisher Scientific) in the ABI QuantStudio 3 system (ThermoFisher Scientific). Primers (Additional file 15) were designed, whenever possible, to span intron sequences or exon-exon boundaries to avoid genomic DNA amplification and tested for specificity and efficiency before experiments. mt-Atp6 and mt-Co2 primers were previously used in . mt-Nd4 primers were previously used in . Beta-actin (Actb) was used as the reference housekeeping gene. Data analysis was carried out with the R statistical programming language, using the 2^-ddCt method . For real-time PCR on retinal tissue, 3 biological replicates (individual mice from which retinas were isolated) were used per time point per strain. Two-way ANOVA (for effect of genotype and age) with Tukey’s range test was carried out on the resulting relative expression values with the R package rstatix. For real-time PCR on 661W cells, three replicate wells were transfected with the same siRNA (or left untreated) in each of the three individual experiments. Each experiment’s mean was used to calculate the mean of means and standard error of the mean. One-way ANOVA with Tukey’s range test was carried out on the three experimental means to determine the statistical significance of the results.
Immunofluorescence on retinal cryosections
Eyes were marked dorsally via cauterization, and 12 μm sections were prepared as previously described . Sections were blocked either in blocking solution (PB-Salt (0.1 M phosphate buffer with 0.8% NaCl and 0.02% KCl) with 0.3% Triton-X (Sigma-Aldrich) and 3% goat serum (Sigma-Aldrich) or 2% horse serum (Sigma-Aldrich)), followed by primary antibody incubation overnight at 4 °C. After three washes with PB-Salt, secondary antibody incubation was performed for an hour at room temperature. Nuclear counterstaining was performed with 4′,6-diamidine-2′-phenylindole di-hydrochloride (DAPI, Thermo Fisher Scientific). Sections were mounted using Mowiol (Sigma-Aldrich) and imaged with a Leica DMI6000B inverted fluorescent microscope (Leica Microsystems, Wetzlar, Germany). Primary and secondary antibodies used in the paper are listed in Additional file 16. Z-stacks were merged with max projection in ImageJ. Images were optimized for histogram minimum and maximum values using Adobe Photoshop.
Availability of data and materials
The scRNA-seq dataset supporting the conclusions of this article is available in the Gene Expression Omnibus repository (GSE183206).
Hartong DT, Berson EL, Dryja TP. Retinitis pigmentosa. Lancet. 2006;368:1795–809.
Bunce C, Xing W, Wormald R. Causes of blind and partial sight certifications in England and Wales: April 2007-March 2008. Eye. 2010;24:1692–9.
RetNet: Summaries. https://sph.uth.edu/retnet/sum-dis.htm. Accessed 4 Feb 2021.
Chang B, Hawes NL, Pardue MT, German AM, Hurd RE, Davisson MT, et al. Two mouse retinal degenerations caused by missense mutations in the beta-subunit of rod cGMP phosphodiesterase gene. Vision Res. 2007;47:624–33.
Wang T, Reingruber J, Woodruff ML, Majumder A, Camarena A, Artemyev NO, et al. The PDE6 mutation in the rd10 retinal degeneration mouse model causes protein mislocalization and instability and promotes cell death through increased ion influx. J Biol Chem. 2018;293:15332–46.
Gargini C, Terzibasi E, Mazzoni F, Strettoi E. Retinal organization in the retinal degeneration 10 (rd10) mutant mouse: a morphological and ERG study. J Comp Neurol. 2007;500:222–38.
Xiong W, MacColl Garfinkel AE, Li Y, Benowitz LI, Cepko CL. NRF2 promotes neuronal survival in neurodegeneration and acute nerve damage. J Clin Invest. 2015;125:1433–45.
Samardzija M, Wariwoda H, Imsand C, Huber P, Heynen SR, Gubler A, et al. Activation of survival pathways in the degenerating retina of rd10 mice. Exp Eye Res. 2012;99:17–26.
Léveillard T, Mohand-Saïd S, Lorentz O, Hicks D, Fintz A-C, Clérin E, et al. Identification and characterization of rod-derived cone viability factor. Nat Genet. 2004;36:755–9.
Aït-Ali N, Fridlich R, Millet-Puel G, Clérin E, Delalande F, Jaillard C, et al. Rod-derived cone viability factor promotes cone survival by stimulating aerobic glycolysis. Cell. 2015;161:817–32.
Punzo C, Kornacker K, Cepko CL. Stimulation of the insulin/mTOR pathway delays cone death in a mouse model of retinitis pigmentosa. Nat Neurosci. 2009;12:44–52.
Shen J, Yang X, Dong A, Petters RM, Peng Y-W, Wong F, et al. Oxidative damage is a potential cause of cone cell death in retinitis pigmentosa. J Cell Physiol. 2005;203:457–64.
Jeon CJ, Strettoi E, Masland RH. The major cell populations of the mouse retina. J Neurosci. 1998;18:8936–46.
Macosko EZ, Basu A, Satija R, Nemesh J, Shekhar K, Goldman M, et al. Highly parallel genome-wide expression profiling of individual cells using nanoliter droplets. Cell. 2015;161:1202–14.
Rohrer B, Pinto FR, Hulse KE, Lohr HR, Zhang L, Almeida JS. Multidestructive pathways triggered in photoreceptor cell death of the rd mouse as determined through gene expression profiling. J Biol Chem. 2004;279:41903–10.
Sharma YV, Cojocaru RI, Ritter LM, Khattree N, Brooks M, Scott A, et al. Protective gene expression changes elicited by an inherited defect in photoreceptor structure. PLoS One. 2012;7:e31371.
Uren PJ, Lee JT, Doroudchi MM, Smith AD, Horsager A. A profile of transcriptomic changes in the rd10 mouse model of retinitis pigmentosa. Mol Vis. 2014;20:1612–28.
Guadagni V, Biagioni M, Novelli E, Aretini P, Mazzanti CM, Strettoi E. Rescuing cones and daylight vision in retinitis pigmentosa mice. FASEB J. 2019;33:10177–92.
Lukowski SW, Lo CY, Sharov AA, Nguyen Q, Fang L, Hung SS, et al. A single-cell transcriptome atlas of the adult human retina. EMBO J. 2019;38(18):e100811.
Menon M, Mohammadi S, Davila-Velderrain J, Goods BA, Cadwell TD, Xing Y, et al. Single-cell transcriptomic atlas of the human retina identifies cell types associated with age-related macular degeneration. Nat Commun. 2019;10:4902.
Cowan CS, Renner M, De Gennaro M, Gross-Scherf B, Goldblum D, Hou Y, et al. Cell types of the human retina and its organoids at single-cell resolution. Cell. 2020;182:1623–1640.e34.
Todd L, Palazzo I, Suarez L, Liu X, Volkov L, Hoang TV, et al. Reactive microglia and IL1β/IL-1R1-signaling mediate neuroprotection in excitotoxin-damaged mouse retina. J Neuroinflammation. 2019;16:118.
Hoang T, Wang J, Boyd P, Wang F, Santiago C, Jiang L, et al. Gene regulatory networks controlling vertebrate retinal regeneration. Science. 2020. https://doi.org/10.1126/science.abb8598.
Clark BS, Stein-O’Brien GL, Shiau F, Cannon GH, Davis-Marcisak E, Sherman T, et al. Single-cell RNA-Seq analysis of retinal development identifies NFI factors as regulating mitotic exit and late-born cell specification. Neuron. 2019;102:1111–1126.e5.
Zheng GXY, Terry JM, Belgrader P, Ryvkin P, Bent ZW, Wilson R, et al. Massively parallel digital transcriptional profiling of single cells. Nat Commun. 2017;8:14049.
Butler A, Hoffman P, Smibert P, Papalexi E, Satija R. Integrating single-cell transcriptomic data across different conditions, technologies, and species. Nat Biotechnol. 2018;36:411–20.
Hoang T, Wang J, Boyd P, Wang F, Santiago C, Jiang L, et al. Cross-species transcriptomic and epigenomic analysis reveals key regulators of injury response and neuronal regeneration in vertebrate retinas. bioRxiv. 2019:717876.
Shekhar K, Lapan SW, Whitney IE, Tran NM, Macosko EZ, Kowalczyk M, et al. Comprehensive classification of retinal bipolar neurons by single-cell transcriptomics. Cell. 2016;166:1308–1323.e30.
Howe KL, Achuthan P, Allen J, Allen J, Alvarez-Jarreta J, Amode MR, et al. Ensembl 2021. Nucleic Acids Res. 2021;49:D884–91.
Power MJ, Rogerson LE, Schubert T, Berens P, Euler T, Paquet-Durand F. Systematic spatiotemporal mapping reveals divergent cell death pathways in three mouse models of hereditary retinal degeneration. J Comp Neurol. 2020;528:1113–39.
Street K, Risso D, Fletcher RB, Das D, Ngai J, Yosef N, et al. Slingshot: cell lineage and pseudotime inference for single-cell transcriptomics. BMC Genomics. 2018;19:477.
Gu Z, Eils R, Schlesner M. Complex heatmaps reveal patterns and correlations in multidimensional genomic data. Bioinformatics. 2016;32:2847–9.
Krizaj D, Demarco SJ, Johnson J, Strehler EE, Copenhagen DR. Cell-specific expression of plasma membrane calcium ATPase isoforms in retinal neurons. J Comp Neurol. 2002;451:1–21.
Liu B, Suyeoka G, Papa S, Franzoso G, Neufeld AH. Growth arrest and DNA damage protein 45b (Gadd45b) protects retinal ganglion cells from injuries. Neurobiol Dis. 2009;33:104–10.
Zhang H, Chen X, Wang J, Guang W, Han W, Zhang H, et al. EGR1 decreases the malignancy of human non-small cell lung carcinoma by regulating KRT18 expression. Sci Rep. 2014;4:5416.
Wang W, Wang Q, Wan D, Sun Y, Wang L, Chen H, et al. Histone HIST1H1C/H1.2 regulates autophagy in the development of diabetic retinopathy. Autophagy. 2017;13:941–54.
Koriyama Y, Sugitani K, Ogai K, Kato S. Heat shock protein 70 induction by valproic acid delays photoreceptor cell death by N-methyl-N-nitrosourea in mice. J Neurochem. 2014;130:707–19.
Samardzija M, Todorova V, Gougoulakis L, Barben M, Nötzli S, Klee K, et al. Light stress affects cones and horizontal cells via rhodopsin-mediated mechanisms. Exp Eye Res. 2019;186:107719.
Weh E, Lutrzykowska Z, Smith A, Hager H, Pawar M, Wubben TJ, et al. Hexokinase 2 is dispensable for photoreceptor development but is required for survival during aging and outer retinal stress. Cell Death Dis. 2020;11:422.
Zhang R, Shen W, Du J, Gillies MC. Selective knockdown of hexokinase 2 in rods leads to age-related photoreceptor degeneration and retinal metabolic remodeling. Cell Death Dis. 2020;11:885.
Chen L, Wu W, Dentchev T, Zeng Y, Wang J, Tsui I, et al. Light damage induced changes in mouse retinal gene expression. Exp Eye Res. 2004;79:239–47.
Joly S, Lange C, Thiersch M, Samardzija M, Grimm C. Leukemia inhibitory factor extends the lifespan of injured photoreceptors in vivo. J Neurosci. 2008;28:13765–74.
Bürgi S, Samardzija M, Grimm C. Endogenous leukemia inhibitory factor protects photoreceptor cells against light-induced degeneration. Mol Vis. 2009;15:1631–7.
Jomary C, Chatelain G, Michel D, Weston A, Neal MJ, Jones SE. Effect of targeted expression of clusterin in photoreceptor cells on retinal development and differentiation. J Cell Sci. 1999;112(Pt 10):1455–64.
Tuson M, Marfany G, Gonzàlez-Duarte R. Mutation of CERKL, a novel human ceramide kinase gene, causes autosomal recessive retinitis pigmentosa (RP26). Am J Hum Genet. 2004;74:128–38.
Jin H, White SR, Shida T, Schulz S, Aguiar M, Gygi SP, et al. The conserved Bardet-Biedl syndrome proteins assemble a coat that traffics membrane proteins to cilia. Cell. 2010;141:1208–19.
Zhang T, Wei Y, Jiang X, Li J, Qiu S, Zhang S. Protection of photoreceptors by intravitreal injection of the Y-27632 Rho-associated protein kinase inhibitor in Royal College of Surgeons rats. Mol Med Rep. 2015;12:3655–61.
Townes-Anderson E, Wang J, Halász É, Sugino I, Pitler A, Whitehead I, et al. Fasudil, a clinically used ROCK inhibitor, stabilizes rod photoreceptor synapses after retinal detachment. Transl Vis Sci Technol. 2017;6:22.
di Rago JP, Macadre C, Lazowska J, Slonimski PP. The C-terminal domain of yeast cytochrome b is essential for a correct assembly of the mitochondrial cytochrome bc1 complex. FEBS Lett. 1993;328:153–8.
Di Gioia SA, Letteboer SJF, Kostic C, Bandah-Rozenfeld D, Hetterschijt L, Sharon D, et al. FAM161A, associated with retinitis pigmentosa, is a component of the cilia-basal body complex and interacts with proteins involved in ciliopathies. Hum Mol Genet. 2012;21:5174–84.
Khan S, Ullah I, Irfanullah, Touseef M, Basit S, Khan MN, et al. Novel homozygous mutations in the genes ARL6 and BBS10 underlying Bardet-Biedl syndrome. Gene. 2013;515:84–8.
Van Schil K, Karlstetter M, Aslanidis A, Dannhausen K, Azam M, Qamar R, et al. Autosomal recessive retinitis pigmentosa with homozygous rhodopsin mutation E150K and non-coding cis-regulatory variants in CRX-binding regions of SAMD7. Sci Rep. 2016;6:21307.
Hlawatsch J, Karlstetter M, Aslanidis A, Lückoff A, Walczak Y, Plank M, et al. Sterile alpha motif containing 7 (samd7) is a novel crx-regulated transcriptional repressor in the retina. PLoS One. 2013;8:e60633.
Omori Y, Kubo S, Kon T, Furuhashi M, Narita H, Kominami T, et al. Samd7 is a cell type-specific PRC1 component essential for establishing retinal rod photoreceptor identity. Proc Natl Acad Sci U S A. 2017;114:E8264–73.
Corton M, Avila-Fernández A, Campello L, Sánchez M, Benavides B, López-Molina MI, et al. Identification of the photoreceptor transcriptional co-repressor SAMD11 as novel cause of autosomal recessive retinitis pigmentosa. Sci Rep. 2016;6:35370.
Van den Berge K, Roux de Bézieux H, Street K, Saelens W, Cannoodt R, Saeys Y, et al. Trajectory-based differential expression analysis for single-cell sequencing data. Nat Commun. 2020;11:1201.
Chen EY, Tan CM, Kou Y, Duan Q, Wang Z, Meirelles GV, et al. Enrichr: interactive and collaborative HTML5 gene list enrichment analysis tool. BMC Bioinformatics. 2013;14:128.
Wheway G, Nazlamova L, Turner D, Cross S. 661W photoreceptor cell line as a cell model for studying retinal ciliopathies. Front Genet. 2019;10:308.
Yu J, de Belle I, Liang H, Adamson ED. Coactivating factors p300 and CBP are transcriptionally crossregulated by Egr1 in prostate cells, leading to divergent responses. Mol Cell. 2004;15:83–94.
Kim J, Kang S-M, Oh SY, Kang SH, Lee I, Hwang JC, et al. Early growth response 1-dependent downregulation of matrix metalloproteinase 9 and mouse double minute 2 attenuates head and neck squamous cell carcinoma metastasis. Cell Physiol Biochem. 2018;50:1869–81.
Wenzel A, Grimm C, Seeliger MW, Jaissle G, Hafezi F, Kretschmer R, et al. Prevention of photoreceptor apoptosis by activation of the glucocorticoid receptor. Invest Ophthalmol Vis Sci. 2001;42:1653–9.
Lee SK, Anzick SL, Choi JE, Bubendorf L, Guan XY, Jung YK, et al. A nuclear factor, ASC-2, as a cancer-amplified transcriptional coactivator essential for ligand-dependent transactivation by nuclear receptors in vivo. J Biol Chem\. 1999;274:34283–93.
Olkku A, Mahonen A. Calreticulin mediated glucocorticoid receptor export is involved in beta-catenin translocation and Wnt signalling inhibition in human osteoblastic cells. Bone. 2009;44:555–65.
Burns K, Duggan B, Atkinson EA, Famulski KS, Nemer M, Bleackley RC, et al. Modulation of gene expression by calreticulin binding to the glucocorticoid receptor. Nature. 1994;367:476–80.
Stubbs FE, Birnie MT, Biddie SC, Lightman SL, Conway-Campbell BL. SKOV3 cells containing a truncated ARID1a protein have a restricted genome-wide response to glucocorticoids. Mol Cell Endocrinol. 2018;461:226–35.
Kolandaivelu S, Chang B, Ramamurthy V. Rod phosphodiesterase-6 (PDE6) catalytic subunits restore cone function in a mouse model lacking cone PDE6 catalytic subunit. J Biol Chem. 2011;286:33252–9.
Siegert S, Cabuy E, Scherf BG, Kohler H, Panda S, Le Y-Z, et al. Transcriptional code and disease map for adult retinal cell types. Nat Neurosci. 2012;15(487–95):S1.
Brennenstuhl C, Tanimoto N, Burkard M, Wagner R, Bolz S, Trifunovic D, et al. Targeted ablation of the Pde6h gene in mice reveals cross-species differences in cone and rod phototransduction protein isoform inventory. J Biol Chem. 2015;290:10242–55.
Roche SL, Wyse-Jackson AC, Byrne AM, Ruiz-Lopez AM, Cotter TG. Alterations to retinal architecture prior to photoreceptor loss in a mouse model of retinitis pigmentosa. Int J Dev Biol. 2016;60:127–39.
Simon P, Schott K, Williams RW, Schaeffel F. Posttranscriptional regulation of the immediate-early gene EGR1 by light in the mouse retina. Eur J Neurosci. 2004;20:3371–7.
Bahrami S, Drabløs F. Gene regulation in the immediate-early response process. Adv Biol Regul. 2016;62:37–49.
Duclot F, Kabbaj M. The role of early growth response 1 (EGR1) in brain plasticity and neuropsychiatric disorders. Front Behav Neurosci. 2017;11:35.
Agca C, Gubler A, Traber G, Beck C, Imsand C, Ail D, et al. p38 MAPK signaling acts upstream of LIF-dependent neuroprotection during photoreceptor degeneration. Cell Death Dis. 2013;4:e785.
Lim CP, Jain N, Cao X. Stress-induced immediate-early gene, egr-1, involves activation of p38/JNK1. Oncogene. 1998;16:2915–26.
Fontainhas AM, Townes-Anderson E. RhoA inactivation prevents photoreceptor axon retraction in an in vitro model of acute retinal detachment. Invest Ophthalmol Vis Sci. 2011;52:579–87.
Mullin M, Lightfoot K, Clarke R, Miller M, Lahesmaa R, Cantrell D. The RhoA transcriptional program in pre-T cells. FEBS Lett. 2007;581:4309–17.
Gehrig A, Langmann T, Horling F, Janssen A, Bonin M, Walter M, et al. Genome-wide expression profiling of the retinoschisin-deficient retina in early postnatal mouse development. Invest Ophthalmol Vis Sci. 2007;48:891–900.
Yasuda M, Tanaka Y, Ryu M, Tsuda S, Nakazawa T. RNA sequence reveals mouse retinal transcriptome changes early after axonal injury. PLoS One. 2014;9:e93258.
Ebner LJA, Samardzija M, Storti F, Todorova V, Karademir D, Behr J, et al. Transcriptomic analysis of the mouse retina after acute and chronic normobaric and hypobaric hypoxia. Sci Rep. 2021;11:16666.
Freeman WM, Bixler GV, Brucklacher RM, Walsh E, Kimball SR, Jefferson LS, et al. Transcriptomic comparison of the retina in two mouse models of diabetes. J Ocul Biol Dis Infor. 2009;2:202–13.
Golob-Schwarzl N, Bettermann K, Mehta AK, Kessler SM, Unterluggauer J, Krassnig S, et al. High Keratin 8/18 ratio predicts aggressive hepatocellular cancer phenotype. Transl Oncol. 2019;12:256–68.
Kim J-W, Yang H-J, Brooks MJ, Zelinger L, Karakülah G, Gotoh N, et al. NRL-regulated transcriptome dynamics of developing rod photoreceptors. Cell Rep. 2016;17:2460–73.
Haider NB, Mollema N, Gaule M, Yuan Y, Sachs AJ, Nystuen AM, et al. Nr2e3-directed transcriptional regulation of genes involved in photoreceptor development and cell-type specific phototransduction. Exp Eye Res. 2009;89:365–72.
Langmann T, Ebert S, Walczak Y, Weigelt K, Ehrengruber MU, Stiewe T, et al. Induction of early growth response-1 mediates microglia activation in vitro but is dispensable in vivo. Neuromolecular Med. 2009;11:87–96.
Bernaudin M, Tang Y, Reilly M, Petit E, Sharp FR. Brain genomic response following hypoxia and re-oxygenation in the neonatal Rat: identification of genes that might contribute to hypoxia-induced ischemic tolerance *. J Biol Chem. 2002;277:39728–38.
Hirata T, Cui YJ, Funakoshi T, Mizukami Y, Ishikawa Y-I, Shibasaki F, et al. The temporal profile of genomic responses and protein synthesis in ischemic tolerance of the rat brain induced by repeated hyperbaric oxygen. Brain Res. 2007;1130:214–22.
Benoit ME, Tenner AJ. Complement protein C1q-mediated neuroprotection is correlated with regulation of neuronal gene and microRNA expression. J Neurosci. 2011;31:3459–69.
Pulido-Salgado M, Vidal-Taboada JM, Saura J. C/EBPβ and C/EBPδ transcription factors: basic biology and roles in the CNS. Prog Neurobiol. 2015;132:1–33.
Lopez de Heredia L, Magoulas C. Lack of the transcription factor C/EBPδ impairs the intrinsic capacity of peripheral neurons for regeneration. Exp Neurol. 2013;239:148–57.
Jones BW, Pfeiffer RL, Ferrell WD, Watt CB, Marmor M, Marc RE. Retinal remodeling in human retinitis pigmentosa. Exp Eye Res. 2016;150:149–65.
Aleman TS, Cideciyan AV, Sumaroka A, Schwartz SB, Roman AJ, Windsor EAM, et al. Inner retinal abnormalities in X-linked retinitis pigmentosa with RPGR mutations. Invest Ophthalmol Vis Sci. 2007;48:4759–65.
Marc RE, Jones BW, Anderson JR, Kinard K, Marshak DW, Wilson JH, et al. Neural reprogramming in retinal degeneration. Invest Ophthalmol Vis Sci. 2007;48:3364–71.
Strettoi E, Pignatelli V, Rossi C, Porciatti V, Falsini B. Remodeling of second-order neurons in the retina of rd/rd mutant mice. Vision Res. 2003;43:867–77.
Phillips MJ, Otteson DC, Sherry DM. Progression of neuronal and synaptic remodeling in the rd10 mouse model of retinitis pigmentosa. J Comp Neurol. 2010;518:2071–89.
Kröger S, Horton SE, Honig LS. The developing avian retina expresses agrin isoforms during synaptogenesis. J Neurobiol. 1996;29:165–82.
Fuerst PG, Rauch SM, Burgess RW. Defects in eye development in transgenic mice overexpressing the heparan sulfate proteoglycan agrin. Dev Biol. 2007;303:165–80.
Gautam M, Noakes PG, Moscoso L, Rupp F, Scheller RH, Merlie JP, et al. Defective neuromuscular synaptogenesis in agrin-deficient mutant mice. Cell. 1996;85:525–35.
Zhang H, Sathyamurthy A, Liu F, Li L, Zhang L, Dong Z, et al. Agrin-Lrp4-Ror2 signaling regulates adult hippocampal neurogenesis in mice. Elife. 2019;8:e45303.
Handara G, Hetsch FJA, Jüttner R, Schick A, Haupt C, Rathjen FG, et al. The role of agrin, Lrp4 and MuSK during dendritic arborization and synaptogenesis in cultured embryonic CNS neurons. Dev Biol. 2019;445:54–67.
MacDonald R, Barbat-Artigas S, Cho C, Peng H, Shang J, Moustaine A, et al. A novel Egr-1-agrin pathway and potential implications for regulation of synaptic physiology and homeostasis at the neuromuscular junction. Front Aging Neurosci. 2017;9:258.
Araki A, Maruyama R, Harada Y, Ishikawa N, Harada T. Analysis of the light-sensitivity of the photoreceptor cells of the ataxia and male sterility (AMS) mouse, an Nna1 mutant. Pathol Int. 2012;62:719–27.
Bosch Grau M, Masson C, Gadadhar S, Rocha C, Tort O, Marques Sousa P, et al. Alterations in the balance of tubulin glycylation and glutamylation in photoreceptors leads to retinal degeneration. J Cell Sci. 2017;130:938–49.
Xu L, Kong L, Wang J, Ash JD. Stimulation of AMPK prevents degeneration of photoreceptors and the retinal pigment epithelium. Proc Natl Acad Sci U S A. 2018;115:10475–80.
Zou T, He J, Chen X, Sun D, Fan X, et al. Rescue of retinal degeneration in rd1 mice by intravitreally injected metformin. Front Mol Neurosci. 2019;12:102.
He L, Wondisford FE. Metformin action: concentrations matter. Cell Metab. 2015;21:159–62.
Zhou G, Myers R, Li Y, Chen Y, Shen X, Fenyk-Melody J, et al. Role of AMP-activated protein kinase in mechanism of metformin action. J Clin Invest. 2001;108:1167–74.
Onichtchouk D, Chen YG, Dosch R, Gawantka V, Delius H, Massagué J, et al. Silencing of TGF-beta signalling by the pseudoreceptor BAMBI. Nature. 1999;401:480–5.
Namwanje M, Brown CW. Activins and inhibins: roles in development, physiology, and disease. Cold Spring Harb Perspect Biol. 2016;8(7):a021881.
Yoshida S, Mears AJ, Friedman JS, Carter T, He S, Oh E, et al. Expression profiling of the developing and mature Nrl-/- mouse retina: identification of retinal disease candidates and transcriptional regulatory targets of Nrl. Hum Mol Genet. 2004;13:1487–503.
Ilicic T, Kim JK, Kolodziejczyk AA, Bagger FO, McCarthy DJ, Marioni JC, et al. Classification of low quality cells from single-cell RNA-seq data. Genome Biol. 2016;17:29.
Gu Z, Gu L, Eils R, Schlesner M, Brors B. circlize implements and enhances circular visualization in R. Bioinformatics. 2014;30:2811–2.
Kuleshov MV, Jones MR, Rouillard AD, Fernandez NF, Duan Q, Wang Z, et al. Enrichr: a comprehensive gene set enrichment analysis web server 2016 update. Nucleic Acids Res. 2016;44:W90–7.
al-Ubaidi MR, Font RL, Quiambao AB, Keener MJ, Liou GI, Overbeek PA, et al. Bilateral retinal and brain tumors in transgenic mice expressing simian virus 40 large T antigen under control of the human interphotoreceptor retinoid-binding protein promoter. J Cell Biol. 1992;119:1681–7.
Wall CE, Whyte J, Suh JM, Fan W, Collins B, Liddle C, et al. High-fat diet and FGF21 cooperatively promote aerobic thermogenesis in mtDNA mutator mice. Proc Natl Acad Sci U S A. 2015;112:8714–9.
El-Merhie N, Baumgart-Vogt E, Pilatz A, Pfreimer S, Pfeiffer B, Pak O, et al. Differential alterations of the mitochondrial morphology and respiratory chain complexes during postnatal development of the mouse lung. Oxid Med Cell Longev. 2017;2017:9169146.
Livak KJ, Schmittgen TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2(-Delta Delta C(T)) method. Methods. 2001;25:402–8.
Barben M, Schori C, Samardzija M, Grimm C. Targeting Hif1a rescues cone degeneration and prevents subretinal neovascularization in a model of chronic hypoxia. Mol Neurodegener. 2018;13:12.
The authors would like to thank Sarah Nötzli and Cornelia Imsand (University of Zürich) for excellent technical support; Dr. Katrin Klee (University of Zürich) for assistance with OCT imaging; Dr. Ge Tan and Dr. Emilio Yángüez (Functional Genomics Center Zürich) for their invaluable guidance in the single-cell transcriptomics experiments and bioinformatics analyses; Dr. Clayton Santiago (Johns Hopkins University), Dr. Paulo Czarnewski, and Dr. Erik Fasterius (National Bioinformatics Infrastructure Sweden) for helpful discussions.
This research was performed with funding by the Swiss National Science Foundation (SNF), grant number 31003A_173008.
Ethics approval and consent to participate
Animal maintenance and experimentation was performed in accordance with the standards of the ARVO Statement for the Use of Animals in Ophthalmic and Vision Research and the regulations of the Veterinary Office of Canton Zürich, Switzerland, under the licenses ZH141/2016 and ZH091/2019. No human participants, tissue, or data were used in this study.
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Validation of the degenerative state of rd10 retinas at P21. Fundus images (left panels) and OCT scans (right panels) of (A) rd10 mouse 1, (B) rd10 mouse 2, and (C) a wild-type mouse at P21. Mice were littermates of rd10 and wild-type mice used for scRNAseq. Red lines in fundus images mark the location of OCT scans. The vertical black bars on OCT images mark the outer retina (outer plexiform layer to photoreceptor segments). Scale bars: 100 μm.
Additional quality control measures on the single-cell transcriptomic dataset. (A) Two-dimensional tSNE plot of the whole dataset before low-quality clusters were filtered. Filtered clusters have the tag “: Omitted” and are labeled with “M” (Müller glia-rods) and “R” (Rods) in the tSNE plot. (B-D) Quality control measures for the omitted rod cluster, compared to the three retained rod clusters C01, C02 and C03. The omitted cluster has lower numbers of total features (B) and transcripts (C) but a higher percentage of mitochondrial transcripts (“mt-genes”) per cell (D). (E) Differential gene expression analysis between omitted and retained rod clusters shows that omitted rods have higher levels of ribosomal transcripts. (F) Differential gene expression analysis between omitted and retained Müller glia cluster shows that omitted Müller glia are positive for rod cell markers (e.g., Sag, Rho) and express lower levels of Müller glia markers (e.g., Dkk3).
Marker gene analysis for the single-cell transcriptomic dataset. To determine which transcripts distinguish each cluster in the dataset for cluster annotation, differential gene expression analysis between each cluster and the remainder of the dataset was performed with the Seurat package.
Early differential gene expression analysis in rods. Wild type rod cluster C01 and rd10 rod cluster C02, which were adjacent in the pseudotime analysis, were compared with differential gene expression analysis to determine transcriptomic differences that occurred with increasing pseudotime values. Avg_logFC values larger than 0 indicate that the gene is upregulated in C02 compared to C01.
Late differential gene expression analysis in rods. rd10 rod cluster C02 and rd10 rod cluster C03, which were adjacent in the pseudotime analysis, were compared with differential gene expression analysis to determine transcriptomic differences that occurred with increasing pseudotime values. Avg_logFC values larger than 0 indicate that the gene is upregulated in C03 compared to C02.
k-means clustering of significantly DE genes among rod clusters according to expression changes along increasing (binned) pseudotime. Rods were grouped into 16 bins along increasing pseudotime. Expression values for all significantly DE transcripts were averaged in each bin, scaled, and the values were used for k-means clustering (k=4) to cluster transcripts into groups of similar behaviour in pseudotime (and thus progressing degeneration). DE: differentially expressed.
Association of gene expression to pseudotime values in the rod dataset. To determine which transcripts were associated with pseudotime, we performed a Wald test using the tradeSeq package (associationTest, no log2 fold-change threshold). The expression pattern of each gene in the entire dataset is modeled along pseudotime, and the Wald statistic is calculated to indicate whether the gene is associated with pseudotime. DE: differentially expressed.
In silico position weight matrix analysis on transcripts differentially expressed in early degeneration. To determine which transcription factors might be associated with significantly differentially expressed transcripts identified in the early differential gene expression analysis (Additional file 2), we used Enrichr to perform an in silico search for sequence motifs representing transcription factor binding sites in the promoter region (-2000 to + 500 bp of the transcriptional start site) of these genes. The Overlap column refers to the number of genes in the input list (differentially regulated genes; numerator) compared to all genes in the database (denominator) that contain the transcription factor binding motif. P-values were calculated with Fisher’s exact test, and adjusted with the Benjamini-Hochberg method. The Combined Score is log(p-value)*z-score, where z-score is the deviation from the expected rank calculated from many random gene lists [57, 112].
In silico position weight matrix analysis on transcripts differentially expressed in late degeneration. To determine which transcription factors might be associated with significantly differentially expressed transcripts identified in the late differential gene expression analysis (Additional file 3), we used Enrichr to perform an in silico search for sequence motifs representing transcription factor binding sites in the promoter region (-2000 to + 500 bp of the transcriptional start site) of these genes. The Overlap column refers to the number of genes in the input list (differentially regulated genes; numerator) compared to all genes in the database (denominator) that contain the transcription factor binding motif. P-values were calculated with Fisher’s exact test, and adjusted with the Benjamini-Hochberg method. The Combined Score is log(p-value)*z-score, where z-score is the deviation from the expected rank calculated from many random gene lists [57, 112].
Expression levels of potential target genes after downregulation of Egr1 in vitro. Expression levels of Egr1, potential EGR1-target genes (Agtpbp1, Gadd45b) and a non-target gene (Apoe) were determined in 661W cells 24 hours after transfection of an siRNA against Egr1 (“siEgr1”, blue) or control siRNA (“siCtrl”, green), or in untransfected controls (“No TF”, red). Shown are the mean of three independent experiments (diamonds) ± standard error of the mean (error bars). Each individual point refers to the mean of an independent experiment. *: p < 0.05; one-way ANOVA with Tukey’s range test. TF: transfection.
Differential gene expression analysis among rd10 and wild-type cones in cluster C04. Differential gene expression analysis was performed to determine transcriptomic changes in rd10 cones compared to wild types at P21 . All transcripts were additionally checked for significant differential expression in early rod degeneration (Additional file 2) or late rod degeneration (Additional file 3). Avg_logFC values larger than 0 indicate that the gene is upregulated in rd10 cones compared to wild types. DE: differentially expressed.
Relative expression levels of selected genes in rd10 and wild-type retinas. Expression levels of selected genes were determined in retinas of rd10 (blue) and wild-type (red) mice at various ages, relative to the average of C57BL/6 P21 retinas. Shown are boxplots with individual data points. N=3 per time point per strain. *: p < 0.05, **: p < 0.01, ***: p < 0.001, ****: p < 0.0001; two-way ANOVA with Tukey’s range test.
Additional stainings for EGR1 in rd10 and wild-type retinas. Co-immunofluorescence for EGR1 (red) and GS (green) as a Müller glia marker was performed in retinas of (A) wild-type control at P21 or (B-C) rd10 mice at P21. (D) Co-immunofluorescence for (red) and OPN1SW (green) as a cone marker in a retinal section of a wild-type control mouse at P35. Arrows point to EGR1-positive Müller glia nuclei close to the outer rim of the INL. The arrowhead indicates an EGR1-positive Müller glia nucleus at its expected position. Scale bars: 50 μm. PS: Photoreceptor segments. ONL: outer nuclear layer. INL: inner nuclear layer. IPL: inner plexiform layer. GCL: ganglion cell layer.
Egr1 expression in the whole dataset. Two-dimensional UMAP plot of the whole dataset, with cells colored according to their normalized expression levels of Egr1. For cluster identities, see B. (B) Heatmap of average normalized expression levels of Egr1 within each cluster, split by genotype. Cluster 28 is composed only of rd10 cells. wt: wild type.
Antibodies used for immunofluorescence.
About this article
Cite this article
Karademir, D., Todorova, V., Ebner, L.J.A. et al. Single-cell RNA sequencing of the retina in a model of retinitis pigmentosa reveals early responses to degeneration in rods and cones. BMC Biol 20, 86 (2022). https://doi.org/10.1186/s12915-022-01280-9
- Single-cell RNA sequencing
- Retinal degeneration
- Retinitis pigmentosa