QTL mapping of human retina DNA methylation identifies 87 gene-epigenome interactions in age-related macular degeneration

Abstract

DNA methylation (DNAm) provides a crucial epigenetic mark linking genetic variations to environmental influence. We analyzed array-based DNAm profiles of 160 human retinas with co-measured RNA-seq and > 8 million genetic variants, uncovering sites of genetic regulation in cis (37,453 mQTLs and 12,505 eQTLs) and 13,747 eQTMs (DNAm loci affecting gene expression), with over one-third specific to the retina. mQTLs and eQTMs show non-random distribution and enrichment of biological processes related to synapse, mitochondria, and catabolism. Summary data-based Mendelian randomization and colocalization analyses identify 87 target genes where methylation and gene-expression changes likely mediate the genotype effect on age-related macular degeneration (AMD). Integrated pathway analysis reveals epigenetic regulation of immune response and metabolism including the glutathione pathway and glycolysis. Our study thus defines key roles of genetic variations driving methylation changes, prioritizes epigenetic control of gene expression, and suggests frameworks for regulation of AMD pathology by genotype-environment interaction in retina.

Keywords: GWAS; aging; eQTL; eQTM; epigenome; gene expression; mQTL; neurodegeneration; retina.

Figure 1. Graphic summary of datasets generated, integrated and analyses performed in the present study and robust identification of retina mQTLs

a: Schematic representation of our genetic, epigenetic, and transcriptomic datasets and methods used in the identification and integration of methylation quantitative trait loci (cis-mQTL), expression quantitative trait loci (cis-eQTL) and expression quantitative trait methylation (cis-eQTM) with AMD GWAS and retina HiC chromatin map. b: Number of CpG sites tested (blue) and identified (grey) on various genomic regions in mQTL analysis. c: QTL enrichment in functional annotations for all retina (red) or retina-specific (orange) cis-mQTLs and all retina cis-eQTLs (blue). m/eQTL fold-enrichment on log2 scale with 95% confidence intervals is shown in descending order based on the retina mQTL fold-enrichment across annotations with >550 variants per QTL type. d: LocusZoom plot showing empirical GWAS association for the top mQTL signal with CpG (cg08027640) for PARK7gene. The diamond indicates the top mVariant for independent mQTL signal. The coloration of the points is determined by their LD with respect to the mQTL in purple. The top GWAS variant in the region is also labeled. The second plot shows −log10(P-values) of the variant association with CpGs in PARK7gene region from mQTL results. The grey and blue diamond’s represent −log10(P values) from mQTLs, respectively. e: Proportion of retina CpGs across different tissues in CpGs of significant mQTLs. f: Proportion of retina mQTLs across different tissues in significant mQTLs.

Figure 2. Characterization and distribution of eQTMs

a: Distribution of the distance between the CpG and the transcription start site (TSS) of the respective gene is plotted against the number of eQTMs. b: Combination chart representing the number of CpG sites tested (pink) and identified (yellow) on various genomic regions in eQTM analysis. c, d, e: DNAm levels are presented on the X-axis and the normalized gene expression levels are shown on the Y-axis. R2 is squared Pearson’s correlation between methylation and gene expression. c: eQTM for CpG cg24846343 located in gene body and GSTT2B on chromosome 22. d: eQTM for CpG cg21653793 located in 5’UTR and ABCA1 on chromosome 9. e: eQTM for CpG cg24307499 located in gene body and NLRP2 on chromosome 19. f: Distribution of eQTMs on different chromosomes based on their size. g: Top panel: Number of CpGs which regulate more than 10 eQTMs and are distributed on various chromosomes. Bottom panel: Cluster of CpGs on chromosomes 16 and 19 on arm p and q.

Figure 3. Associations between DNA methylation and gene expression through genotypes

a: Schematic of bidirectional integrative analysis that integrates summary-level data from independent GWAS with data from mQTL and eQTL. b: Heterogeneity in dependent instruments (HEIDI) model to distinguish pleiotropy from linkage for an observed association between DNAm and gene expression through genotypes. c: Manhattan plots of SMR tests for association between gene expression and DNAm (E2M_SMR). Shown on each y axis are the −log10 (P values) from SMR tests. The red horizontal lines represent the genome-wide significance level (PSMR = 9.29 × 10⁻⁷). d: Manhattan plots of SMR tests for association between DNAm and gene expression (M2E_SMR). The red horizontal lines represent the genome-wide significance level (PSMR = 9.38 × 10⁻⁷). e: Venn-diagram representing common and unique genes identified in E2M_SMR and M2E_SMR associations and bar graph representing the enriched pathways identified in the pathway analysis of common genes. f: Results of variants and SMR associations across DNAm and gene expression (M2E_SMR). The top plot shows −log10(P values) of SNPs from the SMR analysis of DNAm and gene expression (M2E_SMR). The blue diamonds represent −log10(P values) from SMR tests for associations of DNAm and gene expression. The second plot shows −log10(P values) of the SNP association for DNAm probe cg24506221 from the mQTL data. The third plot shows −log10(P values) of the SNP associations for gene expression of GSTM1 from the eQTL data. g: eQTM for CpG cg24506221 located in TSS200 region and GSTM1 on chromosome 1. DNAm levels are presented on the X-axis and the normalized gene expression levels are shown on the Y-axis. R2 is squared Pearson’s correlation between methylation and gene expression.

Figure 4. Associations between DNA methylation, gene expression and AMD GWAS through genotypes

a: Heterogeneity in dependent instruments (HEIDI) model to distinguish pleiotropy from linkage for an observed association between AMD, DNAm or gene expression through genotypes. b: Manhattan plot of SMR tests for association between DNAm and AMD GWAS. Shown on each y axis are the −log10 (P values) from SMR tests. The red horizontal lines represent the genome-wide significance level (PSMR = 5.67 × 10⁻⁷). c: Manhattan plot of SMR tests for association between gene expression and AMD GWAS. The red horizontal lines represent the genome-wide significance level (PSMR = 5.4 × 10−6). d: Results of variants and SMR associations across DNAm and AMD GWAS. The top plot shows −log10(P-values) of SNPs from the DNAm and AMD GWAS. The blue diamonds represent −log10(P-values) from SMR tests for associations of DNAm and AMD GWAS. The second plot shows −log10(P-values) of the SNP association for DNAm probe cg07160278 from the mQTL data. e: Results of variants and SMR associations across gene expression and AMD GWAS. The top plot shows −log10(P-values) of SNPs from the gene expression and AMD GWAS. The blue diamonds represent −log10(P-values) from SMR tests for associations of gene expression and AMD GWAS. The second plot shows −log10(P-values) of the SNP association for DXO gene from the eQTL data. f: LocusZoom plot between the AMD GWAS and the retina mQTL signals for cg11712338 (LINC01004). −log10(P-value) of AMD GWAS with points color coded based on LD (R²) relative to cg11712338 (LINC01004), with highest colocalization posterior probability (CLPP). g: LocusCompare plot comparing −log10(P-value) of AMD GWAS to −log10(P-value) of retina mQTLs acting on cg11712338 (LINC01004).

Figure 5. Colocalization analysis among AMD GWAS, mQTL and eQTL using Coloc and Moloc

a: Schematic representation of colocalization analysis of AMD GWAS, mQTL or eQTL. b: Schematic representation of multiple trait colocalization analysis of AMD GWAS, mQTL and eQTL. c: Circos plot depicting the multiple trait colocalization analysis of AMD GWAS, mQTL and eQTL. The orange lines represent the genes identified in associations between AMD GWAS and DNAm, the cyan lines represent the genes identified in associations between AMD GWAS and gene expression, the purple lines represent the genes identified in associations between AMD GWAS, gene expression and DNAm and the light blue lines represent the genes identified in associations between gene expression and DNAm. The highlighted genes are identified for associations between AMD GWAS, gene expression and DNAm. d: LocusZoom plots of GWAS (genotype and AMD GWAS association), CpG cg05475770 mQTLs (genotype and cg05475770 methylation association) TMEM259eQTLs (genotype and TMEM259 expression association). The y axis shows −log10(P value) of association tests from GWAS, mQTLs and eQTLs. Highlighted variant rs67538026 identified in Moloc analysis of AMD GWAS, mQTLs and eQTLs (GEM). e: LocusZoom plots of GWAS (genotype and AMD GWAS association), CpG cg21565421 mQTLs (genotype and cg21565421 methylation association) MPIeQTLs (genotype and MPI expression association). The y axis shows −log10(P value) of association tests from GWAS, mQTLs and eQTLs. Highlighted variant rs11072507 identified in Moloc analysis of AMD GWAS, mQTLs and eQTLs (GEM). f: Venn-diagram representing common and unique CpGs identified in associations between gene expression and DNAm using moloc (EM_coloc), CpGs identified in associations between gene expression and DNAm using SMR (E2M_SMR) and CpGs identified in associations between DNAm and gene expression using SMR (M2E_SMR). g: Venn-diagram representing common and unique genes identified in associations between gene expression and DNAm using coloc (EM_coloc), genes identified in associations between gene expression and DNAm using SMR (E2M_SMR) and genes identified in associations between DNAm and gene expression using SMR (M2E_SMR).

Figure 6. Hi-C data enables target gene and variant prioritization.

a: Proportion of unique mQTLs, eQTLs and eQTMs overlapping each chromatin compartment b: Number of unique mQTLs, eQTLs and eQTMs overlapping with a loop that are in contact with variants located within ±2.5 kb of the mGene/eGene/Gene TSS were identified as promoter mQTLs/eQTLs/eQTMs while those located >2.5 kb from the mGene/eGene/Gene were identified as distal mQTLs/eQTLs/eQTMs c: Upper panel: LocusZoom plots of CpG cg13422253 ALDH2eQTLs (genotype and ALDH2 expression association) and mQTLs (genotype and cg13422253 methylation association). The y axis shows −log10(P value) of association tests from eQTLs and mQTLs. Highlighted variant rs3858706 identified in E2M_SMR analysis of eQTLs and mQTLs. The lower panel includes associations identified in E2M_SMR for ALDH2 gene. Tracks represent the chromatin loops, E2M_SMR variant, SEs, CREs, H3K27Ac coverage, genes, and TADs. d: Upper panel: LocusZoom plots of CpG cg17020635 mQTLs (genotype and cg17020635 methylation association) and GSTP1 eQTLs (genotype and GSTP1expression association). The y axis shows −log10(P value) of association tests from mQTLs and eQTLs. Highlighted variant rs7108149 identified in M2E_SMR analysis of eQTLs and mQTLs. The lower panel includes associations identified in M2E_SMR for GSTP1 gene. Tracks represent the chromatin loops, M2E_SMR variant, SEs, CREs, H3K27Ac coverage, genes, and TADs.