Large-scale multitrait genome-wide association analyses identify hundreds of glaucoma risk loci

Abstract

Glaucoma, a leading cause of irreversible blindness, is a highly heritable human disease. Previous genome-wide association studies have identified over 100 loci for the most common form, primary open-angle glaucoma. Two key glaucoma-associated traits also show high heritability: intraocular pressure and optic nerve head excavation damage quantified as the vertical cup-to-disc ratio. Here, since much of glaucoma heritability remains unexplained, we conducted a large-scale multitrait genome-wide association study in participants of European ancestry combining primary open-angle glaucoma and its two associated traits (total sample size over 600,000) to substantially improve genetic discovery power (263 loci). We further increased our power by then employing a multiancestry approach, which increased the number of independent risk loci to 312, with the vast majority replicating in a large independent cohort from 23andMe, Inc. (total sample size over 2.8 million; 296 loci replicated at P < 0.05, 240 after Bonferroni correction). Leveraging multiomics datasets, we identified many potential druggable genes, including neuro-protection targets likely to act via the optic nerve, a key advance for glaucoma because all existing drugs only target intraocular pressure. We further used Mendelian randomization and genetic correlation-based approaches to identify novel links to other complex traits, including immune-related diseases such as multiple sclerosis and systemic lupus erythematosus.

Fig. 1. Study design.

MGB Biobank, Mass General Brigham Biobank; DD, disc diameter; LDSC, linkage disequilibrium score regression.

Fig. 2. Manhattan plots displaying POAG GWAS P values.

a, Plot shows 263 loci from POAG MTAG in the European ancestry population. b, Plot shows 312 loci from POAG multiancestry meta-analysis. In these plots, the y axis shows the P values of SNPs in log–log scale. The red horizontal line is the genome-wide significance level at P = 5 × 10⁻⁸. SNPs with P < 1 × 10⁻⁴ are not shown in the plots. Previously unknown loci are highlighted with red dots, and the nearest gene names are in black text. Known SNPs are highlighted with purple dots, and the nearest gene names are in purple text. All tests were two-sided.

Fig. 3. Comparison of the effect sizes for genome-wide significant independent SNPs.

a, Plot shows 261 (two SNPs were unavailable in 23andMe) genome-wide significant independent SNPs identified from POAG MTAG in the European population versus glaucoma GWAS in 23andMe. The Pearson’s coefficient is 0.966 (P = 5.99 × 10⁻¹⁵⁴, n = 261 independent SNPs). b, Plot shows effect sizes for 302 (10 SNPs were unavailable in 23andMe) genome-wide significant independent SNPs identified from the POAG multiancestry meta-analysis versus glaucoma GWAS in 23andMe. The Pearson’s coefficient is 0.958 (P = 1.22 × 10⁻¹⁶⁴, n = 302 independent SNPs). In a and b, previously unknown SNPs are colored in red. c,d, Plots (n = 261 independent SNPs) show POAG MTAG from the European ancestry population versus POAG GWAS from Asian (c) and African (d) ancestry populations. The dots show the effect sizes of SNPs, and the error bars show the 95% confidence interval of the estimations of SNP effect sizes. AFR, African ancestry; ASN, Asian ancestry; EUR, European ancestry.

Fig. 4. Classification of POAG loci into VCDR- or IOP-specific SNPs.

a, Plot shows hierarchical clustering of 263 MTAG POAG loci (a multitrait colocalization approach is shown in Extended Data Fig. 6). Based on the z scores (dashed line shows y = x where the z scores of VCDR and IOP are equal), there is a subset of SNPs that act primarily via IOP (n = 171, dots) and a subset of SNPs acting primarily via VCDR (n = 92, squares). b, Plot displays the effect sizes on VCDR and on POAG for the 92 VCDR-specific SNPs. c, Plot displays the effect sizes on IOP and on POAG for the 171 IOP-specific SNPs. The dots show the effect sizes of SNPs, and the error bars show the 95% confidence interval of the estimations of SNP effect sizes.

Fig. 5. Putative causally associated traits with POAG.

Plots show 14 traits that were associated with POAG from MR (FDR P < 0.05, direction: complex traits → POAG). Different outcome traits are shown in different colors. Different MR methods are displayed in different line types. The dots show the effect sizes of MR estimations, and the error bars show the 95% confidence interval of the estimations. The number of SNPs, effect sizes and P values are presented in Supplementary Table 14. All tests were two-sided. BMI, body mass index; 95% CI, 95% confidence interval.

Extended Data Fig. 1. Comparison of the effect sizes for genome-wide significant independent SNPs by known loci and previously unknown loci.

a, Plot showing effect sizes for known genome-wide significant independent SNPs identified from the POAG multi-trait GWAS in European ancestry versus glaucoma GWAS in 23andMe. b, Plot showing 81 previously unknown genome-wide significant independent SNPs identified from the POAG multi-trait GWAS in European ancestry versus glaucoma GWAS in 23andMe. The Pearson’s coefficient is 0.94 (P = 1.42 × 10⁻³⁸). c, Plot showing effect sizes for known genome-wide significant independent SNPs identified from the POAG multi-ancestry meta-analysis versus glaucoma GWAS in 23andMe. d, Plot showing 109 previously unknown genome-wide significant independent SNPs identified from POAG multi-ancestry meta-analysis versus glaucoma GWAS in 23andMe. The Pearson’s coefficient is 0.939 (P = 1.57 × 10⁻⁴⁸). For the 81 previously unknown genome-wide significant independent SNPs identified from the POAG multi-trait GWAS in European ancestry, the replication rates in an independent cohort using 23andMe were: 38% SNPs (n = 31) passed the genome-wide significance level (P < 5 × 10⁻⁸) in the 23andMe study, 73% SNPs (n = 59) were significant after Bonferroni correction (P < 0.00062), and 96% SNPs (n = 78) reached a nominal significance level (P < 0.05). For the 109 previously unknown genome-wide significant independent SNPs identified from POAG multi-ancestry meta-analysis, the replication rates in an independent cohort using 23andMe were: 38% SNPs (n = 39) passed the genome-wide significance level (P < 5 × 10⁻⁸) in the 23andMe study, 66% SNPs (n = 68) were significant after Bonferroni correction (P < 0.0005), and 96% SNPs (n = 99) reached a nominal significance level (P < 0.05). The dots show the effect sizes of SNPs, and the error bars show the 95% confidence interval of the estimations of SNP effect sizes.

Extended Data Fig. 2. Scatterplot of minor allele frequency (MAF) and effect size (absolute value) for the novel loci in the multi-trait POAG GWAS analysis.

‘1’ corresponds to novel SNPs; ‘0’ corresponds to known SNPs.

Extended Data Fig. 3. Quantile-quantile plots for POAG GWAS.

a, Quantile-quantile plot for multi-trait POAG GWAS in participants of European ancestry. The quantile-quantile plot is based on one million randomly selected SNPs. Linkage disequilibrium (LD) score regression intercept is used to assess the genomic inflation; the intercept is 0.957 (standard error = 0.013, attenuation ratio < 0), and the lambda value is 1.27. b, Quantile-quantile plot for POAG multi-ancestry meta-analysis. The lambda value is 1.28. Because of the multi-ancestry design, LDSC was not performed.

Extended Data Fig. 4. Replication in 23andMe based on quality control annotation.

Comparison of the effect sizes (log odds ratio) for genome-wide significant independent SNPs in our discovery studies and the 23andMe replication cohort. a, Plot shows genome-wide significant independent SNPs identified from the POAG multi-trait GWAS in the European population versus glaucoma GWAS in 23andMe (n = 261 independent SNPs). b, Plot shows genome-wide significant independent SNPs identified from the POAG multi-ancestry meta-analysis versus glaucoma GWAS in 23andMe (n = 302 independent SNPs). The two different colors for ‘not pass’ and ‘pass’ represent the quality control annotation from 23andMe. Most of the 27 SNPs that did not pass the quality control annotated by 23andMe showed high concordance in effect size between the discovery and replication cohorts. The dots show the effect sizes of SNPs, and the error bars show the 95% confidence interval of the estimations of SNPs effect sizes.

Extended Data Fig. 5. Comparison of the effect sizes for 312 POAG genome-wide significant independent SNPs from multi-ancestry meta-analysis against their effect sizes in VCDR and IOP grouped by different P-value thresholds.

a, Plot showing the comparison with VCDR (n = 312 independent SNPs). The x-axis shows the effect sizes in multi-ancestry meta-analysis of POAG. The y-axis shows the effect sizes in VCDR. The SNPs are shown in different colors based on different P values in VCDR (P < 5 × 10⁻⁸; ‘<0.05/312’: 5 × 10⁻⁸ ≤ P < 0.05/312; ‘<0.05’: 0.05/312 ≤ P < 0.05; P ≥ 0.05). b, Plot showing the comparison with IOP (n = 312 independent SNPs). The SNPs are shown in different colors based on different P values in IOP (P < 5 × 10⁻⁸; ‘<0.05/312’: P < 5 × 10⁻⁸ ≤ P < 0.05/312; ‘<0.05’: 0.05/312 ≤ P < 0.05; P ≥ 0.05). The dots show the effect sizes of SNPs, and the error bars show the 95% confidence interval of the estimations of SNPs effect sizes. In our previous work, we have shown that the genetic correlation between VCDR and POAG is 0.5, and many VCDR significant loci are not necessarily associated with POAG. In this figure, many POAG SNPs are not associated with VCDR, and similarly many VCDR SNPs are not associated with POAG.

Extended Data Fig. 6. Classification of POAG loci into VCDR or IOP-specific SNPs based on multi-trait colocalization.

This scatter plot shows the effect sizes of 263 MTAG POAG loci on VCDR (x-axis) and IOP (y-axis). The assigned SNPs (IOP and VCDR in different point shapes) were based on a hierarchical clustering approach. The different point colors for trait combination show the results from multi-trait colocalization. This figure shows a consistent classification of POAG loci into VCDR- and IOP-specific SNPs based on the hierarchical clustering method and the multi-trait colocalization method.

Extended Data Fig. 7. Bivariate genetic correlation analysis identifies 24 traits that are genetically correlated with POAG, VCDR or IOP.

The x-axis shows the genetic correlations and their 95% confidence intervals. The y-axis shows the genetically correlated traits. The dots show the effect size of Mendelian randomization estimations, and the error bars show the 95% confidence interval of the estimations. All tests were two-sided.

Extended Data Fig. 8. Shared genomic region between systemic lupus erythematosus and POAG.

The genomic region (gene ATXN2) has a high posterior probability (PP4 = 0.98) in the Bayesian colocalization analysis. a, LocusZoom plot for systemic lupus erythematosus. b, LocusZoom plot for POAG.

Extended Data Fig. 9. Reverse-directional Mendelian randomization analysis for associated traits with POAG.

Plots show the reverse-directional Mendelian randomization for traits that were associated with POAG (celiac disease was not shown because of the limited number of overlapping SNPs; direction: POAG -> complex traits). Different exposure traits are shown in different colors. Different MR methods are displayed in different line types. In this reverse-directional MR analysis, only the association between VCDR and optic disc area passed multiple testing.