Genetic Control of Chromatin States in Humans Involves Local and Distal Chromosomal Interactions

Abstract

Deciphering the impact of genetic variants on gene regulation is fundamental to understanding human disease. Although gene regulation often involves long-range interactions, it is unknown to what extent non-coding genetic variants influence distal molecular phenotypes. Here, we integrate chromatin profiling for three histone marks in lymphoblastoid cell lines (LCLs) from 75 sequenced individuals with LCL-specific Hi-C and ChIA-PET-based chromatin contact maps to uncover one of the largest collections of local and distal histone quantitative trait loci (hQTLs). Distal QTLs are enriched within topologically associated domains and exhibit largely concordant variation of chromatin state coordinated by proximal and distal non-coding genetic variants. Histone QTLs are enriched for common variants associated with autoimmune diseases and enable identification of putative target genes of disease-associated variants from genome-wide association studies. These analyses provide insights into how genetic variation can affect human disease phenotypes by coordinated changes in chromatin at interacting regulatory elements.

Figure 1. Local QTLs

A) Number of local QTLs (10% FDR) for histone marks, RNA expression (Lappalainen et al., 2013) and DHS sites (Degner et al., 2012). Peaks are classified as promoters (“TSS”), “enhancer (Enh)” and “other” based on chromatin states (Kasowski et al., 2013). B)Chromosome-wide distribution of local QTLs on chromosome 1 for histone marks, DHSs and RNA. The vertical red line marks the location of the local joint QTL described in panels C-E. C) Signal tracks for three histone marks and DHS in a 5kb region around a joint local hQTL/dsQTL/eQTL coinciding with the ZNF695 promoter. The signal is aggregated across individuals by genotype at rs61373194. The position of the QTL SNP is indicated with a vertical dashed red line. D) Boxplots of aggregated signal for gene expression, histone marks and DHS grouped by the genotype of the QTL SNP. The normalized signal corresponds to the signal averaged across the entire peak region as indicated by black dashed lines in panel C. E) Position weight matrix for SPI1. The fourth position corresponds to the location of the motif altering SNP (rs61373194). The genotype with the strongest signal corresponds to a better match to the consensus sequence. Note: the motif is on the (-) strand. Therefore the genotypes are shown for the (-) strand to correspond to the PMW. F) Overlap between local eQTLs and hQTLs (calculated for promoters within 5kb of a histone peak). The overlap is highly significant (Fisher’s exact test); two thirds of all eQTLs are also hQTLs. There are 8,239 promoters that coincide with non-QTL histone peaks.

Figure 2. Evidence of Genetic Coordination Among Distal Functional Elements

A) Schematic of the four possibilities for genetic and physical coordination: the effect of a local QTL on a distal peak can be in the same (+/+ and −/−) or opposite direction (+/− and −/+) and the local QTL can be physically interacting with distal peak or not (bottom vs. top). B) The effect (beta of the regression) of a local QTL on its local peak (y-axis) and distal H3K27ac peak (x-axis) grouped by presence (“interacting”) or absence (“non interacting”) of Hi-C physical links (bottom vs. top). The direction of the effect tends to in the same direction for physically interacting pairs but not for non-interacting pairs (see also Fig. S2F for QTLs distal to the other marks, DHSs, and RNA). C) Enrichment of physical interaction among pairs of local QTLs and all distal genomic features (>50kb) that covary in the same direction with respect to the local QTL SNP. The enrichment is shown for all combinations of histone marks, RNA, and DHSs. (Fisher’s exact test, red bars indicate significance with *p<0.01 and **p<10⁻¹⁰). D) Enrichment for genetic associations (p-value < 10⁻⁶) in physically interacting regions as a function of the minimal distance between the SNP-distal peak pair. Enrichment and confidence intervals are calculated using Fisher’s exact test.

Figure 3. Distal QTLs

A) Number of distal QTLs identified with and without Hi-C. To identify distal QTLs we used SNPs that are a local QTL for any of the histone marks or RNA, split the set of SNP-peak pairs based on whether or not they are physically interacting (Hi-C correlation > 0.4), and calculated the FDR for each set individually (see EEP and Fig. S3A–D). B) Q-Q plot for the different sets of SNP-peak pairs for H3K27ac. The expected p-values were calculated by permuting sample labels (see EEP). Q-Q plots for the other marks as well as a control set using permuted Hi-C links is shown in Fig. S3A–B. C) Distance distribution for local-distal hQTL peak pairs identified using Hi-C interaction data and those identified using a standard hQTL approach. Most distal QTLs are less than 200kb apart from their targets. D) Distribution of absolute effect sizes of local and distal QTLs. Globally, local effects are stronger than the distal effects. E) Heat map of estimated Hi-C interaction counts for a region harboring multiple significant distal and local QTLs. The Hi-C interactions are based on a covariance method (EEP). Hi-C fragments with a correlation score > 0.4 were considered interacting. Local QTLs are indicated on the diagonal and off-diagonal dots indicate distal QTLs (>50kb). Black squares correspond to contact domain calls from (Rao et al., 2014). F) The position of distal QTLs was mirrored (see schematic) to test whether local-distal QTL pairs are more likely to reside within the same TAD. The schematic indicates how QTL SNPs were mirrored relative to the associated feature (histone peak, DHS, or RNA). Shown is the fraction of feature-QTL pairs that share the same TAD (orange) compared to feature-mirrored QTL pairs (grey). Features more frequently share the same TAD with the true QTL than the mirrored QTL. To avoid any bias we used the non-Hi-C aware set of QTLs for this analysis. G) Distribution of fractions of local-distal H3K4me1 QTLs sharing the same TAD for 100 sets of shuffled TADs. Real data corresponding to true TADs is indicated by the red vertical line. For chromatin marks, the fraction of local-distal hQTLs that share the same TAD is significantly higher than for the shuffled TADs (One-sided t-test). To avoid any bias we used the non-Hi-C aware set of QTLs for this analysis. The same analysis for the other marks is shown in Fig. S5D. H) Comparison of replication timing for local-distal QTL pairs. Shown are scatterplots of the replication timing for a H3K27ac peak region (y-axis) and the region of its distal QTL (x-axis; right panel) or mirrored QTL (see schematic in Fig. 4F) (x-axis; left panel). The peak-real QTL pairs show higher correlation than the peak-mirrored QTL pairs (Pearson correlation). I) The cumulative distribution of replication timing difference between peak and QTL regions (orange) and peak and mirrored QTL regions (grey) is shown. The differences for the hQTL- real peak region pairs are significantly smaller than for the corresponding mirrored positions (Wilcoxon rank sum test).

Figure 4. Characterization of Distal QTLs

A) Fraction of chromatin states for local-distal QTL pairs grouped by mark at distal site. Most pairs involve associations between enhancers and enhancers (yellow) or enhancers and TSSs (orange). B) The log odds ratio of observed vs. expected combinations of chromatin states for local and distal QTL pairs is shown for each chromatin mark, RNA and DHS. The number of expected combinations for each mark was obtained by permuting the state labels for the local peaks. (Fishers exact test * = p-value < 0.01) C+D) Number of local enhancer QTLs (C) and local promoter QTLs (D) influencing one or more distal promoters or enhancers (“out-degree”). In-degree plots are shown in Fig. S4. E) Out-degree (as in Figure 4BC) of local enhancer QTLs that are also eQTLs for one or more distal genes. F) In-degree of genes with one or more distal enhancer/promoter. G) Pearson correlation of expression signal among pairs of distal promoters that are associated with the same hQTL (one local, one distal). The correlations are compared to a set of distance matched permuted pairs (permuted genes) and a set of real links with permuted sample names (permuted individuals). Promoters sharing a hQTL are more correlated than permuted data (Wilcoxon rank sum test, p=2.8e-12).

Figure 5. Chromatin Loops and ChIA-PET

A) Example locus showing a genetic variant (rs4405472; dashed vertical line) that is a local hQTL for H3K4me1 and H3K4me3 as well as a distal hQTL for H3K4me3 at a promoter ∼ 200kb away, skipping several genes. The two loci are physically linked as indicated by significant ChIA-PET interaction calls (top) and a Hi-C based chromatin loop call (bottom, (Rao et al., 2014). The green line indicates a chromatin contact domain (i.e. TAD). B)ChIA-PET-links between one or more enhancers/promoters and distal genes based on interaction calls for H3K4me3 and Rad21. C) The log odds ratio of observed vs. expected state combinations for ChIA-PET linked peaks. Expected state combinations are calculated by exhaustively counting all possible combinations of states present at ChIA-PET peaks (Fisher’s Exact test, red=p<0.01) D) Enrichment of genetically associated (p value <10⁻⁶) local-distal hQTL pairs in interacting ChIA-PET fragments (Rad21) or chromatin loops as defined by Rao et al. for each histone mark and RNA. The 5% and 95% confidence intervals are shown in black (*p<0.01; Fisher’s exact test). E) Pairs of promoters that are physically linked (ChIA-PET) show correlation of signal for several features. All features except DHS show significantly higher correlation than a control of permuted, distance-matched links or permuted sample labels (Wilcoxon rank sum test) F) Same as G) but for enhancers linked by ChIA-PET. Interacting enhancers are more correlated than any of the permuted, distance-matched data (p value <10⁻¹⁶ for all histone marks and DHS; Wilcoxon rank sum test).

Figure 6. Transcription Factor Binding Analysis at Local and Distal hQTLs

A) Overlap enrichment of TF binding in H3K27ac peaks that have a QTL. For each TF, we plotted enrichment of having a QTL and being bound by the respective TF. Bars represent the 95% confidence interval. In red are significant enrichments, in gray non-significant ones. B) Overlap enrichment of TF binding in histone mark peaks, focusing on peaks with a local QTL (first three columns) or peaks affecting distal sites (last three columns). The enrichment value is only plotted for TF-histone mark pairs for which the enrichment was significant. Rows are sorted by the fold enrichment for peaks with a local H3K27ac QTL, which is displayed in A). C) Rank enrichment of TF binding in H3K27ac peaks. We plot the fold change enrichment of H3K27ac peaks in TF binding sites at increasing levels of significance for called hQTL peaks (red). The background enrichment was obtained by permuting the p-values between peaks (gray). D) Number of total peaks that can be explained through correlation between a TF motif score and the molecular phenotype. For each molecular phenotype, we computed the correlation between signal and motif score, and defined significant motif disruptions using permutations (5% FDR). The number of peaks with at least one significantly correlated TF motif within 2kb is shown. For each molecular phenotype, peaks were grouped by the number of motif-disrupting SNPs per peak. E) Fraction of H3K27ac hQTLs that are significantly correlated with TF motif disruptions. TF motifs show positive and negative correlation with the local histone mark signal and are sorted by the difference between percent positive and negative correlations. The total number of tested SNP-peak pairs across all H3K27ac peaks is annotated next to the TF name. Only TFs with N>=50 are shown. F) The signal surrounding H3K27ac QTNs was extracted and grouped into 6 clusters (pam clustering, EEP). The aggregate signals for the 6 clusters are shown for the high-, heterozygous- and low-genotypes (blue, purple, red) for H3K27ac, H3K4me3, H3K4me1, and DHS. Nucleosome positioning is indicated by MNase signal extracted from the same regions for a single individual (left to right; signal heat maps are shown in Fig. S6H). As expected, histone signal coincides with MNase signal / nucleosomes, whereas DNase hypersensitivity coincides with nucleosome-free regions. QTNs show concordant effects on the three histone marks.

Figure 7. Effect of regulatory elements and hQTLs on phenotypic diversity

A) Promoter and enhancer histone peaks with an hQTL are significantly enriched for GWAS SNPs compared to histone peaks without an hQTL. Dashed lines indicate the value for GWAS SNPs, while the null distributions indicate values by matching on GWAS SNPs for MAF, LD, and distance to TSS. B) Enrichment of hQTLs in GWAS SNPs. We plot the negative log10 adjusted p-value by log2 odds ratio of the Fisher’s Exact test. C–D) Rank enrichment of hQTLs in GWAS SNPs. At each tier of significance for GWAS variants we plot the fold change enrichment in overlap with hQTLs compared to base overlap of all GWAS variants at any significance level (red). We permute the p-values between the GWAS variants to generate background enrichments (gray). We select one positive enrichment and one negative enrichment example for each mark H3K27ac, H3K4me1 and H3K4me3. E) Regulatory network for Crohn’s disease derived from hQTLs and eQTLs intersected with GWAS SNPs. The network consists of 1) GWAS tag SNPs (gray boxes) connected to QTL SNPs (white nodes) through edges representing LD (gray edges), and 2) QTL SNPs connected to affected regulatory elements (orange triangles) or genes (green nodes) through edges representing either local QTLs (solid lines) or distal QTLs (dotted lines). An orange edge represents an hQTL/dsQTL, a green edge represents an eQTL. Distal QTL edges are labelled with the distance between the QTL SNP and the midpoint of the regulatory element. Regulatory elements are defined as the merged peaks from H3K4me1, H3K4me3, H3K27ac and DHS. If the affected peak in a regulatory element overlaps the transcription start site of a gene, we label the regulatory element with the gene’s name. If multiple QTL SNPs are associated with the same local regulatory element, we label the QTL SNP with the rsID of the most significant SNP; if multiple QTL SNPs are associated with the same distal regulatory element, then for every combination of local-distal association we pick the best-correlated SNP for the local peak. (Here, we only show the parts of the network that contain a gene. For the full network see Figure S7.)