Single-cell multiome and spatial profiling reveals pancreas cell type-specific gene regulatory programs of type 1 diabetes progression

Abstract

Cell type-specific regulatory programs that drive type 1 diabetes (T1D) in the pancreas are poorly understood. Here, we performed single-nucleus multiomics and spatial transcriptomics in up to 32 nondiabetic (ND), autoantibody-positive (AAB⁺), and T1D pancreas donors. Genomic profiles from 853,005 cells mapped to 12 pancreatic cell types, including multiple exocrine subtypes. β, Acinar, and other cell types, and related cellular niches, had altered abundance and gene activity in T1D progression, including distinct pathways altered in AAB⁺ compared to T1D. We identified epigenomic drivers of gene activity in T1D and AAB⁺ which, combined with genetic association, revealed causal pathways of T1D risk including antigen presentation in β cells. Last, single-cell and spatial profiles together revealed widespread changes in cell-cell signaling in T1D including signals affecting β cell regulation. Overall, these results revealed drivers of T1D in the pancreas, which form the basis for therapeutic targets for disease prevention.

Fig. 1.. Cell type–specific map of gene expression in the pancreas.

(A) Design of study profiling human pancreas from ND, ND AAB⁺, and T1D donors using single-cell assays. (B) Uniform Manifold Approximation and Projection (UMAP) plot showing clustering of 276,906 nuclei from single-nuclear RNA-seq of 32 whole pancreas donors from the nPOD biorepository. Clusters are labeled on the basis of cell type and subtype annotations. (C) Dot plot showing the normalized expression levels of selected known marker genes for pancreatic cell types and subtypes. (D) Dot plot of genes with preferential expression across different subtypes of acinar cells (top left) and normalized enrichment score (NES) of pathways enriched in each subtype using fgsea (top right). Box plot showing donor CPM of selected genes with preferential expression in different subtypes of acinar cells. GTPase, guanosine triphosphatase. (E) Representative FOV per condition (ND: top, T1D: bottom) showing (from left to right) immunofluorescence, coarse cell type annotation with the spatial gene panel directly, and finer-grained cell type annotation transferred from the snRNA-seq data. DAPI, 4′,6-diamidino-2-phenylindole. (F) Matrix plots showing the neighborhood enrichment of cell types based on spatial neighbors. Cell type labels are the same as fine-grained annotations in (E). (G) Stacked barplot illustrating the relative abundance of each cell type in each multicellular niche. Cell type labels are the same as fine-grained annotations in (E, left). Dot plot showing the normalized gene expression levels of spatially variable genes across multicellular pancreatic niches (right). (H) Box plot showing normalized cell counts for selected pancreatic cell types and subtypes grouped by T1D status. P values from likelihood ratio test, **FDR < 0.1, *uncorrected P < 0.05. (I) Stacked barplot showing the relative abundance of each multicellular niche per condition. Niches have altered abundance in ND samples denoted as *P < 0.05.

Fig. 2.. Cell type–specific map of accessible chromatin in the pancreas.

(A) UMAP plot showing clustering of 203,348 nuclei from single-nuclear ATAC-seq of 30 pancreas donors from the nPOD biorepository. Clusters are labeled with cell type and subtype identity based on label transfer from the gene expression map. (B) Genome browser showing accessible chromatin signal at the promoter regions of known marker genes for pancreatic cell types. (C) Heatmap showing genome-wide accessibility from chromVAR of sequence motifs for selected transcription factors (TFs) across cell types (left) and box plots showing donor-level accessibility of selected TF motifs across cell types (right). (D) Box plot showing genome-wide accessibility of TF motifs with preferential enrichment in subtypes of acinar cells (left) and log fold change in expression for genes in structural subfamilies for enriched TF motifs, and error bars are SE (right). *FDR < 0.10. (E) Number of cREs identified across all cell types and the percentage of cREs that do not overlap previous catalogs of cREs (27, 53) (top). Example of a pancreatic T cell–specific cRE at the ZNF746 locus. (F) TF sequence motifs enriched in cREs with activity specific to each cell type (left) and bar plots showing −log₁₀ P values of gene sets enriched for proximity to cell type–specific cREs using the Genomic Regions Enrichment of Annotations Tool (GREAT). (G) Example of a cRE active in pancreatic T cells and macrophages that overlaps a candidate T1D risk variant. (H) Box plot showing gene-cRE links per gene per cell type (top) and schematic of TF GRNs (bottom). (I) Matrix plot showing scaled z-score of TF activities for top TFs identified for each cell type using a t test with overestimated variance. (J) Spatial plot of selected TFs showing the TF activity profile (top) and cell type distribution for the respective cell type (bottom).

Fig. 3.. Cell type–specific changes in gene expression in T1D progression.

(A) Number of genes in each cell type with significant (FDR < 0.1) changes in expression in T1D stages compared to nondiabetes (top). Number of pathways enriched in genes with up- and down-regulated expression in each cell type in T1D stages (bottom). The results in all panels include 80 donors (nPOD + HPAP) for endocrine cells and 32 donors (nPOD) for nonendocrine cells. Note that nonendocrine cells were not tested for single ND AAB⁺ association. (B) Volcano plot showing differential expression in β cells in recent-onset T1D compared to ND. (C) Bar plot showing normalized enrichment (NES) of pathways enriched in up- and down-regulated genes in β cells in recent-onset T1D (bottom). MT, Mitochondrial. (D) Scaled expression in spatial profiles of genes with up-regulated expression in T1D in β cells (left). Spatially dependent expression of selected genes up-regulated in T1D in each cell type (right). (E) Pathways with differential expression within spatial niches in T1D compared to ND. EGFR, epidermal growth factor receptor; TNFα, tumor necrosis factor–α. (F) Scatterplot of log fold change in expression of genes in β cells in single or multiple ND AAB⁺ compared to recent-onset (top) and long-duration T1D (bottom). Line shown in each plot is from a linear model of log fold change values, and P values are from Spearman correlation. FC, fold change. (G) Normalized enrichment of pathways in recent-onset T1D and multiple ND AAB⁺. Pathways are colored on the basis of significant enrichment (FDR < 0.1) in either, or both, states. (H) Normalized enrichment of pathways in β cells across each T1D state compared to nondiabetes. (I) Log fold change in expression of selected genes in β cells in each T1D state compared to ND. (J) Normalized enrichment of pathways in other pancreatic cell types in recent-onset T1D or multiple ND AAB⁺. For all panels, *FDR < 0.1, ^•P < 0.05.

Fig. 4.. Epigenomic changes in pancreatic cell types in T1D progression.

(A) Fold enrichment of sequence motifs for TFs enriched in β cell cREs with up-regulated or down-regulated activity in (top) recent-onset T1D or (bottom) ND AAB⁺ (both single and multiple) using snATAC-seq from 30 donors (nPOD). (B) Box plots showing donor-level genome-wide accessibility of selected TF motifs in β cells (left) and α cells (right) from chromVAR across nondiabetes (ND), ND AAB⁺, and recent-onset T1D (T1D). (C) TF GRNs enriched for overlap with genes in biological pathways in β cells altered in T1D progression. (D) Biological pathways in β cells enriched for overlap with the HNF1A GRN (top). β Cell expression of HNF1A in T1D progression, values represent log fold change and error bars are SE (middle). β Cell activity of biological pathways linked to the HNF1A GRN in T1D progression (bottom). (E) TF GRNs enriched for overlap with genes in biological pathways in acinar cells altered in T1D progression. (F) Genome browser views of the TSHR (top) and HLA-A (bottom) loci where β cell cREs with altered activity in T1D were linked to genes with concordant changes in expression in T1D. (G) Genome-wide enrichment of T1D-associated variants in β cell cREs linked to pathways with altered expression in ND AAB⁺. Values represent log enrichment from fgwas, and error bars are 95% confidence interval. (H) Genome browser view of T1D-associated variants and β cell accessible chromatin in ND, ND AAB⁺, and T1D at the IRF1 locus, where candidate T1D variant overlaps a β cell cRE with altered activity in T1D progression. (I) Genome browser view of T1D-associated variants and both β cell and T cell accessible chromatin in ND, ND AAB⁺, and T1D at the STAT4 locus. Candidate T1D variants at this locus overlap T cell cREs.

Fig. 5.. Cell-cell signaling networks altered in T1D progression.

(A) Summary of total interaction strength (top) and number of interactions (middle) for each pancreatic lineage in nondiabetes, ND AAB⁺ (both single and multiple), recent-onset T1D, and long-standing T1D using snRNA-seq from 32 donors (nPOD). Bar plot showing the number of LR interactions per donor FOV in spatial slides, and error bars are SE (bottom). (B) Heatmap showing normalized interaction strength of signals for each cell type among donors which were nondiabetes, ND AAB⁺, recent-onset T1D, and long-standing T1D. Stars represent significance of the difference in interaction strength in each T1D state compared to nondiabetes. **FDR < 0.01 and *FDR < 0.05. (C) Difference in strength of interactions between β cells and other cell types and subtypes in ND AAB⁺, recent-onset T1D, and long-duration T1D. **FDR < 0.01 and *FDR < 0.05. (D) Interaction strength of signals for each cell type summarized by functional categories. **FDR < 0.01 and *FDR<0.05. CAM, cell adhesion molecule. (E) Normalized interaction strength in recent-onset T1D and nondiabetes for ligands with significant change in signal involving β cells. **FDR < 0.01 and *FDR < 0.05. (F) Heatmap per donor showing the interaction score of the top LR interactions from a likelihood ratio test comparing ND and T1D donors. (G) Spatial plots of a representative FOV (T1D: top, ND: bottom) highlighting spots with an interaction between HLA-C and CD8A and the cell types where this interaction occurs. (H) Volcano plot showing genes with up- and down-regulated expression in EndoC-BH1 after BMP5 treatment (left). Pathways enriched in up- and down-regulated genes in BMP5 exposure (right). The experiment was performed using n = 6 biological replicates per treatment. MAPK, mitogen-activated protein kinase. (I) Volcano plot showing genes with up- and down-regulated expression in EndoC-BH1 after progranulin (PGRN) treatment (left). Biological pathways enriched in genes with up- and down-regulated expression after PGRN exposure (right). The experiment was performed using n = 3 biological replicates per treatment. TCA, citric acid cycle.