Whole-exome sequencing study identifies rare variants and genes associated with intraocular pressure and glaucoma

Abstract

Elevated intraocular pressure (IOP) is a major risk factor for glaucoma, the leading cause of irreversible blindness worldwide. IOP is also the only modifiable risk factor for glaucoma. Previous genome-wide association studies have established the contribution of common genetic variants to IOP. The role of rare variants for IOP was unknown. Using whole exome sequencing data from 110,260 participants in the UK Biobank (UKB), we conducted the largest exome-wide association study of IOP to date. In addition to confirming known IOP genes, we identified 40 novel rare-variant genes for IOP, such as BOD1L1, ACAD10 and HLA-B, demonstrating the power of including and aggregating rare variants in gene discovery. About half of these IOP genes are also associated with glaucoma phenotypes in UKB and the FinnGen cohort. Six of these genes, i.e. ADRB1, PTPRB, RPL26, RPL10A, EGLN2, and MTOR, are drug targets that are either established for clinical treatment or in clinical trials. Furthermore, we constructed a rare-variant polygenic risk score and showed its significant association with glaucoma in independent participants (n = 312,825). We demonstrated the value of rare variants to enhance our understanding of the biological mechanisms regulating IOP and uncovered potential therapeutic targets for glaucoma.

Introduction

Elevated intraocular pressure (IOP) is a major risk factor for glaucoma, the leading cause of irreversible blindness worldwide. IOP is also the only modifiable risk factor for glaucoma. Current glaucoma drugs target lowering IOP. Previous studies using directly genotyped and imputed genetic data have uncovered common and some low-frequency variants for IOP^1,2,3,4,5,6. Identifying rare variants that contribute to IOP will help uncover the biological mechanisms regulating this trait and provide improved understanding of IOP regulation and potential therapeutic targets for managing IOP and glaucoma.

Genome-wide association studies (GWAS) have identified over 190 genetic loci associated with IOP^1,2,3,4. These loci have established the contribution of common variants to IOP. The bivariate genetic correlation between IOP and glaucoma was also found to be high (0.49)⁷. While these studies identified numerous loci associated with IOP, these common variants typically show small effect sizes. The role of rare variants for IOP remains to be discovered. Rare variants typically require sequencing and a large sample size to have adequate statistical power.

The UK Biobank (UKB) is a large prospective cohort of half a million adult participants with extensive genetic data linked to physical measurements, health records, family history, and lifestyle information⁸. The recent release of whole-exome sequencing (WES) data now enables the exploration of rare variants for a variety of human traits and diseases and drug targets^9,10, including IOP and glaucoma. Rare variants can have large effect sizes and have demonstrated greater translational potential, e.g., PCSK9 as a target for lowering low-density lipoprotein levels^11,12 and MYOC as a target for gene therapy for treating myocilin-associated glaucoma¹³. These WES variants are also easier to interpret because they directly map to genes.

Using WES data from 110,260 participants in UKB, we conducted an exome-wide association study (ExWAS) to identify rare variants and genes associated with IOP, evaluated their effects on glaucoma in UKB and the FinnGen cohort, and explored potential drug targets of the identified genes. We also constructed a rare-variant polygenic risk score (rvPRS) and tested its association with glaucoma in independent white participants (n = 312,825). To the best of our knowledge, this study represents the largest rare-variant study of IOP to date. Our results uncovered rare variants regulating IOP, and subsequently, furthered our understanding of the biological mechanisms of IOP and potential drug targets for managing glaucoma.

Results

A total of 110,260 UKB participants were included in the IOP WES analysis, of which 98,674 were white. The mean (standard deviation [SD]) of age was 58 (8.1) years and 54% of the participants were female. The average IOP (SD) was 16.0 (3.4; range: 7.0–39.0) mmHg.

We identified 13 rare variants (10 of which are previously unreported) significantly associated with IOP, among which six were identified in white-only (white participants extracted based on a combination of self-reported White ethnicity [UKB data field 21000] and genetic information, see Methods for details) analysis and seven additional ones were identified in pan-ancestry (all ancestry combined) analysis. Table 1 displays the single-variant association results. Our top SNP, rs74315329 (P = 1.22 × 10⁻²⁶) is a well-known stop-gain variant in MYOC, the first gene identified for primary open-angle glaucoma (POAG)¹⁴. Consistently, rs28991009 in ANGPTL7 previously identified in our array-based GWAS² shows significance in this ExWAS using WES data. In white-only results, rs37278669, a nonsynonymous variant (allele frequency [AF] = 0.011%), in BOD1L1, shows a significant association with IOP (P = 5.75 × 10⁻⁹, beta = 4.08) in UKB. BOD1L1 is also significantly associated with the FinnGen phenotype “use of antiglaucoma preparations and miotics” (P = 7.7 × 10⁻⁶). The start-loss variant, rs753877638, in ACAD10 is significantly associated with both IOP (P = 1.30 × 10⁻¹⁰, beta = 8.41, AF = 0.003%) and glaucoma (P = 3.68 × 10⁻⁴) in UKB. In pan-ancestry analysis, rs201956837 in HLA-B is associated with IOP (P = 8.65 × 10⁻⁹, beta = 4.37). Rs201956837 is an intronic variant as well as an upstream transcript variant. The gene HLA-B is highly associated with glaucoma in FinnGen (P = 8.0 × 10⁻⁹). BOD1L1, RALYL, LDB3, ACAD10, CDK11A, and DPF3 are also associated with glaucoma topical treatments (P < 1 × 10⁻⁵, details in Supplementary Data 1). A Manhattan plot of the genome-wide P values for pan-ancestry results is shown in Fig. 1a. The genomic control lambda for white-only and pan-ancestry analyses are 1.01 and 1.02, respectively, which are well under control. The corresponding quantile-quantile plots are shown in Supplementary Fig. 1.

Table 1 Exome-wide significant rare variants for intraocular pressure

Fig. 1: Manhattan plots displaying the −log₁₀(P) for the association between IOP and rare variants and genes.

a Single-variant pan-ancestry results. The dotted horizontal line represents exome-wide significance associations (P < 1 × 10⁻⁸). b Gene-based pan-ancestry results. The dotted horizontal line represents gene-based significance associations (P < 2.5 × 10⁻⁶). Genetic variants or genes are plotted by genomic position. The colors on both plots show the delimitation for chromosome. Gene name is in black text if it has not been previously reported for intraocular pressure. Tests conducted in these analyses were two-sided and no adjustments were made for multiple comparisons.

From SAIGE-GENE analysis, 35 additional genes showed significant associations with IOP and 31 of them not previously published, among which 11 were identified from white-only analysis and 20 additional ones were identified from pan-ancestry analysis. Table 2 displays the gene-based association results. A Manhattan plot of the genome-wide p-values for pan-ancestry results is shown in Fig. 1b. Rare variants in previously known IOP genes, MTOR², EVA1C¹⁵, and CFAP298-TCP10L¹⁵, identified from common-variant investigations show significant gene-based associations with P = 1.08 × 10⁻¹², P = 9.51 × 10⁻¹⁰, and P = 1.34 × 10⁻⁸, respectively. Several of these ExWAS significant IOP genes, such as PTPRB, KIF21A, DNTT, also show a significant association with glaucoma in UKB with P = 3.26 × 10⁻⁵, 0.009, 0.007, respectively. Many of these IOP genes are associated with glaucoma-related traits in FinnGen. For example, CDCA8, HLA-B, RHOC, PPM1J, RPL10A, and TEAD3 are associated with glaucoma (P < 1 × 10⁻⁶). ADRB1, AAK1, IFI27, SYNGR3, and ZNF598 are associated with POAG (P < 1 × 10⁻⁵). Twelve of these genes, including PTPRB, HFM1, TAF1B, AAK1, FOXD1, EHMT1, and DNTT, are associated with glaucoma topical treatments (P < 1 × 10⁻⁵, details in Supplementary Data 1).

Table 2 Genome-wide significant results for intraocular pressure from gene-based analysis

To seek biological support for the identified genes, we evaluated their gene expression using both bulk RNA and single-cell RNA (scRNA) expression datasets. Supplementary Fig. 2 displays the bulk RNA expression information from Genevestigator¹⁶. A number of genes, such as BOD1L1, HLA-B, RPL10A, and RAB4B-EGLN2, are highly expressed in the trabecular meshwork (TM). Several other genes, e.g., ACAD10 and DNTT, show a medium gene expression in TM. Supplementary Figs. 3 and 4 display the scRNA expression information from the Cell atlas of the human ocular anterior segment (OAS)¹⁷ and of aqueous humor outflow pathways (AHOP)¹⁸, respectively. IOP and glaucoma related cell types can include TM fibroblasts, Schlemm canal endothelium (SCE), ciliary muscle (CM), corneal endothelium (CE), and vascular endothelium (VE)^17,18. Most of the identified genes show various levels of expression in these cell types. For example, BOD1L1 is expressed in all the above cell types; HLA-B and PTPRB are expressed in SCE and VE; and ACAD10 is expressed in CM and CE; and RALYL is expressed in CM, CE, and TM fibroblasts, to name a few.

To query potential drug targets, we used the Open Targets online resource. Table 3 displays the current known and proposed drug targets for these IOP rare-variant genes, such as ADRB1 (adrenoceptor beta 1, identified from our gene-based analysis), which is a known drug target for topical beta-adrenergic receptor antagonists, or beta-blockers, known to lower IOP. To the best of our knowledge, our results provided evidence for an unreported association between IOP and ADRB1. ADRB1 is expressed in human TM and ciliary body¹⁹, as well as cardiac tissue (Supplementary Fig. 2). Glaucoma drugs targeting ADRB1 include topical beta-blockers, such as timolol, betaxolol, carteolol, levobunolol, levobetaxolol, and metipranolol. Two older, outdated glaucoma medications include the adrenergic agonists, dipivefrin and epinephrine. In addition, several of these drugs are also used in treating hypertension and cardiovascular disease. PTPRB is highly expressed in the vein and artery endothelium cells (Supplementary Fig. 2). It is a proposed drug target for retinal vein occlusion, diabetic retinopathy and diabetic macular edema. Razuprotafib is a small molecule targeting PTPRB that acts as a negative regulator of Tie2 in diseased vascular endothelium by receptor-type tyrosine-protein phosphatase beta inhibition. EGLN2, a neighboring gene from the readthrough gene RAB4B-EGLN2, has drug trials for roxadustat, daprodustat, and vadadustat, which inhibit a hypoxia-inducible factor prolyl hydroxylase. These drugs target anemia and chronic kidney disease. MTOR is targeted by perhexiline, a drug used for cardiovascular disease that inhibits the serine/threonine-protein kinase mTOR. The MTOR gene is highly expressed in microvessel endothelium cells throughout the eye (Supplementary Fig. 2). RPL26 and RPL10A have three experimental drugs, i.e., ataluren, ELX-02, and MT-3724, two of which (ataluren and ELX-02) work as 80S ribosome modulators while MT-3724 functions as an 80S inhibitor. These drugs are in development for various diseases, such as cystic fibrosis, muscular dystrophy, hemophilia, epilepsy, kidney disease, and leukemia. Drug target genes ADRB1, PTPRB, and RPL10A, among others, were additionally found to have associations with vascular related phenotypes through PheWeb (Supplementary Data 2).

Table 3 Known drug targets for the identified intraocular pressure rare-variant genes

We further constructed a rare-variant polygenic risk score (rvPRS) using the IOP rare variants with P < 5 × 10⁻⁷ from pan-ancestry analysis (Supplementary Table 1) and tested its association with glaucoma in independent UKB white individuals (n = 312,825), who did not participate in the IOP measurements. This rvPRS is significantly associated with glaucoma with odds ratio (OR) per SD = 1.12 and P = 5.13 × 10⁻¹⁶, indicating the relevance of these IOP rare variants in glaucoma. When we used the rare variants identified from white-only subjects, the rvPRS yielded a mitigated association with glaucoma with OR per SD = 1.07 and P = 2.18 × 10⁻⁸. Since the IOP heritability explained by WES rare variants is less than 2% (estimated using GCTA^20,21), the overall prediction improvement over the baseline model (with only age and sex) in terms of the area under the receiver operating characteristic curve (AUC) of the current rvPRS is relatively low, 0.5%, in comparison to more than 5% in AUC improvement from common variants²², which can explain about 40% of the IOP heritability².

Discussion

In this study, we conducted the largest ExWAS of IOP to date using data from UKB. By employing single-variant and gene-based analyses, two complementary frameworks, we have expanded our knowledge of the genetic architecture of IOP, especially of the role of rare variants, beyond previous studies involving microarray data, which mainly covered common variants. In addition to confirming known IOP genes, we identified 40 previously unreported genes for IOP, demonstrating the power of including and aggregating rare variants in gene discovery. About half of these IOP genes are also associated with glaucoma phenotypes, including glaucoma medications, in UKB and FinnGen. Six of these genes are drug targets that are either established for clinical treatment or in clinical trials. Furthermore, we constructed a rvPRS and showed its significant association with glaucoma in independent white subjects.

We also showed that including subjects of all ancestries in a pan-ancestry analysis further improved the statistical power to identify rare variants. It was evident that pan-ancestry analyses identified additional rare variants and genes beyond white-only analyses in both single-variant and gene-based analyses. Testing these IOP variants and genes for their effects in glaucoma-related traits in both UKB and FinnGen and querying for drug targets further increased their translational relevance. Furthermore, the IOP rvPRS constructed using the rare variants identified from pan-ancestry analysis showed an even stronger association signal with glaucoma in independent white subjects than using white-only rare variants.

A concern with multi-ancestry datasets is false-positive signals. Numerous previous GWAS used European subjects only. In some studies, it was further reduced to unrelated European subjects. One way to analyze multi-ethnic GWAS datasets is using meta-analysis^1,23, which is typically used for dealing with common variants. However, rare variants may not have enough carriers in individual ancestral groups, resulting in too few carriers to be analyzed. A pooled approach is an attractive alternative for combining ancestrally diverse populations²⁴, especially for rare variants. Recent advances in statistical genetics tools also made this possible. For example, REGENIE²⁵, a machine learning approach, can avoid the parameter inflation in ultra-rare-variant situations while controlling for both population stratification and sample relatedness. SAIGE²⁶ uses mixed-effects models to adjust for both population stratification and genetic relationship matrix. Using these state-of-the-art methods, we showed the effectiveness of including non-European subjects in pan-ancestry analyses to further increase the study power. With advanced statistical genetics tools that can adjust for both genetic relatedness, principal components (PCs) of genetic ancestry and ancestral clusters, it is feasible to carry out pan-ancestry analyses in a pooled/combined approach. To further validate the robustness of our pan-ancestry results, we analyzed the significant rare variants in several sub-populations (UKB data field 21000), i.e. Black, Asian, and non-white all together, if there were sufficient allele counts for analysis, i.e., minor allele count (MAC) ≥ 5 (REGENIE default), in each sub-population. Supplementary Table 3 and Supplementary Data 3 show the sub-population-specific allele counts and association results, respectively. CDK11A rs556417493 is associated with IOP in Asian participants (AF = 0.08%, beta = 6.51, P = 7.59 × 10⁻⁹). PLK5 rs776910868 is associated with IOP in Black participants (AF = 0.1%, beta = 8.14, P = 8.42 × 10⁻¹⁰). These sub-population-specific results are highly consistent with our pan-ancestry results. ADSS2 rs375507039, PLAU rs367716060, and DPF3 rs933632776 are associated with IOP in non-white participants and could not be analyzed in separate individual sub-populations due to rare allele counts. For the HLA-B variant rs201956837, there is consistent direction of allele effects and effect sizes in white (beta = 4.33, P = 7.76 × 10⁻⁶) and non-white (beta = 4.47, P = 3.02 × 10⁻⁴). Using Fisher’s method to combine the two p-values (white and non-white), we obtained P = 2.35 × 10⁻⁹, which is consistent with our pan-ancestry result, P = 8.65 × 10⁻⁹, for the variant. The consistency of these results with our pan-ancestry results demonstrates the robustness of our analysis. Further research is warranted to maximize the power of pan-ancestry analysis²⁴.

Genetics provides vital information to identify drug targets. The generation of this WES data is sponsored by eight pharmaceutical companies, including Regeneron and AstraZeneca⁹, which clearly shows the value of this dataset to that industry. Drug candidates that have genetics support are twice as likely to be successful than those without genetics support²⁷. Six genes, i.e., ADRB1, PTPRB, RPL26, RPL10A, EGLN2, and MTOR, out of our gene-based analyses have existing therapeutic molecular targets. The most notable one, ADRB1, is the target of cardiovascular and glaucoma drugs, which include the broad class of glaucoma drugs targeting the beta-adrenergic receptor antagonists. For example, timolol was a first-line drug for lowering IOP by blocking the beta-adrenergic receptors in the ciliary body²⁸ to decrease aqueous humor flow²⁹. More recently, timolol has been shown to have an effect on outflow facility³⁰, which also impacts IOP. The other five genes are targets in many clinical trials involving razuprotafib, ataluren, ELX-02, MT-3724, roxadustat, daprodustat, vadadustat, and perhexiline, which provide candidates for drug repurposing for possible glaucoma treatment. For example, razuprotafib has been shown recently as an adjunct to latanoprost for treating glaucoma patients³¹. Razuprotafib also appears to stabilize blood vessels³¹. Roxadustat has proposed pathways affecting blood cell production³². Taken together, many of these drugs appear to be involved with cardiovascular disease and blood flow. Additionally, phenome-wide associations of the identified genes showed numerous significant associations with vascular-related phenotypes (Supplementary Data 2). Validations of cardiovascular relationships and drug targets for these IOP associated genes and recent success with drugs targeting vascular areas for glaucoma as seen with razuprotafib³¹ indicate that it may be possible to repurpose certain drugs that work on cardiovascular disease for glaucoma.

This study is not without limitations. Rare variants have their intrinsic challenges. The rarity of these variants makes their replication far more difficult than common variants. Nevertheless, since IOP is an endophenotype of glaucoma and ~70% glaucoma GWAS hits are also associated with IOP²³, it is reasonable to test these IOP hits for their glaucoma effects⁴ although it should not be assumed that all IOP hits are associated with glaucoma. Furthermore, IOP hits can also provide translation implications for glaucoma management since lowering IOP is currently the sole proven solution for glaucoma treatment. Hence, we checked the significance of these IOP variants and genes on glaucoma-related traits, including glaucoma treatment medication, in both UKB and FinnGen cohorts. In UKB, a combination of self-reported glaucoma and ICD-10/9 codes for glaucoma phenotypes is not homogeneous for specialized glaucoma subtypes, but previous studies have demonstrated the effect of using it for studying POAG genetics^4,33. Despite being the largest WES data currently available, diversity is still low, and European subjects comprise about 94% of the UKB cohort. Other larger ancestral groups are no doubt invaluable and can provide further information for discovery and validation. About half of the IOP rare-variant genes identified were found to be associated with glaucoma-related traits in either UKB or FinnGen. Further studies are required to confirm the remaining ones for their impacts on glaucoma. Furthermore, the best approach for analyzing datasets of ancestrally diverse populations remains an ongoing research topic²⁴, especially for rare variants. We used a combination of self-reported ethnic background and PCs of genetic ancestry to identify white participants who have similar ancestral backgrounds. Despite these efforts, subtle structure may still present. Sub-population labels for our analyses in Black and Asian in Supplementary Tables 3 and Supplementary Data 3 were extracted from self-reported ethnicity (UKB data field 21000). Self-reported ethnic background may be inaccurate for some individuals. However, we used state-of-the-art methods, i.e. REGENIE and SAIGE, which can handle both population structure and cryptic relatedness. There are many different approaches to analyze rare variants, e.g., grouping by predicted loss-of-function (pLOF) variants, missense variants, synonymous variants³⁴, all inclusive³⁵, and sliding window³⁶. To our knowledge, there is no consensus on the best approach to analyze rare variants. It can be phenotype dependent as well. Using other approaches may identify more rare variants and genes associated with IOP. Our rvPRS showed significant association with glaucoma, demonstrating the aggregated effect of the rare variants on glaucoma. However, its discriminatory ability in glaucoma prediction in terms of AUC is still low, which indicates that rare variants from WES may be more useful for biological insight than prediction at present. An exome only comprises about 1% of the human genome. Whole-genome sequencing data should be able to explain more IOP heritability. Hence, the best strategy to incorporate rare variants in PRS construction warrants further studies.

In conclusion, we carried out the largest ExWAS of IOP to date. In addition to showing the efficacy of single-variant and white-only analyses, our study clearly supports using gene-based aggregation and pan-ancestry analyses to further increase the study power. We demonstrated the value of rare variants to enhance our understanding of the biological mechanisms regulating this trait, and uncovered potential therapeutic targets for glaucoma.

Methods

UKB resource

UKB is an ongoing large prospective cohort study. Details regarding this cohort have been described elsewhere^37,38. Briefly, the UKB recruited over 500,000 adult participants (40–70 years of age at enrollment) living in the United Kingdom who were registered with the National Health Service at the study baseline (from 2006 to 2010). Medical information (self-report and electronic health records), family history, lifestyle information, as well as DNA samples, were collected. Ophthalmological data were also collected for a subset of study participants (~118,000). Most participants (~94%) reported their ethnic background as white and the rest originated outside of Europe⁸. UKB was approved by the North West Multi-Centre Research Ethics Committee and written informed consent was obtained from all participants. Our access to the resource was approved by UKB (application number 23424) and we obtained access to fully de-identified data.

FinnGen resource

The Finngen study is a large biobank study focused on the population of Finland³⁹. Over 200,000 participants have been enrolled, genotyped and phenotyped. 500,000 participants are projected to be enrolled by the end of 2023. The study aims to show the power of nationwide biobanks, electronic health records and an isolated population in identifying rare variants associated with different diseases. Data was collected from different Finnish biobanks and digital health care data on Finland citizens starting in 2017. The recruited population has an age average of 63 years and hospital-based recruitment predominates thus far. Phenotypes were built using the International Classification of Disease Ninth and Tenth Revision (ICD-9 and ICD-10) codes. Genotyping was done with a custom Axiom FinnGen1 and legacy arrays and further imputed to 17 million markers based on whole-genome sequences of Finns. Out of 2861 endpoint phenotypes created for this study, 15 are glaucoma related: neovascular glaucoma, primary angle-closure glaucoma, other and unspecified glaucoma, glaucoma, use of antiglaucoma preparations and miotics, juvenile open-angle glaucoma, normotensive glaucoma, glaucoma-related operations, primary open-angle glaucoma (strict), glaucoma (exfoliation), primary open-angle glaucoma, glaucoma secondary to other eye disorders, glaucoma secondary to eye inflammation, glaucoma secondary to eye trauma, and glaucoma suspect. The study used the SAIGE mixed models for their association analyses. The summary statistics are publicly online available (see data availability).

UKB WES and quality control

WES for all UKB participants were generated at the Regeneron Genetic Center^9,10. The sequencing, variant calling, and quality control were detailed previously^9,40. Briefly, sequencing was done on the Illumina NovaSeq 6000 platform using 75 base pair paired-end reads. Variant calling and quality control were performed using the SPB protocol⁴¹. The high-quality WES data have been reported to exceed 20× coverage at 95.8% of targeted bases. We overlapped the data with participants who participated in the ophthalmological measurements and kept all samples that had missing rate <2.5%. We kept autosomal variants with call rate >95% and minor allele count (MAC) ≥ 1 (15.1 million). We annotated these variants using VEP² and annovar⁴².

IOP measurements in UKB

IOP measurements were obtained using the Optical Response Analyzer (Reichert Corp., Philadelphia, PA) and have been described previously⁴³. In brief, both corneal-compensated and Goldman-correlated IOP measurements were collected. We used corneal-compensated IOP for this study since it is less affected by corneal thickness^44,45. The average of both eyes was used for downstream analysis. If only one IOP measurement was obtained, it was used as the final value. Study participants who received eye surgery within 4 weeks prior to the ocular assessment or those with possible eye infections did not receive IOP measurements. Moreover, we excluded study participants with extreme values of IOP, i.e., in the bottom and top 0.3 percentiles, and outliers, including participants who had either eye surgery or used eye drop medications^2,22. Overlapping with the WES data, 98,674 white participants (based on a combination of self-reported White ethnicity [UKB data field 21000] and genetic information; outliers with genetic ancestry at least six SDs from the means of the first two PCs were removed) and 110,260 pan-ancestry (all ancestry combined) participants remained.

Single-variant and gene-based ExWAS analyses

We performed single-variant association analyses using a machine-learning method implemented by REGENIE²⁵, accounting for population stratification and sample relatedness. We analyzed all variants with MAC ≥ 5 (REGENIE default) and minor allele frequency <1% and included age, sex, and the first 10 PCs as covariates. Genetic variants with P < 1 × 10⁻⁸ were declared ExWAS significant⁴⁶. In addition to using European participants, recent ExWAS studies advocate to include participants of all ancestries^47,48. Hence, we performed both white only and pan-ancestry analyses (added an additional covariate for four major ancestral groups, i.e., European, South Asian, East Asian, and African, identified by the K-Means clustering algorithm based on the first 10 PCs of genetic ancestry).

For gene-based association tests, we used SAIGE-GENE²⁶, a generalized mixed model approach that can adjust for both population stratification and genetic relationship. It performs rare-variant collapsing/aggregation tests, such as SKAT-O⁴⁹, burden⁵⁰ and SKAT^51,52. We used predicted loss of function (pLOF) variants as the variants for gene sets. We defined pLOF variants as: stop gained, stop lost, start lost, splice donor, splice acceptor and frameshift based on the VEP⁵³ annotation and gnomAD pLOF variants⁵⁴. We included age, sex, and the first 10 PCs as covariates. Genes with P < 2.5 × 10⁻⁶ were declared significant. We performed both white-only and pan-ancestry analyses (further added dummy variables for the major ancestral groups to the covariates).

Glaucoma lookup in UKB and FinnGen

Since lowering IOP is currently the only glaucoma treatment, we performed a lookup in glaucoma traits in UKB and FinnGen resources for all ExWAS significant IOP rare variants and genes. In the UKB participants, glaucoma cases were identified if they self-reported glaucoma (UKB data fields 6148, 20002) or had an ICD-10 or ICD-9 diagnosis code for glaucoma (UKB data fields 131186, 131188, 41202, 41204, 41076, 41078, 41270), excluding glaucoma secondary to eye trauma, secondary to eye inflammation, secondary to other eye disorders, secondary to drugs, and other glaucoma. The selection of glaucoma based on self-reports and ICD-10 codes has been shown to be effective in previous studies^4,33. Furthermore, the proportion of non-POAG cases in UKB was expected to be small⁵⁵. Controls were identified as those who did not have glaucoma or self-reported eye problems. Overlapping with WES data, 14,378 white cases and 409,571 white controls, 15,606 pan-ancestry cases and 437,417 pan-ancestry controls remained. We further checked our top IOP genes in three FinnGen GWAS summary statistics, i.e., glaucoma, POAG, and use of antiglaucoma preparations and miotics (Freeze 7)³⁹, by querying each of them in their online results (see Data availability).

Phenome-wide associations

For checking broader phenome-wide associations, we used PhenoScanner^56,57 and PheWeb⁵⁸. PhenoScanner consists of over 5000 genetic association datasets from NHGRI-EBI, NHLBI and UKB results. We performed a query of all IOP associated genes to generate associations with glaucoma topical treatment phenotypes (online query default cutoff P < 1.0 × 10⁻⁵, GWAS results source: http://www.nealelab.is/uk-biobank). Supplementary Data 1 shows details of the queried eye-related phenotypes from PhenoScanner, among which there are 15 unique glaucoma topical treatments. It has been reported that dry eye and glaucoma often occur together⁵⁹. Hence, significant artificial tear medication associations were also included. PheWeb uses summary statistics from the UKB to catalog millions of genetic markers across 1,403 binary traits. IOP associated genes queried in PheWeb generated a list of associations sorted by p-value. Out of these associations, phenotype traits related to the eye, cardiovascular, and nervous system were extracted. If no phenotype related to these traits were present, the association with the lowest p-value was reported.

Gene expression

We used Genevestigator¹⁶, a web-based gene expression database, to query bulk RNA information in different human tissues. Expression profiles of the queried genes in eye tissues from 210 human eyes were displayed in box plots showing the level of expression. For scRNAseq expression profiles, we used the Cell atlas of AHOP¹⁸ (queried through Spectacle⁶⁰) and of OAS¹⁷ (queried through the Broad Institute Single Cell Portal [see Web Resources]) online databases. AHOP and OAS data were generated from seven and eleven human samples, respectively^17,18. Each gene was queried to generate a heatmap and a violin plot displaying expression of various cell types related to AHOP and OAS of the eye, respectively.

Drug targets prioritization

To prioritize drug targets for the identified rare-variant genes, we used the Open Targets online resource. For the identified genes, we queried the Open Targets for known drugs, their mechanisms of action (source ChEMBL), and disease information. The druggable genes provide key information on the relevance of these genes on IOP and glaucoma management and potential drugs for repurposing.

rvPRS

From the pan-ancestry single-variant association results, we selected rare variants with P < 5 × 10⁻⁷ excluding intronic and synonymous variants. We assigned weights to these variants based on biological functions similar to that reported by Curtis^35,47. Details of these variants and their weights are shown in the Supplementary Table 1. We then constructed a weighted rvPRS using PLINK similar to our previous approach²², which was calculated as the summation of the number of rare risk alleles weighted by their biological functions. We then tested the association between the standardized rvPRS (subtracted the mean and divided by SD) and glaucoma in independent UKB white participants, who did not participate in the IOP measurements, using logistic regression adjusting for age and sex.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.