Single-Cell Transcriptome Profiling of Pancreatic Islets From Early Diabetic Mice Identifies Anxa10 for Ca²⁺ Allostasis Toward β-Cell Failure

The islets of Langerhans are the pancreatic tissues containing several distinct endocrine cell types, including α-, β-, δ-, pancreatic polypeptide (PP), and ε-cells, in addition to several important nonendocrine cell types, including pancreatic stellate cells, islet resident macrophages, sympathetic and parasympathetic nerves, and vascular cells, such as endothelial cells and pericytes (1). Pancreatic islets have been much investigated because they are critical for glucose homeostasis regulation and are key players in diabetes development.

Type 2 diabetes is triggered by a combination of elevated peripheral insulin resistance and progressive decline in β-cell function and mass (2). The continually rising number of patients with type 2 diabetes is a worldwide health challenge, and the need for early prediction and effective treatment for this disease has become increasingly significant. Several mechanisms have been proposed to elucidate the functional β-cell mass decline in type 2 diabetes, including glucolipotoxicity (3), oxidative stress (4), endoplasmic reticulum (ER) stress (5), inflammation (6), apoptosis, and loss of β-cell identity and function (7). Nevertheless, owing to the cellular heterogeneity of the pancreatic islets, the precise cellular and molecular contexts in the deterioration of insulin secretory function are still unidentified.

Single-cell RNA sequencing (scRNA-seq) is an unbiased and powerful procedure for characterizing different cell types in complex tissues from the perspective of health and disease. scRNA-seq investigations of pancreatic islets have been noted in humans (8–10) and animal models (11–13), highlighting the functional heterogeneity of β-cells, postnatal maturation and dedifferentiation of β-cells, and characterization of specific islet cells in specific conditions. Still, it is uncertain how dysfunctional β-cells advance into type 2 diabetes in the roadmap from healthy to diabetic islets in an obese condition. Moreover, the driver genes have not been recognized in a prediabetic state. Here, we used a mouse model of obesity and diabetes induced by leptin receptor deficiency and conducted robust scRNA-seq to capture and profile cell type–specific gene expression of pancreatic islets in nondiabetic, prediabetic, and diabetic states. We observed a gene related to Ca²⁺ signaling with remarkable correlations to glucolipotoxicity in a prediabetic state and found cellular trajectories and interactions that may influence the transition from a healthy to a dysfunctional β-cell.

All animal husbandry and experiments conformed to the guidelines of the University of Tsukuba’s regulations of animal experiments and were approved by the animal experiment committee of the University of Tsukuba. All experiments were executed following the Animal Research: Reporting of In Vivo Experiments guidelines. We purchased male BKS.Cg-Lepr^db/Lepr^db (db/db) and BKS.Cg-Dock7^m/Dock7^m (^+/+) mice from Charles River Laboratories Japan. All animals were housed in a pathogen-free barrier facility with a 12-h light/dark cycle and were provided free access to normal chow (MF diet; Oriental Yeast Co., Ltd., Tokyo, Japan) and water. Age- and sex-matched littermates were used for all experiments. Mice were sacrificed during the light phase after food deprivation for 4 h. Plasma glucose, insulin, and glycosylated hemoglobin A_1c (HbA_1c) levels were measured as described previously (14).

Seven-week-old male C57BL/6J mice were obtained from CLEA Japan (Tokyo, Japan) and adapted to the environment for 1 week before the experiment. Diabetes was stimulated by intraperitoneal injection of streptozotocin (STZ) 100 mg/kg (Sigma-Aldrich, St. Louis, MO) two times every other day. One week after the last injection, blood glucose levels were checked, and mice with high glucose (≥350 mg/dL) were used in the experiments.

Four-week-old male db/db mice were categorized into three groups and fed the normal chow or a chow diet containing 0.03% empagliflozin (Boehringer Ingelheim, Ingelheim am Rhein, Germany) or verapamil 1 mg/mL (Nacalai Tesque, Kyoto, Japan) in their drinking water. Verapamil and empagliflozin treatments were performed for 2 weeks and 1 week, respectively.

Isolation of mouse islets was performed following the hand-picking protocol described previously (14), with some modifications. Briefly, after clamping the common bile duct at a point close to the duodenum outlet, 2.5 mL Krebs-Ringer bicarbonate HEPES buffer (pH 7.4) comprising 0.5% BSA and 0.4 mg/mL Liberase TL Research Grade (Roche, Basel, Switzerland) were injected into the duct. The swollen pancreas was extracted and incubated at 37°C for 20 min. The pancreas was then dispersed by pipetting, and after washing twice with Krebs-Ringer bicarbonate HEPES buffer, the islets were obtained manually under a stereomicroscope.

For the scRNA-seq, islets from three to four mice per group were pooled and dissociated into single-cell suspension with accutase (Nacalai Tesque) at 37°C for 20 min. Live single cells were isolated by Percoll density gradient separation (25%/65%) at 4°C, and the live-cell suspension was immediately used for the scRNA-seq library using terminator-assisted solid-phase cDNA amplification and sequencing (15). Briefly, reverse-transcribed, exonuclease I–treated BD Rhapsody beads underwent terminator-assisted terminal deoxynucleotidyl transferase reaction, second-strand synthesis reaction, and a first and second round of whole-transcriptome amplification reaction. Size distribution and concentration of amplified cDNA libraries were investigated by using an Agilent Bioanalyzer 2100 system and an Agilent High-Sensitivity DNA kit. Illumina libraries were built from 100-ng amplified cDNA libraries by using NEBNext Ultra II FS DNA Library Prep Kit for Illumina (New England Biolabs). Illumina adapters and unique dual barcodes were included in the hashtag libraries by PCR. Size distribution and concentration of amplified Illumina libraries were examined by using the MultiNA system and KAPA Library Quantification Kit (KAPA Biosystems). Sequencing was conducted with a NovaSeq 6000 sequencer (Illumina, San Diego, CA) following the manufacturer’s instructions. Pooled library concentration was adjusted to 1.75 nmol/L, and 12% PhiX control library (PhiX Sequencing Control V3; Illumina) was spiked into the library.

FASTQ data preprocessing and generation of the single-cell gene expression matrix were performed as described previously (15). To predict background and signal read count distribution, we used the Gaussian mixture model that was previously used to predict the gene expression distribution of scRNA-seq data sets (15,16). Raw sequencing data have been deposited in the Data Bank of Japan database (DRA012714).

The clustering of single cells of each data set was conducted as described previously by using Seurat version 2.3.4 in R version 3.6.3 (17). Each identified cluster was manually annotated by their marker genes that were previously reported as cell subset–defining marker genes. Supplementary Table 1 displays all identified marker genes and cell cluster annotations. Pathway analysis against each marker gene was performed using Ingenuity Pathway Analysis (IPA) (QIAGEN Bioinformatics, Redwood City, CA) to reveal the biological meaning of each marker gene. The technique of constructing gene correlation networks was proposed by Asano et al. (18). We applied several ranking procedures previously proposed for ranking Web pages (PageRank and personalized PageRank) to the networks. We also visualized gene expression and clustering outcomes in a Uniform Manifold Approximation and Projection space using Seurat with default parameters.

Cell-cell interaction among pancreatic islets according to the scRNA-seq data set was performed using the CellChat 1.1.3 package with consideration of population size (population.size = TRUE and raw.use = TRUE in computeCommunProb function) (19). Statistically significant interactions and pathways were chosen by computeCommunProbPathway and aggregate Net functions (thresh = 0.01). Circle plot visualization of the cell-cell interaction network and the strength of the communication between each cell subset were conducted using the netVisual_circle function in CellChat. Outgoing/incoming signaling strength was computed and visualized using the netAnalysis_signalingRole_heatmap function in CellChat.

Isolated islets were incubated overnight at 37°C in RPMI medium (Gibco) comprising 11.1 mmol/L glucose and 2.06 mmol/L l-glutamine and supplemented with 10% FBS, 10 mmol/L HEPES, 1 mmol/L sodium pyruvate, 100 μg/mL streptomycin, and 100 IU/mL penicillin before experiments. Insulin secretion and insulin content of islets were measured as previously described (14).

Recombinant adenovirus expressing full-length mouse annexin A10 (Anxa10) was prepared as detailed elsewhere (20). cDNAs encoding full-length mouse Anxa10 were cloned into the Gateway entry vector pENTR/D-TOPO (Thermo Fisher Scientific, Waltham, MA), and the adenoviral construct was generated by homologous recombination between the entry vector and the pAd/CMV/V5-DEST vector (Thermo Fisher Scientific). Adenoviruses encoding shRNA targeting LacZ and Anxa10 genes were constructed by cloning synthetic DNA into the pENTR/U6 entry vector followed by homologous recombination with the pAd promoterless vector (Thermo Fisher Scientific). The target sequences were as follows: 5′-GCGTGCTGTATTCACACTGAA-3′ for Anxa10 #1, 5′-CCACCATCTTATGATGCTCAT-3′ for Anxa10 #2, and 5′-CTACACAAATCAGCGATTT-3′ for LacZ. Recombinant adenoviruses were propagated in HEK293 cells and purified by CsCl gradient centrifugation as described previously. Isolated Mouse islets were infected with adenoviruses for 2 h in serum-free media and examined at 48–72 h postinfection.

cDNA-encoding mouse Anxa10 was cloned into pcDNA6.2/N-EmGFP plasmid (Invitrogen). MIN6 cells were transfected with pcDNA6.2/N-EmGFP-DEST or pcDNA6.2/N-EmGFP-Anxa10 using Lipofectamine 3000 (Life Technologies) according to the manufacturer’s instructions. At 48 h posttransfection, stably transfected MIN6 cell selection was performed using a selection medium comprising 10 μg/mL puromycin for 2 weeks.

MIN6 cells stably transfected by either pcDNA6.2/N-EmGFP-Anxa10 or pcDNA6.2/N-EmGFP were plated onto a collagen-coated 96-well, black-walled, clear-bottomed plate (Perkin-Elmer Ltd., Beaconsfield, U.K.) and allowed to grow for 24 h. Cells were then loaded with 3-μmol/L Indo-1, AM (AAT Bioquest, Sunnyvale, CA) in a recording medium comprising 2.8 mmol/L glucose, 115 mmol/L NaCl, 5.4 mmol/L KCl, 1.8 mmol/L CaCl₂, 0.8 mmol/L MgCl₂, 20 mmol/L HEPES, 0.04% Pluronic F-127 (Sigma-Aldrich), and 1.25 mmol/L Probenecid (Nacalai Tesque) for 30 min at 37°C. Following the incubation, cells were washed with PBS (-). The medium was then replaced to the new recording medium. Cells were treated with 20 mmol/L glucose, 30 mmol/L KCl, or 1 μmol/L thapsigargin and then chelated with intracellular Ca²⁺ ([Ca²⁺]i) and extracellular Ca²⁺ with 10 mmol/L EGTA with the exception of the thapsigargin treatment. Thapsigargin stimulation was performed in the presence or absence of extracellular Ca²⁺. The fluorescence of Indo-1 was acquired with the multimode plate reader Varioskan LUX (Thermo Fisher Scientific) at 405 and 485 nm with excitation set at 330 nm.

Pancreatic cancer specimens from 17 patients who underwent resection with curative intent from 2017 to 2021 were acquired from the University of Tsukuba Hospital under protocols approved by the ethics committee of the University of Tsukuba (approval no. R02-276). Written informed consent was collected from all patients. Patients receiving chemotherapy and radiation therapy before surgery were not included. Tumor stages were assigned following the Union for International Cancer Control TNM staging. The clinical characteristics of the analyzed patients are shown in Supplementary Table 5.

Immunofluorescence staining using mouse and human pancreatic sections was performed as described previously (14). Specific antibodies used for immunofluorescence were ANXA10 (ab213656; Abcam, Cambridge, MA), ALDH1A3 (1:200 dilution, NBP2-15339; Novus Biologicals, Littleton, CO), glucagon (sc-514592; Santa Cruz Biotechnology, San Francisco, CA), amylase (sc-46657; Santa Cruz Biotechnology), and chymotrypsinogen B (Ctrb) 1/2 (sc-398721; Santa Cruz Biotechnology). Images were obtained using a Leica TCS SP8 confocal laser scanning microscope (Leica Microsystems, Wetzlar, Germany).

Values are expressed as mean ± SEM. Student t tests were used to compare the means of the two groups. One-way or two-way ANOVA was used for multiple group mean comparisons, with Dunnett and Šidák post hoc tests to analyze differences between specific groups. All analyses were performed using GraphPad Prism 7.0 software (GraphPad Software, San Diego, CA).

Data and reagents generated in the current study are available from the corresponding author upon reasonable request.

To investigate cellular heterogeneity and dynamic changes during diabetes progression, we conducted scRNA-seq of pancreatic islet cells isolated from 6-week-old nondiabetic db^+/+ and 6- and 10-week-old db/db mice (Fig. 1A). Body weights, fasted blood glucose levels, fed plasma insulin levels, HbA_1c levels (Supplementary Fig. 1A), and insulin staining of islets (Supplementary Fig. 1B) suggest that the 6- and 10-week-old db/db mice are prediabetic and diabetic mouse models, respectively. We performed bulk RNA-seq of pancreatic islets isolated from 6- and 10-week-old db^+/+ and db/db mice and validated that the 6- and 10-week-old db^+/+ islets had similar gene expression profiles with no differentially expressed genes (DEGs) (Supplementary Fig. 1C–F). Isolated islets were dissociated into single-cell suspension, and by Percoll density gradient centrifugation, viable cells were separated. Then, single-cell cDNA libraries were synthesized using terminator-assisted solid-phase cDNA amplification and sequencing (15). We sequenced 4,956 islet cells at 50,000 reads/cells with quality control metrics of scRNA-seq data (Supplementary Table 1 and Supplementary Fig. 2A–C).

Seurat clustering recognized 20 transcriptionally distinct islet cell clusters (Fig. 1B). Two-dimensional fast Fourier transform–accelerated interpolation-based t-stochastic neighborhood embedding was used to visualize all cell types (20). Following the expression of well-established markers (Supplementary Table 2), we identified cell types including endocrine cells, such as β- (Ins1), α- (Gcg), δ- (Sst), and PP cells (Ppy), in addition to nonendocrine cells, such as endothelial (Plvap), monocyte-derived macrophage (Cd74), stellate (Col1a1), ductal (Krt19), acinar (Prss1), Ki67-positive (Mki67), and unclassified cells (Fig. 1B–D, Supplementary Table 3, and Supplementary Fig. 2D). The β-cell clusters were distributed among three islands: a nondiabetic island containing clusters 6 and 9 from the db^+/+ group; a large prediabetic and diabetic island consisting of clusters 0, 1, 2, 7, and 8 from the db/db group; and an island corresponding to cluster 12 that include β-cells from all groups. Cluster 0 mostly comprises β-cells from 6-week-old db/db mice and is considered prediabetic β-cells. Clusters 1, 2, 7, and 8 primarily consist of β-cells from 10-week-old db/db mice and are considered diabetic β-cells. Furthermore, α-cells were split into three clusters. Clusters 3, 5, and 14 consisted of 10-week-old db/db α-cells, 6-week-old db^+/+ and db/db α-cells, and α-cells from all groups, respectively. A similar clustering pattern was also given by Uniform Manifold Approximation and Projection visualization (Supplementary Fig. 2E and F). The composition of islet cells profiled from each group was similar: β-cells accounted for the major component, and the proportion of β-cells in db/db islets was slightly higher than that of db^+/+ islets (Fig. 1E). Nevertheless, the type and proportion of β-cells were adjusted during diabetes progression (Fig. 1F).

IPA of the DEGs was performed to gain insight into the dynamic changes of β-cells during type 2 diabetes progression (Fig. 2 and Supplementary Tables 3 and 4). Nondiabetic β-cell clusters (clusters 6 and 9) demonstrated significant enrichment of genes associated with insulin secretion and synthesis and nitric oxide signaling (Supplementary Fig. 3A). Prediabetic β-cell cluster (cluster 0) displayed enrichment in ingenuity pathways linked to insulin secretion signaling, unfolded protein response (UPR), oxidative phosphorylation, and mitochondrial dysfunction. The most upregulated DEG in this cluster is Anxa10, a member of the annexin family. To our knowledge, the relationship between Anxa10 and diabetes has not been reported. Among the four diabetic β-cell clusters (clusters 1, 2, 7, and 8), cluster 1 was enriched in ingenuity pathways related to UPR, insulin secretion signaling, and creatine phosphate and proline biosynthesis. Additionally, these β-cells demonstrated the upregulation of aldehyde dehydrogenase 1 family member A3 (Aldh1a3), which has been reported as a dedifferentiation signature gene consisting of β-cells in diabetes (21). IPA pathways in cluster 2 listed UPR and ER stress pathways, as supported by the enrichment of DEGs, including stress-activated nuclear protein transcription regulator 1 (Nupr1) (Supplementary Fig. 3B), serpin family A member 7 (Serpina7), and Aldh1a3. Cluster 7 exhibited DEGs like cholecystokinin (Cck) (Supplementary Fig. 3B), which is upregulated during obesity to enhance β-cell survival (22); death-associated protein like-1 (Dapl1), a suppressor of cell proliferation; and Der1-like domain family member 3 (Derl3), which is upregulated by UPR. Cluster 8 demonstrated enrichment in pathways linked to UPR, oxidative phosphorylation, mitochondrial dysfunction, and ER stress pathways. The DEGs in this cluster included gastrin (Gast) and Cck (Supplementary Fig. 3B), which enhance β-cell proliferation (23), and cocaine- and amphetamine-regulated transcript prepropeptide (Cartpt) (Supplementary Table 3), which is crucial for regulating β-cell function (24).

We also identified transcriptional alterations in type 2 diabetic α-cells (Supplementary Fig. 3C and D and Supplementary Tables 3 and 4). Cluster 5, mainly consisting of nondiabetic α-cells (Fig. 1B), demonstrated substantial enrichment of genes related to translation regulation, including EIF2 and mammalian target of rapamycin (mTOR) signaling. POU class 3 homeobox 4 (Pou3f4), the DEG in cluster 5, encodes a transcription factor that maintains proglucagon gene expression (25). Cluster 14 includes α-cells from nondiabetic and prediabetic groups. N-deacetylase and N-sulfotransferase 4 (Ndst4), encoding the enzyme that facilitates the activation of heparan sulfate sulfation and pancreatic sulfated proteoglycans that control β- and α-cell development (26), were highly expressed in cluster 14. Remarkably, Ndst4 was also listed as one of the DEGs in PP cells (Supplementary Table 3). Cluster 3 primarily consisted of diabetic α-cells. Uroplakin 3A (Upk3a), which was upregulated in α-cells from diet-induced obese mice (27), was highly expressed in type 2 diabetic α-cells. The significantly enriched IPA pathways in this cluster comprised EIF signaling, regulation of eIF4 and p70S6K signaling, and mTOR signaling, indicating the possible participation of this translation initiation pathway in the molecular adaptation of α-cells to diabetic stress.

We determined Anxa10 as a potential biomarker of prediabetic β-cells that might function in the onset and progression of diabetes (Fig. 2 and Supplementary Table 3). We also analyzed mouse transcriptome data sets from the National Center for Biotechnology Information public database and found that Anxa10 expression was increased in islets (Supplementary Fig. 4A) and β-cells (Supplementary Fig. 4B) of db/db mice, as well as in islets from transgenic mice expressing mutant glucokinase (Y214C), which induces β-cell glucotoxicity (Supplementary Fig. 4C). Anxa10 is a member of the annexin family of calcium- and phospholipid-binding proteins (28). Calcium dysregulation can be a central mechanism underlying β-cell dysfunction (29). Additionally, ANXA10 expression is upregulated in various cancers, including pancreatic adenocarcinoma (30).

scRNA-seq showed that Anxa10-positive β-cells were specifically enriched in cluster 0 (Fig. 3A). Conversely, Aldh1a3-positive β-cells were scattered among all db/db β-cells. Consistently, quantitative real-time PCR analysis verified that Anxa10 expression was considerably higher in 6-week-old db/db islets than in 10-week-old db/db islets, while Aldh1a3 expression in db/db islets was consistently elevated during type 2 diabetes progression (Fig. 3B). Immunostaining showed that Anxa10 was robust within islets from 6-week-old db/db mice and gradually declined during 10–12 weeks of age (Fig. 3C). Interestingly, from the inside of the islets, we found the abundance of Anxa10 immunoreactivity was adjusted to the peri-islet in the pancreas of 10–12-week-old db/db mice. Contrarily, Aldh1a3 immunoreactivity was weakly observed in islets from 6-week-old db/db mice and became further elevated during 10–12 weeks of age (Supplementary Fig. 4D). Anxa10 immunoreactivity was also detected in the islets from STZ-induced diabetic mice (Fig. 3D) but not in the islets from ob/ob mice (Supplementary Fig. 4E) and high-fat diet–induced obese mice (Supplementary Fig. 4F). These outcomes indicate that hyperglycemia is crucial for the induction of Anxa10 in the β-cell.

We studied the expression levels of Anxa10 in pancreatic β-cell line MIN6 cells under several culture conditions to identify the inducer of Anxa10 expression in β-cells. The presence of high glucose and high palmitate, mimicking glucolipotoxicity, markedly triggered Anxa10 expression (Fig. 3E). MIN6 cells were also exposed to pharmacological ER stress inducers, indicating that thapsigargin, an inhibitor of the Ca²⁺-ATPase pump that elevates [Ca²⁺]i, dramatically induced Anxa10 expression, whereas tunicamycin, an inhibitor of protein glycosylation, had little influence on Anxa10 expression (Fig. 3E). We treated MIN6 cells with tolbutamide or KCl to determine whether the upregulation of Anxa10 expression is caused by sustainably elevated [Ca²⁺]i concentration and found that elevated Ca²⁺ influx was linked to enhanced Anxa10 expression (Fig. 3E). KCl-induced upregulation of Anxa10 expression was reversed by treatment with the L-type Ca²⁺ channel blocker verapamil. These Anxa10 expression changes were also detected in isolated islets (Fig. 3F). Besides, we treated 4-week-old db/db mice with verapamil or the sodium-glucose cotransporter 2 inhibitor empagliflozin for 2 weeks and observed that inhibition of cytosolic calcium influx by calcium-channel blocker or amelioration of hyperglycemia significantly reduced Anxa10 expression in and around the islets (Fig. 3G and H). Hence, these outcomes show that Anxa10 expression in β-cells is positively controlled by chronic elevation of [Ca²⁺]i concentration stimulated by hyperglycemia, glucolipotoxicity, and ER stress.

To investigate whether ANXA10 expression is also related to human diabetes, multiple transcriptome data sets of pancreatic islets from donors without diabetes and donors with type 2 diabetes were analyzed. In those containing information on this gene, expression levels of ANXA10 in islets were significantly higher in the donors with type 2 diabetes compared with those without diabetes (Fig. 4A). Next, we conducted immunofluorescence staining using human pancreatic sections prepared from samples surgically resected from patients with pancreatic cancer with different diabetic stages. Supplementary Table 5 presents the clinical data. Immunofluorescence analysis showed that acinar cells in nondiabetic and prediabetic groups expressed ANXA10 but less so in the diabetic group (Fig. 4B). In the pancreatic islets, we discovered that the percentage of ANXA10-positive β-cells was considerably higher in the prediabetic group compared with the nondiabetic group (Fig. 4B and C). Contrarily, no significant difference existed in this ratio between the diabetic and nondiabetic groups. This evidence revealed the potential role of ANXA10 in human diabetes pathogenesis.

We conducted the gene correlation network analysis (18) using scRNA-seq data of Anxa10-positive and Anxa10-negative β-cells from 6- and 10-week-old db/db mice to investigate the characteristics of Anxa10-expressing β-cells (Supplementary Fig. 4G and Supplementary Table 6). This evaluation revealed that Anxa10 is positively correlated with oxidative phosphorylation and mitochondrial biogenesis, consistent with the observation that elevation in blood glucose level and mitochondrial glucose metabolism are linked to Anxa10 expression in β-cells. We also studied the subcellular localization of Anxa10 in β-cells and discovered that Anxa10 localizes in the cytosol, nucleus, and microsomal fractions but not in the mitochondria in MIN6 (Supplementary Fig. 4H) and INS-1 832/13 cells (Supplementary Fig. 4I).

To evaluate the role of Anxa10 in β-cells, isolated mouse islets were infected with the Anxa10-expressing adenovirus (Fig. 5A). Anxa10 overexpression did not alter either the basal or the glucose-stimulated insulin secretion (GSIS) but significantly suppressed KCl-stimulated insulin secretion (KSIS) (Fig. 5B) without affecting insulin content (Fig. 5C). Furthermore, we infected MIN6 cells with Anxa10-expressing adenovirus (Supplementary Fig. 5A) and observed that Anxa10 overexpression significantly promoted basal insulin secretion but reduced both GSIS and KSIS (Supplementary Fig. 5B).

We also evaluated the effect of Anxa10 knockdown in islets isolated from 6-week-old db/db mice (Fig. 5D). Knockdown of Anxa10 did not affect GSIS and KSIS but significantly reduced basal insulin secretion in db/db islets (Fig. 5E and F). These results suggest the possibility that elevated basal insulin secretion under hyperglycemic conditions is mediated, at least in part, by altered Anxa10 expression.

Since Ca²⁺ functions significantly in the regulation of insulin release from β-cells (29), Anxa10 might regulate insulin secretion by influencing [Ca²⁺]i levels in β-cells. To evaluate this possibility, we established the MIN6 cell line stably overexpressing Anxa10 and analyzed the effect of Anxa10 overexpression on [Ca²⁺]i concentration changes of MIN6 cells in response to high glucose, KCl, and thapsigargin by applying the Indo-1 Ca²⁺ indicator (Fig. 5G–J and Supplementary Fig. 5C). When the glucose level in the medium increased from 5 mmol/L to 20 mmol/L, a transient steep increase in [Ca²⁺]i followed by a delayed blurred one were detected in control MIN6 cells (Fig. 5G). Overexpression of Anxa10 significantly suppressed high-glucose–induced maximal and delayed [Ca²⁺]i elevation. The increase in [Ca²⁺]i by high K⁺ was considerably attenuated in Anxa10-overexpressing cells compared with control cells (Fig. 5H). In normal extracellular solution, thapsigargin transiently elevated [Ca²⁺]i in both control and Anxa10-overexpressing cells, while the [Ca²⁺]i stimulation by Anxa10 overexpression was significantly less (Fig. 5I). This difference was further increased in Ca²⁺-free extracellular solution: The maximal elevation in [Ca²⁺]i was lower and the Ca²⁺ clearance faster in Anxa10-overexpressing MIN6 cells than in control cells (Fig. 5J). Consistently, the thapsigargin-induced apoptosis rate was substantially lower in Anxa10-overexpressing cells than in control cells (Supplementary Fig. 5D). These outcomes show that Anxa10 alters depolarization-triggered [Ca²⁺]i dynamics by affecting both Ca²⁺ extrusion from the ER and extracellular Ca²⁺ influx via L-type Ca²⁺ channels. Taken together, these results suggest the possibility that Anxa10 is induced by and involved in [Ca²⁺]i handling, which might reciprocally regulate basal and inducible insulin secretions.

We compared single-cell transcriptomes of Anxa10-positive and Anxa10-negative β-cells from 6- and 10-week-old db/db mice to further identify the genes and pathways associated with Anxa10 expression and altered β-cell function (Supplementary Fig. 5E). A volcano plot examination showed that the genes significantly upregulated in Anxa10-positive β-cells, including FXYD domain-containing ion transport regulator 5 (Fxyd5) and tweety family member 1 (Ttyh1). Fxyd5 encodes the γ-subunit of Na⁺, K⁺-ATPase and elevates Na⁺, K⁺-ATPase pump activity (31), and Ttyh1 encodes a Ca²⁺-binding protein localized to the ER (32). The emergence of these genes corroborated our observations that Anxa10 senses and regulates [Ca²⁺]i in β-cells. Additionally, Fxyd5 reportedly moderates cellular junctions, adhesion, and epithelial-mesenchymal transition (EMT) via its capacity to downregulate E-cadherin (33). Hence, we investigated E-cadherin levels in control and Anxa10-overexpressing MIN6 cells. Immunostaining (Supplementary Fig. 5F) and immunoblotting (Supplementary Fig. 5G) revealed that Anxa10 overexpression lowered E-cadherin protein levels in MIN6 cells. Since EMT enhances the loss of β-cell identity (34), Anxa10 might aid the loss of β-cell identity and type 2 diabetes progression in part by inducing EMT through the induction of Fxyd5 expression.

Pancreatic β-cells are heterogeneous in their morphology and function, and scRNA-seq investigations suggested β-cell subpopulations following DEGs (8,10). Nevertheless, the β-cell heterogeneity during type 2 diabetes progression has not yet been fully clarified. We ordered 2,079 β-cells from 6- and 10-week-old db/db mice in pseudotime according to their transcriptional similarity to better understand population heterogeneity of β-cells focusing on early diabetes progression at the single-cell level. Seurat clustering identified nine transcriptionally distinct β-cell subclusters (Fig. 6A–D and Supplementary Tables 7–10). β-Cell subcluster 0 (Beta-0) revealed upregulated Serpina7 expression, a gene enriched in dedifferentiated β-cells (21) and β-cells with a sustained elevation in [Ca²⁺]i (35), and vitamin D binding protein (Gc), a gene partaking in stress-induced β-cell dysfunction (36). Beta-1 expressed high levels of neurogenin 3 (Neurog3), an important transcription factor of the endocrine cell fate decision (37), indicating that these β-cells lose their mature identity and are dedifferentiated. Beta-2 expressed Fxyd5, Csn3, and Gast and was enriched in IPA pathways linked to oxidative phosphorylation, mitochondrial dysfunction, and Ca²⁺ signaling (Supplementary Fig. 6A). Beta-3 was determined by the differential expression of Dapl1 and Ins1. Beta-4 depicted the upregulation of key genes for β-cell function, such as cytochrome c oxidase subunit 6A2 (Cox6a2), Slc2a2, and insulin receptor–related receptor (Insrr), and was enriched in pathways related to insulin secretion signaling and maturity-onset diabetes of the young. Beta-5 displayed the upregulation of an oxidative stress response gene, sortilin-related VPS10 domain-containing receptor 2 (Sorcs2), and enrichment in metabolic pathways partaking in acetyl-CoA and glutamine biosynthesis. Beta-6 expressed high Ppy and peptide YY (Pyy) levels and was enriched in cancer-related pathways. Beta-7 comprised a δ-cell marker gene, somatostatin (Sst), and pancreatic acinar cell marker genes, e.g., trypsin 4 (Try4) and protease, serine 1 (Prss1) (Supplementary Fig. 6B); this subcluster was enriched in acinar cell–related pathways, including SPINK1 pancreatic cancer pathway (Supplementary Fig. 6C). Beta-8 expressed high levels of ER stress marker genes, e.g., tribbles pseudokinase 3 (Trib3) and DNA damage–inducible transcript 3 (Ddit3). Anxa10 was highly expressed in Beta-0, -2, -3, and -8.

Pseudotime ordering of db/db β-cells recognized four pseudotime trajectories during diabetes progression (Fig. 6E and Supplementary Fig. 6D). The trajectory started with Beta-3, transitioned into Beta-0, and then divided into three major arms: insulin secretion signaling, mitochondrial dysfunction, and dedifferentiation. The insulin secretion signaling arm (Beta-3-0-4) indicates transcriptional β-cell changes that improve insulin secretion and β-cell mass expansion to adapt to enhanced insulin demand. We recognized a mitochondrial dysfunction-ER stress arm (Beta-3-0-2-8). The existence of this arm is in line with studies demonstrating that mitochondrial dysfunction and ER stress aid the functional loss of β-cell mass (38). It also suggests that mitochondrial dysfunction may lead to ER stress. We also recognized two dedifferentiation/transdifferentiation trajectories that branched at Beta-1. One is the PP-like arm (Beta-3-0-1-5-6) and the other is the acinar and δ–like arm (Beta-3-0-1-7). Ppy-lineage β-cells allegedly had a lowered glucose-stimulated Ca²⁺ signaling response and were resistant to extended hyperglycemia (39), revealing that transdifferentiated β-cells in Beta-6 might have factors that are alike to Ppy-lineage β-cells. The dedifferentiation acinar and δ–like branch signifies the existence of β-cells that obtained acinar- and δ-cell phenotypes.

Our scRNA-seq and pseudotime evaluation show that obese diabetic conditions cause β-cell to acinar cell transdifferentiation. Indeed, similar to the islet acinar-cell cluster (cluster 17), Beta-7 was enriched in acinar-related pathways (Fig. 7A), such as the SPINK1 pancreatic cancer pathway (Supplementary Fig. 5C). We deeply evaluated pancreatic sections of 6- and 10-week-old db/db mice to further verify the localization of acinar-like cells in db/db islets. Hematoxylin-eosin staining showed acinar-like cells in islets of 10-week-old db/db mice (Fig. 7B). Immunohistochemical analyses of acinar markers showed that amylase-positive and Ctrb-positive cells were specifically detected in islets of 10-week-old db/db mice (Fig. 7C and D). These outcomes suggest that dedifferentiated β-cells could transdifferentiate into acinar-like cells.

The findings shown so far have signified dynamic modifications in β-cell transcriptional signatures during diabetes onset and progression. We used the CellChat algorithm to display the cell-cell communication in islets mediated by ligand-receptor interactions (19). We combined clusters 1, 2, 7, and 8 depicted in Fig. 1B as a diabetic β-cell cluster and clusters 5 and 14 as a nondiabetic α-cell cluster to compare the number and strength of interactions among different cell populations in 6- and 10-week-old db/db islets. The number and strength of cell-cell communication were higher in 6-week-old prediabetic db/db islets than in 10-week-old diabetic db/db islets (Fig. 8A and B and Supplementary Fig. 7A and B). Specifically, endothelial and stellate cells interacted strongly with prediabetic and diabetic β-cell clusters (Fig. 8B). This outcome was in line with the results of our recent research indicating the strong association between poor islet blood flow and T2D progression (40). The endothelial cluster showed a higher number of angiogenesis-related genes, e.g., FMS-like tyrosine kinase 1 (Flt1) and endothelial cell–specific molecule 1 (Esm1), which are stimulated by vascular endothelial growth factor (VEGF) (Supplementary Fig. 7C). Stellate cluster demonstrated higher fibrosis-related genes decorin (Dcn); collagen type I, α-1 (Col1a1); and Col3a1 (Supplementary Fig. 7D). The products of these genes potentially partake in β-cell survival, proliferation, and function (41).

We evaluated the differential outgoing and incoming signaling patterns of each signaling pathway in 6- and 10-week-old db/db islets to identify ligand-receptor interactions contributing to the signaling changes in db/db islets. We compared the information flow for each signaling pathway and discovered midkine and protease-activated receptor (PAR) signaling as the most highly enriched pathway in prediabetic and diabetic islets, respectively (Supplementary Fig. 7E and F). Furthermore, we detected other important ligand-receptor interactions for both outgoing and incoming signaling, such as insulin, Gcg, endothelin (EDN), fibronectin 1, collagen, somatostatin, VEGF, and chemokine. Remarkably, EDN signaling, which can lower islet blood flow (42), actively participated in the outgoing and incoming signaling patterns between α- and endothelial cells (Fig. 8C and D and Supplementary Fig. 7E and F). This association became stronger with diabetes progression. Furthermore, in line with our findings in Figs. 6 and 7, we identified diabetic islet–specific strong interactions between diabetic β- and acinar cells and between diabetic β-cells through the PAR pathway, including trypsins and serine proteases. We also discovered autocrine PAR signaling in diabetic β-cells.

Pancreatic islets comprising different cell populations categorized according to their transcriptome profiles markedly altered diabetes progression. In the current investigation using a robust scRNA-seq analysis combined with a high-sensitivity cDNA amplification technique, we succeeded in characterizing the transcriptome profiles of each islet cell type, their pathological landscape stream, and cell-cell communication network of islet cells linked to diabetes progression.

We identified diabetes stage–dependent β-cell clusters by comparing prediabetic and diabetic islets. Single-cell transcriptome and pseudotime analysis revealed pathological processes, such as insulin secretion, oxidative phosphorylation, mitochondrial dysfunction, ER stress, and dedifferentiation. We then identified pathways, e.g., oxidative phosphorylation, G-protein–coupled receptor (GPCR)–mediated nutrient sensing in enteroendocrine cells, and mitochondrial dysfunction, which can result in oxidative stress, mitochondrial damage, and Ca²⁺ mobilization, as early pathological processes in prediabetic β-cells. Moreover, the pseudotime analysis indicated a hitherto unknown temporal axis of ER stress through elevated oxidative phosphorylation and mitochondrial dysfunction. Mitochondrial dysfunction and ER stress are closely linked events that consequently stimulate β-cell dysfunction (38). Still, the causal nexus of the two in type 2 diabetes was uncertain. Our work implies that ER stress occurs after mitochondrial dysfunction in the context of β-cell dysfunction.

We paid close attention to the early stage of diabetes development and determined Anxa10 as the topmost DEG of prediabetic β-cells. Anxa10 is a member of the annexin family, the largest superfamily of calcium- and phospholipid-binding proteins in eukaryotic cells, and plays a significant role in various cellular and physiological processes, such as calcium signaling, cell motility, proliferation, and differentiation (28). We showed that Anxa10 was upregulated in β-cells in both type 1 and type 2 diabetes models, and Anxa10 expression was enhanced in response to [Ca²⁺]i elevation in mouse islets and MIN6 cells after various metabolic loads. Furthermore, Anxa10 overexpression experiments implied a function of this Ca²⁺- and phospholipid-binding protein in regulating [Ca²⁺]i allostasis and insulin secretion in β-cells. Pathophysiologically, [Ca²⁺]i overload occurs early in the disease course of diabetes (43), and glucolipotoxic conditions lead to a chronic increase in cytosolic Ca²⁺ influx (44). Anxa10 activation reduced KCl-induced Ca²⁺ flux and KSIS. Following these data, Anxa10 induction is a cellular adaptative response of β-cells to glucolipotoxic conditions by suppressing toxic continuous calcium surges, while Anxa10 upregulation in the prediabetic state might facilitate basal insulin secretion and insulin transcription presumably because of the moderately high basal [Ca²⁺]i. Importantly, Anxa10 is linked to HbA_1c in pancreatic islets from individuals with hyperglycemia (45). This prediabetic gene overexpression in cells demonstrated EMT and antiapoptotic properties. Additional investigations are warranted to define the pathological and organelle-specific roles of Anxa10 in β-cells.

Functional β-cell mass reduction is the primary cause of type 2 diabetes. Our scRNA-seq evaluation shows transcriptional changes in β-cells from prediabetes to diabetes, emphasizing β-cell dedifferentiation and transdifferentiation into acinar-like cells (Figs. 6 and 7). The expression of Aldh1a3, a marker for β-cell dedifferentiation, was detected throughout prediabetic and diabetic β-cells, whereas Neurog3, a marker for endocrine progenitor, was expressed in restricted populations of diabetic β-cells (Fig. 6B and C). Our data imply that the dedifferentiation path is divaricated into two transdifferentiation paths from the Neurog3-enriched β-cell subcluster, one being the acinar and δ-cell path and the other the PP cell path. The transdifferentiation of β-cells to α- (7), δ- (46), and PP (39) cells has been noted in both type 1 and type 2 diabetes. Additionally, our study revealed the transdifferentiation of β-cells into acinar cells in diabetes. β-cell to acinar cell transdifferentiation was also detected in severe acute respiratory syndrome coronavirus 2–infected human islets as the pathomechanism of virus-mediated diabetes (47). Our CellChat analysis also showed the crosstalk between acinar-like and diabetic β-cells through the PAR signaling pathway. The PARs are a family of four GPCRs (PAR1–4) that are initiated by proteolytic cleavage and are essential mediators of hemostasis, inflammation, tissue repair, and other physiological and pathophysiological responses (48). In β-cells, thrombin- or trypsin-mediated activation of PAR3 stimulates insulin secretion by mobilizing Ca²⁺ from intracellular stores (49). Additionally, recent study findings indicated that PAR2 signaling is an important mediator of acinar-β-cell crosstalk that affects β-cell homeostasis and viability (50). Furthermore, our CellChat analysis showed that the PAR signaling pathway was an autocrine factor for diabetic β-cells. Hence, we recommend that the outgoing and incoming interaction of PARs in acinar-like cells and diabetic β-cells might be a beneficial adaptation to enhance insulin secretion from β-cells.

Our study has limitations. The db/db mouse strain has limitations as an animal model for understanding human type 2 diabetes because type 2 diabetes in humans is a polygenic disease, unlike the monogenic basis of leptin receptor–deficient db/db mice. We could not include db/db mice >10 weeks of age for scRNA-seq owing to technical difficulties in obtaining a substantial amount of unbroken islets and careful removal of peri-islet acinar cells, thereby preventing an analysis of the islet cell transcriptome during the late progression of type 2 diabetes at multiple time points. Islet scRNA-seq at different time points of db/db mice could make it eligible to fully comprehend type 2 diabetes occurrence and development. Lineage tracing is required to confirm the transdifferentiation of β-cells to acinar cells. Furthermore, while we identified Anxa10 as a new prediabetic gene, the fate and functional consequence of Anxa10-expressing β-cells and Anxa10-expressing peri-islet acinar cells are unclear. Further analyses of the pathophysiological characteristics of Anxa10-expressing cells in diabetic conditions using genetic lineage tracing in vivo is anticipated to provide additional useful information. Also, we could not detect plasma Anxa10 in mice and humans for now, despite our efforts. Additional investigation of clinical patients with diabetes needs to be conducted for the drug discovery and clinical approach.

In conclusion, we determined diabetes-associated genes, islet cell types, and cell-cell communications in the context of β-cell dysfunction. This landscape provides valuable insights for establishing optimal clinical approaches. Significantly, we have identified a new prediabetic marker, Anxa10, that is upregulated by chronically elevated Ca²⁺ in β-cells, which might initially aid the adaptation for metabolic stress under hyperglycemia but consequently cause diabetic pathology. Furthermore, while most research indicates the correlation of β-cell dedifferentiation/transdifferentiation and β-cell mass, our research proposes that β-cell to acinar cell transdifferentiation is related to β-cell dysfunction in the cellular society of diabetic islets.

Acknowledgments. The authors thank Dr. Osamu Shimomura (University of Tsukuba) for helpful advice and technical discussions. The authors also thank Katsuko Okubo and Chizuko Fuku, Department of Endocrinology and Metabolism, Institute of Medicine, University of Tsukubai, for technical assistance and the members of the Shimano laboratory for discussion and helpful comments on the manuscript.

Funding. This work was supported by Grants-in-Aid for Scientific Research on Innovative Areas program (Inflammation Cellular Sociology) JP17H06395 (to H.Sh.) and JP17H06391 (to K.Ma.) from the Ministry of Education, Culture, Sports, Science, and Technology of Japan, as well as JST COI-NEXT grant number JPMJPF2017.

Duality of Interest. T.Og. is currently an employee of Daiichi Sankyo Co., Ltd. No other potential conflicts of interest relevant to this article were reported.

Author Contributions. K.Mo., T.Ma., S.S., T.N., and S.U. performed most of the experiments. K.Mo., T.Ma., S.S., S.U., K.Ma., and H.Sh. designed the study. K.Mo., T.Ma., S.S., S.U., and H.Sh. wrote the manuscript. K.Mo., S.S., T.Og., H.P., Y.A., T.Ok., K.I., and A.O. performed the data analysis. T.T. and T.Od. contributed material support and to the data interpretation of clinical samples. H.O., S.H., T.Mi., Y.T., M.S., H.So., N.Y., and Y.N. contributed to the experimental design, data analysis, and interpretation of data. All authors contributed to data analysis and the final editing of the manuscript. H.Sh. 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.