Retinitis Pigmentosa Results in Neurodegeneration Concomitant With Neuroinflammation, Extracellular Matrix Disorganization and the Upregulation of Neuroprotective Pathways in Glial Cells and Neurons


 BackgroundHereditary retinal degenerations like retinitis pigmentosa (RP) are amongst 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. Here, we in-depth characterized the transgenic VPP mouse model, a genetic model for autosomal dominant RP. MethodsWe 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, immunofluorescent staining and 3D reconstruction. ResultsOne month-old VPP mice showed a significantly higher number of apoptotic photoreceptor cells that resulted in a significantly thinner ONL in three months-old VPP mice, concomitant with an increase in reactivity of microglia and Müller cells. By RNASeq analysis we identified 9,256 dysregulated genes and six significantly associated gene modules in the subsequently performed WGCNA. Gene ontology enrichment showed, amongst others, dysregulation of TGF-β regulated extracellular matrix organization, factors of the (ocular) immune system/response and apoptosis. ConclusionThe predominant effect pointed towards induction of neuroinflammation and the upregulation of neuroprotective pathways like TGF-β, G-protein activated and VEGF signaling that were significantly associated with the VPP transgene-induced photoreceptor degeneration. Thus, modulation of these processes 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 that is 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 amongst the major causes of blindness in younger patients (3,4) caused by the progressive loss of rod and cone photoreceptor cells, respectively (1). Photoreceptors are the light sensitive neurons of the 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 sec, annealing at 68 °C for 45 sec and elongation at 72 °C for 45 sec for 35 cycles.

Microscopy and morphometric analyses (Spider diagram)
The enucleated eyes were xed for 24 h in 2.5% PFA/2.5% glutaraldehyde in sodium cacodylate buffer and processed as described previously (11). 1 µm thick semithin meridional sections were cut and stained after Richardson (12). The sections were analyzed on an Axio Imager Z1 microscope (Carl Zeiss, Jena, Germany) using Zeiss Zen software (Carl Zeiss, Jena, Germany). Thickness of the outer nuclear layer (ONL) was measured and the mean values were plotted as spider diagram as described previously in (13,14).

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 (15,16). The sections were analyzed on an Axio Imager Z1 microscope (Carl Zeiss, Jena, Germany) using Zeiss Zen software (Carl Zeiss, Jena, Germany). TUNEL-positive cells were counted and normalized to the area of the ONL [mm 2 ].

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 (17). All oligonucleotides (supplementary Table 1) were designed to span exon-intron boundaries and purchased from Invitrogen (Carlsbad, California, USA). CFX Manager TM Software and Excel were used to analyse relative mRNA expression levels according to the ∆∆C T -method (18). To perform RNA sequencing total RNA of the retinae was puri ed using the RNeasy Mini Kit by Qiagen (Venlo, 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). 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 Quanti cation 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 oligoattached magnetic beads. Following puri cation, the mRNA was fragmented to an average insert size of 200-400 bases using divalent cations under elevated temperature (94 °C for 4 minutes). Next, the cleaved RNA fragments were reverse transcribed into rst strand cDNA using reverse transcriptase and random hexamer primers. Actinomycin D was added to improve strand speci city by preventing spurious DNAdependent 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 ampli cation, 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 speci c cDNA libraries were created by PCR enrichment. The libraries were quanti ed using the KAPA SYBR FAST ABI Prism Library Quanti cation 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 les were converted into .fastq les with the CASAVA Software v1.8.2.

Bioinformatics
Fastq les were quality controlled with FastQC v0.11.5 (19). All les passed quality control. The reads were aligned against Ensembl Mus musculus GRCm38 version 94 using STAR aligner v2.5.3a (20). One sample (R21753) showed poor read alignments of less than 30% and was removed from further analyses. Reads were quanti ed using salmon v0.8.2 (21). All subsequent analyses were conducted in R v3.5.1. Samples were screened for outliers using PCA and clustering analysis. One sample (R21741) was identi ed as an outlier and removed from further analyses. Thus, the nal sample number was six control and ve VPP retinae. Transcriptional dysregulation was computed using tximport v1.10.0 (22) and DESeq2 v1. 22.1 (23) with genotype as the variable of interest and sex of the mice as a covariate and using ashr (24) 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 (25) 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 (26,27). Heatmaps and k-mer analysis was carried out using ComplexHeatmap v2.3.2 (28). Visualization was carried out using cytoscape v3.7.2 (http://cytoscape.org) with the Reactome FI app v7.2.1 (29). Ontology analysis was carried out using the Enrichr website (30,31). 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 signi cant.

Analysis of photoreceptor degeneration in VPP mice
To assess the degree of apoptosis in retinae of VPP mice, we visualized apoptotic cells in retinae of onemonth-old animals by TUNEL-labelling. Retinae of control animals showed only very few apoptotic cells ( Figure 1A and B). In contrast VPP mice showed a considerable number of TUNEL-positive cells in the outer nuclear layer (ONL) ( Figure 1A and B). Intriguingly, the wide majority of apoptotic cells was located in the central part of the retinae near the optic nerve head (ONH) ( Figure 1B). When quantifying the number of the apoptotic cells on retinal sections through the ONH of one-month-old mice ( Figure 1C) we found 16.28 ± 3.08 apoptotic cells per mm² ONL in retinae of control animals and 173.57 ± 17.84 cells undergoing apoptosis in VPP retinae (p < 0.001). To investigate, whether the apoptosis of photoreceptors affected retinal morphology, we analyzed semithin sections of the eyes of three-month-old control and VPP animals. We did not observe obvious morphological changes in retinae of control animals ( Figure 1D and F). However, VPP mice showed a distinct thinning of the retina that was most pronounced in the central retina ( Figure 1D and F). For quantitative analysis, we measured the thickness of the ONL layer at de ned measuring points along both hemispheres of sagittal semithin-sections though the ONH. The mean values of these individual measurements were plotted in a spider diagram ( Figure 1E) and subsequent statistical analysis con rmed a signi cant thinning of the ONL in VPP mice compared to controls. We furthermore aimed to validate the VPP model by analyzing the mRNA expression levels by qPCR of factors like leukemia inhibitory factor (Lif), broblast growth factor 2 (Fgf2) and endothelin 2 (Edn2) that are well-known to be upregulated in the context of retinal degenerations (15,32,33). Accordingly, we found them signi cantly upregulated in retinae of three-month-old VPP mice (Lif: 22.38 ± 4.13, p = 0.008; Fgf2: 13.17 ± 1.52, p = 0.005; Edn2: 46.33 ± 5.64, p < 0.001) compared to retinae of control littermates (Lif: 1.00 ± 0.18; Fgf2: 1.00 ± 0.14; Edn2: 1.00 ± 0.15) ( Figure 1G).
3.2 Transcriptional alterations in VPP retinae: RNAseq and weighted correlation network analysis (WGCNA) Subsequent, 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 annotation (Mus musculus GRCm38 v. 94) we found 30,796 genes to be expressed in the retina, of which 9,256 were dysregulated (4,636 down-and 4,620 upregulated, Figure 2A, cut off criteria: Benjamin-Hochberg adjusted p-value (p adj ) < 0.05). The top 30 dysregulated genes are shown in supplementary Table 2. Amongst 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) that was reported to be involved in eye development, were signi cantly downregulated. Genes controlling scar formation, such as brinogen-like 2 (Fgl2) as well as apoptosis e.g. caspase1 (Casp1) or Bcl2-interacting killer (apoptosis-inducing) (Bik), were upregulated. Furthermore, we found upregulation of quite a considerable number of genes associated with in ammatory or immune response functions such as C-X-C motif chemokine ligand 13 (Cxcl13), glial brillary acidic protein (Gfap), T-cell receptor T3 gamma chain (Cd3g), chemokine (C-C motif) ligand 5 (Ccl5), C-C motif chemokine ligand 2 (Ccl2) and factors, that are associated with the complement cascade like complement component factor i (C ), complement factor C4B and the Serping1 gene.
Gene ontology enrichment showed, amongst others, involvement of TGF-β regulated extracellular matrix organization, response to cytokine stimuli and apoptosis (Table 1). Photoreceptor loss was indicated by downregulation of rhodopsin signaling pathway components (Table 1). Moreover, we performed weighted gene correlation network analysis (WGCNA) to identify genotype-speci c patterns of dysregulation, upstream regulators and involved signaling pathways. WGCNA clusters co-regulated genes into modules based on their similarity of expression. Since 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 identi ed. 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 identi ed modules (clusters of co-regulated genes as shown by their colors in Figure 2B) and their correlation of expression with the genotype, are illustrated in Figure 2B. The analysis identi ed six signi cantly 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 2D and E and supplementary Figure 1A Table 2). In the Pos2 module, we found signi cant clustering of genes encoding for necroptosis, protein transport, organelle organization, cellular homeostasis, and the ribosome (supplementary Figure 1A and Table 2). The Pos3 module (supplementary Figure 1B and Table  2) showed enrichment for bone morphogenetic protein (BMP) signaling pathway components. In the Neg1 module we observed signi cant clustering of genes regulating cellular component organization and protein kinase activity (Figure 2 E and Table 2). The Neg2 module (supplementary Figure 1C and Table 2) was enriched for genes involved in nucleotide homeostasis and the Neg3 module showed signi cant enrichment for genes encoding for ribosomal proteins and protein transport (supplementary Figure 1D and Table 2).

The glial response to photoreceptor degeneration in VPP mice
In addition to the quantitative information from the RNAseq data and to further validate its results, we performed mRNA in situ hybridization and/or immuno uorescence staining to analyze the cell types expressing transcripts of interest. As we observed a signi cant upregulation of several genes regulating glial reactivity (e.g. glial brillary acidic protein (Gfap), serine protease inhibitor A3N (Serpina3n), lipocalin 2 (Lcn2), we performed immuno uorescence labelling of retinal sections for GFAP. GFAP is a major intermediate lament particularly expressed in astrocytes (34), which is upregulated in response to retinal trauma in astrocytes and Müller cells (34,35). Immuno uorescence staining against GFAP in threemonth-old control mice showed the characteristic staining pattern for astrocytes on top of the retina for both groups (Figure 3 A and 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 ), immuno uorescence staining showed an increased GFAP signal intensity, indicating an enhanced protein expression, on top of the retina and an additional stripe-like staining pattern stretching through the retina, which represents the characteristic morphological appearance of reactive Müller cells. Unlike apoptosis and ONL thinning, we could not detect a difference in GFAP reactivity between the central and the peripheral parts of the retina ( Figure 3A and B). To further validate the RNAseq data, we performed qPCR analyses on retinal samples to determine the relative Gfap expression levels, and also the expression levels of the microglia/macrophage marker ionized calcium-binding adapter molecule 1 (Iba1) and the chemokine (C-C motif) ligand 2 (Ccl2) which is reported to stimulate the migration and reactivity of microglia cells (36,37). QPCR analyses con rmed that these mRNA expression levels were signi cantly upregulated in VPP retinae (Gfap: 8.89 ± 0.63, p < 0.001; Iba1: 6.56 ± 1.03, p = 0.015; Ccl2: 68.74 ± 5.11, p < 0.001) compared to control retinae (Gfap: 1.00 ± 0.19; Iba1: 1.00 ± 0.11; Ccl2: 1.00 ± 0.11) ( Figure 3C). In accordance, our RNAseq data showed comparable increase of Iba1 (5.51-fold, p adj =1.12·10 -23 ) and Ccl2 mRNA (67.51-fold increase, p adj =8.05·10 -14 ; Figure 3F). Furthermore, we used an anti-IBA1 labeling to visualize myeloid cells e.g. microglia and recruited macrophages, in the retinae of control and VPP animals ( Figure 3D). In controls, we observed rami ed IBA1-positive cells in their typical localization on top of the retina and the inner (IPL) and outer plexiform layers (OPL). In contrast, in VPP retinae, IBA1positive cells changed their shape from rami ed 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) ( Figure 3D). Taken together, these results showed a pronounced reactivity of macro-and microglia cells in response to photoreceptor degeneration. To further study the origin of the signi cantly elevated Ccl2 expression in VPP retinae and to supplement the quantitative information from the RNAseq ( Figure 3F) and qPCR data ( Figure 3C), we performed Ccl2 mRNA in situ hybridization on retinal sections combined with immuno uorescence co-labeling of glial cells. Immuno uorescence staining for glutamine synthetase (GS) staining was used to label Müller cells (38) and for GFAP to detect astrocytes and reactive GFAP-positive Müller cells. In control retinae, we observed a rather low number of Ccl2 punctae in the inner nuclear layer (INL), the ONL and a few signals in the retinal ganglion cell layer (GCL) ( Figure 3E and supplementary Figure 4A and B). The number of Ccl2 punctae was markedly increased in the INL and ONL of the VPP retinae ( Figure 3E). When using Ccl2/GFAP/GS co-labeling we observed sparse overlap of Ccl2 in GFAP-positive astrocytes ( Figure 3G and supplementary Figure 4A) and more frequent overlap in GS-positive resting and in GFAP-/GS-positive reactive Müller cells ( Figure 3H and supplementary Figure 4B). However, we also detected Ccl2 mRNA expression in cells other than Müller glia and astrocytes in the neuronal layers of the retina (GCL, INL and ONL) pointing towards additional expression in retinal neurons ( Figure 3E).

.
Thus, we aimed to investigate the impact of these signaling pathways (TGF-β-, G-protein activated-and VEGF-signaling) on the VPP model in detail. Unsupervised hierarchical clustering of the samples that we generated on basis of the Reactome pathway database (29) for the VEGF ( Figure 4A), TGF-β ( Figure 4B) and G-protein mediated signaling pathways ( Figure 4C) 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 2. To further analyze sub-groups of dysregulated pathways, we transformed the Reactome pathways into functional interaction networks (Figure 4 D, E and F). We colored the genes of each network according to their dysregulation: white indicates no signi cant regulation, red genes were up-and blue genes were signi cantly 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 3. This analysis identi ed 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 4E and supplementary Figure 3C) or Tgfbr2 in the TGF-β family signaling network ( Figure 4D and supplementary Figure 3B) and Vegfr2/kinase insert domain receptor (Kdr) in the VEGF signaling network ( Figure 4F and supplementary Figure 3A).
Next, we used RNA/BaseScope® in situ hybridization to supplement the quantitative information from the RNAseq/qPCR data and to identify speci c cell types expressing transcripts of interest. Using a speci c 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 5A). VPP retinae showed signi cantly higher Edn2 expression in the RNAseq data (22.12-fold, p adj =8.68·10 -87 ; Figure 5B) (supplementary Table 2)/qPCR analyses ( Figure 1G) 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 5A). When performing GFAP/GS co-labelling, we observed some Edn2 signals in  Figure 4D). However, we observed the majority of the Edn2 signals in the neurons of the INL and ONL ( Figure 5A). Ednrb mRNA in situ hybridization showed speci c signals in the INL and ONL and some de ned signals in the GCL in control retinae ( Figure 5E). We furthermore observed some Ednrb expression in GFAP-positive astrocytes ( Figure 5G and supplementary Figure 4E) and in GS-positive resting Müller cells ( Figure 5H and supplementary Figure 4F). In VPP retinae, Ednrb was signi cantly upregulated (1,29-fold, p adj =0.007; Figure 5F) in the RNAseq data (supplementary Table 2). In the in situ hybridization, we detected pronounced Ednrb signals in the INL and ONL that overlapped to some extent with GFAP-/GS-positive reactive Müller cells ( Figure 5E and H and supplementary Figure 4F).
Tgfbr2 was signi cantly upregulated (2.23-fold, p adj =2.18·10 -23 ; Figure 6B) in VPP retinae in the RNAseq data (supplementary Table 2). 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 6A and supplementary Figure 4G and H). In VPP retinae, the number of Tgfbr2 punctae increased in the INL and ONL. Immuno uorescence co-labelling con rmed its expression in only some isolated GFAP-positive astrocytes ( Figure 6C and supplementary Figure 4G) and its association with resting, GS-positive and reactive GFAP-/GS-positive Müller cells, respectively ( Figure 6D and supplementary Figure 4H). Yet, we also observed Tgfbr2in 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 its additional expression in neuronal cells ( Figure 6A).
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 6E and supplementary Figure 4J). Moreover, we detected Vegfr2 mRNA signals in the ONL and isolated signals in the GCL that overlapped to some extent with GFAP-positive astrocytes ( Figure 6G and and supplementary Figure 4I). Our RNAseq analysis showed Vegfr2 to be signi cantly upregulated in VPP retinae (2.00-fold, p adj =1.40·10 -41 ; Figure 6F; supplementary Table 2) and accordingly Vegfr2 in situ hybridization showed an increased in expression in the INL and ONL ( Figure 6E). Co-labelling showed its association to and expression in GFAP-/GSpositive reactive Müller cells ( Figure 6H and supplementary Figure 4J). 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 its additional expression in neuronal cells ( Figure 6E).

Discussion
The present data con rm that the VPP model displays the major phenotypic characteristics of the human disease retinitis pigmentosa. Brie y, we demonstrate in comprehensive transcriptome-wide analyses of retinae from three-month-old VPP mice (1) an extensive dysregulation of genes encoding for 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, resting and reactive glia cells, and (4) a dysregulation of extracellular matrix organization, apoptosis and response to cytokine stimuli in WGCNA analyses.
The transcriptional response to photoreceptor degeneration leads to increased expression of genes regulating in ammatory or immune response functions Neuroin ammation is a common hallmark of the pathogenesis of neurodegenerative diseases like Alzheimer's, Parkinson's, multiple sclerosis or retinal degenerations (40,41). Following a neurotoxic event, neuronal stress signals mediate reactivity of microglial cells leading to their proliferation, migration and the secretion of speci c cytokines and chemokines that can exert neurotoxic or neuroprotective effects (41,42). Sustained reactivity of microglia promotes chronic in ammation and may cause irreversible neuronal cell death (41,43,44). Thus, the accumulation of reactive IBA1-positive cells in the OPL and in the sub-neuroretinal space in VPP retinae strongly indicates an ongoing neuroin ammatory process.
Moreover, in the top 30 dysregulated genes, we found a considerable number of genes that are associated with in ammatory 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 neuroin ammation. These ndings are in accordance with previously published data, which describe upregulation of factors like Lif, Ccl2 (Mcp-1), interleukin-1 (Il-1),complement component 1q (C1q) and complement factor H (CFH) in retinae of genetic mouse models of RP (32,(45)(46)(47)(48)(49). As we and others (36,37) have shown Ccl2 is expressed in Müller cells and photoreceptors in the healthy retina and upon retinal damage contributing to the recruitment of microglia/in ltrative macrophages (37,50). However, con icting 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 null background for Ccl2 (VPP;Ccl2 -/-) did not differ to that of transgene VPP mice (32). In contrast, Rutar and colleagues demonstrated that siRNA-mediated knock down of Ccl2 resulted in a signi cantly lower number of apoptotic photoreceptors in rats after lightinduced photoreceptor degeneration and (37,50). Based on our data, it is very likely that in VPP retinae the elevated Ccl2 expression in Müller cells and in the ONL contribute 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 that 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 (47)(48)(49)(51)(52)(53)(54). Activation of the complement system promotes microglia/in ltrative macrophages migration and eventually complement activated lysis (48,51). Still, con icting data exist regarding the exact role of complement system activation and its impact on photoreceptor degeneration. Mice with a de ciency in complement factor D are protected from light-induced photoreceptor degeneration (47), indicating a detrimental role for photoreceptor survival, but the de ciency in complement component 3 (C3) or complement receptor 3 (CR3) in a genetic mouse model of photoreceptor degeneration increases microglia-mediated neurotoxicity to photoreceptors (54). Thus, the detailed function of the complement system and its speci c 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 signi cant upregulation of genes that are associated with in ammatory or immune response functions leads to neuroin ammation 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 (55). Recently, our group and others showed that in response to retinal injury, endothelin 2 (Edn2) is expressed by photoreceptors concomitant with an 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 (14,15,32,33,39). Our RNAseq data, in situ hybridizations, immuno uorescence staining and qPCR analyses con rmed 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 were signi cantly upregulated following induced ocular traumata (39). 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 (45).
The transcriptional response to photoreceptor degeneration leads to upregulation of pro-apoptotic factors and extracellular matrix organization Our transcriptomic analyses also showed an upregulation of factors that are associated with apoptosis (Casp1, Bik) and scar formation such as brinogen like 2 (Fgl2) and Tgf-β1. A well-described characteristic of TGF-β signaling is its contribution to wound healing, tissue brosis and scar formation (56, 57). Accordingly, the TGF-β regulated organization of extracellular matrix was a major hit in our gene ontology analyses and might well be the result of a healing response following photoreceptor degeneration in VPP retinae. We furthermore identi ed an upregulation of 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 (58-60). Moreover, we recently showed that the deletion of TGF-β signaling results in the development of retinal microaneurysms and choroidal neovascularization (13,17), clearly emphasizing its potential to regulate angiogenic processes. We furthermore showed that TGF-β signaling protects retinal neurons from programmed cell death during retinal development (16), thus highlighting its neuroprotective properties (61).
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 (35,62) and eventually results in a glial scar. The identi ed dysregulation of genes involved in neuroin ammation, 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.

Conclusion
The parallel expression of VPP mutant and wildtype rhodopsin (6) results in a signi cant 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 9,000 dysregulated genes in retinae of three-month-old VPP mice. The predominant changes in gene expression point towards induction of apoptosis, scar formation, neuroin ammation, and the upregulation of neuroprotective pathways like TGF-β, G-protein activated (e.g. endothelin) and VEGF signaling in VPP retinae. Using in situ hybridizations combined with cell-type speci c 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. Thus, the modulation of these processes in general or e.g. in glial cells are promising targets for the development of new therapeutic options to delay the degeneration of photoreceptors in diseases like RP.

Consent for publication
Not applicable Availability of supporting data The raw data les of the RNAseq data (suppl. table 2 and 3) are available from the authors upon request.

Competing interests
The authors declare that they have no competing interests.           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 signi cantly dysregulated; red -signi cantly upregulated and blue -signi cantly 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 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 signi cantly dysregulated; red -signi cantly upregulated and blue -signi cantly 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 Figure 5 Upregulation of endothelin signaling in VPP mice A. In situ hybridization for Edn2 (red) and GFAP (green) /GS (purple) immuno uorescence co-labeling in the retinae of three-month-old animals. Nuclei were DAPIstained (blue). Edn2 (red, arrowheads) signals were detectable in the ONL and INL with some sparse signals in the GCL. In the VPP retina, the number of the Edn2 signals (red, arrowheads) was increased in the ONL and INL and the Müller cells were GFAP (green) /GS (purple) -positive, indicating their reactivity.