Genome-Wide Expression Profiling of Patients with Primary Open Angle Glaucoma

Abstract

Purpose.: To identify differentially expressed genes and to elucidate gene interaction networks and molecular pathways possibly contributing to the development of POAG.

Methods.: Genome-wide expression profiling experiments were carried out using ABI high-density oligonucleotide microarrays in leukocytes from 25 POAG patients and 12 age-, ethnicity-, and sex-matched normal controls. Significantly modulated genes were defined as those with a false discovery rate (FDR) <0.01 and an absolute fold change (FC) >1.5. These genes are then mapped to relevant biologic processes and pathways.

Results.: We identified 563 genes that were significantly dysregulated (410 upregulated and 153 downregulated) in POAG compared with normal controls (“POAG gene signature”). These genes were significantly enriched with functions related to, among others, nucleoside, nucleotide, and nucleic acid metabolism, the mitogen-activated protein kinase kinase kinase cascade, apoptosis, protein synthesis, cell cycle, intracellular signaling cascade, and nervous system development and function. Among the most significantly altered canonical pathways in POAG were the ephrin receptor signaling, ubiquitin proteasome pathway, hypoxia signaling, neuregulin, and G-protein coupled receptor signaling. Network analysis revealed potentially critical roles of UBE2, TBP, GNAQ, SUMO1, CREB, p70S6k, IFNG, and CaMKII that are interacting with NF-κB, ubiquitin, proteasome, PI3K/AKT, IL12, and PDGF in the disease pathogenesis.

Conclusions.: Our study revealed blood gene signatures that clearly distinguish POAG patients and normal controls, as well as altered pathways that may shed light on POAG pathogenesis.

Introduction

Glaucoma is one of the leading causes of blindness worldwide ¹ with a prevalence of over 2% in individuals older than 40 years. POAG is a complex, heterogeneous disease that is a major health concern throughout the world. It is estimated that more than 2.25 million Americans aged 40 years and older already have POAG. It is the most common type of glaucoma in Western countries and has risk factors that include elevated intraocular pressure (IOP) and age; but these factors do not predict the presence or degree of visual loss. ² POAG is characterized by the presence of glaucomatous optic neuropathy without an identifiable secondary cause. ^3,4 Abnormally elevated IOP is often associated with POAG and is a major risk factor for this disease. ⁴ A number of studies indicate the familial nature of POAG and support the presence of genetic factors in the pathogenesis of POAG. At least 14 linkage loci have been identified and are designated as GLC1A through GLC1N. ⁴ Several genes have been identified within these loci, including myocilin, optineurin, and WD repeat domain 36 (WDR36). ⁴ Genome Wide Association Studies (GWAS) identified common variants near the CAV1 and CAV2 genes, which are expressed in the trabecular meshwork and retinal ganglion cells that are involved in the pathogenesis of POAG. ⁵ Recently two loci for POAG at TMCO1 and CDKN2B-AS1 were identified using the GWAS technique. ⁶ Despite the great potential for GWAS approach in identifying disease-associated genetic variants in a wide range of human diseases including glaucoma, the technique explains only a limited amount of apparent heritability. ⁷ The human optic nerve is generally not available for study, but gene expression studies in human whole blood have proven valuable in understanding the genetics of a variety of neurologic and ophthalmologic diseases, ^8–13 including POAG. ¹⁴ This study investigated whole genome expression in blood in the hope of identifying specific genes or biological pathways that may contribute to POAG pathogenesis. Interpretation of the significance of gene expression changes in this circumstance is, of course, somewhat complex. Gene expression changes in blood may be present in the entire body (including bone marrow), but affect only the optic nerve. Alternatively, these gene expression changes might affect some other organ and secondarily affect the optic nerve, might affect the optic nerve together with other organs, or might be completely asymptomatic. Nevertheless, complex neurological diseases such as autism ¹⁵ and schizophrenia ¹⁶ have also been successfully studied using similar strategies. 

Materials and Methods

Patients and Control Subjects

Twenty-five Saudi patients were recruited who satisfied clinical criteria for POAG including: appropriate appearance of the optic disc and retinal nerve fiber layer (e.g., progressive thinning or notching of disc rim with nerve fiber layer defects); the presence of characteristic abnormalities in the visual field (e.g., arcuate scotoma, nasal step, or paracentral scotoma) in the absence of other causes of optic nerve damage; age greater than 40 years; open anterior chamber angles bilaterally on gonioscopy; and newly diagnosed POAG patients with no previous treatment or surgical procedure. Exclusion criteria included evidence of secondary glaucoma (e.g., pigmentary dispersion syndrome, narrow anterior chamber angles, or pseudoexfoliation), history of steroid use or ocular trauma or POAG patients already on medications or have had surgical procedure performed. Patients were recruited from the Glaucoma Clinic at the King Khalid Eye Specialist Hospital after signing an informed consent approved by the institutional review board. The mean (SD) age for patients was 53.2 (8.6), age range (41–71) years. Among patients, 15 were male and 10 were female. 

A second group of age- and sex-matched healthy Saudi Arab controls (n = 12) were recruited after a complete ophthalmologic examination proved them to be free from glaucoma. Selection criteria for those subjects were age >40 years, normal IOP, open angles on gonioscopy, and normal optic nerves on examination. All patients and controls were unrelated Saudi Arabs. The mean (SD) age for controls was 55.3 (8.4), age range (42–69) years. Among controls, six were male and six were female. 

All patients and controls were of Saudi ethnicity and this was established in three phases: POAG patients and controls with Saudi nationality as recorded in their medical files and on their medical treatment cards were selected; family and tribal name, which is an indication of the province they came from and helps determine relatedness to other patients; and talking to the patients and asking more questions about their family ancestry. Our group did extensive population genetics work on the Saudi tribes and our database, established during the course of the population-genetics work, can help us determine ethnicity and relatedness accurately and precisely. ^13,17  

Sequence Analysis of MYOC and OPTN

The coding exons, exon-intron boundaries, and promoter regions in the MYOC and OPTN genes were amplified by PCR from genomic DNA for all patients and control subjects and subjected to direct sequencing, as described previously. ¹⁸  

Applied Biosystems Expression Array Analysis

mRNA expression was analyzed using a microarray database (Applied Biosystems Human Genome Survey Microarray [ABI-HGSM] V2.0; Applied Biosystems, Foster, CA) with a microarray analyzer (Applied Biosystems 1700 Chemiluminescent Microarray Analyzer; Applied Biosystems). The ABI-HGSM V2.0 (P/N 4359030) contained 31,700 60-mer oligonucleotide probes representing 29,098 individual human genes. Digoxigenin-UTP labeled cRNA was generated and amplified from 1 μg of total RNA from each sample using a chemiluminescent RT-IVT labeling kit (P/N 4340472; Applied Biosystems) according to the manufacturer's protocol (P/N 4339629). Array hybridization was performed for 16 hours at 55°C. Chemiluminescence detection, image acquisition, and analysis were performed using a chemiluminescence detection kit (P/N 4342142; Applied Biosystems) and a microarray analyzer (P/N 4338036; Applied Biosystems) according to the manufacturer's protocol (P/N 4339629). 

Gene Expression Data Analysis

Images were auto-gridded and chemiluminescent signals were quantified, background subtracted, and spot- and spatially normalized using microarray analyzer software (P/N 4336391; Applied Biosystems). Detection thresholds were used following the manufacturer's recommendations for transcriptome analysis. Detection threshold was set as S/N >3 (a value indicating 99.9% confidence level that the signal is above background level) and quality flag <5000. 

The open source R/Bioconductor packages (Fred Hutchinson Cancer Research Center, Seattle, WA) were employed to analyze data via quantile normalization and to determine significant differences in gene expression levels between POAG patients and controls, ¹⁹ adjusting P values for multiple comparisons by the false discovery rate (FDR) according to the Benjamini-Hochberg procedure, ²⁰ as described previously. ²¹ T-tests were performed on the data to identify differentially expressed genes between POAG (n = 25) and control groups (n = 12) if the probe showed S/N ratio >3 in at least 50% of the samples in either group. Significantly modulated genes were defined as those with FDR <0.01 and absolute fold change (FC) >1.5. Unsupervised two-dimensional hierarchical clustering was performed using Euclidean distance with average linkage clustering. ²² Functional annotation and biological term enrichment analysis was performed by using the Database for Annotation, Visualization and Integrated Discovery (DAVID) Bioinformatics Resources, ²³ Expression Analysis Systematic Explorer (EASE), ²⁴ and the Protein Analysis Through Evolutionary Relationships (PANTHER) Classification System. For each molecular function, biological process, or pathway term, PANTHER calculated the number of genes identified in that category in both a list of differentially regulated genes and a reference list containing all probe sets present on the ABI-HGSM and compared these results using the binomial test to determine if more genes than expected were present in the differentially regulated list. ²⁵ Overrepresentation was defined by P < 0.05. 

Functional, Pathway, and Network Analyses

Functional, pathway, gene ontology, and network analyses were executed using PubGene (In the public domain, http://www.pubgene.org), the Kyoto Encyclopedia of Genes and Genomes (In the public domain, KEGG, http://www.genome.jp/kegg/genes.html), and pathway analysis software (Ingenuity Pathway Analysis [IPA] 6.3; Ingenuity Systems, Mountain View, CA). Gene identifiers for the differentially expressed genes in POAG were mapped to their corresponding gene object in pathway databases. These genes, called “focus genes,” were then used as a starting point for generating biological networks. A score consisting of the negative logarithm of the probability of the focus genes in a network being found together due to random chance was assigned to each network in the dataset to estimate the relevance of the network to the uploaded gene list. Scores ≥2 were considered significant using a 99% confidence level. Significances for biological functions or pathways in the signature genes for these functions or pathways were compared with the ABI-HGSM as a reference set. A right-tailed Fisher's exact test was used to calculate the probability that the biological function or pathway assigned to that data set was explained by chance alone. 

Results

Sequencing the MYOC and OPTN genes in the 25 POAG patients and in controls revealed the presence of one sequence variant in the MYOC gene, g.2259 G > T in exon 3, which resulted in codon change (p.G324V). This sequence variant was detected in six patients (24%) and in three controls (25%). After OPTN sequencing, gene sequence variants g.412 G > A in exon 4 and 469 G > C in exon 5 were each detected in one patient. The g.412 G > A variant was a previously reported polymorphism that did not change an amino acid (p.T34T). ¹⁸ The novel g.469 G > C sequence variant resulted in a codon change (p.Q53H), but was also found with a similar prevalence in control subjects, implying that it also was probably a nondisease-causing polymorphism in this population. 

The ABI-HGSM V2.0 contains 31,700 60-mer oligonucleotide probes representing 29,098 individual human genes. Overall, approximately 17,000 probes were detectable based on signal to noise >3 in >50% of the samples. Using FDR <0.01 and FC >1.5 cut-off criteria to identify the most differentially expressed “POAG gene signature,” 744 probes corresponding to 563 genes were identified as differentially expressed. Among these genes, 410 were upregulated and 153 were downregulated in POAG patients compared with controls (see Supplementary Material and Supplementary Table S1). Unsupervised two-dimensional hierarchical clustering as well as principle component analysis (PCA) clearly separated POAG patients from controls (Figs. 1A, 1B, respectively). 

Figure 1. 

 

Expression profiling analysis. (A) Unsupervised two-dimensional hierarchical clustering separated individuals into either POAG patients or controls. Highly expressed genes are indicated in red, intermediate in black, and weakly expressed in green. (B) Three dominant PCA components containing 74.4% of the variance in the data matrix identified individuals as either POAG patients or controls. (C) Canonical pathway. (D) Functional analysis of differentially expressed genes (up- or downregulated) in POAG. X-axis indicates the significance (−log P value) of the functional/pathway association that is dependent on the number of genes in a class as well as biologic relevance.

IPA analysis of significantly up-/downregulated genes in POAG revealed that ephrin receptor signaling, hypoxia signaling, neuregulin, and G-protein coupled receptor signaling are among the most significantly altered canonical pathways (Fig. 1C). Moreover, the POAG gene signature was significantly enriched with functions related to, among others, DNA replication, recombination and repair, protein synthesis, nervous system development and function, cell death, and cell cycle (Fig. 1D). 

Based on analysis using the PANTHER classification system, genes related to nucleoside, nucleotide, and nucleic acid metabolism, mitogen-activated protein kinase kinase kinase (MAPKKK) cascade, apoptotic processes, DNA metabolism, cell cycle, and intracellular signaling cascade were the most significantly overrepresented (P < 0.02; see Supplementary Material and Supplementary Table S2), results consistent with categories identified by IPA. The most significantly altered pathways included ubiquitin proteasome pathway (P = 3.15×10⁻⁴); insulin/IGF pathway-mitogen activated protein kinase kinase/MAP kinase cascade (P = 2.24× 10⁻³); epidermal growth factor (EGF) receptor (P = 2.88×10⁻³); and fibroblast growth factor (FGF) signaling pathways (P = 4.88×10⁻²; see Supplementary Material and Supplementary Table S3). These analyses shed new light into a large number of biological processes and pathways that may potentially be relevant to POAG. 

To elucidate how significantly dysregulated genes in POAG interact with each other as well as other genes in various pathways, the POAG gene signature was mapped to gene networks using the pathway analysis knowledgebase (Ingenuity Systems). These genes were mapped primarily to top networks related to, among others, post-transcriptional modification, protein synthesis, DNA replication, recombination, and repair, and nervous system development and function. Network analysis also revealed potentially important role of various genes, including UBE2, TBP, GNAQ, GNAO1, CREB, p70S6k, IFNG, and CaMKII, that are interacting directly or indirectly with NF-κB, ubiquitin, proteasome, PI3K/AKT, IL12 and PDGF that may potentially be implicated in POAG pathogenesis (Figs. 2A, 2B). 

Figure 2. 

 

Functional network analysis of POAG gene signature. (A, B) Top two scoring gene interaction networks with high relevancy scores (significance: score >25) for POAG signature genes. A score of 3 indicates that there is 1/1000 (score = −log P value) chance that the focus genes were assigned to a network randomly. Green indicates downregulated and red upregulated genes. Color intensity is correlated with fold change. Straight lines are for direct gene-to-gene interactions and dashed lines are for indirect ones.

Discussion

This study compared whole genome expression profiles in leukocytes from 25 patients with unequivocal POAG but no OPTN or MYOC gene mutations to the expression profile of 12 age-, sex-, and ethnicity-matched controls with no evidence of POAG on examination. A total of 563 genes were significantly dysregulated (410 upregulated and 153 downregulated) with a 1.5-fold expression difference and an adjusted P value <0.01 between POAG patients and controls. These genes were mapped to relevant biologic processes and pathways to gain greater insight into processes that might be activated or repressed and thereby to identify molecular functions and pathways potentially altered in POAG. 

Our data revealed various significantly modulated pathways in POAG patients, including ephrin receptor signaling, ubiquitin proteasome, and EGF receptor signaling. The Eph/ephrin signaling has been demonstrated to have potential role of in glaucoma, in modulation of axonal or glial responses after spinal cord and optic nerve injuries, and associated with retinal ganglion cells (RGC) axon pathology in animal models. ^25–28 Previous studies have also shown that ubiquitin proteasome pathway (UPP) is associated with many age-related dysfunctions of the eye. ^29,30 A recent study has shown that UPP is the major pathway for endogenous optineurin processing. ³¹ Optineurin is a gene linked to normal-tension glaucoma and adult-onset POAG. ³² Finally, the EGF receptor-signaling pathway has been linked to glaucoma through an indirect role related to inducible nitric oxide synthase contained in reactive astrocytes in the optic nerve head. This pathway may contribute to axonal damage if excessive nitric oxide is produced in response to increased IOP through EGF receptor tyrosine kinase, and pharmacological inhibition of inducible nitric oxide synthase may provide neuroprotection in the treatment of glaucoma. ³³  

The functional network and pathway analyses revealed that genes related to antigen presentation, immune response, development, and proliferation were significantly altered in POAG patients (see Supplementary Material and Supplementary Fig. S1). These results are consistent with previous studies suggesting the involvement of immune system regulation in cell fate decisions in glia and retinal ganglion cells that lead to glaucomatous optic nerve degeneration.³⁴ Furthermore, a recent study of global gene expression changes in experimental glaucoma found that the apoptotic, glutamate receptor, and G-coupled signaling pathways were differentially dysregulated in rat RGCs.³⁵ Similar observations were made in this study (see Supplementary Material and Supplementary Table S3, and Fig. 1), where genes related to these three pathways were differentially expressed. Intriguingly, we identified significantly dysregulated genes that are previously associated with visual perception and glaucoma, such as CRYAA, FBLN5, SMAD7, PAX2, MAP2K1, DCN, IGF1R, MC1R, WASL, and OPA1.^36–39 This concordance of results highlights the potential value of studying whole genome expression profile in blood, especially where the target tissue is not available or difficult to obtain. 

A number of genes involved in the hypoxia signaling pathway were also differentially expressed, including CREB3L4, SUMO1, UBE2D3, UBE2D4, UBE2G2, and UBE2S. Low levels of hypoxia have been suggested as a mechanism of damage in glaucoma, ³⁵ and hypoxia has been implicated in loss of RGCs either by apoptosis or necrosis. ⁴⁰ Additionally, free-radicals generated in hypoxic-ischemic conditions result in RGC loss because of an imbalance between antioxidant- and oxidant-generating systems. ⁴¹  

Evaluation of global gene expression patterns in these patients and controls has provided a molecular depiction of POAG and yielded insights into POAG pathogenesis. Nevertheless, there are a number of qualifications to this report. Optic nerve tissue was not available, and the gene-expression profile in these POAG patients was evaluated in blood, as has been done previously in a variety of neurologic and ophthalmologic diseases, ^8–11,13 including POAG. ³⁹ Transcriptional analysis of peripheral blood leucocytes is almost certainly not a perfect model system for studying POAG; however, peripheral blood cells inherit the same genetic information as RGCs, and leukocyte genomic profiling may well reflect pathologically important gene expression changes in RGCs and the optic nerve. 

The number of patients and controls studied is relatively small for a glaucoma report, although not for a gene expression study. ⁴⁰ Under any circumstance, the number of studied individuals was adequate to create a gene-expression profile of POAG patients that was significantly different from that of controls, and strict filtering criteria yielded a discrete number of over- and underexpressed genes involved in an identifiable spectrum of cellular functions. All of the studied individuals are from Saudi Arabia, but there is no reason to believe at this point that the genetic associations of POAG in Saudi Arabia are different from those of other ethnicities. 

There were too many differently expressed genes to carry out real-time PCR confirmation on each of them. Having said that, the microarray analyzer software (Applied Biosystems) employed here contain validated long oligo probes that were tested on numerous tissue samples before the array became commercially available. These longer probes (60-mer) have greater efficiency and specificity, and ABI results have been validated by a number of groups that found concordance between the results from the microarray analyzer (Applied Biosystems) and real-time PCR experiments. ^21,42–44  

This study supports a genetic component in the development of POAG. POAG patients had statistically significant generalized gene expression abnormalities that presumably target metabolic vulnerabilities of the optic nerve at least in part but clearly are not restricted to that structure. The genes and pathways affected are more diverse than expected, and considerable additional study will be necessary to determine how these genetic disturbances may interact to damage the optic nerve gradually over a period of decades. The presence of systemic genetic abnormalities also raises the possibility that novel types of systemic treatment might be effective in preventing or blunting the progression of the disease. 

Supplementary Materials

doc S1 

xlsx S2 

xlsx S3 

xlsx S4 

Footnotes

 Supported by King Faisal Specialist Hospital, Riyadh, Saudi Arabia.

Footnotes

 Disclosure: D. Colak, None; J. Morales, None; T.M. Bosley, None; A. Al-Bakheet, None; B. AlYounes, None; N. Kaya, None; K.K. Abu-Amero, None