Glaucoma Multi-trait genome-wide association study identifies new loci associated with optic

A new avenue of mining published genome-wide association studies includes the joint analysis of related traits. The power of this approach depends on the genetic correlation of traits, which re ﬂ ects the number of pleiotropic loci, i.e. genetic loci in ﬂ uencing multiple traits. Here, we applied new meta-analyses of optic nerve head (ONH) related traits implicated in primary open-angle glaucoma (POAG); intraocular pressure and central corneal thickness using Haplotype reference consortium imputations. We performed a multi-trait analysis of ONH parameters cup area, disc area and vertical cup-disc ratio. We uncover new variants; rs11158547 in PPP1R36-PLEKHG3 and rs1028727 near SERPINE3 at genome-wide signi ﬁ cance that replicate in independent Asian cohorts imputed to 1000 Genomes. At this point, validation of these variants in POAG cohorts is hampered by the high degree of heterogeneity. Our results show that multi-trait analysis is a valid approach to identify novel pleiotropic variants for ONH. 89 strategies. Prioritized genes for ONH traits were highly overlapping and were combined to form a set of 295 unique genes for further functional annotation in FUMA. For IOP and CCT 11 and 116 genes were prioritized respectively. We further investigated the FUMA-mapped genes for enrichment using hypergeo-metric enrichment tests on pre-de ﬁ ned gene-sets derived from MsigDB and WikiPathways. p Values were corrected based on Bonferroni method for the number of tested gene-sets. Statistics and reproducibility . Software used for the data analysis of this study: METAL (https://genome.sph.umich.edu/wiki/METAL), EasyQC (www.genepi-regensburg.de/easyqc), GCTA (http://cnsgenomics.com/software/gcta/), FUMA (https://fuma.ctglab.nl/), LDSC (https://github.com/bulik/ldsc), CPASSOC (http:// hal.case.edu/zhu-web/), and MTAG (https://github.com/omeed-maghzian/mtag). In the single trait analyses and the multi-trait analyses variants surpassing a p value less than 5× 10 − 8 were considered genome-wide signi ﬁ cant and followed up for replication in the Asian replication sample. Replication was de ﬁ ned as variants surpassing nominal signi ﬁ cance level p value< 0.05. Variants considered when located >1 Mb from a previously reported variant.

G laucoma is the most common cause of irreversible blindness in the world 1 . Primary open angle glaucoma (POAG) is the most prevalent type of glaucoma accounting for 74% of all glaucoma cases 2,3 . Intraocular pressure (IOP) and the morphology of the optic nerve head (cup area (CA), disc area (DA), and vertical cup-disc ratio (VCDR)) are important features of the glaucomatous process. For each of these traits, twin studies showed a high heritability (h 2 CA = 0.75, h 2 DA = 0.72, h 2 IOP = 0.55, and h 2 VCDR = 0.48) 4 . Central corneal thickness (CCT) is also a highly heritable trait (h 2 CCT = 0.68-0.95) 5 , which is most likely non-physiologically associated with POAG, but rather biases IOP measurement, the major risk factor of POAG 6,7 . CA, DA, and VCDR, are significantly correlated both at the genetic level 8 . The high genetic correlation found between the optic nerve head (ONH) traits (R g = 0.31-0.83) raises the question whether multi-trait analyses will improve the statistical power of the individual GWAS and will find variants with pleiotropic effects 9 .
For this study, we generated new data on these 5 quantitative traits by imputing 12 European ancestry cohorts from the International Glaucoma Genetic Consortium (IGGC) (n MAX = 31,269) to haplotype reference consortium (HRC) release 1 imputation panel, which includes over 39 million variants 10 . A meta-analysis of these 12 European ancestry studies served as a discovery cohort in the analyses. Replication was performed in five Asian ancestry cohorts that were part of the IGGC. The cohorts of Asian descent were imputed to 1000 Genomes as there is little gain in HRC imputation in this ancestry group because there are no additional Asian samples included in HRC (http://www.haplotype-referenceconsortium.org/participating-cohorts) 11 . We evaluated the added value of multi-trait analyses using two programs: CPASSOC and multi-trait analysis of GWAS (MTAG). Both use aggregated GWAS results. Whereas CPASSOS performs a meta-analysis assuming homogeneous and heterogeneous effects across traits by applying a inter-trait correlation matrix, MTAG basically increases the power of a single trait analyses by incorporating the GWAS findings of correlated traits.
By multi-trait analysis we identified two novel loci associated with the ONH at rs11158547 in-between PPP1R36 and PLEKHG3 and at rs1028727 near SERPINE3 in those of European descent. These loci replicated in the Asian replication sample. Findings for these loci were consistent using a distinct multi-trait approach, MTAG, in both the European and Asian cohorts. This study emphasis that multi-trait analysis in GWAS pleiotropic traits is an effective approach to identify variants harboring correlated traits

Results
Replication of previous CA, DA, VCDR, IOP, and CCT GWAS results. As a validation we first confirmed previously identified loci for CA, DA, VCDR, IOP, and CCT by Springelkamp et al. 8 (n CA = 22,489; n DA = 22,504; n VCDR = 23,899; n IOP = 37,930) and Iglesias et al. 12 (n = 17,803) based on 1000 Genomes imputation. Supplementary Fig. 1 and Supplementary Datas 1 and 2 show the per trait comparison of our meta-analysis of all European ancestry discovery cohorts using the HRC imputation with the results of the meta-analysis by Springelkamp 8 and Iglesias 12 based on 1000 Genomes. Out of 113 (95%), 107 available variants in HRC replicated at a Bonferroni significance level.
IOP and CCT. Next, we conducted a single trait meta-analysis for IOP and for CCT, the two traits that are not likely physiologically related. For IOP, we meta-analyzed a total of 31,269 participants (LDSC intercept = 1.028; SE = 0.0078, Supplementary Data 4) and identified 9 genome-wide significant regions of which two were novel in the HRC-based imputations and had not been uncovered in the IGGC 1000 Genomes analyses before (Table 1, Supplementary Data 5). The lead single-nucleotide polymorphisms (SNPs) in these genomic regions were a common variant rs9853115 near DGKG on 3q27.3 and a rare variant rs150202082 [T] (frequency 0.03) near TCF4 on 18q21.2. rs9853115 failed replication in the Asians (p = 0.9315) and rs150202082 could not be examined since this variant was monoallelic in the Asian individuals. A GWAS by Choquet et al. 15 also identified rs9853115 as new variant associated with IOP in multiethnic cohort of predominately (83%) European ancestry. The same study also identified a novel variant near TCF4, rs11659764, approximately 300 kb upstream of rs150202082 which was in relatively weak LD (r 2 = 0.4). In a recent study from the UKbiobank by Khawaja et al. 16 the same variant near DGKG showed genome-wide significant association with similar effectsize, however, rs150202082 near TCF4 could not be validated in this study.
In the meta-analysis of CCT, a total of 16,204 participants were included (LDSC intercept = 0.989; SE = 0.0082). We identified 31 independent genome-wide significant signals of which three were novel (Table 1 and Supplementary Data 5), including a variant, rs34869, near CDO1. Again, none of the three new variants replicated at a nominal significance level in the Asian samples. Multi-trait analysis by CPASSOC identified four novel variants (Supplementary Datas 8 and 9). In contrast to ONH cross trait analyses, also these could not be replicated in the Asian.
In silico analysis. To investigate the functional and regulatory potential, we annotated the variants in linkage disequilibrium (European LD, r 2 ≥ 0.8) with the lead SNPs at the two new and replicated ONH variants, rs11158547 and rs1028727, using a combination of bioinformatics tools (see Method section). A total of 70 variants in LD with the 2 novel variants were queried. None of the examined variants were predicted to damage protein structure by SIFT, Polyphen, or alternative splicing using Ensembl's Variant Effect Predictor. As all queried variants are noncoding, we reviewed the possible regulatory annotation of these SNPs in experimental epigenetic evidence, including DNAse hypersensitive sites, histone modifications, and transcription factor-binding sites in human cell lines and tissues from the ENCODE 17 and ROADMAP EPIGENOMICS 18 projects, integrated in Haploreg 19 . Annotations of chromatin states indicated that the two novel variants were located in, or in LD with, an active chromatin state region from at least one of the tissues investigated (Supplementary Data 10 and Figs. 3 and 4 for chromatin states in brain tissue). Next, we evaluated the overlaps   of cis-expression quantitative trait loci (eQTL) in several databases (see Methods). In both novel loci PPP1R36-PLEKHG3 and SERPINE3, variants were found to be eQTL's and based on Regulome DB-scores both ONH loci contained variants that were likely to alter binding (Supplementary Data 10).
Gene prioritization, pathway analysis, and gene expression. We explored possible tissue expression and biological functions by pathway analysis for the two novel SNPs. We annotated these SNPs to genes by positional gene mapping, eQTL mapping and chromatin interaction mapping strategies implemented in FUMA 20 (see Method section). For, rs11158547 (PPP1R36-PLEKHG3) 9 genes were assigned to this locus and for rs1028727 (SERPINE3) 21 genes were mapped to this locus (Supplementary Data 11). Pathway analysis based on enrichment of gene-set terms (MsigDB 21 and Wikipathways 22 ) found 5 and 11 Bonferroni significant gene-sets comprising genes mapped to SERPINE3 locus and PPP1R36-PLEKHG3 locus respectively. These pathways were in particular involved in immune response and cancer development (Supplementary Data 12 highlighted gene-sets).
As expression in eye tissues is not available in GTex 23 , we assessed the Ocular Tissue Database 24 . For 26 out of the 30 genes mapped to either rs11158547 or rs1028727 expression data were present in the Ocular Tissue Database (Supplementary Data 11). The highest levels of expression in the optic nerve was found for HSPA2, a gene associated to rs11158547 via lung, tibial nerve and fibroblasts eQTL's and chromatin interaction mapping (Supplementary Fig. 3).
From endophenotypes to glaucoma. We also investigated the translational potential of these two loci to POAG by carrying out a meta-analysis of three POAG studies, NEIGHBOR, Southampton and UK Biobank Eye and Vision Consortium (Ncase = 9450; Ncontrol = 436,824), from European origin. For rs1028727 a negative effect on the VCDR in the present study predicts a decreased risk of POAG which was seen in NEIGH-BOR study and the UKBiobank but not in the Southampton study (Supplementary Data 13). rs11158547 in PPP1R36-PLEKHG3 is predicted to be associated with increased POAG risk based on the positive effect of VCDR. Indeed in all three POAG studies the effect of the SNP is also positive. Pooling the studies based on a fixed effect analysis yields an OR 1.28 (95% CI: 1.15-1.39; p = 4.83 × 10 −8 ) and showed Bonferroni significance (p = 0.025). Given the high degree of heterogeneity of effects at this locus a random effect meta-analysis was carried out which could not confirm this finding in POAG (Supplementary Data 14). Thus the findings were partly but not consistently replicated, awaiting larger and more homogeneous data sets for the final replication.

Discussion
Our results implicate two novel loci, one downstream SERPINE3 and one other in-between PPP1R36 and PLEKHG3, both associated with ONH morphology via CPASSOC multi-trait analysis. SERPINE3 belongs to the clade E family of extracellular serpins. Family members have been described to play a role in other neurodegenerative diseases such as Alzheimer's disease 25 . Recent studies in glaucomatous human postmortem samples and in rat models identified oxidative inactivation of serpins (neuroserpin) as a molecular mechanism of increased plasmin activity leading to neurodegeneration in high ocular pressure conditions 26 . In the trabecular meshwork, serpins (plasminogen activator inhibitor) may mediate the inhibition of matrix metaloprotienase (MMPs) activity induced by transforming Table 1 Genome-wide significant SNPs newly identified for cup area, disc area, vertical cup-disc ratio, intraocular pressure or central corneal thickness in the European HRC discovery.  growth factor-beta enhancement 27,28 . Inactivity of MMPs were found to increase aqueous humor outflow resistance leading to rising IOP 29 . rs11158547 downstream the PLEKHG3 gene, a pleckstrin homology domain containing protein, is also a relatively unknown gene. It contains a guanide nucleotide exchange factor (GEF) domain which is important for Rho-dependent signal transduction 30 . In mice PLEKHG3 knockout is associated with an abnormal anterior chamber depth of the eye (IMPC release 3.2 http://www.mousephenotype.org/data/experiments? geneAccession=MGI:2388284). The other gene close to this locus, PPP1R36, is less likely a candidate gene for ONH. PPP1R36, has been described in connection with autophagy during spermatogenesis. PPP1R36 encodes a regulatory subunit of protein phosphatase 1 which is involved in multiple cellular functions such as metabolism, immune response, apoptosis, meiosis, mitosis, cytoskeletal reorganization, and synthesis. However, its function to the eye has not been characterized. A general insight from this study was that the strength of genetic correlation is an important condition for investigating the pleiotropic effect of genetic variants on traits. The high genetic and phenotypic correlation observed between the ONH traits enabled multi-trait analyses and yielded plausible and replicable results. We validated these novel variants identified by CPASSOC with another multi-trait analysis method MTAG, which also uses summary statistics as input. The main difference between both methods is that CPASSOC test whether a SNP is not associated with any of the traits under the null hypothesis. MTAG, on the contrary produces trait-specific effect estimates for each SNP. The variants discovered in the European CPASSOC analysis that replicated in the Asian CPASSOC analysis also replicated in the MTAG analysis. Sensitivity analysis excluding the Rotterdam studies showed a high correlation of the Shom and Shet statistic with the SHom/SHet statistic from the full analysis including the Rotterdam studies (r = 0.71, p < 2.2 × 10 −16 ; r = 0.68, p < 2.2 × 10 −16 , respectively). The correlation in MTAG for CA, DA, and VCDR was considerably low yet very significant (r CA = 0. 31, p < 2.2 × 10 −16 ; r DA = 0.56, p < 2.2 × 10 −16 , r vcdr = 0.17).
In contrast, the multi-trait analysis between IOP and CCT could not uncover robust new variants. A reason for this observation might be the moderate magnitude of the genetic correlation between IOP and CCT. Also, clinical research has shown that this relation is largely driven by measurement errors in Goldmann applanation tonometry, rather than a pathophysiological process 31,32 .
A potential limitation of this study is the application of a different imputation panel for the discovery and replication phase. The European studies were all imputed to the HRC panel which has a beneficial imputation quality compared to 1000 Genomes. By contract the Asian replications set, was imputed to 1000 Genomes since a recent publication showed that HRC imputations perform less adequate in Asians 11 . Variants associated with cup area and other endophenotypes at genome-wide significance in the European single trait analysis could not be replicated in the Asian replication sample. The use of different imputation panels may be a source bias hampering replication. A theoretical shortcoming, is that the various studies used different methods and equipment to assess ONH parameters among studies. This has most likely reduced the power of the study and has generated most probably false negative rather than false-positive results. To prevent false-positive findings using novel methods, we aimed to replicate the findings of the primary CPASSOC analyses by another analyses using MTAG. The variants discovered in the CPASSOC analysis could also be replicated in MTAG. As both methods are mathematically distinct we concluded that our results are rather robust and independent of the statistical approach. This underscores the strength of the association as it is consistently found by two independent approaches.
We have no doubt that association of the variants to ONH are of interest to the biology community. To evaluate the implications of our findings in the context of glaucoma we studied three independent POAG studies. Up until now we are not able to link our finding to POAG. Although fixed effect meta-analysis showed Bonferroni significance (p = 0.025) for rs11158547 in PPP1R36-PLEKHG3, random effect meta-analysis that takes into account the heterogeneity could not confirm this finding in POAG. The source of the high variability in estimates is unknown and may involve clinical variability and ethnic differences. It is important to realize that the identification of POAG genes is far from complete and work in progress.
In conclusion, we conducted single and multi-trait meta-analysis of five endophenotypes of glaucoma based on HRC imputations in European ancestral populations. The HRC single trait analyses in those of European descent did not yield new loci that could be replicated in Asians. We identified two novel loci for ONH in-between PPP1R36-PLEKHG3 at chromosome 14q23.3 and near SERPINE3 at chromosome 13q14.3 by multi-trait analysis in those of European descent that could be replicated in Asians using CPASSOC. Findings for these loci were consistent using MTAG. The present study underscores that multi-trait analysis in GWAS of true pleiotropic traits in relatively small sample sizes is a powerful approach to identify variants harboring correlated traits. Although these novel loci could not be directly associated with POAG it is likely that the genes in these regions mediate the glaucomatous process by their effect on the optic nerve morphology. For instance the PLEKHG3 gene identified in this study is involved in the Rho signaling cascade, this pathway is   Fig. 3 Regional, chromatin state, and eQTL plot for rs11158547 (PPP1R36-PLEKHG3). Panel a shows the regional assocaiations plots with −log10 p value depicted on the y-axis, genes mapped by either position, eQTL or chromatin interaction are depicted in red on the x-axis; panel b shows 15 core chromatin states of varaints plotted in panel a for 13 brain tissues from Roadmap epidenomes described on the y-axis. E054 ganglion eminence derived primary cultured neurospheres, E053 cortex-derived primary cultured neurospheres, E071 brain hippocampus middle, E074 brain substantia nigra, E068 brain anterior caudate, E069 brain cingulate gyrus, E072 brain inferior temporal lobe, E067 brain angular gyrus, E073 brain dorsolateral prefrontal cortex, E070 brain germinal matrix, E082 fetal brain female, E081 fetal brain male, E125 NH-A astrocytes primary cells. Panel c depicts varaints that overlap eQTLs from selected eQTL databases described in the legend of panel (c).
known to play a crucial role in POAG pathophysiology and is currently targeted for new therapies for POAG 33 . Our bioinformatic analyses suggests that both the PPP1R36-PLEKHG3 and SERPINE3 variants are eQTL's opening avenues to counteract the problem by RNA interference. Further research including exome sequencing and functional studies are needed to further define these genes in the mechanism of POAG.

Methods
Study design. We performed a meta-analysis of European origin GWAS's imputed to HRC reference panel release 1. We analyzed five outcomes: CA, DA, VCDR, IOP, and CCT. The CA phenotype was adjusted for DA in all analyses since these phenotypes are clearly correlated (Pearson correlation coefficient 0.6). Subsequently, we performed multi-trait analysis for CA, DA, VCDR, IOP, and CCT. Replication was carried out for the single trait as well as the multi-trait analysis in a meta-analysis of 5 Asian cohorts imputed to 1000 genomes. We also tested significance of lead SNPs in three independent POAG cohorts.  Fig. 4 Regional, chromatin state and eQTL plot for rs1028727 (SERPINE3). Panel a shows the regional assocaiations plots with −log10 p value depicted on the y-axis, genes mapped by either position, eQTL or chromatin interaction are depicted in red on the x-axis; panel b shows 15 core chromatin states of varaints plotted in panel a for 13 brain tissues from Roadmap epidenomes described on the y-axis. E054 ganglion eminence derived primary cultured neurospheres, E053 cortex-derived primary cultured neurospheres, E071 brain hippocampus middle, E074 brain substantia nigra, E068 brain anterior caudate, E069 brain cingulate gyrus, E072 brain inferior temporal lobe, E067 brain angular gyrus, E073 brain dorsolateral prefrontal cortex, E070 brain germinal matrix, E082 fetal brain female, E081 fetal brain male, E125 NH-A astrocytes primary cells. Panel c depicts varaints that overlap eQTLs from selected eQTL databases described in the legend of panel (c). Association analysis in discovery cohorts. Within each discovery cohort, each genotyped or imputed variant was tested for association with each of the traits, assuming an additive genetic model. The measurements were adjusted for sex, age, and five principal components in all cohorts and if necessary also for cohortspecific covariates (Supplementary Table 8). Family-based studies were adjusted for family structure. Given the clear correlation of CA with DA (Pearson's correlation r = 0.59 in Rotterdam Study I), the CA GWAS was adjusted for DA in all discovery cohorts prior to meta-analysis. Linear regression was employed for studies with unrelated individuals, and linear mixed effects models were used to account for family structure in the family-based studies.
Centralized quality control. Before meta-analysis, a centralized quality control procedure implemented in EasyQC was applied to individual study association summary statistics to identify outlying studies 35 . We included variants with imputation quality ≥ 0.3 (e.g., Minimac R 2 ) and expected minor allele count >6.
Additional checks for quality control were applied on the already filtered datasets including review of P-Z-plots, allele frequency plots and calculation of genomic inflation factor λ.
Meta-analysis of discovery cohorts. The association results of all studies were combined in a fixed effect inverse variance meta-analysis in METAL 36 , since there was no sample overlap or cryptic relatedness as checked by LD score regression (see methods genetic overlap). This tool also applies genomic control by correcting the test statistics to account for small amounts of population stratification or unaccounted relatedness. We also assessed heterogeneity by calculating I 2 values and Cochrans Q-test for heterogeneity as implemented in METAL. After metaanalyses of all available variants, we excluded the variants that were not present in at least three studies. for IOP, and 1288 for CCT. Re-running the meta-analysis excluding TEST-BATS study to show that the significantly younger mean age in this study did not distorted our findings showed nearly perfect correlation between effect estimates from the full analysis and the effect estimates from in the analysis excluding TEST-BATS (CA r = 0.99; DA r = 0.99; VCDR = 0.99; IOP = 0.99). Furthermore, the mean differences between effect estimates found in the full analysis and the effect estimates found in the analysis excluding TEST-BATS were zero (CA mean difference = 0, SD = 0.00; DA mean difference = 0, SD = 0.00; VCDR mean difference = 0, SD = 0.01; IOP mean difference = 0, SD = 0.06), this also suggests that the younger age in the TEST-BATS study has not biased the results.
Selection of independent variants. We examined whether multiple independent variants at a given locus influenced a trait and if they were independent of previous findings, we used the genome-wide complex trait analysis software (GCTA) 37 . This tool performs a stepwise selection procedure to select multiple associated SNPs by a conditional and joint (-CoJo) analysis approach using summary-level statistics from a meta-analysis and LD corrections between SNPs. The three Rotterdam Study cohorts (N = 5815), imputed with the HRC reference panel version 1, were used as the reference to calculate the LD, because it represents the largest discovery studies. LD was calculated between pairwise SNPs, but any SNP further than 10 Mb apart were assumed to not be in LD. All autosomal chromosomes were analyzed, with MAF restricted to ≥0.01 estimated from the three Rotterdam Study cohorts. The independent variants were annotated by Haploreg 19 , see Supplementary  Table 9.
Identification of potential novel variants. Previously, Springelkamp et al. 8 Multi-trait analysis. For multi-trait genome-wide association analysis we applied the CPASSOC package developed by Zhu et al. 13 . We used CPASSOC for two analyses to combine the association results from CA, DA, VCDR, and from IOP and CCT. CPASSOC generates two statistics, SHom and SHet. SHom is similar to the fixed effect meta-analysis method but accounts for the correlation of summary statistics because of the correlated traits. SHom uses the sample size of a trait as a weight instead of variance, so that it is possible to combine traits with different measurement scales. SHet is an extension of SHom, but power can be improved when the genetic effect sizes are different for different traits. To compute statistics SHom and SHet, a correlation matrix is required to account for the correlation among traits or induced by overlapped or related samples from different cohorts.
We followed the approach previously described by Park et al. 42 , to calculate this correlation matrix. Briefly, we used all independent SNPs (r 2 < 0.2) present in datasets that were not associated with any of the traits (−1.96 > Z score < 1.96), and took the Pearson's correlation of their Z-scores 13 . For both tests QQ-plots were created ( Supplementary Fig. 8). Novel loci identified by CPASSOC (p < 5 × 10 −8 ) that were not implicated in the single-trait analysis were validated using a second multi-trait method, MTAG. Similarly, MTAG also utilizes summary statistics as input, but performs LD score regression to estimate the genotypic and phenotypic variance-covariance matrices. In contrast to CPASSOC, MTAG performs association tests for each individual trait by boosting the power of a signal and providing an estimation of the underlying association via the multi-trait variancecovariance structure. We applied MTAG to SNPs MAF > 0.01 for combing the analysis of CA, DA, and VCDR, and the analysis of IOP and CCT. For the European sample we used the 1000Genomes European pre-calculated LD scores and for the Asians the 1000Genomes East-Asian pre-calculated LD scores (https://data. broadinstitute.org/alkesgroup/LDSCORE/). We then validated each of the genomewide significant signals identified by CPASSOC in the MTAG results.
Replication in Asian cohorts imputed to 1000Genomes. All independent SNPs identified with p < 5 × 10 −8 in the discovery stage (single and multi-trait analysis) were carried forward for replication in Asians. For single-trait analyses, we validated these signals in fixed effect meta-analyses previously reported by The genetic overlap between CA, DA, VCDR, IOP, and CCT. To further investigate the genetic overlap among CA, DA, VCDR, CCT, and IOP we used the LD Score regression implemented in LDSC 44 to examine the pattern of genetic correlations. The LD score for each SNP measures the amount of pairwise LD(r 2 ) with other SNPs within 1-cM (centimorgan) windows based linkage disequilibrium. Bivariate LD score regression can estimate the extent to which two phenotypes share genetic variance. Summary statistics of the five meta-analysis were formatted to LDSC input files, we followed quality control as implemented by the LDSC software package (https:// github.com/bulik/ldsc). We used pre-calculated LD scores provided by the developers for each SNP using individuals of European ancestry from the 1000 Genomes project that are suitable for LD score analysis in European populations. SNP heritability estimates for all traits and genetic correlations were then calculated between the traits, see Supplementary Data 3 and Fig. 2.
Bioinformatical annotation. Using the software HaploReg (version 4.1) 19 and RegulomeDB v1.1 45 , we annotated the potential regulatory functions of the replicated GWAS SNPs and their proxies (r 2 > 0.8, 1000 genomes CEU) based on epigenetic signatures. We examined whether these variants (GWAS SNPs and variants in LD with the GWAS SNPs) overlapped with regulatory elements including DNAse hypersensitive sites, histone modifications, and transcription factor-binding sites in human cell lines and tissues from the ENCODE Project and the Epigenetic Roadmap Project. We then used the RegulomeDB score to assess their potential functional consequence, as described previously 46 .
Pathway analysis. We applied FUMA, which uses a three way gene-mapping strategy, to assign genome-wide significant SNPs to genes of interest. For positional mapping, SNPs in LD with the independent SNPs were mapped to genes using a window of 10 kb. eQTL mapping was performed by mapping SNPs to genes up to 1 Mb (cis-eQTL). eQTLs from all tissues available in GTEx v6 23 , Blood eQTL browser 47 , BIOS eQTL browser 48 , and BRAINEAC 49 were selected for the mapping. Chromatin interaction was based on GSE87112 (Hi-C) database as implemented in FUMA. We explored possible biological functions by pathway analysis for all variants that reached genome wide significance in the discovery stage and were nominal significant in the Asian replication set. These 55 associated variants (ONH = 32, IOP = 3, CCT = 20) were assigned to genes by FUMA mapping strategies. Prioritized genes for ONH traits were highly overlapping and were combined to form a set of 295 unique genes for further functional annotation in FUMA. For IOP and CCT 11 and 116 genes were prioritized respectively. We further investigated the FUMA-mapped genes for enrichment using hypergeometric enrichment tests on pre-defined gene-sets derived from MsigDB and WikiPathways. p Values were corrected based on Bonferroni method for the number of tested gene-sets.
Reporting summary. Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Data availability
The genome-wide summary statistics that support the findings of this study will be made available via the NHGRI-EBI GWAS Catalog website (https://www.ebi.ac.uk/gwas/ downloads/summarystatistics) upon publication.