Exploring Early-Stage Retinal Neurodegeneration in Murine Pigmentary Glaucoma: Insights From Gene Networks and miRNA Regulation Analyses

Abstract

Purpose: Glaucoma is a group of heterogeneous optic neuropathies characterized by the progressive degeneration of retinal ganglion cells. However, the underlying mechanisms have not been understood completely. We aimed to elucidate the genetic network associated with the development of pigmentary glaucoma with DBA/2J (D2) mouse model of glaucoma and corresponding genetic control D2-Gpnmb (D2G) mice carrying the wild type (WT) Gpnmb allele.

Methods: Retinas isolated from 13 D2 and 12 D2G mice were subdivided into 2 age groups: pre-onset (1–6 months: samples were collected at approximately 1–2, 2–4, and 5–6 months) and post-onset (7–15 months: samples were collected at approximately 7–9, 10–12, and 13–15 months) glaucoma were compared. Differential gene expression (DEG) analysis and gene-set enrichment analyses were performed. To identify micro-RNAs (miRNAs) that target Gpnmb, miRNA expression levels were correlated with time point matched mRNA expression levels. A weighted gene co-expression network analysis (WGCNA) was performed using the reference BXD mouse population. Quantitative real-time PCR (qRT-PCR) was used to validate Gpnmb and miRNA expression levels.

Results: A total of 314 and 86 DEGs were identified in the pre-onset and post-onset glaucoma groups, respectively. DEGs in the pre-onset glaucoma group were associated with the crystallin gene family, whereas those in the post-onset group were related to innate immune system response. Of 1329 miRNAs predicted to target Gpnmb, 3 miRNAs (miR-125a-3p, miR-3076-5p, and miR-214-5p) were selected. A total of 47 genes demonstrated overlapping with the identified DEGs between D2 and D2G, segregated into their time-relevant stages. Gpnmb was significantly downregulated, whereas 2 out of 3 miRNAs were significantly upregulated (P < 0.05) in D2 mice at both 3-and 10-month time points.

Conclusions: These findings suggest distinct gene-sets involved in pre-and post-glaucoma in the D2 mouse. We identified three miRNAs regulating Gpnmb in the development of murine pigmentary glaucoma.

Glaucoma describes a group of heterogeneous optic neuropathies characterized by the progressive degeneration of retinal ganglion cells (RGCs) and accompanying pathologic outcomes.^1,2 It is considered to be the leading cause of irreversible blindness and a major cause of visual impairment worldwide.^3–6 Despite this status, molecular mediators of spontaneous intraocular pressure (IOP) elevation and progression to blindness remain ill-defined. Glaucoma development is a complex process, and no single factor is sufficient to induce its pathological phenotypes.^7,8 Indeed, even with the presence of well-established risk factors, such as IOP elevation, there is a marked variability in the rate of progression and severity.⁷ A lack of understanding of the genetic and environmental factors to the development of glaucoma has hampered screening efforts.^9,10 This is even more alarming considering that early identification is critical, as glaucoma is largely silent until irreversible visual damage has occurred, and early intervention has been shown to decrease the progression rates.^2,11,12 Furthermore, treatments are limited in controlling IOP elevation, leaving few alternative therapies.^13,14 Thus, continued study of the underlying mechanisms of glaucomatous damage is important for developing effective new screening, prevention, and treatment. 

Rodent models remain one of the most powerful tools for studying the molecular mechanisms of glaucoma.^15,16 Among the various models used, the DBA/2J (D2) mouse is considered a congenital experimental model of glaucoma as it develops a naturally occurring chronic secondary angle-closure glaucoma.^17–20 D2 mice develop a form of glaucoma that results from the abnormal liberation of iris pigment (iris pigment dispersion [IPD]) in the anterior chamber, which obstructs drainage routes for aqueous humor and results in a gradual elevation of IOP. The IPD occurs in D2 mice due to a mutation in the glycoprotein (transmembrane) nmb gene (Gpnmb^R150X). Iris stromal atrophy, a separate phenotype, occurs due to a mutation in the tyrosine protein type 1 gene (Tyrp1^b).^{18,19,21–24} Together, these result in “pigment dispersion syndrome” (PDS) and then pigmentary glaucoma (PG).^22,23,25 The progression of PDS and PG follows a predictable timeline with iris destruction beginning at approximately 4 months of age, followed by IOP elevations due to liberated iris pigment obstruction of the anterior chamber drainage structures at 7 to 9 months, and progressive RGC death after 9 months.^19,26–28 

Despite this predictable pattern, there is still considerable variability in IOP elevation and severity of RGC damage in D2 mice.^19,29 This variation is thought to be due, at least in part, to the environment and other genetic factors that interact with the Tyrp1^b and Gpnmb^R150X mutations. In a previous study, novel double-congenic C57BL/6J (B6) mice were bred to express D2-derived Tyrp1^b and Gpnmb^R150X (B6.D2-Tyrp1^bGpnmb^R150X mice). Although these B6 double-congenic mice showed similar patterns of pigment dispersion to D2 mice, the B6 strains did not develop glaucomatous damage and had a diminished IOP elevation.³⁰ In this study, we compared the retinas from D2 mice (with the Tyrp^b and Gpnmb^R150X mutation) with D2-Gpnmb (D2G) mice that are genetically identical to D2 mice, but contain wild-type Gpnmb. Phenotypically, D2G mice only develop mild IPD due to the presence of Tyrp^b mutation, but neither PDS nor PG at any age.²⁴ Understanding the network in which the Gpnmb^R150X mutation contributes to the development of PDS and PG, is imperative for understanding the pathogenesis of glaucoma. 

Development of glaucoma phenotype in D2 mice includes multiple contributions from genes and factors that are distinct from the two mutations that regulate aqueous outflow in the anterior chamber. A class of small noncoding RNAs (approximately 19–25 nucleotides in length), known as micro-RNAs (miRNAs), can specifically recognize and bind to the 3′ untranslated regions (UTRs) of target genes, and as a result inhibit translation.^31,32 Through this mechanism, a single miRNA can act as a post-transcriptional regulator of multiple mRNA targets and influence cellular functions.³³ We hypothesize that understanding the miRNA expression patterns involved in PDS and PG will help facilitate outlining the orchestrated effects of genes and proteins associated with glaucoma pathogenesis. Retinal expression of miRNAs has been demonstrated to contribute to retinal development and function.^34–37 Additionally, the role of miRNAs in retinal function has been reported in several models of glaucomatous damage using Brown Norway rats,³⁸ Sprague Dawley rats,³⁹ E50K mutant mice,^40,41 and human tissues.⁴² However, investigations concerning miRNA expression in D2 mouse retina have been relatively limited. 

In the current study, we aimed to better characterize the network associated with Gpnmb^R150X that results in D2 glaucoma phenotype. We examined differentially expressed genes (DEGs) and miRNAs between D2 and D2G mice during pre-onset (1–6 months of age) and post-onset (7–15 months of age) glaucoma stages. We identified DEGs comprising the crystallin family and innate immune system family in the pre-onset glaucoma and post-onset glaucoma groups, respectively. Additionally, we used a systems genetics approach utilizing BXD recombinant inbred (RI) mice to construct a proposed miRNA-mediated gene network that may influence the development of the D2 glaucoma phenotype. A genetic reference population is a powerful tool for systems genetics analyses. The current BXD family is the largest mouse genetic reference population panel available and is derived from a cross between B6 and D2 mouse parental strains. The Tyrp1^b and Gpnmb^R150X mutations expressed by D2 mice segregate among BXD progeny, producing a gradient of IPD and PG phenotypes.⁴³ Furthermore, BXD lines have been extensively used in genetic and genomic studies of the eye.^43–49 Here, we use BXD lines to identify potential retinal miRNA regulators of glaucomatous development. We identified several potential miRNA regulators important to pigmentary glaucoma development in the retina of D2 mice. Ultimately, we expect that subsequent investigation of these identified gene families and miRNAs will yield novel insights into the complex underpinnings of glaucoma development and progression. 

Materials and Methods

Animals

The D2 (n = 13) and D2G (n = 12) mice were obtained from the Jackson Laboratory and bred at the University of Tennessee Health Science Center (UTHSC) animal facility. The animals used in this study were housed on a 12:12 light/dark cycle with ad libitum access to water and food. Mice were euthanized at approximately 1 to 2, 2 to 4, 5 to 6, 7 to 9, 10 to 12, and 13 to 15 months of age. At each time point, 2 mice per strain (one male and one female), except for 4 D2 mice at approximately 2 to 4 months and one at approximately 5 to 6 months, were euthanized by cervical dislocation following isoflurane inhalation, after which the retinas were isolated. All animal experiments were conducted in accordance with the procedures approved by the Institutional Animal Care and Use Committee at The University of Tennessee Health Science Center. 

mRNA and miRNA Gene Expression Profiling

Total RNA was isolated and purified from both the retinal tissues (pooled) of each animal using the miRNeasy Mini Kit (Qiagen, Hilden, Germany) following the manufacturer's instructions. Briefly, the tissues were first placed into 2 mL tube with 700 µL TRIzol reagent and homogenized with tissueLyzer II (Retsch, Haan, Germany) for 3 minutes at 200 × r, then transferred into a new tube with 150 µL chloroform followed by centrifugation at 13,000 × r for 15 minutes at −4°C. After phase separation by chloroform, 1.5 × volume of 100% ethanol was added to the aqueous phase. This mixture was loaded into miRNeasy column followed by purification with the RWT and RPE buffers according to the manufacturer's instructions. Total RNA was eluted by adding 50 µL of RNase-free water and stored in a freezer at −80°C. The purity and quantity of RNA was determined by NanoDrop One (Nanodrop, Wilmington, DE, USA), and its integrity was determined with Bioanalyzer 2100 (Agilent, Santa Clara, CA). Extracted RNA samples with OD260/230 greater than 1.80 and RNA integrity number greater than 8 were reversely transcribed into cDNA, and were subsequently hybridized onto the Affymetrix GeneChip Mouse Transcriptome Array 1.0 (Thermo Fisher Scientific, Waltham, MA, USA) for mRNA and Affymetrix Mouse GeneChip miRNA Arrays 4.0 (Thermo Fisher Scientific) for miRNA profiling according to the manufacturer's protocol. 

Data Normalization and Preprocessing

The raw microarray data corresponding to miRNA and mRNA expression were normalized independently using the robust multi-array average⁵⁰ (RMA) method through Affymetrix Expression Console Software. The modified Z-scores method was subsequently applied to re-normalize both the data sets.⁵¹ We calculated the log-base 2 of the normalized values, and the data of each array were converted to their Z-scores followed by multiplying by 2 and adding 8 units to each value. This transformation rescaled all the values to positive and all the data sets to a mean of 8 units with a standard deviation of 2, which resulted in a 2-fold difference in expression corresponding approximately to a one-unit change. In this study, probes with average expression less than 7 were considered to be background noise and removed from the analysis. The normalized mRNA and miRNA expression data can be accessed through our GeneNetwork platform with accession numbers GN803 and GN806, respectively. To access the GN803 (mRNA) dataset on GeneNetwork, the species (mouse [mm 10]), group (glaucoma and aged retina UTHSC), type (retina mRNA), and data set (UTHSC B6D2 Retina Affy MoGene 1.0ST [Sep 16] Gene Level RMA) can be selected. To access GN806 (miRNA), the species (mouse [mm 10]), group (glaucoma and aged retina UTHSC), type (retina mRNA), and data set (UTHSC B6D2 Retina Affy miRNA-4.0 [Nov 16] RMA) can be selected. 

Differentially Expressed Gene and Correlation Analyses

The Limma package from R Bioconductor⁵² was used for determining the statistically significant differences in the expression levels of genes between D2 and D2G for pre-onset (1–6 months of age) and post-onset (7–15 months of age) glaucoma stages, in which sex was included as a covariate. The criteria used for defining DEGs included a P value < 0.05 and fold change > 1.3 between the 2 strains. The gene abundance relationships between Gpnmb and miRNAs were analyzed with Pearson product correlations, and those with a P value < 0.05 were considered statistically significant. 

Gene Enrichment Analysis

The enrichment analysis of DEGs for the identification of overrepresented GO and MPO terms was performed using WebGestalt (http://www.webgestalt.org/) with default parameters.⁵³ Mouse genome was used as reference gene set and a minimum number of genes per category was set to five for both analyses. The Benjamini and Hochberg correction was used to account for multiple testing, and a False Discovery Rate (FDR) < 0.05 indicated significant over-representation. 

Gene Expression Pattern Detection

We used the approach proposed by Dipen P. Sangurdekar, which is an extension of previously published methods,^54,55 to infer the differentially expressed genes and categorize their dynamic expression pattern. This approach has been implemented in “Rnits” statistical package in R environment. Briefly, average expression values per probe per time point for D2 and D2G were subjected to Rnits to build Rnits objects according to the protocol. Then the model was fitted using gene-level summarization and by clustering all genes to identify the differentially expressed genes along with the summary statistics. 

miRNA Target Prediction

The miRNA targets were predicted using miRWalk 2 (http://zmf.umm.uni-heidelberg.de/apps/zmf/mirwalk2/).⁵⁶ We considered the miRNA-target pairs that were supported by one of the following prediction algorithms: miRWalk, miRanda, RNA22, and Targetscan. 

Weighted Gene Co-Expression Network Analysis

We used our previously generated, “Full_HEI_Retina_April_2010_RankInv” retina expression data⁵⁷ to construct gene co-expression networks by using the weighted gene co-expression network analysis (WGCNA) version 1.63 package in R environment (version 3.1.3).⁵⁸ This data set includes expression profiles for 77 BXD RI strains, B6 and D2 parental strains, and reciprocal F1 hybrids, and is publicly accessible through our GeneNetwork website (http://www.genenetwork.org) with the accession number GN267. First, a hierarchical clustering using a set of dissimilarities for the objects being clustered was performed using the function hclust, implemented in WGCNA, to cluster the samples and detect outliers. We considered probes with an average expression >7 to construct the gene co-expression network using WGCNA, as previously described.⁵⁸ Briefly, pair-wise Pearson correlation coefficients were calculated between all probes, then a soft thresholding power (β) of 6 was selected based on approximate scale-free topology (R² > 0.9) to generate a signed weighted adjacency matrix. Further, the adjacency matrix was transformed into Topological Overlap Matrix (TOM), which assesses the transcript interconnectedness. Following this, a dissimilarity measure was calculated. Genes were aggregated into modules by hierarchical clustering based on TOM and further refined using the dynamic tree cut algorithm. 

Quantitative Reverse Transcription PCR

The differences in the expression levels of Gpnmb, miR-125, miR-214, and miR-3076 between D2 and D2G mice were evaluated using quantitative real-time PCR (qRT-PCR). Total RNA was extracted from retinas of D2 and D2G strains at 2 time points (3-4 and 9-11 months, and 3-4 mice per strain per time point). The concentration and integrity of the total RNA was determined with NanoDrop One (Thermo Scientific, Wilmington, DE, USA), and Bioanalyzer 2100 (Agilent, Santa Clara, CA, USA), respectively. Total RNA from each retina was transcribed into cDNA with a reverse transcription Kit (Invitrogen, Carlsbad, CA, USA) following the manufacturer's instructions. Then, the cDNA was used as template to amplify the specific products for individual genes using the SYBR Green-based real-time PCR on a LightCycler 4800 real-time PCR instrument (Roche Applied Science; Indianapolis, IN, USA). The relative expression of each gene was normalized to β-actin using the 2^−ΔΔCt method, and data were presented as mean ± SD based on the average of expression levels calculated. The miRNA expression was detected using the method described previously.⁵⁹ 

Results

Differentially Expressed Genes Between D2 and D2G Mice

We used the Mouse Transcriptome Array 1.0 to quantify the expression levels of retinal genes between D2 and D2G mice during pre- (1–6 months of age) and post-onset (7–15 months of age) glaucoma stages. Retinal samples were collected from mice at approximately 1 to 2, 2 to 4, 5 to 6, 7 to 9, 10 to 12, and 13 to 15 months of age. A total of 397 DEGs were identified (FC > 1.3 and P < 0.05; mean expression >7) between D2 and D2G mice retina. Among those, 314 genes were from the pre-onset, whereas 86 genes were from the post-onset glaucoma group, respectively. Only 3 genes, Sparc, S1pr3, and Gpnmb, were found to be differentially expressed in both pre- and post-onset glaucoma groups. The gene ontology (GO) enrichment analysis revealed that pre-onset glaucoma DEGs are mainly involved in eye structure and function related biological processes (Table 1), with expression levels higher in D2 than in D2G (Table 2), whereas post-onset glaucoma DEGs are significantly enriched in immune system related biological processes (Table 1) with higher expression in D2 than in D2G, except for Gpnmb (Table 3). Furthermore, we looked up the expression levels of Caspase genes, the indicators of degeneration, in our data set and found that Card14 and Casp12 showed differential expressions between D2 and D2G in the age of 13 to 15 months and 10 to 12 months, respectively. In addition, Casp8ap2 showed differential expression in pre-onset group and Casp12 and Casp8 showed differential expression in post-onset group (data not shown). 

Table 1.

 

Top 5 Enriched Biological Processes Represented by DEGs Between D2 and D2G Mice in Each Glaucoma Group

 

Table 2.

 

Eye Structure and Function Related DEGs

 

Table 3.

 

Immune System Process Related DEGs

 

We also performed mammalian phenotype ontology (MPO) enrichment analysis using WebGestalt for the pre- and post-onset glaucoma DEGs, respectively. Our results were in agreement with those obtained from GO analysis. The DEGs in pre-onset glaucoma were significantly involved in a number of annotations, including total cataracts (FDR = 0) and abnormal lens fiber morphology (FDR = 0; Table 4), whereas post-onset glaucoma DEGs were mainly involved in abnormal inflammatory response (FDR = 3.24E-06; see Table 4). 

Table 4.

 

Top 5 Enriched MPOs Represented by DEGs Between D2 and D2G Mice in Each Glaucoma Group

 

In addition, we also detected the gene expression patterns between D2 and D2G with the time-course expression data using Rnits, and identified approximately 1000 differentially expressed genes (P < 0.05). These genes were mainly categorized into 3 expression patterns: (1) highly differentially expressed in pre-onset glaucoma but not in post-onset glaucoma, for example, Fam98b (Fig. 1A); (2) highly differentially expressed in post-onset glaucoma but not in pre-onset glaucoma, for example, Serpine1 (Fig. 1B); and (3) highly differentially expressed in both pre- and post-onset glaucoma, including Brp44l (Fig. 1C). 

Figure 1.

 

Gene expression patterns of three DEGs. (A) Fam98b was differentially expressed in pre-onset glaucoma but not in post-onset glaucoma. (B) Serpine1 was differentially expressed in post-onset glaucoma but not in pre-onset glaucoma. (C) Brp44l was differentially expressed in both pre- and post-onset glaucoma.

miRNAs Targeting Gpnmb

To identify miRNAs that potentially target Gpnmb, we measured miRNA expression levels with time-point matched mRNA expression levels. A total of 1329 miRNAs were predicted to target Gpnmb among the investigated miRNAs (> 2500) spotted on mouse miRNA arrays 4.0. Of the predicted miRNAs, 31 had a significant negative correlation (P < 0.05) with Gpnmb expression, 16 of which either had significant differential expression (P < 0.05) between D2 and D2G or had a fold change > 1.5 in at least one time point. Finally, 3 miRNAs, including miR-3076-5p, and the well-characterized miR-125a-3p and miR-214-5p were predicted to target Gpnmb by at least 2 different prediction methods, and had a P value < 0.05, and fold-change >1.5 (Table 5, Fig. 2). 

Table 5.

 

Summary of Top Three miRNAs Targeting Gpnmb

 

Figure 2.

 

Top three miRNAs significantly negatively correlated with Gpnmb expression. (A) Gpnmb expression levels for D2 and D2G along the time course. (B-D) Correlations between Gpnmb and miR-214-5p, miR-3076-5p, and miR-125a-3p, respectively. The expression values are log2 normalized.

miRNA Mediated Gpnmb-Related Gene Network

To identify Gpnmb-related gene networks, we constructed a gene co-expression network with WGCNA using our previously generated retina expression data in the reference BXD mouse population. One outlier strain, BXD24, was removed based on the sample clustering result. This strain has a mutation in CEP290 that results in a robust photoreceptor degeneration, which skews any data set that includes these mice. Finally, a total of 22922 probes with mean expression level > 7 were used for the network analysis. Soft thresholding power (β) of 6 was used to construct the gene co-expression network to achieve scale free topology (Figs. 3A, 3B). We identified a total of 34 modules (Fig. 3C); merging of these based on the clustering of eigengenes (first principal component of each module) and a height cut off of 0.25 resulted in 23 modules (Figs. 3C, 3D). The module size varied significantly from 73 genes in the “dark olive-green” module to 3887 genes in the “turquoise” module. Gpnmb was part of the “green” module, which contained a total of 705 genes. 

Figure 3.

 

Gene coexpression network analysis of BXD retina expression data using WGCNA. (A) Scale-free fit index (y axis) as a function of the soft-thresholding power (x axis). (B) Mean connectivity (degree, y axis) as a function of the soft-thresholding power (x axis). (C) Hierarchical clustering of modules based on correlations between their eigengenes. (D) Cluster dendrogram and module assignment for mRNA modules.

Among the “green” module genes, 49 genes, including Gpnmb, overlapped with the identified DEGs between D2 and D2G. These genes were identified as either pre- or post-onset glaucoma-related genes. The miRNA-mediated Gpnmb involved gene network, which was constructed using VisANT software⁶⁰ is shown in Figure 4. 

Figure 4.

 

miRNA-mediated Gpnmb gene network. Blue (left) and orange (right) nodes represent pre-onset glaucoma and post-onset glaucoma DEGs, respectively. The network shows targeting of Gpnmb (red node) by three miRNAs (green nodes). Gpnmb further interacts with pre- and post-onset glaucoma genes. The network corresponding to the post-onset glaucoma genes is relatively denser than that of the pre-onset glaucoma genes. A large number of interactions can also be observed between pre- and post-onset glaucoma genes. Each node represents a gene/miRNA, whereas an edge (line) between them represents interaction. The network was constructed using VisANT software.⁶⁰

Validation of Gpnmb and the Targeting miRNA Levels by Quantitative Reverse-Transcription PCR

The expression levels of Gpnmb and its targeting miRNAs (miR-125a-3p, miR-214-5p, and miR-3076-5p) were validated in D2 and D2G at 2 time points of 2 to 4 months (3 months) and 9 to 11 months (10 months) using qRT-PCR. Our results demonstrated significant (P < 0.05) downregulation of Gpnmb (Fig. 5A), whereas miR-125a-3p and miR-214-5p were significantly upregulated at the 3-month (P < 0.01) and the 10-month time points (P < 0.05) in D2 compared to D2G mice. Furthermore, miR-3076-5p was significantly upregulated at the 3-month time point (P < 0.01), however, its overexpression in D2 at the 10-month time point did not attain a statistical significance (Figs. 5B-D). 

Figure 5.

 

The relative expression levels of (A) Gpnmb, (B) miR-125a-3p, (C) miR-214-5p, and (D) miR-3076-5p in retina detected by qRT-PCR. **P < 0.01, and ***P < 0.001.

Discussion

In this study, we characterized the genetic environment associated with pigmentary glaucoma in D2 mice. We compared retinas of D2 with D2G mice and identified DEGs at both pre-and post-onset glaucoma stages. DEGs were predominantly found in the pre-onset glaucoma group. Ontology analysis revealed different biological themes between the differentially expressed genes in the pre- and post-onset glaucoma groups. DEGs identified in the pre-onset glaucoma group were primarily associated with eye structure and function related processes, whereas those identified in the post-onset glaucoma group were primarily associated with immune related biological processes. We identified three miRNAs that target Gpnmb and used the BXD family of mice to propose an miRNA-mediated Gpnmb gene network that influences the development of the D2 pigmentary glaucoma phenotypes. 

To further explore the relationship between these DEGs and glaucomatous damage, we conducted a search of the DEGs identified in our dataset on the Glaucoma Discovery Platform (http://glaucomadb.jax.org/glaucoma).^61,62 The analysis showed that among the 314 DEGs found in the pre-onset glaucoma group (1–6 months), 87 genes were also noted to be differentially expressed in the “no or early axons” damage group. This suggests many additional gene candidates may be involved in the early molecular stages of D2 glaucomatous disease. In the post-onset glaucoma group (6–13 months), out of the 86 DEGs identified, 71 genes were found to be differentially expressed in the “moderate or severe” axons damage group, whereas only 11 genes were noted in the “no or early axons” damage group. These findings indicate a strong relationship between the gene expression changes in the post-onset glaucoma stage and the extent of axon damage. Gene enrichment analyses provide key insight into the importance of the identified DEGs during each stage. 

Early Group Findings

Crystallins are the predominant structural proteins in the vertebrate eye lens that are evolutionarily related to stress proteins.⁶³ They can be subdivided into two major families: alpha and beta gamma. Alpha-crystallins, part of the chaperone-like small heat-shock proteins family, serve to protect against cellular stress.⁶⁴ Beta gamma crystallins are structurally related to microbial stress proteins and are involved in various metabolic and regulatory functions.⁶⁴ Crystallins once represented a group only associated with the lens, but now are part of an expanding family of genes/proteins that are expressed in several cell types associated with various ocular and neural structures.^65–67 

One major function of extralenticular crystallins is their response to cellular stress.^63,68 Alpha crystallins have been previously referred to as “guardians of the retina.”⁶³ Alpha-beta crystallin protects retinal pigment epithelial cells against stress-induced apoptosis.⁶⁹ Critical to photoreceptor outer segment membrane renewal, alpha crystallins are involved in transporting rhodopsin by interacting with post-Golgi membranes. Reduced crystallin expression with age may therefore result in accumulation of damaged proteins that may progress into neurodegeneration. Crystallins also have a role in regulating inflammation in the central nervous system (CNS). Several in vivo and in vitro studies have demonstrated that administration of alpha crystallins reduces inflammation and prevents neurodegeneration in models of optic nerve injury and ischemic optic neuropathy.^70–72 

In the pre-onset glaucoma group, all identified DEG members of the crystallin superfamily (Cryaa, Cryab, Cryba2, Cryba4, Crybb2, Cryga, Crygc-e, and Crygs) were significantly reduced in the glaucomatous retina of D2 mice compared with the D2G group. However, differential expression for the same genes was not observed in the post-onset glaucoma group. Other glaucoma rodent models have described the change in crystallin gene's response to IOP elevation over time. In Winstor rats, 2 weeks after IOP elevation, several members of the crystallin superfamily were downregulated with no difference at 5 weeks.⁷³ In Brown-Norway rats, a similar response pattern was demonstrated with downregulation for 8 days and no difference at 5 weeks.⁷⁴ Our studies emulate the findings in Steele et al. who found an upregulation of transcripts involved in the immune response along with a downregulation of crystallins in retinas from D2 mice at an advanced age.⁷⁵ Interestingly, in a recent large study comparing postmortem human glaucoma and control eyes, crystallins were identified as the most differentially expressed proteins and exclusively downregulated in the glaucoma subset.⁷⁶ A simplistic conclusion from these studies may be that downregulating crystallin gene expression increases RGC susceptibility to glaucomatous nerve damage in an age-dependent manner. A recent study induced experimental glaucoma in Sprague Dawley rats using episcleral vein occlusion and similarly found crystallins of all subclasses to have higher abundance in retinal samples of younger animals in response to IOP increase. The neuroprotective effects of three selected crystallins, CRYBB2, CRYGB, and CRYAB on RGCs isolated from retinal explants of untreated rats in vitro were examined; these crystallins reduced RGC loss in response to elevated pressure.⁷⁷ However, a key aspect to understanding the dynamic nature of crystallin gene expression lies in its connection to Gpnmb function in the D2 model. 

Late Group Findings and Innate Immune Response

The GO enrichment analysis demonstrated that the DEGs identified in the post-onset glaucoma group were primarily associated with the innate immune system, the initial line of defense against foreign organisms. Cells involved in the innate response do not confer long-lasting immunity against the same pathogen to the host. An additional component of the innate response relevant to glaucoma is the complement cascade, which acts to amplify the ongoing response to a pathogen. 

Several of the DEGs expressed included members of the complement cascade. Howell et al. demonstrated that in D2 mice deficient in C1qa, there was a significant decrease in RGC loss and optic nerve degeneration compared to regular D2 mice.⁶¹ C1q is a component of the C1-complex that triggers the classical pathway and is released primarily by microglia. Microglia are innate immune cells with sensor, effector, and phagocytic capacities. They are responsible for the primary active immune response in the CNS. Normally, microglia function to survey and respond to retinal damage by scavenging neuronal debris to limit damage and signaling to other immune effector cells through the secretion of local inflammatory mediators.^78–81 However, continuous activation of retina microglia may be related to RGC degeneration.⁸² Although our analysis demonstrated that DEGs are associated with microglial activation, primarily in the post-onset glaucoma group, previous studies have shown that microglial activation and recruitment are early events in the D2 model.^61,83 Bosco et al. reported that microglial aggregation and activation in the optic nerve head region was solely observed in the D2 mouse and preceded RGC injury.⁸³ Another study by the same group found a significant correlation between early microgliosis at the optic nerve head (ONH) and the late severity of nerve pathology.⁸² However, we found that D2G mice without the Gpnmb^R150X mutation, do not demonstrate early innate immune system changes. 

The presence of Gpnmb has been detected in bone-marrow derived cells of the monocytic lineage and has been proposed as an inflammation suppressor gene.^82,84,85 The resulting truncated form of GPNMB may influence the bone marrow and produce a pro-inflammatory macrophage phenotype. RT-PCR demonstrated significantly lower expression of Gpnmb in both the pre- and post-onset glaucoma groups when comparing D2 to D2G. Overall, differential expression of Gpnmb in both pre- and post-onset glaucoma groups, lower expression of Gpnmb in both the glaucoma stages in D2 mice, and our ontology analysis indicating innate immune system events support an underlying pathophysiological role of Gpnmb and pathologic microglial activation, inducing the glaucomatous nerve damage. 

Key miRNAs Identified and Their Roles in Neurodegenerative and Eye Diseases

The miRNA activity in the CNS has been demonstrated in several neurodegenerative diseases, including amyotrophic lateral sclerosis (ALS), Huntington's Disease (HD), and Alzheimer's Disease (AD).⁸⁶ Expression of miRNAs has been found to be associated with the development of various eye diseases, including ocular diseases, such as age-related macular degeneration⁸⁷ and retinoblastoma.^88–91 Increasing interest is being placed on the role of miRNAs in the pathogenesis of glaucoma.^92,93 We interrogated miRNA-Gpnmb interaction to better understand how miRNAs impact the Gpnmb gene network in glaucoma. We identified 3 miRNAs with significant differential expression between D2 and D2G: miR-125a-3p, miR-3076-5p, and miR-214-5p. Using WGCNA, we identified highly co-expressed gene sets and identified a miRNA-mediated Gpnmb-involved gene network. In fact, two of the identified miRNAs in this study have been evaluated in human vitreous humor samples. Ragusa et al. examined 18 vitreous humor samples from patients with choroidal melanomas, retinal detachment, or macular hole. Of the 94 circulating small RNAs in the vitreous, 10, including hsa-miR-214 and hsa-miR-125a-3p, demonstrated a particular abundance in the vitreous humor, and were either downregulated or not detectable in the serum of healthy donors.⁹⁴ This suggests that there are a specific set of miRNAs secreted in the eyes. 

The miR-125 family has been demonstrated previously to play important roles in cell differentiation, growth and apoptosis.⁹⁵ Our study implicated miR-125a-3p, which has been shown to be involved in the pathogenesis of several neurodegenerative conditions, including multiple sclerosis (MS) and Parkinson's disease.⁹⁶ In MS, miR-125a-3p was found to be upregulated in active lesions from patients as well as in the isolated spinal cord oligodendrocyte precursor cells in mice. In in vitro and in vivo mice models, upregulated miR-125a-3p impaired oligodendrocyte precursor cell differentiation.^97,98 Demyelination may contribute to the pathogenesis of glaucoma and precede axonal degeneration.^99–101 A recent study explored the relationship between the rs12976445 polymorphism in miR-125 and POAG severity.¹⁰² They compared 88 patients with POAG comprising 3 genotype groups (GG, GC, and CC) to evaluate the relationship between each group and the POAG index. It was found that the rs12976445 polymorphism was significantly associated with the risk of POAG via regulation of miR-125a and IL-6R expression. Interestingly, IL-6 has been noted to be involved in a variety of CNS pathologies, including neurodegeneration.^103,104 However, studies demonstrating the effect of IL-6 are variable and have demonstrated both potentially neuroprotective¹⁰⁵ and harmful¹⁰⁶ effects on RGCs. miR-214 may play an important role in neurogenesis. Smith et al. aimed to complete a large-scale expression analysis of different stages of human neurodevelopment using NTera2/D1 (NT2) cells, a “surrogate” of pluripotent embryonic stem cells that undergo a multistep transition from undifferentiated state to terminal neuronal differentiation.¹⁰⁷ During retinoic acid induced transition from progenitors to fully differentiated neural phenotypes, NT2 cells transition through distinct temporal phases beginning with stem cell to precursor expansion, neural commitment, and onset of cell-cycle exit. miR-214 was induced in this beginning phase and was identified as one of 10 known and predicted miRNAs to continue to be expressed in fully differentiated NT2- N neurons and/or NT2-A astrocytes. Metlapally et al. examined scleral miRNA expression profiles of rapidly growing human fetal eyes compared with stable adult donor eyes. They demonstrated increased expression of miR-214 among others in fetal sclera regardless the site of tissue collection.^108,109 

There remain key challenges involved in identifying the role of these miRNAs in regulating Gpnmb. Well-designed time point loss-of-function studies using these miRNAs in the D2 mouse may reveal additional information about their function in regulating Gpnmb. 

In the present study, morphological data were not collected in parallel with gene expression data for each mouse. However, we examined data previously obtained from the D2 colony (at UTHSC; Supplementary Fig. S1) to show phenotype changes over time. Although historic data does provide valuable insight into the mice used in this study, we acknowledge the ecological fallacy that data derived from the colony may not necessarily be applicable to individual mice used in the current study. We aim to validate these findings at the individual level in future work. 

In this study, we investigate the genetic architecture surrounding the Gpnmb^150X mutation and the glaucoma phenotype using D2 mice retinas. We hypothesize that early involvement of inflammation propagates due to diminishing neuroprotective effects from the downregulation of crystallin gene expression in RGCs and upregulation of innate immune system genes later in development. A combination of these events may make RGCs more susceptible to both inflammatory-mediated and IOP-induced damage. Using miRNA target prediction tools, we also identify three high priority miRNA candidates that potentially target Gpnmb. Utilizing systems genetics, we harness the power of WCGNA to elucidate underlying molecular mechanisms by identifying hubs and modules associated with Gpnmb and our miRNA candidates. In future studies, we plan on further exploring this network and how it relates to the regulation and development of the D2 glaucoma phenotype. 

Acknowledgments

Supported by NEI grant EY021200 and Challenge Grant from Research to Prevent Blindness. The funder had no role in the design, execution, or interpretation of the study, or the decision to submit the manuscript for publication. 

The authors of this article have no financial or commercial relationships with any organizations that could have influenced the content of this manuscript. 

Disclosure: Q. Gu, None; A. Kumar, None; M. Hook, None; F. Xu, None; A.K. Bajpai, None; A. Starlard-Davenport, None; J. Yue, None; M.M. Jablonski, None; L. Lu, None