Single-Cell Transcriptomic Profiling in Inherited Retinal Degeneration Reveals Distinct Metabolic Pathways in Rod and Cone Photoreceptors

The cellular mechanisms underlying hereditary photoreceptor degeneration are still poorly understood. The aim of this study was to systematically map the transcriptional changes that occur in the degenerating mouse retina at the single cell level. To this end, we employed single-cell RNA-sequencing (scRNA-seq) and retinal degeneration-1 (rd1) mice to profile the impact of the disease mutation on the diverse retinal cell types during early post-natal development. The transcriptome data allowed to annotate 43,979 individual cells grouped into 20 distinct clusters. We further characterized cluster-specific metabolic and biological changes in individual cell types. Our results highlight Ca2+-signaling as relevant to hereditary photoreceptor degeneration. Although metabolic reprogramming in retina, known as the ‘Warburg effect’, has been documented, further metabolic changes were noticed in rd1 mice. Such metabolic changes in rd1 mutation was likely regulated through mitogen-activated protein kinase (MAPK) pathway. By combining single-cell transcriptomes and immunofluorescence staining, our study revealed cell type-specific changes in gene expression, as well as interplay between Ca2+-induced cell death and metabolic pathways.


Introduction
Retinitis pigmentosa (RP) relates to rare, genetic disorders that cause progressive vision loss [1]. Typically, RP is characterized by a two-step process in which an initial, primary degeneration of rod photoreceptors is followed by a secondary loss of cone photoreceptors, eventually leading to complete blindness [2]. The disease displays a vast genetic heterogeneity, with causative mutations identified in over 85 genes [3]. One of the most extensively studied RP animal models is the retinal degeneration 1 (rd1) mouse. This mouse strain carries a naturally occurring nonsense mutation in the Pde6b gene, encoding for the β subunit of rod cGMP-phosphodiesterase-6 (PDE6). Similar mutations in the human PDE6B gene have been found in RP patients [4].
The mutation-induced loss-of-function of PDE6 leads to photoreceptor cGMP accumulation [5], which activates cyclic nucleotide-gated (CNG) channels. This in turn leads to Ca 2+ influx and, further downstream, activation of Ca 2+ -dependent calpain-type proteases. In parallel, cGMP-dependent activation of protein kinase G (PKG) is associated with activation of histone deacetylase (HDAC) and poly-ADP-ribose-polymerase (PARP) [6]. cGMP-dependent, non-apoptotic cell death appears to be a common mechanism since it was also identified in a diverse set of animal models carrying mutations in a variety of RP genes [6].
Overall, the retina is characterized by high energetic and metabolic demands [7,8], and these demands may be exacerbated by cGMP-dependent activation of CNG-channels and PKG, which in turn may be linked to the depletion of adenosine 5 -triphosphate (ATP) and nicotinamide adenine dinucleotide (NAD + ). Notably, cGMP-dependent Ca 2+influx increases the demand for ATP-dependent Ca 2+ extrusion [9], while PARP uses NAD + to generate poly-ADP-ribose polymer, making it one of the main consumers of NAD + [10]. Thus, ultimately, RD-disease mutations may cause an energetic collapse of the photoreceptor cell [11].
Therefore, insights into the energy metabolism of the retina are crucial for understanding pathology and devising treatments for retinal diseases. Here, we took advantage of single-cell RNA sequencing (scRNA-seq) to investigate the distinct metabolic pathways underlying retinal degeneration in rod and cone photoreceptors. Our analysis suggests that a Ca 2+ -induced activation of the MAPK pathway enhances oxidative phosphorylation (OXPHOS) in rd1 rod photoreceptors. At the same time, rd1 cones may decrease glycolytic activity. These findings may have major ramifications for future therapy developments.

Morphologic Changes in rd1 Retina
The Pde6b mutation in the rd1 mouse initially affects rod photoreceptors in the outer nuclear layer (ONL). Immunofluorescence labeling for PDE6B confirmed the loss of protein expression in the outer segment of rd1 rods when compared to wild-type (WT) (Figure 1a). Concomitantly, a strong reduction of ONL thickness was observed in rd1 retina, while the TUNEL assay revealed a large number of dying cells in the rd1 ONL from post-natal day (P) 11 onwards (Figure 1b; quantifications in Figure 1d,e). Cone photoreceptors, while free of the disease-causing mutation, suffer from a secondary degeneration at a slower rate [12]. Immunofluorescence staining for cone arrestin was used to assess the number of surviving cone cells (Figure 1c; quantifications in Figure 1f) and showed a partial loss of rd1 cones at P17.

Single-Cell RNA Sequencing Analysis Yields 20 Clusters Corresponding to Eight Retinal Cell Types
Using a droplet-based single-cell RNA-seq (scRNA-seq) platform (10× Genomics), single-cell transcriptome analysis was performed on retinal tissues from WT and rd1 mice at postnatal days 11, 13, and 17. Sequencing data were collected from a total of 43,979 cells, of which 23,068 cells (52%) were from WT retinas and 20,911 cells (48%) from rd1 retinas. For unbiased cell classification, we performed a t-stochastic neighbor embedding (tSNE) analysis using Seurat and identified 20 cell clusters, which were further categorized into eight distinct cell types (Figure 2a,c). These were annotated as follows: rod photoreceptors (clusters 0, 1, 2, 3), cone photoreceptors (cluster 9), bipolar cells (clusters 5, 6, 7, 10-14), amacrine cells (clusters 8, 16), Müller cells (clusters 4,19), horizonal cells (cluster 15), microglia cells (cluster 17), and vascular cells (cluster 18). Clusters corresponding to retinal pigment epithelium (RPE) and retinal ganglion cells were not detected in this study, probably due to their relatively low number in the total retinal cell populations. The abundance of WT vs. rd1 cells at each time point is illustrated in Figure 2b, while the relative abundance of cell types in the WT and rd1 cell population is represented in Figure 2d. The percentages of the different cell types for each genotype and time-point are given in Supplementary Table S1. (f) cone cell number. Data from n = 3 animals per group, expressed as mean ± SD. Statistical significance was assessed using Student's t-test; significance levels were: ** p < 0.01. DAPI (grey) was employed as nuclear counterstain. INL, inner nuclear layer.  We then identified the top 35 most regulated genes across the eight most prominent cell types, each of which was characterized by a set of differentially expressed genes (Figure 3). For instance, the genes Rho, Nr2e3, Nrl, Pdc, and Rp1 showed the highest expression on rod, while Opsn1sw, Opn1mw, Pde6h, Arr3, and Gnat2 were expressed in cones. Genes upregulated in Müller cells included Zfp36l1, Dbi, Apoe, Slc1a3, Sparc, and Pcp2. The Pcp2, Trpm1, Isl1, Grm6, and Trnp1 genes were prominent in bipolar cells. The genes Meg3, C1ql1, Snhg11, Tfap2b, and Pcsk1n were enriched in amacrine cells. C1ql1 and Snhg11 were also highly expressed in horizontal cells, as well as Tfap2b, Psk1n, Scl4a3, Calb1, Sept4, and Tpm3. The genes Ctsd, Ccl4, C1qb, and C1qc ranked high in microglial cells. Pcp2, Trpm1, Meg3, and Lgfbp7 were highly expressed in vascular cells. Additional data for the 20 identified clusters, including the top four markers for each cluster, are presented in Supplementary Figure S1. We then identified the top 35 most regulated genes across the eight most prominent cell types, each of which was characterized by a set of differentially expressed genes ( Figure 3). For instance, the genes Rho, Nr2e3, Nrl, Pdc, and Rp1 showed the highest expression on rod, while Opsn1sw, Opn1mw, Pde6h, Arr3, and Gnat2 were expressed in cones. Genes upregulated in Müller cells included Zfp36l1, Dbi, Apoe, Slc1a3, Sparc, and Pcp2. The Pcp2, Trpm1, Isl1, Grm6, and Trnp1 genes were prominent in bipolar cells. The genes Meg3, C1ql1, Snhg11, Tfap2b, and Pcsk1n were enriched in amacrine cells. C1ql1 and Snhg11 were also highly expressed in horizontal cells, as well as Tfap2b, Psk1n, Scl4a3, Calb1, Sept4, and Tpm3. The genes Ctsd, Ccl4, C1qb, and C1qc ranked high in microglial cells. Pcp2, Trpm1, Meg3, and Lgfbp7 were highly expressed in vascular cells. Additional data for the 20 identified clusters, including the top four markers for each cluster, are presented in Supplementary Figure S1.

Differential Gene Expression Analysis in WT and rd1 Retinal Cell Types
An analysis of differentially expressed genes (DEGs) was first carried out on rods and cones at corresponding ages (Supplementary Dataset S1). A volcano plot of DEGs in rd1 retina at P13 showed 151 up-and 111 downregulated genes compared to WT ( Figure  4a). Functional pathways and networks were analyzed for these DEGs with |avg_log2FC| greater than 0.58. This was followed by a gene ontology (GO) enrichment analysis, which indicated that biological processes (BP) related to visual perception, phototransduction, retinal development, and apoptosis were prominently regulated in rd1 rod photoreceptors ( Figure 4a). The 'cellular response to Ca 2+ ions' was also among the top 10 most regulated pathways. Functional enrichment and interactome analysis using the Metascape web portal [13] produced an interactome network in which each enriched term was represented as a node and where connections between nodes were plotted for Kappa similarities over 0.3 across all 262 gene candidates (Supplementary Figure S2). An enrichment of molecular pathways similar to the GO analysis was also found in the KEGG database, where phototransduction, neurotransmitter transport, and synaptic function were among the most regulated processes.
Since cone photoreceptors were affected by the degeneration later than rods (cf. Figure 1c), the corresponding cone pathway and network analysis was performed for the P17 time-point (Figure 4b). A volcano plot for DEGs in P17 cones revealed 219 up-and 154 downregulated genes.
In both GO and KEGG based analyses, the pathways most altered in cones were related to visual perception, retinal development, phototransduction, and synaptic functions. Additionally, tau-protein kinase activity, microtubule nucleation, and protein tetramerization were unique to cone cells as identified by GO pathway analysis. Interestingly, the cGMP-PKG signaling pathway was also related to cone degeneration at P17. Metascape network analysis revealed, among other things, a close interaction between visual perception, retinal development, phototransduction, and photoreceptor cell maintenance (Supplementary Figure S2).

Differential Gene Expression Analysis in WT and rd1 Retinal Cell Types
An analysis of differentially expressed genes (DEGs) was first carried out on rods and cones at corresponding ages (Supplementary Dataset S1). A volcano plot of DEGs in rd1 retina at P13 showed 151 up-and 111 downregulated genes compared to WT (Figure 4a). Functional pathways and networks were analyzed for these DEGs with |avg_log 2 FC| greater than 0.58. This was followed by a gene ontology (GO) enrichment analysis, which indicated that biological processes (BP) related to visual perception, phototransduction, retinal development, and apoptosis were prominently regulated in rd1 rod photoreceptors ( Figure 4a). The 'cellular response to Ca 2+ ions' was also among the top 10 most regulated pathways. Functional enrichment and interactome analysis using the Metascape web portal [13] produced an interactome network in which each enriched term was represented as a node and where connections between nodes were plotted for Kappa similarities over 0.3 across all 262 gene candidates (Supplementary Figure S2). An enrichment of molecular pathways similar to the GO analysis was also found in the KEGG database, where phototransduction, neurotransmitter transport, and synaptic function were among the most regulated processes.
Since cone photoreceptors were affected by the degeneration later than rods (cf. Figure 1c), the corresponding cone pathway and network analysis was performed for the P17 time-point ( Figure 4b). A volcano plot for DEGs in P17 cones revealed 219 up-and 154 downregulated genes.
In both GO and KEGG based analyses, the pathways most altered in cones were related to visual perception, retinal development, phototransduction, and synaptic functions. Additionally, tau-protein kinase activity, microtubule nucleation, and protein tetramerization were unique to cone cells as identified by GO pathway analysis. Interestingly, the cGMP-PKG signaling pathway was also related to cone degeneration at P17. Metascape network analysis revealed, among other things, a close interaction between visual perception, retinal development, phototransduction, and photoreceptor cell maintenance (Supplementary Figure S2). Because of their tight connection to photoreceptors, the pathway analysis was extended to Müller glial cells (Figure 4c). The GO analysis of Müller cells also revealed pathways related to visual perception and retinal development. Furthermore, oxidative phosphorylation, ATP synthesis, and nitric oxide signaling pathways were enriched. A further analysis of pathways changed in rd1 amacrine and horizontal cells is provided in Supplementary Figure S3.
Overall, the analysis of functional pathways and networks indicated that Ca 2+ -signaling and alterations in cellular metabolism might be connected with retinal degeneration and cell death.

Transcriptional Changes Related to Cell Death in rd1 Mutant Retina
Previous research had related rd1 mouse photoreceptor degeneration to non-apoptotic mechanisms that involved the activity of CNG channels, PKG, HDAC, and PARP [6]. We therefore assessed the expression of corresponding genes in rods and cones, starting with that of the various PDE6 genes ( Figure 5). While Pde6b was clearly downregulated in rd1 rods, interestingly, both Pde6a and Pde6g were transiently elevated at P13. Pde6h, which codes for the cone-specific inhibitory subunit of PDE6, was found to be highly up-regulated in cones at P13.
Whatever the case, loss of PDE6 activity produces an accumulation of cGMP within photoreceptors [14,15], keeping CNG channels open. The rod CNG channel is composed of three CNGA1 subunits and one CNGB1 subunit [16]. From our results, both Cnga1 and Cngb1 were significantly upregulated at P11 and P13, but strongly reduced at P17. Cnga3 and Cngb3, which encode the αand β-subunit of the CNG channel in cones, were both down-regulated at P11 and P17.
PKG is another key effector of cGMP-signaling and encoded by the Prkg1 and Prkg2 genes [17]. Prkg1 transcription was upregulated at P17 in both rods and cones, while Prkg2 was not significantly changed.
Downstream of CNG channel and PKG activity, PARP and HDAC have been shown to contribute to rd1 photoreceptor cell death [18]. We found both Parp1 and Parp2 to be downregulated, especially on cones. In rods, Hdac1 was found to be significantly increased at P17, indicating a possible role in the final phase of rod cell death. Both Hdac2 and Hdac3 were essentially downregulated in rd1 rods and cones.
Apart from CNG channels, an important regulator of intracellular Ca 2+ -levels is the Na + /Ca 2+ exchanger-1 (NCX1) encoded for by the Slc8a1 gene. This gene was increasingly upregulated in rods at P13 and P17 but remained unchanged on cones.
The CREB1 and CREM transcription factors together function as a central hub that regulates transcription in response to various stressors, metabolic changes, and developmental signals [19,20]. Importantly, they are targets for PKG phosphorylation and may serve as transducers of cGMP-signaling [21]. A downregulation of Creb1 was seen both in rods and cones at P11 and P17, whereas Crem levels were only increased in rods at P17. We further performed IF on CREB1 and CREM on WT and rd1 retina, at P11, P13, and P17, albeit without finding obvious changes in protein expression (Supplementary Figure S4).
Surprisingly, the gene Aifm1 which encodes the mitochondrial protein apoptosis inducing factor (AIF), an important regulator of programmed cell death [22], did not exhibit any significant transcriptional changes in rods. However, a downregulation of the Aifm1 gene was observed in cone photoreceptors at P17. IF for AIF protein strongly labeled mitochondria-containing structures, such as the photoreceptor inner segments and synapses. However, no significant changes in AIF protein expression were seen between WT and rd1 samples (Supplementary Figure S4).  Hdac1, Hdac2, Hdac3), Na + /Ca 2+ exchanger 1 (NCX1; Slc8a1), cAMP response elementbinding protein 1 (Creb1) and cAMP response element modulator (Crem), and apoptosis inducing factor 1 (Aifm1). The x−axis indicates postnatal day (P), and y−axis depicts average log2 fold change; positive values indicating higher expression in rd1 rods/cones. Bars were color−coded (red for rods; blue for cones) for p−values < 0.05.

Photoreceptors in rd1 Mutant Retina Undergo a Metabolic Switch
Several previous studies indicated a disturbance in energy metabolism during photoreceptor degeneration [23,24], prompting us to examine transcriptional changes on related pathways.

Photoreceptors in rd1 Mutant Retina Undergo a Metabolic Switch
Several previous studies indicated a disturbance in energy metabolism during photoreceptor degeneration [23,24], prompting us to examine transcriptional changes on related pathways.
Using gene set enrichment analysis (GSEA) [25] and the molecular signatures data base (MSigDB) gene sets, we investigated 'glycolysis', 'tricarboxylic acid (TCA) cycle' and 'oxidative phosphorylation (OXPHOS)', as well as the cell death-related pathways 'apoptosis' and 'DNA repair'. Between P11 and P17, all five pathways displayed significant changes in rods and/or cones (Figure 6a). Notably in rods, the energy metabolism-related pathways TCA cycle, OXPHOS, and glycolysis were significantly altered.
To determine how energy metabolism may have changed during rd1 degeneration, we examined sets of individual genes corresponding to each of these three pathways in more detail (Figure 6b). Remarkably, most genes related to glycolysis were strongly upregulated in rods between P11 and P17, while in cones glycolysis genes were downregulated at P17. No clear trend for regulation of TCA cycle-related genes was apparent for either rods or cones. However, genes related to OXPHOS were strongly upregulated in rods at P17, while they were strongly down-regulated in cones at the same time-point. These results indicated a switch in metabolism towards an increased ATP-production in rods. In cones, however, the downregulation of energy metabolism related genes suggested a decreased ATP production at P17, i.e., at the onset of cone degeneration. Normalized average expression of genes involved in glycolysis, TCA cycle, and OXPHOS on photoreceptors were shown in Supplementary Figure S5.  To verify the observed transcriptomic regulation, we performed immunofluorescence staining (IF) for key enzymes of the TCA cycle, OXPHOS, and glycolysis. Citrate synthase (CS) catalyzes the first step of the TCA cycle. As shown by IF, CS was expressed strongly in photoreceptor inner segments and synapses, and in rd1 retina, CS protein expression appeared to weaken as the degeneration progressed from P11 to P17 (Figure 6c). Pyruvate kinase M2 (PKM2) is a glycolytic enzyme that catalyzes the conversion of phosphoenolpyruvate to pyruvate. PKM2 IF was found on photoreceptor inner segments and its expression also seemed to decrease during rd1 degeneration. ATP synthase gamma is a subunit of ATP synthase, a critical enzyme for oxidative phosphorylation and mitochondrial ATP production [26]. It is encoded for by the Atp5c1 gene. ATP synthase gamma was localized prominently on photoreceptor synapses, and in rd1 retina, its expression decreased strongly from P11 to P17.

The MAPK Signaling Pathway Coordinates Energy Metabolism and Cell Death
The mitogen-activated protein kinase (MAPK) pathway is Ca 2+ -dependent [27] and has been implicated in controlling cellular metabolism [28]. Hence, the MAPK pathway could potentially serve as an intermediary between excessive cGMP-signaling and alterations in energy metabolism. We therefore assessed transcriptional changes in MAPK-pathwayrelated genes, identifying significant changes in 46 rod and six cone genes (Figure 7).

Discussion
In this study, we have combined scRNA-seq, immunofluorescence, and cell death detection to gauge the molecular pathways underlying rd1 phenotype. The scRNA-seq dataset annotated 43,979 individual cells grouped into eight distinct retinal cell types and confirmed the key role of cGMP-and Ca 2+ -signaling pathways in rd1 photoreceptor de- Activation of the MAPK signaling pathway typically follows a three-tier kinase module in which a MAP3K phosphorylates and activates a MAP2K, which in turn phosphorylates and activates a MAPK [29]. Among the significantly changed genes, Braf, Taok3, Nlk, Map3k1, Map3k12, Map4k3, and Mapk9 belong to the three layers of activating kinases. Additionally, Fos, Nr4a1, Jun, Jund, and Atf4 are transcription factors involved in MAPK signaling. A significant regulation was also observed for a group of genes that code for voltage-gated Ca 2+ -channel (VGCC) subunits, including Cacna1d, Cacnb2, Cacna2d1, Cacna2d2, and Cacna2d4. This could imply that photoreceptor Ca 2+ -levels may be influenced not only by CNG-channel activity but also by VGCCs as also suggested by recent pharmacological studies [30,31]. Note that of the six genes significantly regulated in cone photoreceptors, five also appear as regulated in rods (Figure 7). At any rate, our data suggests the MAPK pathway as an important regulator of intracellular signaling during rod and cone degeneration and consequently also as a potential target for therapeutic interventions in RD-type diseases.

Discussion
In this study, we have combined scRNA-seq, immunofluorescence, and cell death detection to gauge the molecular pathways underlying rd1 phenotype. The scRNA-seq dataset annotated 43,979 individual cells grouped into eight distinct retinal cell types and confirmed the key role of cGMP-and Ca 2+ -signaling pathways in rd1 photoreceptor degeneration. Molecular profiling of rods and cones indicated a shift from glycolysis toward TCA cycle and OXPHOS activity, likely reflecting increased energy demand. Moreover, our analysis suggested that the MAPK pathway may act as an intermediary between cGMPand Ca 2+ -signaling on the one hand and cellular metabolism on the other hand.

cGMP-and Ca 2+ -Signaling in rd1 Retinal Degeneration
The loss of rod photoreceptors in rd1 mice was previously found to depend on a non-apoptotic cell death mechanism triggered by high intracellular cGMP-levels [6,32]. In photoreceptors, high cGMP likely activates the prototypic effectors PKG and CNG channels, leading, among other things, to increased Na + -and Ca 2+ -influx [33,34]. The continuous depolarization effected by CNG-channel activity causes a sustained activation of VGCCs leading to more Ca 2+ influx [30,35]. The Na + /Ca 2+ exchanger (NCX) type antiporters encoded by the Slc8 gene family utilize the Na + gradient to extrude Ca 2+ [36,37]. The transcriptional upregulation of Slc8a1 may represent an attempt to counterbalance Ca 2+ overload. However, the Na + gradient required for NCX to function is established to a large extent by the ATP-driven Na+/K+ exchanger (NKX), so that NCX-dependent Ca 2+extrusion will likely place an additional burden on photoreceptor energy metabolism [38].
The link between excessive cGMP-signaling, PKG activity, and cell death has been well established for the rd1 mouse and other RD animal models [6,14,39]. However, at present, it is still unclear which PKG isoform may be responsible for photoreceptor death. Here, we observed an upregulation of Prkg1 but not Prkg2 in rod photoreceptors during the critical degeneration phase. Remarkably, in cones increased Prkg1 gene transcription was observed only at P17, i.e., at the beginning of cone degeneration [11,40], indicating that also cone death may be triggered by Prkg1 overactivation.
In retinal degeneration PKG activity is associated with an overactivation of PARP and HDAC [30,41,42]. PARPs are a superfamily of ADP-ribosylating enzymes and are known to play important roles in DNA repair and the maintenance of genome integrity [43][44][45]. The main function of HDACs consists in removing acetyl groups from DNA-binding histone proteins, which is generally associated with a decrease in chromatin accessibility for transcription factors and hence represses gene expression [46]. Although Parp1 and Parp2, as well as Hdac2 and Hdac3, did not show significant changes in P17 rods, it is important to note that activity of PARP and HDAC relies mostly on post-transcriptional regulation and may therefore not be adequately resolved with transcription-based analysis [6,47].

Metabolic Responses of Degeneration Retina
The retina and especially the photoreceptors are characterized by a very high energy expenditure, most of which may be due to the active transport of ions against their concentration and electrical gradients [8]. In the rd1 condition, higher levels of Ca 2+ -influx entail increased ATP-consumption for Ca 2+ extrusion [11,32]. Besides, increased PARP activity consumes large amounts of NAD + [48] and low levels of this critical electron acceptor may further impair mitochondrial ATP production [49].
Paradoxically, even in the presence of oxygen the retina converts a large fraction of its available glucose to lactate rather than oxidizing it completely to carbon dioxide. This phenomenon was first observed by Otto Warburg and is known as 'aerobic glycolysis' or 'Warburg effect' [50]. Per molecule of glucose, aerobic glycolysis produces only two ATP molecules, while the TCA cycle coupled to OXPHOS can generate up to 36 ATP molecules. Thus, from an energetic point of view mitochondrial oxidative metabolism is by for more efficient than cytosolic glycolytic metabolism.
To date, it is not clear why the retina uses the seemingly inefficient aerobic glycolysis to generate ATP, although it has been speculated that glycolytic intermediates, including precursors of nucleic acids, lipids, and amino acids, may be required for retinal anabolic activity [51]. Perhaps even more surprising is that our study suggests that during retinal degeneration glycolytic activity may decrease while OXPHOS appears to increase. Future studies may reveal whether this increases ATP-production or whether it reflects a decrease in the production of glycolytic intermediates. Moreover, our data indicates that this switch in energy metabolism may be regulated by the MAPK pathway.

MAPK Pathway Regulates Crosstalk between Ca 2+ and Energy Metabolism
The MAPK cascade has been shown to play a key role in the regulation of cell proliferation, differentiation, and death [52]. In mammals, MAPKs are divided into three subfamilies: extracellular signal-regulated kinases (ERKs), Jun-N-terminal kinases (JNKs), and p38 kinases [53]. MAPK signaling is activated by a three-tier kinase module that consists of a MAP3K phosphorylating and activating a MAP2K, which in turn activates a MAPK [29]. According to our results, ERK and JNK pathways were upregulated during retinal degeneration. While previously Ca 2+ overload in photoreceptor leads was found to cause activation of Ca 2+ -dependent calpain-type proteases [54,55], high Ca 2+ may additionally activate MAPKs and ERKs directly or indirectly through protein kinase C. The MAPK signaling pathway has recently been found as a key regulator of the Warburg effect in cancer and in metabolic reprogramming [28]. Notably, the MAPK cascade promotes aerobic glycolysis through PKM2 phosphorylation, as well as regulating the expression and activity of transcription factors that directly control the expression of glycolytic enzymes, including Myc-associated factor X (MAX) and c-MYC [56].

Materials and Methods
Animals. C3H HeA Pde6b rd1 / rd1 (rd1) and congenic C3H HeA Pde6b +/+ wild-type (WT) mice were used, two mouse lines that were originally created in the lab of Somes Sanyal at the University of Rotterdam [57]. All efforts were made to minimize the number of animals used and their suffering, notably through the use of in vitro experimentation (see below) and the use of both retinae from the same animal. Animals were housed in the specified pathogen-free facility of the Tübingen Institute for Ophthalmic Research, under standard white cyclic lighting in type-2 long cages with a maximum of five adults per cage. They had free access to food and water and were used irrespective of sex. Protocols compliant with the German law on animal protection were reviewed and approved by the 'Einrichtung fur Tierschutz, Tierärztlichen Dienst und Labortierkunde' of the University of Tübingen (AK 02/19 M, notice acc. to §4 German law on animal protection) and were following the association for research in vision and ophthalmology (ARVO) statement for the use of animals in vision research.
Immunofluorescence. After enucleation, the eyes were fixed at room temperature in 4% paraformaldehyde (PFA) for 45 min and then cryoprotected in 10%, 20%, and 30% sucrose solutions. Then, they were embedded in O.C.T. (Tissue-Tek ® O.C.T. TM compound, Sakura Finetek, Umkirch, Germany), and flash-frozen on liquid nitrogen. The embedded eyes were sectioned to a thickness of 12 µm in a cryotome (Thermo Scientific, CryoStar NX50, Walldorf, Germany) and collected on Superfrost glass slides (R. Langenbrinck GmbH, SuperFrost ® plus, Emmendingen, Germany) and stored at −20 • C until further use. Fixed slides were thawed and dried at 37 • C for 30 min and rehydrated for 10 min in PBS at room temperature (RT; 21 • C). For immunofluorescent labeling, the slides were incubated with blocking solution (10% normal goat serum, 1% bovine serum albumin in 0.3% PBS-Triton X 100) for 1 h at RT. The primary antibodies were diluted (see Table 1) in blocking solution and incubated at 4 • C overnight. The slides were then washed with PBS, three times for 10 min each. Subsequently, a corresponding secondary antibody (see Table 1), diluted in PBS, was applied and incubated for 1 h at RT. Lastly, the slides were washed with PBS and coverslipped using Vectashield with DAPI (Vector, Burlingame, CA, USA). TUNEL assay. Terminal deoxynucleotidyl transferase dUTP nick end labelling (TUNEL) assay was performed using an in situ cell death detection kit (Fluorescein or TMR; Roche Diagnostics GmbH, Mannheim, Germany). Fixed slides were dried at 37 • C for 30 min and rehydrated in phosphate-buffered saline (PBS) solution at RT, for 15 min. Afterward, the slides were treated with 21.5 µg/mL proteinase K (43 U/mL; Cat. No.: 03115879001; Roche Diagnostics, Mannheim, Germany) in TRIS buffer (10 mM TRIS-HCL, pH 7.4) at 37 • C for 5 min. The slides were then washed with TRIS buffer three times for 5 min each. Subsequently, the slides were placed in ethanol-acetic acid mixture (70:30) at −20 • C for 5 min followed by three washes in TRIS buffer and incubation in blocking solution (10% normal goat serum, 1% bovine serum albumin, 1% fish gelatin in 0.1% PBS-Triton X100) for 1 h at RT. Lastly, the slides were placed in the TUNEL solution (labeling with either fluorescein or tetra-methyl-rhodamine) in 37 • C for 1 h and cover slipped using Vectashield with DAPI (Vector; Burlingame, CA, USA).
Single-cell RNA-seq (scRNA-seq) and bioinformatics analysis. scRNA-seq were performed by Jiayin Biotechnology, Ltd. (Shanghai, China). Retinal cellular suspensions (43,979 cells) were loaded on a 10× Genomics Chromium Single Cell instrument (10× Genomics, Shanghai, China) to generate single-cell Gel Beads in Emulsion (GEMs). Barcoded sequencing libraries were conducted following the instruction manual of the Chromium Single Cell 3 Reagent Kits v3 (10× Genomics). Following the library preparation, the sequencing was performed with paired-end sequencing of 150nt each end on one lane of Illumina NovaSeq 6000 (Illumina, San Diego, CA, USA) per sample. Retinal scRNA-seq analyses were performed using the Seurat package in R 14 (R Core Team (2022), R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. URL https://www.R-project.org/ (accessed on 11 May 2021)). Briefly, cells with a significant number of outlier genes (potential polysomes) and a high percentage of mitochondrial genes (potential dead cells) were excluded from further analysis using the "Filter Cells" function [58,59]. Log normalization was used to normalize gene expression. Principal component analysis (PCA) was then performed to reduce the dimensionality of the dataset using t-SNE/UMAP dimensionality reduction. Seurat was used to cluster cells based on the PCA scores. For every single cluster, differentially expressed genes (DEGs) were identified using the "Find All Markers" function in the Seurat package, and the screening threshold was set to |avg_log 2 FC| > 0.58 and p < 0.05. Statistical comparison was performed using the Wilcoxon rank sum test. The identified genes were functionally and taxonomically annotated with major databases, including the Gene Ontology (GO) enrichment analysis, the Kyoto Encyclopedia of Genes and Genomes (KEGG) and Metascape (https://metascape.org/ (accessed on 30 August 2022)), to understand the functional properties and classifications of different genes. The results were shown by the R software (RStudio Team (2022). RStudio: Integrated Development for R. RStudio, PBC, Boston, MA URL http://www.rstudio.com/ (accessed on 30 August 2022)).

Conclusions
The single-cell analysis described in this report provides transcriptional signatures for the diverse populations of neuronal and glial cells in the rd1 mouse model for retinal degeneration. We focused on alterations in metabolic pathways in rod and cone photoreceptors in the critical time period from P11 to P17 during which most of the rod degeneration occurs. The identification of a metabolic switch likely related to altered activity of the MAPK-pathway highlights new possibilities for treatments designed to reprogram cellular metabolism so as to promote cell survival and recovery. Remarkably, the metabolism switch observed in rd1 mutant rods may extend to genetically intact cones, suggesting that very similar treatments could rescue both rods and cones. On a broader level, our study provides further in-depth insights into the pathological mechanisms of rod-cone dystrophy.  Data Availability Statement: Publicly available datasets were analyzed in this study. This data can be found here: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE212183.