Transcriptional Profiling Identifies Upregulation of Neuroprotective Pathways in Retinitis Pigmentosa

Hereditary retinal degenerations like retinitis pigmentosa (RP) are among the leading causes of blindness in younger patients. To enable in vivo investigation of cellular and molecular mechanisms responsible for photoreceptor cell death and to allow testing of therapeutic strategies that could prevent retinal degeneration, animal models have been created. In this study, we deeply characterized the transcriptional profile of mice carrying the transgene rhodopsin V20G/P23H/P27L (VPP), which is a model for autosomal dominant RP. We examined the degree of photoreceptor degeneration and studied the impact of the VPP transgene-induced retinal degeneration on the transcriptome level of the retina using next generation RNA sequencing (RNASeq) analyses followed by weighted correlation network analysis (WGCNA). We furthermore identified cellular subpopulations responsible for some of the observed dysregulations using in situ hybridizations, immunofluorescence staining, and 3D reconstruction. Using RNASeq analysis, we identified 9256 dysregulated genes and six significantly associated gene modules in the subsequently performed WGCNA. Gene ontology enrichment showed, among others, dysregulation of genes involved in TGF-β regulated extracellular matrix organization, the (ocular) immune system/response, and cellular homeostasis. Moreover, heatmaps confirmed clustering of significantly dysregulated genes coding for components of the TGF-β, G-protein activated, and VEGF signaling pathway. 3D reconstructions of immunostained/in situ hybridized sections revealed retinal neurons and Müller cells as the major cellular population expressing representative components of these signaling pathways. The predominant effect of VPP-induced photoreceptor degeneration pointed towards induction of neuroinflammation and the upregulation of neuroprotective pathways like TGF-β, G-protein activated, and VEGF signaling. Thus, modulation of these processes and signaling pathways might represent new therapeutic options to delay the degeneration of photoreceptors in diseases like RP.


Introduction
Retinitis pigmentosa (RP) is a hereditary form of retinal degeneration that results from mutations in any one of more than 70 known susceptibility genes [1,2]. Quite intriguingly, RP is considered as one of the most common hereditary diseases associated with mutations in core splicing proteins resulting in the altered regulation of gene expression [2]. Even though RP is considered a rare genetic disorder, it is still among the major causes

Transcriptional Alterations in VPP Retinae: RNAseq and Weighted Correlation Network Analysis (WGCNA)
Subsequently, we applied next generation RNA sequencing (RNAseq) analyses to investigate the impact of VPP transgene expression and concomitant photoreceptor degeneration on the transcriptome of the retina in three-month-old VPP and control animals. Out of the total of 54,532 genes in the Ensembl gene annotation for mouse (Mus musculus GRCm38 v. 94), we found 30,796 genes to be expressed in the retina, of which 9256 were dysregulated (4636 down-and 4620 upregulated, Figure 1A, cut off criteria: Benjamin-Hochberg adjusted p-value (p adj ) < 0.05). The top 30 dysregulated genes are shown in Supplementary Table S2. Among others, genes regulating processes in neurotransmission like histidine decarboxylase (Hdc); galactosidase beta 1 like 3 (Glb1l3), which is associated with Leber's congenital amaurosis; as well as serine protease 56 (Prss56), which was reported to be involved in eye development, were significantly downregulated. Genes controlling scar formation, such as fibrinogen-like 2 (Fgl2) as well as apoptosis, e.g., caspase1 (Casp1) or Bcl2-interacting killer (apoptosis-inducing) (Bik), were significantly higher expressed in VPP retinae. Furthermore, we found upregulation of quite a considerable number of genes associated with inflammatory or immune response functions such as C-X-C motif chemokine ligand 13 (Cxcl13), glial fibrillary acidic protein (Gfap), T-cell receptor T3 gamma chain (Cd3g), chemokine (C-C motif) ligand 5 (Ccl5), and C-C motif chemokine ligand 2 (Ccl2), as well as factors associated with the complement cascade like complement component factor i (Cfi), complement factor C4B, and the Serping1 gene. Gene ontology enrichment showed, among others, involvement of TGF-β regulated extracellular matrix organization, response to cytokine stimuli, and disturbed cellular homeostasis in VPP retinae (Table 1). Photoreceptor loss was indicated by a lower expression level of rhodopsin signaling pathway components (Table 1). Moreover, we performed weighted gene correlation network analysis (WGCNA) to identify genotype-specific patterns of dysregulation, upstream regulators, and involved signaling pathways. WGCNA clusters co-regulated genes into modules based on their similarity of expression. As this approach does not rely on the traditional dysregulation analysis and the problem of correction for multiple comparisons, more subtle changes and patterns can be identified. In addition, biological key players, e.g., regulatory proteins driving a certain pathway, for a given module can be found by the intra-module analysis. The topology overlay matrix, which represents the co-regulation of expression, for VPP and control animals, as well as the identified modules (clusters of co-regulated genes as shown by their colors in Figure 1B) and their correlation of expression with the genotype, are illustrated in Figure 1B. The analysis identified six significantly associated modules (three positively correlated with the genotype, i.e., higher expressed in the VPP animals (Pos1, 2, 3) and three negatively correlated, i.e., lower expressed in the VPP animals (Neg1, 2, 3) ( Figure 1D,E and Supplementary Figure S2A-D)). The Pos1 module was significantly enriched for genes involved in cellular transport and differentiation ( Figure 1D and Table 2). In the Pos2 module, we found significant clustering of genes regulating necroptosis, protein transport, organelle organization, cellular homeostasis, and the ribosome (Supplementary Figure S2A and Table 2). The Pos3 module (Supplementary Figure S2B and Table 2) showed enrichment for bone morphogenetic protein (BMP) signaling pathway components. In the Neg1 module, we observed significant clustering of genes regulating cellular component organization and protein kinase activity ( Figure 1E and Table 2). The Neg2 module (Supplementary Figure S2C and Table 2) was enriched for genes involved in nucleotide homeostasis and the Neg3 module showed significant enrichment for genes encoding for ribosomal proteins and transport proteins (Supplementary Figure S2D and Table 2). Table 1. Enrichment analysis for dysregulated genes derived from the VPP RNAseq analysis. Enriched pathways and potential upstream regulators were predicted using the indicated databases. For gene ontology enrichment, only the top five non-redundant significantly enriched biological process terms are shown. The numbers following the terms are the combined score as calculated by Enrichr. Only terms with a combined score >5 were considered.
Thus, we aimed to investigate the impact of these signaling pathways (TGF-β-, Gprotein activated-, and VEGF-signaling) in the VPP model in detail. Unsupervised hierarchical clustering of the samples that we generated on basis of the Reactome pathway database [16] for the VEGF (Figure 2A), TGF-β ( Figure 2B) and G-protein mediated signaling pathways ( Figure 2C) demonstrated a perfect separation of the genotypes, highlighting the dysregulation of these pathways in the VPP animals. Furthermore, k-mer analysis (three k-mer groups indicated by numbers on the left side of each heatmap) showed clusters of tightly co-regulated genes. We highlighted some genes of particular interest (e.g., which are known to be involved in neuroprotective or immune modulating processes or to be involved in regulatory functions) in each pathway on the right side of each heatmap. The heatmaps including the full labelling are shown in Supplementary Figure S3. To further analyze sub-groups of dysregulated pathways, we transformed the Reactome pathways into functional interaction networks ( Figure 2D-F). We colored the genes of each network according to their dysregulation: white indicates no significant regulation, red genes were up-and blue genes were significantly down-regulated, respectively. The size of each node corresponds to the log2-fold change of regulation. The fully labeled networks are shown in Supplementary Figure S4. This analysis identified distinct sub-clusters of dysregulated genes, e.g., endothelin 2 (Edn2) and endothelin receptor b (Ednrb) in the G-protein activated signaling pathway network ( Figure 2E and Supplementary Figure S4C) or Tgfbr2 in the TGF-β family signaling network ( Figure 2D and Supplementary Figure S4B) and Vegfr2/kinase insert domain receptor (Kdr) in the VEGF signaling network ( Figure 2F and Supplementary Figure S4A).
In addition to the quantitative information from the RNAseq data derived from whole retinae, and to further validate our results, we next performed mRNA in situ hybridization and/or immunofluorescence staining to identify cell types expressing transcripts of interest. Using a specific probe against Edn2, we detected Edn2 most prominently in the ONL in the control retina. However, we also observed distinct Edn2 signals in the INL and some rather sparse signals in the GCL ( Figure 3A). VPP retinae showed significantly higher Edn2 expression in the RNAseq data (22.12-fold, p adj = 8.68 × 10 −87 ; Figure 3B, Supplementary Table S2, Supplementary Figure S1G) and accordingly showed a marked increase of Edn2 signals in particular in the degenerating ONL and in the INL, but the signals in the GCL remained sparse ( Figure 3A). To identify specific cell types expressing Edn2, we combined in situ hybridization with immunofluorescence staining. We used glial fibrillary acidic protein (GFAP) as a marker for astrocytes and reactive Müller cells. GFAP is a major intermediate filament particularly expressed in astrocytes [17] and upregulated in response to retinal trauma in astrocytes and Müller cells [17,18]. Moreover, immunofluorescence staining for glutamine synthetase (GS) was used to label Müller cells [19]. When performing GFAP/GS/Edn2 co-labelling, we observed the characteristic GFAP staining pattern for astrocytes in the nerve fiber layer for both groups (Supplementary Figure S5A,B). Consistent with the observed upregulation of Gfap in VPP retinae as determined by RNAseq analysis (17.05-fold, p adj = 3.87 × 10 −78 ), immunofluorescence staining showed an increased GFAP signal intensity, indicating an enhanced protein expression, in the nerve fiber layer and an additional stripe-like staining pattern stretching through the retina, which represents the characteristic morphological appearance of reactive Müller cells. Of note, unlike apoptosis and ONL thinning (Supplementary Figure S1A,D,E), which were both more prominent in the central retina, we could not detect a difference in GFAP reactivity between the central (Supplementary Figure S5A  , and G alpha (q) signaling (C). For each heatmap, the genotypes separate perfectly, as indicated by the unsupervised clustering above the heatmaps. Colors (red: upregulated, blue: downregulated) represent the deviation of the mean expression for each gene, independent of genotype. K-mer analysis into three groups revealed clusters of tightly co-regulated genes. Some interesting genes (e.g., neuroprotective or immune modulating function, directly involved in the intracellular signaling) are highlighted on the right. To further visualize sub-groups of pathways that were dysregulated, we converted the Reactome pathways into functional interaction networks. For each network, genes were colored according to their dysregulation state: white-not significantly dysregulated; red-significantly upregulated; and blue-significantly downregulated. The size of the nodes corresponds to the log2-fold change of regulation. The network for TGF-beta signaling is shown in (D), G alpha (q) signaling is shown in (E), and VEGF signaling is shown in (F); R21742-61 = RNAseq sample number. Ednrb mRNA in situ hybridization showed specific signals in the INL and ONL and some defined signals in the GCL in control retinae ( Figure 3E). We furthermore observed some Ednrb expression in GFAP-positive astrocytes ( Figure 3G and Supplementary Figure  S6E) and in GS-positive resting Müller cells ( Figure 3H and Supplementary Figure S6F). In VPP retinae, Ednrb was significantly upregulated (1.29-fold, p adj = 0.007; Figure 3F) in the RNAseq data (Supplementary Table S2). In the 3D reconstruction of the Ednrb/GFAP/GS labelling, we detected pronounced Ednrb signals in the INL and ONL that overlapped to some extent with GFAP-/GS-positive reactive Müller cells ( Figure 3E,H and Supplementary Figure S6F). Yet, Ednrb signals were visible in the INL and ONL that did not overlap with GFAP/ GS, indicating its further expression in neuronal cells.
Tgfbr2 was significantly upregulated (2.23-fold, p adj = 2.18 × 10 −23 ; Figure 4B) in VPP retinae in the RNAseq data (Supplementary Table S2). To further supplement the quantitative information and potentially identify cell types in which it is upregulated, we performed Tgfbr2 in situ hybridization. Control retinae showed distinct signals in the INL and ONL and some scattered punctae in the GCL ( Figure 4A and Supplementary Figure S6G,H). In VPP retinae, the number of Tgfbr2 punctae increased in the INL and ONL. 3D reconstruction of immunofluorescence co-labelling confirmed its expression in only some isolated GFAP-positive astrocytes ( Figure 4C and Supplementary Figure S6G) and its association with resting, GS-positive, and reactive GFAP-/GS-positive Müller cells, respectively ( Figure 4D and Supplementary Figure S6H). Yet, we also observed Tgfbr2 in situ hybridization in the neuronal cell layers of the retina, in particular in the INL and ONL, that did not overlap with GFAP-/GS-positive Müller cells, indicating additional expression in neuronal cells ( Figure 4A).
Vegfr2 mRNA in situ hybridization in control retinae showed numerous signals in the INL that partly overlapped with GS-positive resting Müller cells ( Figure 4E and Supplementary Figure S6J). Moreover, we detected Vegfr2 mRNA in the ONL and isolated signals in the GCL that overlapped to some extent with GFAP-positive astrocytes ( Figure 4E,G and Supplementary Figure S6I). Our RNAseq analysis showed Vegfr2 to be significantly upregulated in VPP retinae (2.00-fold, p adj = 1.40 × 10 −41 ; Figure 4F; Supplementary  Table S2) and, accordingly, Vegfr2 in situ hybridization showed an increase in expression in the INL and ONL ( Figure 4E). 3D reconstruction of co-labelling showed its association with and expression in GFAP-/GS-positive reactive Müller cells ( Figure 4H and Supplementary Figure S6J). Moreover, we detected Vegfr2 signals in the neuronal layers of the retina, again in particular in the INL and ONL, that did not overlap with GFAP-/GS-positive Müller cells, indicating additional expression in neuronal cells ( Figure 4E).

The Glial Response to Photoreceptor Degeneration in VPP Mice
As we observed a significant upregulation of several genes modulating glial reactivity (e.g., glial fibrillary acidic protein (Gfap), serine protease inhibitor A3N (Serpina3n), lipocalin 2 (Lcn2), we performed qPCR analyses on retinal samples to determine the relative Gfap expression levels (control: 1.00 ± 0.19, VPP: 8.89 ± 0.63, p < 0.001) as well as the expression levels of the microglia/macrophage marker ionized calcium-binding adapter molecule 1 (Iba1, control: 1.00 ± 0.11, VPP: 6.56 ± 1.03, p = 0.015) and the chemokine (C-C motif) ligand 2 (Ccl2, control: 1.00 ± 0.11, VPP: 68.74 ± 5.11, p < 0.001), the latter being reported to stimulate the migration and reactivity of microglia cells [20,21] (Supplementary Figure S5D). In accordance, our RNAseq data showed increased Iba1 (5.51-fold, p adj = 1.12 × 10 −23 expression levels), which we further validated using an anti-IBA1 labeling to visualize myeloid cells, e.g., microglia and recruited macrophages, in the retinae of control and VPP animals (Supplementary Figure S5E,F). In controls, we observed ramified IBA1-positive cells in their typical localization in the nerve fiber layer and the inner (IPL) and outer plexiform layers (OPL). In contrast, in VPP retinae, IBA1-positive cells changed their shape from ramified microglia towards amoeboid, reactive microglia in particular in the OPL, and thus in very close association with the degenerating photoreceptors. Moreover, we observed an accumulation of amoeboid-shaped, IBA1-positive cells in the sub-neuroretinal space in close proximity to the retinal pigment epithelium (RPE) (Supplementary Figure S5E,F). Taken together, these results showed a pronounced reactivity of macro-and microglial cells in response to photoreceptor degeneration. To identify the origin of the significantly elevated Ccl2 expression in VPP retinae and to supplement the quantitative information from the RNAseq (Ccl2: 67.51-fold increase, p adj = 8.05 × 10 −14 ; Figure 5B) and qPCR data (Supplementary Figure S5D), we performed Ccl2 mRNA in situ hybridization on retinal sections combined with immunofluorescence co-labeling of glial cells. In control retinae, we observed a rather low number of Ccl2 punctae in the inner nuclear layer (INL) and the ONL, and a few signals in the retinal ganglion cell layer (GCL) ( Figure 5A and Supplementary Figure S6A,B). The number of Ccl2 punctae was markedly increased in the INL and ONL of the VPP retinae ( Figure 5A). When performing 3D reconstruction of Ccl2/GFAP/GS co-labelled sections, we observed sparse overlap of Ccl2 in GFAP-positive astrocytes ( Figure 5C and Supplementary Figure S6A

Discussion
The present data confirm that the VPP model displays the major phenotypic characteristics of the human disease retinitis pigmentosa. Briefly, we demonstrate, in comprehensive transcriptome-wide analyses of retinae from three-month-old VPP mice, (1) an extensive dysregulation of genes modulating apoptosis, processes in scar formation, and components of the (ocular) immune system or response, respectively; (2) a strong genotype-dependent clustering of genes regulating the VEGF, TGF-β, and G-protein activated signaling pathway; (3) the expression of regulatory genes in neurons as well as resting and reactive glia cells; and (4) a dysregulation of extracellular matrix organization and cellular homeostasis in WGCNA analyses.

The Transcriptional Response to Photoreceptor Degeneration Leads to Increased Expression of Genes Regulating Inflammatory or Immune Response Functions
Neuroinflammation is a common hallmark of the pathogenesis of neurodegenerative diseases like Alzheimer's, Parkinson's, multiple sclerosis, or retinal degenerations [22][23][24]. Following a neurotoxic event, neuronal stress signals mediate reactivity of microglial cells, leading to their proliferation, migration, and the secretion of specific cytokines and chemokines that can exert neurotoxic or neuroprotective effects [23,25]. Sustained reactivity of microglia promotes chronic inflammation and may cause irreversible neuronal cell death [23,26,27]. Thus, the accumulation of reactive IBA1-positive cells in the OPL and in the sub-neuroretinal space in VPP retinae strongly indicates an ongoing neuroinflammatory process. Moreover, in the top 30 dysregulated genes, we found a considerable number of genes associated with inflammatory or immune response functions, respectively. Gene ontology enrichment analyses also pointed towards an upregulation of the cellular response to cytokine stimuli, again indicating an ongoing neuroinflammation. These findings are in accordance with previously published data, which describe upregulation of factors like Lif, Ccl2 (Mcp-1), Ccl28, interleukin-1 (Il-1), complement component 1q (C1q), and complement factor H (CFH) in retinae of genetic mouse models of RP [11,[28][29][30][31][32][33]. Quite intriguingly, microarray data from the retina of one RP patient carrying two mutations in the ABCA4 gene showed, among others, an increased expression of complement system genes (complement factor B, complement C2), several cytokines, and cytokine receptors (IL-6, CXCL10, CXCL2), respectively, indicating a neuroinflammatory process in human RP [24]. We as well as others [20,21] have shown that Ccl2 is expressed in Müller cells and photoreceptors in the healthy retina and upregulated upon retinal damage, contributing to the recruitment of microglia/infiltrative macrophages [21,34]. However, conflicting data exist concerning the exact role of Ccl2 in the context of neurodegeneration. Recently, Joly and colleagues showed that the retinal morphology of double mutant mice expressing the VPP transgene on a Ccl2 null background (VPP; Ccl2 −/− ) did not differ from that of transgene VPP mice on a wildtype background [11]. In contrast, Rutar and colleagues demonstrated that siRNA-mediated knock down of Ccl2 resulted in a significantly lower number of apoptotic photoreceptors in rats after light-induced photoreceptor degeneration and [21,34]. Based on our data, we hypothesize that, in VPP retinae, the elevated Ccl2 expression in Müller cells and in the ONL contributes to the attraction/migration and reactivity of microglial cells in this particular region. In our comprehensive analyses, we furthermore detected a considerable number of genes encoding for components of the complement system, which is part of the innate immune system. Various complement factors have been reported to be upregulated in retinae of human patients suffering from RP or in retinae of mice following genetically or light-induced photoreceptor degeneration in mice [30][31][32][35][36][37][38]. Activation of the complement system promotes microglia/infiltrative macrophages migration and eventually complement activated lysis [31,35]. Still, conflicting data exist regarding the exact role of complement system activation and its impact on photoreceptor degeneration. Mice with a deficiency in complement factor D are protected from light-induced photoreceptor degeneration [30], indicating a detrimental role for photoreceptor survival, but the deficiency in complement component 3 (C3) or complement receptor 3 (CR3) in a genetic mouse model of photoreceptor degeneration increases microglia-mediated neurotoxicity to photoreceptors [38]. Thus, the detailed function of the complement system and its specific role in microglia and Müller cells and its contribution to photoreceptor degeneration has yet to be elucidated. Nevertheless, based on our transcriptome-wide data, we conclude that the significant upregulation of genes associated with inflammatory or immune response functions leads to neuroinflammation in VPP retinae, potentially contributing to the degeneration of photoreceptors.

The Transcriptional Response to Photoreceptor Degeneration Leads to the Upregulation of Neuroprotective Factors and Pathways
We additionally analyzed the impact of the VPP model on neuroprotective pathways like TGF-β, G-protein, and VEGF signaling [39]. Recently, our group and others showed that, in response to retinal injury, endothelin 2 (Edn2) is expressed by photoreceptors concomitant with elevated expression of endothelin receptor b (Ednrb) and Gfap, the latter indicating the reactivity of Müller cells, and an increased expression of Lif and Fgf2 [10][11][12]15,40]. Our RNAseq data, in situ hybridizations, immunofluorescence staining, and qPCR analyses confirmed this observation for the VPP model of retinal degeneration. It is of interest to note that our co-labelling experiments showed expression not only of Ednrb, but also of Tgfbr2 and Vegfr2 in resting and reactive Müller cells, clearly indicating the close interplay of neuronal and glial cells. Furthermore, we previously showed that Edn1, Edn2, Ednra, and Ednrb are significantly upregulated following induced ocular traumata [15]. Yet, in VPP retinae, only Edn2 and Ednrb were upregulated. These inconsistent results might well be explained by the different activation patterns of signaling pathways depending on the actual cause of cell death, i.e., light-induced versus genetically-induced cell death [28].

The Transcriptional Response to Photoreceptor Degeneration Leads to Expression of Pro-Apoptotic Factors and Extracellular Matrix Organization
Our transcriptome analyses also showed upregulation of factors associated with apoptosis (Casp1, Bik) and scar formation such as fibrinogen like 2 (Fgl2) and Tgf-β1. A well-described characteristic of TGF-β signaling is its contribution to wound healing, tissue fibrosis, and scar formation [41,42]. Accordingly, "TGF-β regulated organization of extracellular matrix" was a major hit in our gene ontology analysis and might well be the result of a healing response following photoreceptor degeneration in VPP retinae. We furthermore identified the isoform Tgf-β2 as one of the central hub genes in the WGCNA module Pos1. Gene ontology analysis of the Pos1 module showed enrichment for genes involved in cell differentiation. TGF-β signaling modulates manifold processes, e.g., the regulation of early development, cell-cycle control, and cell differentiation [43][44][45]. Moreover, we recently showed that the deletion of TGF-β signaling results in the development of retinal microaneurysms and choroidal neovascularization [13,46], clearly emphasizing its potential to regulate angiogenic processes. Furthermore, TGF-β signaling contributes to extracellular matrix reorganization in ocular diseases such as primary open angle glaucoma [47][48][49]. However, AAV-mediated delivery of TGF-β1 rescued degenerating cones in mouse models of RP [50] and TGF-β signaling protects retinal neurons from programmed cell death during retinal development [14], thus highlighting its neuroprotective properties [51].
The observed gliosis of astrocytes and Müller cells, as indicated, e.g., by elevated Gfap expression levels and the characteristic stripe-like staining pattern of GFAP in retinal sections, is a typical reaction of neuronal tissue to various neurotoxic insults [18,52] and eventually results in a glial scar. Quite intriguingly, recently published data suggest a major and interactive role of glial cells such as astrocytes contributing to neurodegeneration and their cell-specific targeting resulted in accelerated functional recovery compared with untreated animals [53,54].
The identified dysregulation of genes involved in neuroinflammation, neuroprotection, apoptosis, scar formation, and wound healing and the corresponding WGCNA data are not only of interest for researchers working on retina, but might well be of interest for scientists working with other neuronal tissues.
In the future, using material from mouse models of RP or humanized 3D culture models (retinal organoids) combined with advanced transcriptomic techniques such as imaging-based in situ cell-type identification/mapping combined with single-cell RNAsequencing will allow to create a molecularly annotated and spatially dedicated cell atlas of transcriptional changes related to RP [55,56]. Moreover, the use of high-throughput screening and computer-aided drug design can provide novel insights with the overall aim to find new treatment options for neurodegenerative diseases such as RP [57,58].

Mice
The mice were on a 129 SV background and kept in a 12 h light/dark cycle. Mice of both sexes were used for the experiments. Mice, carrying two floxed Tgfbr2 alleles at chromosome 9 [59], thus representing functional wildtype mice, were crossbred with hemizygous VPP mice [6]. The resulting offspring expressed either wildtype rhodopsin (named as control mice = VPP negative animals) or the VPP transgene (referred to as VPP mice), a rhodopsin mutant with point mutations at positions V20G, P27L, and P23H, in addition to wildtype rhodopsin. The VPP mutation results in a progressive retinal neurodegeneration [6]. For genotyping, the following primers were used: 5'-agactgacatggggaggaattcccaga-3 (sense) and 5 -gagctgctcgaagtgactccgacc-3 (antisense). The thermal cycle protocol was denaturation at 94 • C for 30 s, annealing at 68 • C for 45 s, and elongation at 72 • C for 45 s for 35 cycles.

Microscopy and Morphometric Analyses (Spider Diagram)
The enucleated eyes were fixed for 24 h in 2.5% paraformaldehyde (PFA)/2.5% glutaraldehyde in sodium cacodylate buffer and processed as described previously [60]. Then, 1 µm thick semithin meridional sections were cut and stained according to Richardson [61]. The sections were analyzed on an Axio Imager Z1 microscope (Carl Zeiss, Jena, Germany) using Zeiss Zen software (Carl Zeiss, Jena, Germany). The thickness of the outer nuclear layer (ONL) was measured and the mean values were plotted a as spider diagram, as described previously in [13,40]. It is of interest to note that there was no sex-specific difference in the ONL thickness of VPP (male (n = 3) versus female (n = 3)) and control (male (n = 3) versus female (n = 3)) animals.

Apoptosis: TdT-Mediated dUTP-Biotin Nick End Labeling (TUNEL)
TUNEL (DeadEnd Fluorometric TUNEL, Promega, Madison, WI, USA) was used to label apoptotic cells in one-month-old animals, following our previously published protocol [10,14]. The sections were analyzed on an Axio Imager Z1 microscope (Carl Zeiss, Jena, Germany) using Zeiss Zen software (Carl Zeiss, Jena, Germany). Eyes of eleven animals (VPP: n = 11, control: n = 11) were included in the morphometric analyses and their TUNEL-positive cells were counted and normalized to the area of the ONL [mm 2 ]. There was no sex-specific difference in the number of TUNEL positive cells between VPP (male (n = 6) versus female (n = 5)) animals.
To determine Müller cell reactivity of the peripheral and central part of the retina, we measured GFAP fluorescence intensity in the area of the peripheral and central 20% of the retina, which we defined as comparable to the measurements for the spider diagrams [13,40].
We thus used the feature "Intensity Mean Value" of the ZEN software (Carl Zeiss, Jena, Germany) to determine the mean grey value of GFAP and normalized the mean value of the intensity to the total retinal area from the GCL to the outer aspect of the ONL [mm 2 ]. We found no difference in reactivity of Müller cells between the peripheral and central parts of the retinae in VPP animals (n = 3).
To study the co-localization of RNA/BaseScope and GFAP/GS, z-stacks (25-30 images, in total 5.5-7.0 µm thick) were analyzed, transformed into the ortho (orthogonal section) view, and reconstructed as a 3D image using the Zeiss Zen software and the 3Dxl rendering module with the surface function (Carl Zeiss, Jena, Germany).

RNA Isolation and Quantitative Real-Time RT-PCR (qPCR)
TriFast (Peqlab, Erlangen, Germany) was used to isolate total mRNA from retinal tissue and cDNA was synthesized using the iScript cDNA Synthesis Kit (Bio-Rad Laboratories, Inc., Hercules, CA, USA) following the manufacturer's instructions. QPCR analyses were performed on a CFX Realtime PCR Detection System (Bio-Rad Laboratories, Inc.) and as previously described [46]. All oligonucleotides (Supplementary Table S1) were designed to span exon-intron boundaries and purchased from Invitrogen (Carlsbad, CA, USA). CFX Manager TM Software and Excel were used to analyse relative mRNA expression levels according to the ∆∆C T -method [62]. The geometric mean value of the reference genes ubiquitin (Ubc) and guanine nucleotide binding protein subunit beta2 like 1 (Gnb2l1) was used for normalization. To perform RNA sequencing, total RNA of pooled retinae (right and left eye) was purified using the RNeasy Mini Kit by Qiagen (Venlo, The Netherlands).

RNA Sequencing
Library preparation and RNAseq were performed at the service facility 'KFB-Center of Excellence for Fluorescent Bioanalytics' (Regensburg, Germany. www.kfb-regensburg.de, accessed on 11 June 2021). Library preparation and RNAseq were carried out as described in the Illumina TruSeq Stranded mRNA Sample Preparation Guide, the Illumina NextSeq 500 System Guide (Illumina, Inc., San Diego, CA, USA), and the KAPA Library Quantification Kit-Illumina/ABI Prism User Guide (Kapa Biosystems, Inc., Woburn, MA, USA). In brief, 250 ng of total RNA was used for purifying the poly-A containing mRNA molecules using poly-T oligo-attached magnetic beads. Following purification, the mRNA was fragmented to an average insert size of 200-400 bases using divalent cations under elevated temperature (94 • C for 4 min). Next, the cleaved RNA fragments were reverse transcribed into first strand cDNA using reverse transcriptase and random hexamer primers. Actinomycin D was added to improve strand specificity by preventing spurious DNA-dependent synthesis. Blunt-ended second strand cDNA was synthesized using DNA Polymerase I, RNase H, and dUTP nucleotides. The incorporation of dUTP, in place of dTTP, quenched the second strand during the later PCR amplification, because the polymerase does not incorporate past this nucleotide. The resulting cDNA fragments were adenylated at the 3' ends; the indexing adapters were ligated; and, subsequently, specific cDNA libraries were created by PCR enrichment. The libraries were quantified using the KAPA SYBR FAST ABI Prism Library Quantification Kit. Equimolar amounts of each library were sequenced on a NextSeq 500 instrument controlled by the NextSeq Control Software (NCS) v2.2.0, using a 75 Cycles High Output Kit with the single index, paired-end (PE) run parameters. Image analysis and base calling were done by the Real Time Analysis Software (RTA) v2.4.11. The resulting .bcl files were converted into .fastq files with the CASAVA Software v1.8.2.

Bioinformatics
For all samples, at least 30 million reads were analyzed. Fastq files were quality controlled with FastQC v0.11.5 [63]. All files passed quality control. The reads were aligned against Ensembl Mus musculus GRCm38 version 94 using STAR aligner v2.5.3a [64]. One sample (R21753) showed poor read alignments of less than 30% and was removed from further analyses. Reads were quantified using salmon v0.8.2 [65]. All subsequent analyses were conducted in R v3.5.1. Samples were screened for outliers using PCA and clustering analysis. One sample (R21741) was identified as an outlier and removed from further analyses. Thus, the final sample number was six control and five VPP retinae. Transcriptional dysregulation was computed using tximport v1.10.0 [66] and DESeq2 v1.22.1 [67] with genotype as the variable of interest and sex of the mice as a covariate and using ashr [68] as the fold change shrinkage estimator. The Benjamini-Hochberg procedure was used to correct for multiple comparisons (p-adjusted; p adj ). For correlation network analysis, we used the normalized and variance stabilized counts from the DESeq2 analysis. Batch correction for sex was applied with limma v3.38.3 [69] keeping genotype as the variable of interest. The normalized, transformed, and batch corrected counts were used to construct a weighted gene correlation network using WGCNA v1.66 [70,71]. Heatmaps and k-mer analysis were carried out using ComplexHeatmap v2.3.2 [72]. Visualization was carried out using cytoscape v3.7.2 (http://cytoscape.org, accessed on 11 June 2021) with the Reactome FI app v7.2.1 [16]. Ontology analysis was carried out using the Enrichr website [73,74]. Scripts are available upon request.

Statistics
All results that are displayed in bar graphs are expressed as means ± SEM. Data were screened for outliers using the Grubb's outlier test in graph pad prism. Comparisons between the means of two groups were made by a two-tailed Student's t-test. p ≤ 0.05 was considered as statistically significant.

Conclusions
The parallel expression of VPP mutant and wildtype rhodopsin [6] results in a significant increase in apoptosis and thinning of the ONL to half of its thickness in retinae of three-month-old VPP animals. Intriguingly, in our transcriptome-wide analyses, we found more than 9000 dysregulated genes in retinae of VPP mice. The predominant changes in gene expression point towards induction of apoptosis, scar formation, neuroinflammation, and the upregulation of neuroprotective pathways like TGF-β, G-protein activated (e.g., endothelin), and VEGF signaling in VPP retinae. Using in situ hybridization combined with cell-type specific markers, we could show that regulatory factors such as Ccl2, Edn2, Tgfbr2, Ednrb, and Vegfr2 were also expressed in glial cells in addition to neurons. Albeit the relevance of the identified pathways needs further investigation using, e.g., different (cell-type specific) knockout mouse lines, it is tempting to speculate that modulation of neuroinflammation and neuroprotective pathways in general or, e.g., in glial cells, is a promising target for the development of new therapeutic options to delay the degeneration of photoreceptors in diseases like RP.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/ijms22126307/s1, Figure S1: Apoptosis, retinal morphology/morphometry and the expression of neuroprotective factors in VPP and control mice, Figure S2: WGCNA analysis of VPP mice, Figure  S3: Heatmaps of pathway analysis of transcriptomic changes in the VPP mice, Figure S4: Pathway analysis of transcriptomic changes in the VPP mice, Figure S5: The glial response to photoreceptor degeneration in the VPP model, Figure S6: z-stacks and corresponding 3D reconstructions, Table S1: Oligonucleotides for qPCR, Table S2: DEseq2 analysis, Table S3: WGCNA modules.   Table S2  Transforming growth factor -ß receptor type 1 Tgfbr2 Transforming growth factor -ß receptor type 2 Ubc Ubiquitin C VEGF Vascular endothelial growth factor Vegfr1 Vascular endothelial growth factor receptor type 1 Vegfr2 Vascular endothelial growth factor receptor type 2 WGCNA Weighted correlation network analysis