Molecular clustering identifies complement and endothelin induction as early events in a mouse model of glaucoma

Hierarchical clustering allows early disease stages to be identified.

Comparison of eyes grouped by conventional morphological criteria did not allow for sensitive detection of early gene expression changes. We reasoned that grouping the same set of eyes by molecular clustering of gene expression profiles and then comparing these groups would allow more sensitive detection of early changes in glaucoma. Since a critical insult and early damage are known to occur in the optic nerve, clustering was first used to group eyes based on the similarity of gene expression profiles in the ONH. Unbiased hierarchical clustering of the 40 DBA/2J ONH samples and 10 D2-Gpnmb⁺ control samples was performed using the gene expression values of 570 disease-relevant probe sets (Methods and Supplemental Figure 1; supplemental material available online with this article; doi: 10.1172/JCI44646DS1). Hierarchical clustering first linked the most similar eyes, and continued until all eyes were linked in a dendrogram (Figure 1B and Methods). As the least similar eyes were linked at the highest branch points, a cutoff (threshold of relatedness) was applied to separate groups of eyes into independent disease stages (Methods). We selected a cutoff that made biological sense based on conventional morphological criteria, e.g., Gpnmb⁺ controls were grouped into their own stage, and SEV eyes were separate from NOE eyes. We only considered groups containing at least 4 eyes as a stage. Using these criteria, 5 new, molecularly defined, stages of glaucoma were identified (dataset 2, stage 1 to stage 5; Figure 1, B and C). Not surprisingly, 10 of the 40 eyes did not fit into these 5 stages. This reflects the large number of possible molecular states for such a variable disease.

The 5 stages were ordered based on two criteria: (a) Number of DE genes compared with no-glaucoma control: stage 1 had the smallest number of DE genes when stages 1–5 were compared with the Gpnmb⁺ control group (381 DE genes, Figure 1C) and was considered to be the earliest molecular stage. (b) Previously determined morphological damage: although stages 2–4 had a similar number of DE genes (approximately 10,000), only stages 1 and 2 contained eyes with no detectable glaucoma (NOE). In comparison, stage 3 contained 3 eyes with NOE glaucoma and 3 eyes with MOD glaucoma. Stage 4 contained 4 eyes with MOD glaucoma. Therefore, stage 3 was determined to be between stage 2 and stage 4. Stage 5 contained only eyes with SEV glaucoma. As expected, axon number decreased with increasing glaucoma stage (Figure 2A). Also, the number of genes whose expression value was DE by at least 2 SDs from the average of Gpnmb⁺ controls increased from stages 1 to 5 (Spearman’s correlation 0.9, P < 10^–100) (Figure 2B).

Based on optic nerve damage, 3 of the 6 eyes in stage 3 were indistinguishable from no-glaucoma control eyes. The other 3 eyes had MOD glaucoma (Supplemental Figure 2). This indicates that the molecular changes determining the similarity of eyes in stage 3 must have occurred prior to detectable glaucomatous damage. In addition to RGCs, these early molecular changes may be occurring in a variety of other cell types in the ONH, such as astrocytes, microglia, and endothelial cells.

Principal component analysis (PCA) is a powerful method for reducing the dimensionality of data and is used to facilitate microarray experiments (31). Here, PCA was used to independently assess whether the 5 identified molecular stages are distinct. The reduced dimensional position of the 40 DBA/2J and 10 control eyes was plotted using the first two principal component vectors calculated for the 570 disease-relevant probe sets (see Methods). As expected, eyes in each of the 5 molecular stages occupied essentially distinct territories (P < 0.01, Figure 2C). Therefore, PCA supported the hierarchical clustering.

Molecularly defined ONH stages increase sensitivity to detect DE genes in the ONH.

To identify DE genes, molecularly defined stages 1–5 were compared with the Gpnmb⁺ control group (Figure 3). DE genes were defined as being those with a q value (false discovery rate [FDR]) of 0.05 or less. Significantly more DE genes were detected in the ONH for the molecularly defined stages than for the morphological groups (Figure 3A). Considering groups with eyes that had no detectable optic nerve damage and only probe sets that were DE at least 2-fold, 1,385 probe sets were DE in molecularly defined stages 1–3 versus 100 in morphologically defined NOE2 eyes (Figure 3A) and 55 in NOE1 eyes (data not shown). Additionally, many specific genes, such as Timp1, Edn2, Vcan, and Calcb, had a greater fold change in the molecularly as compared with morphologically defined stages (Figure 3B).

ECM and immune changes in the ONH.

We used a variety of software-based annotation tools to provide functional insight into the ONH expression differences detected in this study (see Methods). As it is not possible to show all analyses in this article, we highlight a few enriched gene ontology terms, pathways, and networks that are likely to be important during early stages of glaucoma.

Gene ontology (GO) analysis (performed by DAVID; see Methods) revealed that a significant number of genes with the GO terms “immune response” (Gene Ontology Database 0008283), “leukocyte activation” (GO 0045321), and the related term “chemotaxis” (GO 0006935) were DE in the early molecularly defined stages of disease (Figure 3C). The majority of these changes were not detected in the morphologically based NOE groups. Genes belonging to pathways that are known to affect these same processes were also identified by independent Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis (Figure 3D). Two related KEGG pathways that were significantly overrepresented are “ECM-receptor interactions” (KEGG mmu04512; Figure 4, A and B) and “focal adhesions” (KEGG mmu04510). Neither of these pathways was found to be important in either of the NOE groups or the MOD group (dataset 1). ECM-receptor interactions are mediated by transmembrane molecules such as integrins and proteoglycans and lead to direct or indirect control of cellular activities that include adhesion, migration, and proliferation. These pathway changes may represent early changes to the glial-vascular unit or altered interactions between glia and neurons in the ONH. One member of this pathway, tenascin C (Tnc) is 1.9-fold elevated in stage 2, 4.1-fold elevated in stage 3, and 5.7-fold elevated in stage 4 (Figure 4B) and is localized to the glial lamina region of the ONH (Figure 4C). Showing the human relevance of our findings, expression of TNC was increased in glaucomatous compared with normal human ONHs (32). The importance of ECM-receptor interactions is also supported by computationally distinct k-means clustering, which groups genes that behave similarly across disease stages (Table 1 and Figure 5). K-means clustering also suggested that the complement cascade is an important pathway activated early in the ONH. C1qa is a component of the C1q complex, an early effector in the complement cascade. C1QA is expressed in microglia/macrophages in the ONH during early stages of glaucoma (Figure 6).

Hubs and transcription factors.

Ingenuity pathway analysis (IPA) identifies molecular networks and biological pathways that are overrepresented in a given set of DE genes (see Methods). IPA allows the identification of key genes (or hubs) that lie at the center of any network or subnetwork. Therapeutic modulation of any hub has strong potential to impact a variety of network participants, making hubs attractive targets for future neuroprotective efforts. We separately analyzed the DE genes for stages 1–5 (dataset 2) by IPA. One of the most statistically significant networks included various members of the ECM-receptor interaction pathway and proinflammatory cytokines (Supplemental Figure 3A). Hubs within this network include IL-1β (encoded by Il1b, shown as IL1 in IPA network) and integrin β3 (Itgb3), molecules that are reported to affect ECM composition and metabolism (33–35). The network also includes various other integrins, fibronectin 1 (Fn1), matrix metalloproteinase 19 (Mmp19), tissue inhibitors of matrix metalloproteinases (Timp1 and Timp2), and a number of collagens. Of interest, caspase-1 (Casp1) was upregulated in the earliest glaucoma stage and continued to be increased at later stages (Supplemental Figure 3B). CASP1 has interleukin-converting enzyme activity, and its increase preceded the upregulation of Il1b and Itgb3 genes (Supplemental Figure 3B). Thus, proinflammatory responses are among the earliest changes in the ONH in glaucoma.

Clustering identifies early stages of glaucoma in the retina.

Clustering eyes based on morphological criteria (dataset 1) or the similarity of ONH gene expression profiles (dataset 2) was not sensitive for identifying early changes in the retina. No DE genes were identified in the retina in molecular stages 1 and 2, and only 8 DE genes were identified in stage 3 (Figure 1C). This may be because early changes in the ONH occur prior to changes in the retina, and/or because our sensitivity to detect retinal changes was lower as whole retinas were used (<5% of total retinal cells are RGCs). Alternatively, it is possible that grouping the eyes based on ONH expression patterns is not a sensitive way to determine early glaucoma changes in the retina. It is possible that early expression changes in the retina are completely or partially independent of events in the optic nerve. Thus, we separately clustered the same set of 40 DBA/2J eyes based on the similarity of their retinal gene expression profiles (dataset 3). We used the same methods as were used for the ONH-based clustering and similar criteria for determining stages (see Methods). Retinal clustering divided the eyes into a Gpnmb⁺ control cluster and 4 glaucoma stages (R1–R4) that made biological sense (dataset 3, Figure 7, A and B). Stages R1 and R2 contained only eyes that were previously indistinguishable from each other and controls by conventional morphological criteria. Retinal clustering increased the sensitivity of detecting DE genes at early stages of disease (Figure 7B). Compared with the D2-Gpnmb⁺ control stage, 48 probe sets were DE in stage R1 and 8,664 in stage R2. Many fewer DE genes were detected in the ONH of the eyes in these retinal clusters (59 DE genes, dataset 3, stage R1, compared with 381 DE genes, dataset 2, stage 1; Figure 1C). Since clustering using either the ONH or retina tissues is not sensitive at detecting early DE genes in the other tissue, our analyses suggest that events in the retina and ONH occur asynchronously and are not completely interdependent.

Analyses of the 48 DE probe sets in stage R1 identified the complement cascade as the only significantly overrepresented pathway (P = 5 × 10^–4). The DE genes in stage R1 included C1qa, C1qb, and C1qc. In stage R2, 12 members of the complement cascade were DE (Figure 7C). We have found that C1qa is expressed in RGCs and localizes to synapses in the retina during early stages of glaucoma (36). The complement cascade was also implicated in the ONH by k-means clustering (Table 1 and Figure 5, cluster 6), where it was expressed in microglia (Figure 6).

Early DE genes provide molecular markers to assess glaucoma status.

Genes shown to be DE early in glaucoma may allow the development of marker sets to group eyes with no detectable optic nerve damage into those undergoing early glaucoma and those with no glaucoma. Such markers may also be used to determine the degree of glaucoma or to assess the effectiveness of tested treatments. Since many molecular changes occur in glaucoma and specific changes will vary from eye to eye, a reasonable approach is to assess genes from different pathways and to track the number of DE genes for each of these pathways as a measure of glaucoma status.

To assess the feasibility of this approach, we assessed genes from pathways that were activated early in the ONH or retina in the microarray study (Supplemental Tables 2 and 3). These pathways included ECM-receptor interactions, MAPK signaling, and Toll-like receptor signaling for the ONH tissue and the complement cascade for the retina tissue. We assessed completely new sets of 10-month-old DBA/2J and D2-Gpnmb⁺ control eyes. As we were most interested in eyes with no conventionally detectable glaucoma, we separately assessed the ONH and retina for 21 NOE and 3 MOD DBA/2J eyes, as well as 11 D2-Gpnmb⁺ control eyes. To further refine the location of retinal changes compared with the microarray study, the retinal tissue was enriched for the RGC layer and lacked the inner nuclear layer and photoreceptor layers (see Methods). For each eye, a gene was classified as DE only if the normalized expression level was at least 2 SDs from the average of the 11 control eyes. This approach ensured that only robust expression changes were regarded as DE. For each eye, we determined the total number of DE genes as well as the number of DE genes in each pathway (Supplemental Figure 5). A pathway was considered “activated” in eyes that had at least 2 times the number of DE genes for that pathway relative to any of the controls.

In both ONH and retina, this analysis clearly distinguished the vast majority of eyes with NOE glaucoma from controls (Supplemental Figure 5A). Individual eyes had differing numbers of DE genes, with eyes classified with MOD glaucoma having the greatest total number of DE genes. For the ONH, ECM-receptor interactions, MAPK signaling, and Toll-like receptor signaling pathways were activated in 18 of the 21 eyes with NOE glaucoma and all 3 eyes with MOD glaucoma. For the retina, the complement cascade was activated in 13 of 21 NOE eyes and 3 of 3 eyes with MOD glaucoma. In the NOE eyes, there was no clear relationship between the number of DE genes in the ONH and retina. Most NOE eyes had high numbers of DE genes in the ONH, with lower numbers activated in the retina (16 of 21 having ≥60% the number of genes activated in moderate for ONH compared with only 3 of 21 for retina). However, the reverse was true in some eyes. Most strikingly, eye number 3 had a large number of DE genes in the retina (similar to eyes with MOD glaucoma), but a low number in the ONH (Supplemental Figure 5). This further suggests that early changes in the retina and ONH are not completely dependent on each other and can occur asynchronously. These marker genes (Supplemental Table 2), along with approaches for assessing glaucoma status in individual eyes, are likely to be valuable for future studies.

C1QA-deficient mice are protected from glaucoma.

Changes in the complement cascade have been reported in animal models and human glaucoma. However, it has not been determined whether these changes occur very early and prior to detectable axon loss, or if they are later and possibly secondary (22, 24, 36, 37). Our current study demonstrates that complement cascade changes occur very early and prior to detectable glaucoma by conventional assays, with complement genes being among the first genes to change in the retina (Figure 7C). Given the very early occurrence of these changes and the prior association of complement cascade changes with human glaucoma, we tested the importance of this cascade in DBA/2J mice.

To test the functional importance of the complement cascade in this inherited glaucoma, we generated and analyzed mice with a mutation in the C1qa gene (Figure 8). Clinical examinations of mice from 6.0 to 12.0 months of age assessed the severity of the pigment-dispersing iris disease, which precedes the glaucoma (30, 38). Both C1QA-deficient and wild-type DBA/2J mice developed prominent pigment dispersion and iris atrophy with a similar age of onset, but progression to the most severe iris stage was slightly delayed in C1QA-deficient mice. Subsequent to the iris disease and despite a modest delay in some C1QA-deficient DBA/2J eyes, IOP became elevated in mice of each genotype, and their IOP distributions overlapped extensively (Figure 8A). We previously detected an alteration of IOP similar to that in the C1qa mutants in BAX-deficient DBA/2J mice (39). In these Bax mutants, this IOP change did not profoundly alter the degree of optic nerve damage.

C1qa mutant mice were protected from glaucomatous damage to the optic nerve and retina at both 10.5 and 12.0 months of age (Figure 8B), two key ages in this strain (40). Importantly this protection was very profound at 10.5 months of age, with only 9% of mutant eyes having detectable glaucoma compared with 63% of wild-type eyes (P = 5 × 10^–100, Figure 8B). None of the mutant eyes had SEV glaucoma, while 48% of wild-type eyes were severely affected. At 12.0 months of age, only 25% of mutant eyes had SEV glaucoma compared with 48% of wild-type eyes. The percentage of eyes with no detectable glaucoma was more than doubled in C1qa mutants at 12.0 months (66% versus 30%, P = 2.2 × 10^–9). Protection in C1qa mutant mice was confirmed by axon and soma counts (Figure 8, C–E). Thus, although future experiments are required to completely define the roles of C1QA in ocular drainage and neural tissues during DBA/2J glaucoma, it is likely to have an important role in glaucomatous neurodegeneration. Irrespective of the relative role in different tissues, our findings indicate that complement pathway inhibition has potential as a potent treatment to protect against glaucoma.

Inhibiting the endothelin system protects from glaucoma in DBA/2J mice.

Our current study demonstrates that the Edn2 gene is upregulated in both the retina and the ONH very early in glaucoma (Figure 3B). The endothelin system was previously found to be upregulated in human and animal models of glaucoma (26, 41–43), but our study indicates that the upregulation of Edn2 is a very early event. EDN2 localized to microglia/macrophage-like cells in the nerve fiber layer of DBA/2J retinas (Figure 9, A–D).

Administration of EDN1, which shares receptors with EDN2, is known to kill RGCs. In contrast, EDN2 was suggested to protect photoreceptors from stress (44). To initially assess whether EDN2 damages RGCs, we injected EDN2 peptide into young DBA/2J eyes and found that it damaged both RGCs and the optic nerve (Figure 9, E and F). EDN1 is reported to be harmful to neurons by constricting blood vessels, triggering reactive astrocytosis, and/or inhibiting axon transport, but effects on astrocytes and axon transport may be secondary to altered blood flow (45–47). Since EDN2 is a potent vasoactive peptide, we hypothesized that EDN2 induces vasoconstriction in the ONH and retina during glaucoma. Even if subtle, this vasoconstriction may metabolically stress RGCs and contribute to their demise in glaucoma. In support of this vasoconstriction hypothesis, we identified a decrease in the ratio of lumen area to total blood vessel area in retinas of DBA/2J eyes at early and later stages of glaucoma (see Methods, Figure 9, G and H, and Supplemental Figure 6). This decrease was not uniform in any of the eyes investigated, suggesting local differences that may affect some but not all RGCs. The decreased lumen area is likely a result of vasoconstriction, but an EDN2-induced thickening of the vascular smooth muscle (48) cannot be ruled out. Either way, decreased lumen area is likely to negatively impact perfusion and RGC survival.

To functionally assess the role of the endothelin system in glaucoma, we administered bosentan, an endothelin receptor antagonist, to DBA/2J mice starting at 6.0 months of age (Figure 10). We selected bosentan because it antagonizes both types of endothelin receptors and does not alter blood pressure (see Methods). Also, bosentan has been shown to increase ocular blood flow in human glaucoma patients (49). Administration of bosentan did not alter the onset or progression of either the iris disease or IOP elevation (Figure 10A). However, bosentan significantly reduced glaucoma at both 10.5 and 12.0 months of age (Figure 10, B–E). The protection was especially strong at 10.5 months of age, when 80% of treated eyes had no detectable glaucoma, compared with only 39% of untreated eyes (P = 1.7 × 10^–26). At 12.0 months, 43% of treated eyes had no detectable glaucoma compared with only 31% of untreated eyes (P = 0.01). These experiments implicate early activation of the endothelin system as pathogenic in this inherited model of glaucoma. Endothelin receptor antagonists offer promise as new treatments to alleviate glaucoma.

Mouse strains, breeding, and husbandry.

All experiments were conducted in accordance with the Association for Research in Vision and Ophthalmology statement on the use of animals in ophthalmic research. All mice were housed in a 14-hour light/10-hour dark cycle under previously described conditions (60). The Jackson Laboratory Animal Care and Use Committee approved all of the experiments in this study. All mice in the gene expression profiling study were female and housed under the same conditions in same animal facility, minimizing environmental differences. Our DBA/2J (D2) mouse colony is routinely crossed with DBA/2J mice from The Jackson Laboratory production facility to prevent genetic drift. The D2-Gpnmb⁺ strain (DBA/2J-Gpnmb⁺/Sj), which rarely if ever develops glaucoma in our colony, was used as a control strain for DBA/2J glaucoma studies (50). For the gene expression analyses, more than 500 mice were aged and analyzed to produce sufficient eyes with each glaucoma stage. To generate D2.C1qa^+/– mice, the C1qa^+/– allele (61) was backcrossed to D2 mice for 8 generations. The C1qa mutation was selected due to its established role in retinogeniculate refinement (36). D2.C1qa^+/– mice were assessed with 102 genetic markers evenly distributed across the mouse genome, and only one marker closely linked to the C1qa gene was positive for C57BL/6J (this represents flanking sequence closely linked to the C1qa mutation). The remaining 101 markers were positive for DBA/2J. This indicates that greater than 99% of the genetic background of D2.C1qa^+/– mice is derived from DBA/2J. To generate DBA/2J mice deficient in C1QA and littermate controls, we intercrossed D2.C1qa^+/– mice.

Analysis of glaucomatous damage.

Intracranial portions of optic nerves were processed and analyzed as previously described (38–40). Briefly, optic nerves were fixed in situ for 48 hours, dissected free, processed, and embedded in plastic. One-micrometer-thick cross sections of optic nerve from behind the orbit were cut and stained with paraphenylenediamine (PPD). PPD darkly stains the myelin sheaths and axoplasm of sick or dying axons, but not healthy axons (62). PPD staining allows sensitive determination of disease stage for large numbers of nerves. Three observers blinded to the experimental protocol determined the degree of nerve damage as previously described (38–40, 63). The 3 clearly distinguishable stages of damage were as follows: (a) NOE — less than 5% axons damaged. This amount of damage is seen in age-matched mice of various strains that do not develop glaucoma. It is called “no or early” because some of these eyes are undergoing early molecular stages of disease, but they are not distinguishable from eyes with no glaucoma by conventional morphologic analyses of axon number and nerve damage. (b) MOD — many damaged axons throughout the nerve, with an average of 30% axon loss. (c) SEV — substantial axon loss (>50% lost) and damage. To ensure very distinct stages for this study, we selected severely affected eyes with very few remaining axons (>90% loss). For C1qa and bosentan experiments, differences in optic nerve damage distributions between groups were compared using χ² tests.

Tissue harvesting, RNA isolation, and processing for microarrays.

Following anesthesia, the left eye was dissected free, and the cornea and lens were removed. The optic nerve was cut close to the back of the eye. The ONH was excised using a sharpened 1-mm glass capillary centered over the optic nerve. The punched-out nerve head was immediately placed in RNAlater (Ambion). A small portion of central retina as well as a small portion of choroid and sclera were included in this tissue. This tissue is substantially enriched for the glial lamina. The retina (containing no ONH) was then dissected free from the back of the eye and placed in a separate vial containing RNAlater and stored at 4°C until required. The process was repeated for the right eye. Dissection of each eye was completed in 5 minutes or less (10 minutes per mouse).

For RNA preparation, tissue was homogenized in RLT Buffer (QIAGEN). Total RNA was isolated using the RNeasy Micro Kit (QIAGEN) according to the manufacturer’s protocols, and quality was assessed using a 2100 Bioanalyzer instrument and RNA 6000 Nano LabChip assay (Agilent Technologies). Tissue harvesting and RNA preparation were completed for all mice aged for the gene expression studies. Samples with only the highest-quality RNA were used.

RNA labeling and microarray processing.

Biotin-labeled cDNA was synthesized separately from 20 ng of each total RNA sample per the manufacturer’s protocols (Ovation Biotin System, NuGEN). 2.5 μg of each biotin-labeled and fragmented cDNA sample was then hybridized onto Mouse Genome 430 version 2.0 GeneChip arrays (Affymetrix). Post-hybridization staining and washing were performed according to the manufacturer’s protocols using a Fluidics Station 450 instrument (Affymetrix). Finally, the arrays were scanned with a GeneChip Scanner 3000 laser confocal slide scanner. The images were quantified using GeneChip Operating Software (GCOS) version 1.2. Probe level data were imported into the R statistical software environment ( http://www.r-project.org), and expression values for each probe set on the array were summarized using the robust multichip average (RMA) method with quantile normalization in the R/affy package (64, 65).

A total of 110 arrays were used, assessing 60 ONHs and 50 retinas. The ONHs included 10 D2-Gpnmb⁺ control eyes at 10.5 months, 40 DBA/2J eyes at 10.5 months, and 10 DBA2J eyes at 4.5 months. The 50 retinas included 10 D2-Gpnmb⁺ control eyes at 10.5 months and 40 DBA/2J eyes at 10.5 months. Only 2 arrays failed the quality control; both were 10.5-month DBA/2J MOD retinas. All 110 arrays were processed at the same time. To allow future comparisons, quantile normalization of all of the arrays was performed. All individual gene expression profiles are available in GEO datasets (NCBI, accession number GSE26299).

Statistics for differential gene expression.

Differential gene expression between any two morphologically defined groups or molecularly defined stages of glaucoma was assessed by gene-specific fixed-effect ANOVA methods using the R/maanova package (64). The generalized 1-way ANOVA model can be written as Y_i = μ + GROUP + ε_i, where Y_i represents the log-transformed expression measure of probe set i, μ is the mean for each array, GROUP represents the fixed effect under comparison, and ε_i captures random error absorbing both array and biological components. Given the ANOVA model, all statistical tests were performed using Fs, a modified F statistic incorporating shrinkage estimates of variance components (66). All P values were calculated by permuting model residuals 1,000 times. To address the problem of multiple testing, the permutation P values were adjusted using the FDR method of Storey (67) implemented in R/q value. The resulting q values estimated the proportion of false positives in the list of DE genes and were used to select suitable candidate gene lists for subsequent analysis. Unless otherwise noted, statistical significance was determined at q ≤ 0.05 (FDR <5%). After statistical assessment, probe sets were annotated using Netaffx software (Affymetrix). A wide selection of the DE genes were confirmed to be DE by real-time PCR, including genes with a fold change greater than 2.

Defining early stages of disease using hierarchical clustering.

To generate a list of disease-specific probe sets to cluster the ONH samples, we first made a pairwise comparison between 10 non-glaucomatous D2-Gpnmb⁺ and 40 glaucomatous DBA/2J mice (which have the Gpnmb^R150X mutation) that were 10.5 months of age (Supplemental Figure 1A). As all of these mice were age-, sex-, and strain-matched, genes on this first list were DE due to effects of glaucoma or Gpnmb genotype (glaucoma and genotype list). To avoid DE genes due to Gpnmb genotype, we next generated a genotype-related list of probes that were DE between the 10.5-month-old D2-Gpnmb⁺ and the ten 4.5-month-old DBA/2J ONH samples. Importantly, since neither 4.5-month-old DBA/2J nor D2-Gpnmb⁺ mice have glaucoma in our colony, this second genotype list will not contain genes that are DE due to glaucoma. However, it may contain some genes that differ due to age. The final disease-specific gene list was then generated by subtracting the probe sets on the genotype list from the glaucoma and genotype list (570 DE probe sets, q < 0.001). The 570 probe sets were used to perform hierarchical clustering of ONH samples (using JMP Statistical Discovery version 7.0 software; SAS). Using the average distance metric, hierarchical clustering started with each eye as its own cluster. At each step, the clustering algorithm calculated the distance between each cluster and combined the two clusters that were closest together. This combining continued until all the points were in one final cluster or dendrogram. A threshold or cutoff, based on known differences between control, MOD, and SEV groups, was applied to identify 5 molecular stages that made biological sense. Since the filters that we used to identify disease-specific probe sets were not perfect, all genes other than Gpnmb were included in subsequent group comparisons. PCA was performed using MultiExperiment Viewer (TM4 version 4.5.1; www.http://www.tm4.org/mev) on the 10 D2-Gpnmb⁺ and 40 DBA/2J eyes using the 570 probe sets used for hierarchical clustering.

For the retina, 600 probe sets were DE (q < 0.01) when comparing the 10 D2-Gpnmb⁺ control samples with the 38 DBA/2J retina samples (2 of the 40 retina samples failed quality control). No subtraction was performed on the retina samples, as only Gpnmb was DE when D2-Gpnmb⁺ eyes were compared with either of the NOE groups at 10.5 months. Since Gpnmb clearly differed due to genotype, this gene was not counted as DE due to glaucoma.

A second clustering experiment was performed using probe sets DE between control and MOD eyes to cluster eyes with NOE glaucoma. Providing confirmation of our first clustering experiment, hierarchical clustering placed the NOE eyes in the same stages (data not shown).

Functional analyses and gene set enrichment of DE genes.

The Database for Annotation, Visualization and Integrated Discovery ( DAVID, http://david.abcc.ncifcrf.gov/; ref. 68) was used to add functional annotation to gene lists and provide statistical assessment of the annotations. We focused on the KEGG pathways (69) and GO terms (70). For each list of DE genes, DAVID determines the number of genes in each KEGG pathway and uses a Fisher exact test with Benjamini correction to determine the probability that the number of genes in each pathway would have occurred by chance (68). A corrected P value of 0.05 or less was used to identify significant pathways. Pathways with gene expression changes represented by colors (red: upregulated, green: downregulated) were generated using the KEGG mapper tool ( http://www.genome.jp/kegg/tool/color_pathway.html) (69).

We performed k-means clustering in R/maanova to identify clusters of probe sets whose expression changed in a similar pattern when considering stages 2–5 from the molecularly determined ONH groups. The earliest stage (stage 1) was not included, as there were too few DE probe sets in this stage. The normalized expression values (generated by RMA) of 4,138 probe sets (DE in each of the stages 2, 3, 4, and 5) were clustered into 8 groups (see Supplemental Table 2). Figure of merit analysis was used to determine the most meaningful number of clusters (see Supplemental Figure 4 and ref. 71). We performed k-means clustering using k = 4, 6, 8, 10, and 12. Based on pathway enrichment analysis, 4, 10, and 12 were not appropriate, and there was little difference between the results for 6 and 8. Therefore, we show the results for k = 8.

To assess gene networks and hubs, we used IPA (Ingenuity Systems).

Quantitative real-time PCR.

A completely independent set of eyes from 10-month-old mice were obtained. To avoid RNA degradation, two investigators harvested the eyes of each mouse in parallel. All dissections were performed in an RNase-free environment. Eyes were dissected free (including a significant portion of the myelinated optic nerve) and immediately placed in sterile 1× PBS. The retina with ONH attached was removed from the eye. To separate ONH from retina, an incision was made in the retina as close to the ONH as possible. Choroidal and scleral tissues, as well as the distal portion of the optic nerve, were removed from the ONH (approximately 100 μm of the myelinated nerve was left attached to the ONH). The ONH tissue was then placed immediately into RNA isolation buffer (Ambion) and stored at –20°C. The retina was then placed flat (ganglion cell layer up) onto the sterile membrane in a culture plate insert (Millipore, PR01715) ensuring no bubbles between the retina and membrane. The membrane was cut from the insert and placed flat (retina side up) onto a previously prepared flat surface of frozen OCT (Tissue-Tek). Freezing was instantaneous. The uppermost 80 μm of retina, containing the vast majority of the ganglion cell layer and some inner plexiform layer, was isolated by cryosectioning (4 × 20–μm sections) and immediately placed in RNA isolation buffer and stored at –20°C. For each mouse, the total dissection time never exceeded 15 minutes. RNA was extracted using an RNAeasy Micro Kit (QIAGEN), and RNA quality was assessed using a bioanalyzer (Agilent). Only samples with more than 5 ng of material with no sign of degradation were subjected to further analysis. cDNA synthesis was performed using the WT-Ovation Pico RNA Amplification System (NuGEN). In total, 12 D2-Gpnmb⁺ control and 23 DBA/2J eyes passed quality control.

Expression levels for 91 genes (69 ONH and 22 RGC layer of retina) and 6 normalizers were assessed. The normalizer genes were selected because they had the most consistent expression level, irrespective of damage level, among a previously tested panel of 194 genes (tested in DBA/2J retina and ONH; data not shown). The normalizers were 5730453I16Rik, NM_080456.1, and Slc11a2 for RGC layer and Etfb, Odc1, and Rpn2 for ONH. All primers were designed using Primer3 (version 0.4.0; http://frodo.wi.mit.edu/primer3/) with product sizes of 100–120 bp and optimum melting temperature of 60°C. All sequences are available as supplemental material (Supplemental Table 3). Real-time PCR reactions were prepared using Quant-iT RiboGreen RNA assay kit (Invitrogen) and run on the 7900HT Fast Real-Time PCR System (Applied Biosystems).

For each sample, each gene was interrogated in triplicate. Any technical replicate that was greater than 1.5 cycle threshold (CT) values away from the other 2 replicates was removed. For each gene, the CTs were calculated as the average of the successful replicates. For normalization, ΔCT values were calculated as CT (gene of interest) – geometric mean of CTs for the normalizers. The average and standard deviation of the D2-Gpnmb⁺ control eyes were calculated. The ΔΔCT was defined as the ΔCT (gene of interest) – ΔCT (control average). The fold change was calculated as 2^–ΔΔCT (upregulated) or –2^ΔΔCT (downregulated). A gene was considered DE if the fold change was greater than 2 SDs away from the average fold change for the control eyes. One control eye was removed from the study because it clearly was a far outlier.

RNA in situ hybridization and immunofluorescence.

DBA/2J.Thy1-CFP mice were used to visualize RGCs (15). Immunofluorescence to visualize EDN2 and IBA1 was performed as follows. Eyes from DBA/2J.Thy1-CFP mice were fixed overnight at 4°C in 4% PFA. Retinas were dissected free and incubated in primary antibody diluted in antibody buffer (5% BSA, 1% Triton X-100) for 5 days, with shaking at 4°C. Primary antibodies used were anti-EDN2 (Santa Cruz Biotechnology Inc., 1:100) and anti-IBA1 (Wako, 1:200). Retinas were washed in antibody buffer for at least 6 hours, with shaking at 4°C (changing wash solution at least 5 times). Retinas were incubated in secondary antibodies (Invitrogen, diluted in antibody buffer, 1:1,000) for a further 5 days, with shaking at 4°C. Retinas were rinsed in antibody buffer, washed in 1× PBS, mounted onto a slide in mounting media, and coverslipped. Fluorescence was visualized using a SP5 confocal microscope (Leica).

RNA in situ hybridization was performed as previously described, using 4% PFA perfusion-fixed, 12-μm-thick optic nerve cryostat sections (72). Digoxigenin-labeled (DIG-labeled) riboprobes for C1qa and Tnc were transcribed from cDNA clones (Open Biosystems clone ID: 3592169 and 40099589, respectively). The plasmids were digested with EcoRI (C1qa) and with SpeI (Tnc), and in vitro transcription was performed with T7 polymerase. T3 polymerase was used to generate the control sense probes. For both C1qa and Tnc, the sense probes produced no signal. The detection of hybridized mRNA in sections was performed using the Cy-3 Tyramide Signal Amplification System (PerkinElmer). After in situ hybridization, the sections were incubated in the primary antibodies: mouse anti-phosphorylated neurofilament (2F11; 1:1,000; Dako) and rabbit anti-IBA1 (1:500, Wako). Primary antibodies were diluted in a solution of 10% normal goat serum, 0.5% Triton X-100, and 0.5% BSA in 0.1 M PBS. For the secondary antibodies, goat anti-mouse Alexa Fluor 647 and goat anti-rabbit Alexa Fluor 488 were used at a 1:200 dilution (Invitrogen). The sections were then incubated with DAPI (Invitrogen) and mounted in Fluoromount (Sigma-Aldrich). Fluorescence was visualized using a SP5 confocal microscope.

For all protein and RNA localization, at least 3 sections (ONH) or whole retinas from 10 DBA/2J eyes and 5 D2-Gpnmb⁺ control eyes were assessed.

Vessel measurements and intravitreal injections of EDN2.

Retinas were fixed in 4% PFA overnight and placed flat on a slide and coverslipped. Major arterioles were measured from 24 imaged fields that were evenly distributed around the retina and 1 mm out from the ONH. Images were taken using differential interference contrast (DIC) on an SP5 confocal microscope. For each vessel, the diameter of the both the lumen and the outer vessel wall (taken from widest part of the vessel present in the field) was calculated as the average of 3 measurements. The ratio of the areas of the lumen to the total vessel area was calculated (πr²).

For intravitreal injections of EDN2, 500 μM of EDN2 peptide (Invitrogen, dissolved in 2 μl of sterile 1× PBS vehicle) was injected into the vitreous chamber using a Hamilton syringe with a 35 gauge needle. Vehicle alone was injected into other eyes as a control. After 4 weeks, optic nerves were assessed for optic nerve damage (see Analysis of glaucomatous damage).

Bosentan administration.

Bosentan, an endothelin receptor antagonist, was provided by Actelion Pharmaceuticals. Bosentan was incorporated into standard mouse chow (100 mg/kg, Test Diet). DBA/2J mice were administered bosentan from 6 months of age. Control animals were aged in the same environment on the same mouse chow (minus bosentan). Although bosentan is not reported to affect blood pressure, as it is an inhibitor of both the EDNRA and EDNRB endothelin receptors, some endothelin system inhibitors do affect blood pressure. To ensure that bosentan would not affect blood pressure in DBA/2J mice, a separate cohort of 8 mice was administered bosentan for 2 weeks, and then blood pressures was measured using previously described procedures (73). Blood pressures were not affected by bosentan administration (BP, mmHg ± SEM: 113.5 ± 0.67 treated; 112.7 ± 1.2 control, P = 0.8).

Clinical examination, IOP measurements, and RGC soma assessment.

DBA/2J mice develop a disease of the iris that leads to intraocular pressure elevation and glaucoma. Clinical examinations (30) and IOP measurements (74, 75) were performed as previously described. For clinical examinations, at least 20 eyes from each genotype or treatment group were examined at 6, 8, 10, and 12 months of age. RGC layer soma counts were preformed as described previously (15). Two-tailed Student t tests were performed to determine statistical significance. P values less than 0.05 were considered significant.