3.1 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 one-month-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 defined 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 confirmed a significant 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), fibroblast 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 significantly 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 (padj) < 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 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 upregulated. 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), C-C motif chemokine ligand 2 (Ccl2) and factors, that are associated with the complement cascade like complement component factor i (Cfi), 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-specific 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 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 2B) and their correlation of expression with the genotype, are illustrated in Figure 2B. 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 2D and E and supplementary Figure 1A - D). The Pos1 module was significantly enriched for genes involved in cellular transport and differentiation (Figure 2 D and Table 2). In the Pos2 module, we found significant 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 significant 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 significant enrichment for genes encoding for ribosomal proteins and protein transport (supplementary Figure 1D and Table 2).
3.3 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 immunofluorescence staining to analyze the cell types expressing transcripts of interest. As we observed a significant upregulation of several genes regulating glial reactivity (e.g. glial fibrillary acidic protein (Gfap), serine protease inhibitor A3N (Serpina3n), lipocalin 2 (Lcn2), we performed immunofluorescence labelling of retinal sections for GFAP. GFAP is a major intermediate filament particularly expressed in astrocytes (34), which is upregulated in response to retinal trauma in astrocytes and Müller cells (34, 35). Immunofluorescence staining against GFAP in three-month-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, padj=3.87·10-78), immunofluorescence 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 confirmed that these mRNA expression levels were significantly 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, padj=1.12·10-23) and Ccl2 mRNA (67.51-fold increase, padj=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 ramified 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, 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) (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 significantly 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 immunofluorescence co-labeling of glial cells. Immunofluorescence 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).
3.4 Dysregulation of potentially neuroprotective pathways in VPP retinae: VEGF-, TGF-β- and G-protein mediated signaling
As a follow up to our previously published studies (13, 15, 16) on the neuroprotective properties of signaling pathways such as the transforming growth factor (TGF) –β signaling, G-protein activated signaling and vascular endothelial growth factor (VEGF) mediated signaling, we investigated their potential regulation in the VPP model. Quite intriguingly, our RNAseq data analysis (supplementary Table 2) showed a significant upregulation of genes encoding for members of the G-protein activated signaling family. Here we particularly focused on the endothelin family, as our group and others recently showed that endothelin 2 (Edn2) and endothelin receptor b (Ednrb) are upregulated following photoreceptor damage (15, 32, 33, 39). In accordance the RNAseq data (supplementary Table 2) of the VPP retinae showed a significant increase in Ednrb (1.28-fold, padj = 0.0074) and Edn2 (22.12-fold, padj = 8.68·10-87) expression. We furthermore observed an upregulation of factors involved in TGF–β signaling (e.g. Tgf-β receptor type 1 (Tgfbr1): 1.15-fold, padj = 0.013; Tgfbr2: 2.23-fold, padj = 2.18·10-23; Tgf-β1: 2.24, padj = 1.24·10-12; Tgf-β2: 1.51-fold, padj = 2.20·10-20). Of note, Tgf-β2 was additionally identified as one of the hub gene in the Pos1 module. Vascular endothelial growth factor (VEGF) receptor signaling pathway was identified in the gene ontology enrichment analysis of the significantly upregulated genes and the RNAseq data (supplementary Table 2) showed an increased expression of Vegfr1 (Flt1): 1.25, padj = 4.75·10-7, Vegfr2 (Kdr): 2.14, padj = 1.40·10-41; Vegfb: 1.20, padj = 0.00013; and Vegfc: 1.46, padj = 9.31·10-5.
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 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 3. 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 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 specific 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 5A). VPP retinae showed significantly higher Edn2 expression in the RNAseq data (22.12-fold, padj=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 GFAP-positive astrocytes (Figure 5C and supplementary Figure 4C) and in GS-positive resting and GFAP/GS-positive reactive Müller cells, respectively (Figure 5 D and supplementary 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 specific signals in the INL and ONL and some defined 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 significantly upregulated (1,29-fold, padj=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 significantly upregulated (2.23-fold, padj=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. Immunofluorescence co-labelling confirmed 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 significantly upregulated in VPP retinae (2.00-fold, padj=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-/GS-positive 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).