Genome-wide analysis of multi-ancestry cohorts identifies new loci influencing intraocular pressure and susceptibility to glaucoma.

Elevated intraocular pressure (IOP) is an important risk factor in developing glaucoma, and variability in IOP might herald glaucomatous development or progression. We report the results of a genome-wide association study meta-analysis of 18 population cohorts from the International Glaucoma Genetics Consortium (IGGC), comprising 35,296 multi-ancestry participants for IOP. We confirm genetic association of known loci for IOP and primary open-angle glaucoma (POAG) and identify four new IOP-associated loci located on chromosome 3q25.31 within the FNDC3B gene (P = 4.19 × 10(-8) for rs6445055), two on chromosome 9 (P = 2.80 × 10(-11) for rs2472493 near ABCA1 and P = 6.39 × 10(-11) for rs8176693 within ABO) and one on chromosome 11p11.2 (best P = 1.04 × 10(-11) for rs747782). Separate meta-analyses of 4 independent POAG cohorts, totaling 4,284 cases and 95,560 controls, showed that 3 of these loci for IOP were also associated with POAG.

Primary open angle glaucoma (POAG) is the leading cause of irreversible blindness in the world[1]. The only modifiable risk factor for the development and progression of glaucoma is high intraocular pressure (IOP)[2], and lowering IOP is currently the only therapy that can reduce glaucomatous progression, even in forms of glaucoma that have IOP close to the statistical norm for the population (normal tension glaucoma or NTG)[34]. POAG and IOP are highly heritable; the lifetime risk of developing POAG is 22% among first degree relatives of patients[5], which is approximately 10 times higher than the rest of the population[1]. The IOP heritability is estimated to be approximately 55%[6]. Genetic studies have shown that the genetic risk of POAG and IOP are partly shared; polymorphisms within the TMCO1 gene are associated with both POAG risk [7] and IOP[8]. Studying genetic determinants of IOP is therefore likely to provide critical insights into the genetic architecture of POAG and open new avenues for therapeutic intervention.

In this study we present results from a meta-analysis of genome-wide association studies (GWAS) of IOP from 18 studies participating in the International Glaucoma Genetics Consortium (IGGC) and an assessment of the importance of the genetic findings for susceptibility to POAG (Figure 1). The IOP meta-analysis included 35,296 subjects (7,738 Asians and 27,558 of European descent) drawn from the general populations of seven countries. Demographic characteristics of these population cohorts are given in the Supplementary Table 1. Genotyping assays and imputation to HapMap2 haplotypes were performed at individual sites. Association analyses were performed using an additive model, with IOP as outcome, number of alleles at each polymorphic site as predictors, adjusting for age and sex. IOP levels for participants who were receiving IOP-lowering therapy at the time of the study and whose baseline, pre-treatment levels were not available were imputed as previously described[8]. Subjects who had undergone surgery or had other eye diseases that could affect IOP were removed (Supplementary Note). Secondary analyses were carried out adjusting for central corneal thickness (CCT), which is known to influence IOP measurements [9].

Figure 1

Flowchart of the analyses

Results from a meta-analysis of IOP in participants from 18 general population cohorts were validated in four clinical case control cohorts and checked for transcription regulation activity in three tissues from 856 British subjects.

After applying conventional quality control filters, a fixed effects meta-analysis of 22 autosomes across the cohorts was performed with approximately 2.5 million markers. Within-study genomic inflation factors[10] ranged between 0.992 and 1.043 (Supplementary Table 2 and Supplementary Figure 1), indicating a lack of major population stratification bias within each study. SNPs available in fewer than 16 cohorts or showing large heterogeneity (defined as I2>75% [11]) were removed. We found 145 SNPs (Supplementary Table 3) whose association crossed the conventional genome-wide significance threshold (p<5×10−08) [12]. All of these SNPs clustered around seven separate regions of the genome (Figure 2, Supplementary Figures 2 and 3). Two of the regions associated with IOP in our meta-analysis had previously been implicated in IOP variability: the region near the TMCO1 locus[7,8] (p=2.19×10−09 for rs7555523), and near GAS7 gene[8] (p=1.03×10−11 for rs9913911). A third locus, novel for IOP, was near the CAV1 and CAV2 genes (p=1.87×10−11 for rs10258482) which had previously been associated with POAG[13].

Figure 2

Manhattan plots of the results from the meta-analyses of results from 18 multi-ethnic cohorts from the IGGC

The 22 autosomes are plotted along the x-axis, whereas the values in the y-axis denote the −log10 -transformed p-values from the meta-analysis of association with IOP observed for any of the SNPs. Loci previously associated with IOP or glaucoma have been grayed out.

Novel associations were identified within large linkage disequilibrium (LD) block on chromosome 11 encompassing, among other genes, AGBL2, SPI1 and PTPRJ (best p=1.04×10−11 for rs747782), (Supplementary Figure 1). Two additional loci were mapped on chromosome 9: one at 9q31.1 upstream the ABCA1 gene (p=2.80×10−11 for rs2472493) and the other at 9q34.2 within the ABO blood group gene (p=3.08×10−11 for rs8176743). A fourth region was detected on chr3q25.31, within the FNDC3B gene (p=4.19×10−08 for rs6445055).

Interestingly, of all the loci previously associated with glaucoma or related quantitative traits[14], CDKN2BAS and SIX1/SIX6 are not associated with IOP in the meta-analysis. It is possible these two loci exert their influence on POAG through mechanisms unrelated to intraocular pressure.

Genome-wide significant SNPs from the IOP meta-analysis were then investigated for their effect on the clinical outcome of POAG in four independent cohorts representing a combined 4,284 POAG cases (normal and high tension glaucoma) and 95,560 controls (details about these cohorts in the Supplementary Note). Associations to POAG were found for the novel regions near ABCA1 gene (p=4.15×10−09 for rs2472493), near the FNDC3B gene (p=0.03 for rs6445055) and on the chromosome 11 cluster (p=0.008 for rs12419342). We did not find significant statistical evidence of association of POAG to the ABO locus. The case-control analyses reinforced association evidence at the previously identified loci on TMCO1 (p=1.34×10−16 for rs7555523), CAV1/CAV2 (p=6.27×10−09 for rs10258482) and GAS7 (p=5.22×10−13 for rs12150284). All alleles associated with increased IOP levels also raised glaucoma risk (Table 1).

Table 1

Results for the association with Intraocular Pressure from the general population cohorts for SNPs significant at a multiple testing correction level (p<5×10−08) and their association with POAG in a case-control validation meta-analyses *

SNP					Association with IOP in discovery cohort					Association in POAG Case-Controls
Chr.	Position	id	A1/A2	Nearest Gene	Beta	SE	p-value	Heterogeneityp-value	I2	OR	95% CIs		p value
1	165687205	rs4656461	G/A	TMCO1	0.228	0.039	6.51×10−09	0.46	0.00	1.38	1.28	1.5	2.55×10−15
1	165718979	rs7555523	C/A	TMCO1	0.235	0.039	2.19×10−09	0.55	0.00	1.4	1.3	1.52	1.34×10−16
3	171992387	rs6445055	A/G	FNDC3B	−0.177	0.03	4.19×10−08	0.17	0.24	0.92	0.85	0.99	0.03
7	116150095	rs10258482	A/C	CAV1	0.196	0.029	1.87×10−11	0.81	0.00	1.2	1.13	1.28	6.27×10−09
7	116150952	rs10262524	A/C	CAV1	0.186	0.029	9.69×10−11	0.67	0.00	1.2	1.13	1.28	1.39×10−08
9	107695848	rs2472493	G/A	ABCA1	0.159	0.024	2.80×10−11	4×10−05	0.66	1.24	1.16	1.34	4.15×10−09
9	136131415	rs8176743	T/C	ABO	0.261	0.039	3.08×10−11	0.53	0.00	1.07	0.96	1.19	0.2
11	47468545	rs12419342	C/T	RAPSN	0.153	0.026	4.77×10−09	0.75	0.00	1.09	1.02	1.16	0.008
11	47940925	rs747782	C/T	NUP160,PTPRJ	0.203	0.03	1.04×10−11	0.95	0.00	1.03	0.96	1.11	0.36
11	47969152	rs1681630	T/C	PTPRJ	0.144	0.026	1.69×10−08	0.6	0.00	1.06	0.99	1.12	0.08
11	48004369	rs7946766	T/C	PTPRJ	0.23	0.035	2.71×10−11	0.35	0.09	1.03	0.95	1.12	0.43
17	10031183	rs9913911	G/A	GAS7	−0.179	0.026	1.03×10−11	4×10−04	0.61	0.8	0.75	0.85	2.98×10−13

Abbreviations used: Chr. - chromosome, A1/A2 - reference /alternative alleles, beta- linear regression coefficient (mmHg), SE-standard error of the regression coefficient, OR- Odds Ratios, 95% CIs – 95% Confidence interval for OR.

We then examined whether the effect sizes of SNPs on IOP levels (BetaIOP) were linearly related to their effect sizes on POAG (BetaPOAG) using a causal inference framework as previously described[15]. In a linear regression analysis, we observed a significant association between BetaIOP and BetaPOAG (p=0.03, Supplementary Table 4), suggesting the strength of SNPs’ effect on IOP levels is correlated with their effects on the risk of POAG.

We subsequently investigated the relationship between variants within the seven regions associated with IOP and cis- regulation of mRNA expression in three tissues (adipose, lymphoblastoid cell lines [LCLs] and skin) from a sample of 856 British subjects[16]. The most significant eQTL associations were generally observed in LCLs for most loci, except for CAV1 where effects were strongest in adipose and skin tissues (Table 2). Significant eQTL association was observed for rs4656461 and rs7555523 (p=0.003 and p=0.0001 with TMCO1 and ALDH9A1 transcript expressions in skin and LCLs respectively), rs2024211 (p=5.43×10−16 and 3.84×10−13 with CAV1 in adipose and skin tissues respectively), rs2472493 with ABCA1 (p=3.67×10−05 in LCLs) and rs1681630 with the SPI1 expression on chromosome 11 (p=2.72×10−10 in LCLs) among others (Table 2, Supplementary Table 5A). These SNPs also had the strongest eQTL effects for their respective transcripts (Supplementary Table 5B).

Table 2

Summary of eQTL effects observed in three different tissues extracted from 849 individuals for SNPs associated with IOP*

The SNPs listed are the same as those in Table 1.

SNP					eQTL Effects p-values
Chr	Position	id	A1/A2	Nearest Gene	Adipose	LCL	Skin	Probe ID	Gene
1	165687205	rs4656461	G/A	TMCO1	0.004	0.12	0.003	ILMN_1793829	TMCO1
1	165718979	rs7555523	C/A	TMCO1	0.39	0.0001	0.05	ILMN_1761804	ALDH9A1
3	171992387	rs6445055	A/G	FNDC3B	N.S.	N.S.	N.S.	-	-
7	116150095	rs10258482	A/C	CAV1	N.S.	N.S.	N.S.	-	-
7	116150952	rs10262524	A/C	CAV1	5.79×10−16	8.54×10−05	3.91×10−13	ILMN_1687583	CAV1
9	107695848	rs2472493	G/A	ABCA1	0.19	3.67×10−05	0.36	ILMN_1766054	ABCA1
9	136131415	rs8176743	T/C	ABO	N.S.	N.S.	N.S.	-	-
11	47468545	rs12419342	C/T	RAPSN	0.002	4.32×10−08	0.0003	ILMN_1696463	SPI1
11	47940925	rs747782	C/T	NUP160,PTPRJ	N.S.	N.S.	N.S.	-	-
11	47969152	rs1681630	T/C	PTPRJ	0.006	2.72×10−10	0.002	ILMN_1696463	SPI1
11	48004369	rs7946766	T/C	PTPRJ	0.66	2.02×10−05	0.0066	ILMN_1688627	AGBL2
17	10031183	rs9913911	G/A	GAS7	N.S.	N.S.	N.S.	-	-

Abbreviations used: LCL- lymphoblastoid cell lines, N.S. – no significant association detected

We measured the mRNA expression levels of the identified genes in adult ocular tissues using reverse-transcriptase PCR. We found that most of the identified genes, including TMCO1, FDNC3B, CAV1/CAV2, ABCA1 and GAS7 were expressed in most ocular tissues (Supplementary Table 6). The expression of genes within the chromosome 11 locus showed varied expression levels across ocular tissues.

Gene-based tests or Gene Ontology enrichment analyses failed to identify any new genes or pathways after multiple testing correction (Supplementary Tables 7 and 8).

Altogether, these SNPs explain approximately 1.2% of the IOP heritability in the TwinsUK cohort[17], 1.5% of the IOP phenotypic variability in the Rotterdam study[18] and between 0.6 and 1.2% in Asians. FNDC3B has been associated with central corneal thickness (CCT) [19] and as CCT has a significant effect on IOP measurements[20], we performed an additional meta-analysis of IOP adjusted for age, sex and CCT in a smaller subsample which had CCT measures (19,563 subjects from 13 population cohorts). The association for rs6445055 remained nominally significant, albeit reduced (p= 9.87×10−04, β= −0.121 compared to −0.177 prior to adjustment for CCT). This suggests that this locus has at least some CCT-independent effect over the IOP levels. The association evidence remained consistent, although slightly reduced, for the other loci (Supplementary Table 9).

We report association of variants within the ABCA1 gene and IOP and POAG. A strong eQTL effect was observed in LCLs (p=3.67×10−05) for the most highly associated SNP (rs2472493) in our analyses. The ABCA1 gene is expressed in many tissues[21] and its expression in leukocytes is significantly up-regulated in glaucoma patients[22].

Associations for a number of SNPs within the ABO blood group gene and IOP, although statistically significant and homogeneous across the participating cohorts, were not observed in the glaucoma case-control meta-analysis. This might be due to type I error in the initial meta-analysis or insufficient power to detect a primarily IOP-led effect in cases that included NTG patients, resulting in a type II error in the latter. Four of the nine GWAS polymorphisms associated at genome-wide significance in the ABO locus are non-synonymous variants, determining B blood group[23]. This may be relevant given previous observations that the B blood group is epidemiologically associated with glaucoma, including POAG[24], although the mechanisms remain unclear.

Association was found between IOP and variants lying over a large region on chromosome 11. Of the many genes in that region, eQTL analyses singled out SPI1 and AGBL2 as possible candidates for prioritization in future studies. eQTL analyses also raised the possibility of ALDH9A1 as a candidate for IOP regulation, given its strong expression in ciliary body[25] and location just downstream of the TMCO1-associated variant. The eQTL results also suggest that CAV1 is a stronger candidate than CAV2, although transcription regulation may not be the only mechanism influencing IOP at this locus.

Although IOP and POAG are strongly genetically correlated[26], we explored further their shared genetic backgrounds. Using independent (i.e. not in LD) SNPs with p< 10−06 in the IOP GWAS meta-analysis as described elsewhere[27], we found a statistically significant polygenic overlap between IOP and POAG in the ANZRAG cohort (p=4.33×10−05). The variance explained in POAG was 0.7%, which changed little if less significant SNPs were progressively included in the model (Supplementary Table 10).

There are potential limitations to this study. First, there is variability across the studies in terms of IOP measurement methods, although the differences are likely to be small[28]. In addition, we maximized power to discover genetic variants of small effect size by including multiethnic cohorts, at the risk of introducing heterogeneity into the study. Heterogeneity was however generally low (Table 1) for most of the loci reported, so we consider our results to be conservative. Second, assessment of clinical importance using panels of POAG cases is not equivalent to a formal replication. Even in this case, we expect our results to be over-conservative at the price of reduced sensitivity, which could be a possible reason for non-validation of our association with IOP in the ABO blood group locus. Finally, we based our eQTL analysis on sample tissue availability rather than ideal ocular tissue types. Tissues such as trabecular meshwork would have been preferable, but they are impractical because they are generally less accessible. We tried to circumvent this limitation by studying three different tissues, but caution is required when interpreting eQTL results.

Despite these considerations, our report of seven loci associated with IOP and glaucoma, of which four are novel, is an important step forward to better understanding the mechanisms of IOP regulation-currently the only modifiable risk factor for POAG.

Online Methods

IGGC participants

All studies participating in this meta-analysis are part of the International Glaucoma Genetics Consortium. The discovery cohorts included 27,558 individuals of European ancestry from 14 studies (ALIENOR, BATS, BMES,[29,30] ERF,[31,32] Framingham Family Study,[33] GHS1, GHS2, ORCADES,[34] RAINE, [35-37] RS-I, RS-II, RS-III,[38] TEST[39] and TwinsUK[40]). In addition, 7,738 individuals of Asian ancestry from four cohorts (BES,[41] SCES,[42] SiMES,[43] SINDI,[42]) were included. In addition four case-control population panels were used, all of European ancestry: ANZRAG,[7] MEEI, NEIGHBOR and deCode. General methods, demographics and phenotyping of the study cohorts have previously been described extensively and are provided in the Supplementary Note and Supplementary Table 1. All studies were performed with the approval of their local Medical Ethics Committees and written informed consent was obtained from all participants in accordance with the Declaration of Helsinki.

Phenotype measurements

Eligible participants underwent an ophthalmologic examination including measurements of IOP and, for most but not all studies, measurements of central corneal thickness (CCT). Each participating cohort was phenotyped separately and IOP measurement methods used by each are described in the Supplementary Table 1.

Genotyping & imputation

The study samples were genotyped on either the Illumina (San Diego, CA, USA) or Affymetrix (Santa Clara, CA, USA) platforms. Each study performed single nucleotide polymorphism (SNP) imputation using the genotype data, together with the HapMap Phase II ethnically matched reference panels (CEU, JPT+CHB, or the 4 HapMap populations) on the basis of build 36 databases (release 22 or 24). The Markov Chain Haplotyping software, IMPUTE[44,45] or MACH,[46] were adopted for imputation. A detailed description regarding genotyping platforms and imputation procedures for each study is provided on Supplementary Table 1.

Stringent quality control of genotype data was applied in each cohort. Samples with low call rates (<95%) or with gender discrepancies were excluded. Cryptically related samples and outliers in population structure from principal component analyses were also excluded. SNPs flagged with missingness >5%, gross departure from Hardy Weinberg Equilibrium (P value <10−6) and minor allele frequency (MAF) <1% were removed from further analyses.

Statistical Analysis

For each study, an allele-dosage regression model at each directly genotyped or imputed SNP was conducted to determine its association with IOP. Eyes with prior glaucoma surgery or laser were excluded. For subject receiving IOP-lowering medication, we added 25% to the measured IOP levels to estimate pre-treatment IOP, based on a reported average of 17% to 33% IOP reduction caused by IOP lowering medication in a meta-analysis of clinical trails.[47] The mean of the right and left IOP measurements was used. When data from only one eye were available, the IOP measurement from the available eye was used.

For the analyses, we assumed an additive genetic model where the dosage of each SNP is a continuous variable ranging from 0 to 2 for minor alleles carried. Primary analysis for IOP was adjusted for age and sex. Additional adjustment for principal components was carried out by a few participating cohorts to correct for subtle population substructure.

The per-SNP meta-analyses were performed using the GWAMA software with weighted inverse-variance approach, assuming fixed effects, as for initial discovery purposes the fixed-effects model is preferred for increased statistical power.[48] A Cochran’s Q test and I2 were used to assess heterogeneity across studies[49] . For each participating cohort, only SNPs with sufficient imputation quality scores (proper-info of IMPUTE or R2 of MACH >0.3) were included into the meta-analysis.

Gene-based testing was conducted using VEGAS software[50] on the European ancestry and Asian ancestry meta-analysis results separately. VEGAS incorporates information from the full set of markers within a gene and accounts for LD between markers by using simulations from the multivariate normal distribution. For samples of European descent, we used the HapMap 2 CEU population as the reference to estimate patterns of LD. For Asian ancestry groups, we used the combined HapMap 2 JPT and CHB populations as the reference population to approximate LD patterns. To include gene regulatory regions, SNPs were included if they fell within 50 kb of a gene. We performed meta-analysis on the two sets of gene-based P-values using Fisher’s method.

VEGAS-Pathway analysis[50,51] was carried out using prespecified pathways from Gene Ontology. Pathways of with 10 to 1,000 components were selected, yielding 4,628 pathways. Pathway analysis was based on combining gene-based test results from VEGAS. Pathway P-values were computed by summing χ2 test statistics derived from VEGAS P-values. Empirical VEGAS-Pathway P values for each pathway were computed by comparing the summed χ2 test statistics from real data with those generated in 500,000 simulations where the relevant number (according to the size of the pathway) of randomly drawn χ2 test statistics was summed. To ensure that clusters of genes did not adversely affect results, within each pathway, gene sets were pruned such that each gene was >500 kb away from all other genes in the pathway. Where required, all but one of the clustered genes was dropped at random when genes were clustered. We performed meta-analysis on the two sets of pathway p-values using Fisher’s method.

To investigate shared genetic background by using a large number of autosomal SNPs we made a systematic evaluation of the overlap between the IOP and POAG based on profile scores following approaches previously described[27]. We estimated the relative risk for each SNP of interest based on a discovery set (IOP), with a profile score computed for every individual in a target set of interest (POAG). For each target set individual the profile score was computed as the number of risk alleles weighted by the effect size estimated in the discovery set. The discovery set comprised the European ancestry derived samples from our meta-analysis and the target set was was a set of 590 glaucoma cases and 3956 controls, as previously described[7]. To ensure that there was not a high degree of dependence between the SNPs included in the profile score, we filtered the set of SNPs used in the profile score so that only a set of 149,571 SNPs in low linkage disequilibrium (r2<0.5) were used. We constructed models including more SNPs progressively lowering the threshold of inclusion (i.e. p<0.000001, p<0.00001, p<0.0001, p<0.001, p<0.01, p<0.1, p<0.5). Profiles derived from IOP SNP effects were tested for association with the phenotype (POAG here) using a logistic regression. Variance explained was assessed using Nagelkerke’s pseudoR[2] measure[52].

To assess whether and to what degree IOP levels confer POAG risk, we performed a causal inference analysis using an instrumental variable framework as previously described[15]. In brief, we obtained estimates of effect size (BetaIOP) for the association of a given SNP with IOP from the meta-analysis of the 18 discovery cohorts. For the association of a given SNP with POAG, we obtained estimates of the effect size (BetaPOAG) from the four case-control panels as described above. We selected the the strongest associated SNP from each of the genome-wide significant IOP loci that we identified. To assess whether the strength of SNPs’ association with IOP predicts the risk of POAG, we conducted linear regression analysis using the effect sizes of each SNP for IOP (BetaIOP) as independent variables and the effect sizes of POAG (BetaPOAG) as dependent outcome variables . A total of eight independent IOP-associated SNPs were used for this analysis, including rs7555523 (TMCO1), rs6445055 (FNDC3B), rs10258482 (CAV1), rs2472493 (ABCA1), rs8176743 (ABO), rs747782 (NUP160-PTPRJ), and rs9913911 (GAS7).

Gene expression in human tissues

Adult ocular samples were obtained from normal eyes of an 82-year-old European ancestry female from the North Carolina Eye Bank, Winston-Salem, North Carolina, USA. All adult ocular samples were stored in Qiagen’s RNAlater within 6.5 hours of collection and shipped on dry ice overnight to the lab. Isolated tissues were snap-frozen and stored at −280 °C until RNA extraction. RNA was extracted from each tissue sample independently using the Ambion mirVana total RNA extraction kit. The tissue samples were homogenized in Ambion lysis buffer using an Omni Bead Ruptor Tissue Homogenizer per protocol. Reverse transcription reactions were performed with Invitrogen SuperScript III First-Strand Synthesis kit. The expression of the identified genes were assessed by running 10 ul reactions with Qiagen’s PCR products consisting of 1.26 ul H2O, 1.0 ul 10X buffer, 1.0 ul dNTPs, 0.3 ul MgCl, 2.0 ul Q-Solution, 0.06 ul taq polymerase, 1.0 ul forward primer, 1.0 ul reverse primer and 1.5.0 ul cDNA. The reactions were run on a Eppendorf Mastercycler Pro S thermocycler with touchdown PCR ramping down 1°C per cycle from 72 °C to 55 °C followed by 50 cycles of 94 °C for 30 seconds, 55 °C for 30 seconds and 72 °C for 30 seconds with a final elongation of 7 minutes at 72 °C. All primer sets were designed using Primer3.[53] Products were run on a 2% agarose gel at 70 volts for 35 minutes. Primer sets were run on a custom tissue panel including Clontech’s Human MTC Panel I, Fetal MTC Panel I and an ocular tissue panel.