Human pancreatic islet three-dimensional chromatin architecture provides insights into the genetics of type 2 diabetes

Abstract

Genetic studies promise to provide insight into the molecular mechanisms underlying type 2 diabetes (T2D). Variants associated with T2D are often located in tissue-specific enhancer clusters or super-enhancers. So far, such domains have been defined through clustering of enhancers in linear genome maps rather than in three-dimensional (3D) space. Furthermore, their target genes are often unknown. We have created promoter capture Hi-C maps in human pancreatic islets. This linked diabetes-associated enhancers to their target genes, often located hundreds of kilobases away. It also revealed >1,300 groups of islet enhancers, super-enhancers and active promoters that form 3D hubs, some of which show coordinated glucose-dependent activity. We demonstrate that genetic variation in hubs impacts insulin secretion heritability, and show that hub annotations can be used for polygenic scores that predict T2D risk driven by islet regulatory variants. Human islet 3D chromatin architecture, therefore, provides a framework for interpretation of T2D genome-wide association study (GWAS) signals.

Figure 1. The promoter interactome of human pancreatic islets.

a, Overview of promoter-capture Hi-C (pcHi-C) in human islets. b, Integrative map of the KCNJ11-ABCC8 locus, showing human islet ATAC-seq and ChIP-seq, HindIII bait fragments, and arcs representing high-confidence pcHi-C interactions in human islets and erythroblasts. c, High-resolution annotations of islet open chromatin. ATAC-seq data from 13 islet samples was used to define consistent open chromatin regions, which were classified with k-medians clustering based on epigenomic features. Mediator and H3K27ac binding patterns allowed sub-classification of active enhancer classes I-III. Post-hoc analysis of islet CAGE tags confirmed that transcription start sites are highly enriched in promoters and weakly in class I enhancers. These islet regulome annotations are hereafter Supplementary Data Set 1. d, Average H3K27ac and Mediator signal centered on open chromatin regions for active enhancer subtypes in three human islet (HI) samples and input DNA. e, Overlap of promoter-interacting regions with epigenomic features, expressed as average log2 ratios (and 95% confidence intervals) over the overlaps obtained with 100 sets of distance-matched fragments. Error bars show s.d. across control sets.

Figure 2. Identification of target genes of islet enhancers.

a, We assigned target genes to 39.5% of all 45,683 active enhancers through high-confidence interactions. PAT features allowed imputing the assignment of promoters to another 40% of all active enhancers (see Supplementary Figure 3l,m for further details and evidence that imputed assigments are enriched in sub-threshold interactions). b, Functional correlation of enhancer-gene pairs assigned through high-confidence interactions (n = 18,637 pairs) or imputations (n = 28,695 pairs). Spearman’s Rho values for normalized H3K27ac signal in enhancer-promoter pairs across 14 human islet samples and 51 Roadmap Epigenomics tissues. Control enhancer-gene pairs were enhancers that overlapped a PAT in linear maps but were not assigned to the PAT promoter (n = 9,770 pairs), or other unassigned gene-enhancer pairs from the same TAD (n = 20,186 pairs). c, Concordance of enhancer eQTL-eGene pairs and enhancers-gene pairs assigned through high-confidence interactions (n = 351 pairs) or imputations (n = 293 pairs), relative to distance-matched control regions (n = 579 and 593 pairs, respectively), shown as a fold-change. P values were derived from one-sided Fisher’s exact test. d, Genes assigned to glucose-induced enhancers showed concordant glucose-induced expression. Top: glucose-induced enhancers showed enriched high-confidence (n = 439) or imputed (n = 640) assignments to glucose-induced genes, compared with distance-matched genes from the same TAD. Bottom: glucose-induced enhancers showed no enrichment for assignments to genes that were inhibited by high glucose concentrations (n=196 interacting and n=218 imputed pairs). OR = odds ratio. P values were calculated with Chi-square tests. e, Genes assigned to glucose-induced enhancers through high-confidence interactions (n= 275) or imputations (n=321 pairs) were enriched for glucose-induced promoter H3K27ac, compared with control genes from the same TAD. Box plots represent IQRs, notches are 95% confidence intervals of median, P values are from Wilcoxon’s two-sided signed ranked tests. See also Supplementary Data Set 2.

Figure 3. Identification of gene targets of T2D-relevant enhancers.

a, We assigned gene targets through high-confidence interactions or imputations for 53 (87%) out of 61 T2D-FG associated loci with genetic variants in islet enhancers (Supplementary Table 3). b, Summary of T2D-associated enhancer perturbations presented in this study (see also Supplementary Table 4). NT, not tested. c, Islet pcHi-C analysis defines gene targets of an enhancer bearing T2D-associated variants near CDC123/CAMK1D. The only T2D risk credible set variant that maps to an islet enhancer in the locus (rs11257655, zoomed inset) is assigned to CAMK1D and OPTN (dashed horizontal lines). Islet pcHi-C virtual 4C representations from pooled samples show interactions stemming from both CAMK1D and OPTN promoters towards rs11257655 with ChICAGO >3, but not from CDC123. d, CAMK1D and OPTN mRNA are regulated by the rs11257655-containing enhancer. We deleted the rs11257655-containing enhancer and a nearby control region with a T2D-associated variant (rs33932777) that lacked active chromatin marks in human islets. Cas9 only: n = 6 (2 independent experiments with triplicates). Deletions: n = 8 (2 gRNA pairs in 2 independent experiments with biological duplicates). Bars are means ± s.e.m., normalized by TBP and expressed relative to mean levels of the Cas9 only controls. Statistical significance: two-tailed Student's t-test.

Figure 4. Tissue-specific enhancer hubs regulate key islet genes.

a, Hubs are composed of one or more enhancer-rich PATs (≥ 3 class I enhancers) connected through at least one common interacting enhancer. Turquoise and dashed green lines depict high-confidence and imputed assignments, respectively. Descriptive features of hubs are summarized in Supplementary Figure 8c. b, Islet hubs are enriched in genes showing islet-selective expression. Ratios were calculated relative to all annotated genes. c, Islet hub genes are enriched in annotations important for islet differentiation, function and diabetes. Benjamini-Hochberg adjusted P values from EnrichR are shown (see complete lists in Supplementary Table 5). d, Gene pairs from the same hub show higher RNA correlations across human islet samples and 15 control tissues than gene pairs from the same TAD in which only one gene or neither gene is in a hub. P values were derived with Kruskall-Wallis analysis of variance. e, Enhancer-promoter pairs from the same hub show high H3K27ac correlations across 14 human islet samples and 51 Epigenome RoadMap tissues, compared with pairs from the same TAD in which only one element or neither are in a hub. P values were derived with a Kruskall-Wallis test. f,-g, Culture of 7 human islet donor samples at 4 vs. 11 mM glucose shows concerted changes in H3K27ac in hub enhancers connected with glucose-dependent genes. Hub promoters were ranked by their median fold-change in H3K27ac at high glucose, so that glucose-induced promoters are on the left of the X axis. (f) Median mRNA for genes associated with each hub. (g) Median glucose-dependent fold-change of H3K27ac in enhancers from hubs connecting with each promoter, IQR values in blue shade. In both graphs values are shown as running averages (window = 50). h, Coordinated glucose-induced H3K27ac in enhancers of a hub connected to KIRREL3. Top tracks show RNA and H2K27ac in one representative sample. Bottom insets highlight H2K27ac at 11 mM glucose (red) vs. 4 mM (blue) in regions showing coordinated glucose-induced changes in most hub enhancers, highlighted with black arrows (n = 4 human islet samples). See also Supplementary Table 6, Supplementary Data Set 5.

Figure 5. Tissue-specific topology of the ISL1 enhancer hub.

a, Epigenome annotations and high-confidence pcHi-C interactions from pooled islet samples and total B lymphocytes are shown to illustrate active enhancers, super-enhancers and enhancer clusters distributed across a TAD, while sharing islet-selective 3D interactions with ISL1 and HI-LNC57. b-c, 3D chromatin conformation models of the ISL1 enhancer hub generated from pcHi-C libraries from human islets (b) and total B lymphocytes (c). Images represent the top scoring model from the ensemble of structures that best satisfied spatial restraints. Class I, II and III enhancers are colored in dark to light red and promoters in blue if they are within 200 nm of the ISL1 promoter, or as white spheres if they are further than 200 nm. Note the proximity of lncRNA HI-LNC57 and ISL1 promoters in islets. The models show that active islet regulatory elements interact in a restricted 3D space in islet nuclei. See also Supplementary Figure 10b,c and Supplementary Videos 1 and 2.

Figure 6. The ZBED3 enhancer hub links an enhancer bearing a T2D SNP with multiple target genes.

a, pcHi-C and virtual 4C representations from pooled islet samples for three viewpoints (see also Supplementary Figure 10). The variant with highest posterior probability in this locus (rs7732130) maps to a class I islet enhancer (yellow line, and zoomed inset) that shows interactions with PDE8B (CHiCAGO > 5), and ZBED3, ZBED3-AS1, snoRA47 and S100Z (CHiCAGO > 3, see also Supplementary Figure 11). WDR41 is assigned to rs7732130 by imputation. Dashed horizontal lines show all targets assigned through imputation or high-confidence interactions. b, Analysis of hub and non-hub transcripts after CRISPR activation or inhibition of the transcriptional start site of ZBED3 or the rs7732130-enhancer in EndoC-βH3 cells. Data are presented as means ± s.e.m. of all gRNAs combined per target region (enhancer CRISPRa: 3 gRNAs, CRISPRi: 4 gRNAs, all n=3 independent experiments). Statistical significance: two-tailed Student's t-test.

Figure 7. Islet hub variants impact insulin secretion and provide tissue-specific risk scores.

a, Variant Set Enrichment (VSE) for T2D and FG (n=2,771 variants; Supplementary Table 9) and breast cancer (n=3,048 variants) in high-confidence interacting fragments in islets. Box plots show 500 permutations of matched random haplotype blocks. Red dots indicate significant enrichments (Bonferroni–adjusted P < 0.01). b, T2D and FG GWAS-significant variants are selectively enriched in hub class I islet enhancers. Boxplots show median and IQR. c, Genomic inflation of T2D association P values for non-GWAS significant variants (P > 5x10^-8) from a T2D GWAS meta-analysis (12,931 cases, 57,196 controls) in islet high-confidence interacting regions (magenta), non-interacting islet open chromatin (beige), and all other variants (brown). d, Genomic inflation of T2D association P values for non-GWAS significant variants in hub class I islet enhancers (blue), non-hub islet open chromatin (beige) and all other variants (brown). e, Heritability estimates based on GWAS summary statistics for T2D (12,931 cases, 57,196 controls), insulinogenic index (OGTT, 7,807 individuals), homeostasis model assessment of β-cell function (HOMA-B) and insulin resistance (HOMA-IR) (~80,000 individuals), for indicated islet enhancer domains. Bars show category-specific per-SNP heritability coefficients (τ_c) divided by LD score heritability (h²) of each trait. τ_c coefficients were obtained independently for each trait, controlling for 53 functional annotation categories. Values were multiplied by 10⁷ and shown with s.e.m. f, T2D frequency across 40 bins, each representing 2.5% of individuals in the UK Biobank test dataset (226,777 controls, 6,127 T2D cases) with increasing PRS, calculated with hub (pink dots) or genome-wide variants (light green). g, Odds ratios (OR) for T2D calculated for 2.5% individuals with highest PRS vs. all other individuals, using islet hub (pink) or genome-wide models (green), stratified by BMI and T2D age of onset. Boxplots show ORs for PRS from 100 permutations of pseudo-hubs (IQRs). Z-scores are standard deviations of pseudo-hub averages. See also Supplementary Figure 15 and Supplementary Table 17.