Integrating genetic regulation and single-cell expression with GWAS prioritizes causal genes and cell types for glaucoma

Abstract

Primary open-angle glaucoma (POAG), characterized by retinal ganglion cell death, is a leading cause of irreversible blindness worldwide. However, its molecular and cellular causes are not well understood. Elevated intraocular pressure (IOP) is a major risk factor, but many patients have normal IOP. Colocalization and Mendelian randomization analysis of >240 POAG and IOP genome-wide association study (GWAS) loci and overlapping expression and splicing quantitative trait loci (e/sQTLs) in 49 GTEx tissues and retina prioritizes causal genes for 60% of loci. These genes are enriched in pathways implicated in extracellular matrix organization, cell adhesion, and vascular development. Analysis of single-nucleus RNA-seq of glaucoma-relevant eye tissues reveals that the POAG and IOP colocalizing genes and genome-wide associations are enriched in specific cell types in the aqueous outflow pathways, retina, optic nerve head, peripapillary sclera, and choroid. This study nominates IOP-dependent and independent regulatory mechanisms, genes, and cell types that may contribute to POAG pathogenesis.

Fig. 1. Analysis workflow from POAG and IOP GWAS to causal regulatory mechanisms, genes, pathways, and cell types.

a POAG and IOP associations genome-wide (known and modest associations) were tested for enrichment among expression and splicing quantitative trait loci (e/sQTLs) in GTEx tissues and retina compared to permuted null sets of variants matched on confounding factors, using QTLEnrich (one-sided). In cases where enrichment was found, the lower bound number of e/sQTLs in a given tissue, likely to be true trait associations was estimated using an empirically derived, true positive rate ($\left({\pi }_1\right)$ approach. b Putative causal genes were prioritized per known POAG and IOP genome-wide association study (GWAS) locus by applying two colocalization methods (eCAVIAR, enloc) to all e/sQTLs from 49 GTEx tissues and retina eQTLs that overlapped each locus, followed by two-sample Mendelian Randomization (MR). Overlap of the colocalizing GWAS loci and e/sQTLs with Hi-C (3D chromosome conformation capture), cis-regulatory element (CRE), and super-enhancer (SE) regions from human retina was utilized to further prioritize causal genes. The human and eye images were created with BioRender.com. c All target genes of significantly colocalizing e/sQTLs (e/sGenes) or cell type-specific genes per trait were tested for enrichment in signaling and metabolic pathways (Reactome, KEGG), gene ontologies, and mouse phenotype ontologies using GeneEnrich (one-sided). The POAG cross-ancestry GWAS meta-analysis Manhattan plot was generated using QMplot (https://github.com/ShujiaHuang/qmplot). d Significantly colocalizing e/sGenes were tested for enrichment in specific cell types in single-nucleus RNA-seq data of glaucoma-relevant eye tissues, using ECLIPSER (one-sided). Cell type-specific genes were defined with cell type fold-change>1.3 and FDR < 0.1 per tissue. Cell type-specificity significance per GWAS locus set for a given trait was assessed against a null distribution of loci associated with unrelated, non-ocular traits, using a Bayesian Fisher’s exact test. Genes mapped to GWAS loci with a cell type-specificity score above the 95^th-percentile of null locus scores were proposed as contributing to the trait in the enriched cell type. e Cell type enrichment for the POAG and IOP GWAS was corroborated using two regression-based methods that assess cell type-specificity of trait associations considering all associations genome-wide: stratified-LD score regression (S-LDSC) and MAGMA.

Fig. 2. Enrichment and colocalization analysis of eQTLs and sQTLs with POAG and IOP associations.

Quantile-quantile (Q-Q) plots of POAG cross-ancestry (a) and IOP (b) GWAS -log₁₀ (P-value) compared to expectation for the best eQTL per eGene sets (eVariants with FDR < 0.05) of the most significantly enriched tissues based on adjusted fold-enrichment (colored points), compared to all variants in the GWAS (black points). Grey line represents the diagonal. QTLEnrich (one-sided test; see Methods) was applied to assess GWAS p-value enrichment among the eQTL sets. Bonferroni correction was applied. c Histogram of percent of GWAS loci with ≥1 colocalizing e/sQTL (colocalization posterior probability (CLPP) > 0.01 from eCAVIAR and/or regional colocalization probability (RCP) > 0.1 from enloc) for POAG cross-ancestry, POAG European (EUR) ancestry subset, and IOP European ancestry GWAS meta-analyses. Numbers above the bars represent the number of loci with at least one colocalizing e/sQTL. Red, dark blue, and light blue bars indicate percentage of loci with at least one colocalizing eGene, sGene, or both, respectively. e/sGene, gene with at least one significant e/sQTL. d Scatter plot comparing unique number of e/sGenes that significantly colocalized per GWAS locus versus unique number of e/sGenes tested per locus. Points are color-coded by number of GWAS loci. The black line represents the diagonal. e Violin plots showing the distribution of the unique number of colocalizing eGenes (red), sGenes (dark blue), or both (light blue) per locus for the three GWAS meta-analyses tested. The number of GWAS loci that colocalized with eGenes, sGenes, or both shown in the violin plot are N = 69, 36, 76, respectively, for POAG cross-ancestry, N = 34, 17, 36, respectively, for POAG EUR, and N = 76, 40, 79, respectively, for IOP (Supplementary Data 16). The center line in the box plots contained within each violin plot shows the median, the box edges depict the interquartile range (IQR), and whiskers mark 1.5x the IQR. The violin plot edges represent the minima and maxima values. f, Stacked histogram showing the number of colocalizing e/sGenes per gene biotype for each GWAS. Protein coding (light blue), lincRNA (brown), antisense (grey), pseudogenes (yellow), and other (dark blue).

Fig. 3. Colocalizing e/sQTLs in GTEx tissues and retina with top POAG GWAS loci.

Genes with at least one significant colocalization result are shown for e/sQTLs tested across 49 GTEx tissues and peripheral retina for the top 21 POAG cross-ancestry GWAS loci. GWAS loci were ordered by absolute value of their effect size. Within each locus, genes were ordered based on their chromosome position. Bubble size is proportional to the maximum colocalization posterior probability of all e/sVariants tested for the given gene, QTL type and tissue combination. Points are color-coded by direction of effect (blue if increased expression or splicing increases POAG risk or vice versa; red if increased expression or splicing decreases POAG risk or vice versa). Shape of points indicates colocalization method used: circle (eCAVIAR), triangle (enloc), and square (tested in both methods; results shown for method with maximum posterior probability). Grey or black border denotes variant-gene-tissue-QTL combination that passed quality control (QC) filtering (Methods) and a colocalization posterior probability cutoff above 0.01/0.1 (CLPP/RCP) or 0.5 (higher confidence), respectively. White or black asterisk in the square indicates whether the second method tested passed a posterior probability cutoff 0.01/0.1 (CLPP/RCP) or 0.5, respectively.

Fig. 4. Example of colocalizing e/sQTLs with a top POAG and IOP GWAS locus.

a Colocalization results for all eQTLs (e) and sQTLs (s) across 49 GTEx tissues and retina overlapping the POAG cross-ancestry rs2790053 LD interval, using eCAVIAR and enloc. Genes are ordered by chromosome position. Point size is proportional to maximum colocalization posterior probability (CLPP or RCP) of all e/sVariants tested per gene-QTL-tissue combination. Points are color-coded by direction of effect (blue if increased expression or splicing increases POAG risk or vice versa; red if increased expression or splicing decreases POAG risk or vice versa). Circle: eCAVIAR, triangle: enloc, and square: tested in both methods; results shown for method with maximum posterior probability. Grey or black border denote variant-gene-tissue-QTL combination that passed colocalization posterior probability above 0.01/0.1 (CLPP/RCP) or 0.5, respectively, and QC filtering (Methods). White or black asterisk in the square indicates whether the second method passed CLPP > 0.01/RCP > 0.1 or CLPP/RCP > 0.5, respectively. b LocusZoom plot for TMCO1 sQTL -log₁₀(P-value) in GTEx Cells-Cultured fibroblasts in POAG cross-ancestry GWAS rs2790053 (chr1_165768467_C_G_b38) LD interval. Points color-coded by LD (r²) relative to lead GWAS variant, chr1_165768467_C_G_b38. LocusCompare plots of POAG cross-ancestry GWAS meta-analysis -log₁₀(P-value) versus -log₁₀(P-value) of TMCO1 sQTL (c) or RP11-466F5.8 eQTL (d) in Cells-Cultured fibroblasts. Points color-coded based on LD (r²) relative to e/sVariant with highest CLPP. e Violin plot of normalized intron-excision ratio at chr1:165768269-165768482 (Leafcutter) for TMCO1 in fibroblasts versus genotype of sVariant chr1_165817713_G_A_b38 with highest CLPP for POAG cross-ancestry locus rs2790053. Boxplots represent median and interquartile range (IQR); whiskers mark 1.5x IQR; violin plot edges represent minima and maxima. f TMCO1 gene model (grey boxes: exons) in GTEx fibroblasts showing all intron excision splicing events detected with Leafcutter; zoom-in: alternative splice donor site events on exon 4 whose sQTL colocalized with POAG (red versus blue or green). g Exon boxes in TMCO1 gene model color-coded by exon read counts per base (blue) in GTEx Cells-Cultured fibroblasts; lines connecting exons for all splicing events are color-coded by exon-exon junction read counts (red). Transcript expression shown in Transcripts per Million (TPM) from RSEM. Plots are taken from https://gtexportal.org/.

Fig. 5. Chromatin loops and regulatory elements in human retina support effect of colocalizing e/sQTLs on POAG risk.

a Retina CREs (cyan) derived from epigenetic data overlapping e/sVariants that colocalized with POAG associations in the CDKN2A/B locus. The lead POAG variants from the cross-ancestry (blue line) and European subset (orange line) GWAS are shown in the top track, followed by their linkage disequilibrium (LD) proxy variants (r² > 0.8) in the track below. The significantly colocalizing CDKN2B-AS1 sVariants in Pituitary, CDKN2A eVariants in Brain Cortex, and CDKN2B eVariants in Skeletal Muscle are represented by red lines, and the grey lines represent LD proxy variants to the colocalizing e/sVariants, which are also significant e/sQTLs (FDR < 0.05) for the corresponding gene and tissue. b Retina CREs (cyan) overlapping retina SLC2A12 eVariants that colocalized with the POAG cross-ancestry association. Tracks display the lead POAG cross-ancestry GWAS variant rs2811688 (orange) and its LD proxy variants (grey), followed by the significantly colocalizing SLC2A12 retina eVariants (red) with their LD proxy eVariants, which are also significant eQTLs at FDR < 0.05 (grey). The CRE overlaps the promoter of SLC2A12. c Retina chromatin loops from Hi-C (3D chromosome conformation capture) data, SEs (blue), and CREs (cyan) shown for the RERE POAG locus. Tracks display the lead POAG cross-ancestry GWAS variant rs172531 (orange) and its LD proxy variants (grey), followed by significantly colocalizing RERE eVariants in Nerve Tibial, RERE sVariants in Cell-Cultured Fibroblast cells, and RERE-AS1 eVariants in Adipose Subcutaneous (red), and their LD proxy variants that are also significant e/sQTLs (FDR < 0.05) for the corresponding gene and tissue (grey). Magenta loops have one foot that overlaps or is in LD with the POAG variant and colocalizing e/sQTLs. In all panels, LD proxy variants were computed at r² > 0.8, TADs are represented as solid black lines, and the magenta heatmaps represent Hi-C physical contact maps. CRE Cis-regulatory element, SE Super-enhancer, TAD Topologically associating domain.

Fig. 6. Cell type enrichment of e/sQTL-mapped genes to POAG, IOP and related trait GWAS loci in the ocular anterior segment.

a Significance (circle size, -log10 (P-value)) and fold-enrichment (circle color) of the cell type specificity of e/sGene mapped-GWAS locus sets for POAG cross-ancestry, POAG European (EUR) subset, intraocular pressure (IOP), central cornea thickness, corneal hysteresis, physician-defined vertical-cup-to-disc ratio (VCDR), and machine learning-defined VCDR (VCDR ML), and as positive controls, eye color and cataract, using ECLIPSER (one-sided test), shown for 39 cell types in six ocular anterior segment tissues. Traits (rows) and cell types (columns) were clustered based on hierarchical clustering of Euclidean distance between the cell type-specificity enrichment scores of GWAS locus sets. Red rings: experiment-wide significant (Benjamini-Hochberg (BH) FDR < 0.1). Yellow rings: tissue-wide significant (BH FDR < 0.1). Grey rings: nominal significant (P < 0.05). b, d Cell type specificity fold-enrichment (x-axis) in the anterior segment cell types ranked in descending order for e/sQTL-mapped genes to GWAS loci of the POAG European subset (No. loci (N) = 37) (b) and POAG cross-ancestry (N = 79) (d). Points: fold-enrichment estimates from ECLIPSER. Error bars: 95% confidence intervals. Red: tissue-wide significant (BH FDR < 0.1). Grey: nominal significant (P < 0.05). Blue: non-significant (P ≥ 0.05). The exact p-values (one-sided) for the fold-enrichment of all cell types and GWAS are provided in Supplementary Data 35. c Differential gene expression (log₂(Fold-change), y axis) in the most strongly enriched cell type compared to all other cell types is shown for the set of genes (x axis) driving the enrichment signal of the POAG European GWAS loci in ciliary fibroblasts. The horizontal dashed line represents log₂(Fold-change) of 0.375 (FC = 1.3) and FDR < 0.1 used as the cell type-specificity enrichment cutoff. e Heatmap of fraction of genes that overlap between the e/sGenes driving the enrichment signal for top ranked cell types (P < 0.05) in the anterior segment for POAG cross-ancestry GWAS loci. Numbers refer to fraction of e/sGenes driving the cell type enrichment on each row that overlaps with the genes driving the cell type enrichment on the corresponding column. Hierarchical clustering was performed on both rows and columns using Euclidean distance between fractions. Cell type abbreviations are described in Supplementary Data 35.

Fig. 7. Cell type enrichment of e/sQTL-mapped genes to IOP GWAS loci in the anterior segment of the eye.

a Cell type specificity fold-enrichment (x-axis) in the anterior segment cell types ranked in descending order for e/sQTL-mapped genes to IOP GWAS loci (N = 37). Points: fold-enrichment estimates from ECLIPSER. Error bars: 95% confidence intervals. Red: tissue-wide significant (Benjamini-Hochberg (BH) FDR < 0.1). Grey: nominal significant (P < 0.05). Blue: non-significant (P ≥ 0.05). The exact p-values (one-sided) for the fold-enrichment of all cell types and GWAS are provided in Supplementary Data 35. b Heatmap of fraction of genes that overlap between the e/sGenes driving the enrichment signal for top ranked cell types (P < 0.05) in the anterior segment for IOP GWAS loci. Numbers refer to fraction of e/sGenes driving the cell type enrichment on each row that overlaps with the genes driving the cell type enrichment on the corresponding column. Hierarchical clustering was performed on both rows and columns using the Euclidean distance between fractions. Red asterisks: tissue-wide cell type enrichment BH FDR < 0.1 from ECLIPSER. c Differential gene expression (log₂(Fold-change), y axis) in the most strongly enriched cell type in the anterior segment compared to all other cell types is shown for the set of genes (x axis) driving the enrichment signal of IOP GWAS loci in pericytes (cluster 2). The horizontal dashed line represents log₂(Fold-change) of 0.375 (FC = 1.3) and FDR < 0.1 used as the cell type-specificity enrichment cutoff. Blue asterisks denote genes in IOP loci not associated with POAG risk. d Bubble map displaying the expression of e/sGenes driving the IOP enrichment in pericytes (red box) across all cell types in the anterior segment. The colorbar represents gene expression z-scores computed by comparing each gene’s average expression in a given cell type to its average expression across all cell types divided by the standard deviation of all cell type expression averages. Bubble size is proportional to the percentage of cells expressing the given gene (log(TPK + 1) > 1). Cell type abbreviations are described in Supplementary Data 35.

Fig. 8. Cell type enrichment of e/sQTL-mapped genes to POAG, IOP and related trait loci in retina, optic nerve head and surrounding posterior tissues.

Significance (circle size, -log10(P-value)) and fold-enrichment (circle color) of the cell type-specificity of e/sGene mapped-GWAS locus sets for POAG cross-ancestry, POAG European (EUR) subset, IOP, physician-defined vertical-cup-to-disc ratio (VCDR) or machine learning-defined VCDR (VCDR ML) GWAS, using ECLIPSER (one-sided), shown for all cell types in retina (a) and optic nerve head (ONH), optic nerve (ON), peripapillary sclera, sclera, and choroid (b). Traits (rows) and cell types (columns) were clustered based on hierarchical clustering of Euclidean distance between cell type-specificity enrichment scores of GWAS locus sets. Red rings: experiment-wide significant (Benjamini-Hochberg (BH) FDR < 0.1). Yellow rings: tissue-wide significant (BH FDR < 0.1). Grey rings: nominal significant (P < 0.05). Colored boxes highlight trait-specific enriched cell types. c Cell type-specificity fold-enrichment (x-axis) in ONH and surrounding tissue cell types ranked in descending order for the POAG cross-ancestry GWAS locus set (No. loci=79). Points: fold-enrichment estimates from ECLIPSER. Error bars: 95% confidence intervals. Exact p-values (one-sided) are provided in Supplementary Data 35. d Differential expression (log₂(Fold-change)) in astrocytes (6-Astro) in ONH compared to all other cell types in the posterior tissues is shown for the e/sGenes driving the POAG enrichment signal in astrocytes. Horizontal dashed line represents log₂(Fold-change) of 0.375 (FC = 1.3) and FDR < 0.1, used as the cell type-specificity enrichment cutoff. e Expression profile of e/sGenes driving the POAG cross-ancestry enrichment signal in astrocytes (red box) shown across all ONH and surrounding tissue cell types. Colorbar represents gene expression z-scores computed by comparing each gene’s average expression in a given cell type to its average expression across all types divided by the standard deviation of all cell type expression averages. Bubble size is proportional to percentage of cells expressing the gene (log(TPK + 1) > 1). f LocusCompare plot of -log₁₀(P-value) of POAG cross-ancestry GWAS meta-analysis relative to -log₁₀(P-value) of DGKG retina eQTL, which significantly colocalized with POAG locus rs56233426 (chr3_186411027_G_A). Points are color-coded based on LD (r²) relative to the eVariant with highest eCAVIAR colocalization posterior probability (CLPP = 0.93). Cell type abbreviations are described in Supplementary Data 35.