Genome-wide association study of intraocular pressure identifies the GLCCI1/ICA1 region as a glaucoma susceptibility locus

To discover quantitative trait loci for intraocular pressure, a major risk factor for glaucoma and the only modifiable one, we performed a genome-wide association study on a discovery cohort of 2175 individuals from Sydney, Australia. We found a novel association between intraocular pressure and a common variant at 7p21 near to GLCCI1 and ICA1. The findings in this region were confirmed through two UK replication cohorts totalling 4866 individuals (rs59072263, Pcombined = 1.10 × 10−8). A copy of the G allele at this SNP is associated with an increase in mean IOP of 0.45 mmHg (95%CI = 0.30–0.61 mmHg). These results lend support to the implication of vesicle trafficking and glucocorticoid inducibility pathways in the determination of intraocular pressure and in the pathogenesis of primary open-angle glaucoma.


INTRODUCTION
Elevated intraocular pressure (IOP) is a major risk factor for the development and progression of glaucoma, which is the commonest cause of irreversible blindness worldwide (1). Globally, 60 million people are affected of whom 8.4 million are blind in both eyes (2). Before this end stage there is a greater risk of falls (3), fractures (4) and motor vehicle collisions (5).
The link between high IOP and glaucoma is clearly seen when there is a readily identifiable cause for raised IOP such as trauma, inflammation, dysgenesis or obstruction of the aqueous humour drainage system in the trabecular meshwork of the iridocorneal angle by iris tissue, pigment or other deposits (2). It is also evident even when there is no clear cause or obstruction of the iridocorneal angle. Glaucoma occurring under these latter conditions is known as primary open-angle glaucoma (POAG) and is the most prevalent form of glaucoma, causing the greatest burden of disease (6). The risk of developing POAG is approximately five times higher in subjects with IOPs above the population 95th centile than in subjects with lower IOPs (7); the higher the IOP at screening, the greater the risk of POAG (8). Both mean and maximum IOP have been reported as closely associated with visual field deterioration (8 -15) and optic nerve damage (16). In a group of subjects affected by POAG despite untreated IOPs below the population 95th centile ('normal pressure' glaucoma), worse visual field damage was found in the eye with the higher IOP (17). Furthermore, IOP is currently the only modifiable risk factor for POAG. Therapeutic lowering of intraocular pressure reduces the risk of developing glaucoma (18) and retards disease progression (19)(20)(21). In addition to inherent interest, and its role as a modifiable trait on the causal path to disease, the known heterogeneity of POAG means that a focus on the genetic basis of IOP may be more powerful than studies of POAG directly (22,23).
Heritability estimates for adult IOP range from 0.29 to 0.62 (24)(25)(26)(27)(28). To date, two regions have been reported to be robustly associated with IOP as a result of linkage analyses, one at 10q22 from a Tasmanian glaucoma pedigree (29) and one at 5q22 in a West African cohort (30). A recent Dutch genome-wide association study (GWAS) for IOP (31) also found two associated loci, one of which (TMCO1) overlaps with a previously published POAG GWAS locus (32).
Further elucidation of the genetic basis of IOP would aid in understanding the biology of a major blinding disease (POAG), in risk stratification for development and progression of the disease and in providing therapeutic targets. In order to identify genomic variants that contribute to the determination of IOP, we performed a GWAS as a component of the Wellcome Trust Case Control Consortium 2 (WTCCC2) (http://www.wtccc.org. uk/ccc2/).

RESULTS
A total of 2765 discovery samples from the Blue Mountains Eye Study (BMES) were genotyped on the Illumina 660W-Quad array at the Wellcome Trust Sanger Institute (WTSI). Genotype imputation was performed in all 2302 individuals passing genotype quality control (QC) using IMPUTE2 (33)  For each individual, the mean IOP between their two eyes was calculated. If an individual had undergone surgery or treatment to an eye, data from that eye were excluded. If the difference in IOP between an individual's two eyes was .10 mmHg, the individual was excluded. Eight outlying individuals whose measurements did not conform to a normal distribution were also removed (see Materials and Methods). After genotype and phenotype QC, a total of 2175 samples were available for analysis.
A linear regression fitting an additive model was performed in SNPTEST (https://mathgen.stats.ox.ac.uk/genetics_software/ snptest/snptest.html), where genotype, accounting for imputation uncertainty, was set as the explanatory variable and mean IOP as the response variable ( Fig. 1). Age and sex were included as continuous and categorical covariates, respectively. This analysis will be referred to as the primary scan. The genomic inflation statistic (l) was 1.01 (Supplementary Material, Fig. S1), suggesting that population structure was not a major problem in our discovery analyses.
Replication of the primary scan was attempted at 31 SNPs covering 17 regions (Supplementary Material online, Table S1 and see Materials and Methods), of which two were regions previously found to be associated with IOP (31). Two platforms were used for the first phase of replication; genotypes for four SNPs were taken from the custom built Immunochip (36) (see Materials and Methods), and genotypes for the other 27 SNPs were obtained from a Sequenom plex. Results from discovery and replication were combined using a fixed effects meta-analysis by averaging the estimated effect size parameter across the data sets, weighting by the inverse of the variance in the estimates.
For the first phase of replication, we used individuals from the European Prospective Investigation into Cancer (EPIC) Norfolk cohort (see Materials and Methods). Following QC, data were available for 2833 individuals genotyped on Sequenom and 2493 individuals genotyped on Immunochip (all with age, sex and mean IOP data) with 2461 individuals overlapping between the Sequenom and Immunochip replication cohorts.
The results from our data at SNPs previously found to be associated with glaucoma (32,(37)(38)(39) or IOP (31) are shown in Table 1. Notably, the discovery data in our study contributed to the replication phase in the IOP GWAS (31). An SNP in TMCO1 (rs7555523), which has been associated with IOP (31) showed evidence for association with mean IOP in both our discovery (P ¼ 0.0303) and replication data (P ¼ 1.70 × 10 23 ). At the other region previously associated with IOP (31), on chromosome 17p13 (rs11656696), there was no evidence for association in our discovery data (P ¼ 0.548), whereas there was evidence for association in our replication data (P ¼ 0.016). One of the SNPs previously associated with POAG (39) was also associated with mean IOP in our data, (rs4236601, P ¼ 1.46 × 10 23 ). The remaining previously associated SNPs were not associated (P , 0.1) with IOP in our discovery data.
Following first phase replication, one of the newly identified regions reached a P combined ,5 × 10 28 . This was on chromosome 7p21 and the top SNP (rs59072263, P combined ¼ 2.32 × 10 28 ) is between the genes ICA1 and GLCCI1 ( Fig. 2 and Supplementary Material, Fig. S2). Genotypes at this SNP were imputed in the discovery data and directly genotyped in the replication. To check the reliability of the imputed calls, a subset of 214 discovery samples were directly genotyped at rs59072263. The concordance between the imputed and genotyped calls was 0.96 (eight out of 214 discordant calls) and the correlation was 0.87, suggesting the imputation at this locus is acceptable.
No evidence for replication at P , 0.05 was found at the other 14 regions (Supplementary Material, Table S1). Two SNPs at 11p12 (rs10430986 and rs12222492) showed suggestive evidence for association in the same direction as the discovery data (P ¼ 0.07, P ¼ 0.11), although data were not available for the top SNP in the region . This locus may be worth additional investigation.
A second round of replication of the 7p21 association was next attempted in 2033 individuals from TwinsUK (40), in which the genotypes at rs59072263 were imputed. The SNP had a onesided P-value of 0.053 assuming the same risk allele as the discovery, giving combined evidence across the three cohorts of P ¼ 1.10 × 10 28 (Supplementary Material, Fig. S3), where the G allele is associated with an IOP increase of 0.45 mmHg (95% CI ¼ 0.30 -0.61 mmHg).
Replication of the signal of association at 7p21 was also assessed in the large Dutch meta-analysis (31) of IOP comprising 9680 individuals of European ancestry. The associated SNP rs59072263 was not genotyped in this study, but the best available tag (rs6959703; genotype correlation, r 2 ¼ 0.61) showed a one-sided P-value ¼ 0.017. Imputation of genotypes by the Rotterdam group at rs59072263, using MACH (http://www.sp h.umich.edu/csg/abecasis/MACH/index.html) and the 1000 Genomes data, was also attempted but showed weaker evidence of association (one-sided P-value ¼ 0.178), However, unlike the discovery data and TwinsUK data, which had imputation accuracy ('IMPUTE info' and 'MACH r 2 ', respectively) of 0.95 and 0.89 at rs59072263, the largest component of the Dutch meta-analysis, comprising almost two-thirds of the study, had MACH r 2 , 0.75. Therefore, particularly in light of the clear replication at the best genotyped tag SNP, the lack of replication in the imputed data at P , 0.05 may be due in part to the loss of statistical power which results from the increased noise (imputation uncertainty) in the inferred genotypes (41). We note that the risk allele is the same in all three replication cohorts.
We also undertook a pathway analysis of SNPs with P , 10 24 (see Materials and Methods), but this failed to pinpoint any pathways with a compelling biological connection to IOP biology.

DISCUSSION
We have thus identified a novel association between a region on chromosome 7p21 and mean IOP. Since IOP is on the causal pathway to glaucoma, this association is also relevant for susceptibility to glaucoma. The most associated SNP in this region is flanked by two genes, ICA1 and GLCCI1, both of which have been shown to be expressed in the human eye (information from the EMBL-EBI website, see http://www.ebi.ac.uk/gxa/). ICA1 encodes a protein involved in the regulation of secretory vesicle trafficking (42) and has been implicated as an auto-antigen in insulin-dependent diabetes mellitus (43) and primary Sjögren's syndrome (44). Vesicular metabolism pathways have previously been ascribed a role in the pathobiology of POAG: a casecontrol GWAS of POAG found associations with the CAV1 and CAV2 genes which control vesicle transport in transcytosis mechanisms (39). The latter mediate drainage of aqueous humour from the eye and it is recognized that increased IOP '-' data not available. results from a failure of drainage rather than an overproduction of aqueous humour. The GLCCI1 gene has been ascribed a role in the sensitivity of various tissues to glucocorticoids, initially thymoma cell lines (45) and recently in the degree of response to glucocorticoid therapy in asthma (46). The variant associated with glucocorticoid therapy in asthma (46) was not associated with IOP in our discovery (P ¼ 0.8), and was not in strong linkage disequilibrium (LD) with the IOP associated variant we identified (r 2 ¼ 0.114 in the 1000 Genomes data). Topical glucocorticoid therapy causes raised IOP in susceptible individuals (47) and the likelihood of such an IOP response is greater in POAG patients, POAG suspects and first-degree relatives of POAG patients (48). The BMES has previously reported an association between use of inhaled corticosteroids and a finding of elevated IOP or glaucoma in subjects with a glaucoma family history (49). Furthermore, patients with POAG have elevated blood (50 -54) and aqueous humour (52) levels of the endogenous glucocorticoid cortisol compared with age-matched normal subjects. The first glaucoma gene identified, MYOC (55), was a glaucoma candidate gene because of its glucocorticoid inducibility (56). It is thus plausible that GLCCI1 may influence IOP via the response to endogenous cortisol.
The effect size of the novel associated SNP we report (Table 2) is substantially higher than that of previously reported SNPs for IOP (beta ¼ 0.45 for rs59072263 when compared with 0.19 for rs11656696 and 0.28 for rs7555523). We can get some indication of the possible consequences of the SNP on POAG, based on its measured effect for IOP, from population studies which relate changes in IOP to POAG risk. This approach suggests that each copy of the G allele at rs59072263 gives an expected odds ratio for POAG of 1.08 (57) and, in untreated glaucoma (in our discovery cohort, in common with other cohorts from the developed world, over 50% of glaucoma was undiagnosed at initial ascertainment) (58), each copy of the allele gives an extra 6% likelihood of significant deterioration in vision (59). The novel associated SNP has a higher risk allele frequency (RAF) than previously reported SNPs (RAF ¼ 0.88 for rs59072263 when compared with 0.58 for rs11656696 and 0.12 for rs7555523) so a higher proportion of our cohort will carry at least one potentially deleterious allele.
In summary, we have performed a GWAS which has newly identified a region at 7p21 associated with IOP; the region is between the ICA1 and GLCCI1 genes, both of which may plausibly be invoked in the determination of IOP.

Ethics statement
The Blue Mountains Eye Study received ethical approval by the Western Sydney Area Health Service Human Ethics Committee in 1991. Written, informed consent was obtained from all participants. The EPIC-Norfolk study received approval from the East Norfolk and Waveney Research Governance Committee on 08/03/06 and from the Norfolk1 Research Ethics Committee on 06/02/06.

Discovery samples
The IOP data from the discovery cohort were gathered as part of the Blue Mountains Eye Study (BMES), a survey of vision and common eye diseases in the Blue Mountains region west of Sydney, Australia. The study was approved by the Western Sydney Area Health Service Human Ethics Committee. Written, informed consent was obtained from all participants. The population has been described in detail in a previous report (60). To summarize, a door-to-door census of the study region was carried out based on maps from the Australian Bureau of Statistics. All permanent non-institutionalized residents with birthdates before January 1, 1943 were invited to attend for eye examination. Of the 4433 eligible individuals, 3654 (82.4%) were examined. If the 278 individuals who had died or moved away from the area are excluded, this translates to an 87.9% response rate, which compares well with most population-based research in glaucoma (61)(62)(63). Subjects attending the examination underwent IOP measurement by applanation tonometry using a Goldmann applanation tonometer (Haag-Streit, Bern, Switzerland). Goldmann applanation tonometry is the international standard for accurate determination of IOP in ophthalmic clinical practice. Tonometry could not be performed in 12 subjects (0.3%). In a further nine uniocular subjects, IOP could only be measured unilaterally, in such cases, the value for that single eye was used. Data from subjects receiving past or present treatment designed to lower IOP were excluded from the analysis.

Replication samples
The European Prospective Investigation into Cancer (EPIC) was conceived as a pan-European study of the genetic and environmental determinants of cancer (64). The EPIC-Norfolk cohort The TwinsUK data were obtained from a panel of individuals who were recruited from the UK Adult Twin Registry (40) based at St Thomas' Hospital, London. The subjects were twin volunteers from the general population. Subjects were recruited for studies other than eye studies, and subsequently asked to attend for an eye examination. All subjects provided informed consent, and the study was reviewed by the Local Research Ethics Committee. The methods adopted in this study adhered to the tenets of the Declaration of Helsinki Principles. IOP measurement was performed with the Goldmann applanation tonometer and the ORA (IOPg). Data from subjects receiving past or present treatment designed to lower IOP were excluded from the analysis.

DNA preparation
DNA was extracted from whole blood. Quality was validated using the Sequenom iPLEX assay designed to genotype four gender SNPs and 26 SNPs present on the Illumina Beadchips. DNA concentrations were quantified using a PicoGreen assay (Invitrogen) and an aliquot assayed by agarose gel electrophoresis. A DNA sample was considered to pass quality control if the DNA concentration was ≥50 ng/ml, the DNA was not degraded, the gender assignment from the iPLEX assay matched that provided in the patient data manifest and genotypes were obtained for at least two-thirds of the SNPs on the iPLEX.

Genotyping
Discovery samples were genotyped at the Sanger Institute on the Illumina Infinium platform using the Human660W-Quad, a custom chip designed by WTCCC2 and comprising Human550 supplemented with 60 000 additional probes that were intended to allow the genotyping of common CNVs from the Structural Variation Consortium (66). Replication genotyping for the UK EPIC cohort was carried out at the Sanger Institute using the Illumina Immunochip, a custom chip designed by the Immunochip Consortium and WTCCC2, comprising 196 524 SNPs. For both chips, bead-intensity data were processed and normalized for each sample in BeadStudio; data for successfully genotyped samples were extracted and genotypes called using Illuminus. Replication genotyping on the Sequenom plex was also carried out at the Sanger Institute. Imputation was performed with reference to HapMap release 22 CEU by using IMPUTE2 (33). Genotyping for the TwinsUK cohort was carried out by using Illumina (San Diego, CA, USA) genotyping platforms; the Human Hap 300k Duo and Human Hap610 Quad array. All SNPs passed quality control criteria (Hardy -Weinberg equilibrium P . 0.001, minor allele frequency of at least 0.04, genotyping success rate for the SNP at least 95%). Imputation was performed with reference to the 1000 Genomes panel using minimac (http://genome.sph.umich.edu/wiki/Minimac) and association analysis was performed using mach2qtl.

Quality control
SNPs were excluded if the Fisher information for the allele frequency was not close to unity (information ,0.98) or if the minor allele frequency (MAF) was very low (defined as ,0.01%), or for extreme departures from Hardy -Weinberg equilibrium (HWE P-value ,10 220 ). Full details of the quality control methods employed in the Wellcome Trust Case Control Consortium 2 (WTCCC2) have been published elsewhere (67). Also see Supplementary material online, SNP Imputation.
For quality control of samples, a Bayesian clustering method was used to infer and exclude outlying individuals on the basis of ancestry, call rate, heterozygosity and signal intensity (67). To remove signal intensity outliers observed for raw intensity data, the difference between the A channel intensity and the B channel intensity was averaged over all SNPs on autosomes for each sample. A similar approach was used taking intensity measures from the A channel on the non-pseudo autosomal X chromosomes to identify outliers and infer gender. Samples were removed if their inferred gender was discordant with the recorded gender after cross-checking with original database entries, or if ,90% of the SNPs typed by Sequenom on entry to sample handling (discussed earlier) agreed with the genomewide data. To obtain a set of putatively unrelated individuals, we estimated genome-wide identity by descent (IBD) given identity by state (IBS) for all pairs of individuals using a Hidden Markov Model. One of each pair of related individuals where zero alleles shared IBD (IBD0) was ,95% was removed in the initial association analysis, excluding the member of the pair with the lowest call rate. Where unexpected duplicates were identified, both of the identical samples were removed to avoid incorrect phenotype assignment. Ancestry outliers were removed by projecting the WTCCC2 individuals on to the first two principal components of a PCA of the CEU, YRI and JPT/CHB HapMap individuals. Samples were excluded based on evidence of non-European ancestry.
A total of 463 individuals were excluded from the discovery cohort following these initial quality control checks.

Genomes imputation
We used 120 CEU phased haplotypes from the 1000 Genomes project (34) June 2010 release for the haploid reference panel. After imputation, a total of 7 642 395 SNPs were available for analysis, of which 712 4701 were imputed. SNPs were filtered using an information threshold of 0.7, an allele frequency threshold of 0.1% and Hardy -Weinberg filter of P , 10 220 giving a final total of 6 235 970 SNPs for analysis (5 718 276 imputed, 517 694 genotyped).

Pathway analysis
We undertook a pathway analysis using the approach reported previously in Sawcer et al. (68). Briefly, we took all SNPs with P , 10 24 in the discovery data and we used these to define association regions (with at least one such SNP) and compiled a list of all genes in these regions. This list is matched to the Gene Ontology (GO) pathway database. A Fisher's exact test is used to give a P-value for each pathway.
There were 513 genes in our associated regions of which 394 were in the GO database. The Fisher's exact test highlighted 10 pathways with P , 10 24 namely: anchored to membrane; monooxygenase activity; N,N-dimethylaniline monooxygenase activity; glucuronosyltransferase activity; alkaline phosphatase activity; phospholipid scramblase activity; xenobiotic metabolic process; response to xenobiotic stimulus; cellular response to xenobiotic stimulus; phospholipid scrambling.
Of the pathways identified, those of potentially most interest from a mechanistic viewpoint are 'phospholipid scramblase activity' and 'phospholipid scrambling'. This metabolic activity is a key part of lipid translocation, which drives endo-and exocytosis of microvesicles (69,70). It is thus plausible that this metabolic activity is important in the transcytotic movement of aqueous humour across the trabecular meshwork which is a key determinant of IOP as already mentioned in the manuscript.
We also implemented a number of slight variations on this analysis as in Sawcer et al. (68), for example, only including the nearest gene to an associated SNP. These alternative analyses also failed to pinpoint any pathways with a natural biological connection to IOP biology.

Phenotype
Mean IOP was regressed against age and sex and all eight samples with residuals .10 were removed from the primary analysis. From the replication analysis, 21 individuals with residuals .10 were removed. In the primary analysis, these phenotype outliers were excluded because they can lead to false positive associations. In order to check that we were not missing any true associations by excluding outliers, the analysis was repeated with inclusion of the phenotypic outliers, and a small number of regions (4) were chosen for replication based on this scan. One region for replication was also taken from the analysis fitting a general model with residual outliers included.
Regions were shortlisted for replication if they had at least one SNP with P , 10 25 in the primary scan, or P , 10 25 in the scan including the eight residual outliers and P , 10 24 in the primary scan. Where possible, the SNP with the smallest P-value in each region was replicated, along with the second best hit.

SUPPLEMENTARY MATERIAL
Supplementary Material is available at HMG online.