Progressive dysfunction and failure of insulin-releasing β-cells are a hallmark of type 2 diabetes (T2D). To study mechanisms of β-cell loss in T2D, we performed islet single-cell RNA sequencing of two obese mouse strains differing in their diabetes susceptibility. With mice on a control diet, we identified six β-cell clusters with similar abundance in both strains. However, after feeding of a diabetogenic diet for 2 days, β-cell cluster composition markedly differed between strains. Islets of diabetes-resistant mice developed into a protective β-cell cluster (Beta4), whereas those of diabetes-prone mice progressed toward stress-related clusters with a strikingly different expression pattern. Interestingly, the protective cluster showed indications of reduced β-cell identity, such as downregulation of GLUT2, GLP1R, and MafA, and in vitro knockdown of GLUT2 in β-cells—mimicking its phenotype—decreased stress response and apoptosis. This might explain enhanced β-cell survival of diabetes-resistant islets. In contrast, β-cells of diabetes-prone mice responded with expression changes indicating metabolic pressure and endoplasmic reticulum stress, presumably leading to later β-cell loss. In conclusion, failure of diabetes-prone mice to adapt gene expression toward a more dedifferentiated state in response to rising blood glucose levels leads to β-cell failure and diabetes development.
The development of type 2 diabetes (T2D) is characterized by a gradual decrease of pancreatic islet function accompanied by the disability to secrete sufficient amounts of insulin (1). Currently, no medical treatment can stop or reverse the loss of β-cell mass, resulting in hyperglycemia, severe secondary complications, and premature death. Hence, there is a need to improve our knowledge about the mechanisms driving islet cell failure during the progression from normoglycemia to T2D.
In humans and mice, genetics strongly contribute to islet function and individual T2D risk (2,3). New Zealand Obese (NZO) mice are used as a model for polygenic diabetes with β-cell failure, closely resembling human T2D pathophysiology (4,5). Comparing gene expression signatures in pancreatic islets of diabetes-resistant B6.V-Lepob/ob (OB) and diabetes-prone NZO mice represents an appropriate and valuable strategy for the identification of genetic variants (6,7). Whereas leptin-deficient OB mice are highly insulin resistant but protected from T2D due to adaptive β-cell proliferation, NZO mice exhibit islet cell failure and β-cell apoptosis (8,9). Several diabetes genes have been identified through combining quantitative trait loci (QTL) and islet transcriptome analysis, and their role in β-cell proliferation (Lefty1, Apoa2, Pcp4l1) (7) and cilia function (Kif3a) (10), as well as insulin secretion and β-cell survival (Gjb4) (6), has been demonstrated.
Glucose is a major modulator of β-cell function and fate, and persistent hyperglycemia exerts deleterious effects on insulin secretion and β-cell survival (11,12). The concept of glucotoxicity includes molecular mechanisms such as overstimulation, oxidative stress, and endoplasmic reticulum (ER) stress, which can induce changes in β-cell identity (de- and transdifferentiation) as well as β-cell destruction (13).
Here, we use a synchronized, 2-day glucotoxic stimulus, initiated by dietary carbohydrates after 13 weeks of ketogenic diet, to study early gene expression changes in pancreatic islets of obese mice known to undergo β-cell adaptation (OB mice) or β-cell decompensation (NZO mice) during an extended feeding regimen. Feeding the mice with a fat-enriched carbohydrate-free diet induces severe insulin resistance without affecting islet function, whereas the switch to a carbohydrate-containing diet for 2 days initiates early steps of β-cell dysfunction and failure in NZO islets in a synchronized manner (7,9,10). These early events include decreased expression of Foxo1 and cilia genes and, if extended beyond 2 days, are followed by reduced expression of Pdx1, Nkx6-1, and Mafa; reduced proliferation; and increased apoptosis (10,14). The difference in diabetes susceptibility between NZO and OB mice phenocopies the situation in humans, as there is an overlap of human diabetes susceptibility genes and differentially expressed genes (DEGs) of OB and NZO islets (9).
Recent transcriptomic and functional analyses of single pancreatic endocrine cells have been instrumental in enhancing our understanding of islet cell fate and dysfunction during diabetes (15–18). By using single-cell RNA sequencing (scRNA-seq) (17) combined with diverse bioinformatics analyses, we aim to illuminate the heterogeneity and plasticity of pancreatic islet cells and to decipher the underlying adaptive mechanisms that enable OB islets to escape β-cell loss under metabolic stress conditions.
Research Design and Methods
Five-week-old male NZO (NZO/HIBomDife) and B6.V-Lepob/ob (OB) mice from our own breeding (German Institute of Human Nutrition Potsdam-Rehbrücke [DIfE]; ethics approval no. 2347-33-2019) were fed a carbohydrate-free high-fat diet (−CH) (C 1057/89; Altromin) for 13 weeks (until 18 weeks of age), followed by random allocation into two diet intervention groups: 2 days of either −CH or a carbohydrate-rich high-fat diet (+CH) (self-made with 40% carbohydrates) (7) (Supplementary Fig. 2A).
Pancreata of OB and NZO mice (−CH and +CH groups, n = 3) were prepared for immunohistochemistry as previously described (7). Pancreatic sections were stained for somatostatin (1:1,000, MA5-16987; Invitrogen), insulin (1:50,000, I2018; Sigma-Aldrich), and glucagon (1:2,000, A0565; Dako). Primary antibodies were detected with fluorophore-labeled secondary antibodies at a dilution of 1:200 (anti-rat Alexa Fluor 546, anti-mouse Alexa Fluor 488, anti-rabbit Alexa Fluor 647; Invitrogen) for 1 h at room temperature. Images of 10–15 randomly selected islets per section were documented with a confocal microscope (TCS SP8 X; Leica Microsystems) and analyzed with ImageJ (19). In total, >6,000 islet cells were evaluated per group.
Isolation of Pancreatic Islets
Islet were isolated through standard collagenase digestion (7).
In Vitro Detection of Cellular Stress, Apoptosis, and Proliferation
One day after seeding, INS-1 823/13 cells were transduced with adenovirus expressing scrambled control shRNA (sh-scrmb) or sh-Slc2a2 (targeting both rat and mouse Slc2a2 isoforms; multiplicity of infection [MOI] 5–20) overnight. Afterward, they were incubated with regular INS-1 medium or INS-1 medium containing 30 mmol/L glucose and 0.2 mmol/L palmitate/BSA for 2–3 days. Finally, cells were lysed in radioimmunoprecipitation assay buffer and ER stress and apoptosis markers were detected via immunoblotting (6,10) (Supplementary Table 1).
Islets from C57BL/6J mice were isolated, dispersed via accutase digestion, seeded onto coverslips coated with fibronectin (F1141; Sigma-Aldrich) and extracellular matrix (E1270; Sigma-Aldrich), and cultured in RPMI medium (P04-16500; PAN Biotech UK) for 24 h. Dispersed islets were transduced with sh-scrmb or sh-Slc2a2 (MOI 20; 24 h) before incubation with regular RPMI medium or glucolipotoxic medium (RPMI medium containing 30 mmol/L glucose and 0.2 mmol/L palmitate/BSA) for 2 days. For prevention of overgrowth with fibroblasts, FibrOut (7-15174; CHI Scientific) was added to media. Finally, cells were lysed in QIAzol and RNA isolated and reverse transcribed and ER stress markers detected via quantitative RT-PCR.
Proliferation assays were performed with MIN6 cells cultured in serum-free MIN6 medium with FGF-1 (4034-50; BioVision), FGF-9 (273-F9-025; R&D Systems), or BDNF (B-250; Alomone Labs) for 3 and 24 h. BrdU labeling (100 μmol/L, 2 h) was detected via immunofluorescence staining and analyzed with ImageJ.
Single-Cell Suspension and Single-Cell Libraries
For scRNA-seq, islets from three mice per condition were pooled (34 islets per animal) and dissociated into single-cell suspension. The cell suspension was immediately used for scRNA-seq library preparation with use of the Chromium Single Cell 3′ Reagent Kit, version 3 (PN-1000075; 10× Genomics) according to the manufacturer’s instructions. Libraries were pooled and sequenced on a NovaSeq 6000 (Illumina) with an average read depth of 60,000–90,000 reads/cell.
Preprocessing of scRNA-seq Data and Identification of Doublet-Like Cell Clusters
Preprocessing was performed with the Cell Ranger pipeline (version 3.0.2; 10× Genomics), as well as python3 and Scanpy (version 1.4.4) (20). mm10 (Ensemble release 97) was used as reference genome. Cells with a mitochondrial content <20% and read counts between 1,000 and 150,000 were selected for further analysis. Lastly, we set a threshold of 500 expressed genes for live cells. An overview of cell numbers at each respective filtering step is shown in Table 1.
|.||Start .||Mt in % <20 .||Min count >1,000 .||Max count <150,000 .||Genes/cell >500 .||Doublets .|
|.||Start .||Mt in % <20 .||Min count >1,000 .||Max count <150,000 .||Genes/cell >500 .||Doublets .|
Cells failing quality control were sequentially subtracted from the total number of cells before quality control (Start). Shown are the number of cells remaining after each step. Doublets, exclusion of cell doublets; Genes/cell, number of genes detected per cell; Mt in %, mitochondrial content in percent; Max count, maximum number of reads per cell; Min count, minimum number of reads per cell.
Normalization of scRNA-seq Data
Normalization of scRNA-seq data was performed in independent steps: 1) normalizing the number of reads to the number of expressed genes with the R package scran (22) (version 1.16), 2) log-scaling of data, and 3) via batch correction with scanorama (version 1.4) (23). scanorama used 3,000 highly variable genes, detected per batch and merged afterward (scanpy.pp.highly_variable_genes). If not stated otherwise, subsequent analysis was based on those highly variable genes.
Clustering, Embedding, and Cell Type Annotation
Preliminary clustering of cell types was based on the batch-corrected whole expression profile of single cells as suggested by Büttner et al. (24). For embedding, a neighborhood graph for single cells was calculated, as k-nearest neighbor (KNN) graph, based on the first 50 principal components and considering the nearest 15 neighbors (17). Next, Louvain-based clustering was applied as implemented by Blondel et al. (25) (version 0.6.1). Variant resolutions were used for several different clusters to appropriately detect strong and weak changes in single cell gene expression profiles. Clusters were annotated based on gene expression levels of typical marker genes and hormones (Supplementary Table 2), from Baron et al. (26) with adaptations from Martens et al. (27). Plotting of cell clusters was done via UMAP (28).
Detection of DEGs
DEGs were detected with the MAST R package (29) (version 1.12) and P values corrected for multiple testing through calculation of the false discovery rate (FDR). Genes with an FDR <0.01 and a coefficient >0.2 were defined as differentially expressed. Since input data were already log-transformed, the output was not comparable with classic log (fold change) values. Thresholds similar to those described by others (17) were used for these log-transformed input data.
Gene set enrichment for Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways and Gene Ontology (GO) terms was analyzed with Database for Annotation, Visualization and Integrated Discovery (DAVID) (30).
Calculation of Cell Cycle State
β-Cell Embryonic Scores
β-Cell embryonic scores were calculated (17) based on the GSE87375 data set (36). Expression data from cells between E17.5 and P0 were used to define embryonic β-cells, and P60 was used for mature β-cells. DEGs between those two groups were calculated with log (fold change) > 0.25 and FDR < 0.01. In taking into account only the top 500 genes, scores were assigned with the tl.score_genes scanpy function.
Estimation of β-Cell Trajectories
β-Cell trajectories were calculated with slingshot (version 1.4) (37) using only β-cell–related clusters as input (Beta1, Beta2, Beta3, Beta4, Beta5, BetaP, BG, BDG). To calculate strain-specific effects related to the +CH diet, the data set was split into OB (−CH and +CH) and NZO (−CH and +CH) subsets. The Beta1 cluster was used as root node for both sets to arrange cells in pseudotime and obtain diet effects. For confirmation of the observed pseudotime, randomly selected cells with all DEGs were used to generate heat maps based on the trajectories Beta1-Beta2-Beta3-Beta4 (H1) or Beta1-Beta4-Beta2-Beta3 (H2). Finally, 1,000 permutations were performed, and which trajectory the genes correlated with significantly was determined.
Weighted Gene Correlation Network Analysis
A gene coexpression network was created from normalized and batch-corrected transcriptome data of all Beta4 cells through following the best practice workflow for weighted gene correlation network analysis (WGCNA) (38) (R package, version 1.69) (power = 3, minClusterSize = 30, and deepSplit = 2). Hub genes were defined by geneModuleMembership >0.70 among the top 10 connected genes within a module. Graphical representation was created with the R package igraph (version 1.2.6).
Detection of Differentially Coexpressed Genes Between β-Cell Clusters
The code provided by Tesson et al. (39) was used to detect differential coexpression between the Beta4 and Beta1–3 clusters. For detection of a variety of differences, a power value of 3 and a minclustersize of 30 were used.
Estimation of Cell-Cell Communication Using CellPhoneDB
Cell-cell communication was predicted with use of CellPhoneDB (version 2.0) (40). The normalized and batch-corrected gene expression matrix, as well as the cell type clustering depicted in Fig. 2B, were required to identify the most probable interactions. Interactions were visualized using the R packages ComplexHeatmap (version 2.4.3) and circlize (version 0.4.12).
Sources of Human Genome-Wide Association Study Genes and Mouse Linkage Studies
Detection and Classification of Mouse Single Nucleotide Polymorphisms
Statistics and Plotting
Statistical significance was analyzed with Student t test, one-way ANOVA, or two-way ANOVA through comparison of the test groups with the appropriate control group. Data are presented as mean ± SEM, with significance indicated by asterisks in figures (*P ≤ 0.05; **P ≤ 0.01; ***P ≤ 0.001). Where applicable, the number of replicates (n) is stated in figure legends. If not stated otherwise, enrichment analysis was performed with Fisher exact test. We created plots with R, version 3.6, using the RCircos (47,48), DiagrammeR, and gplots packages.
Data and Resources
All relevant data generated or analyzed during this study are included in this manuscript. Raw scRNA-seq data were deposited under gene accession no. GSE159211 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE159211). Further data sets used or analyzed during the current study are available from the corresponding author on reasonable request. Resources are available on request.
Different Islet Cell Composition in Diabetes-Resistant OB and Diabetes-Susceptible NZO Mice
In contrast to previous bulk RNA-sequencing studies (7,10), we investigated the contribution of different islet cell subtypes from OB and NZO mice fed a −CH or +CH diet for 2 days (10) (Supplementary Fig. 2A). The animals reacted to the diet regimen as expected, with higher blood glucose levels and plasma insulin levels on +CH diet. In line with previous reports (7,9), NZO islets were smaller than their diet-matched OB counterparts and contained significantly fewer β- and more α-cells (Supplementary Results and Supplementary Figs. 2B–F and 3A–D). Subsequently, islets from all four groups (OB −CH, OB +CH, NZO −CH, NZO +CH) were isolated, dispersed, and subjected to scRNA-seq analysis.
In total, 19,440 cells remained following quality control (Table 1). While the mean number of expressed genes increased in OB +CH compared with −CH islets, NZO +CH islets expressed the lowest average number of genes (Supplementary Fig. 4A). Next, cells from all four groups were subjected to unbiased clustering analysis and resulting clusters assigned to cell types based on known marker genes (Fig. 1A and B and Supplementary Table 2). Among the identified populations were six β-cell clusters and one α-, γ-, and δ-cell cluster, plus several polyhormonal (AGD1, AGD2, BDG, BG; abbreviations based on marker gene expression for α- [A], β- [B], γ- [G], and δ- [D] cells) and nonendocrine (acinar and macrophage, endothelial, and stellate cell) clusters (Fig. 1B). As expected, the majority of cells (63–73%) were β-cells (Fig. 1C and Supplementary Table 3). Notably, both OB and NZO mice on −CH diet showed a comparable β-cell cluster composition, with most β-cells belonging to the Beta1 cluster (56.3 and 46.2%, respectively). In addition, NZO −CH islets exhibited a considerable percentage of Beta2 and Beta3 cells and OB −CH islets a higher number of Beta5 β-cells (Fig. 1C and Supplementary Table 3). Beta5 cells expressed a very low number of genes (Supplementary Fig. 4B) and could represent low-complexity libraries.
Drastic changes in cluster distribution were observed in +CH-fed mice, in which β-cell composition changed at the expense of a greatly diminished Beta1 cluster (Fig. 1C). The β-cell cluster signature between OB +CH and NZO +CH diverged; most OB +CH β-cells ended up in Beta4, with a lower percentage in Beta3. In contrast, NZO +CH β-cells were allocated into Beta2 and Beta3 clusters at roughly equal proportions (Fig. 1C). These results suggested that β-cell adaptation versus failure was mainly presented as gene expression changes between Beta1 and Beta2–4 clusters. Because NZO mice proceed to develop stark β-cell defects after 16–32 days on a +CH diet (9), the shift from Beta1 to Beta2/3 likely represents failure to adapt to the diabetogenic stimulus. Alternatively, OB β-cells handle elevated blood glucose levels by switching to a Beta4 cell subtype. Furthermore, a cluster of proliferating β-cells, designated BetaP (where P is proliferation), was identified based on proliferation marker genes (Mki67, Pcna) and G2/M cell cycle phase genes (Supplementary Fig. 4C–E). BetaP was strongly induced by 2-day +CH feeding, particularly in OB +CH (2.6%) compared with NZO +CH (1.6%) islets (Fig. 1C and Supplementary Table 3). Whereas NZO islets contained more α- and slightly more γ-cells than OB islets, the percentage of the δ-cell cluster was highest in OB islets (Fig. 1D). Additionally, OB islets exhibited more polyhormonal cells than NZO islets (Fig. 1E), which could indicate that some OB endocrine cells can proceed toward a developmentally more immature, polyhormonal state. This could protect them from the glucotoxic effects of the +CH diet, as opposed to NZO islets, which appear less capable of doing so. Concerning nonendocrine cells, OB islets contained more acinar cells and macrophages than NZO islets (Fig. 1F).
Predisposition of NZO Beta1 Cells to Metabolic Stress
In line with previous studies, we hypothesized that the changes in NZO islets are maladaptive, leading to functional decline and cell death, while changes in OB islets generally promote cell survival and proliferation. Therefore, we addressed the question of why NZO β-cells on a +CH diet progress toward Beta2/Beta3 rather than the more protective Beta4 cluster. If genetic predisposition of NZO mice played a role, transcriptomic differences between OB and NZO β-cells might be manifest in the Beta1 cluster that is common to both strains on −CH diet. While OB and NZO Beta1 gene expression showed 98.48% (8,490 genes) overlap, 131 genes (1.52%) were differentially expressed (Supplementary Table 4). We recently showed enrichment of cilia-annotated genes in DEGs of OB and NZO islets (10), and 79 of the 131 Beta1 DEGs were annotated to cilia (Fisher exact test P < 0.00001; odds ratio 2.5) (Supplementary Table 4), 33 of which were more abundant in OB and 46 in NZO β-cells. Using linkage analysis from B6-NZO crosses (43), we located 13 of 52 noncilia genes and 26 of the 79 cilia-annotated genes to a metabolically related QTL (Supplementary Table 4), suggesting that these genes are relevant for the divergent development of β-cell clusters during the progression from −CH to +CH (Fig. 2A). The cilia-related genes are involved in different cellular processes (Fig. 2A and Supplementary Table 4); several encode ribosomal proteins (18 Rpl and 13 Rps genes), pointing at altered ribosomal biogenesis and protein synthesis, which, among other processes, is needed during cellular proliferation (49,50). Other genes implicated in β-cell function and survival, such as Isl1 (51–53) and Rgs2 (54), were also more highly expressed in OB Beta1 cells (Supplementary Table 4).
Interestingly, OB Beta1 β-cells expressed more genes linked to positive regulation of insulin secretion (Glp1r, Slc30a8, Isl1, Ucn3), translation (Rpl26, Rps5), oxidation-reduction process (Ndufa7, Ero1lb), and transport (Naca, Slc2a2) (Fig. 2B and Supplementary Table 5), implying a more metabolically active state. At the same time, genes specific for NZO Beta1 β-cells were enriched in pathways linked to protein folding (Canx, Hsp90aa1, Calr) and ER-associated ubiquitin-dependent protein catabolic process (Vcp, Ubxn4, Sec61b) (Fig. 2C and Supplementary Table 5), indicating that these β-cells were already undergoing protein misfolding or oxidative stress. Finally, network analysis by Ingenuity Pathway Analysis demonstrated that targets of the key β-cell transcription factor Pdx1 were more abundant in OB cells (Supplementary Results and Fig. 2D). Altogether, the OB-specific gene expression pattern of Beta1 cells pointed toward active, functional β-cells, whereas that of NZO Beta1 cells suggested an ongoing response to metabolic and ER stress.
β-Cell Clusters of Diabetes-Prone Mice Are Less Capable of Mitigating Effects of Diabetogenic Diet
The data above indicated that β-cells of OB and NZO mice have different capacities to react to a diabetogenic diet. To test the assumption that Beta1 mainly changes into the Beta4 cluster in OB +CH, but into Beta2 and Beta3 in NZO +CH islets (Fig. 1C), we estimated cell-to-cell distances based on similarity of gene expression using slingshot (37). Since these cell-to-cell distances were comparable with pseudotime calculations, we refer to them as such below.
To assess cell fate changes triggered by +CH diet, we defined the predominant cluster in −CH-fed mice (Beta1) as the starting point for the algorithm. Principal pseudotime ordering of the major β-cell clusters was confirmed based on DEGs (Supplementary Fig. 5A), with the most likely lineage progressing from Beta1, via Beta2 and Beta3, to the Beta4 cluster (Supplementary Fig. 5B and C). This lineage was found in both OB and NZO cells. However, while a large proportion of OB β-cells develops into Beta4 (27.7%) with few cells remaining in Beta2/3 (2.7 and 5.1%, respectively), the majority of NZO β-cells fall into clusters Beta2/3 (28.2 and 29.4%) with only few reaching Beta4 (1.9%) (Fig. 3A). Our analysis further predicted that OB Beta4 cells could develop into the proliferation cluster BetaP as well as multihormonal (BG and BDG) or Beta5 cells (Fig. 3A). For NZO islets, percentages of these populations were generally lower, and the algorithm predicted alternative branching points coming from Beta1 and Beta3 instead of Beta4 (Fig. 3A). Together, these factors indicated that similar lineage transitions occurred in β-cells of both strains, but adaptation of NZO +CH β-cell gene expression beyond Beta2/3 was stalled or impeded compared with OB +CH β-cells, a potential underlying cause for the maladaptive response of NZO islets to carbohydrate feeding.
Indeed, dynamic expression changes of several important β-cell and ER stress response genes along the pseudotime lineage indicated that NZO β-cells (Beta2/3) were less adept than OB β-cells (Beta4) at compensating for the increased carbohydrate intake on +CH diet (Supplementary Results, Fig. 3B and C, and Supplementary Fig. 6A–E). For example, Slc2a2/GLUT2 fell drastically during the OB Beta1 to Beta4 cluster transition (Fig. 3D and Supplementary Fig. 6E). Immunohistochemical staining confirmed a substantial, transient reduction of GLUT2 signal in OB +CH versus −CH islets (Fig. 3E), which partially recovered at later time points (16 and 32 days, +CH) (7,9), whereas GLUT2 decreased successively in NZO islets (Fig. 3F). For testing of whether short-term downregulation of GLUT2 provides temporary protection from the toxic effects of glucose overload and thereby protects β-cells from ER stress and apoptosis, its expression was suppressed in INS-1 cells (Fig. 4A–C). While treatment with palmitate and glucose increased the ER stress marker phosphorylated (phospho)-eIF2α in sh-scrmb–transduced cells, this effect was abolished in sh-Slc2a2–infected cells (Fig. 4D), indicating that short-term loss of GLUT2 expression reduces ER stress. In addition, control cells exhibited elevated activation of caspase 3 when treated with glucolipotoxic medium, an effect that was also lower in sh-Slc2a2–transduced cells (Fig. 4E). Because insulinoma cells are partially dedifferentiated and may already express lower levels of Slc2a2 as primary β-cells, we sought to validate our hypothesis of a protective effect due to temporary GLUT2 reduction in dispersed islets using the above conditions. Both treatment with glucolipotoxic medium for 2 days and adenoviral Slc2a2 downregulation resulted in upregulation of Atf4, Hspa5, and Hsp90b1. Thus, lower levels of Slc2a2 correlate with higher levels of these stress response genes (Fig. 4F–I), suggesting that following Slc2a2 downregulation, metabolic stress is resolved successfully, presumably due to lower glucose uptake.
Increased Cell Surface Marker Gene Expression in β-Cells of Diabetes-Resistant Mice on a Diabetogenic Diet
Dysregulated expression of receptors, transporters, channels, and cell adhesion molecules cause β-cell defects (55). We asked whether the maladaptive responses to diabetogenic diet in NZO islets were accompanied by genes coding for cell surface proteins. In the α-, β-, δ-, and γ-cell clusters, 874 of 3,702 predicted islet cell surface genes (55) were detected (mean normalized expression ±0.1). Average expression levels of cell surface protein–encoding transcripts were comparable between β-cell clusters, except for the OB +CH-enriched Beta4, which exhibited much higher expression (Supplementary Fig. 7B). For α- and δ-cells, the highest content of surface marker genes was found in the OB +CH cluster, while there was no difference between groups in γ-cells (Supplementary Fig. 7C).
In our further analysis, we focused on the Beta4 cluster because we hypothesized that its expression pattern participates in the protective adaptations of diabetes-resistant mice. Of 41 cell surface marker genes differentially expressed between Beta1 and Beta4, 36 were more highly and 5 less highly expressed in Beta4 cells (Fig. 5A). The upregulated genes were associated with ion transport (e.g., Atp1b1), N-glycosylation (e.g., Dad1), signal transduction (e.g., Il6ra), cell adhesion (e.g., Bsg), and protein biosynthesis, processing, and transport (e.g., Tmed10). The downregulated genes (Slc2a2, Scl30a8, Glp1r, Trpm5, Atp2a2) play essential roles for β-cell function, confirming a more dedifferentiated phenotype of Beta4 (Supplementary Fig. 7A). A temporary reduction of these genes in the Beta4 cluster might participate in protecting β-cells from glucotoxicity, as shown for Slc2a2/GLUT2 (Fig. 4).
To determine coexpression relationships in Beta4 cells and to identify important gene modules and hub genes, we conducted WGCNA. Overall, 18 modules were identified by the clustering algorithm (Supplementary Table 7). Module M2 contained many genes related to protein processing in ER, protein export, N-glycan biosynthesis, lysosome, and proteasome, with the antiapoptotic oligosaccharyl transferase subunit Dad1 as a hub gene (Fig. 5B and Supplementary Table 7). Among the 70 most connected genes in module M2, 69% were upregulated in Beta4 compared with Beta1 cells (Fig. 5B [circled]) and 13 genes encoded cell-surface proteins (Fig. 5B [highlighted in bold]). The majority of these DEGs were involved in metabolic processes (e.g., Ndufc2), protein biosynthesis (Eef1g), transport (Sec61b) and folding (Dnajc3), and glycosylation (Ddost). Thus, promoting protein biosynthesis and stability through co- or posttranslational N-glycosylation seems to be a key characteristic of the Beta4 cluster. Next, we applied DiffCoEx in combination with WGCNA to detect differentially coexpressed genes (DiffCoEx) unique for one of the 18 Beta4-specific modules. Three modules (M3, M6, M10) exhibited coexpressed genes (Supplementary Fig. 8A), which were specifically and highly connected in Beta4, with limited connections in Beta1–3 (Fig. 5C and Supplementary Fig. 8B and C). In M3, the strongest DiffCoEx were linked to cell proliferation (Mcm10, Bub1, Ercc6l, Cit, Kif11, and Pbk) (Fig. 5C). Candidate genes in M6 were related to DNA repair (Uhrf1, Trip13, Lig1, Dtl, and Mastl) and in M10 to interferon signaling (Ifit1, Ifit3b, Ifit3, and Irf7) (Fig. 5D).
Taken together, these data indicate that OB β-cells—by developing a Beta4 cluster expression pattern—adapt to increased insulin demand and promote survival on a carbohydrate-rich diet via an upregulation of ER chaperones, N-glycosylation, and cell adhesion genes, including many cell surface proteins.
Enhanced Cell-Cell Communication in β-Cells From Diabetes-Resistant Mice
To investigate whether cell-cell communication plays a role in the protective phenotype of OB islets, we analyzed potential ligand-receptor interactions by comparing the distribution of receptor and ligand gene expression within and across cell types using CellPhoneDB (40). A set of ligands (FGF1, FGF9, KL) identified as partners of FGFR1 was more likely to interact specifically in the Beta4 cluster (Fig. 6A), indicating that Beta4 cells could potentially promote survival and proliferation in an autocrine and/or paracrine manner. Furthermore, Beta4 gene expression for receptors activated by CSF1 (Celsr3) and BDNF (Ntrk2, Ncam1) was higher, including a putative receptor-receptor interaction between FGFR1 and NCAM1, implicating processes like cell survival, adhesion, and cell-cell communication (Fig. 6A). We next tested whether FGF1, FGF9, and BDNF have an effect on β-cell proliferation by treating MIN6 cells with these ligands. BDNF (Fig. 6B) increased β-cell proliferation, shown via higher BrdU incorporation into ligand-treated cells after 24 h; for FGF9, this effect did not reach statistical significance (Supplementary Fig. 9A and C). In contrast, FGF1 did not affect cell division (Supplementary Fig. 9B and D). These data suggest that elevated expression of BDNF receptors contributes to prevention of β-cell failure by increasing proliferation in OB +CH-enriched Beta4 cells.
Lastly, genes important for cell-cell communication exhibited different expression patterns between β-cell clusters and other endocrine cell types. Putting a focus on unique ligand-receptor interactions shared only between specific β-cell clusters and δ-cells (Supplementary Table 8), more significant interactions were predicted between δ-cell ligands and their respective receptors in Beta1 and Beta4 cells (Fig. 6C). Of note, SST signaling was likely comparable between β-cell clusters because the most abundant receptor Sstr3 was expressed at similar levels in Beta1–4 (Supplementary Fig. 9E). Finally, only few interactions were seen between δ-cells and the diabetes-prone Beta2 and Beta3 clusters (Fig. 6C), suggesting that impaired cell-cell communication is part of the diabetes susceptibility of NZO islets following a diabetogenic diet. Similarly, more ligand-receptor interactions were observed between Beta4 and the dual hormonal BG cells, whereas fewer were seen between BG and the other β-cell clusters (Supplementary Fig. 9F).
Expression Changes in Mouse Islets Are Associated With Human Diabetes Risk Genes
Finally, we assessed the relevance of our observations from diabetes-susceptible NZO mice for human diabetes (T2D). DEGs from our β-cell clusters were compared with a data set of human genes and loci associated with diabetes susceptibility. Overall, 38 of 665 genome-wide association study (GWAS) genes from the NHGRI-EBI GWAS Catalog (41) (odds ratio 2.1, P = 8.5e−5) corresponded to DEGs detected in our scRNA-seq data of mouse β-cells, which included 21 of 255 genes (odds ratio 3.1, P = 2.13e−5) from the meta-analysis by Mahajan et al. (42) (Fig. 7A). This does not only represent a significant enrichment; the majority of genes (23 of 38) were also located in B6/NZO diabetes-related QTLs (43), further illustrating their relevance for metabolic disease and validity of our analysis. Additionally, eight human diabetes genes (Sctr, Patj, Hnf1a, Plxnd1, Zzef1, Calcoco2, Poc5, Cept120 [identified with GWAS]) feature a missense variant in NZO mice but are not differentially expressed according to the scRNA-seq data set, indicating that NZO β-cells might express proteins with an altered structure or function (Fig. 7A).
Next, we were interested in the classification and distribution of GWAS genes within the β-cell clusters. Among them were six transcriptional regulators, which we hypothesized were promising candidates explaining differences between β-cell clusters and trajectory development. Three were expressed at similarly high levels in Beta1 and Beta4 (Nfat5, Nsd1, Zbtb20), one at higher levels in Beta4 (Ccnd1) and one in BetaP (Hmgb1). Hnf1a was not differentially expressed, but has a missense variant in NZO mice (Fig. 7B), while Zbtb20 was located in a QTL for pancreatic insulin. Together, these transcriptional regulators could in part mediate the differences in OB and NZO β-cell responses to a diabetogenic stimulus.
To assess whether the human diabetes gene variants were conserved in NZO mice, we mapped the 46 diabetes-related GWAS genes to NZO single nucleotide polymorphisms (SNPs). Focusing only on regulatory (down/upstream, 3′-untranslated region [UTR]) and missense variants, we found 13 genes with missense mutations in NZO, of which 9 also have a similar alteration in humans (Fig. 7C). In addition, we found a large number of NZO variants (119 SNPs in 24 genes) in the 3′-UTR, which might affect RNA stability. Altogether, our analysis points to multiple genetic factors, partially shared between mice and humans, as a possible underlying cause for differences in β-cell fates in diabetes-susceptible and -resistant mice, mimicking early T2D in humans.
In the current study, we used diabetes-resistant and -susceptible mouse models and evaluated their islet single-cell transcriptome to illuminate the transition from health to diabetes. We observed islet cell–specific alterations and various trajectories of β-cells in diabetes-prone or -resistant mice in response to a short diabetogenic challenge. This study supported conceptual findings from earlier transcriptome analyses, such as the discovery of diabetes suppressor genes in OB (Lefty, Apoa2, Pcp4l1) and diabetes genes in NZO (Ifi202b, Kif2a) (7,9,10), while enabling us to find new gene expression signatures at much increased resolution.
In our analysis, we mainly focused on β-cells, hypothesizing that OB mice are protected from diabetes due to an altered expression of cilia and antiapoptotic genes and a higher degree of maturation in the absence of carbohydrates. In addition, OB mice adapt to carbohydrate feeding via activation of genes promoting protein biosynthesis and folding, encoding for cell surface proteins and for cell-cell communication. At the same time, genes essential for β-cell function, like Slc2a2, Scl30a8, and Glp1r, were downregulated after feeding of a carbohydrate-containing diet. We believe that—at least via suppression of Slc2a2/GLUT2—β-cells protect themselves from glucotoxic conditions because INS-1 cells showed a lower degree of ER stress (phospho-eIF2α) and a reduced induction of apoptosis (cleaved caspase 3) after Slc2a2 knockdown. In contrast, NZO β-cells showed indications of metabolic pressure, which might participate in the later loss of β-cell function.
Interestingly, the cell cluster pattern of OB and NZO islets on the −CH diet was very similar. Despite higher blood glucose in NZO mice, both strains mainly exhibited β-cells of the Beta1 cluster. However, these cells responded differently to a rise in blood glucose on the +CH diet. While OB Beta1 cells transition to OB +CH-enriched Beta4 cells, NZO Beta1 cells progress into Beta2/3. We propose that the 131 DEGs of OB versus NZO Beta1 cells participate in their divergent development. NZO Beta1 cells showed indications of elevated protein misfolding and oxidative stress, which could expose them to reduced plasticity and prevent an appropriate adaptation to rising blood glucose levels. Conversely, a higher degree of maturity of OB Beta1 cells, reflected by increased expression of Glp1r, Ucn3, Slc2a2, and other PDX1 target genes (shown via IPA), appears to be a key determinant for their diabetes resistance. We propose that if the most differentiated β-cells are able to mount a higher initial response to metabolic changes (e.g., high blood glucose), OB Beta1 cells are better adapted to rising blood glucose levels.
According to pseudotime analysis, Beta1 cells attempt to take a path via Beta2 and Beta3 to become Beta4 cells in both mouse strains. However, the Beta4 cluster is more prominent in OB than in NZO islets, indicating that the transition to Beta4 occurs faster or to a higher degree in OB mice, reducing their diabetes susceptibility. For example, a notable upregulation of N-glycosylation–related genes in the OB-enriched Beta4 cluster might participate in their diabetes resistance, as N-glycosylation defects are linked to T2D pathogenesis (56,57). At the same time, important β-cell genes were downregulated, such as Mafa, Slc30a8, and Slc2a2. These alterations might protect OB β-cells from glucotoxicity, as shown in GLUT2 knockdown experiments where palmitate/glucose treatment showed attenuated effects on ER stress and the induction of apoptosis. This response is expected, since GLUT2 is the main glucose transporter in mouse β-cells and its reduction decreases glucose uptake and metabolism. Importantly, downregulation of GLUT2 in OB +CH mice is transient, since permanent loss contributes to diabetes (58). The downregulation of several typical β-cell genes in OB islets could also indicate that their β-cells transition to a state of reduced β-cell identity, thereby preventing apoptosis (59) and/or allowing proliferation. Recently, Whitticar and Nunemaker proposed an intervention that reduces overactive glucokinase, which restored pulsatility and functionality to β-cells in prediabetic mouse islets (60). Their observations indicate that such strategies mimic a physiologic response occurring in individuals without diabetes, which is or progressively becomes dysfunctional in patients with diabetes.
OB islets carry more cilia and express higher amounts of cilia genes than NZO, and cilia disassemble after 2 days of +CH feeding. This effect is accompanied by an induction of islet cell proliferation (10), a process coinciding with functional immaturity (61). Interestingly, 79 cilia-annotated genes were among the 131 DEGs of the Beta1 cluster in OB. This could indicate that functional cilia dynamics are essential for the capacity of Beta1 cells to develop into the OB-enriched cluster Beta4, of which a small fraction develops into proliferation cluster BetaP. In support of the ability of Beta4 cells to transition to a proliferative state, many MYC target genes related to ribosome function were upregulated. Together, this helps clarify why the proliferative capacity of NZO β-cells was markedly impaired following +CH in comparisons with OB β-cells (7), i.e., due to a lack of or delayed dedifferentiation toward Beta4 and BetaP. The reduction in mature β-cell expression signature in our study is transient and likely allows cells to reenter the cell cycle (transition to BetaP) while maintaining β-cell properties. It differs from pathologic β-cell dedifferentiation, including β- to α-cell conversion, that is observed during T2D progression (62–64).
Finally, the evaluation of cell-to-cell cross talk indicated that Beta1 and Beta4 cells in particular communicate with δ-cells, for instance, via the ephrin ligand EFNA5 (Beta4) and its receptor EPHB2 (δ-cells). Targeting of ephrin receptors by small molecules increased glucose-stimulated insulin secretion (65). Another noteworthy cell-cell signaling event could occur via binding of the ligand semaphorin 4D (SEMA4D) to the plexin receptor PLXNB1, which is important for tissue homeostasis and morphogenesis (66). Furthermore, BDNF receptors are highly expressed in OB-enriched Beta4 cells and treatment of MIN6 β-cells with BDNF increased proliferation, implicating it in the adaptation of OB mice to increased β-cell demand. Altogether, our analyses showed that elevated islet cell-cell communication contributes to diabetes resistance of OB mice.
The NZO mouse is a well-established model of polygenic diabetes (4), sharing many similarities with human T2D. OB mice (on a C57BL/6J background) with their monogenic defect in leptin are comparable with respect to hyperphagia, reduced thermogenesis, and physical inactivity (all slightly more pronounced in OB) (67) but maintain functional islet mass despite insulin resistance and do not develop T2D. Like OB mice, NZO show an element of dysregulated leptin signaling, as they develop peripheral leptin resistance and defective transport of leptin across the blood-brain barrier, suggesting that with progressive age and duration of HFD feeding, defects in leptin signaling in NZO become more comparable with OB mice (68). Altogether, that makes these two strains a good combination for identifying diabetes risk genes. A limitation of the inbred NZO strain is genetic homogeneity, which means that it only reflects a subpopulation of all possible human diabetes susceptibility genes. However, several diabetes-related genes discovered in NZO mice overlapped with human diabetes GWAS genes (9). Furthermore, we detected an enrichment of human T2D genes among our DEGs (Fig. 7), demonstrating the relevance of our mouse model for human T2D.
In conclusion, our intensive scRNA-seq–based analysis of islets from diabetes-resistant and -susceptible mice has helped to clarify the picture of the heterogeneous responses of islet cells exposed to diabetogenic conditions. In diabetes-prone mice, the failure to adapt gene expression leads to higher metabolic stress and β-cell failure. In diabetes-resistant mice, exposure to a diabetogenic diet induces transcriptional processes that might participate in β-cell protection. This includes activation of genes involved in N-glycosylation, cell adhesion, cell surface markers, and indications for elevated cell-cell communication.
This article contains supplementary material online at https://doi.org/10.2337/figshare.20080700.
P.G. and T.S. contributed equally.
Acknowledgments. The skillful technical assistance of Anett Helms (A.S. laboratory) is gratefully acknowledged.
Funding. The study was supported by the German Ministry of Education and Research (BMBF: DZD grant 82DZD00302) and the Brandenburg State.
Duality of Interest. No potential conflicts of interest relevant to this article were reported.
Author Contributions. P.G. performed data analysis. T.S. performed experiments and data analysis and wrote the manuscript. M.Sta. conceived and designed the experiments and participated in the drafting of the manuscript. M.Ste. performed experiments. M.B. contributed to the acquisition and analysis of data. P.M.S. performed experiments. M.J. performed data analysis. E.Z. performed data analysis. H.A. performed experiments and data analysis. S.R.B. critically edited and revised the manuscript. F.J.T. performed data analysis. H.L. performed study conception and critically edited and revised the manuscript. A.S. performed study conception and design and wrote the manuscript. All authors subsequently edited the manuscript and approved its final version for publication. A.S. is the guarantor of this work and, as such, had full access to all the data in the study and takes responsibility for the integrity of the data and the accuracy of the data analysis.
Prior Presentation. Parts of this study were presented orally at the 2022 Diabetes Congress of the German Diabetes Association (DDG), Berlin, Germany, 25–28 May 2022.