Loss of mature β-cell function and identity, or β-cell dedifferentiation, is seen in both type 1 and type 2 diabetes. Two competing models explain β-cell dedifferentiation in diabetes. In the first model, β-cells dedifferentiate in the reverse order of their developmental ontogeny. This model predicts that dedifferentiated β-cells resemble β-cell progenitors. In the second model, β-cell dedifferentiation depends on the type of diabetogenic stress. This model, which we call the “Anna Karenina” model, predicts that in each type of diabetes, β-cells dedifferentiate in their own way, depending on how their mature identity is disrupted by any particular diabetogenic stress. We directly tested the two models using a β-cell–specific lineage-tracing system coupled with RNA sequencing in mice. We constructed a multidimensional map of β-cell transcriptional trajectories during the normal course of β-cell postnatal development and during their dedifferentiation in models of both type 1 diabetes (NOD) and type 2 diabetes (BTBR-Lepob/ob). Using this unbiased approach, we show here that despite some similarities between immature and dedifferentiated β-cells, β-cell dedifferentiation in the two mouse models is not a reversal of developmental ontogeny and is different between different types of diabetes.
Introduction
Insulin-secreting pancreatic β-cells are essential for maintaining blood glucose homeostasis, and their loss or dysfunction underlies all types of diabetes. In type 1 diabetes (T1D), β-cells are targeted by an autoimmune attack. In type 2 diabetes (T2D), β-cells fail due to work overload and a toxic metabolic environment brought about by obesity and peripheral insulin resistance. In recent years, it has become clear that not all β-cells are permanently lost in either type of diabetes. Instead, chronically stressed β-cells lose their functionally mature phenotype and shift to a dysfunctional state in a process called dedifferentiation. Such β-cell dedifferentiation is seen in humans (1–6) as well as in murine models of both T1D and T2D (7,8). The progression to overt diabetes can be prevented if diabetogenic β-cell stress is alleviated in time, before the functionally mature β-cell mass is lost in T2D (9,10). Furthermore, “sleeping” dedifferentiated β-cells that evade immune destruction may be “awakened” to remature and restore some of the lost functional β-cells in recent-onset T1D (3). Thus, drugs that work by directly reversing or preventing β-cell dedifferentiation are critically needed (11,12).
The term “β-cell dedifferentiation” to describe the loss of mature β-cell phenotype was first coined more than two decades ago (13,14). However, what exactly constitutes dedifferentiated β-cells remains debated (15). Previously, it was proposed that β-cells in diabetes dedifferentiate in the reverse order of their normal developmental ontogeny (8). This model predicts that dedifferentiated β-cells resemble β-cell progenitors (Fig. 1, top). An alternative model suggests that β-cell dedifferentiation is a stress type–specific process caused by disruption of specific gene regulatory networks by the diabetogenic environment, thus resulting in a stress type–specific loss of functional maturity, without assuming a “true” cell progenitor identity (16). This model, which we call the Anna Karenina model (based on the opening sentence in Tolstoy’s novel by the same name, “All happy families resemble one another, each unhappy family is unhappy in its own way” [17]), predicts that, in each type of diabetes, β-cells will lose their mature phenotype in a unique manner, depending on how their genetic network is perturbed by a particular diabetogenic environment (Fig. 1, bottom). We reason that distinguishing between the reversal of ontogeny model and the stress type–specific model will help with directing new approaches to restoring mature β-cell function to dedifferentiated β-cells in diabetes.
Here, we test the two models of β-cell dedifferentiation in diabetes. Specifically, we test whether under different types of diabetic stress, dedifferentiated β-cells resemble one progenitor state or whether each type of diabetes produces β-cells that are dedifferentiated in their own way. We do so by elucidating how the transcriptional landscape of β-cells changes during their maturation in normal development, and their dedifferentiation in different types of diabetes, using a β-cell–specific lineage-tracing system in mice. This approach enables us to follow β-cells both during the normal course of their development and during their dedifferentiation in diabetes and allows for direct, unbiased comparison between the gain of β-cell maturation in development and the ways it is lost upon different types of diabetogenic insult.
Research Design and Methods
Mice
All animal experiments were conducted in accordance with the University of Wisconsin-Madison institutional animal care and use committee guidelines under protocol no. M005221. BTBR-Lepob/ob, NOD, and ICR (“WT”) mice were obtained from The Jackson Laboratory and Envigo. Insulin2-Cre;Rosa26-lox-stop-lox-H2BmCherry mice were previously reported (18). Blood glucose and weight were measured in nonfasted animals with the OneTouch Ultra2 glucometer (LifeScan, Milpitas, CA) at the animal facility before islet collection. Mice with blood glucose >300 mg/dL were considered diabetic. Islet isolation was performed as previously described (18,19). Isolated islets were dissociated with 0.25% trypsin-EDTA before sorting through BD FACSAria II for mCherry+ cells.
RNA Sequencing
RNA was isolated from FACS-sorted lineage-traced β-cells from ICR embryos, neonates, and adult mice; NOD adult (diabetic and nondiabetic) mice; and BTBR-Ob/Ob and BTBR-Ob/+ adult mice using phenol chloroform extraction (TRIzol) and QIAGEN RNeasy Plus Mini Kit. The purity and integrity were verified via the NanoDrop One Spectrophotometer and Agilent 2100 Bioanalyzer. DNA libraries were generated with the SMART-Seq v4 Ultra Low Input RNA Kit for Sequencing (Takara Bio, Mountain View, CA) for cDNA synthesis and the Nextera XT DNA Library Preparation Kit (Illumina, San Diego, CA) for cDNA dual indexing. Full-length cDNA fragments were generated from 1–10 ng total RNA by Takara Bio SMART (Switching Mechanism at 5′ end of RNA Template) technology. cDNA fragments were fragmented and dual indexed in a single step via the Nextera kit’s simultaneous transposon and tagmentation step. Quality and quantity of finished libraries were assessed with an Agilent DNA series chip assay (Agilent Technologies, Santa Clara, CA) and a Qubit HS kit (Invitrogen, Carlsbad, CA), respectively. Libraries were standardized to 2 nmol/L. cDNA fragments were sequenced for 1 × 100 base pairs on the Illumina NovaSeq 6000 at sequencing depth 25–30 million reads per sample. Two samples (1 BTBR-Ob/Ob and 1 BTBR-Ob/+) were sequenced for 2 × 150 base pairs on Illumina NovaSeq 6000 with use of standard sequencing by synthesis chemistry at the same sequencing depth. Images were analyzed using the standard Illumina Pipeline, version 1.8.2, and bcl2fastq v2.20.0.422 was used for base calling. Quality control of both single-end and paired-end raw sequencing data was conducted using FastQC-0.11.7 and MultiQC-1.9. All samples passed quality control, as they had uniformly high base quality and sequence quality. Mild adapter contaminations were detected, and we decided not to perform adapter trimming. Raw sequencing data were aligned to mm10 reference genome using Bowtie 1.2.2 under default settings. A gene-by-sample count matrix was estimated using RSEM 1.3.0 under default settings. After combination of the two batches into one data set, genes with average expression >1 were filtered out. Median-by-ratio normalization was conducted on combined data to account for sequencing depth artifact and batch effects. This results in a normalized expression matrix with 63 samples and 16,455 genes.
Hierarchical Clustering and Principal Component Analysis
The normalized expression matrix was log2-transformed and further adjusted for potential batch effects by removeBatchEffect() in the limma R package (v3.44.3). The 15% most variable genes were identified using varFilter() in the genefilter R package (v1.70). Hierarchical clustering was performed on these highly variable genes using Spearman correlational distance and Ward’s linkage method in the cluster R package (v2.1) and visualized on dendextend R packages (v1.14). For principal component analysis (PCA), eigenvectors were calculated using the prcomp() function, and three-dimensional visualization was generated by the Plotly R package (v4.9.2.1).
Differential Expression and Fold Change
Genes having nonzero expression in three or more samples and at least five reads total were retained for differential expression (DE) analysis. DESeq2 (v1.28.1) was used to identify DE genes. Specifically, we applied DESeq2 to obtain gene-specific P values, which were converted to q values using the Benjamini and Hochberg method. A gene was considered DE with q value <0.05 and shrunken log2 fold change >1.5, which was estimated with the adaptive shrinkage estimator from the ashr package inside lfcShrink() in the DESeq2 package. Visualization was done in BioVenn. Gene ontology enrichment analysis was performed with Database for Annotation, Visualization and Integrated Discovery (DAVID), v6.8 (https://david.ncifcrf.gov). Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment was deduced with gauge R package (v2.40.1) with a 0.1 q value cutoff.
Data and Resource Availability
All sequencing data are available in the Gene Expression Omnibus (GEO) repository under accession no. GSE168743. Resources are available from the corresponding author upon reasonable request.
Results
Transcriptional Relationships Between β-Cell Maturation in Postnatal Development and Their Dedifferentiation in Different Types of Diabetes
To test the transcriptional relationship between β-cell maturation and their dedifferentiation in different types of diabetes, we used our previously reported murine β-cell–specific lineage-tracing system (18,20). This system is made by crossing mice transgenic for Insulin2-Cre with mice carrying a floxed reporter of histone H2B fused to mCherry (Rosa26-lox-stop-lox-H2BmCherry). In this system, any cell that had ever expressed the Insulin gene is permanently marked with nuclear mCherry. This reporter mouse line thus enables us to isolate and investigate β-cells through development and functional maturation, as well as through the progression of diabetes, using a single-platform method. We crossed this system into the nonobese diabetic (NOD) model of autoimmune T1D and into the BTBR-Lepob/ob (BTBR-Ob/Ob) model of obesity-related T2D. We FACS purified lineage-traced β-cells from healthy mice during postnatal development, through adulthood, and during the progression to diabetes in the different models. We next subjected the samples to whole-genome RNA sequencing (RNA-seq). We thus generated gene expression data from four time points during β-cell development and maturation (embryonic day [E] 18.5, postnatal day [P]1, P7, and P10), as well as healthy adult mice and diabetic mice (defined as having non-fasting blood glucose levels >300 mg/dL).
We performed unsupervised bottom-up hierarchical clustering of the samples based on the top 15% most variable genes using Spearman correlation as the distance metric. This method identified three large clusters, which we termed “development,” “healthy adult,” and “diabetic” (Fig. 2). Importantly, WT samples (ICR genetic background) and nondiabetic Ob/+ samples (BTBR genetic background) clustered together, without apparent separation between them, confirming that our method correctly distinguishes between the disease conditions and not between genetic backgrounds (see Supplementary Fig. 1 for metadata on all samples, including blood glucose, weight, sex, age, sequencing batch, and genetic background). Interestingly, three of the NOD nondiabetic samples clustered together with the healthy adult samples, and four of the NOD nondiabetic samples clustered with the diabetes samples, suggesting that transcriptional changes related to β-cell stress can be detected before the increase in blood glucose in these mice.
Ontogeny of β-Cell Maturation and Dedifferentiation
To distinguish, in an unbiased manner, between the reversal of ontogeny model and the stress type–specific model of β-cell dedifferentiation in diabetes, we generated a multidimensional trajectory map of the transcriptional states of β-cells as they mature during development and as they lose their mature identity in each of the two types of diabetes (Fig. 3). We reasoned that if the reversal of ontogeny model is correct, then diabetic β-cells are expected to cluster along the developmental trajectory. On the other hand, if the stress type–specific model is correct, then diabetic β-cells will not cluster with any progenitor stage. PCA of the top 15% most variable genes among the groups was used to generate a three-dimensional spatial distribution map of the samples. We found that the first three principal components captured 46.5% of the variation between the samples—PC1 (26.8% of the variation), PC2 (13.5% of the variation), and PC3 (6.2% of the variation)—clearly separating the healthy adult samples, the development samples, and the diabetic samples into three distinct clusters (Fig. 3, left). The P10 samples, which are included in the development group based on their unsupervised hierarchical clustering with this group, are located in the PCA between the earlier developmental samples and the healthy adult samples, confirming the continuity of the normal maturation trajectory. Further separation was seen between the NOD diabetic (T1D) and the BTBR-Ob/Ob diabetic (T2D) samples (Fig. 3, right), although this separation may be amplified by differences in genetic background between the two diabetes models. Again, the NOD nondiabetic samples were divided between the NOD diabetic and the healthy adult samples, indicating that loss of β-cell maturation in NOD mice precedes the onset of overt diabetes. Thus, our analyses using two independent unsupervised mathematical methods suggest that β-cells in the above two diabetes models lose their mature identity but do not return to any developmentally relevant stage.
Gene-Specific Expression Changes in β-Cell Maturation and Dedifferentiation
To validate our unbiased clustering results, we directly examined the expression of a broad list of published markers of mature β-cell identity (19,21–29), “β-cell disallowed” genes (30–32), markers of immature β-cells and non-insulin-expressing β-cell precursors (6,8,9,19,33–36), and islet hormones (Fig. 4). Several markers of immature β-cells and β-cell progenitor genes (MafB, Nnat, Sox17, Fev, Pck1, and Myc), as well as most “β-cell disallowed” genes (Ldha, Hk1, Mylk, Igfbp4, Ndrg2, Pcolce, and Slc16a2), showed downregulation during normal β-cell maturation but were not reexpressed in either type of diabetes. One “disallowed” gene, Ly6a, was reexpressed in the T1D group, and one disallowed gene, Aldh1a3, was reexpressed in the T2D group. Most mature β-cell genes were already present at high levels in the development group (which included semimature P10 pups [37]) and were downregulated in the T2D group and, to a lesser extent, in the T1D group. Of the known β-cell maturation markers, only Ucn3 was upregulated during β-cell maturation and downregulated in both types of diabetes, confirming previous reports by us and others that Ucn3 is one of the most sensitive markers for the fully mature β-cell state (18,19,38,39). Conversely, Dlk1, a marker for immature β-cells (19,40), is downregulated in maturation and is reexpressed in both types of diabetes, while Gast appears to be downregulated in maturation and to be reexpressed specifically in T2D, as was previously reported (33). We did not see reexpression of Neurog3 or any of the other markers of early β-cell precursors in either type of diabetes.
Further comparisons of all genes expressed at >1.5-fold and <1.5-fold with an adjusted P < 0.05 (q < 0.05) in each nonmature condition (development, T1D, and T2D) compared with all samples that clustered together in the healthy adult group showed little overlap among the nonmature groups, confirming our observation that dedifferentiated β-cells in either of the diabetes groups do not revert to a developmentally relevant transcriptional state (Fig. 5). For an additional validation of our findings, we compared the list of differentially expressed genes in our study to those of Neelankal John et al. (41), who reported on RNA-seq analysis of islets from diabetic and nondiabetic LepRdb/db mice. Despite the differences in genetic backgrounds (BTBR vs. C57BL/KsJ), T2D model (Lepob/ob vs. LepRdb/db), and sequencing approach (sorted lineage-traced β-cells vs. whole islets), we found remarkable similarities in gene expression between the two models (Supplementary Fig. 2). A full list of genes in each group, GO terms and biological processes, and KEGG pathways enriched in each of the groups are presented in Supplementary Tables 1–4. Genes that significantly trend between the three NOD groups (NOD nondiabetic samples that clustered with the healthy samples, NOD non-diabetic samples that clustered with the diabetic samples, and the NOD diabetic samples) are shown in Supplementary Fig. 3. These gene-specific analyses confirm distinct gene signatures for β-cell maturation during normal postnatal development and their dedifferentiation in each type of diabetes. We concluded that β-cell dedifferentiation in diabetic NOD and BTBR-Ob/Ob mice is not a reversal of developmental ontogeny and is different for each type of diabetes.
Discussion
Preventing or reversing β-cell dedifferentiation is a promising approach to restoring glycemic control in people with diabetes, in particular T2D. To this end, it is essential to understand the genetic mechanisms leading to β-cell dedifferentiation under different diabetogenic conditions. Several studies over the last decade proposed that β-cells in diabetes dedifferentiate in reverse order of their normal developmental ontogeny. This was shown by the loss of mature β-cell markers in diabetic β-cells, concomitant with the reexpression of several β-cell progenitor genes, such as Neurog3, Sox9, Myc, and in some cases even Nanog and Oct4 (8,9,34,35). Other studies, however, reported the loss of mature β-cell markers in diabetic β-cells without reexpression of progenitor-stage transcription factors (16,18,41) or found that dedifferentiated β-cells resemble immature (neonatal) β-cells to some extent but are not Neurog3-expressing progenitors (12,33). It thus remains debated whether dedifferentiated β-cells in diabetes revert to a progenitor-like state, or whether they lose their mature identity without reverting to any ontogeny-relevant stage and, if so, whether different types of diabetogenic stresses push β-cells to different dedifferentiated trajectories. Distinguishing between these models is important for identifying genetic program-specific intervention points to prevent or reverse β-cell dedifferentiation in diabetes.
We set out to distinguish between the different models of β-cell dedifferentiation in an unbiased manner, using unsupervised analysis of the transcriptional landscapes of both β-cell maturation during development and their dedifferentiation in two mouse models of diabetes, namely, NOD (a model for T1D) and BTBR-Ob/Ob (a model for T2D). We used the same lineage-tracing reporter system to isolate β-cells both during development and during the progression to diabetes in the two different models of the disease. This allowed us to compare β-cell maturation and their dedifferentiation in diabetes using one unperturbed system. We reasoned that superimposing the complete transcriptional states of the diabetic samples on the developmental ontogeny transcriptional map will directly resolve between the two models: if β-cells in diabetes revert to any development-relevant transcriptional state, then the dedifferentiated samples will cluster along the developmental trajectory. On the other hand, if β-cell dedifferentiation in diabetes is not a reversal of ontogeny, despite upregulation of some genes that are expressed also in progenitors, then the dedifferentiated β-cells will not cluster with any development-relevant stage. We report that, despite some similarities between immature and dedifferentiated β-cells, such as reduced expression of several maturation markers and increased expression of some disallowed genes, β-cells dedifferentiation, at least in the two mouse models tested here, is not a reversal of developmental ontogeny and may be different between T1D and T2D.
It is worth noting that our analyses here focused on late (postnatal) β-cell maturation. It is possible that if we compared our diabetic β-cells to earlier (embryonic) progenitors, we would have found that there may be different entry points to β-cell dedifferentiation, but the trajectories eventually converge to a stage resembling embryonic β-cell precursors. However, we did not see reexpression of any known marker of early β-cell precursors, including Neurog3, in either of the diabetes groups, even in samples from mice with severe diabetes for several weeks. Furthermore, our approach using an Insulin promoter–based genetic lineage-tracing system instead of a transient Insulin promoter–driven fluorescent reporter to isolate the cells means that we would have observed dedifferentiated β-cells even if they were to drift into a non-insulin-expressing precursor state. This lineage-tracing system would have also detected β-cell transdifferentiation into other endocrine and nonendocrine cell types should that have been a substantial phenomenon. On the other hand, our use of a constitutive (noninducible) lineage-tracing approach means that our system detects any new cells that turn on insulin expression during or after the progression to diabetes—not only preexisting Insulin-expressing cells. Similarly, our bulk RNA-seq approach may miss unique subpopulations of cells that may respond differently to the diabetogenic stresses and thus would miss any potential mosaicism of β-cell dedifferentiation in the islet under the various conditions. Another aspect of our approach that may confound our results is the possibility of contamination from small numbers of non-β cells, which are hard to sieve out using our bulk RNA-seq approach, such as mature acinar cells (few in early postnatal pancreata and increasing in adults), immune cells (higher probability of contamination in T1D samples), or adipose cells (more abundant in samples from obese diabetic mice). Indeed, several genes associated with such contamination were detected in our comparisons. While such contamination may possibly skew our unsupervised clustering analyses to some extent, our FACS sorting using lineage-traced β-cells and our bulk RNA-seq approaches compensate for such rare events because of the analyses being done on relatively purified β-cell populations and the depth of sequencing, which is not possible with single-cell RNA-seq. Most importantly, our gene-specific analyses using a large list of known markers of β-cell development and maturation provide independent confirmation that dedifferentiated β-cells in diabetic NOD and BTBR-Ob/Ob mice lose their mature β-cell identity but do not return to any developmentally-relevant state. That said, our results do not dispute that β-cells in other models not tested here, such as FoxO1-null mice (8) and mice subjected to a fasting-mimicking diet (42), could return to a Neurog3-expressing progenitor state. With Neurog3 being a master regulator of an embryonic proto-endocrine transcriptional program (43,44), it is conceivable its reexpression in these unique models may force a more developmentally relevant cell identity that is not seen in β-cells from diabetic NOD or BTBR-Ob/Ob mice. It is also important to remember that the NOD and BTBR-Ob/Ob mouse models represent a very narrow view on human T1D and T2D, respectively, and the results obtained from them do not reflect the diversity of stressors and genetic components of the human disease. Finally, while our results show that dedifferentiated β-cells in diabetic NOD and BTBR-Ob/Ob mice are not reverting to any relevant progenitor state, they do not tell us yet what these cells do become. It would be interesting to see whether these dedifferentiated β-cells reflect dysfunctional “empty” or “sleeping” β-cells reported recently in human diabetes (3).
We propose that at least in the case of diabetic NOD and BTBR-Ob/Ob mice, each type of diabetes produces β-cells that are dedifferentiated in their own way, supporting the Anna Karenina model of β-cell dedifferentiation. We hope that these results will provide a valuable resource in the efforts of finding genetic and pharmacological intervention points for preventing and possibly reversing β-cell dedifferentiation in diabetes.
This article contains supplementary material online at https://doi.org/10.2337/figshare.14770101.
Article Information
Acknowledgments. The authors thank N. Sharon (Technion - Israel Institute of Technology), D. Ben-Zvi (The Hebrew University of Jerusalem), A. Helman (The Hebrew University of Jerusalem), and A. Attie (University of Wisconsin-Madison) for valuable comments and discussion. The authors are grateful to all present and former members of the Blum laboratory, especially Melissa Adams, Jennifer Gilbert, Bayley Waters, Emily Cade, Ron Fleminger, and Emily Maritato, for help on this project. The authors are also grateful to the University of Wisconsin-Madison Biotechnology Center Gene Expression Core for RNA-seq.
Funding. This work was funded by National Institute of Diabetes and Digestive and Kidney Diseases grant 1R56-DK-115837, JDRF grant 2-SRA-2018-621-S-B (to B.B.), and NIH Small Instrument Grant 1S10RR025483-01 to the University of Wisconsin-Madison School of Medicine and Public Health Flow Cytometry Core. M.N.B. acknowledges the support of a postdoctoral fellowship from the Morgridge Institute for Research.
Duality of Interest. No potential conflicts of interest relevant to this article were reported.
Author Contributions. S.D.N. and B.B. contributed to conceptualization. B.B., C.K., and S.D.N. contributed to methodology. S.D.N. contributed to data acquisition. S.D.N, M.N.B., Z.N., and J.B. contributed to formal analysis. S.D.N. and B.B. contributed to writing the original draft. All authors contributed to writing, review, and editing. B.B. contributed to funding acquisition. C.K. and B.B. contributed to supervision. B.B. 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. A non–peer-reviewed version of this article was submitted to the bioRxiv preprint server (https://www.biorxiv.org/content/10.1101/2021.02.16.431507v1) on 17 February 2021.