Bioinformatic Prioritization and Functional Annotation of GWAS-Based Candidate Genes for Primary Open-Angle Glaucoma

Background: Primary open-angle glaucoma (POAG) is the most prevalent glaucoma subtype, but its exact etiology is still unknown. In this study, we aimed to prioritize the most likely ‘causal’ genes and identify functional characteristics and underlying biological pathways of POAG candidate genes. Methods: We used the results of a large POAG genome-wide association analysis study from GERA and UK Biobank cohorts. First, we performed systematic gene-prioritization analyses based on: (i) nearest genes; (ii) nonsynonymous single-nucleotide polymorphisms; (iii) co-regulation analysis; (iv) transcriptome-wide association studies; and (v) epigenomic data. Next, we performed functional enrichment analyses to find overrepresented functional pathways and tissues. Results: We identified 142 prioritized genes, of which 64 were novel for POAG. BICC1, AFAP1, and ABCA1 were the most highly prioritized genes based on four or more lines of evidence. The most significant pathways were related to extracellular matrix turnover, transforming growth factor-β, blood vessel development, and retinoic acid receptor signaling. Ocular tissues such as sclera and trabecular meshwork showed enrichment in prioritized gene expression (>1.5 fold). We found pleiotropy of POAG with intraocular pressure and optic-disc parameters, as well as genetic correlation with hypertension and diabetes-related eye disease. Conclusions: Our findings contribute to a better understanding of the molecular mechanisms underlying glaucoma pathogenesis and have prioritized many novel candidate genes for functional follow-up studies.


Introduction
The term 'glaucoma' refers to a group of ocular disorders characterized by the loss of retinal ganglion cells and the degeneration of their axons [1]. Primary open-angle glaucoma (POAG) is the most common form of glaucoma. While the exact cause of POAG is still unknown, there is clear evidence that age, sex, and intraocular pressure (IOP) are important risk factors. However, genetic factors also play a significant role [1]. Indeed, early evidence from twin and family studies revealed substantial glaucoma heritability [2]. Later on, linkage studies enabled researchers to map the chromosomal locations of a number of rare pathological variants in genes (MYOC, OPTN, and WDR36) that co-segregated with the disease in families [3,4]. More recently, genome-wide association studies (GWASs) have identified a large number of common genomic variants associated with POAG in unrelated individuals. For example, a recent meta-GWAS in 12,315 POAG cases and 227,987 controls, from the Genetic Epidemiology Research on Adult Health and Aging (GERA) cohort and the UK Biobank, identified or replicated more than 70 genomic regions associated with the disease, of which 14 were novel genetic loci [5]. However, a major drawback of GWASs is that the identified variants merely tag genomic regions without providing definitive information on the likely causal genes and functional mechanisms underlying the statistical associations with a particular disease or phenotype [6].
Often, GWASs simply report the nearest genes for each locus. However, in many cases, there are several co-inherited variants in strong linkage disequilibrium (LD) with one another, making it difficult to distinguish the causal variants underlying the statistical association [7]. Moreover, several examples have clearly shown that functionally relevant genes are sometimes located at large distances from the significant GWAS loci [8,9]. This calls for further post-GWAS bioinformatics follow-up studies mapping the identified GWAS associations to the likely causal genes using more robust evidence rather than physical distance to the GWAS loci. Only a limited number of previous studies have attempted to determine the functional characteristics of GWAS-based POAG candidate genes [10,11]. However, a systematic post-GWAS approach that ranks the candidate genes in order of their relevance/causality based on different sources of biological evidence has not been performed yet. Therefore, deep investigations of disease-related genes' expressions across different tissue types are required to better understand disease mechanisms and to avoid off-target reactions in subsequent therapeutic approaches.
The relationship between a GWAS association signal, the so-called single-nucleotide polymorphism (SNP), and the related causal gene (variant) is often not clear. An associated SNP may: (i) alter the amino acid coding (i.e., a nonsynonymous SNP (nsSNP)), changing the protein structure and potentially the function of a gene directly, or (ii) may indirectly exert its phenotypic effect through influencing regulatory sequences and, hence, the expression of a gene. Therefore, bioinformatic post-GWAS pipelines will typically check whether GWAS signals will be in high linkage disequilibrium with nonsynonymous SNPs within nearby genes and use publicly available expression quantitative trait loci (eQTL) resources from relevant tissues to check whether associated SNPs in the identified loci overlap with genomic loci related to altered gene expression [12]. Despite being helpful to some extent, finding such an overlap can also be misleading, as it does not provide conclusive evidence that the target gene expression is also associated with the phenotype. There are a number of ways to provide additional information, including transcriptome-wide association studies (TWAS).
TWAS is a popular approach to test gene expression-phenotype associations. However, this approach does not infer causality, as it is highly vulnerable to environmental confounder effects. Moreover, it requires individual-level phenotype and expression data from the same population, which is hardly ever available on a large scale, limiting the power of detecting true associations. A recently developed new TWAS approach uses previously trained transcriptome models to predict the genetic component of gene expression levels, a method known as transcriptome imputation [13]. This method allows the use of expression and phenotype data from different samples and it concurrently minimizes environmental confounding effects through the use of genetically predicted gene expressions. A recent implementation of this TWAS approach, called MetaXcan, only requires summary statistics to test gene expression associations with the outcome [14,15]. Another summary statisticsbased TWAS approach is summary data-based Mendelian randomization (SMR) [16]. This method uses genetic variants as instrumental variables to test for the causative effects of gene expressions (exposure) on disease outcome. We take advantage of the two latter approaches, MetaXcan and SMR, in the current work. Furthermore, given that genetic variants associated with complex diseases, such as POAG, are predominantly located The resulting 50 independent gSNPs (i.e., after clumping analysis; Table 1) were used as input for performing in silico sequencing and in silico look-ups of pleiotropic associations with other traits in the GWAS catalog. Using a 1Mb region on either side of the independent top gSNPs, a cut-off value for LD was set at r 2 > 0.50 and the analysis was restricted to the European population. Annotation of the gSNPs together with their linked SNPs was carried out using the ANNOVAR software (version 23rd March 2019) [26]. The possible damaging effects of nonsynonymous SNPs on protein structure and function were predicted using the Sorting Intolerant From Tolerant (SIFT) [27] PROVEAN [28] and Polymorphism Phenotyping (PolyPhen) [29] scoring tools. Furthermore, the pleiotropic effect of POAG-associated loci was assessed using the GWAS catalog database (version 17 March 2019) [30].

Co-Regulated Genes within the POAG-Associated Loci
We ran Data-driven Expression Prioritized Integration for Complex Traits (DEPICT) [33] using the full set of POAG GWAS summary statistics ( Figure 1). The significance threshold for index SNPs was set to p < 1 × 10 −5 and clumping was based on r 2 < 0.05 as the LD metric and a physical distance of 500 kb. DEPICT systematically prioritizes the most likely causal genes, gene sets and tissue enrichments based on gene function predictions, even for uncharacterized genes [33]. It performs functional predictions using 14,461 'reconstituted gene sets' based on a curated data set of 77,840 human expression microarrays. DEPICT Figure 1. Summary of the analysis pipeline: phase one shows the gene prioritization pipeline which was performed based on (i) 50 independent gSNPs (via nearest genes and coding consequence) and (ii) the full set of POAG GWAS summary statistics, using co-regulation, gene expression, and epigenetic regulation approaches. It also presents a summary of the quality control thresholds, reference panel used, data used for model predictions, method of data analysis, and the final number of genes prioritized for each part separately and for the entire pipeline. Similarly, phase two summarizes the functional and tissue enrichment analyses, which were performed using the 142 prioritized genes and the full set of POAG GWAS summary statistics.

Co-Regulated Genes within the POAG-Associated Loci
We ran Data-driven Expression Prioritized Integration for Complex Traits (DEPICT) [33] using the full set of POAG GWAS summary statistics ( Figure 1). The significance threshold for index SNPs was set to p < 1 × 10 −5 and clumping was based on r 2 < 0.05 as the LD metric and a physical distance of 500 kb. DEPICT systematically prioritizes the most likely causal genes, gene sets and tissue enrichments based on gene function predictions, even for uncharacterized genes [33]. It performs functional predictions using 14,461 'reconstituted gene sets' based on a curated data set of 77,840 human expression microarrays. DEPICT first identifies all genes at the trait-associated loci and then estimates their cofunctionality to prioritize the most likely causal genes [33]. Summary of the analysis pipeline: phase one shows the gene prioritization pipeline which was performed based on (i) 50 independent gSNPs (via nearest genes and coding consequence) and (ii) the full set of POAG GWAS summary statistics, using co-regulation, gene expression, and epigenetic regulation approaches. It also presents a summary of the quality control thresholds, reference panel used, data used for model predictions, method of data analysis, and the final number of genes prioritized for each part separately and for the entire pipeline. Similarly, phase two summarizes the functional and tissue enrichment analyses, which were performed using the 142 prioritized genes and the full set of POAG GWAS summary statistics.

MetaXcan
In order to identify genes whose expression levels are associated with POAG, free of non-genetic confounders, we performed PrediXcan [13] analysis using the summary data-based pipeline named MetaXcan [15] (Figure 1). This pipeline integrates the summary statistics from genetically predicted transcriptome models with GWAS results to test for the association between the genetic component of gene expression levels and the outcome. We only included GWAS variants with minor allele frequency (MAF) ≥ 1% and, given their high LD ratio, excluding all variants and probes within the major histocompatibility complex (MHC) region. Our analysis was based on the transcriptome model of whole blood from the DGN cohort [34]. The significance level of association p-value was set to 4.39 × 10 −6 regarding 11,397 genes (or variant sets) being tested in the DGN dataset (0.05/11,397 genes). Results were then retained if their prediction performance p-value was smaller than a Bonferroni corrected level (0.05/nsign; with nsign being the number of significant genes).
We then tested whether the associated variants for gene expression (eQTLs) and glaucoma are colocalized using the COLOC R package to address LD contamination concerns in significant genes [35]. The test is based on approximate Bayes factors on five hypotheses of (i) no causal variant (H0), (ii) causal variant for glaucoma only (H1), (iii) causal variant for gene expression only (H2), (iv) two distinct causal variants (H3) and (v) single causal variant for gene expression and glaucoma (H4). We took a probability of H3 < 0.5 and H4 > 0.5 as acceptable evidence of colocalized signals to filter out MetaXcan TWAS association results, which can be due to LD, as suggested by Barbeira et al. [15].

Summary Data-Based Mendelian Randomization (SMR) Based on Gene Expression Data
We performed SMR analysis combined with the Heterogeneity Independent Instruments (HEIDI) test (Figure 1), which jointly uses GWAS summary statistics and eQTL data from independent studies; providing more power to detect causal associations [16,36]. In SMR, the top cis-eQTL for each gene is used as 'instrumental variable' and 'gene expression' is considered the 'exposure' for the phenotype. This Mendelian randomization framework enables a test for the causal effect of the genetic variant (i.e., the eQTL) on the phenotype through gene expression [16]. Since the observation of a significant SMR association may be due to two distinct but linked causal variants, one affecting gene expression and the other influencing the phenotype, we also conducted the HEIDI test to ensure that the trait-gene expression associations are driven by the same genetic variant, not confounded by linkage [16].
We used two gene expression datasets: (i) blood eQTL summary data from the eQTLGen consortium (n~32,000) [37] and (ii) a large set of eQTL data from a meta-analysis of 10 brain regions (n effective~1194) [31]. We used the 1000 Genomes Project phase 3 data panel [24] for LD calculations. The following exclusion criteria were applied after which we retained remaining eligible SNPs for the analysis: (i) rare genetic variants with MAF below 1%, (ii) variant and probes within the MHC region, and (iii) variants with inconsistent alleles or MAF differences > 0.20 amongst pairs of three datasets (GWAS, 1000 Genomes Project, and eQTL datasets). Bonferroni's corrected significance level of p < 3.26 × 10 −6 was set for the blood (0.05/15,352 genes) and p < 6.79 × 10 −6 for the brain (0.05/7361 genes) SMR results. Similarly, a Bonferroni corrected significance level of p ≥ 2.78 × 10 −3 (0.05/18 SMR significant probes) in blood and p ≥ 5.0 × 10 −3 (0.05/10 SMR significant probes) in the brain dataset was set for HEIDI test.
Given the likely role of epigenetics in complex disease, we also conducted epigenomewide studies on POAG with SMR using both blood and brain mQTL data (MSMR, see below). That is, the SMR framework affords a test for the causal effect of the genetic variant (i.e., the mQTL) on the phenotype through methylation. We then mapped the significant methylation sites to their nearest genes using the FDb.InfiniumMethylation.hg19 R package [38] (Figure 1).
To further investigate the underlying mediating mechanism from DNA to POAG, we performed 3xSMR analyses on GWAS of POAG using blood and brain data, separately. The 3xSMR analysis was based on three sets of SMR results: (i) GWAS vs. mQTL (MSMR); in which methylation was the exposure and POAG was the outcome, (ii) mQTL vs. eQTL (MESMR); in which methylation was the exposure and gene expression was the outcome, and (iii) GWAS vs. eQTL (ESMR); in which gene expression was the exposure and POAG was the outcome, as described in the section above. The results of these three sets of analyses are integrated all into one causal model [31,36]. This enabled us to identify the associations between DNA, methylation, and gene expression, which consequently may led to the development of POAG (i.e., DNA→Methylation→Gene expression→POAG) ( Figure 2). Jaffe et al. [41]. We then performed MESMR using the same brain mQTL data against brain eQTL data reported by Qi et al. [31]. In the final step, we completed the 3xSMR analysis using the brain ESMR results from the previously performed analyses described above.
We used the same quality control applied in TWASs, i.e., GWAS variants with MAF below 1%, variants and probes within the MHC region, and variants with inconsistent alleles or MAF differences > 0.20 amongst pairs of four input datasets (GWAS, 1000 Genomes Project, eQTL and mQTL datasets) were excluded. We used the blood mQTL summary data from a meta-analysis (n = 1980) recently reported by Wu et al. [36] and performed MSMR analysis to find the likely causal DNA methylation sites for POAG. Then, to map those methylation sites to the subsequently regulated genes, we repeated SMR analysis on mQTL data [32] using blood eQTL data from the eQTLGen consortium (n = 32,000) [37] (MESMR). Finally, we used the blood ESMR results from the previous steps to retain genes with significant causal associations for all three steps (3xSMR).
For MSMR in the brain, we used the brain mQTL summary data from the Qi et al. meta-analysis (estimated effective n = 1160) [31] of ROSMAP [39], Hannon et al. [40] and Jaffe et al. [41]. We then performed MESMR using the same brain mQTL data against brain eQTL data reported by Qi et al. [31]. In the final step, we completed the 3xSMR analysis using the brain ESMR results from the previously performed analyses described above.

Functional and Tissue Enrichment Analysis
The biological pathways through which POAG-associated loci act were examined using two approaches. First, we used DEPICT to identify the enriched gene sets, alongside the corresponding functions, and tissue types for genes in POAG-associated loci. DEPICT uses the same data as for functional predictions and gene prioritization (see phase one), for gene set enrichment analysis, and a set of 37,427 human gene expression microarrays for tissue enrichment analysis of 209 tissue/cell type annotations. For gene set enrichment analysis, we applied the Affinity Propagation Clustering algorithm (APCluster R package [42]), as previously suggested and implemented by Ligthart et al. [43]. The clustering was conducted based on pairwise correlation of gene sets. Second, the prioritized gene lists based on phase one (i.e., nearest genes, genes with nsSNPs linked to POAG loci, and significant genes from DEPICT, ESMR, MSMR, and MetaXcan analyses) were merged and used as the input to run functional enrichment analysis with the GeneMANIA algorithm [44] as previously described [23] and with the Ingenuity Pathway Analysis (IPA; QIAGEN) software ( Figure 1). IPA core-analysis yields the relationships, canonical pathways, diseases, and functions most relevant to the uploaded set of prioritized genes. In addition, we also performed a sensitivity functional enrichment analysis in GeneMANIA, using a subset of genes that showed two or more lines of evidence out of all five approaches that we used (see methods). Furthermore, we used the OTDB described by Wagner et al. [45] to assess whether prioritized genes are overrepresented in ocular tissues. The database contains microarray gene expression values of >20,000 genes in ten human ocular tissues (choroid, ciliary body, cornea, iris, lens, optic nerve, optic nerve head, retina, sclera, and trabecular meshwork [TM]). We performed right-tailed Fisher's exact tests [46] to evaluate the enrichment of prioritized genes amongst the genes with gene expression values in the top 25% [47].

Genetic Correlation of POAG with Other Traits
We applied the (bivariate) LD score regression method [48] to estimate the genetic correlation of POAG with 597 UK Biobank traits with available GWAS summary statistics. POAG GWAS summary statistics data (~7.7 million SNPs) were uploaded to LDHub, an online tool dedicated to estimating the genetic correlation of traits of interest [49].

In Silico Sequencing
In silico sequencing of 50 independent lead SNPs (gSNPs) from Choquet et al. [5] returned 3250 and 1493 SNPs that are, respectively, in moderate (r 2 > 0.50) and high (r 2 > 0.80) LD (Supplementary Table S1). One hundred and ten of the aforementioned 1493 SNPs were in complete LD (r 2 = 1), indicating that these linked SNPs represent the same association signal. Annotation of these 3250 linked SNPs using the ANNOVAR software [26] detected nine nsSNPs. Two of the nine nsSNPs (rs3753841 and rs2274224) were in perfect LD (r 2 = 1) with the corresponding linked gSNPs (rs993471 and rs3891783, respectively; Supplementary Table S1). These nine nsSNPs were located in seven genes, of which one (ACP2) was novel for POAG and four (ACP2, SH2B3, SIX6, and C14orf39) did not overlap with the nearest gene list of the 50 gSNPs (Table 2).

MetaXcan
After the quality control steps, 11,397 variant sets were used to impute gene expression levels using the Depression Gene Network (DGN) transcriptome model. From the total of 14 significant genes, whose predicted gene expressions were associated with POAG, two genes with prediction performance p ≥ 3.57 × 10 −3 (0.05/14) were filtered out and 12 significant genes remained, of which eight showed acceptable evidence of co-localization signals (Table 3). Three colocalized genes identified in MetaXcan analysis (NR1H3, LTBP3, and EHBP1L1) were also significant in SMR analysis (see below) of blood eQTL ( Table 2).

SMR Based on Gene-Expression Data
After applying the quality control criteria, 15,352 genes were retained for the SMR analysis. Using blood gene expression profile as a source, our analysis yielded 18 genes previously implicated/associated with glaucoma. Thirteen out of these 18 have no significant evidence of linkage confounding (based on heterogeneity in dependent instruments [HEIDI] p ≥ 2.78 × 10 −3 ; Table 4). Three (BICC1, LTBP3, and ABCA1) out of these 13 significant genes were also identified by at least two additional prioritization methods. Nine out of the 13 identified genes were novel for POAG (Table 2).          Similarly, in the brain eQTL, 7361 genes were eligible for SMR analysis. Our analysis returned 10 genes whose expression profile was significantly associated with glaucoma; eight of which have no significant evidence of linkage confounding (based on heterogeneity test, HEIDI p ≥ 5 × 10 −3 ) ( Table 4). Two of these eight genes, TXNRD2 and CDKN2B-AS1, overlapped with nearest genes and one (CDKN2B) with genes prioritized in the DEPICT analysis. Interestingly, the significant association of three HEIDI-passed genes (LRRC37A2, LRRC37A4P, and RP11-707O23.5) was also confirmed by SMR analysis of blood eQTL.

SMR Based on Methylation Data and 3xSMR
Methylation QTL (mQTL)-based SMR analysis (MSMR) of blood and brain data yielded, respectively, 14 (SMR p < 5.60 × 10 −7 ) and 16 (SMR p < 5.39 × 10 −7 ) methylation sites significantly associated with glaucoma, with no evidence of LD contamination (based on heterogeneity test; HEIDI p ≥ 1.28 × 10 −3 and HEIDI p ≥ 1.39 × 10 −3 , respectively; Table 5). Collectively, our SMR analysis of POAG GWAS and mQTL of blood and brain identified 27 novel CpG sites in 15 genes. Of these, two genes are new (AFAP1-AS1 and TBKBP1), i.e., were not identified by the previous prioritization methods. Amongst significant MSMR genes, AFAP1-AS1 and NR1H3 are novel for POAG ( Table 2). The integrative analysis of mQTL and eQTL data in the blood (MESMR) resulted in 32,420 DNA methylation (DNAm) sites significantly associated (SMRp < 1.85 × 10 −8 ; regarding Bonferroni correction based on 2,697,257 tests) with expression levels of 10,680 genes not rejected by the HEIDI test (p ≥ 3.26 × 10 −7 ). These results were used to link glaucomaassociated gene expressions to glaucoma-associated methylation levels, and identified the BICC1 gene with its genetic regulation in glaucoma to be explained by a likely causal chain (DNA→Methylation→Expression→Glaucoma). The integrative analysis of mQTL and eQTL data in the brain (MESMR) resulted in 10,685 DNAm sites significantly associated (SMRp < 1.44 × 10 −8 ; regarding Bonferroni correction based on 3,482,629 tests) with expression levels of 3305 genes not rejected by the HEIDI test (p ≥ 2.16 × 10 −6 ). We did not detect a likely causal chain from DNA to glaucoma through methylation and expression in the brain.

Integration of Results (Phase One)
Taken all these results together, AFAP1 and BICC1 were simultaneously highlighted in five of the six approaches and 34 of the 142 prioritized genes were supported by at least two lines of evidence (Table 2), suggesting that these genes are involved in glaucoma development by multiple mechanisms including coding, transcription, and manipulation of targeted gene expression. On the other hand, less than half of the nearest genes (n = 23) overlapped with prioritized genes identified in the pipeline ( Table 2). That is, 92 likely causal genes identified by different prioritization approaches were not the genes nearest to the lead gSNPs.

Functional and Tissue Enrichment Analysis
We describe in this section the results of different functional and tissue enrichment analyses. The former include results from GeneMANIA Gene Ontology (GO) enrichment analysis, ingenuity pathways analysis (IPA) functional enrichment analysis of its manually curated biological pathways, DEPICT gene set enrichment (GSE) analysis, and affinity propagation clustering (APC) of its enriched gene sets. The latter, i.e., tissue prioritization, is based on DEPICT tissue enrichment analysis with its expression data from 209 cell types/tissues, and evaluating gene expression of the prioritized genes in 10 ocular tissues. The results of these analyses are described in order below.  GeneMANIA functional enrichment analysis revealed 246 significantly enriched gene ontology (GO) terms for 142 prioritized genes (Supplementary Table S3A). Sensitivity analysis in the 34 genes with two or more lines of evidence yielded 174 significantly enriched GO terms (Supplementary Table S3B). Along with several closely related pathways, most of the most significant GO terms were related to the extracellular matrix (ECM) (GO:0031012, q = 6.95 × 10 −56 ), transforming growth factor-β (GO:0007179, q = 4.20 × 10 −11 ), cardiovascular traits (e.g., blood vessel development, GO:0001568, q-value=2.40 × 10 −7 ), heart development (GO:0007507, q = 2.52 × 10 −6 ), Wnt signaling (GO:0016055, q = 2.52 × 10 −6 ), retinoic acid receptor binding (GO:0042974, q = 1.18 × 10 −4 ), and eye development (GO:001654, q = 3.09 × 10 −3 ; Supplementary Table S3A,B). Several of these pathways, including heart development, regulation of protein phosphorylation, embryo development, and Wnt-protein signaling, were in line with the most prioritized terms from DEPICT. Similarly, IPA identified 64 canonical signaling pathways (Supplementary Table S3C). Based on the high percentage of focus molecules in our datasets, RAR activation (p = 1.32 × 10 −5 ), leptin signaling (p = 7.94 × 10 −4 ), and aryl hydrocarbon receptor signaling (p = 1.09 × 10 −3 ) were the most strongly enriched canonical signaling pathways constructed in IPA.
Using the full set of POAG GWAS-summary statistics [5], DEPICT's functional enrichment (FE) analysis identified 269 biological pathways enriched by glaucoma-associated loci (FDR < 0.05; Supplementary Table S2B). These pathways were based on five annotation categories (see Figure 3), which include Gene Ontology (GO), Kyoto Encyclopedia of Gene and Genomes (KE), REACTOME (RE), Mouse Phenotypes (MP), and Protein-Protein Interactions (PI). Cardiovascular-related terms, e.g., abnormal aorta morphology from the MP pathway resource, and Wnt and TGF-β signaling terms from the KE pathway resource, were significantly enriched at FDR < 0.05 (Figure 3), suggesting that dysregulation in any of these pathways may contribute to POAG development. For the sake of briefness, we only visualized the top five gene sets per each annotation category ( Figure 3). Furthermore, affinity propagation clustering of enriched gene sets yielded 37 gene set clusters at FDR < 0.05, including artery morphology, vasculature development and cell motility regulation. The list of all cluster centers (nodes) and pairwise Pearson correlation (edges) between the nodes is summarized in Figure 4. In the DEPICT TE (tissue/cell type) enrichment analysis, 20 tissues were prioritized at FDR < 0.05 including tissues from the urogenital system and more specifically female reproductive organs, aortic and heart valves, and arteries (Supplementary Table S2C). Based on the assumed relevance for POAG, here we only show results for sense organ tissues (SO), exclusively those from the eye, as well as from the top ten nervous system tissues ( Figure 3).    Our ocular tissue database (OTDB) investigations using Fisher's exact test showed a statistically significant overrepresentation of prioritized genes in four (sclera, TM, ciliary body, and choroid) of the 10 ocular tissues (p-value < 0.005; i.e., 0.05/10) ( Figure 5). More specifically, prioritized genes were overrepresented in the sclera, TM, and ciliary body by >1.5 fold.

In Silico Pleiotropy Look-Up and Genetic Correlation
In silico pleiotropy analysis using GWAS catalog with European subpopulation yielded 139 linked SNPs (LD, r 2 > 0.80), that were previously associated with complex disease types including 41 with glaucoma itself. A complete list of all the traits identified by the in silico pleiotropic look-up along with the linked SNPs (r 2 > 0.50) and their nearest genes is summarized in Supplementary Table S1. In Figure 6, we more specifically present the number of shared highly linked SNPs (r 2 > 0.80) between glaucoma and six phenotypic categories that have been hypothesized in the literature to have a relationship with glaucoma, i.e., traits and diseases related to the eye (e.g., intraocular pressure) [50], anthropometry (e.g., body mass index) [22], blood pressure [19], lipids and type 2 diabetes [22], neurodegenerative disorder (e.g., Alzheimer's disease) [21], and psychiatry (e.g., depression) [51]. Compared to other traits and diseases, POAG showed the highest genetic overlap with IOP and vertical cup-to-disc ratio, both of which are glaucoma endophenotypes with considerable heritability [2]. This can be considered as a proof of concept and internal validation of the pipeline.

In Silico Pleiotropy Look-Up and Genetic Correlation
In silico pleiotropy analysis using GWAS catalog with European subpopu yielded 139 linked SNPs (LD, r 2 > 0.80), that were previously associated with com disease types including 41 with glaucoma itself. A complete list of all the traits iden by the in silico pleiotropic look-up along with the linked SNPs (r 2 > 0.50) and their ne genes is summarized in Supplementary Table S1. In Figure 6, we more specifically pr the number of shared highly linked SNPs (r 2 > 0.80) between glaucoma and six pheno categories that have been hypothesized in the literature to have a relationship with coma, i.e., traits and diseases related to the eye (e.g., intraocular pressure) [50], an pometry (e.g., body mass index) [22], blood pressure [19], lipids and type 2 diabetes neurodegenerative disorder (e.g., Alzheimer's disease) [21], and psychiatry (e.g., de sion) [51]. Compared to other traits and diseases, POAG showed the highest genetic lap with IOP and vertical cup-to-disc ratio, both of which are glaucoma endopheno with considerable heritability [2]. This can be considered as a proof of concept and in We also visualized genetic correlations (rg) of POAG with traits from six phenotypic categories (Figure 7). POAG showed significant correlations with cardiometabolic disease (hypertension, rg = 0.08, p = 0.011), blood pressure (systolic blood pressure, rg = 0.06, p = 0.039; high blood pressure, rg = 0.08, p = 0.012) and eye diseases (diabetes-related eye disease, rg = 0.20, p = 0.025, other eye problems, rg = 0.44, p = 4.44 × 10 −9 , and senile cataract, rg = 0.27, p = 0.022; Figure 7). No significant genetic correlations were observed between glaucoma and Alzheimer's disease or depression. Except for other eye problems, correlations were no longer significant after correcting for multiple testing of 597 UK Biobank traits.

Integration of Results (Phase Two)
Some signaling pathways, e.g., developmental (Wnt/β-catenin signaling, p = 2.82 × 10 −3 ), retinoic acid receptor activation (p = 1.32 × 10 −5 ), and cardiac-related (Cardiac Hypertrophy Signaling, p = 1.10 × 10 −2 ), were in common with biological processes enriched in DEPICT and GeneMANIA. In line with functional enrichment results, tissue investigations also highlighted the eye, as well as cardiovascular tissues, to be the appropriate contexts of POAG genes. In completion, pleiotropy and genetic correlation analyses also revealed shared genetic loci between glaucoma and IOP, cup-to-disc ratio, as well as blood pressure.

Figure 7.
Dot and whisker plot: the genetic correlations and 95% confidence intervals between POAG and six phenotypic categories were calculated using the LD score regression method.

Integration of Results (Phase Two)
Some signaling pathways, e.g., developmental (Wnt/β-catenin signaling, p = 2.82 × 10 −3 ), retinoic acid receptor activation (p = 1.32 × 10 −5 ), and cardiac-related (Cardiac Hypertrophy Signaling, p = 1.10 × 10 −2 ), were in common with biological processes enriched in DEPICT and GeneMANIA. In line with functional enrichment results, tissue investigations also highlighted the eye, as well as cardiovascular tissues, to be the appropriate contexts of POAG genes. In completion, pleiotropy and genetic correlation analyses also revealed shared genetic loci between glaucoma and IOP, cup-to-disc ratio, as well as blood pressure.

Discussion
We aimed to prioritize the most likely causal genes and identify the underlying biological processes involved in POAG through post-GWAS bioinformatics analyses. Our systematic post-GWAS approach spotted 142 genes as the most likely causal and/or relevant genes for POAG, 64 of which were novel. Among the prioritized genes were seven genes with nsSNPs linked to POAG genomic loci (LD r 2 > 0.5), and 34 (23.9%) of the prioritized genes were supported by at least two lines of evidence. With considerable overlap, DEPICT and GeneMANIA functional enrichment analysis revealed 269 and 246 signaling pathways associated with glaucoma, respectively. ECM was the most significant and repeatedly implicated pathway in GeneMANIA. Along with several closely related pathways, TGF-β signaling, blood vessel development, heart development, and retinoic acid receptor signaling were also significantly overrepresented. Furthermore, tissues from the female reproductive as well as the cardiovascular system were significantly enriched

Discussion
We aimed to prioritize the most likely causal genes and identify the underlying biological processes involved in POAG through post-GWAS bioinformatics analyses. Our systematic post-GWAS approach spotted 142 genes as the most likely causal and/or relevant genes for POAG, 64 of which were novel. Among the prioritized genes were seven genes with nsSNPs linked to POAG genomic loci (LD r 2 > 0.5), and 34 (23.9%) of the prioritized genes were supported by at least two lines of evidence. With considerable overlap, DEPICT and GeneMANIA functional enrichment analysis revealed 269 and 246 signaling pathways associated with glaucoma, respectively. ECM was the most significant and repeatedly implicated pathway in GeneMANIA. Along with several closely related pathways, TGF-β signaling, blood vessel development, heart development, and retinoic acid receptor signaling were also significantly overrepresented. Furthermore, tissues from the female reproductive as well as the cardiovascular system were significantly enriched for POAG genes. POAG-prioritized genes were overrepresented in four ocular tissues: sclera, ciliary body, TM, and choroid.

Highlights of Separate Analyses
Six of the 64 novel POAG genes (NR1H3, ACP2, EHBP1L1, LRRC37A2, LRRC37A4P, and RP11-707O23.5) presented two or more lines of evidence (Table 2). NR1H3 is one of the top genes associated with IOP [52], a prominent risk factor for POAG, whereas ACP2 is found to be overexpressed in the cerebellum and brain stem in neuronal ceroid lipofuscinoses (CLN3) mice [53]. CLN3 disease is an inherited disorder that affects the nervous system, and children with this disorder are characterized by progressive neurological degeneration and vision loss. Using data derived from GWAS of large consortia, previous reports showed the association of EHBP1L1 with myopia and BP (diastolic BP and mean arterial pressure) [54,55]. Both myopia [56] and BP [19] are risk factors of glaucoma. Further studies are required to uncover how these novel genes are relevant to POAG.
Of the nine linked nsSNPs identified, two (rs3753841 and rs2274224) candidate variants were perfectly linked (LD r 2 = 1) to POAG lead gSNPs, and rs3753841 mapping to COL11A1 gene was predicted to be 'deleterious' and 'possibly damaging' by SIFT [27] and Polyphen [29], respectively. Indeed, SNP rs3753841 has been reported previously to be associated with glaucoma [57,58]. This result provides evidence that rs3753841 alters an amino acid sequence in the collagen α chain precursor protein, a major component of the ECM, and contributes to susceptibility of glaucoma. This may be explained through the contribution of the modified collagen to IOP induction by generating outflow resistance in aqueous humor [59].
Furthermore, nsSNP rs8940 has previously been annotated in the coding region of CAV2 in Australian and Swedish glaucoma patients, but its association with POAG was not significant after adjusting for the effects of a correlated (LD r 2 = 0.63) SNP rs4236601 [60]. CAV2 codes for caveolin protein family members, which is a specialized plasma membrane raft forming flask-shaped invaginations. It is involved in cell proliferation, transcellular transport, membrane lipid homeostasis, mechanotransduction, and signal transduction [61]. Caveolin 1 and caveolin 2 inhibit endothelial nitric oxide synthase enzyme activity within the caveolae, and alteration of this pathway has been associated with abnormal nitric oxide generation and TM function [62]. Caveolae are available in different retinal cell types including retinal vascular, retinal pigment epithelium, and Müller glia cells [61]. Cav-1, the principal protein of caveolae, modulates neuroprotective responses, and ablation of Cav-1 in mice and zebrafish was associated with defects in retinal pigment epithelium differentiation and STAT3 activation in the retina [63,64].
Of the 50 nearest genes identified in GWAS, only 23 (46%) genes were detected in the gene prioritization analyses, suggesting that not all nearest genes have a functional impact on protein coding, gene expression, or the regulation of transcription. This highlights the limited mapping resolution of GWAS results due to the complicated linkage disequilibrium structure of the genome, i.e., GWAS lead SNPs are not always the causal functional variants. This has been confirmed by a previous study which reported about 80% of the common GWAS variants are within 33.5 Kbp of the underlying causal variants [65]. Furthermore, our SMR analysis using methylation data also detected significant loci, strengthening the hypothesis that while SNPs located in the coding region may act by altering amino acid sequence and protein functions, genetic variants in the non-coding region may cause diseases through other mechanisms, including methylation and regulation of gene expression [66].
Whilst the role of epigenetics is well studied in mental disorders [67] and cancer [68], exploration of the role of DNA methylation in eye diseases is limited [69]. One study reported the association of a CpG site at CDKN2B in normal-tension glaucoma [70], another study found CpG sites at SKI and GTF2H4 in the retinal pigment epithelium of age-related macular degeneration patients [71]. Our combined SMR analysis of POAG GWAS data and mQTL of blood and brain data identified 27 novel DNA methylation sites for POAG in 15 genes, highlighting the role of epigenetics in gene expression and glaucoma. Moreover, our 3xSMR analysis, which was used for linking glaucoma-associated gene expression to glaucoma-associated methylation levels, enabled us to find two methylation sites (cg05938607 and cg12342675) at the BICC1 gene, confirming its likely causal role in methylation mediated genetic regulation in glaucoma (DNA→methylation→expression→glaucoma; Figure 2).

Post-GWAS Analyses Yielding Most Reliable Data
Three genetic loci (AFAP1, BICC1, and ABCA1) were identified in four or more approaches in our pipeline, confirming earlier evidence that these genes have a likely role in glaucoma pathogenesis. Actin filament associated protein 1 (AFAP1), also known as AFAP110 or AFAP-110, encodes the nonreceptor tyrosine kinase (Src) binding partner protein and affects actin filament organization in response to cellular signals [72]. Activation of Src signaling, in turn, has been implicated in the attenuation of ECM degradation via the inhibition of plasminogen activator expression, and the application of dasatinib, a potent Src signaling inhibitor, in rats improved the TGF-β2-induced adhesive and contractile characteristics of the TM and also attenuated ECM deposition [73]. BICC1, Bicaudal-C (BicC) family RNA binding protein 1, is involved in gene expression during embryonic development and is a negative regulator of the Wnt signaling pathway [74]. Previous studies reported the association of BICC1 with POAG [5,75], major depressive disorder [76], and high myopia [77]. The GTEx, CAGE, and FANTOM5 projects showed high average RNA expression levels in the eye [78]. BICC1 and AFAP1 were the two most highly prioritized genes; we further searched DEPICT gene set enrichment results in order to find the most likely relevant pathways through which these two genes act. Our co-regulation analyses predicted the highest functional similarity of BICC1 with an ECM-related cluster of gene sets. Within the cluster, the strongest correlation of BICC1 points towards abnormal tendon morphology, with corneal thinning also among the gene sets in the cluster ( Figure 8A). In addition to cell motility, our analysis predicts AFAP1's potential involvement in vasculature development as well as neuron differentiation ( Figure 8B).
ATP-binding cassette transporter A1 (ABCA1) protein is a cholesterol and phospholipid transporter across the cell membrane, and mutation of this gene is associated with cancer and lipoprotein metabolism abnormality [79,80]. Although the association of AFAP1 and ABCA1 genes with cancer has been widely studied [81,82], there is a paucity of data showing the underlying mechanism in glaucoma. Cui et al. used reverse transcription polymerase chain reaction and immunolabeling approaches to examine the expression of AFAP1, ABCA1, and GMDS in human ocular tissues. Both AFAP1 and ABCA1 showed significant gene expression in the retina, retinal ganglion, optic nerve, and TM cells [83]. Our TWAS analyses showed that in contrast to the potential protective effect of increased AFAP1 expression, ABCA1 expression level is positively correlated with the risk of POAG. In line with this result, there is supporting evidence that ABCA1 is involved in the apoptosis of retinal ganglion cells in rat glaucoma models [84].

Functional Assessments
POAG is a multifactorial disease with multiple pathways involving different tissues. In normal human development, complex signaling pathways control cells' proliferation, differentiation and fate, motility, and death; thus hypo/hyperactivation of either of these pathways may lead to the development of chronic diseases including POAG [85]. Functional enrichment analysis showed that the 142 prioritized causal genes act through several pathways, but most importantly, candidate genes were overrepresented in ECM, TGF-β signaling, blood vessel development, heart development, angiogenesis, and the retinoic acid receptor signaling pathway. Overrepresentation of the ECM and TGF-β pathways is in line with previous in vitro studies which reported TGF-β-induced morphological changes in the ocular structures, and possibly the development of POAG, using human TM cells as well as optic nerve astrocytes [86].
The ECM is a group of more than 300 complex and dynamic proteins regulating many biological processes and structures. The composition of ECM is unique in each organ, and it is continuously remodeled to regulate tissue homeostasis [87]. Previous studies have suggested that the cross-linking and deposition of ECM proteins in TM are responsible for aqueous humor outflow resistance and IOP elevation in glaucoma patients [86]. Additionally, glaucoma patients have elevated levels of TGF-β in their aqueous humor, and high TGF-β has been shown to increase the synthesis of ECM in human TM cells [86,88]. Few studies have described the mechanisms and cascade of events in how ECM leads to TM dysfunction and increased IOP. Prior experimental studies reported that the induction of ECM cross-linking LOX genes by three TGF-β isoform proteins change ECM stiffness in human glaucomatous TM cells. The authors inferred that TGF-β-mediated overexpression of LOX activity is partially responsible for the increased aqueous humor outflow resistance and IOP elevation [89].  TGF-β is a polypeptide whose signaling has a pleiotropic effect in several cell functions depending on the cellular context. For example, TGF-β receptor signaling, mediated via canonical SMAD pathways, controls the expression of hundreds of genes, while via noncanonical pathways, it regulates cell polarity and microRNA maturation [90,91]. Although TGF-β in normal tissues has pleiotropic roles in multiple biological processes and does not lead to fibrosis, it has a major impact on fibrosis and scarring following glaucoma surgery [92]. Of the three TGF-β isoforms (TGF-β1, TGF-β2, and TGF-β3), increased levels of TGF-β2 in the aqueous humor was associated with fibrosis of TM in POAG patients [93]. In addition, a recent study reported that the TGF-β-induced increase in collagen expression by TM cells is linked with phosphorylation of PTEN, and the manipulation of PTEN activity has been suggested as having therapeutic potential to prevent fibrosis of TM in POAG patients [94]. Other studies demonstrated that TGF-β2 increases expression of PAI-1 and prevents the activation of matrix metalloproteinases (enzymes responsible for the degradation of most ECM proteins) via both tissue-type (tPA) and urokinasetype urokinase (uPA) plasmin activators, thereby increasing the cross-linking of ECM components in TM cells via transglutaminase [95,96]. Another in vitro study demonstrated a cascade of events initiated by TGF-β2 could lead to ECM deposition. TGF-β2 injected into human TM cultured cells significantly increased fibronectin (ECM glycoprotein) levels. In non-glaucomatous eyes, TM cells secrete BMP-4 and antagonize the TGF-β2 pathway. In contrast, in POAG eyes, Gremlin (glycoprotein) is overexpressed in the TM cells, and this inhibits the antagonistic effect of BMP-4 on TGF-β2, thereby leading to increased ECM crosslinking and elevated IOP [97].
Overrepresentation of prioritized genes in cardiovascular-related pathways is consistent with the fact that abnormal blood vessel growth, through the activation of vascular endothelial growth factor (VEGF), has been implicated in neovascular glaucoma, agerelated macular degeneration, and diabetic retinopathy [98]. Similarly, insufficient blood supply to the optic nerve head, either due to elevated IOP or low BP, is a prominent risk factor put forward to explain the development and progression of glaucoma through epidemiological studies [99]. Wnt signaling, combined with several signaling pathways, controls crucial tissues and organs during embryonic development [100]. Accordingly, a recent review showed that Wnt signaling has a vital role in angiogenesis and vascular morphogenesis in eye development, and increased activation of this signal was implicated in diabetic retinopathy, wet age-related macular degeneration, corneal neovascularization, and retinopathy of prematurity [101].
Our functional analysis implied the significance of retinoic acid (RA) and its signaling pathways in glaucoma pathogenesis. RA, a derivative of vitamin A, is essential for regulating genetic transcription that controls a wide range of biological processes including development and homeostasis, as well as cellular apoptosis and survival [102]. The identification of this pathway may strengthen the previous evidence that reported the involvement of CYP26A1 (which was included among the list of prioritized genes) and CYP26C1 in the regulation of RA metabolism, eye development, and maturation of visual function [103]. Moreover, RAs and their receptors (RARs and RXRs) have an important role in eye development and physiology [104], and in a recent study, Prat et al. demonstrated that myocilin (a pathogenic genetic biomarker for glaucoma) expression is regulated by RA via the RARE-DR2 promotor [105]. Furthermore, intraocular injection of RA receptor agonists significantly reduced the number of injured axons of the retinal ganglion cells in frogs, while the injection of antagonists and also the inhibition of RA synthesis with disulfiram had the opposite effect [106]. The role of RA signaling in the pathogenesis of POAG via myocilin transcription, or through its direct effect on retinal ganglion cells, may provide clues for new therapeutic approaches in glaucoma. Our findings highlight the urogenital system and, more specifically, female reproductive organs as tissues in which POAG genes are highly expressed. This is in line with the previous evidence of a protective role for female hormones in glaucoma [107][108][109]. Aortic and heart valves, which were among the most enriched tissues, support the results of our functional enrichment analysis regarding vasculature development and artery morphology, as well as previous studies [10]. These observations may also partly explain the shared genetics that we observed for glaucoma and blood pressure. Additional research is required to elucidate the relevance of digestive tract systems and to determine if this is only explained by the potential role of the gut microbiota in glaucoma progression, through the gut-retina axis [110]. Furthermore, the expression of POAG genes in sclera and choroid is in agreement with previous studies that showed an association between ocular rigidity (structural stiffness of the ocular tissues, including sclera, choroid, and cornea) and glaucoma [111,112]. For example, the biomechanical response of the optic nerve head to IOP stress is determined by the mechanical properties of the sclera, i.e., the sclera plays an important role in transmitting forces created by changes in IOP to the optic nerve head, thereby linking the pressure-glaucoma damage pathway [113,114].
In silico pleiotropy look-ups in the GWAS catalog confirmed the already-known genetic overlap between POAG and its related endophenotypes. For example, sixty-five linked glaucoma SNPs showed a pleiotropic effect on IOP, which is a modifiable risk factor for glaucoma ( Figure 6). In addition, several glaucoma-linked SNPs were associated with optic disc parameters, including optic cup area, optic disc area, and vertical cup-to-disc ratio, a clinical metric of glaucoma [115].
Previous epidemiological studies have also suggested a close relationship between glaucoma and the two most common neurodegenerative diseases (Alzheimer's and Parkinson's disease), suggesting that glaucoma is the consequence of an age-related neurodegenerative process that affects the visual system [21,116]. Additionally, clinical studies have demonstrated that extracellular amyloid β (Aβ) accumulation, hyperphosphorylated tau protein (pTau) aggregation, and glial reaction are common pathological mechanisms in both glaucoma and Alzheimer's diseases [117,118]. However, our pleiotropic look-ups and bivariate LD score regression analysis demonstrated no significant genetic correlation, highlighting that non-genetic factors (environmental factors), or undiscovered confounders, might account for the reported relationships. Our functional enrichment analysis outcomes agreed with previous comparable studies that reported the overrepresentation of ECM, TGF-β, lipid metabolism, developmental, and cardiovascular-related pathways in glaucoma, but could not confirm the importance of inflammation-related factors (e.g., nuclear factor-κB (NFκB), and tumor necrosis factor α (TNF-α)) [10,11]. The most plausible explanation for the apparent discrepancy, regarding the inflammation-related pathway, is the difference in the gene prioritization methodology used. We included genes that may have a role through different mechanisms (coding consequence, gene expression, co-regulation, and epigenetic regulation), whereas the data used as the input for the functional enrichment analysis in Janssen et al. [11] and Danford et al. [10] were significant genes curated from previous linkage, candidate-gene, and genome-wide association studies. IPA [10,11] and ConsensuspathDB [10] were used for functional annotation, whereas in addition to IPA, we used DEPICT as well as the GeneMANIA algorithm, along with its corresponding composite networks based on about 600 million interactions, which may also partly explain the difference in the identified functional pathways.

Strengths and Limitations
Despite the very large sample sizes of our input data sets and different statistical approaches used to find likely causal genes and pathways, there are still some limitations regarding our work. More specifically, we had to use expression data from the blood and brain to increase power for TWAS analyses, because only very limited sample sizes are available for eye tissues. Inevitably, as in any other research, our work is built on the basis of the previous studies in terms of data or methods, which may have biased our findings towards previous knowledge. However, we used agnostic methods and included data from large-scale studies, offering ample opportunities to uncover new aspects of the disease, while concurrently avoiding false interpretations by using conservative significance cut-off values by adjusting for multiple hypothesis testing. We acknowledge that while we were conducting this research, 19 (NR1H3, ACP2, EHBP1L1, ANGPTL2, CALD1, CLIC5, COL16A1,  H1F0, KLF5, KREMEN1, LOC100147773, MAPT, MIR4778, MNAT1, NPEPPS, SLC2A12,  THSD4, TRIB2, TRIOBP) of the 64 novel prioritized genes were detected in two recent POAG GWAS publications [119,120].

Conclusions
We shed light on the underlying biology of glaucoma in terms of likely involved genes and pathways using the most comprehensive recent GWAS on POAG. Our prioritized genes, particularly those with multiple lines of evidence, are eligible for further in vitro and in vivo studies. Such studies, together with improved knowledge on highlighted biological pathways including ECM organization in TM, retinoic acid signaling as well as blood vessel development, may ultimately yield new and more effective therapies. Author Contributions: N.G.A. designed the study, performed pleiotropy and genetic correlation analyses, gene expression assessment and enrichment analysis of the ocular tissues, functional analyses using IPA, the interpretation of the results, and drafted the paper. Z.K. conducted in silico sequencing, DEPICT, TWASs in blood and brain, MSMR analyses in blood and brain, functional analysis using GeneMANIA, and was a major contributor in writing the manuscript. S.P. contributed to the initial review of published glaucoma papers and collecting GWAS loci. A.V., N.J., A.A.B., and H.S. performed critical reviews and supervision over the analyses and the draft. All authors have read and agreed to the published version of the manuscript.
Funding: This project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 661883. Additional funding has been provided by the Rotterdamse Stichting Blindenbelangen under grant number B20150036. The funders had no role in study design, data collection, and analysis, decision to publish, or preparation of the manuscript.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The POAG GWAS dataset analyzed during the current study is available in the GWAS catalog repository (http://ftp.ebi.ac.uk/pub/databases/gwas/summary_ statistics/ under Study accession GCST006065, accessed on 21 April 2019). Whole blood and brain eQTL datasets for TWAS analyses are available from eQTLGen website (https://www.eqtlgen.org/, accessed on 10 December 2019) and SMR data resource (https://cnsgenomics.com/software/smr/, accessed on 19 December 2019), respectively. Gene expression data for ocular tissues are available at Ocular Tissue DataBase (OTDB) (https://genome.uiowa.edu/otdb/, accessed on 23 January 2020). All results generated during this study are included in this published article and its Supplementary Information Files.