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.
     __________________________________________________________________

   Pancreatic cell types show changes in gene regulation and activity
   associated with T1D progression and risk.

INTRODUCTION

   Type 1 diabetes (T1D) is a complex endocrine disorder characterized by
   autoimmune destruction of β cells in the pancreatic islets, leading to
   lifelong dependence on insulin therapy. The destruction of β cells in
   T1D is caused by interactions with multiple cell types in and
   surrounding the islet microenvironment including infiltrating immune
   cells, other endocrine cells, and endothelial cells ([69]1–[70]3). Cell
   types in the pancreas outside the local islet environment, such as
   exocrine acinar and ductal cells, are also increasingly implicated in
   T1D pathogenesis ([71]4, [72]5). β Cells themselves likely contribute
   to the development of T1D as well through response to environmental
   factors, external signaling to immune, β, and other cell types, and
   cellular survival ([73]6). The sequence of events in the pancreas that
   drives initiation of β cell autoimmunity and progression through stages
   of T1D, however, as well as the role of each pancreatic cell type in
   these processes, remains poorly understood.

   Seroconversion to autoantibody positivity (AAB^+) against islet
   proteins (i.e., self-antigens) precedes T1D onset in nearly all cases
   and is used as a clinical biomarker of T1D progression ([74]7, [75]8).
   Individuals at T1D diagnosis can present with a differing number and
   type of autoantibodies, which are associated with varying rates of
   disease incidence; for example, the presence of a single islet AAB has
   a relatively low lifetime risk of T1D, whereas individuals with
   multiple AAB have disease rates more than 90% ([76]9–[77]11). As
   clinical presentation of T1D does not occur until a large fraction of β
   cells has been destroyed, there is a window of time between
   seroconversion and T1D onset where disease processes can potentially be
   halted or reversed ([78]7). Even after onset of T1D, residual β cell
   mass could potentially be modulated therapeutically to restore
   metabolic function ([79]12). Defining changes in disease-relevant cell
   types across the stages of T1D progression would both improve our
   understanding of the mechanisms of T1D as well as reveal potential
   targets to prevent or reverse disease. Furthermore, an improved
   understanding of key changes associated with progression would also
   help identify biomarkers of T1D, which are particularly needed in the
   early stages of disease to identify progressors and candidates for
   therapeutic intervention ([80]13).

   Single-cell technology, the focus of this work, enables profiling of
   individual cells within the pancreas ([81]5, [82]14). Previous
   single-cell studies of the pancreas in T1D have been limited in that
   they focused primarily on gene expression profiling of dispersed cells
   ([83]4, [84]15), which does not provide information on the spatial
   localization of cellular transcriptomes within the pancreas nor the
   genomic elements driving changes in gene activity. Recent developments
   in spatial transcriptomics enables profiling cells in their native
   location ([85]16), which enables understanding cell type–specific
   changes in the context of specific cellular neighborhoods and niches in
   the pancreas. This is particularly important in the context of T1D
   which has extensive heterogeneity in disease processes within the
   pancreas ([86]17). In addition, single-cell epigenome profiling, for
   example, using single-nucleus Assay for Transposase-Accessible
   Chromatin using Sequencing (snATAC-seq) or single-cell multiome [paired
   single-nucleus RNA sequencing (snRNA-seq) + snATAC-seq], can reveal
   transcriptional regulators, cis-regulatory elements (cREs), and gene
   regulatory networks (GRNs) driving altered gene expression in T1D
   ([87]14, [88]15). Critically, GRNs and cREs can be intersected with
   T1D-associated variation to infer cell type–specific regulatory
   programs that may play a causal role in driving disease ([89]5,
   [90]18).

   Previous single-cell studies have also been limited in the extent to
   which they have captured key windows of T1D progression and
   pathogenesis ([91]4, [92]15). Specifically, nondiabetic (ND) AAB^+
   donors in these efforts were largely those with single glutamic acid
   decarboxylase (GAD) autoantibodies ([93]4), which have a relatively
   lower risk of developing T1D compared to multiple AAB^+ donors and do
   not reflect the full arc of progression to T1D ([94]19). Furthermore,
   many of the T1D donors in these studies had long-standing T1D where
   disease processes are potentially more dormant, whereas profiling
   donors who had more recently developed T1D may give greater insight
   into active disease processes. Third, these studies profiled purified
   islets, and profiling tissue sampled from the whole pancreas may offer
   greater insight into genomic changes during T1D progression in exocrine
   cells and other cell types outside of the islet microenvironment.

   In this study, we addressed these limitations by performing single-cell
   gene expression and epigenome profiling in whole pancreas from 32 ND,
   ND single and multiple AAB^+, recent-onset T1D, and long-standing T1D
   organ donors, as well as spatial transcriptomics in a subset of ND and
   recent-onset T1D donors. We determined changes in pancreatic cell type
   abundance, cellular pathways, GRNs, and cell-cell signaling across
   these stages of T1D progression and pathogenesis and, using T1D
   association data, identified pathways and gene networks that may play a
   causal role in the development of T1D.

RESULTS

A comprehensive, multimodal, spatially resolved map of pancreatic cell types

   We obtained pancreatic tissue from 32 donors in the Network for
   Pancreatic Organ Donors with Diabetes (nPOD) biorepository including 11
   ND, 9 ND autoantibody positive (ND AAB^+), and 12 T1D which we
   separated into 7 recent onset (<1 year from diagnosis) and 5 longer
   duration (>5 years from diagnosis) (table S1). Within the ND AAB^+
   group, most organ donors, by our study design, had multiple
   autoantibodies (multiple ND AAB^+). For all samples, we performed
   snRNA-seq and snATAC-seq assays, and for eight of the samples, we
   performed single-nucleus multiome (joint snRNA-seq and snATAC-seq in
   the same nucleus) assays (table S1). In addition, for six of the
   samples, we performed spatial transcriptomic assays using the CosMx
   Spatial Molecular Imager ([95]Fig. 1A).

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

   [96]Fig. 1.
   [97]Open in a new tab

   (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.

   After extensive barcode quality control and filtering steps (see
   Materials and Methods), we integrated data from all 32 donors using
   Harmony ([98]20) and clustered 276,906 gene expression profiles
   ([99]Fig. 1B and figs. S1 and S2). We annotated the resulting 18
   clusters based on the expression of known cell type marker genes which
   revealed 12 pancreatic cell types including exocrine (acinar and
   ductal), endocrine (α, β, and δ), immune (T cell, B cell, macrophage,
   and mast), stellate, endothelial, and Schwann cells ([100]Fig. 1, B and
   C, and table S2). Cell type clusters had broadly consistent
   representation across donors and donor characteristics (figs. S2 and
   S3). We aggregated expression profiles for all cells in the cell type
   and derived normalized expression levels of each gene using counts per
   million (CPM) (data S1). For each cell type, we further identified
   genes with expression levels specific [false discovery rate
   (FDR) < 0.1] to the cell type which revealed both known and previously
   unreported sets of genes with cell type–specific expression (table S3);
   for example, well-known genes with expression specific to β cells
   included INS, IAPP, and G6PC2 as well as others with no now known role
   in β cell function (e.g., PLCH2, NRG2, RBFOX3, and MTUS2).

   Several cell types displayed multiple subclusters representing both
   known cell subtypes, such as active and quiescent stellate cells, blood
   vessel cells (BVECs), lymphatic endothelial cells (LECs), and
   MUC5b^+ ductal cells, as well as several potential subtypes of acinar
   cells ([101]Fig. 1, B and C). As the genomic properties of these
   subtypes have not been completely described previously, we derived sets
   of marker genes for each subtype (see Materials and Methods and table
   S3). For BVECs and LECs, in addition to reported marker genes PLVAP
   (BVECs) and FLT4 (LECs), we observed specific up-regulation of genes in
   each subtype such as INHBB, BMP6, FCN3, and PCAT19 in BVECs and EFNA5,
   COLEC12, and MYCT1 in LECs. In MUCB5^+ ductal cells, there was
   up-regulation of ERN2, CYP2C18, MYO7B, and DMBT1 compared to the
   primary subtype of ductal cells. For acinar cells, the primary cluster,
   which we annotated as “basal” acinar cells, was enriched for genes and
   pathways involved in digestive enzyme production and secretion. Other
   clusters included “high-enzyme” acinar cells with higher expression of
   enzymes such as chymotrypsin (CTRB1/2), trypsinogen (PRSS1 and PRSS2),
   lipase (PNLIP), carboxyl ester lipase (CEL), chymotrypsin-like elastase
   (CELA3A/B), and increased oxidative phosphorylation and translation,
   “signaling” acinar cells with increased signaling and stress-response
   activity, and “signaling/differentiation” acinar cells with increased
   signaling, metallothionein (MT1/MT2), and identity and differentiation
   genes (REG1A/B and PTF1A) ([102]Fig. 1D).

   To next characterize the spatial organization of pancreatic cell types,
   we performed RNA in situ hybridization of 1010 genes with CosMx from a
   subset of donors including three ND and three recent-onset T1D (tables
   S1 and S4). We imaged a total of 82.6M transcripts from 71 fields of
   view (FOVs) in pancreatic sections from three ND (32 FOVs) and three
   recent-onset T1D donors (39 FOVs) (fig. S4A) and assigned transcripts
   to 392,248 cells overall (80 median genes and 200 median transcripts
   per cell), using the CosMx default segmentation. We performed
   unsupervised clustering of cellular gene expression profiles, which
   revealed nine distinct clusters including exocrine, endocrine,
   endothelial, immune, and mast cells (fig. S4B). We next mapped
   finer-grained cell type annotations from snRNA-seq using moscot (fig.
   S4, B and C) ([103]21), which revealed 14 cell types and subtypes that
   were confirmed on the basis of marker gene expression ([104]Fig. 1E and
   fig. S4). Spatial neighborhood enrichment using squidpy ([105]22)
   revealed expected cell type clustering including acinar subtypes,
   ductal subtypes, endocrine cells (β, α, and δ), and connective cells
   (e.g., endothelial, immune, and stellate) ([106]Fig. 1F).

   Next, we sought to determine whether spatial neighborhoods form
   recurrent niches across the pancreas, by defining niches involving a
   cell type using a gene-gene covariance matrix ([107]23) in a spatial
   neighborhood of 30 cells. We recovered six niches in total,
   characterized by cell type abundance (see Materials and Methods),
   including three exocrine (acinar-enriched, ductal-enriched, and MUC5b
   ductal-enriched) niches, one endocrine niche, one niche including both
   endocrine and exocrine cells (endo-exo), and one niche consisting of
   connective cells ([108]Fig. 1G). To characterize each niche, we
   identified spatially variable genes (Moran’s I > 0.2, P < 0.05) that
   captured gene signatures specific to the niche ([109]Fig. 1G). In the
   acinar-enriched niche, marker genes from the basal and high-enzyme cell
   types showed strong spatial clustering (PRSS2 and REG1A). In
   comparison, the ductal-enriched niche had more spatial association with
   signaling and signaling/differentiation acinar cells (MT1X, SOD2, and
   MT2A). In the MUC5b ductal–enriched niche, spatially variable genes
   were strongly associated with immune interactions (HSPA1A, HLA-A, and
   B2M). In addition, the endocrine niche had highly distinct patterns
   which highlighted multiple endocrine-specific genes (e.g., INS, GCG,
   SST, and IAPP) ([110]Fig. 1G).

   Last, we determined whether there were changes in abundance of cell
   types and subtypes in T1D progression based on snRNA-seq (see Materials
   and Methods). There was a significant decrease (likelihood ratio test,
   FDR < 0.1) in β cells ([111]Fig. 1H and table S5) although we still
   observed residual β cells in T1D particularly in recent-onset
   (ND = 1.5%, recent-onset T1D = 0.93%). We also observed a significant
   decrease (FDR < 0.1) in δ cells in T1D and increased abundance of
   immune populations in ND AAB^+ and recent-onset T1D. There was also
   nominal evidence (P < 0.05) for altered abundance of specific subtypes
   including high-enzyme acinar (P = 0.037) and MUC5b^+ ductal cells
   (P = 0.049). We next asked whether there were corresponding changes in
   the abundance of specific niches in T1D in spatial profiles. We
   quantified the pairwise similarity between ND and T1D spatial graphs
   using Wasserstein distance (fig. S4D) ([112]24), which revealed
   substantial changes in the underlying structure of endocrine cells (α
   and β) in T1D. We also observed significant changes in the abundance of
   the endocrine and MUC5b^+ ductal niche in T1D (P < 0.05) ([113]Fig.
   1I).

Comprehensive map of pancreatic cell type–accessible chromatin

   To understand how the epigenome may drive changes in cell type–specific
   gene expression in T1D, we next created a matched map of accessible
   chromatin in pancreatic cell types. Of the 32 nPOD donors with
   snATAC-seq assays, 30 passed quality control (see Materials and
   Methods). We filtered, integrated, and clustered accessible chromatin
   profiles from the 30 donors and annotated cell type identity by label
   transfer of the gene expression map using Seurat (see Materials and
   Methods) ([114]25). After filtering nuclei with low transfer
   predictions (<0.5), there were 203,348 chromatin profiles mapping to
   the same cell types and subtypes ([115]Fig. 2A and figs. S6 and S7). We
   estimated that label transfer was >97% accurate at the cell type level
   by comparing the predicted and actual identity of accessible chromatin
   profiles in single-cell multiome data. We also confirmed that predicted
   cell types were accessible at the promoter regions of key marker genes
   ([116]Fig. 2B). The proportions of each cell type were highly
   correlated between expression and chromatin maps (r = 0.98, P = 1.7 ×
   10^−13; fig. S8).

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

   [117]Fig. 2.
   [118]Open in a new tab

   (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 ([119]27, [120]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[10] 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).

   We identified transcription factor (TF) binding motifs preferentially
   enriched in each pancreatic cell type and subtype using chromVAR
   ([121]26). At the cell type level, enriched sequence motifs revealed
   key regulators of each cell type; for example, β cells and other
   endocrine cells were enriched for RFX and FOXA motifs; ductal cells for
   HNF1, ONECUT, and TEAD motifs; endothelial cells for ETV, FLI, and
   GABPA motifs; and T cells for RUNX, ETV, and ETS motifs, among others
   ([122]Fig. 2C and table S6). Motif enrichments also highlighted
   regulators that distinguished accessible chromatin profiles of cell
   types within specific lineages; for example, NEUROD1 and NR3C1 had
   stronger enrichment in β compared to other endocrine cells ([123]Fig.
   2C). Acinar cells showed distinct sets of enriched TF motifs across
   different subclusters, including signaling acinar cells which were more
   enriched for FOS/JUN, activating transcription factor (ATF), and
   Nuclear Factor Erythroid (NFE) motifs ([124]Fig. 2, C and D, and table
   S6). In high-enzyme acinar cells, the strongest enrichments were for
   TFs such as Zinc finger E-box-binding (ZEB), Snail transcription factor
   (SNAI1-3), and Transcription Factor (TCF3-4), which were also the most
   enriched motifs in acinar cells overall compared to other cell types
   ([125]Fig. 2, C and D, and table S6). As structurally related TFs often
   have similar motifs, we linked TF motifs enriched in subclusters to
   specific TFs in the same structural subfamily with concordant
   expression patterns. For example, FOSL2 and JUNB/D, as well as ATF3,
   NFE2L2, and BACH1/2, were increased in signaling acinar cells, and TCF3
   had increased expression in high-enzyme acinar cells ([126]Fig. 2D).

   For each cell type and subtype, we next defined candidate cREs. We
   derived “pseudo”-bulk accessible chromatin profiles by aggregating
   reads from all cells for that cell type or subtype and identified cREs
   by performing peak calling with Model-based Analysis of ChIP-seq
   (MACS2). In total, there were 368,688 cREs across cell types and an
   average of 94.3k cis-regulatory elements (cREs) per cell type (data
   S2). Among cREs in our study, 9.4 and 7.4% were unique compared to a
   pan-tissue ([127]27) and pancreas-specific ([128]5) cRE catalog,
   respectively, such as a T cell cRE directly upstream of ZNF746
   ([129]Fig. 2E). We identified cREs with cell type–specific activity by
   comparing accessible chromatin profiles across cell types (data S3; see
   Materials and Methods). Cell type–specific cREs were enriched for
   sequence motifs of key cell type TFs as well as proximity to genes
   involved in cell type–specific function (tables S7 and S8). For
   example, β cell–specific cREs were significantly enriched (FDR < 0.1)
   for proximity to insulin secretion–related pathways and RFX, FOXA,
   NEUROD, and NKX6.1 TF motifs, whereas endothelial-specific cREs were
   significantly enriched for proximity to angiogenesis, blood vessel
   morphogenesis, and vasculature pathways and FLI, ETS, and ETV TF motifs
   (tables S7 and S8). We also identified cREs specific to several of the
   subtypes within acinar cells; for example, signaling acinar-specific
   cREs were enriched for JUN, FOS, and ATF motifs.

   Because of the scarcity of immune populations in the pancreas, the
   epigenome of resident and infiltrating pancreatic immune cells has not
   been extensively described. In our study, we identified multiple immune
   cell types including T cells, macrophages, B cells, and mast cells,
   although available cell numbers only enabled defining cREs in T cells
   and macrophages. T cell–specific cREs were significantly enriched for
   proximity to genes involved in T cell activation, T cell receptor
   complex, and cytokine receptor activity and motifs for ETS, ETV, and
   RUNX TFs, while macrophage-specific cREs were enriched for
   immune-related processes and PU.1 and SPIB motifs ([130]Fig. 2F).
   Compared to a previous study which profiled several whole pancreas
   donors, more than double the number of cREs were identified in each
   cell type (T cells: 58.8k versus 24.5k; macrophages: 114.3k versus
   55.7k). The increased number of cREs improved annotation of
   T1D-associated variants; for example, candidate T1D variant rs947474
   [posterior probability of association (PPA) = 0.88] ([131]5) overlapped
   a pancreatic T cell and macrophage cRE not identified in these cell
   types in the pancreas previously and not active in other pancreatic
   cell types ([132]Fig. 2G).

   We next predicted networks of genes regulated by TF activity in each
   pancreatic cell type (see Materials and Methods). We linked cREs to
   target genes in each cell type using the activity-by-contact (ABC)
   method, which revealed an average of 46,474 cRE-target gene links per
   cell type, as well as based on promoter proximity (data S4). Using ABC
   and promoter proximity, genes were linked to, on average, 2.8 cREs per
   cell type ([133]Fig. 2H). We identified genes which had highly cell
   type–specific cRE links (see Materials and Methods), and genes with
   highly cell type–specific cRE links included key marker genes such as
   INS in β cells, GCG in α cells, IL-2, IFNGR1, and GZMA in T cells, and
   MARCOS in macrophages. In each cell type, we next constructed GRNs for
   366 TFs by combining (i) cRE-target gene links, (ii) TF motif
   predictions in cREs, and (ii) TF and target gene expression levels
   ([134]Fig. 2H; see Materials and Methods and data S5). We annotated
   likely cellular functions of TF GRNs by identifying biological pathways
   with gene sets that significantly overlapped TF GRNs. There were
   thousands of significant relationships linking TF GRNs to biological
   pathways across all cell types (Fisher’s test, FDR < 0.1) (data S6),
   which annotated many known regulators of pathway activity as well as
   many putative functions of TFs.

   Last, we used spatial transcriptomic data in combination with cell
   type–specific TF GRNs to infer TF activity within cell types and
   subtypes in the pancreas. Briefly, we used a univariate linear model to
   predict the observed gene expression based on TF-gene interaction
   weights, from which we scored TFs as active or inactive in each cell
   type ([135]28). We identified TFs with endocrine-specific activity in
   line with the previously described regulators of endocrine cell
   activity, such as NEUROD1, as well as high activity of PAX6 in β cells,
   where it is a key regulator of β cell identity, function, and survival
   ([136]Fig. 2, I and J) ([137]29). Among other cell types, we inferred
   high activity for BHLHA15/MIST1 in acinar cells, where it may play a
   role in the maintenance of pancreatic acinar identity ([138]30), and
   highly specific activity for MEOX2 in endothelial cells and RUNX3 in T
   cells ([139]Fig. 2, I and J). Integrating GRNs with spatial
   transcriptomic profiles thus confirmed the specificity of key TFs
   regulating pancreatic cell types, including for TFs not measured on the
   spatial panel directly.

Pancreatic cell type gene expression in T1D progression

   Changes in genome-wide gene activity within each pancreatic cell type
   during progression to T1D are poorly understood. We therefore
   identified genes and biological pathways in each cell type with altered
   activity in ND AAB^+ and T1D. To increase our power to detect changes
   in endocrine cell types, we also used single-cell RNA-seq from purified
   islets of 48 donors from the Human Pancreas Analysis Program (HPAP)
   ([140]4, [141]15, [142]31) including 27 ND, 11 ND AAB^+ (9 single, 2
   multiple), and 10 T1D (6 recent, 4 long-standing). For each cell type
   and subtype, we derived gene counts per sample, tested for differential
   expression in single and multiple ND AAB^+ and recent and long-standing
   T1D compared to nondiabetes, and considered genes significant at
   FDR < 0.1 (see Materials and Methods). We further performed gene set
   enrichment of differential expression results for each cell type and
   subtype and identified pathways with significant (FDR < 0.1) changes in
   activity in each condition (see Materials and Methods).

   Marked gene expression changes were observed in β cells in T1D
   ([143]Fig. 3A). In recent-onset T1D, 704 genes in β cells had a
   significant change (FDR < 0.1) in expression, where the most
   up-regulated genes included major histocompatibility complex (MHC)
   class I and related (CD74 and B2M) genes, cytokines, and
   cytokine-induced genes (IL15, GBP2, and IFIT3), cytokine-responsive TFs
   (STAT1/4 and IRF1), and components of the 20S proteosome ([144]Fig. 3B,
   fig. S9, and data S7). We also observed up-regulation of MHC class II
   genes in T1D, particularly HLA-DPB1, as well as MHC class 1 modulators
   such as NLRC5 (fig. S9). At the pathway level, there was up-regulation
   of antigen processing and presentation, interferon signaling,
   interleukin signaling, and Janus kinase–signal transducers and
   activators of transcription (JAK-STAT) signaling and down-regulation of
   oxidative phosphorylation, translation, mitochondrial function,
   mitosis, mRNA processing, protein folding and localization, endoplasmic
   reticulum (ER)–Golgi transport, and autophagy ([145]Fig. 3C and table
   S9). We examined whether specific pathways up-regulated in T1D showed
   heterogeneity in expression across single β cells and pathways
   including antigen presentation had evidence for bimodal expression
   while others such as interferon and JAK-STAT signaling did not (fig.
   S10). Further, a higher proportion of β cells from ND AAB^+ and T1D
   donors expressed antigen presentation pathways compared to ND (fig.
   S10).

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

   [146]Fig. 3.
   [147]Open in a new tab

   (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.

   Compared to recent-onset T1D, the most significant changes in gene
   expression generally differed in long-standing T1D, including
   down-regulation of INS and key genes involved in β cell function such
   as GLIS3 and G6PC2 (data S7). In addition, antigen presentation and
   class I MHC genes were less pronounced, specific IRF TFs had higher
   expression, and class II MHC genes had stronger up-regulation in
   long-standing T1D (fig. S9). At the pathway level, interferon signaling
   was significantly enriched (FDR < 0.1) in long-standing T1D, although
   there were overall few pathways with significant changes in activity
   potentially due to the low number of β cells at this stage (table S9).
   We identified more nominal enrichment for multiple pathways including
   up-regulated antigen presentation and cytokine receptor interactions
   and down-regulated autophagy, insulin processing, β cell regulation,
   and β cell development (table S9).

   Given marked changes in gene expression in β cells in recent-onset T1D,
   we further characterized whether these pathways had altered activity
   within specific localizations in the pancreas. Of the genes with
   altered expression in β cells in recent onset T1D and present in the
   spatial gene panel, almost all (95%) were up-regulated in T1D in
   spatial profiles ([148]Fig. 3D and fig. S11). Furthermore, multiple
   up-regulated genes in T1D such as MHC class I genes (e.g., HLA-A and
   B2M) showed spatially dependent expression patterns (Moran’s I > 0.2)
   within endocrine, immune, and ductal cells ([149]Fig. 3D). We further
   characterized pathways in recent-onset T1D with expression profiles
   dependent on specific niches and altered in T1D progression. We
   identified pathways in the PROGENy database in LIANA+ ([150]32) to
   predict pathways preferentially active in a niche using a multivariate
   linear model. We identified multiple pathways with niche-dependent
   expression, including hypoxia in the endocrine niche ([151]Fig. 3E).
   When further assessing T1D-specific changes, pathways related to
   hypoxia and inflammation such as tumor necrosis factor–α and JAK-STAT
   were differentially active in T1D (fig. S11).

   In contrast to T1D, few individual genes had significant changes in
   expression in β cells in either single or multiple ND AAB^+ ([152]Fig.
   3A). We determined whether more subtle changes might be occurring at
   these stages. Genes altered in recent-onset T1D had significantly
   correlated effects in multiple ND AAB^+, although not in single ND
   AAB^+ ([153]Fig. 3F). At the pathway level, antigen processing and
   presentation were up-regulated in both single and multiple ND AAB^+,
   and interferon signaling was up-regulated in multiple ND AAB^+
   ([154]Fig. 3, G and H, and table S9). Among key genes in these
   pathways, MHC class I genes and interferon signaling IRF TFs were
   up-regulated in multiple but not in single ND AAB^+ ([155]Fig. 3I and
   data S7). We also identified pathways altered specifically in single
   and multiple ND AAB^+ and not in recent-onset T1D; for example, heat
   stress response was up-regulated in single and multiple ND AAB^+,
   extracellular matrix (ECM) organization, cytokine-cytokine
   interactions, and G protein–coupled receptor (GPCR) ligand binding were
   all down-regulated in multiple ND AAB+, and transforming growth
   factor–β (TGF-β) signaling was down-regulated in single ND AAB^+
   ([156]Fig. 3H and table S9). In addition, class II MHC antigen
   presentation was strongly up-regulated in multiple ND AAB^+, but not in
   single ND AAB^+, including class II MHC genes HLA-DBP1 and HLA-DRB1
   (data S7 and table S9). These results highlight that single and
   multiple ND AAB^+ have both shared and distinct genomic changes in β
   cells compared to T1D.

   Changes have been reported in the exocrine pancreas in T1D and at-risk
   individuals ([157]33), and in our study, we observed marked changes in
   exocrine cell gene expression in T1D progression. In basal acinar
   cells, there were 255 genes with altered expression in recent-onset
   T1D, almost all of which (95%) had decreased expression ([158]Fig. 3A
   and data S7). Basal acinar and other acinar subtypes showed
   down-regulation of numerous pathways in recent-onset T1D including
   those related to signaling, stimulus response, metabolism, and protein
   transport ([159]Fig. 3A and table S9). In multiple ND AAB^+, the
   high-enzyme acinar subtype showed significantly higher expression of
   antigen presentation, interferon signaling, and immune-related
   pathways, as well as increased activity of amino acid metabolism, which
   is necessary for enzyme production, carbohydrate and glucose
   metabolism, transcriptional activity, and respiration ([160]Fig. 3J and
   table S9). We also observed down-regulation of genes in ductal cells in
   T1D associated with small molecule transport, stimulus response,
   cytokine signaling, and RNA processing but no evidence for changes in
   ND AAB^+ (table S9).

   Other cell types in islets and the surrounding microenvironment also
   had significant changes in activity across entire pathways during
   progression to T1D. In α cells, antigen presentation, interferon
   signaling, and other pathways were significantly increased in T1D with
   less pronounced effects for these pathways in multiple ND AAB^+ and, in
   contrast to β cells, little change in single ND AAB^+ ([161]Fig. 3J and
   table S9). δ Cells showed more prominent changes in multiple ND AAB^+,
   including significantly increased hypoxia and heat stress response and
   cell cycle–related pathways and decreased signaling pathways, as well
   as in single ND AAB^+ ([162]Fig. 3J and table S9). In endothelial
   cells, we observed increased interleukin-2 (IL-2) and JAK-STAT
   signaling as well as Stem Cell Factor (SCF)/KIT signaling, which
   promotes angiogenesis ([163]34, [164]35), in recent-onset T1D
   ([165]Fig. 3J and table S9). In activated stellate cells, there was
   decreased expression of translation and RNA processing in ND AAB^+ and
   down-regulation of many pathways in recent-onset T1D (table S9). We
   observed few significant changes in gene or pathway activity in immune
   (T cell and macrophage) cells, although this could be due to the small
   number of cells profiled for these cell types.

   Together, these results reveal key genes and pathways altered in
   pancreatic cell types in ND AAB^+ and T1D donors with both shared and
   distinct changes in ND AAB^+ compared to T1D, which in ND AAB^+
   included antigen presentation, interferon signaling, ECM-related and
   stress response pathways in β cells, and metabolism and immune
   signaling in acinar cells.

Changes in the pancreatic cell type–specific epigenome in T1D progression

   We next examined to what extent altered gene and pathway activity in
   pancreatic cell types in T1D progression is driven by changes in the
   epigenome using snATAC-seq profiles from 30 ND, ND AAB^+, and T1D
   donors from nPOD. First, we identified cREs in each cell type with
   altered activity in T1D progression using a linear mixed model to
   account for pseudo-replication (see Materials and Methods). We observed
   significant changes (FDR < 0.1) in cRE activity in ND AAB^+ (single and
   multiple) and T1D for most pancreatic cell types (data S8). β Cell cREs
   with increased activity in recent-onset T1D were significantly enriched
   (FDR < 0.1) for sequence motifs of steroid hormone receptors (NC3C1 and
   NR3C2), NF-Y (NFYA, NFYB, and NFYC), interferon response factors (IRF2
   and IRF7), and stress-response TFs (ATF4, STAT1, and CEBPG) among
   others ([166]Fig. 4A and table S10). Conversely, cREs with decreased
   activity in T1D were significantly enriched for sequence motifs of TFs
   involved in core β cell functions, such as HNF1 and RFX, with many β
   cell identity TFs (NKX6.1 and PDX1) and other TF families including
   FOXA and MEF showing more nominal enrichment ([167]Fig. 4A and table
   S10). We also identified sequence motifs enriched in β cell cREs
   altered in ND AAB^+, including IRF, TCF, and STAT TF motifs in cREs
   with increased activity and MEF, RFX, and NFAT TFs in cREs with
   decreased activity, although other T1D-associated motifs such as HNF1
   showed no change in ND AAB^+ ([168]Fig. 4A and table S10). Sequence
   motifs were also enriched cREs altered in T1D progression for other
   pancreatic cell types, such as MEF and RFX TF motifs in α cells, RUNX
   TF motifs in activated stellate cells, STAT TF motifs in endothelial
   cells, and FOS/JUN motifs in ductal cells.

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

   [169]Fig. 4.
   [170]Open in a new tab

   (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.

   We determined next whether TF motifs enriched in T1D-associated cREs in
   pancreatic cell types had broader, genome-wide changes in activity in
   T1D progression by modeling sequence motif accessibility across
   individual cells using chromVAR (see Materials and Methods) ([171]26).
   In β cells, we observed consistent changes in the genome-wide
   accessibility of specific sequence motifs in T1D progression, including
   increasing accessibility of IRF motifs and decreasing accessibility of
   RFX, FOXA, and MEF motifs from ND AAB^+ to T1D states ([172]Fig. 4B and
   table S11). In other cases, sequence motifs had different patterns in
   ND AAB^+ and T1D, such as decreased accessibility of HNF1 and increased
   accessibility of PAR-related and hormone receptor TFs in T1D only and
   opposed accessibility of SIX TFs in ND AAB^+ and T1D. While α cells
   showed similar increases in accessibility of hormone receptor, stress
   response, and PAR-related TFs in T1D progression as in β cells, there
   were also several marked differences such as increased accessibility of
   MEF and RFX motifs in ND AAB^+ and recent-onset T1D, respectively
   ([173]Fig. 4B and table S11).

   We used TF GRNs to determine which TFs drive changes in pathway
   activity in T1D progression. In β cells, pathways altered in ND AAB^+
   and T1D had highly specific links to TF GRNs, suggesting key regulators
   of pathway activity in T1D progression ([174]Fig. 4C and data S6). For
   example, pathways up-regulated in β cells in T1D and ND AAB^+ such as
   interferon signaling were linked to GRNs for IRF TF motifs, and antigen
   processing and presentation were linked to NFY, IRF, and nuclear factor
   κB (NF-κB) TF GRNs, while down-regulated pathways in T1D such as ER and
   Golgi-related processes were linked to CREB3L1, XBP1, and other TF
   motifs ([175]Fig. 4C). We also identified TF GRNs linked to pathways
   altered specifically in ND AAB^+, such as heat stress–related pathways
   and heat shock factor (HSF) TF GRNs, ECM-related pathways and ETS, ELK
   and ELF TFs, and GPCR signaling pathways and RFX and FOXA GRNs
   ([176]Fig. 4C and data S6). While we observed a strong change in HNF1
   motif accessibility, as well as HNF1A expression, in β cells in T1D
   ([177]Fig. 4, B and D), no pathways linked to the HNF1 GRN had a
   significant change in expression in T1D. However, there was a more
   nominal change in β cell development and function pathways linked to
   the HNF1 GRN in T1D ([178]Fig. 4D), supporting that reduced HNF1
   activity likely underlies altered β cell function in T1D, as has been
   shown in the context of type 2 diabetes ([179]36).

   Similarly, in other pancreatic cell types, TF GRNs were linked to
   pathways with altered activity in ND AAB^+ or T1D. For example, in
   enzyme-high acinar cells, metabolic pathways altered in ND AAB^+ were
   linked to GRNs for specific TFs such as glucose metabolism and HNF1,
   amino acid metabolism and STAT1, and oxidative phosphorylation and MEF
   and FOS TF GRNs ([180]Fig. 4E and data S6). In activated stellate
   cells, fibrin-related pathways up-regulated in ND AAB^+ were
   significantly linked to ELK, HOX, CEBP, and other TF GRNs. In
   endothelial cells, IL-2 and JAK-STAT signaling pathways up-regulated in
   T1D were strongly linked to NF-κB (REL and RELA) and IRF TF GRNs, and
   SCF/KIT signaling was also linked to HOX family TF GRNs, among others.
   We further explored changes in TF activity inferred from spatial gene
   expression profiles of TF GRNs across cell types, which revealed
   increased activity of immune regulation, inflammation and signaling TFs
   (e.g., STAT3, RBPJ, FOSL2, and JUND), and reduced activity of
   endocrine-related TFs (e.g., PAX6, GLI3, MAFA, INSM1, and NEUROD1), in
   T1D compared to nondisease (fig. S12).

   We next annotated specific β cell cREs altered in T1D progression with
   putative target genes and assessed changes in regulatory programs at
   specific loci. There were 114 β cell cREs with altered activity in T1D
   progression linked to genes with significant changes in expression. For
   example, a β cell cRE on chromosome 14 in the first intron of TSHR had
   increased accessibility in recent-onset T1D and was linked to TSHR,
   which had among the largest increases in expression in recent-onset T1D
   ([181]Fig. 4F). We identified similar cREs up-regulated in recent-onset
   T1D linked to genes with highly up-regulated expression including HLA-A
   ([182]Fig. 4F), as well as CD74, GAD1, IL15, and STAT1/4. In other
   cases, we observed epigenomic changes in β cells that may precede
   changes in expression of cognate target genes. For example, a cRE
   upstream of IAPP had reduced accessibility in early T1D, although IAPP
   itself only had a significant decrease in expression in longer-duration
   T1D.

   Given pathways and transcriptional regulators with altered cell type
   activity in T1D progression, we determined whether any changes before
   T1D onset showed evidence for a role in genetic risk of T1D. We tested
   for enrichment of cREs linked to genes in each pathway for
   T1D-associated variants genome-wide (excluding the MHC locus) using
   fgwas (see Materials and Methods) ([183]5, [184]37). In β cells,
   several pathways altered in ND AAB^+ were enriched for T1D-associated
   variants including antigen processing and presentation (log
   enrich = 4.48), class II MHC antigen presentation (log enrich = 4.74),
   and interferon signaling (log enrich = 6.00) as well as several
   extracellular interaction-related processes (focal adhesion and laminin
   interactions) and GPCR signaling ([185]Fig. 4G). By comparison,
   multiple other pathways previously implicated in driving T1D risk in β
   cells, such as apoptosis, autophagy, mitophagy, and senescence, showed
   limited to no enrichment ([186]Fig. 4G). Among other cell types, we
   found evidence for enrichment of immune, metabolism, and
   transcription-related pathways in high-enzyme as well as basal acinar
   cells (fig. S13).

   We further identified specific T1D risk loci that may alter regulatory
   activity of disease-enriched pathways in key cell types such as β
   cells, T cells, and other immune populations and exocrine cells during
   T1D progression. We identified candidate causal variants at known T1D
   loci by overlapping cREs altered in T1D progression with published
   fine-mapping data ([187]5). In β cells, multiple candidate causal
   variants at the IRF1 locus overlapped cREs with increased activity in
   T1D including at the promoter and downstream of IRF1 ([188]Fig. 4H and
   table S12). There was increased β cell expression of IRF1 through
   stages of T1D progression, and IRF1 is a driver of β cell interferon
   responses, which is a pathway broadly enriched for T1D-associated
   variants ([189]Fig. 4H). Conversely, at the STAT4 locus, we identified
   cREs with increased activity in β cells as well as T cells, although
   candidate causal variants for T1D at the STAT4 locus only overlapped
   cREs active in T cells ([190]Fig. 4I). This finding supports that while
   increased STAT4 activity in β cells is observed in T1D, the STAT4 locus
   more likely affects T1D risk through altered T cell function.

   Proinflammatory cytokine exposure has been used as a model to study
   genetic risk of T1D in β cells ([191]18, [192]38), but the extent to
   which cytokine-induced changes in β cells captures transcriptional
   regulators and sites of T1D risk in actual T1D progression is not well
   established. Several TF motifs enriched in cytokine-responsive cCREs
   mirrored those enriched in cCREs up-regulated in T1D including IRF and
   STAT TFs, while other TFs such as steroid hormones may be more specific
   to T1D progression. We next compared the overlap of T1D risk loci with
   cCREs altered in T1D progression and cytokine-responsive cCREs from a
   prior study ([193]38). Risk variants at a subset of T1D loci overlapped
   both sets of cCREs, including the IRF1 locus. These results suggest
   that cytokine exposure partially recapitulates epigenomic changes in β
   cells during T1D progression and may help model disease mechanisms at
   specific loci. Together, these results reveal transcriptional
   regulators and networks altered in T1D progression, including those
   regulating pathways that likely play a causal role in T1D such as
   antigen presentation, interferon signaling, and extracellular
   interactions in β cells.

Changes in pancreatic cell-cell signaling in T1D progression

   External signaling between cell types is a key driver of changes in
   cell type–specific regulation and function, and therefore, we lastly
   identified cell-cell signaling interactions in the pancreas altered in
   T1D progression. We first inferred cell-cell interactions using
   snRNA-seq data for 32 nondiabetes, ND AAB^+ (both single and multiple),
   and T1D donors in nPOD using 1939 ligand-receptor (LR) pairs in
   CellChat (see Materials and Methods) ([194]39), which revealed 87,650
   interactions significant (FDR < 0.1) in at least one condition (table
   S13). Grouping ligands into functional categories revealed classes of
   outgoing signals preferentially produced by each cell type; for
   example, hormones, neuropeptides, and cell adhesion molecules from
   endocrine cells and enzymes from exocrine cells (fig. S14 and table
   S14).

   We identified cell-cell interactions with changes in activity in T1D
   progression using a permutation test and considered changes significant
   at FDR < 0.1 (see Materials and Methods). Overall, there was a
   reduction in the number and strength of interactions in recent-onset
   and long-standing T1D compared to nondiabetes, which was largely driven
   by exocrine cells ([195]Fig. 5A). In both ND AAB^+ and recent-onset
   T1D, there was increased strength of interactions involving endocrine
   cells and other cell types, although the total number of interactions
   was reduced ([196]Fig. 5A). We further identified cell-cell
   interactions among cells in spatial niches and determined changes in
   T1D using spatial transcriptome profiles. We identified spatially
   coexpressed LR interactions by Moran’s bivariant extension in SpatialDM
   ([197]40) using LR pairs from CellChat ([198]39). We compared the
   number of detected interactions considering each FOV as technical
   replicates of a donor and observed significant heterogeneity across
   donors and, like dispersed cell data, fewer interactions in T1D
   compared to ND donors (H statistic = 19.6, P = 0.0015) ([199]Fig. 5A).

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

   [200]Fig. 5.
   [201]Open in a new tab

   (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.

   Among specific cell types, endocrine cells displayed significant
   increases in both outgoing and incoming signaling in recent-onset T1D
   ([202]Fig. 5, B and C). We also observed significant increases in
   incoming signaling to endothelial, ductal, and activated stellate
   cells, as well as nominal changes in basal and high-enzyme acinar,
   immune, and stellate cells, in recent-onset T1D. Summarizing signaling
   by functional category revealed broad classes of cell type–specific
   signals altered in T1D; for example, β and other endocrine cells had
   increased signaling from cell adhesion molecules, whereas T cells had
   increased antigen presentation and interleukin signaling ([203]Fig.
   5D). We further examined changes in signaling between specific pairs of
   cell types in T1D progression (table S15). Significant changes
   (FDR < 0.1) in recent-onset T1D included increased incoming and
   outgoing signaling involving β cells, including between β cells
   themselves ([204]Fig. 5C), as well as increased signaling for α cells,
   outgoing signaling from high-enzyme acinar cells, and incoming
   signaling to endothelial cells.

   Given the importance of external signaling to β cells in T1D, we
   focused specifically on signals involving β cells. In recent-onset T1D,
   autocrine/paracrine signals incoming to β cells with significant
   changes in activity included cell adhesion molecules NRXN1, CADM1, and
   NEGR1 from all endocrine cell types and the secreted factor BMP5 from β
   cells ([205]Fig. 5E). In addition, high-enzyme acinar cells had
   increased signaling of trypsinogen ([206]Fig. 5E), and stellate cells
   had increased signaling of ECM and cell adhesion molecules to β cells.
   Among immune cells, signals with significant changes in signaling to β
   cells included GZMA and CCL5 from T cells and VSIR, PGRN, and LGALS9
   from macrophages ([207]Fig. 5E). In return, β cells had increased
   signaling of IL7 and MHC class I genes HLA-A and HLA-C to T cells, as
   well as increased signaling of BMP5, EFNA5, DLK1, and ANGPTL2 to
   macrophages. Notably, multiple β-immune cell signals altered in a T1D
   map to T1D risk loci (e.g., DLK1, HLA-A, HLA-C, and IL7R) ([208]5).

   We next identified differential interactions (P < 0.05) in spatial
   profiles by performing a likelihood ratio test, which provided support
   for many T1D-associated interactions identified in dispersed cell data.
   For example, interactions involving HLA class I (e.g., HLA-C), APP,
   SPP1, and BMP5, as well as ECM-related interactions, were altered in
   T1D ([209]Fig. 5, F and G). We also identified additional interactions
   enriched in T1D donors, for example, between migration inhibitory
   factor MIF and its transmembrane receptor CD74, consistent with
   previous studies ([210]41), and involving several chemokines. Next, we
   identified spatially coexpressed LR pairs using the Moran’s I score in
   Liana+ ([211]32). We obtained the top interactions associated with each
   niche using non-negative matrix factorization (see Materials and
   Methods). In T1D, an interaction between APP and CD74 was enriched in
   the endocrine niche, where APP is involved in inflammation and could
   promote immune responses in T1D (fig. S15). Conversely, interactions
   involving INS, IGF1R, INSR, and CALM1, among others, were enriched in
   the endocrine niche from nondisease donors (fig. S15).

   Several ligands that signal to β cells in T1D progression including
   BMP5 and PGRN have not been previously implicated in T1D, and cellular
   responses to these ligands may contribute to changes in gene activity
   observed in β cells in T1D. We thus determined the effects of in vitro
   exposure to these ligands on gene expression using the β cell model
   EndoC-BH1. Exposing β cells to BMP5 in culture revealed 1926 genes with
   significant change (FDR < 0.1) in expression, where the most
   up-regulated genes were ID1-4 and SMAD6-7, known targets of BMP that
   regulate proliferation and differentiation, and the β cell identity
   gene MAFA ([212]Fig. 5H and table S16). More broadly, BMP5 exposure
   up-regulated pathways (FDR < 0.1) related to TGF-β signaling,
   glycolysis, secretion, and lipid metabolism and down-regulated pathways
   such as antigen presentation and chemokine signaling ([213]Fig. 5H).
   Second, PGRN encodes secreted proteins produced by macrophages and
   ductal cells. Upon exposure to granulin, 491 genes had a significant
   change (FDR < 0.1) in expression including up-regulation of β cell
   function and insulin secretion genes MAFA, ISL1, SOX4, CRY2, and
   down-regulation of apoptosis-related genes PEA15, PDCD5, and CCAR1
   ([214]Fig. 5I and table S17). More broadly, granulin up-regulated
   cholesterol and glycerolipid metabolism and down-regulated interleukin
   signaling and inflammation, transcription and translation, and cell
   death. Overall, these results support that both BMP5 and PGRN suppress
   pathways up-regulated in β cells in T1D and therefore may play a
   protective role in β cells during T1D progression. These results
   together reveal broad changes in predicted cell-cell signaling in T1D
   progression most prominently among endocrine cells and niches but also
   involving other cell types, including signals altered in T1D that
   modulate T1D-relevant regulatory programs in β cells.

DISCUSSION

   Single-cell and spatial profiling of human pancreas donors revealed
   extensive changes in the abundance, regulation, and signaling of
   specific cell types in T1D progression, including processes that play a
   likely causal role in driving disease. In β cells, class I and class II
   MHC antigen presentation and interferon signaling pathways, TF
   regulators of these pathways, and cREs linked to genes in these
   pathways all had up-regulated activity in recent-onset T1D and ND
   AAB^+. Antigen presentation was altered as early as single ND AAB^+
   donors, suggesting that aberrant antigen presentation in β cells may be
   an initial triggering event in T1D. A larger proportion of β cells from
   ND AAB^+ and T1D donors expressed antigen presentation pathways,
   although whether subsets of these cells drive initiation or
   exacerbation of immune responses requires further investigation.
   Antigen presentation and interferon signaling pathways in β cells were
   broadly enriched for T1D-associated variants, and specific risk loci
   for T1D were linked to key genes in these pathways such as IRF1. In
   contrast, we found limited evidence that pathways directly related to
   apoptosis, as well as other processes implicated in T1D in β cells such
   as autophagy, senescence, and mitophagy, harbor T1D risk and are thus
   more likely consequences of disease. It has been long hypothesized that
   β cells affect genetic risk of T1D through increased cell death
   ([215]41–[216]47). Our results support that β cells may primarily
   contribute to T1D risk via the initiation or exacerbation of immune
   responses, which necessitates different cellular models and phenotypic
   readouts to understand their role in disease.

   In addition to shared pathways, gene activity in β cells and other
   pancreatic cell types had distinct changes in ND AAB^+ compared to
   recent-onset T1D, revealing that genomic profiles before T1D onset are
   only partially intermediate to those in T1D. In addition, the lack of
   individual genes with highly significant changes expression in ND AAB^+
   suggests that changes at these stages are likely more subtle, in
   contrast to previous reports ([217]4). Several pathways in β cells were
   altered specifically in single and multiple ND AAB^+ such as heat shock
   response and ECM organization. Heat shock responses are activated by a
   variety of stressors, promote antigen presentation in β cells, and can
   act as chaperones for antigens and thus may plausibly contribute to the
   initiation of autoimmunity ([218]48, [219]49). The breakdown of ECM is
   also an important process in T1D, as both a precursor to immune
   invasion and by affecting intrinsic β cell function ([220]50). We
   observed a similar pattern of both shared and distinct changes in the
   epigenome of β cells in ND AAB^+ compared to T1D, including increased
   NEUROD1 activity and decreased SIX TF activity. There were also shared
   and distinct features in T1D based on the duration of disease; for
   example, a more pronounced reduction in β cell function in
   long-standing T1D.

   In contrast to β cells, changes in pathway activity in α cells were
   largely restricted to multiple ND AAB+ and T1D, including antigen
   presentation and interferon response pathways and transcriptional
   regulators of these pathways. This supports that immune responses are
   more pronounced within β cells compared to α cells particularly in the
   early stages of T1D, which may reflect differences in immune targeting
   as well as the intrinsic properties of each cell type. The latter is
   supported by in vitro studies showing pronounced responses of β cells
   to external stressors ([221]18). A previous study revealed changes in α
   cell function and gene expression in single ND AAB^+ using data from
   HPAP ([222]51), although there were overall few genes with altered
   expression which supports our findings that genomic changes in α cells
   before T1D are likely subtle. In addition, several TF families such as
   RFX and MEF2 had different patterns of accessibility between α cells
   and β cells in T1D, further highlighting the unique responses of each
   cell type to disease. Conversely to α cells, δ cells had altered
   activity of multiple stress and inflammatory response pathways both in
   single and multiple ND AAB^+, as well as decreased abundance in T1D,
   suggesting that they may play an as-of-yet unappreciated role in T1D
   progression.

   Given that we profiled pancreatic tissue samples and not purified
   islets, our study was uniquely placed to reveal changes in the exocrine
   pancreas compared to previous single-cell studies ([223]4, [224]15). We
   identified multiple clusters of acinar and ductal cell types which had
   distinct genomic profiles and may represent heterogeneous subtypes of
   these cell types. In acinar cells, subclusters were broadly related to
   enzyme production and signaling responses, and previous reports
   highlighted similar heterogeneity in secretory and idling acinar cells
   ([225]52). Similar hormone-producing and signaling states have been
   reported in endocrine cells ([226]53) and thus may represent a common
   property of secretory cells. Resolving exocrine subclusters revealed
   genomic changes within specific exocrine subtypes in T1D. Enzyme-high
   acinar and MUC5B^+ ductal cells were more abundant in ND AAB^+ donors,
   and acinar subtypes had altered immune signaling, metabolism, and
   transcriptional pathways, as well as increased signaling to β cells, in
   T1D progression. Specific pathways within acinar cells altered in T1D
   progression also harbored T1D-associated variants, further supporting a
   role for exocrine pancreas in T1D risk ([227]4, [228]5) and providing
   new in-roads to determine how cellular processes in acinar cells
   contribute causally to T1D.

   Signaling relationships between pancreatic cell types revealed incoming
   and outgoing external signals during progression to T1D. Cell-cell
   signaling between immune and β cells highlighted known signals in T1D
   ([229]54, [230]55), as well as potential mechanisms of genes implicated
   in T1D risk such as DLK1 and IL7 ([231]5). Additional signals incoming
   to β cells in T1D such as BMP5 and PGRN have no prior known role in
   disease. BMP5 has increased autocrine/paracrine signaling in T1D and in
   vitro suppressed antigen presentation- and chemokine-related genes and
   enhanced expression of several genes linked to β cell proliferation and
   function. Other BMP proteins have been shown to both enhance and
   inhibit β cell function, maturity, and proliferation ([232]56–[233]58),
   where the direction of effect depends on the level of BMP signaling.
   PGRN suppresses class I MHC expression and T cell infiltration of
   ductal adenocarcinoma cells in the context of pancreatic cancer and has
   been shown to promote proliferation in mouse models of β cells
   ([234]59, [235]60). Signaling pathways altered in T1D, particularly
   those involved in T1D risk, may represent therapeutic areas for
   preserving β cell function to prevent or reverse T1D.

   There were several limitations of the data and analyses in our study to
   highlight. Current single-cell and spatial assays provide sparse
   profiles per cell which can affect downstream analyses. For example,
   although we used data integration methods designed to overcome sparsity
   in spatial data ([236]21), there are still challenges in annotating
   cell types and subtypes among cells with sparse profiles. In addition,
   analyses such as differential expression, TF activity, and cell-cell
   communication likely underrepresent true biological signals that are
   lowly expressed and have more limited detection in single-cell and
   spatial assays. The data from single-cell and spatial assays also
   represent a limited sample of the whole pancreas. For example, our
   single-cell assays profiled the head of the pancreas, and whether cell
   type–specific changes in T1D in this region reflect those in the body
   and tail is unknown. Similarly, spatial assays profiled FOV selected on
   the basis of specific criteria which limited our ability to define
   changes in T1D more broadly across the pancreas. Overcoming these
   limitations with current technology could consist of repeated assays
   per donor covering many different regions and FOV. Last, cell-cell
   interactions in our study represent predictions based on gene
   expression, which may not necessarily reflect ligand and receptor
   protein abundance nor capture direct LR interactions, and additional
   experiments are needed to validate these interactions.

   Last, one key area for future studies to address is the still
   relatively limited profiling of donors across all dimensions of T1D
   pathogenesis and progression. For example, we profiled few donors from
   nPOD that were single ND AAB^+, which prohibited understanding changes
   in nonendocrine cells in the earliest stages of disease initiation. In
   addition, while we grouped ND donors by the number of autoantibodies,
   there is further granularity in how T1D stages can be defined. For
   example, stage 2 of T1D is marked by both multiple autoantibody
   positivity and reduced β cell function ([237]61), and refined
   characterization of T1D stages may help reveal drivers of dysglycemia
   in T1D progression. The extent of differences in genomic profiles of
   heterogeneous subgroups of T1D, for example, based on age of onset, HLA
   background, or other variables ([238]62–[239]64), is also largely
   unknown. Continued collection of single-cell and spatial data from much
   larger sample sizes will be instrumental in enabling all these
   analyses. Last, studies pairing pancreas data with other T1D-relevant
   tissues from matched donors, for example, pancreatic lymph nodes, will
   help understand immune infiltration and cross-tissue signaling in
   driving T1D in the pancreas.

   In summary, our study revealed gene regulatory changes in pancreatic
   cell types in T1D progression and highlighted pathways, regulatory
   networks, and signals that may play a causal role in T1D; efforts that
   inform both new directions for mechanistic studies and targets for
   therapies to prevent or reverse T1D. We provide these data and maps in
   visualization tools at [240]http://t1d-pancreas.isletgenomics.org to
   further enhance their utility to the research community. More broadly,
   our study highlights the utility of single-cell multiomic and spatial
   analysis to reveal insight into cellular processes underlying
   progression to complex disease.

MATERIALS AND METHODS

Sample selection

   Whole pancreas tissue was obtained from the nPOD biorepository
   according to federal guidelines with informed consent obtained from
   each donor’s legal representative. Studies were considered exempt and
   approved by the Institutional Review Board of the University of
   California San Diego. We selected 7 T1D donors with more recent onset
   (<1 year from diagnosis) and 5 T1D donors with longer duration (>5
   years from diagnosis), along with 11 age- and sex-matched ND
   individuals. We also selected nine ND donors with T1D autoantibodies
   (ND TD AAB^+), most of which had multiple antibodies although one donor
   was single GAD^+. In total, 32 donors were obtained for genomic
   profiling. For all 32 donors, frozen tissue samples were obtained from
   head of the pancreas (table S1).

Single-cell assays

Tissue homogenization

   For each sample, we homogenized roughly 40 mg of flash-frozen pancreas
   tissue using mortar and pestle on liquid nitrogen, and ground tissue
   was used as input for the different single-nucleus assays.

Generation of snATAC-seq data

   Ground pancreas tissue was resuspended in 1 ml of nuclei
   permeabilization buffer [10 mM tris-HCl (pH 7.5), 10 mM NaCl, 3 mM
   MgCl[2], 0.1% Tween 20 (Sigma-Aldrich), 0.1% IGEPAL-CA630
   (Sigma-Aldrich), 0.01% digitonin (Promega), and 1% fatty acid–free
   bovine serum albumin (BSA) (Proliant, 68700) in molecular biology–grade
   water]. Nuclei suspension was filtered with a 30-μm filter (CellTrics,
   Sysmex) and then incubated for 5 min at 4°C on a rotator. Nuclei were
   pelleted with a swinging bucket centrifuge (500g, 5 min, 4°C;
   Eppendorf, 5920 R) and washed with 1 ml of wash buffer [10 mM tris-HCl
   (pH 7.5), 10 mM NaCl, 3 mM MgCl[2], 0.1% Tween 20, and 1% BSA
   (Proliant, 68700) in molecular biology–grade water]. Nuclei were
   pelleted and resuspended in 10 μl of 1× nuclei buffer (10x Genomics).
   Nuclei were counted using a hemocytometer, and 15,360 nuclei were used
   for tagmentation. snATAC-seq libraries were generated using the
   Chromium Single Cell ATAC Library & Gel Bead Kit v1.1 (10x Genomics,
   1000175), Chromium Chip H Single Cell ATAC Kit (10x Genomics, 1000161),
   and indexes (Single Index Kit N Set A, 1000212) following the
   manufacturer’s instructions. Final libraries were quantified using a
   Qubit fluorometer (Life Technologies), and the nucleosomal pattern was
   verified using a TapeStation (High Sensitivity D1000, Agilent).
   Libraries were sequenced on NextSeq 500, HiSeq 4000, and NovaSeq 6000
   sequencers (Illumina) with the following read lengths
   (Read1 + Index1 + Index2 + Read2): 50 + 8 + 16 + 50. Libraries were
   sequenced to an average depth of 333M reads (table S1).

Generation of snRNA-seq data

   Ground pancreas tissue was suspended in 500 μl of nuclei buffer: 0.1%
   Triton X-100 (Sigma-Aldrich, T8787), 1× EDTA free protease inhibitor
   (Roche or Pierce), 1 mM dithiothreitol (DTT), and ribonuclease (RNase)
   inhibitor (0.2 U/μl; Promega, N211B), and 2% BSA (Sigma-Aldrich,
   SRE0036) in phosphate-buffered saline (PBS). Sample was incubated on a
   rotator for 5 min at 4°C and then pelleted with a swinging bucket
   centrifuge (500g, 5 min, 4°C; 5920 R, Eppendorf). The supernatant was
   removed, and the pellet was resuspended in 400 μl of sort buffer [1 mM
   EDTA and RNase inhibitor (0.2 U/μl) in 2% BSA (Sigma-Aldrich, SRE0036)
   in PBS] and stained with DRAQ7 (1:100; Cell Signaling Technology,
   7406). A total of 75,000 nuclei were sorted using an SH800 sorter
   (Sony) into 50 μl of collection buffer [RNase inhibitor (1 U/μl) and 5%
   BSA (Sigma-Aldrich, SRE0036) in PBS]. Sorted nuclei were then
   centrifuged at 1000g for 15 min (Eppendorf, 5920R; 4°C, ramp speed of
   3/3), and the supernatant was removed. Nuclei were resuspended in 18 to
   25 μl of reaction buffer [RNase inhibitor (0.2 U/μl) and 1% BSA
   (Sigma-Aldrich, SRE0036) in PBS] and counted using a hemocytometer. A
   total of 16,500 nuclei were loaded onto a Chromium controller (10x
   Genomics). Libraries were generated using the 10x Genomics, Chromium
   Next GEM Single Cell 3′ GEM, Library & Gel Bead Kit v3.1 (10x Genomics,
   1000121), Chromium Next GEM Chip G Single Cell Kit (10x Genomics,
   1000120), and indexes (Single Index Kit T Set A, 10x Genomics, 1000213
   or Dual Index Kit TT Set A, 10x Genomics, 1000215) according to the
   manufacturer’s specifications. cDNA was amplified for 12 polymerase
   chain reaction (PCR) cycles. SPRISelect reagent (Beckman Coulter) was
   used for size selection and cleanup steps. Final library concentration
   was assessed by the Qubit dsDNA HS Assay Kit (Thermo Fisher
   Scientific), and fragment size was checked using a TapeStation (High
   Sensitivity D1000, Agilent). Libraries were sequenced using a NextSeq
   500 and a Novaseq6000 (Illumina) with the following read lengths
   (Read1 + Index1 + Index2 + Read2): 28 + 8 + 0 + 90 (single index) or
   28 + 10 + 10 + 90 (dual index). Libraries were sequenced to an average
   depth of 206M reads and captured 1580 genes per nucleus (table S1).

Generation of joint single-nucleus RNA and ATAC-seq data (Multiome)

   Ground tissue was resuspended in 1 ml of wash buffer [10 mM tris-HCl
   (pH 7.4), 10 mM NaCl, 3 mM MgCl[2], 0.1% Tween 20 (Sigma-Aldrich), 1%
   fatty acid–free BSA (Proliant, 68700), 1 mM DTT (Sigma-Aldrich), 1×
   protease inhibitors (Thermo Fisher Scientific, [241]PIA32965), and
   RNasin (1 U μl^−1; Promega, N2515) in molecular biology–grade water].
   Nuclei suspension was filtered with a 30-μm filter (CellTrics, Sysmex)
   and pelleted with a swinging bucket centrifuge (500g, 5 min, 4°C;
   Eppendorf, 5920 R). Nuclei were resuspended in 400 μl of sort buffer
   [1% fatty acid–free BSA, 1× protease inhibitors (Thermo Fisher
   Scientific, [242]PIA32965), and RNasin (1 U μl^−1; Promega, N2515) in
   PBS] and stained with 7-aminoactinomycin D (1 μM; Thermo Fisher
   Scientific, A1310). A total of 120,000 nuclei were sorted using an
   SH800 sorter (Sony) into 87.5 μl of collection buffer [RNasin (1 U
   μl^−1; Promega, N2515) and 5% fatty acid–free BSA (Proliant, 68700) in
   PBS]. Nuclei suspension was mixed in a ratio of 4:1 with 5×
   permeabilization buffer [50 mM tris-HCl (pH 7.4), 50 mM NaCl, 15 mM
   MgCl[2], 0.5% Tween 20 (Sigma-Aldrich), 0.5% IGEPAL-CA630
   (Sigma-Aldrich), 0.05% digitonin (Promega), 5% fatty acid–free BSA
   (Proliant, 68700), 5 mM DTT (Sigma-Aldrich), 5× protease inhibitors
   (Thermo Fisher Scientific, [243]PIA32965), and RNasin (1 U μl^−1;
   Promega, N2515) in molecular biology–grade water] and incubated on ice
   for 1 min before pelleting with a swinging bucket centrifuge (500g, 5
   min, 4°C; Eppendorf, 5920 R). The supernatant was gently removed,
   and ~10 μl was left behind to increase nuclei recovery. A total of 650
   μl of wash buffer [10 mM tris-HCl (pH 7.4), 10 mM NaCl, 3 mM MgCl[2],
   0.1% Tween 20 (Sigma-Aldrich), 1% fatty acid–free BSA (Proliant,
   68700), 1 mM DTT (Sigma-Aldrich), 1× protease inhibitors (Thermo Fisher
   Scientific, [244]PIA32965), and RNasin (1 U μl^−1; Promega, N2515) in
   molecular biology–grade water] was added with minimal disturbance of
   the pellet, and samples were centrifuged again with a swinging bucket
   centrifuge (500g, 5 min, 4°C; Eppendorf, 5920 R). The supernatant was
   gently removed without disturbing the pellet, leaving ~2 to 3 μl
   behind. Approximately 7 to 10 μl of 1× nuclei buffer (10x Genomics)
   were added, and nuclei were gently resuspended. Nuclei were counted
   using a hemocytometer, and 18,300 nuclei were used as input for
   tagmentation. Single-cell multiome ATAC and gene expression libraries
   were generated following the manufacturer’s instructions (Chromium Next
   GEM Single Cell Multiome ATAC + Gene Expression Reagent Bundle, 10x
   Genomics, 1000283; Chromium Next GEM Chip J Single Cell Kit, 10x
   Genomics, 1000234; Dual Index Kit TT Set A, 10x Genomics, 1000215;
   Single Index Kit N Set A, 10x Genomics, 1000212) with the following PCR
   cycles: 7 cycles for preamplification, 8 cycles for ATAC index PCR, 7
   cycles for cDNA amplification, and 12 cycles for RNA index PCR. Final
   libraries were quantified using a Qubit fluorometer (Life
   Technologies), and the size distribution was checked using a
   TapeStation (High Sensitivity D1000, Agilent). Libraries were sequenced
   on NextSeq 500 and NovaSeq 6000 sequencers (Illumina) with the
   following read lengths (Read1 + Index1 + Index2 + Read2): ATAC (NovaSeq
   6000), 50 + 8 + 24 + 50; ATAC (NextSeq 500 with custom recipe),
   50 + 8 + 16 + 50; RNA (NextSeq 500, NovaSeq 6000), 28 + 10 + 10 + 90.
   Libraries were sequenced to an average depth of 104M reads for RNA-seq
   and 247M reads for ATAC-seq and captured 1789 genes per nucleus (table
   S1).

Quality control and filtering

   Single nuclei ATAC data were processed and aligned to reference genome
   hg38, and duplicate reads were removed using CellRanger ATAC (v1.1.0).
   Chromatin accessibility for each sample was quantified in 5-kb genome
   windows as previously described ([245]53). Nuclei with less than 1000
   unique ATAC-seq fragments were removed. Initial quality control was
   performed to retain cells in each sample using the following metrics
   unique usable reads > 5000, fraction promoters used > 0.01, and
   transcription start site (TSS enrichment) (TSSe) > 0.3 using scanPy
   v1.8.0. Doublets were removed using Amulet v1.0 per sample ([246]65).
   After quality control, snATAC-seq profiles from two donors (6366 and
   6459) were excluded from further analysis due to low barcode numbers
   and poor clustering.

   Single-nuclei RNA samples were processed using CellRanger (v6.0.1) with
   reference genome hg38 ([247]66). Individual samples were processed for
   quality initially by removing nuclei with less than 500 expressed
   genes. Doublets were detected for each sample using DoubletFinder
   (v2.0.3) using an expected doublet rate of 4% for all samples
   ([248]67). In effort to reduce ambient RNA contamination largely driven
   by acinar cells, we used SoupX (v1.6.2) and selected acinar marker
   genes, REG1A and PRSS1, to estimate contamination rates ([249]68). Gene
   expression count matrices were then corrected for this predicted
   contamination, and these correct counts were used for both clustering
   and downstream analysis.

   Paired multiome data were processed and aligned, and multiplet reads
   were removed using cellranger arc (v2.0.0) with the reference genome
   hg38. Individual sample quality control was done using both modalities
   to remove low-quality nuclei without a minimum of 500 expressed genes
   and 1000 ATAC-seq fragments. Ambient RNA contamination was removed
   using SoupX (v1.6.2) using the same parameters as previously described.
   Doublets were detected and removed for both modalities using
   DoubletFinder (v2.0.3) and Amulet (v1.0), with the same parameters as
   above for single modality data ([250]65, [251]67).

Generation of spatial transcriptomic data

   Pancreatic tissue from six nPOD organ donors—three with T1D (6228,
   6247, and 6456) and three ND (6431, 6339, and 6229), matched by age and
   sex—was selected for spatial transcriptomic profiling on the CosMx
   platform (NanoString, Seattle, WA). For each donor, five consecutive
   Formalin-Fixed, Paraffin-Embedded (FFPE) tissue sections from the
   pancreatic body region were cut at a thickness of 4 μm. Sections #1,
   #2, #4, and #5 were mounted on the back of VWR Superfrost Plus Micro
   Slides, centered within the scanning area. After sectioning, the slides
   were air-dried overnight at room temperature, sealed, and immediately
   shipped with desiccant and ice packs to the NanoString facility
   (Seattle, Washington), where they were processed within 2 weeks of
   receipt. Section #3 was triple-stained for CD3, insulin, and glucagon
   using chromogen-based immunohistochemical staining using the Mach2
   Double Stain 1/Mach2 Double Stain 2 HRP-AP Polymer Detection Kit
   according to the manufacturer’s instructions (Biocare Medical, Pacheco,
   CA), and chromogens used included Betazoid Diaminobenzidine (DAB)
   (CD3), Warp Red (insulin), and Ferangi Blue (glucagon; all from Biocare
   Medical). Slides were then counterstained with hematoxylin. After
   staining, the slide was digitized at ×20 magnification using an Aperio
   CS2 slide scanner (Leica Biosystems Inc., Wetzlar, Germany), and this
   image served as a reference for FOV selection during CosMx data
   processing. The FOVs were selected by prioritizing specific features
   such as insulitic islets, islets with few insulin-positive cells,
   insulin-negative islets, and areas of inflammation in acinar tissue.
   The gene panel used for spatial imaging included 1010 genes, including
   1000 genes from the Human Universal Cell Characterization RNA Panel and
   10 additional custom genes selected for this project. The imaging
   experiments using CosMX were performed at NanoString (Seattle, WA).
   Cell segmentation was performed by NanoString using Giotto ([252]69),
   which included using immunofluorescence for glucagon to mark islets,
   CD3 or CD45 to mark immune cells, and PanCK for ductal cells +
   4′,6-diamidino-2-phenylindole.

Quality control and transcriptomic clustering of segmented cells

   For downstream analysis of spatial transcriptomes, we used the Python
   toolkits Scanpy ([253]70) and Squidpy ([254]22). For each slide, we
   imported matrices containing the gene expression, metadata, and
   positions of segmented cells. We defined a unique cell name and created
   a merged anndata object with data from all the slides. We adopted a
   standard filtering strategy, removing cells with less than 10 detected
   genes and removing genes detected in less than 300 cells. We then
   normalized the counts per cell, such that every cell has the same total
   count after normalization (1 × 10^6), and we log-transformed the
   counts.

Clustering

Gene expression

   After individual sample quality control, high-quality barcodes from
   single modality snRNA-seq and the RNA modality of our multiome data
   were clustered for 40 samples (32 snRNA and 8 snRNA multiome) using
   Seurat (v4.3) ([255]71). Quality control metrics such as high
   mitochondrial percentage (>1%), high number of genes detected (>4000
   genes), and high number of RNA counts (>7500) were used to remove
   low-quality barcodes. A combined clustering was created using principal
   components (PCs) from gene expression. We used Harmony ([256]20)
   (v1.0.3) to correct the PCs for batch effects across samples, sex, and
   sequencing technology. Clusters were removed with low number of cells
   (<10 cells) and with quality metrics such as the number of detected
   genes and RNA counts lower than other clusters. Additional doublet
   cells were removed on the basis of the expression of 2+ canonical
   markers from unrelated cell types. The final clustering resolution of
   0.5 was determined empirically based on maximizing the recovery of
   known pancreatic cell types and subtypes as distinct clusters.

   We leveraged gene expression profiles specific to the wide array of
   pancreatic cells from previous work to broadly label each snRNA-seq
   cluster as one of the following types: α (GCG), β (INS), endothelial
   (PLVAP), lymphatic endothelial (FLT4), ductal (CFTR), acinar (REG1A),
   stellate (PDGFRB), and variety of immune cells including T cells
   (CD3D), macrophages (C1QC), and mast cells (KIT) (table S2). Using cell
   type markers previously used to annotate cell type and subtype
   populations such as activated stellate (COL6A3) and quiescent stellate
   (SPARCL1), we were able to annotate these clusters. We identified
   previously characterized ductal subtype MUC5b ductal cells from the
   presence of known marker genes such as MUC5B, TIFF3, and CRISP3
   ([257]52).

   Marker genes of acinar subclusters were identified using DESeq2
   ([258]72) (v1.34), followed by gene set enrichment of subcluster marker
   genes in Kyoto Encyclopedia of Genes and Genomes (KEGG)
   ([259]73–[260]75) and REACTOME ([261]76) pathways using the fgsea
   package (v1.20) in R. Briefly, this was done by first creating two sets
   of sample pseudo-bulk count matrices of SoupX corrected gene expression
   for each cell type, one set which has the summation of count per sample
   per gene for that cell type and another with the summation of counts
   per sample per gene for all other cell types. We then performed DESeq
   for each cell type by concatenating these two matrices as our input and
   using cell type as the outcome variable with sample ID as a covariate.

Accessible chromatin

   We first merged 40 samples (32 snATAC samples and 8 multiome snATAC
   samples) from 29 donors using read counts in 5-kb windows using Signac
   ([262]77) (v1.9.0). We then performed latent semantic indexing of the
   combined snATAC data using Signac ([263]77). Harmony (v1.0.3) was used
   to correct for batch effects using the covariates sample, sex, and
   sequencing technology ([264]20). Clustering was performed on the
   batch-corrected PCs using graph-based Leiden clustering. We removed
   nuclei with a TSSe score < 2 and removed clusters with less than 10
   cells or with overall lower-quality metrics, such as fraction of read
   in peaks, number of ATAC fragments per barcode, and fraction of reads
   in promoters compared to other clusters. After an initial window-based
   clustering, we called peaks using MACS2 ([265]78) (v2.2.7.1)
   (parameters: -q 0.05 --nomodel --keep-dup all) on each cluster and then
   repeated the entire clustering process using a consensus set of peaks
   merged across clusters. Additional doublets were manually removed based
   off the presence of promoter accessibility of other cell type marker
   genes. This was done using nine known marker genes (INS, GCG, REG1A,
   REG2B, CTRB2, PRSS1, PRSS2, CFTR, and C1QC); promoter region was
   considered 2 kb upstream of the TSS. Data were clustered again after
   the removal of doublets. The clustering resolution of 0.5 was defined
   empirically based on maximizing the recovery of known pancreatic cell
   types and subtypes as distinct clusters. To label cell types, we first
   assigned gene names to peaks that overlapped 2 kb upstream of TSS and
   gene body using the gene activity function in Signac and then
   determined gene activity in established marker genes for each cell type
   and subtype.

   We next performed label transfer on the snATAC object using our gene
   expression map as reference and the peak-based chromatin data as query
   in Signac. Because of the size of the chromatin data, before label
   transfer, we randomly split the barcodes in the object into smaller
   subsets. We used the 2k most highly variable features from the gene
   expression map to derive transfer anchors using canonical correlation
   analysis. These anchors were then used to transfer to the chromatin map
   using the TransferData function in Seurat (v4.3). After each subset
   object was done with label transfer, we merged the objects and
   reclustered all the chromatin data together using the same methods
   described above. Last, we removed cells with low prediction scores
   (max.predicted.score < 0.5), and all cells passing this threshold were
   labeled with the predicted cell type annotation. For acinar cells, we
   summed the prediction scores of all acinar subtypes then filtered by a
   combined acinar max.predicted.score < 0.5.

   To determine the accuracy of label transfer, we used single-cell
   multiome data where the identity of the accessible chromatin profile is
   known from the paired gene expression profile. Since the gene
   expression and chromatin profiles for these nuclei were analyzed
   separately, we could use them as an independent check. We identified
   multiome barcodes present in both chromatin and gene expression maps
   and then calculated the percentage of accessible chromatin barcodes
   with matching cell type assignments in label transfer and from the
   paired gene expression profile. Because of the limited transferring of
   subtypes in the chromatin modality, we calculated a percentage at both
   the cell type and subtype levels.

Clustering of segmented cells and cell type annotation in spatial data

   To cluster the segmented cells, we first integrated the samples using
   scVI v1.1.2 ([266]79). We performed integration by condition using the
   slide as a categorical covariate. We then used the latent
   representation to create a shared nearest-neighbor graph and compute
   Uniform Manifold Approximation and Projection (UMAP) for
   two-dimensional visualization. We performed hierarchical clustering on
   the scVI latent space at resolutions of 0.5 and 0.7, and we identified
   15 and 16 transcriptomic clusters for ND and T1D, respectively. To
   annotate cell types, we identified marker genes enriched in each
   cluster for knowledge-based cell type annotation. We detected endocrine
   cells by hormone expression, β (INS and IAPP) and α (GCG and TTR); we
   also identified exocrine cells positively expressing epithelial marker
   EPCAM, ductal (SOX9 and KRT19), and acinar (EGF, DLL1, and JAG1); we
   further annotated endothelial cells (PECAM1 and VWF), fibroblasts (VIM
   and COL1A1), immune cells (CD4 and CD8A), and mast cells (CPA3 and
   TPSAB1).

Cell type label transfer from reference snRNA-seq data

   To achieve a finer annotation on the spatial context, we transferred
   the cell type labels from the dissociated reference to the spatial data
   using spatial mapping function from moscot v0.3.5 ([267]21). First, we
   performed pseudo-bulking of dissociated data using decoupler v1.6.0
   ([268]28). We found the optimal combination of parameters for the
   spatial mapping task by hyperparameter tunning per FOV, and we used
   cosine distance between the modalities. For the annotation mapping, we
   selected the label of the annotated cell with the highest matching
   probability.

Identification of spatial cellular neighborhoods

   Cellular neighborhoods in the spatial context were computed per FOV
   using the squidpy ([269]22) function spatial_neighbors, where we used
   generic coordinates and considered 30 nearest neighbors.

Identification and annotation of multicellular spatial niches

   To identify multicellular niches, we computed the covet representation
   implemented in envi v0.3.0 ([270]23) per FOV. We used the default
   parameters, which included 64 genes to represent the covariance matrix.
   We then created a shared nearest-neighbor graph using the covet
   representation and performed unsupervised Leiden clustering with a
   resolution of 0.2. To annotate the clusters, we evaluated the relative
   cell type abundance in each group per FOV and performed hierarchical
   clustering. We aggregated “acinar basal,” “acinar high-enzyme (enz),”
   “acinar signal,” and “acinar signaling (sig)/differentiation (diff)”
   subtypes in the acinar niche, “ductal” and “MUC5b ductal” subtypes in
   the ductal niche, “α,” “β,” and “δ” subtypes in the endocrine niche,
   and “act stellate,” “Q. stellate,” “endothelial,” “macrophage,” and “T
   cells” in the connective tissue niche.

Downstream analysis

Final peak calling and signal tracks

   Cell type–specific set of chromatin peaks were derived using MACS2
   ([271]78) v2.2.7.1 on the final cell type annotations of our chromatin
   map using the following parameters -q 0.05 --nomodel --keep-dup all.
   These peak calls were used to accessible chromatin signal tracks in
   UCSC genome browser ([272]80).

Marker CREs

   Cell type–specific cREs were derived for each cell type and subtype. We
   first created a set of union peaks across the whole dataset. This was
   achieved by limiting peak size for all called peaks to 300 bp by
   centering any peaks larger than 300 bp at their summit and extending
   coordinates 150 bp in either direction. We grouped peaks based on
   overlap to create clusters of peaks. Within each cluster, the peak with
   the highest read count at its summit was identified as the reference
   peak for the region. We then generated a list of peaks that did not
   overlap any of the reference peaks and iteratively identified
   additional reference peaks again until no peaks remained.

   We used this set of union peaks to calculate two sets of sample level
   pseudo-bulk matrices per cell type as follows: First, we aggregated the
   number of ATAC fragments within peaks per donor per cell type and then,
   for each cell type, created a second matrix with the summation of
   fragments from all other cell types. Normalized count matrices were
   generated by dividing the number of fragments within a peak by the
   total number of fragments for that sample in that cell type and then
   multiplying by a scaling factor (1 × 10^6). Cell type–specific
   regulatory elements were then determined for each cell type by
   comparing the normalized count matrix for a given cell type with the
   normalized count matrix of all other cell types summed together. To
   test enrichment of a given peak for each cell type, we performed a
   logistic regression model using sample ID as a covariate and corrected
   for multiple tests using the Benjamini-Hochberg correction method
   (FDR < 0.1). We limited the marker cREs per cell type to the top 5000
   cREs ranked by fold change. We performed motif enrichment of marker
   cREs for each cell type compared to a background of all cREs in the
   cell type using HOMER ([273]81) v5.0.1 and retained enriched motifs at
   FDR < 0.01. We also tested for gene set enrichment in marker cREs using
   GREAT ([274]82) (v4.0.4).

Normalized gene expression levels

   We derived normalized gene expression profiles for each cell type by
   creating aggregate count matrices by donor per cell type. Counts were
   normalized per million (CPM) by dividing the counts for each gene by
   the total counts per donor and multiplying by 1 × 10^6.

Cell type proportion changes

   We first scaled the counts for each cell type in a sample to 10,000
   cells per sample. For several cell types, we excluded samples with
   abnormally high counts (6278 for β and δ; 6393 for T cells and B cells;
   6375 for MUC5b^+ ductal cells). We then created a linear model of the
   log-transformed counts as a function of disease status (ND, ND AAB^+,
   recent-onset T1D, and long-standing T1D), age, sex, and body mass index
   (BMI), as well as a linear model without disease status. We performed
   comparison of the models using a likelihood ratio test in package
   lmtest in R and considered P values from the test significant at 0.05.

Differential gene expression

   To determine disease-related changes in gene expression, we performed
   differential analysis using DESeq2 ([275]72) v1.34. Using snRNA-seq
   data, we derived pseudo-bulk count matrices for each cell type by
   aggregating all barcodes of a donor for each gene on a per cell type
   basis. We created the count matrices from the SoupX ([276]68)–corrected
   expression counts and then rounded counts in the matrix to the nearest
   integer. We included sex, age, and BMI, as well as proportion of β
   cells, as covariates in the model. For endocrine cell types, we
   included expression counts from scRNA-seq of 48 donors from the HPAP
   consortium ([277]31) derived from a previously created single-cell map
   ([278]15) and included an additional covariate in the model for cohort.
   For a given cell type, we only used samples with at least 20 cells,
   except for long-standing T1D β cells where we included all samples, and
   donor 6234 was excluded from analyses due to aberrant numbers of
   endocrine cells. For nonendocrine cell types, genes were tested if
   detected in at least two samples per condition and had >10 counts
   across all tested conditions, while for endocrine cell types, genes
   were tested if present in at least half of the donors per condition. We
   further excluded genes for each cell type that are established marker
   genes for a different cell type. Multiple test correction was performed
   using the Benjamini-Hochberg correction, and we considered genes
   significant at FDR < 0.1.

Differential cRE accessibility

   Using cell type–specific peak calls from MACS2 ([279]78) v2.2.7.1 per
   cell type, we created peak by barcode fragment count matrices all
   snATAC-seq donors for each disease condition. Peaks with low
   accessibility were removed from analysis, as determined by the average
   accessibility of peak across all samples less than median accessibility
   of all peaks across all samples. In addition, for each cell type,
   samples were removed with <10 barcodes in that cell type. Last, cell
   types with less than 10 cells were not used in this analysis. We tested
   each disease condition against ND using glmer ([280]83) in R using the
   logistic regression model [peak accessibility ~ disease + scale
   (FRiP) + scale (count) + (1|sample)] using a binary peak count matrix.
   We used the fixed covariates of fraction reads in peak (FRiP) and ATAC
   fragment count (count) to account for sequencing depth variation and
   used sample ID as a random effect to adjust for sample variation. We
   considered sample as a random effect to mediate pseudo-replication of
   barcodes from the same donors. Cell types with more than 30,000 cells
   were subsampled down to 10,000 for this analysis. Fold change was
   calculated by dividing the average accessibility of peaks between
   conditions. Multiple test correction was performed using the
   Benjamini-Hochberg method, and we considered cREs significant at
   FDR < 0.1.

Pathway enrichment during T1D using gene expression input

   To test for pathways enriched by disease, we performed gene set
   enrichment analysis (GSEA) ([281]84, [282]85). Using the results from
   our differential expression analysis, input genes were ranked using the
   following formula [−log[10](P value) × log[2] fold change], and the
   fgsea package (v1.20) in R was run using both KEGG ([283]73–[284]75)
   and REACTOME ([285]86–[286]92) databases [parameters: eps = 0.0,
   minSize = 0, maxSize = 1000]. Pathways were considered significant
   using FDR < 0.1.

Motif enrichment

   We used chromVAR ([287]26) to measure z-scored motif accessibility in
   snATAC-seq data. To do so, we prepared peak count data for input to
   chromVAR by converting the fixed peak sparse count matrix into a
   SummarizedExperiment and estimated Guanine/Cytosine (GC) content bias
   using chromVAR built in method ([288]26). Human TF motifs from JASPAR
   2022 ([289]93) were accessed using the JASPAR2022 Bioconductor package
   in R, and motifs were annotated to peaks using motifmatchr (v.1.21.0)
   in R. The SummarizedExperiment and motif annotations were used as
   inputs into chromVAR computeDeviations function to derive GC
   bias–corrected motif accessibility z-scores.

Motifs enriched in cell types

   TF motifs were filtered for those with an accessibility of >1.2 based
   on chromVAR built-in computeVariability function. Cell types with <50
   cells were excluded. Cell type motif accessibility z-scores were
   averaged and plotted with pheatmap (v1.0.13) and RColorBrewer (v1.1-3)
   in R.

Motifs enriched in acinar subtypes

   After subsetting the motif matrix to barcodes from acinar cells, we
   averaged motif accessibility of each acinar subtype per sample and then
   tested each motif using a two-way analysis of variance (ANOVA) across
   acinar subtypes also including a donor variable. We then calculated FDR
   from the P values using the q value package in R. To identify which
   specific subtype a significant motif was most enriched in, motifs were
   further tested using a two-way ANOVA comparing motif accessibility
   within the subtype to the average motif accessibility for the other
   acinar subtypes together also including a donor variable. P values for
   each motif were corrected by the Holm’s method. Motifs were annotated
   to subclusters based on being significant in the pan-subtype ANOVA,
   significant in the post hoc ANOVA with Holm’s correction, and having
   the highest average deviation score in the given cluster.

Motif differential accessibility

   To identify motifs with differential accessibility across disease
   states, we used a linear mixed model using the lmerTest (v3.1-3)
   package in R. We identified motifs in a cell type enriched in cREs with
   altered activity in ND AAB^+ or T1D. For these motifs, accessibility
   was modeled by barcode using encoded variables to contrast
   autoantibody, recent-onset and long-duration T1D against ND controls.
   Scaled fractions of reads in peaks and scaled number of counts were
   used as fixed effect covariates, and a random effect for sample was
   used to account for pseudo-replication. Samples with <10 cells in the
   cell type were excluded, and cell types with <50 cells or disease
   states with <20 cells and three samples were not tested. We obtained P
   values from the resulting models. Motif accessibility was averaged by
   sample and disease state to make box plots. Motif accessibility per
   condition was generated by averaging sample average motif
   accessibility, and volcano plots were generated by comparing difference
   in motif accessibility versus −log[10] q values, with a difference
   threshold of 0.25 and q value of 0.05 for dashed lines and coloring and
   labeling of samples.

Motif enrichment in differential accessible CREs

   To identify TF motifs enriched in cRE differential accessibility in
   each cell type, we used HOMER ([290]81) (v5.0.1). For each cell type,
   we identified cREs with nominal association (uncorrected P < 0.05) and
   split cREs by fold change as input and user HOMER function
   findMotifsGenome with a background of all cREs for the cell type with a
   size parameter of 200 and a masked version of the human genome hg38.
   Multiple test correction was done using the Benjamini-Hochberg method,
   and significant motifs were considered at FDR < 0.1.

ABC analysis

   To link cREs to target genes, we used ABC ([291]94) v0.2. This was done
   by first creating .bam files for each cell containing only barcodes
   from the accessible chromatin map. Since the HiC reference panel used
   was in hg19 genome build, cell type bams and peaks were converted to
   hg19 using CrossMap ([292]95) v0.6.3, and we called peaks for each cell
   type with MACS2 v2.2.7.1 using this genome build. To further improve
   enhancer activity prediction, we used publicly available Histone H3
   lysine 27 acetylation (H3K27ac) chromatin immunoprecipitation
   sequencing data for acinar, ductal, α, β, and δ cells ([293]96). We
   predicted candidate regions and enhancer activity for each cell type
   using the following flags: --peakExtendFromSummit 250,
   --nStrongestPeaks 150000 and a list of genes with nonzero expression
   (CPM > 0) in that cell type. After ABC analysis, links were converted
   back to hg38 using CrossMap. We identified genes with cell
   type–specific cRE link profiles by calculating the proportion of the
   total number of ABC links for that gene by cell type and calculating
   Shannon entropy based on the proportion.

Constructing TF GRNs

   To determine GRNs, we constructed units of TFs linked to cCREs linked
   to genes. We first used a position frequency matrix (PFMatrixList
   object) of TF DNA binding preferences from the JASPAR 2022 database