Identification of beta cells in T1D, AAb+, and ND donors from HPAP scRNAseq data

A: Integrated UMAP of islet cell types recovered from HPAP scRNAseq profiling of T1D, AAb+, and ND donors. B: Feature plot showing SoupX-corrected and normalized INS expression among islet cells. C: Dot plot of scaled marker gene expression among islet cell populations. A full table of DEGs by islet cell type is available in Supplementary Table 2. D: Beta cell counts per donor. E: Frequency of beta cells per donor as a fraction of total endocrine cells (beta cell count divided by the sum of alpha, beta, delta and gamma cells per donor). F: Volcano plot of DEGs between T1D and ND beta cells calculated using edgeR with donor age, sex, BMI and 10X Chromium kit chemistry as covariates. A pseudobulk method was used by aggregating counts in beta cells by donor (Methods). The Benjamini-Hochberg (BH) procedure was used to calculate FDR. The dashed line indicates FDR = 0.05.

Identification of a T1D enriched beta cell population through GRN inference-based clustering

A: Volcano plot of differentially active regulons calculated using pySCENIC with RegDiffusion (Methods) between T1D and ND beta cells. Regulon scores were averaged by donor. Differential activity testing was performed using a linear model with donor age, sex, BMI and 10X Chromium kit chemistry as covariates. BH-calculated FDR values are shown. The dashed line indicates FDR = 0.05. B: UMAP of beta cells clustered on regulon scores. C: Frequency of each beta cell cluster by donor. D: Frequency of each beta cell cluster by clinical status. Each point is one donor. Disease categories were compared using the Wilcoxon Rank Sum Test. BH-calculated FDR values are shown. NS: not significant. Cluster 4 was not calculated since cells were mainly from a single donor (HPAP-022). E: Violin plots of non-rank normalized regulon activity scores between cluster 3 (red) and non-cluster 3 (blue) beta cells. All panels are differentially expressed (Methods) with adjusted p values (Bonferroni) < 0.05. F-H: Feature plots of rank-normalized regulon activity scores. I: Principal component 4 score by clinical status. Each point is the average beta cell PC4 score for one donor. Disease categories were compared using the Wilcoxon Rank Sum Test. Uncorrected p values are shown. The top 6 positive and negative PC4 contributing features were CEBPD, ETS2, JUNB, ETV6, IRF1, BCL6 (positive) and MAFA, MAFB, XBP1, CHD2, ZBTB20, JUND (negative).

Transcriptional profile of cluster 3 beta cells

A: Heatmap of average expression for each beta cell cluster of the top 35 positive differentially expressed genes in C3 (by adjusted p value). Values were scaled for visualization. Cluster 4 was omitted (prior to scaling) since cells were mainly from a single donor (HPAP-022). B-E: Gene set enrichment plots for selected pathways (full list of gene sets in Supplementary Table 3). Genes were ranked by fold change in C3 beta cells (vs non-C3 beta cells), with filtering of lowly expressed genes (Methods, full table of GSEA output in Supplementary Table 10). BH-calculated FDR values are shown. NES: normalized enrichment score. F-H: Violin plots showing expression of selected transcripts for C3 and non-C3 beta cells. All panels are differentially expressed (Methods) with adjusted p values (Bonferroni) < 0.05 unless indicated by NS (not significant). I-J: Frequency of beta cells expressing TFs (I) and class I HLA genes (J) between C3 and non-C3 beta cells. For J, an inset for HLA-E is shown. All adjusted p values (Bonferroni) < 0.05 unless indicated by NS (not significant). P values derived from DEG testing. Full C3 differential expression tabular data is available in Supplementary Table 9.

In vitro cytokine stimulated beta cells are similar to cluster 3 beta cells

Beta cells treated in vitro with various stressors were identified from scRNAseq data from Maestas et al. (Methods, Supplementary Fig. 10A-H). A: In vitro beta cell expression of C3 up and down gene sets. Gene sets were synthesized by taking the top 50 non-ribosomal/mitochondrial up- and down-regulated genes (by p value) in HPAP C3 beta cells (Methods). GSVA was used to calculate pathway scores for aggregated (pseudobulk) counts for each sample. Each point is a GSVA score for one sample, colored by donor (CTRL: control, CM: cytokine mix with IFNG+IL1B+TNFA, BFA: Brefeldin A, TG: Thapsigargin). All have BH-calculated FDR < 0.05 versus CTRL unless indicated by NS (not significant). B-C: Paired expression plots for indicated transcripts in beta cells treated with cytokines (B) or BFA (C) versus control beta cells. Log2(counts per million + 1) of aggregated (pseudobulk) counts per sample are shown. Each point is the expression for one sample, colored by donor. Differential expression testing was performed using edgeR with donor as a covariate (Methods). BH-calculated FDR values are shown. Full tabular data for Maestas et al. differential gene expression is found in Supplementary Table 13. ***: FDR < 0.001, *: FDR < 0.05, NS: not significant. D: Violin plots of in vitro cytokine-treated and control beta cells (pooled from all donors). All panels are differentially expressed (Methods) with adjusted p values (Bonferroni) < 0.05. Shared ISGs are those differentially expressed in both in vivo C3 vs C 3 beta cells (HPAP) and in vitro CM-treated vs control beta cells (Maestas et al.). Non-shared ISGs are those only differentially expressed by in vitro CM-treated vs control beta cells (Maestas et al.). Full tabular data for HPAP C3 DEGs are found in Supplementary Table 9 and for Maestas et al. in Supplementary Table 13. E: Changes in beta cell motif accessibility versus transcript abundance between CM and CTRL treatment. Beta cells were isolated from Maestas et al. snMultiome data (Supplementary Fig. 11). Differential motif accessibility was calculated using chromVAR (Supplementary Table 15). Difference in mean chromVAR z-score for CM beta cells minus mean CTRL z-score is shown on the y-axis. Log2 fold change between CM and CTRL beta cells is shown on the x-axis, calculated from scRNAseq pseudobulk DEG testing (panel B). TFs with significant changes in both transcript abundance and motif accessibility are colored red. Linear regression between differential transcripts and differential motif accessibility (blue line) and Pearson correlation coefficient are shown. F: Top hits of de novo HOMER motif analysis of differentially accessible beta cell snATACseq peaks after IFNG+IL1B+TNFA stimulation from Benaglio et al. (Methods).

Identification of a homologous alpha cell population to cluster 3 beta cells

A: UMAP of alpha cells from the HPAP dataset clustered on regulon scores. B: Feature plots of rank-normalized regulon activity scores. C: Gene expression log2 fold change values for C3 beta cells versus non-C3 beta cells (x axis) plotted against log2 fold change values for C3 alpha cells versus non-C3 alpha cells (y axis). Genes were filtered to those expressed by ≥ 1% of cells in at least one test group within both comparisons. Full tabular results are in Supplementary Tables 9&19. Transcripts of regulons used to define C3 beta cells are shown in red. D: Frequency of C3 alpha cells by clinical status. Each point is one donor. Disease categories were compared using the Wilcoxon Rank Sum Test. NS: not significant. E: Frequency of C3 beta cells among all beta cells plotted against frequency of C3 alpha cells among alpha cells for each donor with a boxed inset. The Spearman correlation coefficient (ρ) is shown (p = 1.49E10-7).

Technical differences in HPAP islet scRNAseq dataset

A: UMAP plot of HPAP islet cells after dimensionality reduction (PCA on normalized transcript counts) without integration methods. Cells are colored by 10X Chromium kit chemistry. B: Feature plot of INS expression to highlight beta cells.

Identification of beta cells in T1D, AAb+, and ND donors from HPAP scRNAseq data

A: Integrated UMAP of HPAP islet cell types split by clinical status. B-C: Stacked bar plot showing frequency (B) and count (C) of each islet cell type per donor. Beta cell frequencies/counts shown here are those following post-islet clustering refinement and low feature count filtering of the beta cell population (Methods). D: Integrated UMAP of islet cell types by 10X Chromium reagent kit chemistry. E-F: Feature plots of islet expression of IAPP (E) and MAFA (F).

Intermediate beta cell UCell clustering

A: Intermediate beta cell clustering based on UCell scores (Methods). Gene sets for UCell score calculation are in Supplementary Table 3. Cells in cluster U3 (circled) were identified as likely contaminants and removed. B-G: Feature plots showing normalized expression of indicated transcripts. H: Volcano plot of DEGs between T1D and AAb+ beta cells calculated using edgeR with donor age, sex, BMI and 10X Chromium kit chemistry as covariates (as in Fig. 1F). A pseudobulk method was used by aggregating counts in beta cells by donor (Methods). BH-calculated FDR values are shown. The dashed line indicates FDR = 0.05.

Identification of a T1D enriched beta cell population through GRN inference-based clustering

A: Volcano plot of differentially active regulons between T1D and AAb+ beta cells (as in Fig. 2A). Regulon scores were averaged by donor. Differential activity testing was performed using a linear model with donor age, sex, BMI and 10X Chromium kit chemistry as covariates. BH-calculated FDR values are shown. The dashed line indicates FDR = 0.05. B: Beta cell UMAP by individual donor. C-D: Beta cell UMAP split by clinical status (C) and 10X Chromium reagent kit chemistry (D). E: Count of each beta cell cluster per donor. F: Frequency of C3 beta cells (among all beta cells) by clinical status. Each dot is one donor, shaded by donor age. G: Frequency of C3 beta cells (among all beta cells) versus donor age. Each dot is one donor, colored by clinical status. H-M: Feature plots of rank-normalized regulon activity scores.

Beta cell cluster 3 is stable across clustering resolutions

Clustree plot showing clustering solutions at different resolutions (Zappia & Oshlack, 2018). A resolution of 0.3 was used for the UMAP in Fig. 2B. Edges from parent nodes contributing less than 10% of cells to the child node are not shown.

Differential regulon activity between T1D and ND/AAb+ beta cells within each beta cell subcluster

Volcano plots of differentially active regulons between T1D and ND/AAb+ beta cells for each beta cell cluster. Regulon scores were averaged by donor for each cluster. Donors with less than 5 beta cells in a particular cluster were excluded for comparisons within that cluster. Differential activity testing was performed using a linear model with donor age, sex, BMI and 10X Chromium kit chemistry as covariates. BH-calculated FDR values are shown on the y axis. Dashed lines indicate FDR = 0.05.

Stability of RegDiffusion model for SCENIC GRN inference

Regulons were scored after independently re-training the RegDiffusion model 9 additional times on the same beta cell population as used for the original run (original run used for all analyses in the dashed box, 10 total trainings). Pairwise Pearson correlation matrices for beta cell activity of the top 6 up and top 3 down C3-enriched regulons are shown. Greyed rows/columns indicate that the regulon was not detected in a given run. The matrix diagonal is also colored grey.

Transcriptional profile of cluster 3 beta cells

A: Frequency of beta cells expressing indicated ISGs between C3 and non-C3 beta cells. B: Frequency of beta cells expressing indicated oxidative stress/DNA damage genes between C3 and non-C3 beta cells. C-H: Violin plots showing expression of indicated transcripts for C3 and non-C3 beta cells. I: Frequency of beta cells expressing lineage-defining beta cell TFs between C3 and non-C3 beta cells. For all panels, all adjusted p values (Bonferroni) < 0.05 unless indicated by NS (not significant). P values derived from DEG testing (Methods). Full cluster 3 differential expression tabular data is available in Supplementary Table 9. J: Percent of cells expressing indicated transcripts for the 3 T1D samples with the most C3 beta cells. Trends are shown (not significance) due to low cell counts in individual donors.

DEGs between T1D and ND/AAb+ beta cells within each beta cell subcluster

Volcano plots of DEGs between T1D and ND/AAb+ beta cells within each cluster calculated using edgeR with donor age, sex, BMI and 10X Chromium kit chemistry as covariates. A pseudobulk method was used by aggregating counts by donor. Donors with less than 5 beta cells in a particular cluster were excluded for comparisons within that cluster. BH-calculated FDR values are shown on the y axis. Dashed lines indicate FDR = 0.05.

Beta cell identification and gene expression from Maestas et al. in vitro scRNAseq data

A: Integrated UMAP from reanalysis of in vitro islet scRNAseq data Maestas et al. B-C: UMAP showing cells colored by donor (B) and treatment condition (C; CTRL: control, CM: cytokine mix with IFNG+IL1B+TNFA, BFA: Brefeldin A, TG: Thapsigargin). D: Dot plot of scaled islet marker gene expression within each cluster. E: Feature plots of beta cell marker transcripts. F: The number of unique genes expressed per cell within each cluster. G: UMAP showing cells identified as beta cells from Maestas et al. scRNAseq data (highlighted in red, beta cells defined as cells belonging to clusters M1, M4, M7, M8, or M15). H: Beta cell counts recovered for each sample (donor + treatment combination). I: Violin plots of in vitro cytokine-treated and control beta cells (pooled from all donors). All panels are differentially expressed with adjusted p values (Bonferroni) < 0.05. Shared ISGs are those differentially expressed in both in vivo C3 vs non-C3 beta cells (HPAP) and in vitro cytokine-treated vs control beta cells (Maestas et al.). Non-shared ISGs are those only differentially expressed by in vitro cytokine-treated vs control beta cells (Maestas et al.). Full tabular data for HPAP C3 DEGs are found in Supplementary Table 9 and for Maestas et al. in Supplementary Table 13.

Beta cell identification from Maestas et al. in vitro snMultiome data

A: Integrated ATAC and RNA UMAP from Maestas et al. islet snMultiome data. B: UMAP showing cells colored by treatment condition (CTRL: control, CM: cytokine mix with IFNG+IL1B+TNFA). C: Feature plots of beta cell marker transcripts. D: Dot plot of scaled islet marker gene expression within each cluster. E: The number of unique genes expressed per cell within each cluster. F: UMAP showing cells identified as beta cells from Maestas et al. snMultiome data (highlighted in red, beta cells defined as cells belonging to clusters A0, A1, or A3).

Intermediate alpha cell UCell clustering

A: Intermediate alpha cell clustering based on UCell scores (Methods). Gene sets for UCell score calculation are in Supplementary Table 3. Cells in cluster U2 (circled) were identified as likely contaminants and removed. B-H: Feature plots showing normalized expression of indicated transcripts.

Identification of an alpha cell population homologous to cluster 3 beta cells

A-B: Alpha cell UMAP split by clinical status (A) and 10X Chromium reagent kit chemistry (B). C: Alpha cell UMAP by individual donor. D: Count of each alpha cell cluster per donor. E: Frequency of each alpha cell cluster per donor. F-G: Feature plots of rank-normalized regulon activity scores. H-I: Gene set enrichment plots for selected pathways (full list of gene sets in Supplementary Table 3). Genes were ranked by fold change in C3 alpha cells (vs non-C3 alpha cells), with filtering of lowly expressed genes (Methods, full table of GSEA output in Supplementary Table 20). BH-calculated FDR values are shown. NES: normalized enrichment score. J: Frequency of C3 alpha cells (among all alpha cells) versus donor age. Each dot is one donor, colored by clinical status.