Abstract
Purpose.:
Intraocular pressure (IOP) is an important clinical parameter in the evaluation of ocular health. Elevated IOP is a major risk factor for primary open-angle glaucoma (POAG). The goal of this study was to identify rare and less common variants that influence IOP.
Methods.:
We performed an exome array analysis in a subset of 1660 individuals from a population-based cohort, the Beaver Dam Eye Study. Associations with IOP were tested on 45,849 single nucleotide variants and 12,390 autosomal genes across the genome.
Results.:
Intraocular pressure was suggestively associated with novel variants located in FAR2 at 12p11.22 (rs4931170, P = 1.2 × 10−5), in GGA3 at 17q25.1 (rs52809447, P = 6.7 × 10−5), and in PKDREJ at 22q13.31 (rs7291444, P = 7.4 × 10−5). Gene-based analysis found suggestive associations between IOP and the genes HAP1, MTBP, FREM3, and PHF12. We successfully replicated the associations with GAS7 (P = 7.4 × 10−3) for IOP, and also identified a previously reported POAG locus in the CAV1/CAV2 region to be associated with IOP (P = 3.3 × 10−3). This association was confirmed in a meta-analysis with three published genome-wide association studies (Pcombined = 4.0 × 10−11).
Conclusions.:
Our results suggest that novel genetic variants and genes with multiple, less common variants may play a role in the control of IOP. The implication of the caveolin genes, CAV1/CAV2, as a common genetic factor influencing both IOP variations and POAG may provide new insights of the underlying mechanism leading to glaucoma and glaucomatous visual field loss.
Intraocular pressure (IOP) is an important physiologic characteristic in maintaining the structure and function of the eye. Normal IOP is regulated by the balance between the production and the drainage of aqueous humor through the trabecular meshwork in the anterior chamber angle of human eyes. Elevated IOP contributes to the accelerated death of retinal ganglion cells and the damage of the optic nerve, and ultimately results in progressive visual field loss and irreversible blindness in glaucoma patients.
1 It is a major risk factor of glaucomatous optic neuropathy, and the only target for current glaucoma therapy. Clinical trials have shown that IOP-lowering treatments in eyes with ocular hypertension are effective both in preventing glaucoma onset and delaying its progression.
2,3
Intraocular pressure is a heritable polygenic trait with environmental influences.
4 The heritability for IOP is estimated to range from 0.30 to 0.42 in populations of European ancestry,
4–7 suggesting genetics may play an important role in IOP. Genome-wide studies have located regions on chromosome 10q22 and chromosome 19p with significant genetic linkage to IOP in family-based studies.
8,9 In addition, genome-wide association studies (GWAS) have identified multiple common single nucleotide polymorphisms (SNPs) associated with IOP,
10–13 some of which are in the overlapping GWAS regions for primary open-angle glaucoma (POAG).
14,15 The shared association variants for IOP and glaucoma suggest a possible common pathway for this endophenotype and disease.
The development of exome arrays that target protein-coding variation offers new opportunities to assess the role of rare and low-frequency coding variants in human complex traits. The primary goal of this study was to identify rare (minor allele frequency (MAF) < 1.0%) and less common (1.0% ≤ MAF ≤ 5%)
16 variants associated with IOP and to replicate previously reported GWAS loci for IOP in a population-based cohort of European ancestry from the Beaver Dam Eye Study (BDES).
The BDES is a population-based cohort study established in Beaver Dam, Wisconsin, in 1987.
17 From 1987 to 1988, the BDES completed a census of the city and township of Beaver Dam and identified 3715 households with at least one occupant aged from 43 to 84 years. Of the 5925 individuals eligible to enroll, 4926 (83.14%) individuals underwent an ocular examination and a personal history questionnaire during 1988–1990.
18 Pedigrees have been constructed for 2783 participants with known familial relationships in the catchment area of the study.
4 To improve the sampling efficiency while maintaining the full spectrum of the phenotypes, we selected individuals from the full BDES cohort and included those with extreme baseline IOP or refractive error measurements. This sampling resulted in a random sample with the full spectrum of IOP phenotype represented, since refractive error extremes were also included (
Supplementary Fig. S1). The study followed the recommendations of the Declaration of Helsinki. Informed consent was obtained from all study participants and the institutional review board at the University of Wisconsin approved all protocols.
Intraocular pressure was measured with a Goldmann applanation tonometer.
18 Unreliable measurements, as assessed by the trained observers, were excluded from the analysis. Blood pressure was recorded as the mean of the second and third measurements according to the Hypertension Detection and Follow-up Program protocol.
19 A detailed medical history, including information about hypertension, diabetes, and other medical conditions, and history of medication use were obtained from each participant.
The higher IOP measurement of the two eyes at the baseline visit was used in these analyses. The distribution of intraocular pressure was approximately normal but slightly right skewed (
Supplementary Fig. S1). To reduce the impact of extreme outliers, we winsorized the data.
20 All individuals with the trait more extreme than 5 standard deviations from the mean were assigned the values at 5 standard deviations from the mean (
Supplementary Fig. S1). Before the association analysis, the winsorized trait was first adjusted in a multiple linear regression model for potential confounding variables including age, sex, systolic blood pressure, and IOP treatment. The residual IOP values were used for the association analysis. In this study, we present the adjusted IOP values (sum of the residuals and the intercept from the regression model) with mean centered age and systolic blood pressure.
A total of 1908 BDES participants underwent genotyping with the Illumina HumanExome BeadChip (San Diego, CA, USA). Genotyping was performed at the Genetic Resources Core Facility at Johns Hopkins Institute of Genetic Medicine. Genotyping calling was done with Illumina's GenTrain clustering algorithm in GenomeStudio.
Standard quality control measures were used. Samples with sex inconsistencies that could not be easily resolved (
n = 15) were excluded. All samples had a call rate > 98%. We assessed the cryptic relatedness on the basis of pairwise identity by descent (IBD) sharing and removed one individual from each pair of samples that represented unexpected duplicates or first- and second-degree relatives (IBD sharing > 20%,
n = 126). Samples with unavailable IOP measurement of either eye or missing values of any covariates were also excluded (
n = 107). Of the 242,901 genotyped variants, we excluded nonautosomal variants (
n = 5465), variants with genotyping call rate < 98% (
n = 6038), as well as monomorphic variants (
n = 130,783). The concordance rate with 44 HapMap controls was 99.82% and the genotype concordance among 21 masked duplicate sample pairs was 99.99%. Since departures from Hardy-Weinberg equilibrium (HWE) can be a result of evolutionary forces, genotyping errors, population admixture, and marker-trait associations,
21 we evaluated HWE on the final set of variants, and reported the
P value estimated from the HWE exact test for significant associations. Our final data set includes 1660 individuals and 100,615 polymorphic variants.
The BDES population self-reported as primarily European-American. To evaluate the genetic ancestry we performed principal component analysis (PCA) by using SMARTPCA in EIGENSTRAT
22 with six HapMap phase III populations (
Supplementary Fig. S2). There was no increase in genomic inflation (
λ = 1.007) (
Supplementary Fig. S3).
We tested single variants for association with IOP under an additive linear regression model in PLINK.
23 We considered rare (0.3% ≤ MAF < 1%,
n = 10,744) as well as low-frequency and common variants (MAF ≥ 1%,
n = 35,093) in the association analysis. However, given the sample size of this study, variants with a MAF < 1% may be prone to bias. Variants with a MAF ≥ 1% correspond to at least 33 copies of the minor allele in the analyzed cohort, which should allow stable estimates of sampling errors and provides adequate statistical power to detect low-frequency variants with effect size (
β) greater than 3.5 at
α = 0.0001 (estimated in QUANTO
24).
Linkage disequilibrium (LD)–based pruning identified 33,056 variants with pairwise correlation < 0.20. We set the significance threshold across the exome to P < 1.5 × 10−6, corresponding to a Bonferroni correction for these 33,056 independent tests.
Genome-wide association studies have identified seven loci associated with IOP and/or POAG. For IOP these include
TMCO1 (1q22-q25),
10,11,14 GLCCI1/ICA1 (7p21),
12 MVB12B (9q33.3),
13 and
GAS7 (17p13.1).
10,11 For POAG these include
CAV1/CAV2 (7q31),
15,25,26 CDKN2B-AS1 (9p21.3),
14,27–29 and
SIX1/SIX6 (14q23).
27,29 For replication, we considered (1) the exact variant, and (2) variants in strong LD (
r2 > 0.70) with the reported variant in the same GWAS loci. For variants in previously reported IOP loci, associations were considered as replicated if the
P value was <0.01. For variants in previously reported POAG loci, we evaluated their associations with IOP in our study and conducted a meta-analysis on those with
P < 0.01, combining our results with published studies for IOP in a fixed-effect model using an inverse sample-size weighted approach implemented in METAL.
30
We performed gene-based analysis by using SKAT-O implemented in R.
31–33 The residual IOP values were regressed on rare and less common variants (MAF ≤ 5%) in a gene region, allowing them to have different directions and magnitude of effects. Gene regions with at least two variants in desired MAF range were included. We defined the cumulative minor allele frequency (CMAF) as the sum of the MAF of all variants in a gene region. To ensure accurate statistical inference for this sample size, we focused this analysis on gene regions with CMAF ≥ 1.0%, although we present associations for all CMAF. A total of 12,390 genes (75,936 variants) across the autosomal chromosomes were tested in gene-based analysis. To optimize the statistical power, we applied the default weight
Display Formula
in SKAT-O, which up-weights the rare variants. The significance threshold across all identified gene regions was set to a Bonferroni corrected
P < 4.0 × 10
−6 (
n = 12,390 tests).
Annotation was performed by using SeattleSeq Annotation Server 138 version 9.01 under GRCh37/hg19 (
http://snp.gs.washington.edu/SeattleSeqAnnotation138/; provided in the public domain by the University of Washington, Seattle, WA, USA). Except for 57 insertions and deletions, all variants that passed quality control filters were annotated.
The characteristics of the BDES participants for the analyzed cohort and for the entire population are summarized in
Table 1. Females constituted 57.3% of the analyzed cohort and the mean age was 60 years for both men and women. The distribution of IOP in our cohort had a mean 15.6 mm Hg with standard deviation 4.1 mm Hg, which is comparable to other large population-based cohorts of European descent and to the full BDES study. The mean IOP was 15.9 mm Hg in the full BDES study (
Table 1); it was 16 and 15.4 mm Hg in the Blue Mountains Eye Study
34 and the NEIGHBOR consortium,
35 respectively. Intraocular pressure was normally distributed with a similar mean and median measurement (15.6 mm Hg versus 15.0 mm Hg, respectively). After winsorization and covariate adjustment, the distribution of adjusted IOP remained approximately normal (
Supplementary Fig. S1). In all, 8.3% of the individuals had elevated IOP, defined as >21 mm Hg, and 2.5% had reported IOP medical treatment. There were 45 individuals with self-reported glaucoma diagnosis at either eye.
We tested 35,093 low-frequency and common variants (MAF ≥ 1%) individually for association with IOP (
Fig. 1) and no variants met the threshold for significance (
P < 1.5 × 10
−6) (
Table 2). However, we did identify suggestive associations. A common variant (rs4931170, MAF = 38.3%) on chromosome 12p11.22 in the intronic region of the fatty acyl-CoA reductase 2 gene (
FAR2) was associated with elevated IOP. Each copy of the minor allele G of rs4931170 led to 0.58–mm Hg (95% confidence interval [CI] = 0.33–0.83,
P = 1.2 × 10
−5) increase in IOP (
Fig. 2A).
Additionally, a low-frequency coding variant (rs52809447, MAF = 1.08%) that results in the substitution of glutamic acid with glycine (p.Glu97Gly) on chromosome 17q25.1 in the exonic region of Golgi-associated, gamma adaptin ear containing, ARF-binding protein 3 gene (
GGA3) was associated with IOP. Individuals with one copy of the minor allele C have on average an IOP 2.43 mm Hg (95% CI = 1.23–3.63,
P = 6.7 × 10
−5) higher than noncarriers (
Fig. 2B). For this low-frequency variant there were no individuals with two copies of the C allele in this study.
The third suggestive locus for IOP was a common variant on chromosome 22q13.31 (rs7291444, MAF = 15.8%), substituting a threonine with a proline (p.Thr992Pro) in the polycystin (PKD) family receptor for egg jelly gene (
PKDREJ). Each copy of the minor allele G of rs7291444 was associated with 0.69–mm Hg (95% CI = 0.36–1.02,
P = 7.4 × 10
−5) increase in IOP (
Fig. 2C).
Analysis of rare variants (0.3% ≤ MAF < 1.0%) identified eight variants in six genes that had large effects on IOP (
β > 3.0) with suggestive significance (
Table 2). Variants in three of these gene regions,
MTBP (rs61753798),
FREM3 (rs72938299, rs17017968, and rs7679078), and
RNASET2 (rs41269593), were also identified in gene-based analyses.
Except for the
MVB12B and
GLCCI1/ICA1 region, our exome array contained either the exact variant or at least one variant in strong LD for five published GWAS loci:
TMCO1 and
GAS7 for IOP, and
CAV1/CAV2,
CDKN2B-AS1, and
SIX1/SIX6 for POAG (
Table 3).
We successfully replicated the previously reported association in growth arrest–specific 7 gene (
GAS7) region for IOP.
GAS7 on chromosome 17p13.1 was first identified as a susceptibility locus for IOP in a meta-analysis of a large discovery cohort from the Netherlands.
10 The minor allele A of the lead variant rs11656696 in
GAS7 region was associated with lower IOP (MAF = 43%,
β = −0.26,
P = 9.8 × 10
−9).
10 An independent meta-analysis later identified an association in the
GAS7 region on a different variant rs12150284 (
β = −0.49,
P = 2.4 × 10
−6), a locus 2.6-kb upstream from rs11656696 (
r2 = 0.62).
11 Neither of these variants was genotyped in our data set. In our study, rs9897130 was the top
GAS7-associated SNP (
Fig. 3A), a locus that is highly correlated with rs12150284 (
r2 = 0.70) and is also associated with lower IOP (MAF = 48.4%,
β = −0.33,
P = 7.4 × 10
−3). In these independent studies, all three variants in the
GAS7 region showed associations with lower IOP and similar magnitudes of effects.
Among the variants in published GWAS loci for POAG, rs4236601 on chromosome 7q31, which contains the caveolin genes
CAV1 and
CAV2, showed a significant association with IOP (
Table 3). This locus was first discovered with POAG in an Icelandic population (odds ratio [OR] = 1.36, 95% CI = 1.23–1.50,
P = 5.0 × 10
−10)
15 and was later replicated in two independent studies of POAG in populations of European ancestry.
25,26 In these POAG studies, cases consist of both normal-tension and high-tension POAG patients, but Wiggs et al.
25 have shown no association among high-tension POAG individuals for rs4236601. In our study the minor allele A of variant rs4236601 was associated with elevated IOP (
Fig. 3B;
β = 0.42,
P = 3.3 × 10
−3) and is consistent in direction with these other studies.
Our finding is in agreement with three GWAS of IOP for rs4236601 in European ancestry populations (van Koolwijk et al.
10: MAF = 29%,
β = 0.19,
P = 1.1 × 10
−4; Ozel et al.
11:
β = 0.38,
P = 2.8 × 10
−4; Strange et al.
12:
P = 1.5 × 10
−3). We performed a meta-analysis combining our results with these three studies and identified an overall genome-wide significant association for IOP at rs4236601 (
Table 4,
Pcombined = 4.0 × 10
−11) in the
CAV1/CAV2 gene region (
Supplementary Fig. S4).
Gene-based tests offer an alternative to single-variant analysis, which is often underpowered to detect associations with rare variants. Among the 6900 genes with a CMAF ≥ 1%, no gene regions were significant (
P < 4.0 ×10
−6). Suggestive associations between IOP and four gene regions were identified (
Table 5): the huntingtin-associated protein 1 gene (
HAP1) on chromosome 17q21.2-q21.3 (CMAF = 2.11%,
P = 9.3 × 10
−6), the MDM2-binding protein gene (
MTBP) on chromosome 8q24.12 (CMAF = 1.26%,
P = 2.6 × 10
−5), the FRAS1-related extracellular matrix 3 gene (
FREM3) on chromosome 4q31.21 (CMAF = 4.49%,
P = 9.4 × 10
−5), and the PHD finger protein 12 gene (
PHF12) on chromosome 17q11.2 (CMAF = 1.14%,
P = 9.6 × 10
−5).
For the 5490 gene regions with a CMAF < 1.0% (
Table 5), the microfibrillar-associated protein 2 gene (
MFAP2) on chromosome 1p36.1–p35 reached significance (CMAF = 0.6%,
P = 6.2 × 10
−8). Two singleton variants were identified in this gene: one missense mutation (p.Val130Ile) in an individual with an IOP of 36 mm Hg and the other missense mutation (p.Gln33Pro) in an individual with an IOP of 10 mm Hg. As with
MFAP2, other reported gene regions with a CMAF < 1.0% contained mostly singleton and doubleton variants that were seen in only a few individuals in our study (
Table 5).
In this study, we performed an exome array analysis on 1660 BDES participants to identify rare and less common variants influencing intraocular pressure. We successfully replicated the associations of GAS7 in our analysis. Our meta-analysis identified a novel significant association of IOP with CAV1/CAV2, a region previously identified to be a POAG locus. We also reported suggestive associations for three novel variants in FAR2, GGA3, and PKDREJ and four gene regions (HAP1, MTBP, FREM3, and PHF12) with rare and low-frequency variants.
Our results lend support to a role for
CAV1/CAV2 in the pathogenesis of elevated IOP. Proteins encoded by caveolin-1 (
CAV1) and caveolin-2 (
CAV2) genes are major components of the caveolae plasma membranes. Caveolin-1 is expressed in several retinal cell types, including photoreceptor, retinal vascular endothelia cells, Müller glia, and retinal pigment epithelium cells, and has been linked to ocular pathologic processes including autoimmune uveitis, diabetic retinopathy, and POAG.
36 The role of caveolin-1 in the retina is largely unknown. Genetic ablation of caveolin-1 might cause retinal functional deficits due to disruptions in microenvironmental homeostasis
36 and blood–retina barrier breakdown.
37 A cytotoxic agent that increases IOP and aqueous outflow resistance in mice can also increase the expression of caveolin-1 in the treated human trabecular meshwork cells.
38 In a study by Elliot and colleagues (Elliott MH, et al.
IOVS 2014;55:ARVO E-Abstract 2888),
Cav-1 knockout mice experience prolonged elevation of IOP and significant reduction in pressure-dependent outflow, accompanied with morphologic change of caveolae in endothelial cells, suggesting CAV-1 may have a role in IOP regulation. Further investigations are necessary to elucidate the potential mechanisms of caveolin-1 and caveolin-2 in controlling IOP and its relationship with POAG. But this finding of a gene significantly associated with both IOP, an endophenotype, and POAG may yield new understandings of the underlying mechanisms leading to glaucoma and glaucomatous visual field loss.
Our suggestive SNP-associated genes may also have biologic relevance.
FAR2 belongs to the short-chain dehydrogenase/reductase superfamily whose protein products convert fatty acyl-CoA into fatty alcohols in wax biosynthesis, and the expression of FAR2 was found to be highest in the eyelid.
39 The products of the
GGA family, including
GGA3, regulate the protein trafficking between the trans-Golgi network and the lysosome, and have a significant role in Alzheimer disease pathogenesis, contributing to the increased accumulation of the amyloid-β protein (Aβ).
40
The gene-based tests are helpful in evaluating a gene whose causal variants may have different directionality and frequency. Interestingly, another gene identified in our gene-based analysis,
HAP1, which encodes a Huntington's disease–associated protein, may also play a role in the regulation of Aβ protein levels in neurons by controlling intracellular trafficking.
41 FREM3 is a member of an extracellular matrix protein family and has a strong implication in Fraser syndrome. Its analogous protein in mice presents high expression levels in a variety of retinal cells during eye development,
42 lending biologic plausibility to the associations in our study. The
MFAP2 gene association is most relevant. As a component of elastic microfibrils, MFAP2 is involved in the pathogenesis of exfoliation syndrome and exfoliative glaucoma.
43 However, the cumulative minor allele frequency for
MFAP2 is low (<1%) and we are cautious in interpreting the significance of these associations in the BDES cohort. These gene-based associations reflect extremely rare variants in this population (singletons and doubletons), and require much larger sample sizes, as well as exome or whole-genome sequencing, to fully capture the dispersion of variants.
The lack of significant associations in our study may be due to insufficient power. At a significance level of 0.0001, the single-variant analysis had >80% power to detect large effects (
β > 3.5) for variants with 1% MAF, but was underpowered for low-frequency variants with moderate or small effect sizes. The design of the exome array, which is constrained to protein-altering variants, has resulted in a limited number of polymorphic variants when the sample size is moderate. In our analyzed cohort of 1660 individuals, almost 70% of the variants in the genotyping assay were monomorphic or singletons. The remaining variants across the whole genome led to a mean coverage of six variants in each gene region. Although statistical approaches such as SKAT and SKAT-O are designed to detect the effects of low-frequency and rare variants,
32 the power of these approaches are diminished when the small number of variants in the sample population cannot provide adequate gene coverage.
Additional studies and data across populations are warranted and encouraged for a comprehensive assessment of the rare and low-frequency variants associated with IOP. These studies, coupled with existing genome-wide association studies of common variants, may start to elucidate the causal alleles and mechanisms for IOP and its association with POAG.
We thank the participants of the Beaver Dam Eye Study.
Supported by the National Eye Institute of the National Institutes of Health under Grants R01EY021531 and U10006594.
Disclosure: F. Chen, None; A.P. Klein, None; B.E.K. Klein, None; K.E. Lee, None; B. Truitt, None; R. Klein, None; S.K. Iyengar, None; P. Duggal, None