Heterogeneous enhancer states orchestrate β cell responses to metabolic stress | Nature Communications

Nature Communications volume 15, Article number: 9361 (2024) Cite this article

Metrics details

Obesity-induced β cell dysfunction contributes to the onset of type 2 diabetes. Nevertheless, elucidating epigenetic mechanisms underlying islet dysfunction at single cell level remains challenging. Here we profile single-nuclei RNA along with enhancer marks H3K4me1 or H3K27ac in islets from lean or obese mice. Our study identifies distinct gene signatures and enhancer states correlating with β cell dysfunction trajectory. Intriguingly, while many metabolic stress-induced genes exhibit concordant changes in both H3K4me1 and H3K27ac at their enhancers, expression changes of specific subsets are solely attributable to either H3K4me1 or H3K27ac dynamics. Remarkably, a subset of H3K4me1+H3K27ac- primed enhancers prevalent in lean β cells and occupied by FoxA2 are largely absent after metabolic stress. Lastly, cell-cell communication analysis identified the nerve growth factor (NGF) as protective paracrine signaling for β cells through repressing ER stress. In summary, our findings define the heterogeneous enhancer responses to metabolic challenges in individual β cells.

The rise of the obesity pandemic and the concurrent global increase in the incidence of type 2 diabetes (T2D) present a significant challenge to public health. Dysfunction of β cells serves as an early hallmark of prediabetes and is essential for the progression of T2D1. Overnutrition leads to metabolic and inflammatory stress in β cells, thus impairing their insulin secretion capacity2,3. Extensive research has shed light on the cellular and molecular mechanisms underlying β cell dysfunction. Significant differences in gene expression, glucose-stimulated insulin secretion, and stress responses exist between normal and dysfunctional β cells4,5,6,7,8. Moreover, in both normal and diabetic islets, considerable heterogeneity at transcription, protein, and functional levels has been observed across all major endocrine cell types9,10,11,12. Recent advances in single-cell RNA sequencing (scRNA-seq), single-cell Assay for Transposase-Accessible Chromatin sequencing (scATAC-seq), and single-cell multiome have established transcriptional and epigenetic signatures that define the heterogeneous states of β-cell13,14,15,16,17,18,19,20,21,22,23. In addition, by integrating chromatin accessibility, transcription factor motif enrichment, and gene expression at the single-cell level, these studies were able to infer the β cell-specific gene regulatory network (GRN)24 and link T2D GWAS variants to β cell function13.

Gene expression is governed by transcription factors (TFs) and epigenetic cofactors that bind to both proximal and distal enhancers, collectively referred to as cis-regulatory regions (CREs)24. The heterogeneity in gene expression and CRE activity has been instrumental in defining subpopulations in islet endocrine cells, particularly in β cells13,14. In mature β cells, the heterogeneous gene expression corresponds to functional heterogeneity at single-cell level22,25,26,27. While sometimes described in simplistic binary terms as open or closed, CREs are, in fact, demarcated by a spectrum of histone modifications28,29. The chromatin states defined by the histone code are linked to diverse biological functions, providing a nuanced perspective on gene regulation30. Functional studies in islets have revealed that the perturbation of epigenetic modulators for different enhancer marks showed distinct phenotypes, suggesting both collaborative and divergent regulatory networks31,32,33,34. Nevertheless, the precise way in which specific enhancer states contribute to gene expression dynamics in islet cells at the single-cell level remains to be fully elucidated. Furthermore, the role of enhancer activity in generating heterogeneous gene expression, especially in the context of β cell dysfunction, is not well defined.

The recent development of multi-omic platforms such as Paired-Tag35 and CoTECH36, based on the combinatorial indexing platform, enables joint profiling of histone modifications and RNA from single nuclei at an unprecedented scale. Here, we use Paired-Tag to investigate the multifaceted cellular trajectories defined by the integration of transcriptome and epigenome of H3K27ac or H3K4me1, in islets of mice under normal chow (NC) or high-fat diet (HFD). Our results demonstrated consistent trajectories of islet dysfunction defined by both histone marks. The trajectories from healthy to dysfunctional β cells are defined by elevated metabolic stress, as well as reduced β cell identity and function at both transcription and enhancer level. Interestingly, the degree of coupling between H3K27ac/H3K4me1 and gene expression is highly locus-specific. The complex relationship can display either concordant or discordant patterns between H3K27ac and H3K4me1 in different gene subsets. Through gene regulatory network (GRN) analysis and pharmacological perturbation, we revealed distinct regulatory modes for inflammation- and ER stress-induced genes. Surprisingly, we found that a large portion of H3K4me1+/H3K27ac- primed enhancers in NC islets were lost after HFD. These enhancers are occupied by the pioneer factor FoxA2. The loss of primed enhancers is linked to the suppression of β cell function genes by HFD. Lastly, we identified intra-islet cell-cell communication that regulates HFD-induced transcriptome changes in β cells and identified NGF as a protective paracrine factor for β cells. Importantly, the administration of NGF curbs excessive ER stress-induced β cell dysfunction and ameliorates hyperglycemia in diabetic Akita mice. Taken together, our data presented a dynamic and comprehensive regulatory roadmap of β cell dysfunction under metabolic stress, revealed the essential and often discordant H3K27ac/H3K4me1 enhancers activities, and identified NGF-TrkA (Tropomyosin receptor kinase A) signaling that can be modulated to ameliorate β cell dysfunction.

To identify the gene expression dynamics and their relationship with enhancer states in islet dysfunction, we employed Paired-Tag35 to simultaneously quantify the enhancer marks (H3K27ac or H3K4me1), and RNA transcription from single nuclei. It is well documented that HFD induces islet dysfunction due to both glucolipotoxicity and insulin resistance in mice2. We harvested islets from male and female C57BL/6J mice under normal chow (NC) or high-fat diet (HFD) (8 weeks or 16 weeks), which caused increased body weight, mild hyperglycemia (Supplementary Fig. 1), and β cell proliferation and stress37,38,39. Using a “split-and-pool” strategy to index both RNA and enhancers from the same nuclei35, this protocol enables us to establish multi-omic datasets of “RNA+H3K4Me1” and “RNA+H3K27ac” (Fig. 1A).

A Schematic experimental design. Created in BioRender. Wei, Z. (2022) BioRender.com/n56g182. B Uniform manifold approximation and projection (UMAP) embedding showing the clustering of 56,748 single nuclei by single nuclei RNA profile. C Row-normalized expression of top enriched genes for each islet cell type (MPs: macrophages, FBs: Fibroblasts, ECs: Endothelial cells). Color bars in ID represent individual biological samples. D RNA velocity flow projected in the UMAP space. E RNA velocity-based latent time of the nuclei shown in D. F Smoothed gene expression trends of the top 1,000 genes whose expression values correlate best with β cell state probabilities (from beta-hi to beta-low), sorted according to peak in pseudotime. G The RNA velocity (upper panel) and the expression (lower panel) for representative genes (Iapp, Slc2a2, Hspa5, Ddit3) in each cluster by scVelo. Positive velocity: up-regulated; negative velocity: down-regulated. H UMAP of 3 states of mouse β cells (high, low, and proliferating (Mki67+)) density from NC or HFD mice, showing the relative proportion of beta-low cells are increased in HFD islets. I The percentage of cells in each state (high, low, and Mki67+) in NC and HFD β cells. J Venn diagram of HFD-repressed (upper panel) and HFD-induced (lower panel) genes in high, low, and Mki67+ β cells. K Heatmap showing the odds ratio of the top GO terms calculated by enrichR124 for up and downregulated DEGs within each cell type (*adjusted p < 0.05). Top panel: GO terms for genes repressed by HFD. Lower: GO terms for genes upregulated by HFD. Adjusted p-values were calculated using Fisher-exact test and Benjamini-Hochberg method for correction.

After quality control (see methods and Supplementary Fig. 2) for both RNA and DNA reads, we first integrated the single nuclei RNA data across all time points and both sexes. We identified all major endocrine and nonendocrine islet cell types (Fig. 1B), consistent with multiple published scRNA-seq and snRNA-seq13,14,15,16,17,18,19,20,21,22,23. We were able to clearly define various known cell populations (α, β, δ, PP, ductal, endothelial, fibroblast, and macrophage) based on cell type-specific expression signatures (Fig. 1B, C, Supplementary Fig. 3A). Next, we used scVelo40 and CellRank41 to unbiasedly determine the cell kinetics from healthy (NC) state to dysfunctional (HFD) state in the integrated β cell population. Interestingly, we found clear directions of RNA velocities pointing to a subset of β cells (Fig. 1D), which is consistent with the pattern of latent time (Fig. 1E). Multiple scRNA-seq studies in mouse and human islets have identified subpopulations of β cells13,14,17,42,43. Based on clusters defined by Seurat44 (Supplementary Fig. 3B) and driver genes resolved along latent time (Fig. 1F), we defined two major β cell subpopulations in our datasets (β-hi, beta-low), and an additional minor population of proliferating β cells (Mki67-β) (Fig. 1D). The top driver genes (Fig. 1F, G) enriched in beta-hi states include hormonal genes (Ins2, Iapp, Chga), β cell-essential transcription factors (Rfx6, Foxo1), and genes essential for nutrient metabolism and insulin secretion (Rptor, Syt7, Pak3, Dcx). Conversely, the beta-low state is characterized by genes related to the unfolded protein response (UPR) (Pdia3, Ddit3, Atf3, Herpud1, Hspa5, Eif2ak3), ER-associated degradation (Sel1l, Herpud1) and inflammatory stress (Nfkbiz) (Fig. 1F). These results confirmed that the cellular trajectories defined here reflect a progression from healthy β cells to stressed and dysfunctional cells. Though we defined the two states based on integrated RNA, the beta-low and beta-hi states do not simply reflect diet conditions (NC or HFD). In fact, both NC and HFD β cells contain substantial beta-low and beta-hi state cells (Fig. 1H). However, the ratio of beta-hi to beta-low changed dramatically in HFD islets (Fig. 1H, I), suggesting that HFD induces a progression from beta-hi state to beta-low state. RNA velocity results suggest that beta-low cells have the most evident dynamics (Fig. 1D). Indeed, when we quantified HFD-induced gene expression changes in β-hi, beta-low, and Mki67+ β cells, the number of differentially expressed genes was much larger in beta-low cells (Fig. 1J, Supplementary Data 1). Gene ontology analysis further confirms that HFD-repressed genes are related to insulin secretion, while HFD-induced genes are enriched in excessive ER stress-induced unfolded protein response (UPR) (Fig. 1K, Supplementary Fig. 3C). Lastly, we found a substantial overlap of dynamic genes defined here with the published GWAS-associated T2D genes or previous single-cell studies in human T2D islets16,45 (Supplementary Fig. 4), suggesting that the dysfunction trajectories of mouse and human β cells may be evolutionarily conserved. Together, these results demonstrate that at the transcriptome level: 1) β cells are heterogeneous, 2) HFD induces a transition from beta-hi state to beta-low state, and 3) the main signature of this dysfunction trajectory is the downregulation of β cell function genes and the elevation of ER stress and proinflammatory genes induced by HFD.

Our analysis also revealed a minor β cells population (Mki67 + -β) that highly expresses proliferation markers such as Mki67, Top2a, Cenpp, and Smc4. (Fig. 1C, Supplementary Fig. 5A, B). These cells showed substantial differences in gene expression when compared with either beta-hi or beta-low cells (Supplementary Fig. 5C, D). Gene ontology enrichment analysis revealed that in contrast to beta-hi or beta-low cells, the Mki67+ β cell-enriched genes are uniquely associated with cell cycle, mitotic chromosome condensation, and separation (Supplementary Fig. 5E). Motif analysis showed that the H3K27ac and H3K4me1 enhancers in Mki67 + -β cells are specifically enriched in multiple TF binding motifs that are known to be related to β cell proliferation, such as FoxO146, PITX247, STAT548, and VDR49 (Supplementary Fig. 5F).

Lastly, we examine sexual dimorphism in gene expression in mouse islets. Sex differences have profound effects on islet function and gene expression50. Our analysis revealed that cells do not separate into sex-specific clusters (Supplementary Fig. 6A), indicating that overall gene signatures for each cell type are similar in both sexes. When comparing the male and female β cells at 0, 8, and 16 weeks, we found a total of 84 genes that are differentially expressed in male or female β cells (adj. p < 0.05, fold>1.25). As expected, this list includes X-linked (Tsix, Xist, Gm26992) or Y-linked (Uty, Eif2s3y) genes (Supplementary Fig. 6B–E). For the other autosomal genes, Ins1 and Ins2 are expressed at a higher level in male β cells, consistent with previous scRNA-seq results51.

To further investigate enhancer dynamics in islets, we integrated multimodal single-cell datasets of H3K27ac+RNA or H3K4me1+RNA using Seurat (Fig. 2A, B). The integrated datasets successfully classify similar cell types as defined solely by RNA (Fig. 2A, B), and the majority of cells retained their identities defined by RNA (Supplementary Fig. 7A, B). To explore the link between gene expression and cell type-specific enhancers, we first pooled H3K27ac or H3K4me1 reads by individual endocrine cell types. As anticipated, both H3K27ac and H3K4me1 displayed highly specific enrichment at loci of cell type-specific genes, such as Gcg, Ins1, Sst, Mki67, and Ppy (Fig. 2C). To further quantify the coupling between H3K27ac/H3K4me1 and target gene expression, we used the FigR package52 to compute the DORC (domain of regulatory chromatin) value, a measurement of accumulated H3K27ac or H3K4me1 signals corresponding to each gene locus. Globally, the DORCs of both H3K27ac and H3K4me1 clustered coherently with gene expression (Supplementary Fig. 7C, D). Moreover, for the top 50 cell type-specific genes defined by RNA, both H3K27ac and H3K4me1 DORCs of these loci are also highly cell-type specific (Fig. 2E, F, Supplementary Data 2 & 3). Beyond identifying cell-type specific enhancers, we were also able to quantify the strength of correlations between individual H3K4me1+ or H3K27ac+ enhancers and target expression at both bulk and single-cell level (Fig. 2G, H). Our previous analysis pinpointed a set of heterogeneous genes that align with the latent time of β cell dysfunction trajectory (Fig. 1F). Nevertheless, it is not clear whether these heterogeneous gene expression patterns correlate with local enrichment of H3K27ac or H3K4me1. To identify genes displaying concordant dynamics of RNA and H3K27ac/H3K4me1 levels, we generated pseudotime profiles of H3K27ac+RNA or H3K4me1+RNA multiome from healthy (beta-hi) to dysfunctional states (beta-low) using ArchR53 (Fig. 2I, J). We subsequently ranked individual cells based on dysfunction pseudotime and identified genes of which the dynamics of RNA and H3K27ac/H3K4me1 are coupled (Fig. 2K, L). Among the genes identified, those elevated during β cell dysfunction are notably associated with excessive ER stress and inflammatory responses (Fosb, Pdia3, Ddit3, Hspa5, Hsp90b1, Manf, etc.), while the β cell function genes (Slc2a2, Tcf7l2, Vdr, Abcc8, Glp1r, Neurod1, etc.) are repressed (Fig. 2K, L, Supplementary Fig. 8). Taken together, these findings identified genes of which the H3K27ac and/or H3K4me1 enhancer activities are coupled with gene expression.

UMAP embedding from single-nucleus joint profiles of transcriptome and histone modification (H3K27ac (A) or H3K4me1 (B)). Total number of nuclei: 37,847 in A and 18,919 in B. Genome tracks of aggregated normalized signals of H3K27ac (C) and H3K4me1 (D) of individual endocrine cell types (α, β, δ, Mki67-beta, and PP cells) at Gcg, Ins1, Sst, Mki67, and Ppy loci.Heatmaps of RNA expression (left) and the accumulated H3K27ac (E) or H3K4me1 (F) signals (right). The top 50 differentially enriched genes for each cell type were selected. Correlations between individual H3K4me1 (G) or H3K27ac (H) enhancer activity and the expression of Glrx3, a gene involved in glutathione metabolism and a marker for stressed β cells125. Top 3 panels: genome tracks of pooled H3K4me1 (G) or H3K27ac (H) signals in α, β, and δ cells. Mid panel: H3K4me1 (G) or H3K27ac (H) peaks. Bottom panel: spider plot represents the strength of linkage between H3K4me1 (G) or H3K27ac (H) and Glrx3 expression. (I-J) UMAP shows the trajectory of β cell dysfunction from β-high to beta-low states in H3K27ac(I) and H3K4me1(J). Pseudo-time values were overlaid on the UMAP embedding; the smoothed line and arrow represent the visualization of the trajectory path from the spline fit. K Heatmaps showing the enrichment patterns of H3K27ac (left) and expression (right) of 1,657 genes of which the H3K27ac and RNA are correlated along the pseudo-time order. L Heatmaps showing the enrichment patterns of H3K4me1 (left) and expression (right) of 964 genes of which the H3K4me1 and RNA are correlated along the pseudo-time order.

To infer GRN based on H3K27ac and H3K4me1, we first employed ArchR53 to compute the enrichment of TF motifs in individual β cells. Consistent with the heterogeneous transcriptome and H3K27ac/H3K4me1 in single cells, the motif activities for individual TFs involved in β cell function (FoxA2) or stress response (FOS, IRF1, NFKB2) are highly heterogeneous with relative enrichment in beta-hi or beta-low subpopulations (Fig. 3A). To discern putative TFs that may regulate β cell dysfunction, we correlated TF expression with motif activity along the trajectories of β cell dysfunction (Fig. 3B, C, Supplementary Data 4). This approach identified several known TFs associated with identity and function of β cells (FOXA2, FOXO1, HNF1β, NR4A1), nuclear hormone signaling (RARα/β, RXRα/γ, ERα, PPARδ, RORα), and UPR/inflammatory stress (FOS/FOSB, IRF1/7, ATF3/6, RELA/B) (Fig. 3B, C, Supplementary Data 4). To distinguish activators from repressors, we utilized FigR52 to rank the mean regulation scores of all TFs, revealing those that predominantly act as activators (IRF1, FOXA3, RARB) versus those function mainly as repressors (JUN, NR3C1, TBX1) at H3K27ac or H3K4me1 enhancers (Fig. 3D, E). Next, we aggregated all GRNs based on the clustering of TFs and target genes. The analysis uncovers that inflammation modulators (e.g., RELA, RELB, IRF1, NFKB1, STAT1) tend to coactivate downstream proinflammatory targets (Fig. 3F, G). Conversely, TFs known to repress inflammatory responses, such as glucocorticoid receptor (NR3C1) and vitamin D receptor (VDR), cluster together and counteract stress-induced targets (Fig. 3F, G). The resemblance between H3K27ac or H3K4me1 GRNs highlights that the transcriptional responses to inflammation in HFD islets are coupled with concordant H3K27ac and H3K4me1 dynamics. Furthermore, GRNs reveals putative H3K27ac- or H3K4me1-associated TF regulators for individual targets (Supplementary Fig. 9). For instance, NR4A1, an orphan nuclear receptor that regulates β cell proliferation54, emerged as a prominent activator for Itgb1, a gene required for β cell expansion and maturation55,56,57(Supplementary Fig. 9A, B). Meanwhile, other activators of Itgb1 were specific to either H3K27ac or H3K4me1 GRNs (Supplementary Fig. 9A, B). Taken together, the inferred GRNs from H3K27ac or H3K4me1 demonstrate that many pro- and anti-inflammatory TFs synergize to co-regulate downstream transcriptional responses in metabolically stressed β cells.

A Feature plot visualization of Z-score normalized TF motif activities of H3K27ac (FOS, FoxA2, upper 2 panels), and H3K4me1 (IRF1, NFKB2, lower two panels) showing the single cell heterogeneity of TF activities in β cells. B Heatmaps showing 76 TFs with a highly correlating TF expression (left) and H3K27ac motif activities (right) over pseudo-time. C Heatmaps showing 113 TFs with a highly correlating TF expression (left) and H3K4me1 motif activities (right) over pseudo-time. Mean regulation score (signed, −log10 scale) shows TF activators(left) and TF repressors(right) based on DORCs of H3K27ac (D) or H3K4me1 (E). Heatmap of DORC regulation scores for all significant TF-targets pairs in H3K27ac (F) and H3K4me1 (H–J) qPCR of expression of selected IL1β-induced NF-κB targets in MIN6 cells treated with or without SAHA (5 uM) (H), GNE 781 (1 uM) (I), or GSK-LSD1 (5 uM) (J). Gene expression was normalized to the expression of Gapdh (For Tgfbr1 (H), Ntrk2 (J), n = 3 biological independent experiments; for rest genes, n = 4 biological independent experiments). Data are shown as the mean + S.E.M. Analysis was done using two-tailed Student’s t-test. *p < 0.05, **p < 0.01, and ***p < 0.001. Source data and exact p values are provided as a Source Data file.

The GRNs of NF-κB pathway TFs (NFKB1/2, RELA/B) suggest a robust correlation between H3K27ac/H3K4me1 enhancers and inflammation-driven dysfunction of β cells. However, whether the H3K27ac/H3K4me1 enhancer activities are required for individual NF-κB downstream gene expression in β cells remains undetermined. Based on the in silico H3K27ac/H3K4me1 GRNs (Fig. 3F-G), we selected a number of NFKB1/RELA target genes of which the expression is significantly linked to H3K27ac and/or H3K4me1, and tested the causality between H3K27ac/H3K4me1 enhancer activities and gene expression by perturbing the epigenetic modulator of H3K27ac/H3K4me1 (CBP, HDACs, or LSD1) with HDACs inhibitor SAHA58, CBP inhibitor GNE78159, and LSD1 inhibitor GSK-LSD160. In MIN-6, a mouse β cell line, the global H3K27ac is significantly increased upon SAHA treatment, and decreased by GNE-781 (Supplementary Fig. 9E, left). Global H3K4me1 level is increased by GSK-LSD1 (Supplementary Fig. 9E, right). Next, we applied IL-1β to activate NF-κB signaling, and found that many of our selected genes, such as Itgb1, Tgfbr1, Jak2, Il1r2, and Ezh2, are significantly elevated when treated with SAHA (Fig. 3H). Conversely, when using GNE781 to inhibit CBP and specifically reduce H3K27ac, we found that the expression of NF-κB targets such as Traf3ip2, Nr1h4, and Hdac9, are significantly reduced (Fig. 3I). Lastly, we boosted H3K4me1 by inhibiting H3K4me1 demethylase LSD1 and found that the expression of NF-κB targets (Ntrk2, Nr1h4, Itgb1, and Trafip2) are elevated (Fig. 3J). Collectively, results from these pharmacological perturbations provide strong evidence that the NF-κB mediated inflammatory responses within β cells depend on the dynamics of H3K27ac and H3K4me1 enhancers.

In a simplified model of gene regulation, the activation of a gene from a silent state is concomitantly “coupled” with the increase in both H3K27ac and H3K4me1 enhancer activities. However, gene expression heterogeneity in islets is mostly quantitative rather than binary. Whether the heterogeneity of gene expression is associated with simultaneous changes in H3K27ac and H3K4me1 at a single-cell level requires further investigation. We first plotted the correlation coefficient of RNA-H3K27ac and RNA-H3K4me1 for each individual gene (Fig. 4A, Supplementary Data 5). The H3K27ac-RNA coefficient correlates with the H3K4me1-RNA coefficient (r = 0.51) (Fig. 4A, Supplementary Data 5). This is best demonstrated by a group of 182 “concordant genes” (coefficient r > 0.7 for both H3K27ac and H3K4me1, green dots in Fig. 4A), of which both H3K27ac and H3K4me1 are highly correlated with gene expression dynamics. For the majority of genes (Ddit3, Manf, Fos, Atf6, and G6pc2) in this category (Fig. 4A, Supplementary Data 5), the expression and H3K27ac/H3K4me1 levels are elevated in dysfunctional β cells (Fig. 4B, C). In contrast, we also found a substantial number of outliers that are “discordant” between H3K27ac and H3K4me1 dynamics (group 2 and 3 in Fig. 4A). For group 2 genes (95 in total), transcription is coupled with H3K27ac but not H3K4me1 (Fig. 4A, red dots). For group 3 genes (607 in total), transcription is coupled with H3K4me1 but not H3K27ac (Fig. 4A, yellow dots). Interestingly, the 3 groups of genes are associated with distinct gene ontology (GO). Group 1 genes are highly induced in a stressed state (Fig. 4B, C) and are mostly related to UPR/ER stress (Fig. 4D). In contrast, group 2 genes are enriched in the calcium signaling pathway (Fig. 4D), whereas group 3 genes are enriched in secretion, inflammatory response, and to a less extent, ER stress pathways (Fig. 4D). For individual discordant genes, this divergence can simply be demonstrated by comparing the trajectories of H3K27ac and H3K4me1 side-by-side (Fig. 4E). For group 2 genes such as Sik3, Ate1, and Sec24d, H3K4me1 levels decrease along with the dysfunction trajectory, which correlates with the RNA dynamics. The H3K27ac level, on the other hand, remains largely constant and therefore are decoupled from the RNA changes (Fig. 4E). For group 3 genes such as Dusp10, Lrp11, or Kcmb2, it is the H3K27ac, rather than H3K4me1, that is coupled with gene expression (Fig. 4E). As the majority of H3K27ac and H3K4me1 peaks are located in the proximal/distal enhancer regions, the contribution of H3K27ac at promoters to this discordance is likely limited (Supplementary Fig. 10). The concordance and discordance of two enhancer marks with RNA dynamics may be attributed to the dynamic equilibrium between histone acetyltransferases (HATs) to histone deacetylases (HDACs) for H3K27ac, and histone methyltransferases (HMTs) to lysine demethylases (KDMs) for H3K4me1, and is likely governed by the preference of epigenetic cofactors by local TFs. Indeed, while some TF modulators are enriched in both H3K27ac and H3K4me1 enhancers (ATF3, ATF6, ZFP148), a large number of TFs are exclusively associated with either H3K27ac (YY1, FOS, FOSB, NR4A1, RREB1, HNF4α) or H3K4me1 (FoxA2, FoxO1, HNF1β, MAFA, PDX1) (Fig. 4F). Next, we performed RNA sequencing (RNA-seq) on islets from 16 week HFD mice treated with vehicle, HDAC inhibitor (SAHA), which increases H3K27ac61,62, or CBP/p300 inhibitor GNE-781 and C646, which mildly decrease the global H3K27ac63,64, to determine transcriptional changes. In MIN6 cells, the global H3K27ac is significantly increased upon SAHA treatment, and decreased by GNE-781 and C646 (Supplementary Fig. 9E, left). In primary islets, exposure to SAHA induced 2,255 genes, while GNE-781 or C646 repress 3,262 genes and 1,303 genes, respectively. When we compared the perturbed genes with 855 genes of which RNA dynamics is coupled with H3K27ac (group 1 and group 3 in Fig. 4A), we found that 94 H3K27ac-RNA coupled genes are upregulated by SAHA, whereas 320 H3K27ac-RNA coupled genes are downregulated by GNE-781 or C646 (Fig. 4G). Gene ontology of these perturbed genes (412 genes in total) shows the association with response to ER stress and insulin secretion (Fig. 4H). The substantial overlap indicates that for a majority of genes in group 1 and 3, H3K27ac directly contributes to the heterogeneous gene expression. Next, we evaluated the H3K27ac perturbation on β cell heterogeneity by using a well-characterized cell surface gene, CD63, as a surrogate marker for “healthy” β cells15. In our dataset, CD63 is expressed at a higher level in beta-hi cells compared to beta-low cells (Supplementary Fig. 11A). Interestingly, treatment of GNE-781, C646, or SAHA all reduced the CD63hi β cells% compared to vehicle (Supplementary Fig. 11B), suggesting that any perturbation of the “optimal” H3K27ac level may shifted the heterogeneity of β cells to a more dysfunctional state. Lastly, we also found a significant overlap between genes in Fig. 4A and H3K27ac/H3K4me1-associated genes in human islets/progenitors defined by bulk ChIP-seq in previous studies65,66 (Supplementary Fig. 12), suggesting that the H3K4me1/H3K27ac-to-RNA linkage may be conserved in human β cells. Taken together, these data suggest that the concordance between H3K27ac and H3K4me1 is highly variable across individual genes. The coupling (or decoupling) of H3K27ac/H3K4me1 to transcription appears to be determined by the preferential (or lack of) engagement of HATs/HDACs or HMTs/KDMs by local TFs.

A Scatter plot showing the correlation coefficients of H3K27ac-RNA (y-axis) or H3K4me1-RNA(x-axis) for all genes. Group 1 (green): Genes showing a concordant pattern of H3K27ac-RNA and H3K4me1-RNA correlations. Group 2 (yellow): Genes showing only significant correlations between H3K4me1 and RNA. Group 3 (red): Genes showing only significant correlations between H3K27ac and RNA. Group 4 (grey): Other genes that do not fit into the above correlation patterns. B Pseudo time heatmap of the H3K4me1 gene score (left) and gene expression (right) for all group 1 genes. C Pseudo time heatmap of the H3K27ac gene score (left) and gene expression (right) for all group 1 genes. D Gene ontology for group 1-3 genes showing enrichment in distinct enriched categories. Adjusted p-values were calculated using Benjamini-Hochberg correction. E H3K4me1 (left) and H3K27ac (right) gene scores of representative genes in group 2 (Sik3, Ate1, Sec24d) and group 3 (Dusp10, Kcnmb2) showing the discordant patterns of H3K4me1 and H3K27ac over the β cell dysfunction trajectory. Dots represent gene scores of an individual pseudotime-ordered histone modifications profile. The smoothed line and arrow represent the visualization of the trajectory path from the spline fit. F Heatmap showing the correlation between TF expression and motif activity (left: H3K4me1, right: H3K27ac). G Pie chart showing the H3K27ac-coupled genes (group 1 and 3 from A) that are significantly suppressed by CBP inhibitor GNE-781 or P300 inhibitor C646 (red), or induced by HDACs inhibitor SAHA (green), in primary mouse islets treated with SAHA (5 uM), GNE-781 (1 uM) or C646 (5 uM) for 24 h. H Gene ontology categories of 414 overlapped genes in G.

When generating the multiome UMAPs of H3K27ac and H3K4me1 using the same pipeline, we noticed that β cells in H3K4me1-RNA UMAPs are more diffused compared to H3K27ac (Fig. 2A, B). To investigate the underlying cause of this discrepancy we labeled all β cells by time points (NC, HFD 8 weeks, or HFD 16 weeks) in H3K4me1-RNA UMAP. This revealed an obvious separation between NC and HFD (8 or 16 weeks) (Fig. 5A). To further investigate the HFD-induced differences in H3K4me1, we first pooled H3K4me1 and H3K27ac signals by each point. We found that many H3K4me1 enhancers in NC are lost after HFD and these enhancers are generally devoid of H3K27ac signals (Fig. 5B). In development, the H3K4me1+/H3K27ac- enhancers are usually classified as poised or primed enhancers67,68 (Fig. 5C). These enhancers are usually occupied by pioneer factors such as FoxA1/269 in progenitors, which open up the chromatin and eventually lead to activation in mature β cells70. However, the prevalence and function of primed enhancers in terminal differentiated mature β cells are less clear. We defined primed enhancers in NC islets and found that about 14% of all H3K4me1 peaks are H3K4me1 + /H3K27ac-. The pattern of H3K4me1 and H3K27ac in these primed enhancers is consistent with published ChIP-seq datasets (Supplementary Fig. 13). Interestingly, the majority of these enhancers lose H3K4me1 upon HFD (Fig. 5E), suggesting an overhaul of primed enhancers by HFD. The most enriched motifs in the primed enhancers are for Fox family TFs (FoxA1/2/3, FoxM1) (Fig. 5F), consistent with their known role as pioneer factors. In postnatal β cells, FoxA2 is essential for regulating key β cell identity genes71,72,73. The expression pattern of FoxA2 is associated with its motif enrichment in H3K4me1+ enhancers, but not in H3K27ac+ enhancers (Fig. 4F). Indeed, by comparing with previously published FoxA2 cistrome in mouse islets71, we found that the majority of the 14,116 NC primed enhancers are occupied by FoxA2 (Fig. 5G). Primed enhancers in development are usually associated with inactive genes. However, many primed enhancers defined in this study reside within gene loci which also harbors classic active enhancers (H3K4me1+/H3K27ac+) in cis (Fig. 5B). In this case, the question is whether the loss of primed enhancers could negatively influence transcription. To investigate this possibility, we first identified primed peaks that are coupled with RNA and found that ~14% (2,044) of all primed enhancers are in this class (Fig. 5H), of which 44% are occupied by FoxA2 (Fig. 5H). GO analysis of the putative targets of the 2,044 enhancers showed enrichment of TCA cycle and peptide hormone stimulus (Fig. 5I), including genes that are known to be linked to β cell function in diabetes such as Hnf1a74,75 and Arap176,77 (Fig. 5J, K). Taken together, these data revealed a surprising finding that the overnutrition stress caused a genome-wide loss of H3K4me1+/H3K27ac-, FoxA2 occupied, primed enhancers. The reprogramming of these primed enhancers contributes to HFD-repressed islet function genes.

A The integrated H3K4me1-RNA UMAP of β cells at NC (0 weeks) and HFD (8 and 16 weeks) shows a distinct pattern of NC β cells. B The pseudo-bulked H3K4me1(red) and H3K27ac(blue) at Fam210b, Ptpn1, and Slco3a1 show the complete loss of H3K4me1+H3K27ac- primed enhancers (grey) in HFD islets. C Schematic demonstration of the conventional model of primed (H3K4me1+H3K27ac-) and active (H3K4me1+H3K27ac+) enhancers. D Pie chart of the percentage of primed enhancers in all H3K4me1 enhancers. E Heatmap of all K4 + K27- primed enhancers in NC and HFD (8 and 16 weeks), showing a majority of the primed enhancers in NC are lost in HFD. The proportion of primer peaks relative to the total number of H3K4me1 peaks. F Motif enrichment of the primed enhancer in E. P-value were calculated using HOMER’s default method. G Heatmaps of pseudo-bulked H3K4me1, H3K27ac, and mouse islets FoxA2 ChIP-seq (from GSM 1306337) show the extensive occupancy of FoxA2 at primed enhancers. H Pie chart of the percentage of primed enhancers (2,044 peaks, 14%) that correlate with gene expression (left) and the percentage (44%) of 2,044 peaks that are occupied by FoxA2. I Top enriched GO terms for genes associated with the 2,044 peaks in H. Linkage plot for Hnf1a (J) and Arap1 (K) loci. Top left: pseudo-bulked H3K4me1 in NC and HFD β cells, grey regions indicate primed enhancers associated with Hnf1a gene expression. Top right: violin plot of Hnf1a expression in NC and HFD β cells. Middle: FoxA2 ChIP-seq, FoxA2 motif, and H3K4me1 peaks. Bottom: spider plots showing significant enhancer(H3K4me1)-RNA links.

The extensive transcriptome and epigenome differences between NC and HFD β cells are induced by dramatic changes in nutrient availability and the microenvironment of the islets. While the profound impacts on islet function by hyperglycemia, hyperlipidemia, and hyperinsulinemia in T2D have been extensively studied, many intra-islet paracrine pathways remain less explored. We utilized NicheNet78 to infer the ligand-receptor relationship and specifically look for ligands from non-β cells that account for the changes of target gene expression between beta-hi and beta-low cells. We found multiple ligands that are known to be essential for β cell functions, such as GCG, IAPP, and FGF1 (Fig. 6A, Supplementary Fig. 14). We decided to focus on the less characterized ligands that regulate ER stress targets, as elevated ER stress is a major phenotype in beta-low cells. MANF, a secretion factor known to repress ER stress in β cells79, showed up in our list (Fig. 6A). Interestingly, we also found that NGF has a high predicted ligand-target regulatory potential (Fig. 6A). NGF is known for its variable functions in different tissues, including regulating neuron survival and growth80. Previous studies showed that NGF induced F-actin reorganization to regulate insulin secretion in β cells81,82,83. However, its potential role in antagonizing β cell ER stress and dysfunction has not been characterized. Expression pattern of Ngf suggests that it is mostly produced in ductal and fibroblasts (Fig. 6B, left), which is in contrast with the widely expressed Manf (Fig. 6B, right). This is further confirmed by an integrated islet scRNA-seq atlas57 (Fig. 6C). Although the NicheNet does not differentiate the direction of gene changes, previous studies in neuron cell lines such as PC12 suggest potential repression of ER stress by NGF84. NGF receptor TrkA is expressed in mouse islets81. Published RNA-seq of purified α, β, and δ cells from DiGruccio MR et al.85 also showed that Ntrk1 (TrkA) is enriched in β cells but not α or δ cells (Supplementary Fig. 15). NGF can be detected in lysate from pancreas, and the protein concentration is reduced in HFD-fed mice (Fig. 6D), which is consistent with scRNA-seq results (Fig. 6C). Exposure to NGF (100 ng/ml81,86) led to decreased expression of ER-stress related genes (Ddit3, Ppp1r15a, and Atf3) in islets from HFD mice (Fig. 6E). To further confirm that NGF directly repress excessive ER stress in β cells, we treated MIN6 cells with thapsigargin (TG) to induce dramatic transient ER stress. Consistent with previous reports, MANF is able to significantly reduce several UPR genes such as Atf3 (Fig. 6F). Importantly, we found that NGF is able to repress the expression of various TG-induced ER stress markers such as Atf3, Ddit3 (Chop), Hspa5, Herpub1, Ppp1r15a, etc. (Fig. 6G). This is consistent with the protein level changes of CHOP (Fig. 6H). The effect of NGF can be reversed by TrkA antibody (Fig. 6H), suggesting that the beneficial effect of NGF is mediated by the TrkA receptor.

A Heatmap showing NicheNet predicted ligand-target pairs between beta cells and other cell types. B UMAP shows restricted expression patterns of Ngf and Manf. C Published single-cell gene expression data of Manf in β cells (left), Ngf in ductal cells, and fibroblasts (right two) from normal and diabetic mouse islets. D Pancreatic NGF content in male mice (B6) fed with normal chow (NC) and HFD for 16 weeks (n = 4 (NC) or 5(HFD) biologically independent samples). Data shown as mean ± S.E.M. Analysis was done using two-tailed Student’s t-test. E RT-qPCR of Ddit3, Ppp1r15a, and Atf3 expression in HFD-fed mice islets treated with or without NGF. Data are shown as mean ± S.E.M. (n = 6 biologically independent experiments). Analysis was done using two-tailed Student’s t-test. *p < 0.05, **p < 0.01, and ***p < 0.001. F RT-qPCR of Atf3 expression in TG-stressed MIN6 cells treated with or without MANF (100 ng/ml). Data are shown as mean ± S.E.M. (n = 3 biologically independent experiments (MANF + TG), n = 4 biological independent experiments (Ctrl, MANF, TG)). Analysis was done using two-tailed Student’s t-test. *p < 0.05, **p < 0.01, and ***p < 0.001. G RT-qPCR of Atf3, Ddit3, Hspa5, Herpub1, and Ppp1r15a in TG-stressed MIN6 cells treated with or without NGF (100 ng/ml). Data shown as mean ± S.E.M. (n = 4 biological independent experiments). Analysis was done using two-tailed student’s t-test. *p < 0.05, **p < 0.01, and ***p < 0.001. H Western blotting (top) and relative band intensity (bottom) showing CHOP (Ddit3) protein levels in TG-stressed MIN6 cells treated with or without NGF. I Schematic of experimental design for NGF administration in Akita mice. Created in BioRender. Wu, J. (2023) BioRender.com/a64z386. J 16-hour fasting blood glucose level of Vehicle or NGF-treated Akita mice. (n = 5 (control) or 5 (NGF) animals per group). K Ad lib serum insulin level measured by ELISA. (n = 5 (control) or 6(NGF) animals per group). L Whole pancreas insulin content measured by ELISA. All data are shown as mean ± S.E.M. (n = 5 (control) or 4 (NGF) animals per group). Analysis was done using two-tailed Student’s t-test. *p < 0.05, **p < 0.01. Source data and exact p values are provided as a Source Data file.

In mouse islets, CHOP is a central mediator for β cell apoptosis. The fact that NGF is able to reduce CHOP protein level in MIN6 cells suggests that administration of NGF may be able to rescue β cell dysfunction induced by UPR stress. To test this, we administered recombinant NGF to Ins2WT/C96Y (Akita) mice, a genetic mouse models for severe diabetes87. These mice carry Cys96Tyr mutation in Ins2 gene, which caused a conformational change of insulin, resulting in severe misfolding stress, and eventually led to a significant reduction of total insulin (both WT and mutant) level and β cell apoptosis87. Genetic deletion of CHOP ameliorates the hyperglycemia of Akita mice88, suggesting repression of CHOP by NGF could reduce the β cell stress in this model. We injected recombinant mouse NGF intraperitoneally at a dose of 1 mg/kg89,90 daily (Fig. 6I). After 2 weeks of treatment, we found a significant reduction of fasting blood glucose in treated mice (Fig. 6J). In addition, we found that the serum insulin, as well as the total insulin content in pancreata, are increased in NGF-treated mice (Fig. 6K, L). Altogether, these data suggest that NGF is able to antagonize excessive ER stress in β cells and ameliorate β cells dysfunction.

A number of recent studies of scRNA-seq, scATAC-seq, and multiome have revealed the heterogeneity states of β cells at transcription and chromatin accessibility levels13,14,15,16,17,18,19,20,21,22,23. However, due to technical limitations, the contribution of individual cis-element states, defined by combinatory histone codes and DNA modifications, to the target gene expression has not been characterized. In this study, we take advantage of Paired-Tag, which is able to quantify histone modification instead of accessibility, to define the contribution of two related yet distinct enhancer marks, H3K4me1 and H3K27ac, to heterogenous gene expression in beta cell dysfunction. Though both marks are enriched at enhancers, each of them is separately regulated by a distinct set of epigenetic machinery. While our study has confirmed several principal discoveries of previous single-cell islet studies, it also reveals unexpected discrepancies of H3K27ac and H3K4me1 in regulating gene expression dynamics in beta-cell dysfunction, as further discussed below.

GRNs are used to study the interplay between TFs, transcription, and chromatin environment. Recent scRNA-seq and scATAC-seq datasets enable inference of GRNs in diabetic islets and differentiating β cells13,91,92, which reveals extensive links between CREs and target gene expression. In this study, we infer the GRNs from single nuclei H3K27ac/H3K4me1 and RNA. GRNs here reflect the specific linkage between the H3K27ac/H3K4me1 activity and gene expression. Although active enhancers are enriched for both H3K4me1 and H3K27ac67,68, the two modifications have distinct localization patterns and dynamics, which are determined by TF binding and the balance between HAT/HDACs and HMT/KDMs, respectively. Therefore, it is not surprising that for many HFD-induced genes, such as Atf3 and Ddit3, both H3K27ac and H3K4me1 levels are coupled with transcription dynamics along the trajectories of β cell dysfunction. However, we also identified a substantial number of genes of which the heterogeneous gene expression is only coupled with either H3K27ac or H3K4me1. Although a few published bulk ChIP-seq studies93,94 have revealed the discrepancy of H3K4me1 and H3K27ac at the bulk level, the data in our study suggest that the discordance may contribute to the heterogeneous gene expression rather than the “all-or-none” activation in the development setting. This discrepancy may be a result of the local TFs that preferentially engage in HAT/HDACs or HMTs/KDMs. For example, FoxA2 expression is highly correlative with its motif enrichment in H3K4me1 CREs but not H3K27ac CREs (Fig. 4F), suggesting that the change of FoxA2 binding is more correlated with local H3K4me1 levels, which is consistent with a previous study showing its ability to deposit H3K4me1 during pancreatic differentiation95.

The observed correlation between H3K27ac/H3K4me1 to transcription does not automatically establish causality between enhancer activity and transcription. Various research findings suggest that alterations in H3K27ac or H3K4me1 often lead to changes in selective targets rather than globally interfering with gene expression31,59,65. In islet cells, the correlation between CRE activities and gene expression varies widely among genes. Notably, genes with the most significant RNA-H3K27ac/H3K4me1 coupling are functionally distinct. This is supported by the enrichment of stress response and β cell function pathways in concordant and discordant genes (Fig. 3H–J, Fig. 4D), as well as the clustering of TF regulators (Fig. 3F, G). Furthermore, transient perturbation of global H3K27ac/H3K4me1 disrupted target stress-related expression (Fig. 3H–J, Fig. 4G), suggesting that the dynamic NF-κB and ER stress signaling in dysfunctional β cells depend on epigenetic heterogeneity. Together, these data support the model that the enhancer activity in individual cells is a significant contributor to the heterogeneity of β cell dysfunction in metabolic stress.

It is worth noting that the levels of H3K27ac and H3K4me1 at individual loci reflect the local equilibrium between HAT and HDACs, or HMT and KDMs. The transcriptional output is ultimately determined by the epigenetic complexes and TFs at specific loci in individual cells and to a lesser extent, by H3K27 acetylation itself96. However, given that the dynamics of H3K27ac and H3K4me1 mirror the dynamic balances of p300/CBP/HDACs59 and MLL3/4/LSD197, these modification levels can serve as a reliable indicator for CREs activity.

Primed enhancers are typically described in the context of developmental processes. The “priming” of these enhancers is initiated by pioneer factors such as FoxA1/2 and GATA, which localize at inaccessible nucleosomes and initiate chromatin remodeling, recruitment of secondary TFs, and eventually gene activation69. The initial stage of the primed enhancers is usually positive for H3K4me1 but not H3K27ac69. Unlike their conventional role in gene activation, the dynamics of these primed enhancers in mature islets, especially under nutritional stress, are less defined. Here we found that about 14% of H3K4me1 regions in NC islets are primed enhancers and are enriched for FoxA2 binding. Unexpectedly, these primed enhancers are largely lost under HFD stress. FoxA2 is also identified as a top TF modulator in H3K4me1-RNA multiome, but not in H3K27ac-RNA dataset, suggesting its role in regulating the heterogeneity of the primed enhancers. Nuclear FoxA2 level is significantly reduced after HFD98, which may contribute to the reduction of primed enhancers. An alternative explanation would be the redistribution of FoxA2 binding after nutritional stress. The exact mechanisms of FoxA2 in reprogramming H3K4me1 enhancers in HFD β cells will require further investigation.

What is the functional consequence of losing these primed enhancers in response to HFD? The majority of primed enhancers in islets do not significantly correlate with gene expression. The loss of H3K4me1 at these sites may signify a reallocation of the HMT to stress-induced enhancers and therefore indirectly contribute to stress-induced transcription. However, a substantial subset (14%) of primed enhancers is correlated with gene expression, specifically enriching gene categories tied to essential metabolic functions and β cell identity within the islets (Fig. 5I-K). The loss of these primed enhancers may present a distinct mechanism of cellular dysfunction.

Excessive ER stress, induced by inflammation and glucolipotoxicity, causes β cell dysfunction and contributes to the development of T2D99,100,101,102. In this study, we demonstrated that the increased expression of UPR genes in response to excessive ER stress is a consistent theme in the multiomic trajectories leading β cells to a dysfunctional state. While the adaptive UPR is beneficial to the β cell adaptation to nutrition stress, maladaptive UPR eventually leads to the reduction of β cell function genes99,100,101,102, as shown in the most stressed/dysfunctional cells in our multiomic trajectories. Chronic excessive ER stress can be reversible, offering a potential avenue for pharmacological intervention. However, currently, there is a lack of pharmacological approaches to effectively mitigate damage from excessive ER stress in β cells. Through an investigation of cell-cell communication, we identified several secretory factors, including nerve growth factor (NGF), that repress ER stress markers in β cells. While NGF is well-characterized in neurite growth, our data indicate it can also act on β cells to specifically repress ER stress genes. NGF expression in fibroblasts and ductal cells within islets suggests a local paracrine signaling mechanism that may protect β cells from excessive ER stress. We demonstrated that recombinant NGF effectively blocks ER stress in HFD islets or in MIN-6 cells stressed by synthetic inducers thapsigargin (TG), and this action required the TrkA receptor (Fig. 6G, H, Supplementary Data 6). Most importantly, when treated with recombinant NGF, Akita mice which are highly diabetic due to excessive ER stress, showed reduced blood glucose, increased serum insulin, and increased pancreatic insulin content (Fig. 6J–L). Mechanistically, this may partly be due to NGF’s ability to suppress CHOP (Ddit3), which is a key factor inducing hyperglycemia in the Akita mouse model88. While classic paracrine hormones such as glucagon, somatostatin (SST), and urocortin-3 (UCN3) are extensively studied103, other cell-cell communications specific to β cell dysfunction are less explored. Given NGF’s potent ability to repress ER stress in β cells, further investigation into the mechanistic insights and potential pharmacological modulation of this paracrine pathway may prove valuable.

Defining the heterogenous beta cell states has been challenging, as previous literature suggests numerous ways of biological distinct subpopulations exist in healthy and dysfunctional islets. Our initial classification by Seurat reveals up to 6 potential subpopulations in the main beta cell clusters (Supplementary Fig. 3B). However, the challenge is to determine the biological characteristics for individual clusters. In addition, it is challenging to define the trajectories of 6 clusters in beta cell dysfunction. Therefore, we considered the definition in multiple published single cell studies13,14,15,16,17,45, and merged 6 populations into two major β cell states. It should be noted that this dichotomy is very likely an oversimplification of the complex beta cell heterogeneity. Future studies, especially with spatial-omics, single-cell electrophysiology, and lineage tracing, will be needed to further define the potential substates of beta cells.

In conclusion, we presented a comprehensive analysis of enhancer dynamics in islets from lean and obese mice. We identified both shared and distinct gene regulatory mechanisms by H3K27ac and H3K4me1 enhancers. Our results suggest that β cells exhibit substantial heterogeneity not only in their transcriptional profile but also at the enhancer level. The heterogeneity is partly driven by the inflammatory and ER stress, contributing to the dysfunction of islets. Interestingly, although both H3K27ac and H3K4me1 are enhancer marks and correlated to gene activation, our study reveals unique dynamics and dependencies for these two marks in regulating specific gene subsets. Surprisingly, the primed enhancers in lean islets are globally lost under HFD stress. The loss is connected to the decreased expression of genes vital for β cell function, indicating a significant role for primed enhancers and pioneer factors in β cell dysfunction. Lastly, we combined cell-cell communication analysis and pharmacological perturbation to identify paracrine factors regulating HFD-induced stress in β cells, and discovered that NGF is able to potently suppress the ER stress in β cells and ameliorate hyperglycemia in Akita mice. NGF, previously uncharacterized in islets, could be a target for therapeutic interventions to alleviate ER stress-related β cell dysfunction. Taken together, our results suggest distinct epigenetic regulatory modes, defined by combinations of TFs and epigenetic modulators, orchestrate β cell dysfunction. Future investigations will be needed to further dissect the heterogeneous enhancer states in cellular dysfunction.

The Institutional Animal Care and Use Committee (IACUC) at Mayo Clinic Arizona deemed that this research complies with all relevant ethical regulations (Protocol no. A00005200-20-R23).

SAHA (T1583), GNE781 (T15405), GSK-LSD1 (T2315), Thapsigargin (TQ0302) were from Targetmol (Boston, MA, USA). Native mouse NGF protein (N-100) was from Alomone Labs. Murine IL-1β (211-11B) and human MANF recombinant protein (450-06) were purchased from PeproTech (Cranbury, NJ, USA).

Rabbit monoclonal anti-H3K4me1 (D1A9), rabbit monoclonal anti-H3K27ac (D5E4), mouse monoclonal anti-CHOP (L63F7), rabbit monoclonal anti-α/β-Tubulin, HRP-linked anti-rabbit IgG (#7074), and HRP-linked anti-mouse IgG (#7076) antibodies were from Cell Signaling Technology (Danvers, MA, USA). Mouse monoclonal anti-TrkA antibody (clone MNAC13) was from Sigma (St Louis, MO, USA).

MIN6 cells were maintained in Dulbecco’s modified Eagle’s medium (DMEM) (Corning, 10-013-CV) containing 1 mM sodium pyruvate (Corning, 25-000-CI), 15% fetal bovine serum (Corning, 35-015-CV), 1% penicillin-streptomycin (Gibco, 15140-122), and 55 uM 2-mercaptoethanol (Gibco, 21-985-023). For experiments in Figs. 3 and 4, MIN6 cells were pretreated with GNE 781 (1 uM), SAHA (5 uM), or GSK-LSD1 (5 uM) respectively for 24 h. For IL-1β stimulation, MIN6 cells were pretreated with GNE 781 (1 uM), SAHA (5 uM), or GSK-LSD1 (5 uM) respectively for 24 h, followed by treatment of IL-1β (5 ng/ml) for 5 h. For experiments in Fig. 6, MIN6 cells were pretreated with NGF (100 ng/ml) for 3 h or MANF (200 ng/ml) for 24 h, followed by treatment of TG (0.5 uM) for 12 h.

All mice were housed in a 12-h light-dark cycle at an ambient temperature of 23 °C. Mouse experiments were approved by the Institutional Animal Care and Use Committee (IACUC) at Mayo Clinic Arizona. C57BL/6J mice were from Jackson laboratory (cat. no. 000664). For Paired-tag, ten female or male mice (n = 10) on postnatal week 12 were used for islet isolation. For high-fat diet treatment, mice were fed with 60% fat diet paste (BioServ, S1850) for 8 weeks and 16 weeks before tissue was collected. Akita mice (male) were from Jackson Laboratory (Cat. no. 003548). Akita mice (male) on postnatal day 84 (week 12) were randomly grouped and received NGF (1 mg/kg/mouse) or vehicle (Control) (PBS) daily via intraperitoneal injections for 16 days. For fasting blood glucose, Akita mice were fasted for 16 h before monitored blood glucose levels. Serum insulin levels were measured via the Mouse Ultrasensitive Insulin ELISA kit (ALPCO, Catalog 80-INSMSU-E10). For pancreatic insulin content, insulins were extracted via acid-ethanol extraction from the whole pancreas isolated from Akita mice, and measured by Mouse Ultrasensitive Insulin ELISA kit (ALPCO, Catalog 80-INSMSU-E10). Insulin content was normalized by the total protein amount measured by BCA protein assay kit (Pierce, 23225).

Mouse pancreatic islets were isolated as previously described104. 0.5 mg/ml collagenase P (Sigma, 11249002001) diluted in HBSS buffer (Fisher, 14170161) and injected through the common bile duct. The perfused pancreas was collected and incubated at 37 °C for 21 min. Digested exocrine cells and intact islets were separated via centrifugation over Histopaque-1077 (Sigma, 10771), and intact islets were manually selected.

We followed the published protocol35,105 with minor modifications. Briefly, after nuclei isolation with nuclei isolation buffer: 0.2% IGEPAL CA-630 (Sigma #I8896), 5% BSA (Roche #03117057001) and 1 mM DTT (Sigma # D9779) in PBS (Thermo Fisher Scientific #10010-23), supplemented with 1× Proteinase Inhibitor (Roche #05056489001), 0.5 U/µL SUPERaseIn (Invitrogen # AM2694), and 0.5 U/µL RNase OUT (Invitrogen #10777019), each 300,000 of nuclei were aliquot into the 12 1.5 mL low-bind tubes. Nuclei were spin-down and resuspended in 30 µL MED#1 buffer: 20 mM HEPES (Gibco #15630106), 300 mM NaCl (Sigma #S7653), 0.5 mM Spermidine (Sigma #85558), 1× Proteinase Inhibitor, 0.5 U/µL SUPERase In, 0.5 U/µL RNase OUT, 0.01% IGEPAL-CA630, 0.01% Digitonin (Sigma #D141), 2 mM EDTA (Invitrogen #AM9261)) and keep on ice.

2 µg of H3K4me1 antibody (Cell Signaling #5326S) or H3K27ac antibody (Cell Signaling #8173S) were added into 12 of 200 µL tubes containing 20 µL MED#1 buffer, and pA-Tn5 protein (purified in house) with 12 DNA barcodes were added and incubated at room temperature with gentle rotation for 1 h. The 12 tubes of antibody-pA-Tn5 mix were then mixed with each tube of nuclei, respectively and the incubation was carried out in 4 °C with gentle rotation overnight. The nuclei were then spun-down and washed two times with MED#2 buffer (20 mM HEPES, 300 mM NaCl, 0.5 mM Spermidine, 1× Proteinase Inhibitor, 0.5 U/µL SUPERase In, 0.5 U/µL RNase OUT, 0.01% IGEPAL-CA630, 0.01% Digitonin) and resuspended in 50 µL MED#2 buffer. Tagmentation reaction were activated out by adding 2 µL of 250 mM MgCl2 (Sigma #63069), carried out in a ThermoMixer set at 37 °C, 550 rpm for 60 min and quenched by adding 16.5 µL of 40 mM EDTA. Nuclei were then spin-down and reverse transcription was carried out with Maxima H minus reverse transcriptase (Thermo #EP0751). Nuclei were then barcoded by ligation-based combinatorial barcoding with T4 DNA Ligase (NEB #M0202L), aliquoted into 3.5k nuclei sub-libraries, and lysed. Library preparation were then carried out and sequenced with read cycles 150 (read1) + 6 (index1) + 6 (index2) + 150 (read2) on a NovaSeq 6000 platform.

Isolated C57BL/6J HFD-fed mice islets were treated with SAHA (5 uM), C646 (5 uM), or GNE 781(1 uM) for 24 h. Treated islets were washed with Versene (Gibco, 15040-066) and dissociated using TrypLE Express (Gibco, 12605-010) for 10 min at 37 °C. Dead cells were gated and removed using DAPI staining (1:3000, Cayman Chemical) on a BD LSRFortessa flow cytometer with FacsDIVA software (version 9.0). Beta cells were isolated based on their granularity and autofluorescence of FAD/FMN (FITC) at wavelength 488 nm. Cell populations were identified with CD63-APC monoclonal antibody (1:50, eBioscience) where 50% of CTL beta cells set the relative benchmark for differentiating high from low. Data was analyzed using FlowJo, BD Biosciences).

ELISA were performed per manufacturer’s instructions: mouse insulin (Alpco, 80-INSMS-E01), mouse NGF(RayBiotech, ELM-bNGF-1).

Total RNA was extracted as previously described106 with slight modification. Briefly, MIN6 cells (AddexBio, C0018008) in 96-well plates were washed once with PBS (Fisher, SH30256.LS), and incubated with 20 ul/well of lysis buffer (10 mM Tris-HCl, 0.5% NP40, 150 nM NaCl) for 5 min. 12 ul cell lysate was used for reverse transcription. Reverse transcription was performed with AzuraFle cDNA Synthesis Kit (Azura, AZ-1998). Real-time quantitative PCR with reverse transcription (RT-qPCR) was performed using Forget-Me-Not EvaGreen qPCR Master Mix (Biotium, 31041). Data was normalized to the levels of the internal control Gapdh and showed as the mean with SEM error bars. Statistical comparisons were made using Student’s t-test.

MIN6 cells seeded in 6-well plate were washed with ice-cold PBS and lysed using ice-cold RIPA buffer (Sigma, R0278) with protease inhibitor (Sigma, 11836145001). Lysates were mixed with sample loading buffer and denatured at 95 °C for 10 min. Samples were analyzed with SDS-PAGE on 4-20% Mini-PROTEAN TGX precast protein gels (Bio-Rad, 456-1033) with equal amount of protein loading. Proteins were visualized by transferring onto PVDF membrane and incubating with anti-CHOP or anti-α/β-Tubulin primary antibodies and horseradish peroxidase (HRP)-linked secondary antibodies.

Preprocessing of data were carried out with the scripts available from Paired-Tag35. Briefly, cellular barcodes were extracted from Read2 and assigned to each sample barcodes (12 initial tubes for tagmentation and reverse transcription) and combination of ligated barcodes.

Cellular barcodes and the linker sequences are read by Read 2. The first bases of barcode (BC) no. 1, BC no. 2, and BC no. 3 should locate within the 84th−87th, 47th−50th and 10th−13th bases of Read 2. We identify the positions of barcodes by matching the linker sequences adjacent to the cellular barcodes. A bowtie reference index was generated with all possible cellular barcode combinations (96 × 96 × 12) and barcode sequences were mapped to the cellular barcode reference using bowtie107 with the parameters: -v 1 -m 1 --norc (reads with more than 1 barcode mismatch and can be assigned to more than 1 cell were discarded). NextEra adaptor sequences were trimmed from 3′ of DNA and RNA libraries, Poly-dT sequences were further trimmed from 3′ of RNA libraries, and low-quality reads were excluded from further analysis.

Cleaned reads were first mapped to a mouse GRCm38 reference genome with STAR108 (v.2.7.3a) for RNA or bowtie2109 for DNA. Duplicates were removed based on the mapped position, cellular barcode, PCR index, and UMI. Before generating the cell-counts matrices, DNA bam files were further filtered by removing high-pileup positions (cutoff = 10) regardless of cellular barcode, PCR index, and UMI. RNA alignment files were converted to a matrix with cells as columns and genes as rows. DNA alignment files were converted to a matrix with cells as columns and 5-kb bins as rows. In addition, the annotated DNA bam files were converted to fragment files using scATACpro’s simply_bam2frags_singleEnd.pl script, which are bed-like matrices containing chromosome, start, stop, cell ID, and number of fragments contained in the region110.

Cells with fewer than 200 features in RNA were removed. Seurat (v4.3.0)111 R package was used to perform filtering, normalization, dimensionality reduction, clustering, and differential expression analysis. The following criteria were applied to each sample to remove low-quality cells: gene number between 200 and 3000, and mitochondrial content <10%. Doublets were predicted by the DoubletFinder (v2.0.3) algorithm112. After filtering, a total of 56,748 cells were left. For integration, the 3000 most conserved genes in all samples were identified with SelectIntegrationFeatures function. Subsequently, the PrepSCTIntegration function was used to identify conserved features for integrating the datasets. Integration anchors were identified using FindIntegrationAnchors function with a scTransform normalization method. Seurat’s IntegrateData function was implemented to integrate the datasets. DEG between conditions within each cluster were identified using the FindMarkers function. The first 20 principal components were chosen to reduce dimensions. The first two uniform manifold approximation and projection (UMAP) dimensions were used to visualize cell clusters. Cell clusters were identified using FindClusters (20 principal components and 0.3 resolution), and the cell clusters expressing the same marker genes will be merged. The β cell cluster was used for next-step pseudo- or latent time analysis. Density plot, produced by the Nebulosa113 package.

Comparisons between selected cell populations were performed using NicheNet78 functions and standard workflow.

RNA velocity analysis was based on the initial processing of the Paired-tag RNA sequencing data using scVelo (v0.2.5)40 and further analysis using CellRank41. We used Velocyto114 to generate a loom file from Paired-tag output data (‘velocyto run’, see also http://velocyto.org/velocyto.py/tutorial/cli.html) that differentiates between spliced, unspliced, and ambiguous gene counts. This loom file was used for pre-processing and clustering via Seurat(v4.3.0)111. scVelo(v0.2.5) was used to estimate trajectories based on spliced vs unspliced RNA for each cell and gene. We used scVelo to compute RNA velocity, then set up CellRank’s VelocityKernel and computed a transition matrix based on RNA velocity. We have combined the VelocityKernel with the ConnectivityKernel to emphasize gene expression similarity and visualized the transition matrix in a low-dimensional embedding (UMAP). CellRank41 computed terminal states, the number of terminal states was defined by using the cr.tl.terminal_states() function. Then, it defined initial states with the cr.tl.initial_states() function. Next, fate maps were computed to define the likelihood that an individual cell reaches a certain fate. Subsequently, the individual fate maps were aggregated to cluster-level fate maps by calculating latent time and running PAGA115. Gene expression trends along pseudo- or latent time were calculated based on a general additive model (GAM) with the cr.pl.gene_trends() function. Heatmaps of top 1000 genes whose expression levels correlated best with the probability of reaching the terminal stage were plotted using the cr.pl.heatmap() function.

Paired-Tag data for mouse islet cells was processed using Signac and Seurat. We retained cells with total RNA counts >200 and total single-cell DNA profiles bin counts >500. Functions within Seurat were used for analysis. RNA counts were normalized with SCTransform with mitochondrial percent per cell regressed out. Principal component analysis (PCA) was performed on RNA, and UMAP was run on the first 30 principal components (PCs). The optimum number of PCs was determined to be 30 PCs using an elbow plot. The single-cell DNA profiles peaks counts were normalized with term-frequency inverse-document-frequency (TFIDF). Dimension reduction was performed with singular value decomposition (SVD) of the normalized single-cell DNA peaks matrix. The single-cell DNA profiles UMAP was created using the second through the 30th LSI components. The weighted nearest neighbor (WNN) graph was determined with Seurat’s FindMultiModalNeighbors to represent a weighted combination of both modalities. The first 30 dimensions of the Paired-Tag RNA reduction and the second through the 30th dimensions from the Paired-Tag DNA profiles reduction were used to create the graph.

Differentially expressed genes (DEGs) were determined for HFD versus control for each cell type. Within each cell type, the gene expression data was log-normalized with a scale factor of 1×105. Differential expression was assessed using MAST116 for genes present in at least 25% of either HFD or control cells. Week and sex were included as covariates in the MAST model. Genes with a Bonferroni-adjusted p value < 0.01 and an absolution log2(fold change)>0.25 were determined to be significant.

The R package enrichR117,118 was used for all gene set enrichment analyses. Sets of DEGs were used as input to look for enrichment in GO Biological Process 2021, GO Molecular Function 2021, GO Cellular Component 2021, and KEGG 2021 databases. Terms with an adjusted p-value < 0.05 were considered to be enriched.

To infer peak-gene putative regulatory links, we used FigR 0.1.052, restricted to the normalized and corrected counts of the 10 Paired-tag samples with both RNA and DNA profiles available. We take a default (10 kb) window around each gene’s TSS, and compute the Spearman correlation across all cells between their peak accessibility counts (mean-centered) and the normalized RNA expression. For each peak-gene pair correlation, we use (default n = 100) background peaks matched for GC content and accessibility to correlate to the same gene, so that we can test for significance (permutation p-value). Downstream analyses were performed using gene-peak pairs with a p-value < 0.05. To infer regulatory interactions, we used the getDORCScores and runFigRGRN functions, on all genes with more than 10 linked peaks. TF-gene pairs with an absolute regulation score greater than the default (>1) are considered putative interaction.

Known and de novo motif analyses of H3K4me1/H3K27ac primed enhancer peaks were performed with HOMER119.

We performed cellular trajectory analyses as previously described using ArchR53. We ran the createArrowFiles and ArchRProject functions to create the object in R required for analysis, followed by addPeakSet and addPeakMatrix for quantification. This was performed using one core by setting addArchRThreads(threads = 1). We defined the trajectory backbone of cell groups or clusters based on the beta cell states, such as the order of beta-hi, and beta-low cells in the beta cell trajectory. We then created a trajectory using the ‘addTrajectory’ function and plotted the pseudotime values on UMAP embedding using the ‘plotTrajectory’ function. We next plotted pseudotime heatmaps of TFs, gene scores, and peak accessibility using the ‘plotTrajectoryHeatmap’ function. To identify positive TF regulators, we performed an integrative analysis of gene scores and motif accessibility across pseudotime using the ‘correlateTrajectories’ and plotTrajectoryHeatmap functions.

FoxA2 ChIP-seq data was obtained from the GEO (GSM1306337). H3K4me1 ChIP-seq data was obtained from the NCBI SRA (SRR9840908, SRR9840909). H3K27ac ChIP-seq data of mouse pancreatic islets was obtained from the NCBI SRA (SRR9840900, SRR9840901). H3K27ac ChIP-seq data in EndoC-βH1 cells was obtained from the NCBI SRA (SRR9336435, SRR9336436, SRR9336437, SRR9336438, SRR9336439, SRR9336440, SRR9336441, SRR9336442). Bowtie2 (v2.3.3.1) (-very-sensitive) was used to map ChIP-seq reads to the mouse reference genome GRCm38 or human reference genome GRCh38. Duplicate reads were filtered out using the MarkDuplicate function from Picardtools v.2.17.0 (http://broadinstitute.github.io/picard/). Reads per kilobase and million mapped read (RPKM)-normalized bigWig files were generated with bamCoverage from deepTools v3.3.2.

The GWAS data was obtained from the NHGRI-EBI GWAS Catalog120 using “gwasrapidd”121(v0.99.17) R package. We queried the GWAS Catalog for studies on type 2 diabetes mellitus using the get_studies() function. This search returned 226 GWAS studies related to type 2 diabetes mellitus. Function get_associations() was used to get association IDs, p-values, and odds ratio for every variant. Out of the 5438 associations found in the GWAS Catalog, 4871 meet the p-value threshold of 1 × 10−6, leading to the identification of 678 genes. Next, we compiled a list of genes associated with human islet function from two recently published scRNAs-eq studies16,45.

RNA-seq data from mouse islets and PP2 cells was obtained from the NCBI SRA (SRR6159442, SRR6159443, SRR6159444, SRR6159445, SRR6159446, SRR10708482). The reads were then mapped against GRCm38 or GRCh38 reference genome using STAR (v 2.7.9a)108. Gene counts were normalized using the variance stabilizing transformation (VST) function of DESeq2 (v1.38.3) R package122. Heatmaps of gene expression were generated from the Z-scores of VST gene counts using ComplexHeatmap (v2.15.4)123.

Sample sizes were not predetermined by statistical methods. Investigators were not blinded to allocation during experiments and outcome assessment. GraphPad Prism was used for analysis and graph plotting. The statistical methods are listed in figure legends.

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.

The sequencing data obtained in this study have been deposited at the NCBI Gene Expression Omnibus (GEO) (http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE241639) under accession number GSE241639. This paper also analyses existing, publicly available data. FoxA2 ChIP-seq data was obtained from the GEO (GSM1306337). H3K4me1 ChIP-seq data was obtained from the NCBI SRA (SRR9840908, SRR9840909). H3K27ac ChIP-seq data of mouse pancreatic islets was obtained from the NCBI SRA (SRR9840900, SRR9840901). H3K27me3 ChIP-seq data was obtained from the NCBI SRA (SRR6728238). H3K27ac ChIP-seq data in EndoC-βH1 cells was obtained from the NCBI SRA (SRR9336435, SRR9336436, SRR9336437, SRR9336438, SRR9336439, SRR9336440, SRR9336441, SRR9336442). RNA-seq data from mouse islets and PP2 cells was obtained from the NCBI SRA (SRR6159442, SRR6159443, SRR6159444, SRR6159445, SRR6159446, SRR10708482). Ntrk1 RNA-seq data was obtained from the GEO (GSE80673). Source data are provided with this paper.

No new software was developed during this study. Code used for processing mouse Paired-tag data sets is available at https://github.com/256wangliu/Multiomics_diabetes.

Ashcroft, F. M. & Rorsman, P. Diabetes mellitus and the beta cell: the last ten years. Cell 148, 1160–1171 (2012).

Article CAS PubMed PubMed Central Google Scholar

Prentki, M., Peyot, M. L., Masiello, P. & Madiraju, S. R. M. Nutrient-induced metabolic stress, adaptation, detoxification, and toxicity in the pancreatic beta-cell. Diabetes 69, 279–290 (2020).

Article CAS PubMed Google Scholar

Donath, M. Y., Dalmas, E., Sauter, N. S. & Boni-Schnetzler, M. Inflammation in obesity and diabetes: islet dysfunction and therapeutic opportunity. Cell Metab. 17, 860–872 (2013).

Article CAS PubMed Google Scholar

Haythorne, E. et al. Diabetes causes marked inhibition of mitochondrial metabolism in pancreatic beta-cells. Nat. Commun. 10, 2474 (2019).

Article ADS PubMed PubMed Central Google Scholar

Kim, H. & Kulkarni, R. N. Epigenetics in beta-cell adaptation and type 2 diabetes. Curr. Opin. Pharmacol. 55, 125–131 (2020).

Article CAS PubMed PubMed Central Google Scholar

Lee, Y. S., Wollam, J. & Olefsky, J. M. An integrated view of immunometabolism. Cell 172, 22–40 (2018).

Article CAS PubMed PubMed Central Google Scholar

Butler, A. E. et al. Beta-cell deficit and increased beta-cell apoptosis in humans with type 2 diabetes. Diabetes 52, 102–110 (2003).

Article CAS PubMed Google Scholar

Weir, G. C. & Bonner-Weir, S. Five stages of evolving beta-cell dysfunction during progression to diabetes. Diabetes 53, S16–S21 (2004).

Article CAS PubMed Google Scholar

Roscioni, S. S., Migliorini, A., Gegg, M. & Lickert, H. Impact of islet architecture on beta-cell heterogeneity, plasticity and function. Nat. Rev. 12, 695–709 (2016).

CAS Google Scholar

Benninger, R. K. P. & Kravets, V. The physiological role of beta-cell heterogeneity in pancreatic islet function. Nat. Rev. 18, 9–22 (2022).

CAS Google Scholar

Camunas-Soler, J. et al. Patch-Seq links single-cell transcriptomes to human islet dysfunction in diabetes. Cell Metab. 31, 1017–1031.e1014 (2020).

Article CAS PubMed PubMed Central Google Scholar

Bader, E. et al. Identification of proliferative and mature beta-cells in the islets of Langerhans. Nature 535, 430–434 (2016).

Article ADS CAS PubMed Google Scholar

Wang, G. et al. Integrating genetics with single-cell multiomic measurements across disease states identifies mechanisms of beta cell dysfunction in type 2 diabetes. Nat. Genet. 55, 984–994 (2023).

Article CAS PubMed PubMed Central Google Scholar

Chiou, J. et al. Single-cell chromatin accessibility identifies pancreatic islet cell type- and state-specific regulatory programs of diabetes risk. Nat. Genet 53, 455–466 (2021).

Article CAS PubMed PubMed Central Google Scholar

Rubio-Navarro, A. et al. A beta cell subset with enhanced insulin secretion and glucose metabolism is reduced in type 2 diabetes. Nat. Cell Biol. 25, 565–578 (2023).

Article CAS PubMed PubMed Central Google Scholar

Fang, Z. et al. Single-cell heterogeneity analysis and CRISPR screen identify key beta-cell-specific disease genes. Cell Rep. 26, 3132–3144.e3137 (2019).

Article CAS PubMed PubMed Central Google Scholar

Dror, E. et al. Epigenetic dosage identifies two major and functionally distinct beta cell subtypes. Cell Metab. 35, 821–836.e827 (2023).

Article CAS PubMed PubMed Central Google Scholar

Shapira, S. N., Naji, A., Atkinson, M. A., Powers, A. C. & Kaestner, K. H. Understanding islet dysfunction in type 2 diabetes through multidimensional pancreatic phenotyping: The Human Pancreas Analysis Program. Cell Metab. 34, 1906–1913 (2022).

Article CAS PubMed PubMed Central Google Scholar

Qiu, W. L. et al. Deciphering pancreatic islet beta cell and alpha cell maturation pathways and characteristic features at the single-cell level. Cell Metab. 25, 1194–1205 e1194 (2017).

Article CAS PubMed Google Scholar

Baron, M. et al. A single-cell transcriptomic map of the human and mouse pancreas reveals inter- and intra-cell population structure. Cell Syst. 3, 346–360.e344 (2016).

Article CAS PubMed PubMed Central Google Scholar

Li, J. et al. Single-cell transcriptomes reveal characteristic features of human pancreatic islet cell types. EMBO Rep. 17, 178–187 (2016).

Article CAS PubMed Google Scholar

Shrestha, S. et al. Combinatorial transcription factor profiles predict mature and functional human islet alpha and beta cells. JCI Insight 6, https://doi.org/10.1172/jci.insight.151621 (2021).

Augsornworawat, P. et al. Single-nucleus multi-omics of human stem cell-derived islets identifies deficiencies in lineage specification. Nat. Cell Biol. 25, 904–916 (2023).

Article CAS PubMed PubMed Central Google Scholar

Badia, I. M. P. et al. Gene regulatory network inference in the era of single-cell multi-omics. Nature reviews https://doi.org/10.1038/s41576-023-00618-5 (2023).

Liu, J. S. & Hebrok, M. All mixed up: defining roles for beta-cell subtypes in mature islets. Genes Dev. 31, 228–240 (2017).

Article CAS PubMed PubMed Central Google Scholar

Cefalu, W. T. et al. Heterogeneity of Diabetes: beta-Cells, Phenotypes, and Precision Medicine: Proceedings of an International Symposium of the Canadian Institutes of Health Research’s Institute of Nutrition, Metabolism and Diabetes and the U.S. National Institutes of Health’s National Institute of Diabetes and Digestive and Kidney Diseases. Diabetes https://doi.org/10.2337/db21-0777 (2021).

Dai, X. Q. et al. Heterogenous impairment of alpha cell function in type 2 diabetes is linked to cell maturation state. Cell Metab. 34, 256–268.e255 (2022).

Article CAS PubMed PubMed Central Google Scholar

Kim, S. & Wysocka, J. Deciphering the multi-scale, quantitative cis-regulatory code. Mol. Cell 83, 373–392 (2023).

Article CAS PubMed Google Scholar

Consortium, E. P. et al. Expanded encyclopaedias of DNA elements in the human and mouse genomes. Nature 583, 699–710 (2020).

Article ADS Google Scholar

Calo, E. & Wysocka, J. Modification of enhancer chromatin: what, how, and why? Mol. cell 49, 825–837 (2013).

Article CAS PubMed Google Scholar

Wortham, M. et al. Nutrient regulation of the islet epigenome controls adaptive insulin secretion. J. Clin. Investig., 133, https://doi.org/10.1172/JCI165208 (2023).

Scoville, D. W. et al. MLL3 and MLL4 Methyltransferases bind to the MAFA and MAFB transcription factors to regulate islet beta-cell function. Diabetes 64, 3772–3783 (2015).

Article CAS PubMed PubMed Central Google Scholar

Wong, C. K. et al. The p300 and CBP transcriptional coactivators are required for beta-cell and alpha-cell proliferation. Diabetes 67, 412–422 (2018).

Article CAS PubMed Google Scholar

Zhang, L. et al. CBP/p300 HAT maintains the gene network critical for beta cell identity and functional maturity. Cell Death Dis. 12, 476 (2021).

Article CAS PubMed PubMed Central Google Scholar

Zhu, C. et al. Joint profiling of histone modifications and transcriptome in single cells from mouse brain. Nat. Methods 18, 283–292 (2021).

Article CAS PubMed PubMed Central Google Scholar

Xiong, H., Luo, Y., Wang, Q., Yu, X. & He, A. Single-cell joint detection of chromatin occupancy and transcriptome enables higher-dimensional epigenomic reconstructions. Nat. methods 18, 652–660 (2021).

Article CAS PubMed Google Scholar

Mosser, R. E. et al. High-fat diet-induced beta-cell proliferation occurs prior to insulin resistance in C57Bl/6J male mice. Am. J. Physiol. 308, E573–E582 (2015).

CAS Google Scholar

Hull, R. L. et al. Dietary-fat-induced obesity in mice results in beta cell hyperplasia but not increased insulin release: evidence for specificity of impaired beta cell adaptation. Diabetologia 48, 1350–1358 (2005).

Article CAS PubMed Google Scholar

Kleinert, M. et al. Animal models of obesity and diabetes mellitus. Nat. Rev. 14, 140–162 (2018).

Google Scholar

Bergen, V., Lange, M., Peidli, S., Wolf, F. A. & Theis, F. J. Generalizing RNA velocity to transient cell states through dynamical modeling. Nat. Biotechnol. 38, 1408–1414 (2020).

Article CAS PubMed Google Scholar

Lange, M. et al. CellRank for directed single-cell fate mapping. Nat. Methods 19, 159–170 (2022).

Article CAS PubMed PubMed Central Google Scholar

Miranda, M. A., Macias-Velasco, J. F., Schmidt, H. & Lawson, H. A. Integrated transcriptomics contrasts fatty acid metabolism with hypoxia response in beta-cell subpopulations associated with glycemic control. BMC Genom. 24, 156 (2023).

Article CAS Google Scholar

Tritschler, S. et al. A transcriptional cross species map of pancreatic islet cells. Mol. Metab. 66, 101595 (2022).

Article CAS PubMed PubMed Central Google Scholar

Hao, Y. et al. Integrated analysis of multimodal single-cell data. Cell 184, 3573–3587.e3529 (2021).

Article CAS PubMed PubMed Central Google Scholar

Weng, C. et al. Single cell multiomic analysis reveals diabetes-associated beta-cell heterogeneity driven by HNF1A. Nat. Commun. 14, 5400 (2023).

Article ADS CAS PubMed PubMed Central Google Scholar

Zhang, T. et al. FoxO1 plays an important role in regulating beta-cell compensation for insulin resistance in male mice. Endocrinology 157, 1055–1070 (2016).

Article CAS PubMed PubMed Central Google Scholar

Rulifson, I. C. et al. Wnt signaling regulates pancreatic beta cell proliferation. Proc. Natl Acad. Sci. USA 104, 6247–6252 (2007).

Article ADS CAS PubMed PubMed Central Google Scholar

Jackerott, M. et al. STAT5 activity in pancreatic beta-cells influences the severity of diabetes in animal models of type 1 and 2 diabetes. Diabetes 55, 2705–2712 (2006).

Article CAS PubMed Google Scholar

Wei, Z. et al. Vitamin D switches BAF complexes to protect beta cells. Cell 173, 1135–1149.e1115 (2018).

Article CAS PubMed PubMed Central Google Scholar

Gannon, M., Kulkarni, R. N., Tse, H. M. & Mauvais-Jarvis, F. Sex differences underlying pancreatic islet biology and its dysfunction. Mol. Metab. 15, 82–91 (2018).

Article CAS PubMed PubMed Central Google Scholar

Liu, G. et al. Single-cell RNA sequencing reveals sexually dimorphic transcriptome and Type 2 Diabetes genes in mouse islet beta cells. Genom. Proteom. Bioinforma. 19, 408–422 (2021).

Article Google Scholar

Kartha, V. K. et al. Functional inference of gene regulation using single-cell multi-omics. Cell Genom 2, https://doi.org/10.1016/j.xgen.2022.100166 (2022).

Granja, J. M. et al. ArchR is a scalable software package for integrative single-cell chromatin accessibility analysis. Nat. Genet. 53, 403–411 (2021).

Article CAS PubMed PubMed Central Google Scholar

Tessem, J. S. et al. Nkx6.1 regulates islet beta-cell proliferation via Nr4a1 and Nr4a3 nuclear receptors. Proc. Natl Acad. Sci. USA 111, 5242–5247 (2014).

Article ADS CAS PubMed PubMed Central Google Scholar

Win, P. W., Oakie, A., Li, J. & Wang, R. Beta-cell beta1 integrin deficiency affects in utero development of islet growth and vascularization. Cell tissue Res. 381, 163–175 (2020).

Article CAS PubMed Google Scholar

Diaferia, G. R. et al. beta1 integrin is a crucial regulator of pancreatic beta-cell expansion. Development 140, 3360–3372 (2013).

Article CAS PubMed PubMed Central Google Scholar

Tixi, W. et al. Coordination between ECM and cell-cell adhesion regulates the development of islet aggregation, architecture, and functional maturation. bioRxiv, 2022.2004.2027.489466 https://doi.org/10.1101/2022.04.27.489466 (2023).

Kelly, W. K. & Marks, P. A. Drug insight: Histone deacetylase inhibitors–development of the new targeted anticancer agent suberoylanilide hydroxamic acid. Nat. Clin. Pr. Oncol. 2, 150–157 (2005).

Article CAS Google Scholar

Raisner, R. et al. Enhancer activity requires CBP/P300 Bromodomain-dependent Histone H3K27 acetylation. Cell Rep. 24, 1722–1729 (2018).

Article CAS PubMed Google Scholar

Egolf, S. et al. LSD1 inhibition promotes epithelial differentiation through derepression of fate-determining transcription factors. Cell Rep. 28, 1981–1992.e1987 (2019).

Article CAS PubMed PubMed Central Google Scholar

Richon, V. M., Sandhoff, T. W., Rifkind, R. A. & Marks, P. A. Histone deacetylase inhibitor selectively induces p21WAF1 expression and gene-associated histone acetylation. Proc. Natl Acad. Sci. USA 97, 10014–10019 (2000).

Article ADS CAS PubMed PubMed Central Google Scholar

Wang, M., Chen, Z. & Zhang, Y. CBP/p300 and HDAC activities regulate H3K27 acetylation dynamics and zygotic genome activation in mouse preimplantation embryos. EMBO J. 41, e112012 (2022).

Article CAS PubMed PubMed Central Google Scholar

Bowers, E. M. et al. Virtual ligand screening of the p300/CBP histone acetyltransferase: identification of a selective small molecule inhibitor. Chem. Biol. 17, 471–482 (2010).

Article CAS PubMed PubMed Central Google Scholar

Romero, F. A. et al. GNE-781, A highly advanced potent and selective bromodomain inhibitor of cyclic adenosine monophosphate response element binding protein, binding protein (CBP). J. Med Chem. 60, 9162–9183 (2017).

Article CAS PubMed Google Scholar

Vinckier, N. K. et al. LSD1-mediated enhancer silencing attenuates retinoic acid signalling during pancreatic endocrine cell development. Nat. Commun. 11, 2082 (2020).

Article ADS CAS PubMed PubMed Central Google Scholar

Ramos-Rodriguez, M. et al. The impact of proinflammatory cytokines on the beta-cell regulatory landscape provides insights into the genetics of type 1 diabetes. Nat. Genet. 51, 1588–1595 (2019).

Article CAS PubMed PubMed Central Google Scholar

Creyghton, M. P. et al. Histone H3K27ac separates active from poised enhancers and predicts developmental state. Proc. Natl Acad. Sci. USA 107, 21931–21936 (2010).

Article ADS CAS PubMed PubMed Central Google Scholar

Rada-Iglesias, A. et al. A unique chromatin signature uncovers early developmental enhancers in humans. Nature 470, 279–283 (2011).

Article ADS CAS PubMed Google Scholar

Balsalobre, A. & Drouin, J. Pioneer factors as master regulators of the epigenome and cell fate. Nat. Rev. 23, 449–464 (2022).

Article CAS Google Scholar

Gao, N. et al. Dynamic regulation of Pdx1 enhancers by Foxa1 and Foxa2 is essential for pancreas development. Genes Dev. 22, 3435–3448 (2008).

Article CAS PubMed PubMed Central Google Scholar

Bastidas-Ponce, A. et al. Foxa2 and Pdx1 cooperatively regulate postnatal maturation of pancreatic beta-cells. Mol. Metab. 6, 524–534 (2017).

Article CAS PubMed PubMed Central Google Scholar

Gao, N. et al. Foxa1 and Foxa2 maintain the metabolic and secretory features of the mature beta-cell. Mol. Endocrinol. 24, 1594–1604 (2010).

Jia, S. et al. Insm1 cooperates with Neurod1 and Foxa2 to maintain mature pancreatic beta-cell function. EMBO J. 34, 1417–1433 (2015).

Article MathSciNet CAS PubMed PubMed Central Google Scholar

Haliyur, R. et al. Human islets expressing HNF1A variant have defective beta cell transcriptional regulatory networks. J. Clin. Investig. 129, 246–251 (2019).

Article PubMed Google Scholar

DeForest, N. et al. Human gain-of-function variants in HNF1A confer protection from diabetes but independently increase hepatic secretion of atherogenic lipoproteins. Cell Genom. 3, 100339 (2023).

Article CAS PubMed PubMed Central Google Scholar

Kulzer, J. R. et al. A common functional regulatory variant at a type 2 diabetes locus upregulates ARAP1 expression in the pancreatic beta cell. Am. J. Hum. Genet. 94, 186–197 (2014).

Article CAS PubMed PubMed Central Google Scholar

van de Bunt, M. et al. Transcript expression data from human islets links regulatory signals from genome-wide association studies for Type 2 diabetes and glycemic traits to their downstream effectors. PLoS Genet. 11, e1005694 (2015).

Article PubMed PubMed Central Google Scholar

Browaeys, R., Saelens, W. & Saeys, Y. NicheNet: modeling intercellular communication by linking ligands to target genes. Nat. Methods 17, 159–162 (2020).

Article CAS PubMed Google Scholar

Park, S. J. et al. Discovery of endoplasmic reticulum calcium stabilizers to rescue ER-stressed podocytes in nephrotic syndrome. Proc. Natl Acad. Sci. USA 116, 14154–14163 (2019).

Article ADS CAS PubMed PubMed Central Google Scholar

Aloe, L., Rocco, M. L., Balzamino, B. O. & Micera, A. Nerve growth factor: a focus on neuroscience and therapy. Curr. Neuropharmacol. 13, 294–303 (2015).

Article CAS PubMed PubMed Central Google Scholar

Houtz, J., Borden, P., Ceasrine, A., Minichiello, L. & Kuruvilla, R. Neurotrophin signaling is required for glucose-induced insulin secretion. Dev. Cell 39, 329–345 (2016).

Article CAS PubMed PubMed Central Google Scholar

Pingitore, A. et al. Fine tuning of insulin secretion by release of nerve growth factor from mouse and human islet beta-cells. Mol. Cell. Endocrinol. 436, 23–32 (2016).

Article CAS PubMed Google Scholar

Samario-Roman, J. et al. NGF and its role in immunoendocrine communication during metabolic syndrome. Int. J. Mol. Sci. 24, https://doi.org/10.3390/ijms24031957 (2023).

Shimoke, K. et al. NGF-induced phosphatidylinositol 3-kinase signaling pathway prevents thapsigargin-triggered ER stress-mediated apoptosis in PC12 cells. Neurosci. Lett. 389, 124–128 (2005).

Article CAS PubMed Google Scholar

DiGruccio, M. R. et al. Comprehensive alpha, beta and delta cell transcriptomes reveal that ghrelin selectively activates delta cells and promotes somatostatin release from pancreatic islets. Mol. Metab. 5, 449–458 (2016).

Article CAS PubMed PubMed Central Google Scholar

Frodin, M. et al. Glucose, other secretagogues, and nerve growth factor stimulate mitogen-activated protein kinase in the insulin-secreting beta-cell line, INS-1. J. Biol. Chem. 270, 7882–7889 (1995).

Article CAS PubMed Google Scholar

Wang, J. et al. A mutation in the insulin 2 gene induces diabetes with severe pancreatic beta-cell dysfunction in the Mody mouse. J. Clin. Investig. 103, 27–37 (1999).

Article CAS PubMed PubMed Central Google Scholar

Oyadomari, S. et al. Targeted disruption of the Chop gene delays endoplasmic reticulum stress-mediated diabetes. J. Clin. Investig. 109, 525–532 (2002).

Article CAS PubMed PubMed Central Google Scholar

Tomlinson, R. E. et al. NGF-TrkA signaling in sensory nerves is required for skeletal adaptation to mechanical loads in mice. Proc. Natl Acad. Sci. USA 114, E3632–E3641 (2017).

Article CAS PubMed PubMed Central Google Scholar

Xu, D. et al. Efficient delivery of nerve growth factors to the central nervous system for neural regeneration. Adv. Mater. 31, e1900727 (2019).

Zhu, H. et al. Understanding cell fate acquisition in stem-cell-derived pancreatic islets using single-cell multiome-inferred regulomes. Dev. Cell 58, 727–743.e711 (2023).

Article CAS PubMed PubMed Central Google Scholar

Weng, C. et al. Single-cell lineage analysis reveals extensive multimodal transcriptional control during directed beta-cell differentiation. Nat. Metab. 2, 1443–1458 (2020).

Article CAS PubMed PubMed Central Google Scholar

Kubo, N. et al. H3K4me1 facilitates promoter-enhancer interactions and gene activation during embryonic stem cell differentiation. Mol. Cell 84, 1742–1752.e1745 (2024).

Article CAS PubMed Google Scholar

Kang, Y., Kim, Y. W., Kang, J. & Kim, A. Histone H3K4me1 and H3K27ac play roles in nucleosome eviction and eRNA transcription, respectively, at enhancers. FASEB J. 35, e21781 (2021).

Article CAS PubMed Google Scholar

Lee, K. et al. FOXA2 is required for enhancer priming during pancreatic differentiation. Cell Rep. 28, 382–393.e387 (2019).

Article CAS PubMed PubMed Central Google Scholar

Millan-Zambrano, G., Burton, A., Bannister, A. J. & Schneider, R. Histone post-translational modifications - cause and consequence of genome function. Nat. Rev. 23, 563–580 (2022).

Article CAS Google Scholar

Cao, K. et al. An Mll4/COMPASS-Lsd1 epigenetic axis governs enhancer function and pluripotency transition in embryonic stem cells. Sci. Adv. 4, eaap8747 (2018).

Article ADS PubMed PubMed Central Google Scholar

Blandino-Rosano, M. et al. Raptor levels are critical for beta-cell adaptation to a high-fat diet in male mice. Mol. Metab. 75, 101769 (2023).

Article CAS PubMed PubMed Central Google Scholar

Arunagiri, A. et al. Proinsulin misfolding is an early event in the progression to type 2 diabetes. eLife 8, https://doi.org/10.7554/eLife.44532 (2019).

Back, S. H. & Kaufman, R. J. Endoplasmic reticulum stress and type 2 diabetes. Annu. Rev. Biochem. 81, 767–793 (2012).

Article CAS PubMed PubMed Central Google Scholar

Engin, F., Nguyen, T., Yermalovich, A. & Hotamisligil, G. S. Aberrant islet unfolded protein response in type 2 diabetes. Sci. Rep. 4, 4054 (2014).

Article PubMed PubMed Central Google Scholar

Ozcan, U. et al. Endoplasmic reticulum stress links obesity, insulin action, and type 2 diabetes. Science 306, 457–461 (2004).

Huising, M. O. Paracrine regulation of insulin secretion. Diabetologia 63, 2057–2063 (2020).

Article CAS PubMed PubMed Central Google Scholar

Yoshihara, E. et al. Disruption of TBP-2 ameliorates insulin sensitivity and secretion without affecting obesity. Nat. Commun. 1, 127 (2010).

Article ADS PubMed Google Scholar

Zhang, Y. et al. Single-cell epigenome analysis reveals age-associated decay of heterochromatin domains in excitatory neurons in the mouse brain. Cell Res. 32, 1008–1021 (2022).

Article CAS PubMed PubMed Central Google Scholar

Shatzkes, K., Teferedegne, B. & Murata, H. A simple, inexpensive method for preparing cell lysates suitable for downstream reverse transcription quantitative PCR. Sci. Rep. 4, 4659 (2014).

Article PubMed PubMed Central Google Scholar

Langmead, B., Trapnell, C., Pop, M. & Salzberg, S. L. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol. 10, R25 (2009).

Article PubMed PubMed Central Google Scholar

Dobin, A. et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics 29, 15–21 (2013).

Article CAS PubMed Google Scholar

Langmead, B. & Salzberg, S. L. Fast gapped-read alignment with Bowtie 2. Nat. Methods 9, 357–359 (2012).

Article CAS PubMed PubMed Central Google Scholar

Yu, W., Uzun, Y., Zhu, Q., Chen, C. & Tan, K. scATAC-pro: a comprehensive workbench for single-cell chromatin accessibility sequencing data. Genome Biol. 21, 94 (2020).

Article CAS PubMed PubMed Central Google Scholar

Stuart, T. et al. Comprehensive integration of single-cell data. Cell 177, 1888–1902.e1821 (2019).

Article CAS PubMed PubMed Central Google Scholar

McGinnis, C. S., Murrow, L. M. & Gartner, Z. J. DoubletFinder: Doublet detection in single-cell RNA sequencing data using artificial nearest neighbors. Cell Syst. 8, 329–337.e324 (2019).

Article CAS PubMed PubMed Central Google Scholar

Alquicira-Hernandez, J. & Powell, J. E. Nebulosa recovers single-cell gene expression signals by kernel density estimation. Bioinformatics 37, 2485–2487 (2021).

Article CAS PubMed Google Scholar

La Manno, G. et al. RNA velocity of single cells. Nature 560, 494–498 (2018).

Article ADS PubMed PubMed Central Google Scholar

Wolf, F. A. et al. PAGA: graph abstraction reconciles clustering with trajectory inference through a topology preserving map of single cells. Genome Biol. 20, 59 (2019).

Article PubMed PubMed Central Google Scholar

Finak, G. et al. MAST: a flexible statistical framework for assessing transcriptional changes and characterizing heterogeneity in single-cell RNA sequencing data. Genome Biol. 16, 278 (2015).

Article PubMed PubMed Central Google Scholar

Chen, E. Y. et al. Enrichr: interactive and collaborative HTML5 gene list enrichment analysis tool. BMC Bioinforma. 14, 128 (2013).

Article Google Scholar

Kuleshov, M. V. et al. Enrichr: a comprehensive gene set enrichment analysis web server 2016 update. Nucleic Acids Res. 44, W90–W97 (2016).

Article CAS PubMed PubMed Central Google Scholar

Heinz, S. et al. Simple combinations of lineage-determining transcription factors prime cis-regulatory elements required for macrophage and B cell identities. Mol. Cell 38, 576–589 (2010).

Article CAS PubMed PubMed Central Google Scholar

Buniello, A. et al. The NHGRI-EBI GWAS Catalog of published genome-wide association studies, targeted arrays and summary statistics 2019. Nucleic Acids Res. 47, D1005–D1012 (2019).

Article CAS PubMed Google Scholar

Magno, R. & Maia, A. T. gwasrapidd: an R package to query, download and wrangle GWAS catalog data. Bioinformatics 36, 649–650 (2020).

Article CAS PubMed Google Scholar

Love, M. I., Huber, W. & Anders, S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 15, 550 (2014).

Article PubMed PubMed Central Google Scholar

Gu, Z., Eils, R. & Schlesner, M. Complex heatmaps reveal patterns and correlations in multidimensional genomic data. Bioinformatics 32, 2847–2849 (2016).

Article CAS PubMed Google Scholar

Xie, Z. et al. Gene set knowledge discovery with Enrichr. Curr. Protoc. 1, e90 (2021).

Article CAS PubMed PubMed Central Google Scholar

Xin, Y. et al. Pseudotime ordering of single human beta-cells reveals states of insulin production and unfolded protein response. Diabetes 67, 1783–1794 (2018).

Article CAS PubMed Google Scholar

Download references

This work was supported by funding from the NIH (R01 DK132651, K01 DK120808, UL1 TR002377), Roubos Family Fund in Research, and the Mayo NORP program to Z.W.

These authors contributed equally: Liu Wang, Jie Wu.

Department of Physiology and Biomedical Engineering, Mayo Clinic College of Medicine, Scottsdale, AZ, USA

Liu Wang, Jie Wu, Madeline Sramek, S. M. Bukola Obayomi & Zong Wei

Department of Genetics and Genome Sciences, Case Western Reserve University, Cleveland, OH, USA

Peidong Gao & Yan Li

Department of Physiology and Biomedical Engineering and Division of Endocrinology, Diabetes and Metabolism, Mayo Clinic College of Medicine, Rochester, MN, USA

Aleksey V. Matveyenko

Division of Endocrinology, Mayo Clinic, Scottsdale, AZ, USA

Zong Wei

You can also search for this author in PubMed Google Scholar

You can also search for this author in PubMed Google Scholar

You can also search for this author in PubMed Google Scholar

You can also search for this author in PubMed Google Scholar

You can also search for this author in PubMed Google Scholar

You can also search for this author in PubMed Google Scholar

You can also search for this author in PubMed Google Scholar

You can also search for this author in PubMed Google Scholar

L.W., J.W., and Z.W. designed the study. J.W. performed the Paired-Tag, in vitro, and mouse studies. L.W. performed all the computational analysis. M.S. performed immunostaining. S.B.O. and P.G. provided technical assistance. Y.L. and A.M. advised on the manuscript. L.W., J.W., and Z.W. wrote the manuscript with input from all authors.

Correspondence to Zong Wei.

The authors declare no competing interests.

Nature Communications thanks David Hodson, Annette Schürmann, and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. A peer review file is available.

Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Open Access This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/.

Reprints and permissions

Wang, L., Wu, J., Sramek, M. et al. Heterogeneous enhancer states orchestrate β cell responses to metabolic stress. Nat Commun 15, 9361 (2024). https://doi.org/10.1038/s41467-024-53717-0

Download citation

Received: 23 October 2023

Accepted: 18 October 2024

Published: 30 October 2024

DOI: https://doi.org/10.1038/s41467-024-53717-0

Anyone you share the following link with will be able to read this content:

Sorry, a shareable link is not currently available for this article.

Provided by the Springer Nature SharedIt content-sharing initiative