The impact of proinflammatory cytokines on the β-cell regulatory landscape provides insights into the genetics of type 1 diabetes

Abstract

The early stages of type 1 diabetes (T1D) are characterized by local autoimmune inflammation and progressive loss of insulin-producing pancreatic β cells. Here we show that exposure to proinflammatory cytokines reveals a marked plasticity of the β-cell regulatory landscape. We expand the repertoire of human islet regulatory elements by mapping stimulus-responsive enhancers linked to changes in the β-cell transcriptome, proteome and three-dimensional chromatin structure. Our data indicate that the β-cell response to cytokines is mediated by the induction of new regulatory regions as well as the activation of primed regulatory elements prebound by islet-specific transcription factors. We find that T1D-associated loci are enriched with newly mapped cis-regulatory regions and identify T1D-associated variants disrupting cytokine-responsive enhancer activity in human β cells. Our study illustrates how β cells respond to a proinflammatory environment and implicate a role for stimulus response islet enhancers in T1D.

Extended Data Fig. 1. Chromatin characterization of human pancreatic β cells exposed to pro-inflammatory cytokines.

a, Pearson correlation values between replicates in different assays and conditions (see Supplementary Note 2). b, Volcano plots of ATAC-seq (left) and H3K27ac ChIP-seq (right) changes obtained after exposure of EndoC-βH1 to IFN-γ and IL-1β; green and red dots correspond to sites with ∣log₂ FC∣>1 and FDR adjusted P<0.05 as calculated by fitting a negative binomial model in DESeq2. Chromatin changes are classified as “gained” and “lost” chromatin sites whereas non-significant changes are defined as “stable”. c, Chromatin accessibility and H3K27ac enrichment changes observed in EndoC-βH1 are largely replicated in human pancreatic islets as illustrated by the distribution of log₂ fold change at regions as classified in b in EndoC-βH1. Dotted lines indicate log₂ fold change thresholds (∣log₂ FC ∣>1). Box plot limits show upper and lower quartiles, whiskers extend to 1.5 times the interquartile range and the notch represents the confidence interval around the median. d, Hierarchical clustering using normalized ATAC-seq and H3K27ac ChIP-seq read counts at EndoC-βH1 IREs shows that samples cluster according to treatment, suggesting that the differences caused by the proinflammatory cytokines are greater than those derived by the sample heterogeneity. HI=Human pancreatic islets, EndoC=EndoC-βH1 e, Distribution of distances to nearest TSS for the different types of regulatory elements, showing that IREs, compared with stable regulatory elements (SREs), are preferentially located distally to TSS. f, Mean sequence conservation score of IREs and a randomized set of IREs in placental mammals. Peaks were extended from the center 1kb to each direction and mean score was calculated in 50bp windows. g, Sequence composition analysis of IREs (n=3,009) illustrating the top identified de novo motifs. Colors for matched genes correspond to RNA-seq (name) or protein (underlined) status (red=down-regulated, blue=equal-regulated, green=up-regulated, black/no line= not expressed/detected).

Extended Data Fig. 2. Exposure to pro-inflammatory cytokines drives changes in the transcriptome and proteome of pancreatic β cells.

a, Volcano plot of RNA-seq genes, showing up-regulated genes (green) and down-regulated genes (red) upon exposure of EndoC-βH1 to cytokines. Vertical lines indicate the log₂ fold change threshold (∣log₂ FC∣>1) and horizontal line indicates the FDR adjusted P cutoff for significance (FDR adjusted P<0.05) calculated by fitting a negative binomial model in DESeq2. b, Distribution of RNA-seq counts in human islet samples in the genes previously classified as up, down or equal-regulated in EndoC-βH1 cells. Boxplot limits show upper and lower quartiles, whiskers extend to 1.5 times the interquartile range and the notch represents the confidence interval around the median. c, Volcano plot for multiplex proteomics, showing in green the up-regulated proteins and in red the down-regulated, which have a Q-value<0.1 and ∣log₂ FC∣>0.58. Vertical lines indicate the log₂ fold change thresholds. d, Protein-protein Interaction (PPI) network generated from up-regulated proteins after cytokine exposure. Node color indicates belonging to same interacting community and background corresponds to specific pathway enrichment. e, Proportion of up, equal or down-regulated proteins encoded by genes located <15 kb from IREs or SREs. *** Chi-squared test P<0.001. f, An additive effect on gene up-regulation was observed for multiple IREs located at <40kb of a gene. Box plot limits show upper and lower quartiles, whiskers extend to 1.5 times the interquartile range and the notch represents the confidence interval around the median. ANOVA P<2.2×10⁻¹⁶. g, View of the LY6E locus, whose expression is induced after cytokine exposure and is coupled with chromatin changes in the vicinity.

Extended Data Fig. 3. Characterization of β-cell IREs.

a, Genes associated to different classes of IREs (classified as in Fig. 2a) show cytokine-induced expression in EndoC-βH1. CYT=cytokine exposed, CTRL=control. Boxplot limits show upper and lower quartiles, whiskers extend to 1.5 times the interquartile range and notch represents the median confidence interval. ***Wilcoxon test P<0.001. b, Sequence conservation score of IREs and a corresponding randomized set used as control. c, Distribution of distances to nearest TSS of the different classes of open chromatin sites. Line indicates the threshold used to classify them as “promoters”. d, Number of IREs overlapping regions annotated as “Strong” or “Weak” enhancers by ENCODE ChromHMM. *Chi-squared P<2×10^-16. e, f, Top hits for de novo motif analysis in opening (e) and primed enhancers (f). Colors for matched genes correspond to RNA-seq (name) or protein (underlined) status (red=down-regulated, blue=equal-regulated, green=up-regulated, black/no-line=not-expressed/detected). g, Diagram showing the percentage of colocalization between the TF binding sites identified by de novo motif analysis in SRE and primed enhancers (i.e. excluding sites <2Kb from a TSS). Label size indicates number of regions containing the TF binding sites and line width/intensity percentage of regions in which two motifs colocalize. h, Odds-ratio for finding a motif pair in the same enhancer in primed vs. SRE. Only significant pairs (FDR-adjusted Fisher Exact test P<0.001) are shown. Immune and islet-specific TF motifs colocalize more often in primed compared to SRE chromatin sites. i, Percentage of overlap between EndoC-βH1 different classes of open chromatin and islet-specific TFs obtained by ChIP-seq in untreated human islets. j, Volcano plot showing differentially methylated sites (depicted in red) in EndoC-βH1 exposed or not to cytokines. Dotted lines indicate the threshold for methylation differences or significance using limma moderated t-test. k, Distribution of demethylated and stable CpGs according to different classes of open chromatin.

Extended Data Fig. 4. Deconstructing cytokine induced cis-regulatory networks in β cells.

a, Gene Regulatory Network (GRN) derived from IREs and their putative target genes. Squares represent the IREs inferred TF binding sites (motifs logos and TF matches are shown on the right side) and the ellipses represent their putative target genes (see methods). The size of the squares reflects the number of connections (edge count) while the gene node size reflects the log₂ fold change of RNA expression after cytokine exposure. The resulting GRN is an interconnected scale-free network composed of 648 nodes and 3,589 edges. Genes regulated exclusively by primed IREs are represented in blue while green depicts opening IREs regulated genes. Red denotes genes regulated by both types of IREs. In each of these three groups the representation of the hierarchy is based on the principle of network centrality where authoritative nodes are located more proximal to the core. b, Comparison between the degree distribution of the observed GRN (black triangles) and a random generated network (blue squares) having the same number of nodes and edges. The bell-shaped degree distribution of random graph denotes a statistically homogeneity in the degree of small and large nodes. In contrast, the observed network showed a high frequency of small degree nodes and a low frequency of highly connected nodes as is typical of a scale-free network. c, Bar plot of gene ontology biological process enrichment analysis. Gene-ontology analysis was performed using all target genes in the GRN. Functional enrichment analysis was performed by Metascape (http://metascape.org). Only terms with P<0.001 and with at least 3 enriched genes were considered as significant. Color is proportional to their P values.

Extended Data Fig. 5. 3D chromatin changes induced by exposure of human islets to pro-inflammatory cytokines.

a, Violin plots showing the distribution of CHiCAGO scores of contacts, detected by pc-HiC experiments in untreated human islets, between stable and induced enhancers and their target genes. SREs engage chromatin contacts with higher interaction scores compared to those detected for IREs. *** Wilcoxon test P<0.001. b, c, d, Views of the 3D chromatin contacts of CIITA (b), SOCS1 (c) and RSAD2 (d) promoters obtained by UMI-4C performed in islets exposed or not to pro-inflammatory cytokines. In yellow we highlight those IREs that gain contacts with the up-regulated gene promoter. A heatmap under the 4C track represents the log₁₀ odds ratio (OR) of the UMI-4C contacts difference in cytokine vs. control and a small black diamond on top of the contact heatmap indicates a significant difference in contacts between cytokine-treated and control samples (Chi-squared P<0.05). ATAC-seq peaks are represented by rectangles, shaded from gray to green proportionally to the cytokine-induced H2K27ac log₂ fold change observed at that site.

Extended Data Fig. 6. Cytokine-induced islet regulatory elements are enriched in T1D associated variants.

a, EndoC-βH1 cytokine-induced regulatory elements (IREs) overlap more often than expected T1D associated variants while the opposite is true for T2D. EndoC-βH1 cytokine-invariant regulatory elements (SREs) are instead enriched for T2D, but not T1D associated variants. Each dot denotes the Varian Set Enrichment (VSE) score in IREs or SREs regions. Boxplot shows the enrichment distribution of the matched null permutated data sets. Red dots indicate that the difference is statistically significant as determined by VSE (Bonferroni adjusted P < 0.05). Box plot limits show upper and lower quartiles, whiskers extend to 1.5 times the interquartile range and the notch represents the confidence interval around the median. b-f, Representative regional plots of different T1D risk loci containing T1D variants overlapping IREs and up-regulated genes. R² values are based on 1KG phase 3 EUR and the leading SNPs in the locus is represented by a diamond. If different leading variants are present in the same locus, their proxies are depicted in different colors. Yellow squares highlight those variants that overlap a human islet IRE. IREs are depicted as boxes, with the filling color corresponding to the H3K27ac log₂ fold change. g, The IRE bearing the T1D associated variant rs78037977 is marked by the ENCODE ChromHMM classification as a “strong enhancer” (orange) in other non β-cell lines (left). ENCODE ChromHMM classification in non β-cell lines for the IRE bearing the T1D associated variant rs193778. h, i, Allele-specific luciferase experiments for rs78037977 (h) and rs193778 (i) in untreated EndoC-βH1. ANOVA followed by Bonferroni correction * P < 0.05; ** P < 0.01; *** P < 0.001. Bars represent mean ± sd.

Extended Data Fig. 7. Human islets and EndoC-βH1 Glucose-Stimulated Insulin Secretion (GSIS).

GSIS was assessed, in pancreatic human islets (a) and EndoC-βH1 cells (b). In the case of EndoC-βH1 cells, the experiments were performed upon exposure or not to IFNγ (1000 U/ml) + IL1β (50U/ml) for 48h. Data are mean plus range of four to eight independent experiments, and are expressed as the ratio between glucose stimulated and basal insulin secretion. *P < 0.05, **P < 0.01, ***P < 0.001, for the indicated comparisons (paired t test (a) or ANOVA followed by Bonferroni correction (b)). NT=Non treated.

Extended Data Fig. 8. ATAC-seq quality control.

a, Agilent TapeStation profiles obtained by chromatin tagmentation of human islets and EndoC-βH1 samples showing the laddering pattern of ATAC-seq libraries. The band sizes correspond to the expected nucleosomal pattern. *Notice that samples HI-19 CTRL and CYT were used as examples to illustrate the expected fragment distribution pattern in ATAC-seq experiments in Raurell-Vila et al.. b, Summary of per-replicate sequencing metrics, showing total library sizes, percentage of aligned reads, percentage of mitochondrial aligned reads, normalized strand cross-correlation coefficient (NSC, values significantly lower than 1.1 (<1.05) tend to have low signal to noise or few peaks) and relative strand cross-correlation coefficient (RSC, values significantly lower than 1 (<0.8) tend to have low signal to noise). c, TSS enrichment over a 4kb window centered on genes TSS compared to a set of genes randomized along the genome, showing the expected pattern of open chromatin centered on the TSS. d, Percentage of total reads found at called open chromatin peaks classified as distal (>2kb from TSS) or promoters (≤2kb from TSS) compared to a randomized set of peaks. e, UCSC views at islet-specific loci (NKX6.1, PDX1 and NEUROD1) showing the high reproducibility of ATAC-seq profiles among independent replicates.

Figure 1.. Pro-inflammatory cytokines exposure causes profound remodeling of the β-cells regulatory landscape.

a, Summary of the experimental design. The number of EndoC-βH1 (EC) and human pancreatic islets (HI) samples used in different assays is indicated. b, Correlation between chromatin accessibility and H3K27ac deposition, each dot corresponds to a chromatin site. The point fill refers to the ATAC-seq and the border to the H3K27ac classification (gained=green; lost=red; stable=grey). The dotted box depicts the regulatory elements referred as induced regulatory elements (IREs) and the lighter shade of green depicts a subtype named neo IREs (see text). c, Correlation between changes in RNA expression and protein abundance in EndoC-βH1 cells. Point fill and border indicate the classification of RNA-seq and protein respectively (up-regulated=green; down-regulated=red; equal-regulated=grey). d, Genes proximal to IREs (see methods) show cytokine-induced expression in EndoC-βH1 exposed or not to pro-inflammatory treatment. CYT=cytokine exposed, CTRL=control. ***Two-sided Wilcoxon test P<0.001. Box plot limits show upper and lower quartiles, whiskers extend to 1.5 times the interquartile range. e, Translation of proteins encoded by IRE-associated genes is induced by cytokine exposure in EndoC-βH1. This is shown by the significantly different (Two-sided Wilcoxon test P<2×10⁻¹⁶) log₂ FC distribution of protein abundance obtained after cytokine exposure, for proteins encoded by genes associated with IREs or stable regulatory elements (SREs). f, Representative view of the IFN-inducible guanylate binding proteins GBP4 and GBP5, illustrating their up-regulation upon cytokine exposure and the nearby induction of IREs characterized by gains in chromatin accessibility and enrichment in H3K27ac (green boxes).

Figure 2.. The β cell response to pro-inflammatory cytokine unveils neo and primed regulatory elements.

a, Classification of ATAC-seq open chromatin sites upon exposure of human β cells to IFN-γ and IL-1β. IREs=induced regulatory elements, SREs=stable regulatory elements. b, View of the SOCS1 locus, a gene strongly induced upon pro-inflammatory cytokine exposure. We here depict representative examples of primed (blue box) and neo IREs (green boxes). c, Box plot distribution of ATAC-seq and H3K27ac normalized tag counts at different classes of IREs. Box plot limits show upper and lower quartiles, whiskers extend to 1.5 times the interquartile range, individual data points represent outliers and the notch represents the confidence interval around the median. d, Islet-specific TF occupancy at neo, primed and stable regulatory elements. Read density for PDX1, NKX2.2, FOXA2, NKX6.1 and MAFB was calculated in 10bp bins in 1kb windows centered on the regulatory element. Lines represent means, while the grey shades depict the standard deviation. e, Footprint analysis of ISRE motifs in neo (left) and primed regulatory elements (right) in cells exposed or not to IFN-γ and IL-1β (control = blue; cytokines = orange). f, Violin plots showing the distribution of DNA methylation β-values in neo and primed enhancers, exposed or not to pro-inflammatory cytokines. *** Two-sided Wilcoxon test P<0.001. g, Model showing two types of IREs driving the response to pro-inflammatory cytokines in human β cells.

Figure 3.. Cytokine exposure induces changes in human islet 3D chromatin structure.

a, b, View of the UMI-4C chromatin contacts of TNFSF10 (a) and GBP1 (b) promoters, before and after exposure to pro-inflammatory cytokines. Yellow boxes indicate IREs that gain contacts with the up-regulated gene promoters. A heatmap under the 4C track represents the log₁₀ odds ratio (OR) of the UMI-4C contacts difference in cytokine vs. control. Small black diamonds on top of the contact heatmap indicate a significant difference between cytokine-treated and control samples 3D chromatin contacts (Chi-squared P<0.05). ATAC-seq peaks are represented by rectangles shaded from gray to green proportionally to the cytokine-induced H2K27ac log₂ fold change observed at that site (RE=regulatory elements track). c, Distribution of the UMI-4C contacts log₂ fold changes (cytokines vs. control) at the different types of islets open chromatin sites classified as in Fig. 2a. The data, obtained by analyzing viewpoints centered at the promoter of cytokine-induced genes, show that the chromatin structural changes are preferentially happening at IREs. Box plot limits show upper and lower quartiles, whiskers extend to 1.5 times the interquartile range, individual data points represent outliers and the notch represents the confidence interval around the median. Two-sided Wilcoxon test * P<0.05, ** P<0.01, *** P<0.001.

Figure 4.. Cytokine-induced islet regulatory elements map to T1D associated regions and guide the identification of functional risk variants.

a, Islet IREs are enriched for T1D but not T2D risk variants while the opposite is true for islet SREs. Significant Varian Set Enrichment (VSE) scores are depicted in red (Bonferroni adjusted P<0.05). Box plot limits show upper and lower quartiles, whiskers extend to 1.5 times the interquartile range and notch represents the median confidence interval for distributions of matched null sets (500 permutations). b, rs78037977 overlaps an IRE bound by islet-specific TFs under basal conditions. c, Luciferase assays in EndoC-βH1 exposed to cytokines show that the sequence exerts enhancer activity which is reduced in the T1D associated allele (G). d, UMI-4C in human islets show that the IRE containing rs78037977 engages multiple distal chromatin contacts. e, Zoom-in at one induced chromatin contact, mapping to the up-regulated TNFSF18 gene. f, Variant rs193778 maps to a phylogenetically conserved IRE. g, Luciferase assays in EndoC-βH1 exposed to cytokines show significantly increased enhancer activity of the risk (G) allele compared to the non-risk (A) allele. h, UMI-4C in islets using the promoter of DEXI as viewpoint, show a chromatin contact with the IRE bearing the T1D risk variant. ATAC-seq peaks are represented, in d and h, by rectangles shaded proportionally to the H2K27ac log₂ fold change (RE=regulatory elements). Statistical significance in c and g was determined by one-way ANOVA tests followed by Bonferroni’s correction. *P<0.05, **P<0.01, ***P<0.001. Bars represent mean ± sd.