Characterizing pancreatic β-cell heterogeneity in the streptozotocin model by single-cell transcriptomic analysis.

OBJECTIVES: The streptozotocin (STZ) model is widely used in diabetes research. However, the cellular and molecular states of pancreatic endocrine cells in this model remain unclear. This study explored the molecular characteristics of islet cells treated with STZ and re-evaluated β-cell dysfunction and regeneration in the STZ model.
METHODS: We performed single-cell RNA sequencing of pancreatic endocrine cells from STZ-treated mice. High-quality sequencing data from 2,999 cells were used to identify clusters via Louvain clustering analysis. Principal component analysis (PCA), t-distributed stochastic neighbor embedding (t-SNE), uniform manifold approximation and projection (UMAP), force-directed layout (FDL), and differential expression analysis were used to define the heterogeneity and transcriptomic changes in islet cells. In addition, qPCR and immunofluorescence staining were used to confirm findings from the sequencing data.
RESULTS: Untreated β-cells were divided into two populations at the transcriptomic level, a large high-Glut2 expression (Glut2high) population and a small low-Glut2 expression (Glut2low) population. At the transcriptomic level, Glut2low β-cells in adult mice did not represent a developmentally immature state, although a fraction of genes associated with β-cell maturation and function were downregulated in Glut2low cells. After a single high-dose STZ treatment, most Glut2high cells were killed, but Glut2low cells survived and over time changed to a distinct cell state. We did not observe conversion of Glut2low to Glut2high β-cells up to 9 months after STZ treatment. In addition, we did not detect transcriptomic changes in the non-β endocrine cells or a direct trans-differentiation pathway from the α-cell lineage to the β-cell lineage in the STZ model.
CONCLUSIONS: We identified the heterogeneity of β-cells in both physiological and pathological conditions. However, we did not observe conversion of Glut2low to Glut2high β-cells, transcriptomic changes in the non-β endocrine cells, or direct trans-differentiation from the α-cell lineage to the β-cell lineage in the STZ model. Our results clearly define the states of islet cells treated with STZ and allow us to re-evaluate the STZ model widely used in diabetes studies.

Introduction

Diabetes mellitus is a metabolic disorder resulting from loss or dysfunction of insulin-producing β-cells in the pancreas or lack of insulin effect. Endogenous repair is a promising strategy to restore β-cell mass and relieve diabetes. Previous studies described various β-cell regeneration pathways, such as replication of pre-existent β-cells [1], intra-islet cell trans-differentiation [[2], [3], [4]], and β-cell neogenesis [5,6]. Studies of in vivo regeneration of β-cells require appropriate animal models [7,8].

The streptozotocin (STZ)-induced diabetic model is widely used in diabetes research. A single high dose of STZ results in near-total ablation of β-cells in rodents, leading to insulin deficiency, hyperglycemia, polydipsia, and polyuria, all of which mimic human type-1 diabetes [9]. However, rats that received low-dose STZ treatments combined with high-fat diets have been used as a model for type 2 diabetes [10,11]. Due to its ease of use and ability to mimic diabetic phenotypes, the STZ model has been used by several groups to study β-cell regeneration after acute β-cell injury [4,[12], [13], [14]]. One group observed δ-cell-dependent β-cell regeneration and diabetes recovery in juvenile mice after a single high dose of STZ [4]. Cheng et al. identified an intermediate cell type (vimentin+/MafB+) as β-cell progenitors following STZ-induced extreme β-cell ablation [12]. Grossman et al. studied the effect of insulin treatment in the STZ model and found that it could restore β-cell mass as well as normoglycemia after 120 days of glycemic control in female mice [13]. However, the molecular characteristics of each islet lineage after STZ treatment at the single-cell level remain undescribed, which hinders our understanding of the cell fate plasticity of endocrine lineages and β-cell state and regeneration in the STZ model.

In this study, we performed single-cell RNA sequencing (scRNA-seq) to investigate the single-cell transcriptomic profiles of all endocrine lineages over 9 months after a single high dose of STZ. We analyzed the characteristics of islet endocrine cells after STZ treatment to describe their cell heterogeneity, maturation, functional state, and trans-differentiation.

Materials and methods

Mice

Ngn3-Cre; Rosa-RFP [15] and Ins1-RFP [50] transgenic mice were used in this study. Male mice were weaned to regular chow at postnatal day 21 (P21). All of the procedures were approved by the Institutional Animal Care and Use Committee of Peking University.

In vivo STZ treatment

STZ was dissolved in 50 mM sodium citrate buffer (pH 4.5) and injected within 15 min of preparation. Two-month-old male mice were injected intraperitoneally with a single high dose of STZ (200 mg/kg, V900890, Sigma) after an overnight fast. The control mice were injected with an equal volume of sodium citrate buffer.

Insulin treatment and blood glucose monitoring

To prevent death in long-term experiments, the mice received subcutaneous insulin implants (LinBit Pr-1-B, LinShin Canada) when hyperglycemic (>400 mg/dL). The first insulin implant was administered 2 days after STZ injection. Non-fasted blood glucose levels were determined from the tail vein using a Contour Next blood glucose meter.

Insulin secretion assay

Eight-week-old male C57BL/6 mice were injected with STZ as previously described. Three or six days after STZ treatment, blood was collected from STZ-treated and control mice after fasting for 6 hours. Fasting serum insulin levels were determined using a Rat/Mouse Insulin ELISA Kit (EZRMI-13K, Millipore) according to the manufacturer's instructions.

Fluorescence-activated cell sorting (FACS)

The pancreas was perfused and digested using 0.5 mg/mL collagenase P (11213873001, Roche) as previously reported [16]. After centrifugation and removal of the supernatant, pancreatic tissue was dissociated into single cells via 0.25% trypsin-EDTA treatment for 5 min at 37 °C, and digestion was terminated with 0.4 volumes of fetal bovine serum (FBS). The cells were sorted using a BD FACS Aria SORP cell sorter.

Single-cell RNA-seq

Library preparation was conducted following the modified STRT-seq protocol [17,18]. Briefly, after cell sorting, single Ngn3-Cre; Rosa-RFP cells were transferred quickly into lysis buffer by a mouth pipette. Next, we reverse-transcribed the RNA and amplified the resulting cDNA over 18 cycles of PCR. The quality of the cDNA was assessed by qPCR of the housekeeping gene Gapdh. cDNA of good-quality single cells was pooled together and purified. Biotinylated pre-indexed primers were added to the cDNA during an additional 4 cycles of PCR. Approximately 300 ng of purified PCR product was sheared to approximately 300 bp by sonication. Sequencing libraries were prepared using a Kapa Hyper Prep Kit (KK8502, Kapa Biosystems).

Bulk cell reverse transcription quantitative PCR (RT-qPCR)

Ins1-RFP+ cells from the control and STZ-treated D12 Ins1-RFP mice were obtained by FACS. Total RNA from sorted cells was prepared using an RNAprep Pure Kit (for micro samples) (DP420, Tiangen). First-strand cDNA was synthesized using HiScript II Q RT SuperMix (R223, Vazyme). Real-time PCR was performed with 2X M5 HiPer SYBR Premix EsTaq (MF787-01, Mei5 Biotechnology).

The following primer sequences were used for RT-qPCR:

Actb: forward 5′-TCCTGAGCGCAAGTACTCTGT-3′,

reverse 5′-CTGATCCACATCTGCTGGAAG-3′;

Gapdh: forward 5′-ATGGTGAAGGTCGGTGTGAAC-3′,

reverse 5′-GCCTTGACTGTGCCGTTGAAT-3′;

Glut2: forward 5′-GTGCTGCTGGATAAATTCGCC-3′,

reverse 5′-ATTGCAGACCCAGTTGCTGA-3′;

Ucn3: forward 5′-AAGCCTCTCCCACAAGTTCTA-3′,

reverse 5′-GAGGTGCGTTTGGTTGTCATC-3′;

Ins1: forward 5′-GCAAGCAGGTCATTGTTTCAAC-3′,

reverse 5′-AAGCCTGGGTGGGTTTGG-3′;

Nkx6.1: forward: 5′-AACACACCAGACCCACGTTCT-3′,

reverse 5′-ATCCCCAGAGAATAGGCCAAG-3′;

Pdx1: forward 5′-AACTTGAGCGTTCCAATACGGA-3′,

reverse 5′-CAGCCGCCTTTCGTTATTCTTA-3′;

Vdr: forward 5′-GGCTTCCACTTCAACGCTATG-3′,

reverse 5′-ATGCTCCGCCTGAAGAAAC-3′;

Etv1: forward 5′-TGCCTAGCTGCCACTCCAT-3′,

reverse 5′-GCCTGGGACCCTTCTCAAAC-3′;

Klf10: forward 5′-AGAAGAACCCACGGAAAT-3′,

reverse 5′-GAGGAAGGCACAGCAAAG-3′;

Tspan8: forward 5′-GGAGTTCCGTTTACCCAAAGA-3′,

reverse 5′-CTCAATAACGGCCAGTCCAA-3′;

Ins2: forward 5′-TGGCTTCTTCTACACACCCA-3′,

reverse 5′-TCTAGTTGCAGTAGTTCTCCA-3′.

Immunofluorescence staining

Tissues were fixed in 4% paraformaldehyde overnight at 4 °C. For frozen sections, the tissues were cryoprotected in 30% sucrose, embedded in Optimal Cutting Temperature (OCT) media (Thermo Fisher, 6502), and frozen on dry ice. For paraffin sections, the tissues were dehydrated with 30–100% ethanol and embedded in paraffin. The sections (6 μm thick) were probed with primary antibodies against insulin (1:500, Abcam, ab7842), NKX6.1 (1:100, DSHB, F55A12-c), GLUT2 (1:50, Santa Cruz Biotechnology, sc-518022), and glucagon (1:200, Millipore, AB932) followed by incubation with Alexa Fluor 594 AffiniPure Donkey Anti-Guinea Pig IgG (H + L) (Jackson ImmunoResearch, 706-585-148) and Alexa Fluor 488 Donkey Anti-Mouse IgG (H + L) antibody (Thermo Fisher, A21202) as the secondary antibodies. Images were acquired using a Zeiss Axioimager M2 fluorescence microscope or Zeiss LSM 710 NLO and DuoScan System.

Quantification of gene expression from single-cell RNA-seq

Modified STRT-seq libraries were sequenced on an Illumina HiSeq 4000 System to generate 150-bp paired-end reads. Read1 included transcript sequences, and Read2 included an 8-bp cell-specific barcode and an 8-bp unique molecular identifier (UMI) sequence. The fastq file for each single-cell sample was split based on 8-bp cell-specific barcode sequences. Then 8-bp UMIs were added into the corresponding Read1. The polyA sequences in Read1 were trimmed. We aligned preprocessed Read1 reads to the mouse genome (mm10) with TopHat (v2.1.0) [19]. Mapped reads were annotated to genes using featureCounts (v1.5.3) [20]. The gene ID was appended as the XT tag in the bam file. After sorting and indexing the bam file with SAMtools (v.1.3.1) [21], we quantified the UMIs of each gene with UMI-tools (v0.5.0) using the following parameters: “count --per-gene --gene-tag = XT --method unique -I indexed_sorted_bam -S out_file” [22] and generated the UMI matrix. To normalize the sequencing depth, the UMI matrix was transformed to transcripts per 0.1 million (TP0.1M) (one UMI was considered as one transcript). To eliminate the influence of individual genes on the transcriptome, we ignored extremely highly expressed genes (Ins1, Ins2, Gcg, Sst, Ppy, Pyy, Iapp, and Yam1) in the calculations of the total transcripts. Single-cell samples were discarded if their mapped read count was less than 50,000 or if less than 1,500 genes were detected.

Cell type identification

To eliminate the index switching effect of the HiSeq 4000 System (5–10% of sample indexes were incorrectly assigned) [23], we subtracted 10% of the maximum TP0.1M value for each gene. Then we calculated a log2 (subtracted TP0.1M + 1) matrix to be used as the input for Seurat (v3.0.2) [24]. The top 2,000 variable genes were identified using the FindVariableFeatures function in the Seurat package. Cell cycle-related genes (including cell cycle regulators, such as cyclins and CDKs) were identified with hierarchical clustering and were excluded in downstream analyses. Principal component analysis (PCA) was performed with filtered variable genes using RunPCA in the Seurat package. The first 20 principal components (PCs) were selected for t-distributed stochastic neighbor embedding (t-SNE), uniform manifold approximation and projection (UMAP), and Louvain clustering using RunTSNE, RunUMAP, and FindNeighbors/FindClusters in the Seurat package, respectively. Cell types were assigned according to the marker gene expressions: β-cells (Ins1+; Ins2+), α-cells (Gcg+), δ-cells (Sst+), PP-cells (Ppy+), ductal cells (Spp1+), acinar cells (Cel+), and immune cells (Ptprc+).

β-cell heterogeneity analysis

The log2 (TP0.1M + 1) matrix was imported into Seurat. The top 1,000 highly variable genes were identified, and we retained genes detected in at least five samples. Cell cycle-related genes and index switching effect-related genes (including genes highly expressed in non-β-cells, such as the exocrine enzymes Cel and Ctrb1) were identified using hierarchical clustering and were excluded from further analyses. After PCA, we observed that the major variation between the control and STZ-treated β-cells was represented on PC1. Therefore, we performed hierarchical clustering on the top 200 high-loading genes of PC1, and divided the β-cells into three subpopulations. We performed differential expression analyses (group 1 vs group 2/3 and group 2 vs group 3) using FindMarkers in Seurat with the parameter “logfc.threshold = 0.6”. To reduce noise, we only retained differentially expressed genes co-expressed with at least 10 other differentially expressed genes (ρp > 0.2) [25]. We also performed hierarchical clustering on the retained differentially expressed genes to identify and exclude index switching effect-related genes. Gene ontology (GO) enrichment analyses were conducted using GOstats (v2.50.0) [26]. To evaluate β-cell state transitions, we performed PCA on the β-cells with filtered differentially expressed genes, and fit the principal curve on the PCA plot using princurve (v2.1.4) [27]. The lambda values from the princurve results were extracted and used as the pseudo-time values.

Comparison of immature β-cells and Glut2low cells

The P3, P12, and P21 β-cell identification pipeline was similar to the pipeline described in Section 2.10. PCA, UMAP, and forced-directed layout (FDL) were performed on P3, P12, P21, and adult β-cells with the top 2,000 highly variable genes (excluding cell cycle-related genes and index switching effect-related genes). To define subgroups in the P3, P12, and P21 β-cells, we performed hierarchical clustering on P3, P12, P21, and adult β-cells with genes differentially expressed among three groups of adult β-cells (Figure 2D). If a P3, P12, or P21 β-cell clustered with group 2 adult β-cells, we considered this β-cell to be a group 2 β-cell. To identify the maturation-related genes, we performed differential expression analysis between P3 and adult group 1 β-cells and filtered genes using the same criteria described in Section 2.11. We performed differential expression analysis between the P3 and adult group 1 β-cells using FindMarkers in Seurat with the parameter logfc.threshold = 0.6 and only retained genes co-expressed with at least 10 other differentially expressed genes (ρp > 0.2). Then we performed hierarchical clustering on the retained genes to identify and exclude index switching effect-related genes. The overlap between the maturation-related genes and STZ-related genes (genes in Figure 2D) is exhibited in the Venn diagram (Figure 3E).

Figure 2

β-cell heterogeneity.A and B: PCA and UMAP plots of the β-cells. Each dot represents a single cell sample. The colors denote the group (left), treatment condition (middle), and cell source (right). The circled numbers indicate the cell source labeled in Figure 1A. C: Bar plot showing the proportion of the three cell groups in the β-cells from the control and STZ-treated mice. D: Heatmap showing the expression patterns of genes differentially expressed among the three groups. Each column represents a single cell sample, and each row represents one gene. The gene count of each gene cluster is labeled on the left of the heatmap. Transcription factors (TFs) included in each gene cluster are labeled on the right of the heatmap. E: Selected GO terms enriched in the three clusters in (D). F: Expression levels of differently expressed genes are shown on the PCA plot (top) and violin plot (middle). The black line within each violin plot indicates the median of expression level. The expression levels of the indicated genes were verified via bulk cell RT-qPCR and normalized to β-actin (Actb) expression (bottom). Data are shown as means ± SEM. Student t-test: ∗ p-value < 0.05; ∗∗ p-value < 0.01; ∗∗∗ p-value < 0.001. G: Immunofluorescent staining of insulin and NKX6.1 in frozen sections from the control (top) and STZ-treated D1 (bottom) pancreatic tissues. The arrowheads indicate NKX6.1low β-cells in the STZ pancreas. Scale bars: 20 μm. H: The violin plots showing the pseudo-time value for each group and cell source. The circled numbers indicate the cell source labeled in Figure 1A. The black line within each violin plot indicates the median of the pseudo-time values.

Figure 3

cells represent a status distinct from immature β-cells.A, B, and C: PCA (A), UMAP (B), and FDL (C) plots of all of the Ngn3-Cre; Rosa-RFP+ β-cells. Each dot represents a single cell sample. The cells are colored by cell stage, cell type, treatment condition, and expression level of Glut2, respectively. Cell counts are indicated in brackets. The black and red curves indicate maturation and STZ-related trajectory, respectively. D: Immunofluorescent staining of insulin and GLUT2 from P3, P12, P21, and P60 pancreatic tissues. The arrows indicate insulin+GLUT2low β-cells. Scale bars: 20 μm. E: Venn plot showing the overlap between the maturation-related genes and the genes differentially expressed among groups 1–3 in Figure 2D. F–H: Heatmap showing the expression patterns of 865 maturation exclusively related genes (F), 134 maturation and STZ-related genes (G), and 317 STZ exclusively related genes (H) in Figure 3E among the Ngn3-Cre; Rosa-RFP+ β-cells. Each column represents a single cell sample, and each row represents one gene. The gene count of each gene cluster is labeled on the left of the heatmap.

Figure 1

Identification of cell types.A: Overview of 2,400 Ngn3-Cre; Rosa-RFP+ cells analyzed by scRNA-seq in this study. The cell # column shows the cell counts from the indicated treatment conditions. D (M): days (months) after STZ-treatment. Insulin was subcutaneously implanted at D2 after STZ treatment. N: no insulin implant. ∗:no insulin implant during M7-M9. B–C: The UMAP and t-SNE plots show the distribution of cells with different types and treatment conditions. Each dot represents a single cell sample. Cell counts are indicated in brackets. The circled numbers indicate the cell source labeled in (A). D: Expression levels of marker genes in each cell type. The point color from blue to red represents the expression level from low to high. The point size represents the detected ratio of the marker gene in the corresponding cell type.

For FDL analysis, we used FindNeighbors in the Seurat package to generate the shared nearest neighbor (SNN) matrix based on the PCA results. Then we built a graph using the graph.adjacency function in the igraph package (v1.2.4.1) [28] according to the SNN matrix. Finally, we generated the FDL results using the layout_with_fr function in the igraph package.

Heterogeneity analyses of α/δ/PP-cells

We applied the same pipeline in a β-cell heterogeneity analysis to explore the heterogeneity of α-cells, δ-cells, and PP-cells, respectively. We selected the top 1,000 highly variable genes, filtered out genes detected in less than 5 samples, and excluded cell cycle-related and index switching effect-related genes. The retained genes were used for PCA, UMAP, and hierarchical clustering. To evaluate whether the subgroups also existed in the α-cells, δ-cells, and PP-cells, we calculated the average silhouette widths [29] using the silhouette function in the cluster package (v2.0.5) [30] based on the hierarchical clustering. A relatively greater average silhouette width indicated an unambiguous clustering result. If a peak of the average silhouette width on N clusters was observed, the cells were expected to cluster into N subpopulations. In contrast, if no distinct peak of average silhouette width was observed, the cells were expected to be homogeneous.

Results

scRNA-seq of the STZ-treated pancreatic endocrine cells

To explore the transcriptomic states of individual pancreatic endocrine cells in the STZ-treated mice, we isolated cells from the STZ-treated mice via FACS and performedscRNA-seq using a modified STRT-seq (mSTRT-seq) approach [17,18]. We used the Ngn3-Cre; Rosa-RFP mice to enrich RFP+ pancreatic endocrine cells (Figure 1A–D). The mice were treated with a single high dose (200 mg/kg) of STZ and were assessed from day 6 (D6) to month 9 (M9) (Figure 1A). At D1 after STZ treatment, we observed an impaired islet structure, with many INS+ cells diminished (Supplementary Figure S1A). At D3 and D6, body weight significantly decreased, and the treated mice exhibited severe hyperglycemia and significantly reduced insulin concentrations compared with the control groups (Supplementary Figure S1B). These results indicated that the β-cells were significantly damaged in our STZ model. In this study, only mice with hyperglycemia at D2 after STZ treatment were used for subsequent experiments (Supplementary Figure S1C). To prevent the mice from dying from extreme hyperglycemia and study the cell states over a longer time period, we subcutaneously implanted an insulin capsule at D2 after STZ treatment [13]. We measured non-fasting blood glucose at 4 p.m. every three days and implanted insulin pellets when hyperglycemia (>400 mg/dL) was detected. The blood glucose level remained below 400 mg/dL after insulin implantation, but the mice without implantation remained hyperglycemic (>600 mg/dL) and subsequently died (Supplementary Figure S1D). However, two mice that did not receive an insulin supply but survived at M2 after STZ treatment were included in this study (Figure 1A). At D12, M5, and M9 after STZ treatment, we collected the pancreatic tissues to examine the morphology of the islets and observed few β-cells in the damaged islets (Supplementary Figure S1E and F).

In total, 2,400 Ngn3-Cre; Rosa-RFP+ mSTRT-seq samples passed quality control. On average, 9 × 104 unique molecular identifiers (UMIs) and 4,500 genes were detected in each single cell (Supplementary Figure S2A and B and Supplementary Table 1). Using UMAP [31] and t-SNE with Louvain cell clustering algorithms [24], we divided the cells into seven cell types based on their marker gene expressions: β-cells (Ins1+ and Ins2+), α-cells (Gcg+), δ-cells (Sst+), PP-cells (Ppy+), ductal cells (Spp1+), acinar cells (Cel+), and immune cells (Ptprc+) (Figure 1B–D and Supplementary Figure S2C). Ngn3 is also expressed in a fraction of the multipotent pancreatic progenitor cells and trunk cells (progenitors of ductal cells and endocrine progenitors) [32], which can also produce an exocrine population. These non-endocrine cells could be sorted from the single-cell suspension of the STZ-treated pancreas. For this reason, we obtained many non-endocrine cells from the Ngn3-Cre; Rosa-RFP mice. Notably, we collected at least 150 cells for scRNA-seq from each of the STZ-treated animals with an insulin supply. Because the cells from a range of time points were clustered together on the UMAP and t-SNE plots and did not show variations at the transcriptomic level (Figure 1B and C), we concluded that the cell states after STZ treatment were stable over time and that inter-individual variations of the transcriptomic profiles between different time points could be ignored in this analysis.

β-cell heterogeneity in the STZ model

Next we focused our analyses on β-cells. PCA and UMAP revealed three β-cell groups (Figure 2A and B). Group 1 included cells from the non-STZ-treated mice, group 2 contained cells from both the STZ-treated and non-STZ-treated mice, and group 3 exclusively included cells from the STZ-treated mice (Figure 2C). We identified 451 differentially expressed genes, including 23 transcription factors (TFs), among these three β-cell groups (Figure 2D and Supplementary Table 2). Cluster a genes were downregulated in group 2/3 cells and included key TFs important for maintaining β-cell identity, such as Nkx6.1, Pdx1, and Vdr [[33], [34], [35]], and genes associated with β-cell function, such as Glut2 (Slc2a2), Ucn3, and Ins1 (Figure 2F). Gene ontology (GO) enrichment analysis confirmed that cluster a genes were enriched for terms related to β-cell functions (Figure 2E and Supplementary Table 3). Cluster b genes were upregulated in the group 2/3 cells and enriched for the GO term of negative regulation of response to external stimulus (Figure 2E and Supplementary Table 4). The GO analysis also showed that fructose metabolism was upregulated in the group 3 β-cells (Figure 2E and Supplementary Table 5). Several genes, such as Fbp1, Fbp2, and Pdk1, have been reported to participate in fructose metabolism, resulting in reduced insulin secretion (Supplementary Tables 2 and 5) [[36], [37], [38]]. Our scRNA-seq analysis revealed that these genes were upregulated in the group 3 β-cells, suggesting that the group 3 β-cells were in a dysfunctional state regarding secreting insulin to regulate glucose homeostasis. Etv1 and Rorc, negative regulators of insulin secretion [39,40], were included in cluster b/c (Figure 2D), consistent with the lower expression of Ins1 in group 2/3 (Figure 2F), and suggesting weak insulin secretion in the group 2/3 cells. We selectively validated the downregulation of Glut2, Ucn3, Ins1, Nkx6.1, Pdx1, and Vdr and the upregulation of Etv1, Klf10, and Tspan8 in the sorted STZ-treated β-cells from Ins1-RFP transgenic mice (D12) [50] compared to the control using RT-qPCR (Figure 2F) and also confirmed the reduced expression of NKX6.1 in the STZ-treated mice via immunofluorescence staining (Figure 2G).

To analyze whether the group 2/3 β-cells could be restored to functional group 1 β-cells, we aligned the cells according to their positions along the group 1 to group 3 trajectory (Figure 2H). We found that the group 1 β-cells were quickly eliminated after STZ treatment, and the group 3 β-cells appeared over time (Figure 2H). However, only a few group 1 cells were detected in the STZ pancreas, independent of the presence of an insulin implant (Figure 2H and Supplementary Figure S3A). This finding indicated that the group 2 cells transformed into group 3 β-cells but did not restore their group 1 β-cell identity over a long period of time following STZ treatment. Taken together, our analyses revealed the heterogeneity of both the normal (group 1 and group 2 cells) and STZ-treated β-cells (group 2 and group 3 cells). In the STZ-treated animals, the group 1 β-cells were selectively killed, whereas the group 2 β-cells escaped STZ-mediated cell death but correspondingly altered their gene expression and morphed into a pathological state represented as group 3 β-cells.

It has been reported that a small fraction of β-cells locates to the islet edge and lacks expression of Glut2 and Ucn3 [41]. We re-analyzed the transcriptional profiles of the Ucn3- (Glut2low) and Ucn3+ (Glut2high) β-cells from this previous study [41]. We found that the 229 genes highly expressed in the group 1 cells and 163 genes highly expressed in the group 2 β-cells (no STZ treatment) were also generally highly expressed in the Ucn3+ (Glut2high) and Ucn3- (Glut2low) cells, respectively (Figure 2D and Supplementary Figure S3B and C). Therefore, based on this global similarity of gene expression patterns, we presumed that the group 2 β-cells included “virgin” β-cells identified by van der Meulen et al. [41].

Glut2low β-cells in the adult mice represented a state distinct from immature β-cells

Glut2low β-cells in the adult mice were considered immature β-cells because they displayed a low expression of the maturation marker gene Ucn3 (Figure 2F) [41,42]. To evaluate the maturation status of these Glut2low β-cells, we performed mSTRT-seq on the Ngn3-Cre; Rosa-RFP cells from the postnatal day 3 (P3), P12, and P21 mice (Supplementary Figure S4). Combining these datasets with adult Ngn3-Cre; Rosa-RFP datasets, we performed PCA, UMAP, and FDL analyses (Figure 3A–C). On the PCA plot, most of the β-cells were arranged along a path from P3 to adult Glut2high that reflected the β-cell maturation process (Figure 3A). Curiously, we observed a few Glut2low cells from the P3, P12, and P21 mice (Figure 3A), but at P12 and P21, not P3, the Glut2low cells were separated from the Glut2high cells on the developmental trajectory (Figure 3A). By performing immunostaining, we observed GLUT2low β-cells at the islet periphery at various developmental stages (Figure 3D). The distribution pattern of the GLUT2low β-cells was consistent with Van der Meulen's findings [41], which further indicated that the Glut2low β-cells (group 2) in the healthy mice identified by scRNA-seq were “virgin” β-cells. These findings clearly showed that Glut2low cells developed during the β-cell maturation process. However, STZ-treated adult β-cells (consisting of group 2/3 cells) were not located along this maturation path (as indicated by the black curve in Figure 3A). The UMAP and FDL analyses showed a similar result (Figure 3B and C). Moreover, the genes that were variably expressed during the maturation process displayed limited overlap with the genes variably expressed between the Glut2high and Glut2low β-cells from our STZ study (Figure 3E–H and Supplementary Table 6). The overlapping portion of the Venn diagram, including genes such as Ucn3 and Glut2 (Supplementary Table 6), was related to both β-cell maturation and β-1/2/3 cell transition processes.

The expression profiles of the non-β endocrine cells were not affected in the STZ model

To investigate whether non-β endocrine cells were affected in the STZ model, we extended our analyses of the α-cells, δ-cells, and PP-cells. Although a small fraction of cells expressed at least two of the Ins1/Ins2, Gcg, Sst, and Ppy genes, these multi-hormone expressing cells were defined as a specific endocrine cell type based on their transcriptomic profiles (Figure 1B and C). We found that α-cells, δ-cells, and PP-cells from the control and STZ-treated mice overlapped on the PCA and UMAP plots and did not form distinct clusters associated with STZ treatment (Figure 4A and B). Additionally, hierarchical clustering based on highly variable genes could not distinguish non-β endocrine cells between the control and STZ-treated mice, indicating that the transcriptional statuses of the α-cells, δ-cells, and PP-cells were not affected by STZ treatment, unlike the β-cells (Figure 4C). Furthermore, the silhouette analysis [29] (Materials and Methods) suggested that in the control and STZ-treated mice, the α-cells, δ-cells, and PP-cells were homogeneous, whereas the β-cells from the non-treated control mice by themselves were heterogeneous (Supplementary Figure S5A and B), and after a high dose of STZ, the β-cells existed as three subpopulations (Figure 4D). Therefore, we concluded that STZ treatment did not affect the transcriptomic profiles of the non-β endocrine cells.

Figure 4

Transcriptomic profile of non-β endocrine cells in the STZ model.A and B: PCA (A) and UMAP (B) of α-, δ-, and PP-cells from the Ngn3-Cre; Rosa-RFP pancreas. Each dot represents a single cell sample. The colors denote the treatment condition. C: Hierarchical clustering of highly variable genes on the α-, δ-, PP-, and β-cells. Each column represents a single cell sample, and each row represents one gene. D: Silhouette analyses of the hierarchical clustering of the α-, δ, PP-, and β-cells in (C). A relatively higher average silhouette width indicates an unambiguous clustering result. If a peak of the average silhouette width on N clusters is observed, the cells are expected to be clustered into N subpopulations. In contrast, if no distinct peak of the average silhouette width is observed, the cells are expected to be homogeneous. E-G: Three dimensions (3D) PCA (E), UMAP (F), and FDL (G) of the endocrine cells from the Ngn3-Cre; Rosa-RFP pancreas. Each dot represents a single cell. The colors denote cell type and treatment condition. Expression levels of Ins1, Ins2, Gcg, Sst, Ppy, and Glut2 are shown on the right.

It has been reported that, under normal physiological conditions, Glut2lowUcn3low “virgin” β-cells represent an intermediate stage in α to β trans-differentiation [41]. We therefore performed PCA, UMAP, and FDL analyses on endocrine lineages (Figure 4E–G and Supplementary Figure S5C and D). However, we did not detect a clear STZ-induced trans-differentiation trajectory from the α-cells to β-cells on the PCA, UMAP, and FDL plots (Figure 4E–G and Supplementary Figure S5C and D). Moreover, we performed hierarchical clustering analysis between the group 1 β-cells (Glut2high), group 2 β-cells (Glut2low), and the α-cells, but did not find that the group 2 β-cells were more similar to the α-cells at the transcriptomic level (Supplementary Figure S5E and F). These analyses suggest that trans-differentiation from the α-cells to Glut2low β-cells may not be a dominant mechanism in the maintenance of the small population of β-cells after STZ treatment.

Discussion

The STZ model is commonly used to study type-1 and type-2 diabetes [10,11]. However, precise characterization of the β-cell states in the STZ model has been challenging due to the rarity of β-cells and the model's impaired islet structure. In this study, we performed scRNA-seq on endocrine cells from reporter mice over a long period of time (D6 to M9) after STZ treatment. We found that the adult pancreas included two groups of β-cells that differed in their Glut2 expression (Glut2high and Glut2low). In our previous single-cell study, we detected several Glut2low β-cells [43]. However, due to a low number of sequenced cells, this small population of Glut2low β-cells could not be identified as a cell cluster in that study. We also confirmed the existence of Glut2low Ucn3low β-cells under normal physiological conditions, a population of cells previously shown to locate to the edge of the islets [41]. Although some genes related to β-cell functions were downregulated in the Glut2low Ucn3low cells, our transcriptomic analysis suggests that the Glut2low Ucn3low β-cells induced by STZ did not reflect immature β-cells, as they differentially expressed ∼300 genes that were not related to maturation (Figure 3E). Glut2 mediates STZ membrane transport, as a consequence, most of the Glut2high β-cells were damaged by STZ, but the Glut2low cells and possibly a small fraction of the Glut2high β-cells survived. However, these surviving cells gradually altered their states after STZ treatment despite exogenous insulin application. This result may suggest that normal functional β-cells or normal islet structures are necessary to maintain the status of Glut2low cells. In this study, we did not detect restoration of the Glut2low to Glut2high cells even 5–9 months after STZ-treatment.

Non-β endocrine cells have been reported to convert into insulin-secreting β-like cells under certain conditions [2,4,41,44,45]. Genetic inactivation of Dnmt1 and Arx in adult mice [44] or Pdx1-enforced expression in endocrine progenitors [45] could convert pancreatic α-cells into β-cells. Thorel et al. found that α-cells could transform into β-cells after near-total β-cell ablation in a diphtheria-toxin-induced model [2]. During trans-differentiation, α-cells gradually upregulated β-cell markers, such as Pdx1 and Nkx6.1, and eventually expressed insulin. Van der Meulen et al. also observed α- to β-cell trans-differentiation via a “virgin” β-cell state under normal conditions [41]. Moreover, an anti-malarial drug artemisinins was shown to drive the in vivo α-to β-cell fate transition by enhancing GABA signaling [46], although the main conclusion of that study has been challenged by others [47,48]. These findings encouraged us to investigate the process of cell fate trans-differentiation in the STZ model. If trans-differentiation from the non-β endocrine cells to the β-cells occurs in the STZ model, we would expect to observe heterogeneity in the non-β endocrine cell populations, which may reflect their state changes and reveal a defined cell fate transition trajectory via scRNA-seq analysis. However, we did not detect heterogeneity in the non-β endocrine cell populations or the previously described cell fate transition trajectory in this study (Figure 4). We acknowledge that the number of cells analyzed in this study may not be sufficient to detect cell trans-differentiation events in a small proportion of cells. A larger scale single-cell analysis combined with a lineage-tracing approach is required to investigate whether trans-differentiation events occur in the STZ model.

Many genes related to β-cell functions and diabetes development were differentially expressed in the STZ-treated β-cells compared to the functional Glut2high β-cells, indicating group 2/3 (Glut2low) cells were not fully functional in responding to glucose and insulin secretion. Group 3 cell upregulated genes Fbp1 and Fbp2 are primarily expressed in gluconeogenic organs, such as liver and muscle. Fbp1 negatively regulates insulin secretion, and Fbp2 is proposed to regulate indirect glycogen synthesis from lactate. Both genes were found to be upregulated in diabetic animal models [36] and diabetic subjects [38]. FBP inhibitors have been designed to improve insulin secretion, and several inhibitors, such as MB05032, have been reported to increase insulin secretion in the MIN6 β-cell line and in isolated mouse islets [49]. The transcriptomic signatures we observed in the group 3 cells may reveal other therapeutic targets to improve β-cell function. In addition, the group 2/3 (Glut2low) cells may represent a platform for screening factors and drugs that can reverse β-cell dysfunction via restoration of Glut2low to Glut2high cells as well as to study the molecular mechanisms related to this process. Our single-cell analyses have deepened our understanding of β-cell heterogeneity, the transcriptomic profiles of β-cell subsets, and the applicability of the STZ diabetes model.

Data availability

The RNA-seq data were submitted to the GEO repository under accession number GSE137909.

Funding

This study was supported by the National Key Research and Development Program of China (2019YFA0801500), the Ministry of Science and Technology of the People's Republic of China (973 Program 2015CB942800), the National Natural Science Foundation of China (31521004, 91753138), funding to C.-R.X., and funding provided by Boehringer Ingelheim Pharma GmbH.

Author contributions

C.-R.X. conceived the project. W.Z., J.Y., J.D.H., and C.-R.X. designed the study. Y.F., X.-X.Y., Y.Z., M.-Y.H., L.Y., and L.-C.L. conducted the research. W.-L.Q., W.Z., M.F., J.Y., J.D.H., and C.-R.X. analyzed the data. Y.F., W.-L.Q., J.D.H., and C.-R.X. wrote the manuscript. C.-R.X. is the guarantor of this work and, as such, had full access to the study data and takes responsibility for the integrity of the data and the accuracy of the data analysis.