Integrative analysis of liver-specific non-coding regulatory SNPs associated with the risk of coronary artery disease

Abstract

Genetic factors underlying coronary artery disease (CAD) have been widely studied using genome-wide association studies (GWASs). However, the functional understanding of the CAD loci has been limited by the fact that a majority of GWAS variants are located within non-coding regions with no functional role. High cholesterol and dysregulation of the liver metabolism such as non-alcoholic fatty liver disease confer an increased risk of CAD. Here, we studied the function of non-coding single-nucleotide polymorphisms in CAD GWAS loci located within liver-specific enhancer elements by identifying their potential target genes using liver cis-eQTL analysis and promoter Capture Hi-C in HepG2 cells. Altogether, 734 target genes were identified of which 121 exhibited correlations to liver-related traits. To identify potentially causal regulatory SNPs, the allele-specific enhancer activity was analyzed by (1) sequence-based computational predictions, (2) quantification of allele-specific transcription factor binding, and (3) STARR-seq massively parallel reporter assay. Altogether, our analysis identified 1,277 unique SNPs that display allele-specific regulatory activity. Among these, susceptibility enhancers near important cholesterol homeostasis genes (APOB, APOC1, APOE, and LIPA) were identified, suggesting that altered gene regulatory activity could represent another way by which genetic variation regulates serum lipoprotein levels. Using CRISPR-based perturbation, we demonstrate how the deletion/activation of a single enhancer leads to changes in the expression of many target genes located in a shared chromatin interaction domain. Our integrative genomics approach represents a comprehensive effort in identifying putative causal regulatory regions and target genes that could predispose to clinical manifestation of CAD by affecting liver function.

Keywords: CRISPR; GWAS; SNP; STARR-seq; cholesterol; coronary artery disease; enhancer; functional genomics; hepatocyte; liver.

Figure 1

CAD/MI SNPs are enriched in regulatory regions of hepatocytes (A) Flower plot depicting the percentage of CAD/MI GWAS SNPs that are also associated with type 2 diabetes (T2D), triglycerides (TG), high-density lipoproteins (HDL), total cholesterol, low-density lipoprotein (LDL), body mass index (BMI), basal metabolic rate (BMR), blood pressure (BP), and nonalcoholic fatty liver disease (NAFLD) in the UK Biobank. (B) Enrichment analysis of non-promoter regions in hepatocyte chromosomal interactions for enhancer- (H3K4me1-3, and H3K27ac) and repressor- (H3K27me3, H3K9me3) associated histone marks, and DNase I hypersensitive sites (DNaseHS). (C) Radar chart showing the enrichment of GWAS variants within HepG2 chromatin interactions. (D) Washu genome browser shot showing the location of SORT1 (chr1:109782257–109979272), H3K27ac ChIP-seq track for liver and HepG2, CAD-risk SNPs that fall within the looping ends, and PCHi-C interactions in HepG2 cells. Interacting restriction fragments are represented as boxes connected by a line on the HEPG2 PCHi-C track. (E) Gene ontology analysis of the target genes identified from PCHi-C data.

Figure 2

Identification of target genes regulated by CAD SNP harboring enhancers using PCHi-C and cis-eQTL analysis (A) Venn diagram showing the intersection of cis-eQTL genes and PCHi-C target genes. The names of the 25 most common genes are shown. (B) Expression of CAD GWAS SNP target genes (cis-eQTL and PCHi-C) in different cell types of the human liver based on scRNA-seq. A zoomed-in heatmap is shown for hepatocyte-specific genes.

Figure 3

Identification of regulatory regions influenced by inflammatory conditions in hepatocytes (A) Changes in H3K27ac signal and gene expression within enhancer hubs upon inflammatory stimulus. Hub promoters were ranked by their median fold change (FC) in H3K27ac upon cytokine treatment (2–23 h), so that inflammatory-induced promoters are on the left of the x axis. Similarly, median fold change in mRNA expression is shown for genes associated with each hub. (B) Gene pairs located within the same super enhancer or enhancer hub region show higher expression correlation across inflammatory treatment conditions than gene pairs from the random regions. p values were derived using the Kruskall-Wallis analysis of variance. (C and D) Examples of coordinated inflammation-induced H3K27ac at chromatin hubs encompassing IL1A, IL1B, MERTK, and CEBPD.

Figure 4

Association of target gene expression with liver-related traits Circos plot (circular manhattan plot) showing the association of target gene expressions with the KOBS data. Shown are (A) lipid traits (cholesterol, HDL, TG, LDL, and FFA) and (B) glucose and insulin levels in the KOBS cohort. Asterisk (^∗) denotes genes with a significant association also in the Hybrid Mouse Diversity Panel (HMDP). The height of the plot indicates the significance of the association (p value).

Figure 5

Analysis of allele-specific binding activity of CAD SNPs (A) Volcano plot of DeepSEA-predicted probability differences between the reference and alternative alleles for FOSL2 and CEBPB motifs that overlap CAD/MI GWAS alleles (lead+proxies). Two example SNPs, rs17293632 and rs9591145, predicted to result in high-significance allele-specific binding of the TF are shown. E-value is defined as the expected proportion of SNPs with larger predicted effect (from reference allele to alternative allele) for a given chromatin feature. (B) Heatmap of selected SNPs demonstrating allele-specific TF binding based on BaalChIP-analysis. See Figure S5 for a complete list. (C) A dot plot demonstrating negative correlation between the percentage of SNPs and the percentage of TFs demonstrating a change in TF binding. The x axis represents percent of SNPs that demonstrate significant allele-specific binding (>2-fold) for a given TF among all studied CAD/MI SNPs. The higher value indicates that among all SNPs demonstrating an allele-specific binding in BaalChIP, a specific TF is more likely to exhibit allele-specific binding. The y axis represents the percent from all studied TFs exhibiting an allele-specific binding across CAD/MI SNPs that exhibit allele-specific binding by at least one TF. The higher value indicates that the allele-specific binding of a specific TF could be more likely to translate to an allele-specific binding by the other studied TFs. (D) EMSA for the rs9591145 SNP showing that the “T” allele significantly gains binding affinity of CEBPB compared to the “G” allele.

Figure 6

Allele-specific activity of enhancers investigated by STARR-seq (A) Dot plot depicting the STARR-seq input DNA library ref-allele proportions in relation to the experimentally quantified RNA library ref-allele proportions in HepG2 cells. SNPs with significant allele-specific effects are highlighted in blue. Selected genes associated with the studied enhancer regions are shown. (B–D) Heatmap of selected enhancers showing a significant allele-specific difference compared to the most common haplotype (1) observed in the European population for which the target genes were defined by (B) PCHi-C, (C) proximity, and (D) cis-eQTL analysis. Asterisk (^∗) represents the SNP that demonstrated allele-specific binding in BaalChIP while hatch mark (#) denotes SNPs that were predicted to disrupt TF binding by DeepSEA and double s (§) is an experimentally validated location from Figure 7.

Figure 7

CRISPR-mediated genetic perturbations of enhancer hubs (A and B) Washu genome browser shots of two enhancer hubs containing (A) TSPAN14 and (B) SFXN2 where the CRISPR-Cas9 system was used to delete the enhancer region harboring a CAD GWAS variant. GRO-seq showing enhancer RNA (eRNA) for HEPG2 comes from GSE92375 (GSM2428726). (C and D) Analysis of the effect of enhancer deletion on gene expression within the TSPAN14 and SFXN2 hubs in HepG2 cells. qPCR was performed for genes located in the same hub as well as for genes in adjacent hubs. (E) Analysis of the effect of CRISPRa-mediated activation of enhancer variants in six selected chromatin hubs. For locus information, see Figure S9. Gene expression data are presented as the mean ± SEM of three independent experiments. The statistical significance was evaluated using a two-tailed Student’s t test or Mann-Whitney U test. For all bar plots, significance is denoted with asterisk. ^∗p < 0.05, ^∗∗p < 0.005, and ^∗∗∗p < 0.0005.