Genetic effects on liver chromatin accessibility identify disease regulatory variants

Abstract

Identifying the molecular mechanisms by which genome-wide association study (GWAS) loci influence traits remains challenging. Chromatin accessibility quantitative trait loci (caQTLs) help identify GWAS loci that may alter GWAS traits by modulating chromatin structure, but caQTLs have been identified in a limited set of human tissues. Here we mapped caQTLs in human liver tissue in 20 liver samples and identified 3,123 caQTLs. The caQTL variants are enriched in liver tissue promoter and enhancer states and frequently disrupt binding motifs of transcription factors expressed in liver. We predicted target genes for 861 caQTL peaks using proximity, chromatin interactions, correlation with promoter accessibility or gene expression, and colocalization with expression QTLs. Using GWAS signals for 19 liver function and/or cardiometabolic traits, we identified 110 colocalized caQTLs and GWAS signals, 56 of which contained a predicted caPeak target gene. At the LITAF LDL-cholesterol GWAS locus, we validated that a caQTL variant showed allelic differences in protein binding and transcriptional activity. These caQTLs contribute to the epigenomic characterization of human liver and help identify molecular mechanisms and genes at GWAS loci.

Keywords: ATAC-seq; GWAS; caQTL; cardiometabolic traits; chromatin accessibility; eQTL; transcription factor motif.

Figure 1

Joint profiling of gene expression and chromatin accessibility in human liver tissue (A) RNA-seq and ATAC-seq was performed in liver samples from 20 donors. (B) Distribution of consensus ATAC peak widths in base pairs. (C) Percent of consensus ATAC peaks by chromatin state in liver tissue from the Roadmap Epigenomics Project. All peaks, gray; 50,000 most accessible consensus peaks, black; quiescent represents unannotated regions. (D) Heritability enrichment of GWAS variants for multiple traits in all 223,265 liver ATAC peaks using stratified LD score regression. Points represent fold enrichment (proportion of heritability divided by proportion of SNPs within ATAC peaks) and error bars represent standard error. Significant enrichment (enrichment_p < 0.05), black; non-significant enrichment (enrichment_p > 0.05), gray. (E) Comparison of the distribution of expression between genes with and without an ATAC peak overlapping the transcription start site (TSS).

Figure 2

Identification and characterization of caQTLs (A) caQTLs identified using variants within 100 kb or 1 kb of peak centers. (B) Comparison of effect sizes between caQTLs and simple allelic imbalance (Pearson’s R = 0.75). The red line is the one-to-one line for caQTL effect sizes. (C) Comparison of effect sizes between caQTLs and H3K27ac QTLs (Pearson’s R = 0.40). The red line is the one-to-one line for caQTL effect sizes. (D) Comparison of the number of caPeaks and non-caPeaks assigned to each chromatin state in liver tissue from the Roadmap Epigenomics Project. caPeaks, purple; non-caPeaks, gray; quiescent represents unannotated regions. (E) Enrichment of caQTL variants in liver chromatin states. Error bars represent 95% confidence intervals. ^∗ indicates significant enrichment (p < 0.0071).

Figure 3

Disruption of TF binding motifs by caQTL variants (A) Allele affinities for TF binding and chromatin accessibility for variants within caPeaks and in strong LD with the caQTL lead variant (r² > 0.8). (B) Association of caQTL status with motif disruption status. Only the 109 TFs with at least 20 motifs disrupted by caQTL variants were included in the analysis, and only the 29 significant associations (p < 4.6 × 10⁻⁴) are shown. Error bars indicate 95% confidence intervals. (C) Percent of disrupted motifs for which the allele with higher chromatin accessibility matched the motif better. Percents are shown for the 29 TFs that had at least 20 motifs disrupted by caQTL variants. Black line, percent for all disrupted motifs across all tested TFs; red line, average percent across the 29 TFs.

Figure 4

Prediction of target genes for caPeaks using four approaches (A) Illustrations of four approaches to predict caPeak target genes. (B) Hi-C chromatin contact shown as an arc between caPeak191932 and the SNX10 promoter. Selected ATAC-seq signal tracks are shown for each caQTL genotype of rs12534816. More accessible homozygotes, purple; heterozygotes, black. (C) Genome browser image showing the correlation across rs12740374 genotypes of caPeak9372 and a peak at a SORT1 promoter. The purple arrow indicates the caPeak and the gray arrow indicates the promoter peak. (D) The same peak correlation with points representing normalized peak counts of individual samples colored by rs12740374 genotype. (E and F) SORT1 eQTL associations at the signal colocalized with the caQTL for caPeak9372 (E) and caQTL associations with caPeak9372 (F). In both plots, the caQTL lead variant within 1 kb of the peak center is indicated by a purple diamond and LD is based on 1000G phase 3 Europeans. (G) Comparison of directions of effect among all colocalized caQTL and eQTL signals. The A allele represents the more accessible allele than C, and more red marks indicate higher gene expression. (H) UpSet plot comparing the number of shared and unique caPeak-gene links identified by the four approaches. It is not possible for a caPeak-gene pair to be predicted using all four methods because if a caPeak is TSS proximal, it cannot form a Hi-C loop with the same gene and it cannot be a distal caPeak correlated with the promoter peak for the same gene.

Figure 5

A plausible regulatory mechanism at the EFHD1 locus for plasma liver enzyme levels (A–C) GWAS association with plasma levels of the liver enzyme alanine transaminase in Japanese individuals (A), eQTL association for EFHD1 (B), and caQTL associations for caPeak119621 (C). For all three plots, the caQTL lead variant within 1 kb of the peak center is indicated by a purple diamond and LD is based on 1000G phase 3 East Asians (A) or Europeans (B and C). Additional plots are shown in Figure S7. (D) Hi-C chromatin contact shown as an arc between caPeak119621 and the EFHD1 promoter. Selected ATAC-seq signal tracks are shown for each rs13395911 genotype. More accessible homozygotes, purple; heterozygotes, black; less accessible homozygote, gray. (E) Transcription factor ChIP-seq peaks in liver tissue from ENCODE that overlap caPeak119621. (F) Sequence logo plot for the FOXA2 motif s disrupted by caQTL variant rs13395911 (arrow). The motif match is shown on the negative strand, and variant alleles in (D) and (E) are shown on the positive strand.

Figure 6

Identification of a putative functional variant at the LITAF locus for LDL cholesterol (A and B) eQTL association for LITAF (A) and caQTL associations for caPeak75869 (B) at an LDL cholesterol GWAS signal. In both plots, the caQTL lead variant within 1 kb of the peak center is indicated by a purple diamond, and LD is based on 1000G phase 3 Europeans. Additional plots are shown in Figure S8. (C) Hi-C chromatin contact between caPeak75869 and the LITAF promoter. Selected ATAC signal tracks are shown for each rs57792815 genotype. More accessible homozygotes, purple; heterozygotes, black; less accessible homozygotes, gray. (D) Transcriptional activity of a 666-bp DNA element spanning caPeak75869 and containing rs3784924, rs11644920, and rs57792815 in HepG2 hepatocytes, THP-1 monocytes, THP-1 differentiated macrophages, and LPS-stimulated THP-1 macrophages. The DNA element was tested in the forward orientation relative to the genome (reverse orientation in Figure S8G). V, empty vector; H1, haplotype 1 of more accessible alleles rs3784924-A, rs11644920-A, and rs57792815-T; H2, haplotype 2 of less accessible alleles rs3784924-G, rs11644920-T, and rs57792815-C. Symbols represent 4–5 independent clones for each haplotype tested in duplicate wells; bars indicate mean ± standard deviation; p values from t tests of allelic differences. (E) EMSA using HepG2 nuclear extract (NE) shows allelic differences in protein binding for rs11644920. rs3784924 and rs57792815 are shown in Figure S8H. Green arrow, band represents T-allele-specific binding; black arrows, T-allele-preferential binding; white arrow, non-specific binding. Competition probes were unlabeled and included in 10-fold excess. (F) TF ChIP-seq peaks in liver tissue from ENCODE that overlap caPeak75869.