We used summary statistics from our recent GWAS meta-analysis (11). Briefly, we conducted a meta-analysis, including 538,515 individuals of European ancestry (59,944 cataract cases and 478,571 controls) from the GERA (25) and UK Biobank (26, 27) cohorts. The meta-analysis was performed using the R package “meta” (28) and fixed-effects summary estimates were calculated for an additive model. In total, 9,056,148 single nucleotide variations (SNVs) passing quality control were used for the TWAS analyses.

For the GERA cohort, all study procedures were approved by the Institutional Review Board of the Kaiser Permanente Northern California, and written informed consent was obtained from all participants. For the UK Biobank, this research has been conducted using the UK Biobank Resource project #14105.

Local (cis) eQTL datasets for 54 tissue types were downloaded from the FUSION website. These reference data were sourced from the Genotype Tissue Expression (GTEx) Project v8 (N=49 tissue reference panels) (20), the CommonMind Consortium (CMC) (N=2 tissue reference panels) (29), the Metabolic Syndrome in Men Study (METSIM) (N=1 tissue reference panel) (30), the Netherlands Twin Registry (NTR) (N=1 tissue reference panel) (31), and the Cardiovascular Risk in Young Finns Study (YFS) (N=1 tissue reference panel) (32). Supplementary Table S1 reports the datasets sources, number of individuals, and number of imputable genes for each tissue reference panel.

We conducted a TWAS of cataract using FUSION (18), which computes predictive models for eQTLs from reference data, and tests the association between predicted gene expression with a trait from GWAS summary statistics. As previously done (33), we performed tissue-specific TWAS analyses using FUSION default settings and the three following data inputs: 1) the above-mentioned GWAS summary statistics for cataract; 2) FUSION gene expression predictive models for 54 reference tissues; and 3) 1000 Genomes (European ancestry) Phase 3 data from the 1000 Genomes Project (34) as a reference panel for linkage disequilibrium (LD). Model weights for tissue-specific gene expression regressed on SNVs were computed from best linear unbiased predictor (BLUP), Bayesian sparse linear mixed model (BSLMM), least absolute shrinkage and selection operator (LASSO), and elastic net regression, as well as from the model with the top associated SNV.

A total of 325,513 gene-tissue-pairs (representing 37,920 unique genes across 54 tissue reference panels) were tested for associations between imputed gene expression with cataract susceptibility. Associations with a Bonferroni significance p-value less than 1.54 x 10^-7 (=0.05/325,513) were considered significant. Novel TWAS genes were defined as those located over 1 Mb apart from any previously reported cataract GWAS loci (i.e., no prior GWAS SNVs within 1 Mb from the start or end of the gene).

To assess whether GWAS SNVs colocalized with eQTLs, we conducted a Bayesian colocalization analysis using the COLOCv3.2.1 software, which is implemented in FUSION using marginal expression weights, for Bonferroni-significant TWAS associations (35). Thus, we tested the hypothesis that a single variant in each TWAS-significant model was associated with both cataract (from the GWAS) and imputed gene expression. Bayesian posterior probability greater than 0.9 was considered supporting evidence for colocalization.

To determine if the TWAS associations were conditionally independent of the GWAS hits, we conducted conditional analyses by adjusting transcriptome-wide associations for SNV-level effects from GWAS. Specifically, we used the COJO software program to adjust the GWAS summary statistics (the meta-analyzed results from the GERA and UK Biobank European samples) by the most statistically significant risk variants within 1 Mb of each TWAS gene (36). Using the marginal TWAS associations from the single-tissue analysis, we performed a FUSION joint analysis for cataract-associated genes located on the same chromosome region within each reference panel.

To identify tissues potentially relevant to cataract, we assigned the 54 tissue reference panels to 12 anatomical categories as per Strunz et al. (2020) (37): adipose (n = 3 reference panels), brain (n = 15), cardiovascular (n = 9), female reproductive (n = 3), gastrointestinal (n = 7), gland (n = 11), lung (n = 1), skeletal muscle (n = 1), skin (n = 2), tibial nerve (n = 1), and transformed fibroblasts (n = 1). Supplementary Table S1 lists the tissue reference panels and their corresponding anatomical categories used for this analysis. We assessed the frequency of Bonferroni-significant TWAS genes in each anatomical category. Because more Bonferroni-significant TWAS genes are expected from eQTL reference panels with more tissue donors and more imputable genes, we used the hypergeometric test to estimate the probability of observing at least as many TWAS-significant genes from all the gene-tissue pairs that we tested in each anatomical category.

We also conducted sex-specific TWAS analyses using sex-specific GWAS summary statistics (i.e., women and men analyzed separately) and tissue reference panels (i.e., ovary, uterus, and vagina eQTLs for women; and prostate and testis eQTLs for men).

We conducted an omnibus test in FUSION for associations with cataract across multiple tissues. Specifically, TWAS associations from all 54 tissue reference panels were jointly analyzed accounting for correlation between expression weights across tissues. Two filters were applied to the omnibus test results to consider a multi-tissue gene expression test significant: 1) using a Bonferroni correction, we divided the α of 0.05 by the effective number of genes tested (n = 13,328), and retained genes with omnibus test p-values less than this value (P < 3.75x10^-6); and 2) genes with a minimum tissue-specific p-value suggestive of a significant association (P < 1x10^-5) were retained as described by Barbeira et al (17).

The iSyTE 2.0 database, which contains meta-analyzed mouse lens gene expression data across different stages, was used to examine the expression of the novel genes identified in the current study in the lens tissue (21–23). Mouse orthologs of the human candidate genes for the novel cataract genes identified in the current TWAS analyses were examined in iSyTE, which contains meta-analyzed lens transcriptome data generated on microarrays or RNA-sequencing (RNA-seq) (21, 23). Mouse whole lens tissue gene expression datasets at embryonic day (E) stages E10.5, E11.5, E12.5, E16.5, E17.5, E19.5, and postnatal (P) day stages P0, P2, and P56, in addition to isolated lens epithelium at P28 were available on the Affymetrix 430 2.0 platform (GeneChip Mouse Genome 430 2.0 Array and/or 430A 2.0 Array) and were used in this analysis. Further, mouse whole lens tissue gene expression datasets at stages P4, P8, P12, P20, P30, P42, P52, and P60 were available on the Illumina platform (BeadChip MouseWG-6 v2.0 Expression arrays), and were used in this analysis. We also examined RNA-seq data generated on mouse whole lens tissue at E10.5, E12.5, E14.5 and E16.5. Additionally, because lens-enriched expression of a candidate gene has proven to be an effective predictor of its role in the lens (21, 22), the lens-enrichment of these candidate genes was also investigated at these stages. “Lens-enriched expression” is a measure of expression of a candidate gene in the lens compared to that in mouse whole embryonic body (WB) as described (21–23, 38, 39). Microarray expression data is publicly available on several gene-specific perturbation mouse models that exhibit lens defects or cataract, as described (21). We analyzed these datasets to examine potential changes in expression of the novel cataract candidate genes, as done before (11). Additionally, to gain insights into expression of candidate genes specifically in lens epithelial or fiber cells, we examined previously described RNA-seq data from isolated epithelium and fiber cells (40, 41). Gene expression analysis was performed as previously described (11, 23, 42). The University of Delaware Institutional Animal Care and Use Committee (IACUC) reviewed and approved the animal protocol.

We found that 202 gene-tissue pairs reached the Bonferroni significance threshold for their associations between imputed gene expression with cataract susceptibility ( Supplementary Table S2 ). While increased predicted expression was associated with cataract risk for 79 of the Bonferroni-significant gene-tissue pairs (e.g., IGHMBP2-colon sigmoid, z = 8.25), decreased predicted expression was associated with cataract for 123 Bonferroni-significant gene-tissue pairs (e.g., MRPL21 - whole blood, z = -7.64) ( Figures 1 , 2 ).

These 202 gene-tissue pairs were represented by 27 unique genes across 54 tissue reference panels ( Figures 1 , 2 ). Importantly, 2 of the 27 unique genes did not overlap previously identified cataract GWAS loci: INSRR on chromosome 1, and CEP95 on chromosome 17 ( Table 1 ). Furthermore, out of the 27 genes, we found that 9 (33.3%) were located in the 11q13.3 genomic region which was previously identified in our GWAS (11). These include: TESMIN, AP000808.1, MRPL21, IGHMBP2, MRGPRD, MRGPRF, AP003071.4, MRGPRF-AS1, and TPCN2 ( Figure 1 ).

Interestingly, 12 genes were Bonferroni-significant in only one tissue reference panel; these included 10 genes within previously reported cataract-associated loci: ARL4D (thyroid); CDC42BPA (heart atrial appendage); CDKN2A (brain cortex); GSTM2 (whole blood); OCA2 (brain); PKD2L1 (lung); SEMA4D (spleen); ST6GALNAC4 (blood); TESMIN (adipose subcutaneous); and TPCN2 (brain cerebellar hemisphere); and 2 genes newly identified in the current study: CEP95 (skin sun exposed lower leg), and INSRR (brain nucleus accumbens basal ganglia) Supplementary Table S2 .

To assess whether common genetic variants underly eQTL and GWAS associations with cataract, we conducted a colocalization analysis for the 202 Bonferroni-significant gene-tissue pairs. We found that 128 (63.4%) of the Bonferroni-significant gene-tissue pairs had a colocalized variant associated with both cataract risk (from GWAS) and predicted gene expression based on our TWAS results (column COLOC.PP4 in Supplementary Table S2 ).

To identify TWAS signals for cataract independent of GWAS genome-wide significant risk variants, we repeated the FUSION analysis with GWAS summary statistics conditioned on the top GWAS SNV in each of the 202 Bonferroni-significant gene-tissue pairs. We found that all gene-tissue pairs reached nominal significance (P < 0.05) ( Supplementary Table S2 ). Furthermore, we assessed joint TWAS associations in tissue reference panels with more than one Bonferroni-significant gene on the same chromosome region ( Supplementary Table S3 ). Of the four pairwise joint models including six unique genes, all the associations were attenuated but retained marginal significance (P<0.05). All six of these genes (IGHMBP2, TPCN2, MRPL21, MRGPRF-AS1, AP000808.1, and MRGPRD) are located within the 11q13.3 chromosome region, which has been previously identified as a GWAS susceptibility locus for cataract (11, 12).

Because cataracts are more common in women (4) and genetic susceptibility loci specific to women have been previously identified (11), we evaluated sex-specific TWAS associations. We used sex-specific GWAS summary statistics and tissue reference panels, i.e., GWAS summary statistics from women for TWAS of ovary, uterus, and vagina eQTLs; and GWAS summary statistics from men for TWAS of prostate and testis eQTLs. We found that 22 of the sex-specific tests reached the Bonferroni significance level that we applied to the main analysis (P<1.54x10^-7), including 3, 5, and 4 genes for ovary, uterus, and vagina, respectively, and 7 and 3 genes for prostate and testis, respectively ( Supplementary Table S4 ). Of those 22 sex-specific associations, 9 unique genes were identified, all of these genes (ITPKB, AC104162.1, AP000808.1, MRPL21, IGHMBP2, CAPRIN2, CLEC18A, LINC01229 and AC003681.1) were located nearby previously identified GWAS loci for cataract (11, 12). For instance, while differential gene expression of MRPL21 at 11q13.3 was associated with cataract in the 5 sex-specific tissue reference panels (i.e., ovary, uterus, vagina, prostate, and testis), differential gene expression of ITPKB was associated with cataract in vagina only.

Across the 54 tissue reference panels, the greatest number of Bonferroni-significant gene-tissue pairs was observed for the GTEx adipose subcutaneous, artery tibial, and thyroid datasets (seven gene-tissue pairs for each dataset), followed by the GTEx artery aorta, heart atrial appendage, esophagus muscularis, and skin sun exposed datasets (six gene-tissue pairs for each dataset) ( Supplementary Table S1 ). To identify tissues potentially relevant to cataract, we assigned tissue reference panels to anatomical categories as described above in the Methods ( Supplementary Table S1 ). In the tissue-specific TWAS results, gastrointestinal tissues represented one of the highest proportion of the 202 Bonferroni-significant genes (43 genes; 21.3%) ( Supplementary Figure S2 ). Interestingly, this gastrointestinal tissue type was the only anatomical category significantly enriched in our results, after accounting for the number of gene-tissue pairs tested per anatomical category (p-value from the hypergeometric test = 0.0055) ( Supplementary Table S5 ).

The multi-tissue TWAS using the omnibus test in FUSION revealed 86 genes for which imputed expression was associated with cataract susceptibility (Bonferroni p-value <0.05/13,328 effective gene tests ≈ 3.75x10^-6 and minimum tissue-specific p-value < 1 x 10^-5) ( Supplementary Table S6 ). Interestingly, 14 of the 86 multi-tissue associated genes were also associated with cataract in tissue-specific models, including INSRR (chr1), CDC42BPA (chr1), CDKN2A (chr9), CDKN2B (chr9), SEMA4D (chr9), ST6GALNAC4 (chr9), PKD2L1 (chr10), TESMIN (chr11), IGHMBP2 (chr11), BMP4 (chr14), OCA2 (chr15), TBX6 (chr16), RHOT1 (chr17) andARL4D (chr17). In addition, 18 of the 86 multi-tissue associated genes were located outside of previously described risk loci ( Table 1 ). These included: CROCC (chr1), ATP13A2 (chr1), CAPZB (chr1), SH3YL1 (chr2), EFEMP1 (chr2), CREB1 (chr2), ANKH (chr5), NIPSNAP2 (chr7), HERC4 (chr10), PICALM (chr11), COQ6 (chr14), ESRRB (chr14), SENP3 (chr17), GOSR2 (chr17), DDX5 (chr17), AC007773.1 (chr19), EIF2S2 (chr20), and ASIP (chr20).

We identified the mouse orthologs for 19 of the 20 novel genes as follows. For NIPSNAP2, in the Affy and Illumina microarray data, the gene alias for mouse Gbas was used. For AC007773.1, ZNF507 (mouse ortholog, Zfp507) and DPY19L3 (mouse ortholog, Dpy19l3) were considered as candidate genes. We first examined the expression of these genes in the lens tissue across various stages using the iSyTE microarray database (21, 22). While majority of the genes were found to be expressed, several exhibited robust expression ( Figure 3A ). For example, Atp13a2, Capzb, Crocc, Efemp1, Gbas, Gosr2, Picalm, Senp3 and Zfp507 had high expression in Affymetrix datasets. When examined for “enriched expression” in the lens, several candidates (e.g., Atp13a2, Capzb, Cep95, Crocc, Dpy19l3, Efemp1, Esrrb, Gbas, Gosr2, Insrr, Picalm, and Senp3) were identified ( Figure 3B ). Moreover, RNA-seq data from whole lens tissue confirmed 10 of the mouse orthologs to have expression or enriched expression in the lens ( Supplementary Figures S3, 4 ). Further, expression data from isolated lens epithelium and fiber cells at different time-points, spanning embryonic (E14.5 through newborn) through aging stages (3 months through age 2 years) showed that all the novel cataract candidate genes with mouse orthologs exhibit robust expression in the epithelium and/or fiber cells ( Supplementary Figure S5 ). Interestingly, this cell-specific data also shows that three candidate genes are specifically enriched in the postnatal epithelium (e.g., Efemp1, Esrrb, Insrr) and a subset of these exhibit progressively high expression with aging in the epithelium (e.g., Efemp1, Insrr). Finally, all novel candidate genes exhibited differential expression in at least one gene-perturbation mouse models of lens defects/cataract ( Supplementary Figure S6 ).