Genome-wide meta-analysis identifies 127 open-angle glaucoma loci with consistent effect across ancestries

Primary open-angle glaucoma (POAG), is a heritable common cause of blindness world-wide. To identify risk loci, we conduct a large multi-ethnic meta-analysis of genome-wide association studies on a total of 34,179 cases and 349,321 controls, identifying 44 previously unreported risk loci and confirming 83 loci that were previously known. The majority of loci have broadly consistent effects across European, Asian and African ancestries. Cross-ancestry data improve fine-mapping of causal variants for several loci. Integration of multiple lines of genetic evidence support the functional relevance of the identified POAG risk loci and highlight potential contributions of several genes to POAG pathogenesis, including SVEP1, RERE, VCAM1, ZNF638, CLIC5, SLC2A12, YAP1, MXRA5, and SMAD6. Several drug compounds targeting POAG risk genes may be potential glaucoma therapeutic candidates.

P rimary open-angle glaucoma (POAG) is the leading cause of irreversible blindness globally 1,2 . The disease is characterized by progressive optic nerve degeneration that is usually accompanied by elevated intraocular pressure (IOP). Neuroprotective therapies are not available and current treatments are limited to lowering IOP, which can slow disease progression at early disease stages; however, over 50% of glaucoma is not diagnosed until irreversible optic nerve damage has occurred 2,3 .
The majority of known risk loci for POAG have been identified through GWAS in participants of European descent, followed by replication in other ethnic populations. However, previous observational studies have shown that individuals of African ancestry, followed by Latinos and Asians, have higher POAG disease burden compared to those with European ancestry 3,[16][17][18] , suggesting important differences in genetic risk and highlighting the need to compare the genetic architecture of these ethnic groups.
In this study, we report the results of a POAG multi-ethnic meta-analysis on 34,179 cases and 349,321 controls, identifying 127 risk loci (44 not previously reported at genome-wide significance levels for POAG). The identified risk loci have broadly consistent effects across European, Asian, and African ancestries. We show that combining GWAS data across ancestries improves fine-mapping of the most likely causal variants for some loci. By integrating multiple lines of genetic evidence we identify the most likely causal genes, some of which might contribute to glaucoma pathogenesis through biological mechanisms related to extracellular matrix cell adhesion, intracellular chloride channels, adipose metabolism, and YAP/HIPPO signaling.

Results
Discovery of previously unreported POAG risk loci in Europeans. We performed a four-stage meta-analysis (Fig. 1). In the first stage, we conducted a fixed-effect meta-analysis of 16,677 POAG cases and 199,580 controls of European descent. The participating studies are detailed in Supplementary Data 1. We identified 66 independent genome-wide significant (P < 5e-08) single nucleotide polymorphisms (SNPs) (Supplementary Data 2), of which 16 were not previously identified (i.e., uncorrelated with previously reported SNPs). There was no evidence of inflation due to population structure (linkage disequilibrium (LD) score regression 19 intercept 1.03, se = 0.01, Supplementary Fig. 1A).
Significant POAG risk loci in Asians and Africans. In the second stage, we completed a fixed-effect meta-analysis of 6935 POAG cases and 39,588 controls of Asian descent, and a separate fixed-effect meta-analysis of 3281 POAG cases and 2791 controls of African ancestry (Supplementary Data 1). Ten loci were significantly associated with POAG (P < 5e-8) in the meta-analysis of Asian studies (Supplementary Data 3), all of which are known POAG loci, and at least nominally (P < 0.05) associated with the European meta-analysis. While only one of these loci had a P < 0.05 in Africans, eight had consistent direction of effects. For the African meta-analysis, one locus (rs16944405 within IQGAP1) reached the genome-wide significance level (P = 3e-08). This locus has not been previously reported for POAG, and in this study, was not associated with POAG in Europeans (P = 0.315) and Asians (P = 0.075) (Supplementary Data 3). The LD score regression intercept was 0.99 (se = 0.009) for Asians and 0.95 (se = 0.006) for Africans, suggesting that these results are not influenced by population structure.
Consistent genetic effect across ancestries. As part of the second stage, we replicated the stage 1 European Caucasian findings in an independent dataset comprising 7,286 self-reported cases and 107,362 controls of European descent from the UK Biobank study (UKBB) (Supplementary Data 1) 20 , as well as in the metaanalyzed Asian and African datasets described above. We replicated the European Caucasian results in each ancestry group separately (Supplementary Data 2 and Supplementary Fig. 2), followed by combining the three replication datasets in a fixed-   effect meta-analysis (17,502 cases and 149,741 controls) to maximize statistical power. Of the 66 significant loci in the European Caucasian meta-analysis, 56 were nominally significant (P < 0.05), and 37 after Bonferroni correction (P < 7.6e-04) in the meta-analyzed replication cohort. The effect sizes had a Pearson correlation coefficient (r) = 0.82 (Fig. 2).
There was moderately high cross-ancestry concordance both for genome-wide significant loci and across the genome. For the genome-wide significant SNPs, the European SNP effects were correlated (Pearson correlation coefficient (r) = 0. 68 Fig. 2B, C). Of the 68 SNPs available in the Asian meta-analysis, 60 (88%) showed the same direction of effect as European Caucasians, and of the 66 SNPs available in the African meta-analysis, 55 (83%) showed the same direction as European Caucasians. The genetic correlation across the genome estimated using the approach implemented in Popcorn v0.9.9 21 was even higher: r = 0.85 (95% CIs 0.70-1.00) for European-Asian and r = 0.75 (95% CIs −0.93 to 2.43) for European-African. Although the concordance amongst the top SNPs was clear for the European-African comparison, larger sample sizes will be required to narrow the CIs on the European-African genome-wide correlation estimate.
Discovery of previously unknown POAG risk loci in the crossancestry meta-analysis. In the third stage, given the large genetic correlation between ancestries, we performed a fixed-effect metaanalysis of the results from stage 1 and 2 (34,179 cases vs. 349,321 controls) and identified 127 independent genome-wide significant loci, located at least >1 Mb apart ( Fig. 3 and Supplementary  Fig. 3). Of these, 44 loci were not previously associated with POAG at genome-wide significance levels (Supplementary Data 4). All loci identified in the European meta-analysis were also significant at the genome-wide level in the combined ancestry meta-analysis except for three loci, two of which were not previously identified (OVOL2 and MICAL3), and one previously reported (EGLN3/SPTSSA). Of note, four of the risk loci (MXRA5-PRKX, GPM6B, NDP-EFHC2, and TDGF1P3-CHRDL1) are on the X chromosome, representing the first POAG risk loci on a sex chromosome. We also identified an association of a human leukocyte antigen (HLA) gene (HLA-G/HLA-H) with POAG. All the lead SNPs have MAF > 0.01 in Europeans, except for two variants: rs74315329 (MAF = 0.0026 in 1000G Europeans) a well-known nonsense variant in MYOC 22,23 , and rs190157577 (MAF = 0.0013 in 1000G Europeans) an intronic variant in LINC02141/LOC105371299.
The POAG risk loci were strongly replicated in 23andMe. In the fourth stage, we validated the association of the genome-wide significant SNPs from stage 3 in a dataset comprising 43,254 participants with self-reported POAG (defined as those who reported having glaucoma excluding angle-closure glaucoma or other types of glaucoma) and 1,471,118 controls from 23andMe, Inc. Of the 127 loci, the association results for 125 SNPs were available in 23andMe, 120 of which (96%) were replicated at P < 0.05, and 106 (85%) after Bonferroni correction for 125 independent tests (Supplementary Data 4). The correlation of the effect size was r = 0.98 (95% CIs 0.977-0.989). In total, the genome-wide significant loci in this study (N = 127) collectively explain 9.4% of the POAG familial risk. The previously known loci (N = 83) explain 7.5% of the familial risk, and the previously unreported loci (N = 44) explain an additional 1.9%.
Most of the risk loci associated with POAG involve known glaucoma-related endophenotypes. Several highly heritable endophenotypes are related to POAG risk including IOP, structural variation of the optic nerve characterized as vertical cup-todisc ratio (VCDR) and variation in thickness of the retina cell layers including the retinal nerve fiber layer (RNFL) and the ganglion cell inner plexiform layer (GCIPL) 24 .
All POAG risk variants identified to date have also been associated with either IOP or with VCDR or both. We investigated the association of the POAG loci identified in this study with IOP and VCDR, using previous GWAS data for IOP (N = 133,492) 11 and for VCDR (N = 90,939) 15 . Figure 4a shows that the majority of loci (89 of 123; four were unavailable for IOP) are also associated with IOP (red and green dots on Fig. 4a). For the 34 loci with unclear effect on IOP (purple and blue dots on Fig. 4a, full data in Supplementary Data 5), we plotted the POAG effect sizes against VCDR effect sizes and determined that 24 of the 32 SNPs (2 SNPs unavailable for VCDR), have a clear effect on VCDR (purple dots on Fig. 4b). Eight of the POAG loci did not appear to have a clear effect on IOP or VCDR, although a small effect on glaucoma via a small change in IOP or VCDR could not be ruled out. The overall correlation of effect sizes between all POAG risk loci and IOP was 0.53, and between POAG and VCDR was 0.31, in line with previously published genetic correlation estimates 25 . To better visualize clustering of the POAG SNPs based on their effect on IOP/VCDR, we created Of the 127 POAG genome-wide significant lead SNPs, 116 were available in a GWAS of GCIPL thickness (N = 31,536; Supplementary Data 6). Of these, 14 loci were associated with GCIPL thickness at nominal significance threshold (P < 0.05) and four (PLEKHA7, MAPT, LINC01214-TSC22D2, and POU6F2) after Bonferroni correction for the number of tests (P < 0.05/116). Similarly, 13 loci were nominally associated with RNFL, and three (PLEKHA7, MAPT, and SIX6) after Bonferroni correction. These results suggest that these loci may impact glaucoma pathogenesis through modulation of retinal thickness.
Given that three POAG risk loci that we identified (loci containing MAPT, CADM2, and APP) have also been implicated in Alzheimer's disease (AD) and dementia [26][27][28] , we asked whether the same causal variants may underlie these loci and whether there is evidence for genetic sharing in any of the other genome-wide significant POAG loci. Applying the Bayesianbased colocalization method, eCAVIAR 29 , to the cross-ancestry and European POAG GWAS meta-analyses and the publicly available AD GWAS data 30 30,31 , there was weak support for colocalization at six loci (CLPP = 0.01-0.14; four loci from the cross-ancestry and three from the European POAG meta-analysis with one overlapping locus; see Supplementary Data 8), though none of these POAG loci reached genome-wide significance in the AD GWAS (AD variant P-values on the order of 10 −4 to 0.05). We note that the colocalization results with this larger AD GWAS meta-analysis 30,31 might be slightly inflated due to the large overlap of UK biobank samples between the POAG and AD meta-analyses. We further estimated the genome-wide genetic correlation between POAG and AD using LD score regression (LDSR) 32 that adjusts for sample overlap, on the two AD GWASs; 30,31 the correlation estimates were 0.03 (95% CIs: −0.11-0.16; P = 0.7) and 0.14 (95% CIs: 0.003-0.28; P = 0.049), respectively.
We next investigated genetic correlation between POAG and a range of other traits using bivariate LDSR 32 through the LD Hub platform (http://ldsc.broadinstitute.org/ldhub/). Only glaucoma, self-report glaucoma, and "Other eye problems" were significantly associated after adjustment for multiple testing for 758 traits (P < 6.6e-05; Supplementary Data 9). Some other traits in UKBB such as myopia (short-sightedness); systolic blood pressure; seeing a psychiatrist for nerves, anxiety, tension or depression; and suffering from nerves, showed some evidence for association at P < 0.003 (Supplementary Data 9).
Cross-ancestry fine-mapping. Incorporating GWAS data across European, Asian, and African ancestries allowed us to improve fine-mapping of the most likely causal variants. For 10 loci (including previously unidentified loci GJA1/HSF2, SEPT7, and MXRA5/PRKX), the posterior probability of finding a causal SNP Fig. 3 Manhattan plots for the cross-ancestry meta-analysis. Each dot represents a SNP, the x-axis shows the chromosomes where each SNP is located, and the y-axis shows −log10 P-value of the association of each SNP with POAG in the cross-ancestry meta-analysis (34,179 cases vs. 349,321 controls). The red horizontal line shows the genome-wide significant threshold (P-value = 5e-8; −log10 P-value = 7.30). The nearest gene to the most significant SNP in each locus has been labeled.
in Europeans improved after including Asian and African data (improvements from posterior probabilities <0.9 to >0.9 or from <0.8 to >0.8; Supplementary Data 10). For eight loci (of which THRB and SMAD6 were not known), although the posterior probability of a SNP being causal in Europeans was high (>0.9 at least for one SNP), there was still a slight improvement after including the other ancestries (Supplementary Data 10). In contrast, the cross-ancestry data made fine-mapping worse for three loci where the posterior probabilities in Europeans were >0.8 but declined to <0.8 after incorporating data from the other ancestries. For the rest of the loci, the posterior probabilities did not change significantly after including Asian and African data. Overall, the best causal SNPs in Europeans changed for 52 of 127 loci after including data from the other ancestries (Supplementary Data 10). For the remaining 75 loci, at least one SNP remained the best causal SNP in both fine-mapping using European data alone, as well as across ancestries, and 23 of the 127 lead SNPs identified in the meta-analysis remained the best causal SNP in cross-ancestry fine-mapping.
Gene-based and pathway-based results. We performed genebased and pathway-based tests using MAGMA v1.07b 33 for each ancestry separately, followed by combining P-values across ancestries using Fisher's combined probability test 34 . We identified 205 genes that passed the gene-based Bonferroni-corrected threshold (P < 0.05/20174), corresponding to an additional seven independent risk loci located at least >1 Mb apart from the risk loci identified in the single-variant-based test. Supplementary Data 11 presents significant genes within these seven loci. Expression of the risk genes identified in MAGMA gene-based analysis were significantly enriched in artery and nerve tissues, reflecting the widely recognized neuronal and vascular character of glaucoma ( Supplementary Fig. 5A, B).
Pathway analysis identified 21 significant gene-sets surviving the Bonferroni-corrected threshold (P < 0.05/10678). These included previously identified pathways such as collagen formation and vascular development 10,11 , and highlighted additional pathways involved in lipid binding and transportation such as apolipoprotein binding and negative regulation of lipid storage (Supplementary Data 12). Genes involved in these pathways that demonstrated suggestive association (P < 5e-05) with POAG in the MAGMA gene-based test are summarized in Supplementary Data 13.
Functional relevance of the identified POAG risk loci. We used multiple lines of genetic evidence to investigate the functional relevance of the identified risk loci, and to prioritize causal variants and target genes. A summary of these results for the previously unknown loci is provided in Table 1, with additional details presented in Supplementary Data 14. The following paragraphs describe these findings in further detail.
First, the relevance of the identified risk loci was investigated by examining their roles in regulation of gene expression, as well as chromatin interactions. Approximately 76% (96 out of 127) of the lead SNPs or those in high LD (r 2 > 0.8) with the lead SNPs have also been reported to be significant expression quantitative trait loci (eQTLs) (FDR < 0.05) in various tissues (Supplementary Data 14 and 15). Moreover, the identified risk loci have 34,724 unique significant (FDR < 1e-6) chromatin interactions in various tissues/cell lines involving 4882 genes ( Supplementary Fig. 6). Of these, 425 genes overlap with the eQTL genes (286 genes if only considering eQTL results for the genome-wide significant SNPs). In addition, the lead SNPs or SNPs with LD r 2 > 0.8 with the lead SNPs for 124 risk loci identified in this study overlap with one end of these chromatin interactions (Supplementary Data 14).
Additional support for the pathogenicity of the POAG risk loci comes from the predicted pathogenicity scores: 20 SNPs had CADD scores >12.37, suggesting that these SNPs have deleterious effects (Supplementary Data 16) 35 . Overall, three lead SNPs are protein-altering variants and 12 lead SNP are in high LD (r 2 > 0.8) with a protein-altering variant (Supplementary Data 14), suggesting pathogenic effects through protein-coding roles of these variants (e.g., rs61751937 a missense variant in SVEP1). In addition, 24 SNPs had RegulomeDB 36 scores ≤3, supporting regulatory roles for these SNPs (Supplementary Data 16).
To investigate which risk loci are more likely to affect POAG by modulating gene expression, we used two transcriptome-wide association study (TWAS) approaches: Summary Mendelian randomization (SMR) 37 and MetaXcan 38 . For MetaXcan, we used our POAG cross-ancestry meta-analysis statistics, RNA-seq and genotype data from peripheral retina (EyeGEx) 39 and 44 GTEX tissues. Following Bonferroni correction for the maximum number of genes tested (N = 7209) in 45 tissues (P < 1.5e-07), we identified 100 significant genes, which were selected as the most likely causal genes based on the integration of eQTL data (also see below and Supplementary Data 14 and 17). Of these significant genes, three (AKR1A1, DDIT4L/LAMTOR3, and C4orf29) were located >1 Mb apart from (the other loci) identified using single-  rs172531  Intron  RERE  NO  RPL7P11  RERE  RERE  NO  NO  RERE  NO  rs941125  Intron  GLIS1  NO  SLC25A3P1  GLIS1  NO  NO  NO  GLIS1  NO  rs4076000  Intergenic  ELOCP18/RPE65  NO  CTBP2P8  WLS  NO  NO  NO  NO  NO  rs12566440  Intron  GPR88  NO  VCAM1  VCAM1  VCAM1  NO  NO  GPR88  NO  rs4542196  Intron  DDR2  NO  DDR2  DDR2  NO  NO  NO  NO  DDR2  rs12623251  Intergenic  TRIB2  YES  AC013471.1  TRIB2  NO  NO  NO  TRIB2  NO  rs6713914  Intron  LOC105374754  NO  AC009970.1  NO  NO  NO  YES  NO  NO  rs12613800  Intron  ZNF638  NO  RNU6-105P  ZNF638  ZNF638  ZNF638  NO  ZNF638  NO  rs9852634  Intron  THRB  YES  RPL31P20  NO  NO  NO  NO  NO  NO  rs1500708  Intron  ARHGEF3  YES  ARHGEF3-AS1  a This column indicates whether the GWAS lead SNP was identified to be the most probable causal SNP across ancestries in the multi-ethnic fine-mapping analyses performed using the software PAINTOR. b Where the GWAS lead SNP or SNPs in LD r 2 > 0.8 with the lead SNP overlap with one end of a chromatin interaction, the most significant gene involved in this interaction has been shown. c Where the GWAS lead SNP or SNPs in LD r 2 > 0.8 with the lead SNP is an eQTL in any GTEx tissue or the other eQTL databases used in this study (see the Methods), the most significant target gene has been shown. d Where the GWAS lead SNP or SNPs in LD r 2 > 0.8 with the lead SNP is an eQLT in retina in EyeGEx study, the most significant target gene has been shown. e Where the GWAS lead SNP is a protein-altering variant or in LD r 2 > 0.8 with a protein-altering variant, the corresponding gene has been shown. f This column indicates whether the GWAS lead SNP has a CADD score > 12.37. g This column shows the loci in which one or multiple genes were significant in any of the gene-based tests (MAGMA, MetaXcan, SMR, and FOCUS) used in this study. "Multiple" indicates that several genes are involved. These genes have been named in Supplementary Data 14. h This column shows the loci for which the reported nearest genes or the significant genes identified in the gene-based tests above are members of at least one significant POAG pathway identified in this study. i "Multiple" indicates several genes.
variant and other gene-based tests performed in this study and were not previously reported for POAG (Supplementary Data 17). In a post hoc analysis looking solely at retina, two additional genes (CNTF and MPHOSPH9) were significant (given Bonferroni correction threshold for 6508 genes in retinal tissue). Additionally, we integrated our GWAS meta-analysis summary statistics with eQTL data from blood (CAGE eQTL summary data, N = 2765) and retina (EyeGEx eQTL data, N = 406) using SMR. Given that these eQTL data were obtained from people of European descent, we restricted this analysis to our European meta-analysis to ensure that different gene expression and LD structure patterns between ancestries did not influence the SMR findings. In retina and blood, 16 genes passed the SMR significance threshold corrected for the maximum number of 8516 genes tested in two tissues (P < 2.9e-06), of which eight had a P > 0.05 in the heterogeneity in dependent instruments (HEIDI) test 37 implemented in SMR, suggesting that the same association signals drive both gene expression and POAG risk, at these loci (Supplementary Data 18). Although the majority of the risk loci identified through the MetaXcan and SMR approaches were also identified in the meta-analysis, these analyses help with prioritizing the most likely functionally relevant genes. To further identify the most plausible causal genes based on gene expression data, we used the approach implemented in FOCUS v0.5 40 , a probabilistic framework that assigns a posterior probability for each gene causally driving TWAS associations in multiple tissues (Supplementary Data 19 summarizes the genes with a posterior probability >0.6).
Integrating data from several lines of evidence described above, as well as the cross-ancestry fine-mapping and genetic pathways, provided support for specific genes potentially influencing POAG risk particularly RERE, VCAM1, ZNF638, CLIC5, SLC2A12, YAP1, MXRA5, and SMAD6 (Table 1 and Supplementary Data 14). For example, rs3777588, a lead GWAS SNP in this study, is an intronic variant within CLIC5. The CADD score for this variant is 16.58, providing support for the pathogenicity of this variant (determined by a CADD score > 12.37). The lead SNP is also an eQLTL for CLIC5 in both GTEx and retina, and genebased analysis by incorporating eQTL data supported the involvement of CLIC5 in POAG risk. This approach also helped to select best genes near the lead SNPs or to shift the focus from the nearest genes to genes further away. For example, rs12846405, a lead GWAS SNP in this study, is an intergenic variant located between MXRA5 and PRKX. Based on integration of eQTL data and gene-based analysis, MXRA5 was prioritized as the most likely causal gene in this locus.
Differential expression of the previously unknown risk loci in eye tissues. We investigated the expression of the genes nearest to the lead SNP for the novel loci in 21 healthy eye tissues 11 (Supplementary Data 20 and Supplementary Fig. 7). Clustering analysis shows that the majority of these genes were expressed in eye tissues ( Supplementary Fig. 7A). We examined the differential expression of the previously unreported genes in ocular tissues likely to be involved in POAG pathogenesis, namely trabecular meshwork, ciliary body and optic nerve head, and found 36/51 (71%) of the genes differentially expressed in these tissues compared to the other eye tissues tested in this study (Supplementary Data 20 and Supplementary Fig. 7B).
Drug targets. At least 16 of the POAG risk genes (nearest to the lead SNPs) are targeted by existing drugs, some of which are already in use/clinical trials for several eye or systemic diseases (Supplementary Data 21). The functional relevance of 14 of these 16 drug target genes is supported by the bioinformatic functional analyses we used in this study (i.e., eQTL, chromatin interaction, etc; Supplementary Data 21). We discuss the relevance of some of these drugs in the discussion section below.
Sex-stratified meta-analysis. We identified a very high genetic correlation (rg = 0.99, se = 0.06) between POAG in men versus women (European stage 1 and UKBB self-reports combined). We also performed cross-ancestry, sex-stratified meta-analyses using a subset of the overall study with sex-stratified GWAS available (Supplementary Data 1; Supplementary Data 22 and 23; Supplementary Fig. 1C, D). Only one signal near DNAH6 appeared to have a female-specific effect (2:84828363[CA], OR = 1.6, P = 3.28e-09 for women; OR = 1.05, P = 0.56 for men).
Subtype-stratified meta-analysis. Based on IOP levels, POAG can be classified into two major subtypes: high-tension glaucoma (HTG) in which IOP is increased (>21 mmHg), and normal tension glaucoma (NTG) in which IOP remains within the normal range. We performed cross-ancestry subtype-specific metaanalyses using 3247 cases and 47,997 controls for NTG (Normal tension glaucoma defined as glaucoma with IOP < 21 mmHg), and 5144 cases and 47,997 controls for HTG (high-tension glaucoma with IOP > 21 mmHg (Supplementary Data 24 and 25 and Supplementary Fig. 1E, F). All NTG and HTG loci were also significantly associated with the overall POAG meta-analysis except for one locus near FLNB that was significant for NTG (lead SNP rs12494328[A], OR = 1.18, P = 1.7e-08), but did not reach the significance threshold for POAG overall (P = 7.5e-07). However, this SNP was significant in the 23andMe replication study (rs12494328[A], OR = 1.06, P = 1.35-e12), and has previously been associated with optic nerve head changes 41 . Overall, all NTG loci were at least nominally associated (P < 0.05) with HTG (and vice versa) except for rs1812974 (top SNP near ARHGEF12). Although this SNP had the same direction of effect for NTG, the effect was significantly larger for HTG than NTG (P = 0.007). Similarly, several other loci had significantly larger effects on one subtype (e.g., CDKN2B-AS1, FLNB, and C14orf39 had larger effects on NTG than HTG) (Supplementary Data 24 and 25). Overall, the genetic correlation between NTG and HTG was estimated to be 0.58 (se = 0.08) using LD score regression and the meta-analysis summary data from Europeans.

Discussion
In this large multi-ethnic meta-analysis for POAG, we identified 127 risk loci for POAG, of which 44 were not previously identified. We also identified additional risk loci using gene-based tests and highlighted genetic pathways involved in the pathogenesis of POAG. We observed relatively consistent genetic effects for POAG across ancestries. The risk loci include genes that are highly expressed in relevant eye tissues, nerves, arteries, as well as tissues enriched with these components. Functional relevance of the identified risk loci were further supported by eQTL and chromatin interaction data.
We identified a significant correlation between the POAG effect sizes of genome-wide significant SNPs, as well as all the SNPs throughout the genome, across Europeans, Asians, and Africans. Although previous studies have suggested that the genetic architecture of POAG might differ between Africans and Europeans 42 , we observed a moderate correlation (r~0. 45) between effect sizes of the POAG risk loci in Europeans and Africans ( Supplementary Fig. 2C), and the correlation was higher between Europeans and Asians (r~0.7). Although the overall correlation is moderately high across ancestries, there are genomic regions where the LD pattern differs by ancestry and our fine-mapping approach showed that incorporating GWAS data across ancestries improved the probability of finding a causal variant for 18 loci in this study, including known (e.g., AFAP1 and RELN loci) and unknown (e.g., GJA1/HSF2 and SEPT7) loci. However, the most probable causal variants in Europeans remained the same for~60% (75 out of 127) of the risk loci even after incorporating Asian and African GWASs. Overall, due to the relatively lower statistical power of our African studies, the fine-mapping results in this study were not strongly influenced by African GWASs, emphasizing that larger African POAG GWASs are required for better cross-ancestry fine-mapping in the future.
We also identified a previously unknown association of a human leukocyte antigen (HLA-G/HLA-H) with POAG. The HLA system is a gene complex encoding the major histocompatibility complex proteins in humans. These cell-surface proteins are responsible for the regulation of the human immune system. The most significant SNP in this region (rs407238) has been associated at the genome-wide significance level for other traits such as Celiac disease, intestinal malabsorption, disorders of iron metabolism, multiple blood traits, hyperthyroidism, multiple sclerosis, hip circumference, and weight (https://genetics. opentargets.org/variant/6_29839124_C_G). The mechanism of action of the lead SNP appears to be via IOP (UKBB IOP GWAS P = 8.8e-06).
The gene-sets enriched for the risk loci identified in this study indicate two major pathogenic mechanisms for POAG: (1) vascular system defects, mainly the molecular mechanisms contributing to blood vessel morphogenesis, vasculature development, and regulation of endothelial cell proliferation, and; (2) lipid binding and transportation-mainly the molecular mechanisms involved in intracellular lipid transport, apolipoprotein binding, negative regulation of lipid storage, and positive regulation of cholesterol efflux. Involvement of the vascular system in the pathogenesis of POAG is further supported by our results showing enrichment of the expression of the POAG risk genes in arteries and vessels. Molecular targets in these pathways can be potential candidates for treatment of POAG. Two of the significant genesets in this study (phagocytosis engulfment and negative regulation of macrophage derived foam cell differentiation) suggest an important role of immune system defects in increasing the risk of POAG.
Integrating several lines of genetic evidence provided support for specific genes within the identified risk loci that could influence risk through known and unknown processes. MXRA5 and SMAD6 are both involved in transforming growth factor (TGF) beta-mediated extracelluar matrix remodeling 43,44 , a process known to contribute to POAG risk 45 . Additionally, a SVEP1 missense allele was associated with POAG risk (rs61751937). SVEP1 encodes an extracellular matrix protein that is essential for lymphangiogenesis in mice, through interaction with ANGPT2 (the product of another POAG risk gene identified in this study), and modulation of expression of TEK and FOXC2 in knockout mice 46 . Lyphangiogenesis has an important role in the development of Schlemm's canal required for outflow of fluid from the eye 47,48 , and two other genes necessary for lyphangiogenesis and Schlemm's canal development (TEK, ANGPT1) cause childhood glaucoma 49,50 . Interestingly, SVEP1 was shown to be a modifier of TEK-related primary congenital glaucoma 51 . VCAM1 is an extracellular matrix cell adhesion molecule involved in angiogenesis and possibly regulation of fluid flow from the eye 52 . RERE mutations are a cause of neurodevelopmental disorders that can involve the eye 53 , providing further evidence for a role of ocular development in adult glaucoma 54 . While RERE has also been associated with VCDR 55 , this is the first association with POAG.
Genes involved in biological processes not previously known to contribute to glaucoma have also been implicated by this study. CLIC5 encodes a chloride channel that functions in mitochondria 56 and could have a role in ocular fluid dynamics. ZNF638 is a zinc finger protein that regulates adipose differentiation 57 and has been implicated in the genetic regulation of height 58 . SLC2A12 is a glucose transporter that is also involved in fat metabolism 59 . YAP1 is an oncogene that is a main effector of the HIPPO tumor suppressor pathway and apoptosis inhibitor 60 , processes that could influence retinal ganglion cell survival in glaucoma. In mice, heterozygous deletion of Yap1 leads to complex ocular abnormalities, including microphthalmia, corneal fibrosis, anterior segment dysgenesis, and cataract 61 .
Several proteins encoded by genes within the identified POAG risk loci are targets of some currently approved drugs. For instance, COL4A1 is targeted by ocriplasmin, a collagen hydrolytic enzyme that is currently used to treat vitreomacular adhesion (adherence of vitreous to retina). This drug can degrade the structural proteins including those located at the vitreoretinal surface 62 . Clinical trials are in progress to evaluate ocriplasmin therapy for several eye conditions including macular degeneration, diabetic macular edema, macular hole, and retinal vein occlusion. Some other drug candidates targeting proteins encoded by POAG-associated loci are also currently under consideration for treating dementia and cardiovascular diseases including acitretin, a retinoid receptor agonist targeting RARB, which has been considered for treatment of AD in ongoing clinical trials. Also, dipyridamole, a 3′,5′-cyclic phosphodiesterase inhibitor, targets PDE7B and current clinical trials are testing therapies based on dipyridamole for diseases such as stroke, coronary heart disease, ischemia reperfusion injury, and internal carotid artery stenosis. Given that our pathway analyses highlighted the involvement of vasculature development and blood vessel morphogenesis in POAG pathogenesis, dipyridamole could be a potential therapy for POAG through modulation of blood flow, which can be defective in POAG 63,64 . Further studies to confirm the functionality of these POAG risk genes in vivo and in vitro may support the suitability of repurposing these drugs as alternative treatments for POAG. Moreover, comprehensive fine-mapping is required to identify the most likely causal genes that can be targeted by currently approved drugs.
This study has several strengths and limitations. The main strength includes identification of risk loci contributing to the development of POAG across ethnic groups, an advance over prior POAG GWAS that have mainly focused on individuals from a single ancestry group. We showed that combining GWAS data across ancestries increases the power of gene mapping for POAG. Another strength is the integration of GWAS, gene expression, and chromatin interaction data to investigate the functional relevance of the identified loci, as well as to identify the most plausible risk genes.
A limitation of this study is that although the majority of the cases were clinically confirmed POAG, our data included >7200 glaucoma cases from the UKBB obtained through self-reports. However, we observed a very high concordance between the GWAS results for clinically validated cases versus self-report. Additionally, the vast majority of our results replicated in selfreport data from 23andMe. The second limitation of this study is its relatively low statistical power for the subtype-specific analyses (especially for the NTG subset), limiting the ability of this study to identify subtype-specific loci. Larger NTG GWASs are required to dissect the genetic heterogeneity between POAG subtypes. Third, although where possible each participating study adjusted for the effect of age in their association testing prior to the metaanalysis, in a subset of studies, cases and controls were not matched for age and future studies should fully investigate the effect of the identified risk loci across different age strata, particularly for loci where certain alleles are strongly associated with other age-related conditions. Finally, although we investigated the functional relevance of the identified risk loci using bioinformatic analyses, we did not confirm their functionality in vitro and in vivo. Further studies to investigate the biological roles of these risk loci with respect to POAG pathogenesis in relevant eye tissues will further shed light on the molecular etiology of POAG.
In conclusion, this study identified a strong cross-ancestry genetic correlation for POAG between Europeans, Asians, and Africans, and identified 127 genome-wide significant loci by combining GWAS results across these ancestries. The crossancestry data improved fine-mapping of causal variants. By integrating multiple lines of genetic evidence, we implicate previously unknown biological processes that might contribute to glaucoma pathogenesis including intracellular chloride channels, adipose metabolism and YAP/HIPPO signaling.

Methods
Study design and participants. We obtained 34,179 POAG cases and 349,321 controls including participants of European, Asian, and African descent from 21 independent studies across the world. Number of cases and controls, and distribution of age and sex for each study are summarized in Supplementary Data 1. The phenotype definition and additional details such as genotyping platforms for each study are provided in Supplementary Information. For most of the studies, we restricted glaucoma to POAG based on the ICD9/ICD10 criteria. However, considering that POAG constitutes the majority of glaucoma cases in Europeans 65 , we also included 7286 glaucoma self-reports from UK Biobank to replicate findings from the ICD9/ICD10 POAG meta-analysis in Europeans and to maximize the statistical power of the final stage meta-analysis (please see below). Informed consent was obtained from all the participants, and ethics approval was obtained from the ethics committee of all the participating institutions.
We performed a four-stage meta-analysis to combine GWAS data from the participating studies. In the first stage, we conducted a meta-analysis of the POAG GWAS in Europeans (16,677 POAG cases and 199,580 controls). In the second stage, we performed independent meta-analyses of POAG GWAS in Asians (6935 cases and 39,588 controls) and in Africans (3281 cases and 2791 controls) (Supplementary Data 1). As part of the second stage, the Asian and African metadata, as well as data from a GWAS of 7286 self-report glaucoma cases and 107,362 controls of European descent from UKBB were used to validate the findings from the European Caucasian meta-analysis. The UKBB self-report GWAS was completely independent of the UKBB IC9/ICD10 POAG GWAS; all the UKBB POAG cases and controls from the first stage, as well as their relatives (Pi hat >0.2), were removed from the self-report GWAS dataset. In the third stage, we combined the results from stage 1 and 2 to increase our statistical power to identify POAG risk loci across ancestries. In the fourth stage we replicated the stage 3 findings in a dataset from 23andMe.
To investigate sex-specific loci for POAG, we also conducted a meta-analysis of POAG in males and females separately. For this analysis, we had GWAS data from a subset of the overall POAG meta-analysis, including 10,775 cases and 123,644 controls for males, and 10,977 cases and 144,606 controls for females (Supplementary Data 1). Similarly, to identify risk loci for the HTG and NTG subtypes, we performed a subtype-specific meta-analysis using 3247 NTG cases and 47,997 controls, and 5144 HTG cases and 47,997 controls.
Quality control (QC) and imputation. Study-specific QC and imputation details have been provided in Supplementary Information. Overall, SNPs with >5% missing genotypes, minor allele frequency (MAF) < 0.01, and evidence of significant deviations from Hardy-Weinberg equilibrium (HWE) were excluded. In addition, individuals with >5% missing genotypes, one of each pair of related individuals (detected based on a p-hat > 0.2 from identity by descent calculated from autosomal markers), and ancestry outliers from each study (detected based on principal component analysis including study participants and reference samples of known ancestry) were excluded from further analysis (for more details please see Supplementary Information).
Imputation for studies involving participants of European descent was performed in Minimac3 using the Haplotype Reference Consortium (HRC) r1.1 as reference panel through the Michigan Imputation Server 66 . However, for a study of Finnish population (FinnGen study), whole-genome sequence data from 3775 Finnish samples were used as reference panel for better population-specific haplotype matching, which results in a more accurate imputation. For studies involving Asian and African participants, imputation was performed using the 1000 Genomes samples of relevant ancestry. SNPs with MAF > 0.001 and imputation quality scores (INFO or r 2 ) > 0.3 were taken forward for association analysis.
Association testing. Association testing was performed assuming an additive genetic model using dosage scores from imputation, adjusting for age, sex, and study-specific principal components as covariates, using software such as PLINK2 66,67 , SNPTEST v2.5.1 68,69 ,SAIGE v0.36.3 70 , EPACTS v3.2.6 (https:// genome.sph.umich.edu/wiki/EPACTS), and Rvtests v2.0.3 71 . For studies with a large number of related individuals, mixed-model association testing was performed to account for relatedness between people. For the X chromosome analysis, we used the following approach to allow for dosage compensation: females were coded as 0 (homozygous for non-effect allele), 1 (heterozygous), and 2 (homozygous for effect allele) while males were coded as 0 (no effect allele) and 2 (one effect allele). The covariates were the same as for the association testing for the autosomes except removing sex.
To confirm the validity of combining GWAS results across populations comprising different ancestries, we estimated the genome-wide genetic correlation for POAG between the populations of European, Asian, and African descent participating in this study. For this purpose, we used Popcorn v0.9.9 21 , a toolset that provides estimates of genetic correlation while accounting for different genetic effects and LD structure between ancestries. For this analysis, the LD scores for each ancestry population were estimated using the 1000G populations (Europeans, Asians, and Africans), and SNPs were filtered based on the default MAF = 0.05 in Popcorn.
We performed within and between ancestry meta-analyses using a fixed-effects inverse-variance weighting approach in METAL (the version released on 2011-03-25) 72 using SNP effect point estimates and their standard errors. The presence of heterogeneity between SNP effect estimates across studies were investigated using the Cochran's Q test implemented in METAL. To identify multiple independent risk variants within the same locus using GWAS summary statistics obtained from the meta-analysis, we used the Conditional and Joint (COJO) analysis implemented in GCTA v1.26 73 . Q-Q and Manhattan plots were created in R-3.2.2, and regional association plots in LocusZoom v1.4 74 .
We used the univariate LD score regression 19 intercept for each study separately as well as for the meta-analyzed results to ensure that the test statistics did not include model or structural biases such as population stratification, cryptic relatedness, and model misspecification. To investigate the genetic correlation between POAG and AD, we used bivariate LD score regression 32 using two large AD GWAS meta-analyses 30,31 . To investigate the genetic correlation between POAG and the other traits, we used bivariate LD score regression through the LD Hub platform (http://ldsc.broadinstitute.org/ldhub/).
The association of the POAG risk loci identified in this study with its major endophenotypes, IOP and VCDR, was investigated using summary statistics from a recent GWAS meta-analysis for IOP (N = 133,492) 11 and VCDR (N = 90,939) 15 .
The variance in the POAG familial risk explained by the loci identified in this study (N = 127) was calculated based on where p i and β i refer to the MAF and the magnitude of association of the i-th SNP, respectively, and log (λP) is the familial relative risk obtained from observational studies. The estimates for p i and β i in this study were obtained from UKBB and European POAG metaanalysis, respectively, and log (λP) was 9.2 estimated in a previous study 75 .
23andMe replication. We validated the genome-wide significant risk loci from our cross-ancestry meta-analysis (127 independent SNPs) and subtype analyses (7 independent SNPs) in a subset of 23andMe research participants of European descent comprising 43,254 POAG cases and 1,471,118 controls. POAG cases were defined as those who reported glaucoma, excluding those who reported angleclosure glaucoma or other types of glaucoma. Controls did not report any glaucoma. Association testing was performed using logistic regression assuming an additive genetic model, adjusting for age, sex, top five principal components, and genotyping platform as covariates.
Cross-ancestry fine-mapping. We used PAINTOR v3.0 76,77 to perform a crossancestry fine-mapping for the 127 risk loci identified in this study. For this analysis, the GWAS summary statistics for 1 Mb either side of the lead risk SNPs were extracted from European (including UKBB self-reports), Asian, and African metaanalyses, separately. To account for different LD patterns between ancestries, we created ancestry-specific LD matrices between SNPs using 1000G phase 3 as a reference panel. We allowed for the presence of two causal SNPs per locus. To investigate any advantage of fine-mapping across ancestries, we compared the posterior probabilities of the prioritized causal SNPs in Europeans separately, as well as across ancestries, without including any annotation data.
Gene-based and pathway-based tests. Gene-based and gene-set (pathway) based tests were performed using the approach implemented in MAGMA v1.07b 33 . We performed this analysis for each ancestry separately, and the P-values were then combined across ancestries using Fisher's combined probability test. The significance threshold for gene-based test was set to P < 0.05/20174, and for pathwaybased tests to P < 0.05/10678, accounting for the maximum number of independent genes/pathways tested. In addition, MAGMA was used to investigate the enrichment of the expression of the significant risk genes in GETX v6 tissues (P < 1e-03 accounting for 53 tissues tested).
To identify loci with effect on POAG risk due to modulation of gene expression, we also used alternative gene-based tests that integrate GWAS summary statistics with eQTL data throughout the genome (TWAS-based approaches). For this purpose, we used MetaXcan v0.3.3 38 , SMR v0.69 37 , and FOCUS v0.5 40 . MetaXcan uses GWAS summary statistics to impute the genetic component of gene expression in different tissues based on a reference eQTL panel. We used the EyeGEx eQTL data from retina 39 , as well as 44 GTEX tissues as reference eQTL data for this study. The GWAS input for this analysis included the summary statistics obtained from the cross-ancestry meta-analysis. We set the significance threshold to Bonferroni-corrected P < 1.5e-07, accounting for the maximum number of genes tested (N = 7209) in 45 tissues. To investigate the probability of each gene causally driving TWAS associations, we used the approach implemented in FOCUS 40 , a probabilistic framework that assigns a posterior probability to each gene. We used the gene expression reference weights provided in FOCUS, which combines the expression weights obtained from GTEX tissues with the weights provided in FUSION 78 obtained from several tissues, including adipose, peripheral blood, whole blood, and brain.
SMR uses a Mendelian Randomization framework to identify genes whose expression is likely modulated through the same variants associated with the outcome of interest (POAG). For the SMR analysis, we used the following eQTL data: CAGE eQTL summary data from blood (N = 2765) and EyeGEx eQTL data from retina (N = 406). The SMR significance threshold was set to P < 2.9e-06, accounting for the maximum number of 8516 genes tested in two tissues. A heterogeneity P > 0.05 from the HEIDI test implemented in SMR implies that we cannot reject the null hypothesis that a single causal variant is likely to affect both gene expression and POAG risk for these loci.
Gene expression. RNA was extracted from the corneal layers, trabecular meshwork, ciliary body, retina, and optic nerve tissues from 21 healthy donor eyes of 21 individuals. After quality control, the RNA was sequenced using Illumina NextSeq® 500 (San Diego, USA). Trimgalore (v0.4.0) was used to trim low-quality bases (Phred score < 28) and reads shorter than 20 bases after trimming were discarded. Data analysis was done with edgeR (version 3.22.5) 79 . Only genes expressed a minimum of 10 times (1.5 counts per million) in at least five dissected tissues were kept. The RNA count libraries were normalized using trimmed mean of M-values method 80 . Two-group differential expression analysis was done via negative binomial generalized linear model in edgeR 81 . The RNA expression in ciliary body, trabecular meshwork, and optic nerve head, which are involved in aqueous production, drainage and principal site of glaucoma injury, respectively, was compared to the remaining eye tissues.
Drug targets. We used Open Targets 82 to search for drugs currently in use or in clinical trials for treating other ocular or systemic diseases that target the POAG risk genes identified in this study. These drugs can be potentially repurposed as alternative treatments for POAG, owing to in vivo and in vitro confirmation of the functionality of the target genes in the pathogenesis of POAG.
Bioinformatic functional analyses. The bioinformatic functional analysis to investigate the functional relevance of the identified risk loci for POAG were performed through the FUMA platform v1.3.5 83 using the following dataset/toolsets: GTEX eQTL v6 84 ; Blood eQTL browser 85  Colocalization analysis. To test whether the POAG loci tag a shared causal variant or haplotype with Alzheimer's disease (AD), given that several loci overlap, we applied the bayesian-based colocalization method, eCAVIAR v2.0.0 29 to all 123 and 66 autosomal POAG significant variants from the cross-ancestry and European subset POAG GWAS meta-analyses, respectively, and the AD GWAS metaanalysis from Kunkle et al. 30,31 and Jansen et al. 31 that lack summary statistics for chromosome X. Loci were defined by identifying the outermost variant on either side of each lead POAG variant corresponding to r 2 > 0.1 relative to the POAG variant with an added 50 kb on either side. We used the default parameters of eCAVIAR (http://genetics.cs.ucla.edu/caviar/manual.html), assuming up to two causal variants per locus. We used the European derived samples from 1000 Genome Project Phase 3 as the reference panel to compute the LD matrix between all variant pairs within each locus needed in eCAVIAR, since both AD GWAS meta-analyses consisted of samples from European ancestry. We used a colocalization posterior probability (CLPP) above 0.01 as the significance cutoff as recommended by the tool 29 .