Targeted pharmacological therapy restores β-cell function for diabetes remission

Dedifferentiation of insulin-secreting β cells in the islets of Langerhans has been proposed to be a major mechanism of β-cell dysfunction. Whether dedifferentiated β cells can be targeted by pharmacological intervention for diabetes remission, and ways in which this could be accomplished, are unknown as yet. Here we report the use of streptozotocin-induced diabetes to study β-cell dedifferentiation in mice. Single-cell RNA sequencing (scRNA-seq) of islets identified markers and pathways associated with β-cell dedifferentiation and dysfunction. Single and combinatorial pharmacology further show that insulin treatment triggers insulin receptor pathway activation in β cells and restores maturation and function for diabetes remission. Additional β-cell selective delivery of oestrogen by Glucagon-like peptide-1 (GLP-1–oestrogen conjugate) decreases daily insulin requirements by 60%, triggers oestrogen-specific activation of the endoplasmic-reticulum-associated protein degradation system, and further increases β-cell survival and regeneration. GLP-1–oestrogen also protects human β cells against cytokine-induced dysfunction. This study not only describes mechanisms of β-cell dedifferentiation and regeneration, but also reveals pharmacological entry points to target dedifferentiated β cells for diabetes remission. β-cell dedifferentiation has emerged as a contributing mechanism in type 1 and type 2 diabetes development. Here Sachs et al. show that a pharmacological treatment that combines insulin and a GLP-1–oestrogen conjugate reverses dedifferentiation, and improves β-cell function and hyperglycaemia in diabetic mice.

T he progressive loss or dysfunction of insulin-producing β-cell mass ultimately leads to type 1 diabetes (T1D) or type 2 diabetes (T2D), respectively 1 . Current pharmacological treatments do not stop the decline of β-cell function and number that leads to glucose excursions and, eventually, to devastating microvascular and macrovascular complications. Hence, the ideal treatment should be initiated when the first diabetic symptoms appear, to protect or regenerate glucose-sensing and insulin-secreting β cells for optimal blood glucose regulation and to prevent secondary complications. Recently, T1D progression has been halted by anti-CD3 immunotherapy for 2 years (ref. 2 ), but it will be important to test whether additional β-cell regenerative therapy can further or permanently delay the onset of diabetes. Intensive insulin therapy at disease onset has been shown to partially restore β-cell function, which slows the disease progression in patients with T1D and T2D [3][4][5] . Despite both therapies providing similar glycaemic control, early intensive insulin treatment, compared with oral anti-diabetic drugs, in patients with T2D better preserves β-cell function, which suggests additional glucose-independent beneficial effects of insulin therapy 6,7 . Thus, an understanding of the mechanisms of β-cell dysfunction and pharmacological replenishment is urgently required to stop or reverse diabetes progression and to improve the therapeutic options for patients. Dedifferentiation of β cells has been observed in genetic mouse models of T1D and T2D as well as in patients with diabetes, and is characterised by the loss of the expression of key maturation marker genes (for example, Slc2a2 (also known as Glut2) and Ucn3) and by impaired insulin secretion, and thereby contributes to β-cell dysfunction and hyperglycaemia [8][9][10] . To investigate whether dysfunctional β cells under hyperglycaemic conditions can be targeted pharmacologically to restore β-cell function, we explored the multiple-low-dose model of streptozotocininduced diabetes (mSTZ) in mice. STZ specifically ablates β cells, but when STZ is injected in multiple low doses, some residual β cells can survive 11 . Furthermore, the absence of genetic lesions and GLP-1-oestrogen and insulin polypharmacotherapy. Oestrogen and glucagon-like peptide 1 (GLP-1) have been repeatedly implicated in the treatment of diabetes due to insulinotropic and β-cell protective effects in preclinical studies 16,17 . Chemically optimised GLP-1 analogues profoundly improve glucose and body weight management in obese people and individuals with T2D 18 . However, severe gynaecological, oncogenic and mitogenic side effects preclude chronic oestrogen use as a drug for diabetes and, as yet, GLP-1 analogues have failed to preserve β-cell function and mass in obese and diabetic humans 17 .
To circumvent the gynaecological, oncogenic and mitogenic actions of oestrogen, we recently designed and evaluated a stable GLP-1-oestrogen conjugate, which reversed the metabolic syndrome in diet-induced obese male and female mice 19 . Here, we used the GLP-1-oestrogen conjugate to test whether the specific delivery of oestrogen into the GLP-1 receptor protein (GLP-1R)-expressing β cells could restore β-cell functionality. GLP-1-oestrogen treatment for 100 d was more efficacious in decreasing fasting glucose (Fig. 1b) and increasing fasting C-peptide (Fig. 1c) and insulin levels (Extended Data Fig. 2a) than either of the monoagonists (oestrogen or GLP-1 alone). Moreover, only GLP-1-oestrogen treatment improved pancreatic islet architecture (Fig. 1d) and increased β-cell number, as compared with those in mSTZ-diabetic mice treated with vehicle (Fig. 1e,f). These effects were independent of body weight loss (Extended Data Fig. 2b).
Polypharmacotherapy holds the potential to simultaneously activate redundant or additive pathways to enhance efficacy, and to enable reduced dosing of the individual components and consequently reduce the risk of unwanted side effects 20 . We tested the combination of insulin and GLP-1-oestrogen to investigate a triple pharmacological approach to enhance the efficacy of both compounds, and particularly to lessen the amount of insulin required. The combination therapy normalised glycaemia (Fig. 1b) and increased C-peptide levels ( Fig. 1c), as compared with those in mSTZ-diabetic mice treated with vehicle. Furthermore, the combination therapy displayed a superior effect compared to treatment with insulin alone to limit weight gain (Extended Data Fig. 2b), normalise glucose tolerance (Extended Data Fig. 2c), increase pancreatic insulin content (Extended Data Fig. 2d) and increase β-cell number (Fig. 1e,f). Importantly, we were able to reduce the insulin dose by 60% (10 nmol kg −1 ) compared with that in the insulin monotherapy (25 nmol kg −1 ) and still achieve superior therapeutic outcomes, which reduces the risk of hypoglycaemia and unintended weight gain.

GLP-1-oestrogen targets β cells.
We next wanted to confirm the absence of systemic toxicity that was related to the oestrogen component of the GLP-1-oestrogen conjugate, a pre-requisite for clinical use. To that end, we investigated whether GLP-1-oestrogen (doses up to 10× higher than are generally used in mouse experiments, and at least 1,000× the plasma oestradiol exposure as compared with that in women on hormone replacement therapy (see Methods)) increased uterus weight in ovariectomised (OVX) rats after two weeks of treatment (Extended Data Fig. 3a and Methods).  In contrast to treatment with oestrogen alone, no treatment-related effect was observed with the GLP-1-oestrogen conjugate (Extended Data Fig. 3b), which is consistent with previously reported results in OVX mice 19 .
To confirm β-cell-specific targeting of the GLP-1-oestrogen conjugate, we used a double knock-in fluorescent reporter mouse model (Foxa2-Venus Fusion (FVF) × Pdx1-BFP (blue fluorescent protein) fusion (PBF); FVFPBF DHom ) 21 , which allows αand β-cell sorting (Extended Data Fig. 3c,d). Male FVFPBF DHom mice develop maturity onset diabetes of the young owing to reduced Pdx1 levels in islets, accompanied by hyperglycaemia, reduced β-cell number and impaired islet architecture at weaning age 21 . In this genetic diabetes model, none of the therapies used in this study improved glycaemia after four weeks of treatment (Extended Data Fig. 3e,f). These results suggest that β-cell function cannot be restored pharmacologically in the presence of the genetically induced β-cell lesions in the FVFPBF DHom mice. However, the GLP-1-oestrogen conjugate, but not the monoagonists, specifically increased β-cell granularity of FVFPBF DHom mice in comparison to control mice, which shows that the conjugate selectively targets β cells (Extended Data Fig. 3g).

GLP-1-oestrogen improves human β-cell function.
To provide human relevance to the findings, we tested GLP-1-oestrogen and the monoagonists in human micro-islets in the absence or presence of β-cell stressors (cytokine cocktail, see Methods). After acute compound exposure, GLP-1-oestrogen was more potent than either of the individual components in increasing glucose-stimulated insulin secretion (GSIS) from human micro-islets in comparison to that in untreated and monoagonist-treated micro-islets (Fig. 3a). We next exposed the human micro-islets to cytokines to determine whether the beneficial effects of GLP-1-oestrogen protect against stress-induced impairment of β-cell functionality (Fig. 3b). Seven-day treatment with GLP-1-oestrogen enhanced GSIS, and exceeded the effects that were seen with either of the individual components (Fig. 3b). Moreover, only GLP-1-oestrogen treatment increased the total insulin content of cytokine-exposed human micro-islets (Fig. 3c). This was independent of changes in the total ATP content (Extended Data Fig. 4a) and in caspase luciferase activity (Extended Data Fig. 4b), which suggests that the compound treatment improved functionality, but did not improve cell survival of human micro-islets. These results show that GLP-1-oestrogen is superior to both of the monoagonists in improving β-cell function in homoeostasis, and in acting upon cytokine stress in both mice and humans. β-cell heterogeneity in homoeostatic and healthy mice. To elucidate the molecular mechanisms that underlie β-cell failure in mSTZ diabetes, and β-cell recovery after the different therapeutic approaches, we performed scRNA-seq of isolated islets from mice that responded to treatment (Extended Data Fig. 5).
In normal islet homoeostasis, we identified the four main endocrine cell subtypes, α cells, β cells, δ cells and pancreatic polypeptide (PP) cells, by unbiased graph-based clustering. Clusters were annotated on the basis of predominant endocrine hormone expression of glucagon (Gcg), insulin (Ins), somatostatin (Sst) and PP (Extended Data Fig. 6a,b and Methods). For each of the four endocrine subtypes we identified a specific marker gene signature (Supplementary  Table 1 and Methods). Refined clustering of insulin-positive cells revealed the presence of two main β-cell subpopulations, β1 and β2 (Extended Data Fig. 6c and Methods). Single-cell trajectory inference suggests a continuum of transcriptional states, rather than discrete phenotypes within β cells and a transition between β1 and β2 subpopulations (Extended Data Fig. 6d and Methods). We found a progressive increase in expression of β-cell maturation marker genes (for example Ins1, Ins2 and Ucn3 (ref. 22 )), and genes of the secretion machinery (G6pc2, Sytl4 and Slc2a2), as well as a concomitant decrease in expression of the β-cell immaturity gene Mafb 23 , and pan-endocrine lineage marker genes (Chga and Chgb) along the pseudotime trajectory from β2to β1-cells (Extended Data Fig. 6d and Supplementary Table 2). The expression of transcription factor genes that are associated with β-cell identity (for example Pdx1, Nkx6.1 and NeuroD1) was unchanged (Extended Data Fig. 6d). The β2-cell cluster was characterised by the downregulation of genes that are involved in insulin secretion, oxidative phosphorylation and cell-cycle inhibition, as well as by the upregulation of genes that are involved in cAMP and WNT signalling (Extended Data Fig. 6e), both of which observations are suggestive of a more immature and/ or proliferative state for β2 cells 24 . We consistently observed upregulation of cell-cycle-associated genes, such as Ki67 and Cdk1 in the immature β2 subpopulation and, in accordance with this result, 16 out of 403 of the β2-cells, but only 2 out of 5,319 of the mature β1-cells, were classified as being involved in cell cycling (Extended Data Fig. 6f (see Methods)). Taken together, these results confirm the co-existence of mature (β1) and immature and/or proliferative β cells (β2) in healthy mouse islets 24,25 . In addition, we found subpopulations of polyhormonal cells that could be distinguished from doublets and ambient RNA, both of which are common problems with the scRNA-seq technology (Extended Data Fig. 6g,h and Methods).
β-cell dedifferentiation at the single-cell level. To obtain a deeper understanding of the cell autonomous and non-cell autonomous effects that underlie chemical β-cell ablation in the islet cell niche, we performed a scRNA-seq survey of diabetic islets after 100 d of persistent hyperglycaemia. Unsupervised clustering and embedding of the scRNA-seq data revealed altered endocrine subtype composition and cell-intrinsic gene expression profiles that were indicated by a shifted cell cluster location in the Uniform Manifold Approximation and Projection (UMAP) space of mSTZ-diabetic compared to healthy mice ( Fig. 4a and Supplementary Table 3). In particular, there was a threefold decrease in the proportion of β cells from mSTZ-treated mice (β-mSTZ), and β-mSTZ formed a cluster that was clearly distinct from healthy β cells (Fig. 4a).
In contrast to the results for β cells, few transcriptional changes were detected for α cells, δ cells and PP cells (Extended Data Fig. 7a-f and Supplementary Table 3). We next sought to describe the progression from healthy to dysfunctional β cells and the associated gene expression changes, using single-cell trajectory inference. Cells were ordered on the basis of a cell-to-cell distance metric that was calculated using the concept of diffusion pseudotime (see Methods). We identified a cellular trajectory in which cells transitioned from mature to immature to β-mSTZ cells (Fig. 4b).
Remaining β-mSTZ cells expressed low Ins1 and/or Ins2 messenger RNA and also showed sustained low expression of β-cell identity transcription factor genes, such as Pdx1, Nkx2.2, Nkx6.1, Pax6, Isl1 and NeuroD1 (refs. [26][27][28][29][30][31] (Fig. 4b). Pathways that are associated with β-cell maturity and functionality included genes with downregulated expression, whereas genes in ER stress and oxidative phosphorylation pathways were upregulated in the remaining β-mSTZ cells; together these observations indicate an ER stress response and β-cell dysfunction ( Fig. 4c and Supplementary Table 3). Along this trajectory, expression of key markers of β-cell maturity and functionality gradually decreased concomitantly with an increase in the very few known markers of β-cell immaturity and dedifferentiation (for example Aldh1a3 (ref. 32 ) and gastrin 33 ) (Fig. 4b,d). Strikingly, our single-cell analysis uncovered a large number of upregulated genes and pathways in β-mSTZ cells that are not expressed at all or are expressed only subtly in mature, functional murine β cells (Fig. 4e). We confirmed the increased expression of, for example, Cck and Slc5a10 proteins by immunohistochemistry in mSTZ diabetic mice (Fig. 4f). These identified targets may serve as biomarkers for dysfunctional β cells, and may have the potential to be part of druggable pathways to restore β-cell function. Some of these were also identified recently in β cells and pancreata of T1D and T2D human specimens 34,35 (Extended Data Fig. 7g,h). There is some debate as to whether β-cell dedifferentiation resembles reversal to a pluripotent (Oct3/4, Nanog, Sox2) or endocrine progenitor state (Neurogenin 3 or Neurog 3, hereafter called Ngn3), whether it is part of normal phenotypic variation described as β-cell heterogeneity, or whether it resembles a glucotoxic-induced reversible state 36 . Although dedifferentiated β cells had been characterised previously by upregulation of pluripotency or endocrine transcription factors 8,9,15 , in our study expression levels of Sox9, Pou5f1 (Oct3/4), Myc and Ngn3 genes were unaltered in mSTZ-treated β cells (Extended Data Fig. 7i).
To further characterise the maturation state of mSTZ-derived β cells, we compared our data set to β-cell expression profiles during embryonic (E17.5) to postnatal development (P60) 37 . We assessed transcriptional similarity of β-cell subpopulations using Partitionbased graph abstraction (PAGA) after data integration and computation of a common embedding. The enforcement of integration using only genes associated to β-cell maturation in the reference data allowed us to match maturation states independent of differences in other biological processes between the two data sets (see Methods). PAGA uses a statistical model to measure the relatedness of groups of single cells (see Methods). We found that dedifferentiated β-mSTZ cells were more strongly connected to early time points of the maturation data set, whereas β2 cells clustered to intermediate time points and β1 cells to late time points (Fig. 5a). Remarkably, the inferred cellular trajectory from β-mSTZ to β2 to β1 cells aligned with the trajectory of β-cell maturation from embryonic (E17.5) to mature β cells (P60) of the reference data set (Fig. 5b). An increase in known β-cell maturation and a decrease of immaturity markers along this trajectory further indicated that the cells follow a similar differentiation or dedifferentiation programme (Fig. 5c,d). To validate these findings, we scored each β cell using the gene sets that are characteristic for the start (E17.5/P0) and end Secretion (mean ± s.e.m.) of donor 1 = 0.13 ± 0.02 ng ml −1 , donor 2 = 0.56 ± 0.04 ng ml −1 and donor 3 = 0.51 ± 0.06 ng ml −1 after chronic vehicle (no stress) exposure. Box plot of all data points. Line indicates the median. # P indicates P value to vehicle under the no stress condition. *P indicates P value to vehicle cytokine exposure. Data from no stress versus cytokine stress were analysed by an unpaired two-sided t-test (t = 4.776, d.f. = 33). Otherwise, data were analysed by a one-way ANOVA with donor as random effect followed by Tukey post-hoc (F low dose (4, 82) = 6.35; F medium dose (4, 81) = 7.65; F high dose (4, 81) = 21.91). c, Total insulin content of human micro-islets after 7-d exposure to cytokine stress and effect of oestrogen, GLP-1 or GLP-1-oestrogen treatment. n = 6 micro-islets of n = 3 human donors for each condition. Insulin content (mean ± s.e.m.) of donor 1 = 41.11 ± 3.73 ng per islet, donor 2 = 30.86 ± 3.36 ng per islet and donor 3 = 82.73 ± 3.99 ng per islet after chronic vehicle (no stress) exposure. Box plot of all data points. Line indicates the median. # P indicates P value to vehicle no stress condition. *P indicates P value to vehicle cytokine exposure. Data from no stress versus cytokine stress were analysed by an unpaired two-sided t-test (t = 7.429, d.f. = 32). Otherwise, data were analysed by a one-way ANOVA with donor as random effect followed by Tukey post-hoc (F low dose (4, 80) = 4.12; F medium dose (4, 78) = 3.01; F high dose (4, 79) = 3.31). Additional information is available in Source data.
(P60) point of the developmental trajectory (Supplementary Table  4 and Methods). Healthy mature β cells (β1) scored high for maturity genes, whereas healthy immature cells (β2) and dedifferentiated cells (β-mSTZ) scored higher for the embryonic immaturity gene set (Fig. 5e). Taken together, these results imply that during the transition from healthy β1 to β2 to dedifferentiated β-mSTZ, β cells revert, at least in part, back to a more immature and, further, to an embryonic-like state. To separate the altered maturation state from other processes induced in dedifferentiated β cells, we compared differentially regulated gene ontologies and pathways between embryonic (E17.5/P0) and mature (P60) β cells and between β-mSTZ cells and healthy control β cells (see Methods). According to the trajectory and PAGA analysis, these β-cell states of the reference data set correspond best to dedifferentiated β-mSTZ cells and mature, healthy β cells, respectively. Embryonic and mSTZ-diabetic β cells shared downregulation of molecular processes connected to β-cell function and maturity (for example, insulin secretion and FoxO signalling) compared to healthy mature β cells (Extended Data Fig. 8a and Supplementary  Table 4), whereas genes involved in oxidative phosphorylation and gene and protein transcription were upregulated (Extended Data Fig. 8b and Supplementary Table 4). Specific to embryonic β cells was a downregulation of lipid and carbohydrate metabolism and an upregulation of WNT signalling (Extended Data Fig. 8a,b). This corresponds to known mechanisms of β-cell maturation during embryogenesis 37 . Interestingly, in mSTZ-diabetic β cells, but not the embryonic β cells, we found an upregulation of pathways and ontologies associated with ER stress and a decreased expression of genes involved in insulin and MAPK signalling in comparison to healthy mature β cells (Extended Data Fig. 8a,b).
Thus, β-cell dedifferentiation involves partial reversal to an embryonic or immature β-cell programme and upregulation of an ER stress response and altered signalling state. These results suggest that surviving β cells are dedifferentiated, which shows that mSTZinduced diabetes is a good model in which to study mechanisms of β-cell dedifferentiation and redifferentiation in the absence of genetic lesions.

Mechanisms of β-cell redifferentiation.
In line with the pharmacological data, single-cell analysis of the different treatments revealed that β cells of mice treated with vehicle (Extended Data Fig. 9a), oestrogen (Extended Data Fig. 9b) and GLP-1 (Extended Data Fig. 9c) remained dedifferentiated. By contrast, we observed an increased fraction of immature β2-cells from GLP-1-oestrogentreated mice (Extended Data Fig. 9d). In PEG-insulin-treated mice (Extended Data Fig. 9e) and GLP-1-oestrogen plus PEG-insulincotreated mice (Extended Data Fig. 9f), almost no dedifferentiated β cells remained and most cells clustered with immature β2-cells.
To further assess the transcriptional state of β cells from the treated mice, we calculated a cell-to-cell distance so that cells could be ordered along the cellular trajectory from dedifferentiated to healthy β cells (see Methods). On this one-dimensional axis, β cells of mice treated with PEG-insulin or with the combination of PEG-insulin and GLP-1-oestrogen were located closest to β cells from healthy mice (Fig. 6a-c). This transcriptional similarity to healthy β cells was further supported by PAGA (see Methods). In the PAGA graph,   Table 3). Cells were pooled from mice treated with no STZ-vehicle (n = 3) and mSTZ-vehicle (n = 3). We used limma-trend to find differentially expressed genes (see Methods). Gene enrichment was done with EnrichR using Fisher's exact test to identify regulated ontologies or pathways (see Methods). d, Immunohistochemical analysis of Glut2, Ucn3, Aldh1a3 and gastrin in β cells of mSTZ mice and healthy mice at study end. Images are representative of mice treated with no STZ-vehicle (n = 3) and mSTZ-vehicle (n = 3). Scale bars, 50 μm; scale bar zoom-in, 20 µm. e, Gene expression along the trajectory from β1 to β-mSTZ (as in c) of 29 genes specifically expressed in β cells from mSTZ-treated mice (expression in <5% of no-STZ-β cells and in >25% of β-mSTZ cells, see Methods). Cellular locations of associated proteins are indicated. f, Immunohistochemical analysis of Slc5a10 and Cck in β cells of mSTZ and healthy mice at study end. Images are representative of mice treated with no STZ-vehicle (n = 3) and mSTZ-vehicle (n = 3). Scale bars, 50 μm; scale bar zoom-in, 20 µm.
β cells of PEG-insulin-treated and GLP-1-oestrogen plus PEGinsulin-cotreated mice showed the strongest connection to healthy β cells (Fig. 6d). The observed overall re-establishment of the healthy β-cell expression profiles was substantiated by an increased expression of β-cell maturity markers and decreased expression of immaturity and dedifferentiation markers along the inferred trajectory ( Fig. 6e-g). Moreover, Ucn3 expression recovered during the pharmacological treatment (Extended Data Fig. 10a). This shows

Ucn3
Maturity score Embyronic or immaturity score that the maturation state before treatment was different from that achieved after treatment. Hence, upon PEG-insulin or PEG-insulin plus GLP-1-oestrogen treatment, β cells adopt a molecular immature yet functional phenotype that is sufficient for blood glucose normalisation and diabetes remission. β cells from mice treated with PEG-insulin or the combination of PEG-insulin and GLP-1-oestrogen were grouped in distinct β-cell subpopulations, albeit at a similar maturation state (Fig. 6b). This implies the existence of a compound-specific mechanism-of-action (MOA) that underlies the recovery of β-cell function. To investigate the distinct MOAs of the different treatments, we identified the β-cell-specific transcriptional signature of treated and mSTZderived β cells (Supplementary Table 5). An increased expression of genes involved in β-cell functionality and maturation (Fig. 6h) was common to both treatments and may result from improved blood glucose levels and/or stimulation of shared pathways of insulin and GLP-1-oestrogen signalling. GLP-1, oestrogen and insulin receptor activation regulate MAPK and FoxO signalling [38][39][40] , both of which were increased after PEG-insulin, GLP-1-oestrogen and PEGinsulin and GLP-1-oestrogen co-therapy (Fig. 6h). Although we cannot dissect the signalling contribution of each individual receptor, we believe that the polypharmacological approach might potentiate the simultaneous activation of commonly regulated pathways.
Unexpectedly, treatment with PEG-insulin elicited a β-cellspecific stimulation of the insulin signalling cascade as well as stimulation of the recently characterised RNA polymerase II mediated pathway 41 (Fig. 6h). Hence, our data suggest that direct effects of insulin on β cells contribute to the improvement of β-cell function and recovery, as was proposed for T2D 42,43 .
Our goal was to use GLP-1-oestrogen to selectively deliver oestrogen into β cells. We consistently found β-cell-specific induction of the ER-associated degradation (ERAD) pathway and induction of transfer RNA signalling in PEG-insulin and GLP-1-oestrogen cotreated mice (Fig. 6h). ERAD mitigates ER stress, which, when unresolved, contributes to functional β-cell mass loss in T1D and T2D 44 . Chemical and genetic disturbances of ERAD impair β-cell function 44 . We found an increased proinsulin-to-C-peptide ratio in mSTZ-diabetic mice, which is used as an ER-stress surrogate in diabetes 45 (Extended Data Fig. 10b). GLP-1-oestrogen, PEG-insulin and GLP-1-oestrogen and PEG-insulin co-therapy normalised this ratio (Extended Data Fig. 10b). Recently, it has been reported that oestrogen, via nuclear oestrogen receptor alpha signalling, stabilizes the ERAD proteins Sel1l and Hrd1in β cells, an effect associated with diabetes amelioration in Akita mice 46 . We observed increased costaining for insulin and Sel1l in GLP-1-oestrogen and PEG-insulin cotreated islets only 25 d after treatment initiation (Extended Data Fig. 10c). In addition to upregulation of Sel1l and Hrd1 (also known as Syvn1), ERAD-associated gene expression (for example, expression of Sdf2l1, Herpud1, Dnajb11, Dnajb9, Derl3 and Hspa5 genes) was specifically increased in β cells of GLP-1-oestrogen and PEG-insulin cotreated mice (Extended Data Fig. 10d). The Sdf2l1, Herpud1 and Hspa5 genes encode ERAD-proteins that have beneficial effects on β cells that allow correct insulin folding and/or function [47][48][49] . Dnajb9 and Dnajb11 are chaperone proteins that might aid correct insulin protein folding. Derl3 is required for ER-associated degradation 50 . The specific role of Derl3 in β cells is unknown, but interestingly, Derl3 expression was shown to protect cardiomyocytes against ER-stress-induced death by enhancing ERAD 50 . These results support a role for ERAD activation by GLP-1-oestrogen and PEG-insulin cotreatment that induces a treatment-specific molecular profile for the protection and regeneration of β cells (Fig. 6b).
In a similar manner, tRNA signalling is known to be an intracellular target of oestrogen, and an increased abundance of tRNA has been associated with proliferating cells 51,52 . Indeed, from the single-cell data, we observed the highest fraction of proliferative β cells in the PEG-insulin and GLP-1-oestrogen conjugate cotreated mice (Fig. 7a). Moreover, we found increased β-cell proliferation in GLP-1-oestrogen and PEG-insulin cotreated mice that was already evident after 25 d of treatment, and that, importantly, was not evident in the single-treatment groups (Fig. 7b). Stressed β cells, such as those under chronic hyperglycaemic conditions, lack an adequate response to GLP-1 therapy probably as a result of decreased expression of GLP-1R [53][54][55] . Here, we reasoned that the restoration of glycaemia in mSTZ mice, notably through additional chronic insulin therapy, may increase the expression of GLP-1R. Indeed, we found progressively increased levels of GLP-1R on the surfaces of β cells of mice with improved glucose levels (Fig. 7c). This may have facilitated the enhanced delivery, uptake and action of oestrogen and GLP-1 in β cells, especially for those cells cotreated with GLP-1oestrogen and insulin.
To examine whether other endocrine cells have contributed to the regeneration of functional β cells, we explored cluster relations and possible cellular transitions using PAGA and RNA velocity estimation (see Methods). For this we included β-mSTZ cells as the origin (starting point) of treated cells and we investigated where cells moved from that baseline (Fig. 8a). We found no direct connection or cell movement from other (non-β) cell populations towards redifferentiated β cells. We next examined the RNA velocity and potential fate of treated endocrine cells (Fig. 8b). The velocity of some of the immature β cells observed after GLP-1-estrogen, PEG-insulin and combined treatment pointed towards mature β cells of healthy mice, thus further substantiating β-cell redifferentiation. Moreover, the scRNA-seq data suggested that neogenesis was not increased after 100 d of treatment, as the expression levels of Ngn3 mRNA remained unchanged in endocrine cell subtypes (Fig. 8c). We also found no indication, from Ngn3 immunostaining  Table 5). We used limma-trend to find differentially expressed genes (see Methods). Gene enrichment was done with EnrichR using Fisher's exact test to identify regulated ontologies or pathways (see Methods). Cells of n = 3 mice per treatment were pooled. in tissue sections of earlier time points, that overt neogenesis contributed to β-cell regeneration (Fig. 8d). Together, these results suggest that redifferentiation of β cells along dedifferentiation and redifferentiation trajectories is the main mechanism that underlies the re-establishment of functional β cells in the mSTZ model by treatment with GLP-1-oestrogen or PEG-insulin, or by their cotreatment. By using a combination of low dose insulin with GLP-1-oestrogen treatment, we were able to trigger a β-cell-specific transcriptional response that was characterised by increased β-cell proliferation and enhanced functionality.

Discussion
Herein we have established the mSTZ model of diabetes as a model to study β-cell dysfunction and dedifferentiation. Single-cell profiling of remaining β cells identified many markers of β-cell dedifferentiation that were undescribed previously and that code for surface molecules, receptors and secreted proteins. These may be used as biomarkers or may allow the detection, isolation and characterisation of dedifferentiated β cells. This could reveal pathomechanisms of T1D and T2D and may have the potential to identify unique diagnostic markers and therapeutic targets. Using scRNA-seq we were able to delineate a β-cell fate trajectory in which cells transitioned from mature to immature to dedifferentiated β cells, which implies that β cells can be characterised by a continuum of transcriptional states that reflect discrete phenotypes. Inference of cell transitions using the RNA velocity concept further suggested that there  was no ongoing transdifferentiation from other non-β and nonendocrine cells towards dedifferentiated β cells. Upregulation of the endocrine master regulator Ngn3 might depend on the severity of hyperglycaemia; high glucose levels ( >33 mM) have been shown to induce Ngn3 expression 9,15 , whereas lower levels ( <25 mM) had no effect 14,56 . We show that β-cell dedifferentiation in mSTZ-diabetic mice is independent of induction of Ngn3 + endocrine progenitors and that the transcriptional state of dedifferentiated β cells is more similar to late embryonic or early postnatal β cells.
Recently, the Kushner laboratory has provided evidence that some of the remaining insulin in the blood stream of patients with long-term T1D 57 originates from dedifferentiated β cells and/or from polyhormonal non-β cells that function as 'insulin microsecretors' 58 . Similarly, the Korsgren laboratory found histological evidence for β-cell dedifferentiation at T1D onset 59 . Therefore, triggering the redifferentiation of dedifferentiated β cells seems to be an intuitive approach for the treatment of diabetes that does not involve β-cell proliferation or neogenesis per se 60 . Preclinical as well as clinical findings from patients with type 1 and 2 diabetes suggest that a transient recovery of β-cell dysfunction occurs upon glycaemia normalisation by intensive insulin treatment by either β-cell rest or redifferentiation 3,15 . By using scRNAseq we can dissect endocrine subtype-specific treatment responses and show that insulin treatment triggers transcriptional changes in β cells that are connected to insulin and/or IRS signalling. This supports the idea that in addition to lowering the glucotoxic stress on β cells, direct insulin-or IGF-signalling improves β-cell health and performance and can redifferentiate β-cell mass in diabetic models 43 . Importantly, the redifferentiated β cells that were induced by insulin therapy were functional and responded to physiological stimuli, as indicated by increased fasting plasma C-peptide levels.
Moreover, we show that targeted delivery of oestrogen using GLP-1 as a peptide carrier and intensive insulin co-therapy by a distinct MOA alleviates hyperglycaemia, increases fasting C-peptide levels and redifferentiates β cells, whereas it reduces daily insulin requirements by 60% and limits weight gain in mice. The enhanced restoration of GLP-1R expression in dedifferentiated β cells by the GLP-1-oestrogen and insulin cotreatment renders them susceptible to targeted delivery of oestrogen. As has been proposed previously in the Akita mouse model 46 , we observed that stimulating the ERAD pathway by GLP-1-oestrogen beneficially influences β-cell physiology in rodent models of diabetes. In future studies, it might be of specific interest to test GLP-1-oestrogen with and without PEG-insulin co-therapy in genetically perturbed mouse models of ERAD. Finan et al. showed previously that peptide-based targeting prevented adverse side effects of oestrogen, such as uterus and tumour growth 19 . There was also no measurable oestrogen-induced increase in bone content owing to the effect of the absence or limited expression of GLP-1R on off-target tissues and cells 19 . Here, we extended the safety profile of GLP-1-oestrogen and demonstrated that GLP-1-oestrogen did not stimulate uterine tissue growth in OVX rats. This study further verified that there is insufficient free, systemic oestrogen to drive toxicity as well as the selectivity and specificity of GLP-1-mediated oestrogen targeting. The strategy to use GLP-1 as a carrier may be adopted to selectively target any other small molecule or biologic to β cells. The prerequisite for the transport of the molecule of interest into the target cell (that is, stressed and dedifferentiated β cells) is adequate GLP-1R expression. Under hyperglycaemic conditions, adjunctive treatments that reduce the glycaemic burden, such as chronic insulin therapy as demonstrated here, can facilitate the restoration of GLP-1R expression in stressed β cells. Notably, chronic PEG-insulin treatment also increased functional β-cell number and the scRNA-seq data suggest a direct effect of insulin on β cells. Thus, combinatorial pharmacological treatments that include insulin might have additional beneficial effects on β-cell survival, protection, proliferation and function. In summary, our work has identified mechanisms and pathways of β-cell dedifferentiation and opens avenues for pharmacological targeting of these dedifferentiated cells for diabetes remission.
Methods mSTZ treatment. mSTZ (Sigma-Aldrich cat. no. S0130) was injected intraperitoneally in 8-week-old male C57BLJ/6 mice (n = 125) at 50 mg kg −1 for five consecutive days following the mSTZ model to induce diabetes. A subset of agematched male mice was injected with ice-cold citrate buffer (pH 4.5) as control animals (n = 20). C57BLJ/6 mice were obtained from Janvier Labs. Ten days after the last mSTZ injection, fasting blood glucose was taken, as well as fasting plasma, to determine fasting insulin and C-peptide levels. We included hyperglycaemic mice with fasting blood glucose levels of >190 mg dl −1 (n = 116). We estimated β-cell function and mass of mSTZ-treated mice by combining fasting blood glucose levels, the homeostatic model assessment (HOMA)-β-score and the ratio of fasting C-peptide to blood glucose levels. Among mSTZ-treated mice, animals with fasting glucose levels of >25th percentile, a HOMA-β-score of <25th percentile and a C-peptide to blood glucose ratio of <25th percentile were excluded from the study (n = 9).
Pharmacological study in mSTZ mice. mSTZ diabetic mice were randomised and evenly distributed to different treatments according to fasting blood glucose levels. Ten days after the last STZ injection, mice were allocated to different treatments of daily subcutaneous injection with vehicle (PBS; n = 17, not mSTZ-treated), vehicle (PBS; n = 17, mSTZ-treated), a GLP-1 analogue (n = 16), oestrogen (n = 14), GLP-1-oestrogen (n = 28, of which n = 11 mice were switched to vehicle (PBS) treatment after 12 weeks of GLP-1-oestrogen treatment), PEG-insulin (n = 13), or GLP-1-oestrogen and PEG-insulin (n = 16) at the doses indicated in Fig. 1a for 100 d. Mice were housed up to four per cage on a 12:12-h light-dark cycle at 22 °C and given free access to a normal chow diet (Altromin 1314) and water. Compounds were administered in a vehicle of PBS (Gibco) and were given by daily subcutaneous injections at the indicated doses at a volume of 5 μl g −1 body weight for the indicated durations as indicated in the figure legends. The investigators were not blinded to group allocation during the in vivo experiments or to the assessment of longitudinal endpoints. All rodent studies were approved by and performed according to the guidelines of the Animal Use and Care Committee of Bavaria, Germany.
Study in FVFPBF DHom mice. Male 8-week old FVFPBF DHom mice with fasting blood glucose > 250 mg dl −1 were randomised to vehicle (n = 7), oestrogen (n = 5), GLP-1 (n = 9) or GLP-1-oestrogen (n = 11) treatment according to their fasting blood glucose levels. FVFPBF DHom mice were obtained from breeding. Mice were treated daily for four weeks with subcutaneous injections. Fasting blood glucose was measured after a 6-h fast. We single-or group-housed the mice on a 12-h light, 12-h dark cycle at 22 °C and provided free access to food and water. This study was approved by and performed according to the guidelines of the Animal Use and Care Committee of Bavaria, Germany.
Uterotrophic assessment in OVX rats. The study was designed in accordance with the Endocrine Disruptor Screening Programme Test Guidelines OPPTS 890.1600: Uterotrophic Assay, a standardised in vivo screening test intended to evaluate the ability of a chemical to elicit biological activities consistent with agonists of natural oestrogens (for example, 17ß-oestradiol). The test measures the increase in uterine weight or uterotrophic response in comparison with non-treated controls 61 . The study was further designed according to accepted pharmacological principles and followed Good Laboratory Practice (conducted by Envigo). A total of 44 OVX female Sprague-Dawley rats (Charles River) were supplied for the study, of which 40 animals were allocated to treatment groups (randomised by body weight to ensure equal group mean starting body weight) and the remaining 4 animals were allocated as spares. On the day of dosing (following 14-22 d of acclimatisation), rats were approximately 9-12 weeks of age, and weighed 217 g to 348 g. Four groups of 8 rats each were treated for 14 consecutive days with  To synthesise the derivatized oestrogen for construction of the conjugate, we reacted oestradiol 17-acetate (Sigma) and a tenfold excess of ethyl 2-bromoacetate in dioxane in the presence of K 2 CO 3 under reflux conditions and agitation. After the removal of dioxane in vacuo, we resuspended the intermediate product in dioxane with 1 N NaOH with heat, followed by a re-acidification by HCl in dichloromethane. After the removal of DCM in vacuo, the crude product was resuspended in aqueous solvent containing at least 20% MeOH, 20% acetronitrile (ACN, VWR Chemicals) and 1% AcOH. The crude extract was purified by reversed-phase HPLC.
The derivatized oestrogen was added covalently to the side-chain amine on the C-terminal lysine amide residue. We used a C-terminal N-methyltrityl-l-lysine (Lys(Mtt)-OH) residue, whose side chain was orthogonally deprotected with four 10-min treatments with 1.5% TFA, 2% TIS and 1% anisole in DMF. The oestrogen was attached through this side-chain amine of the C-terminal lysine after treatment with a threefold excess of the purified derivatized oestrogen and DIC/HOBt-Cl in N-methyl-2-pyrrolidone. After TFA cleavage as described above, the crude extract was resuspended in aqueous buffer containing 20% ACN and 0.1 M NH 4 HCO 3 . The conjugate was purified as described above.

PEG-insulin. Pegylated insulin (PEG-insulin) was prepared by insulin N-terminal
amine reductive amination with 20 K methoxy PEG propionaldehyde. In brief, human insulin was dissolved in 50 mM sodium acetate buffer (pH 5.0) and 50% ACN. A 30-fold excess of sodium cyanoborohydride and a 1.5-fold excess of methoxy PEG propionaldehyde (M-ALD-20K, JenKem Technology) was added to the buffer containing insulin for 3 h at room temperature with stirring. Purification by reverse phase chromatography on a C-8 column in 0.1% TFA acetonitrile solvents yielded the final product at greater than 95% purity.
Oestrogen (17β-oestradiol, Sigma-Aldrich) was dissolved in 100% ethanol (Sigma-Aldrich) at a concentration of 1 mg ml −1 and was diluted with PBS to the required concentration.
Blood parameters. Blood was collected from tail veins after a 4-h fast, using EDTAcoated microvette tubes (Sarstedt). Blood was immediately chilled on ice. Plasma was separated by centrifugation at 5,000g at 4 °C for 10 min using a micro centrifuge and was stored at −20 °C until further usage. Plasma insulin and C-peptide (Crystal Chem) and Proinsulin (Alpco) were quantified by ELISA (enzyme-linked immunosorbent assay) following the manufacturer's instructions. Four-hour fasting blood glucose levels were determined using a handheld glucometer (FreeStyle).
Pancreas dissection. Adult pancreata were dissected and fixed in 4% PFA in PBS for 24 h at 4 °C. The tissues were cryoprotected in a sequential gradient of 7.5%, 15% and 30% sucrose-PBS solutions at room temperature (2 h incubation for each solution). Next, the pancreata were incubated in 30% sucrose and tissue embedding medium (Leica) (1:1) at 4 °C overnight. Afterwards, they were embedded in a cryoblock using tissue-freezing medium (Leica), frozen in dry ice and stored at −80 °C. Sections of 20 μm thickness were cut from each sample, mounted on a glass slide (Thermo Fisher Scientific) and dried for 10 min at room temperature before use or storage at −20 °C.
Automatic tissue analysis. Stained tissue sections were scanned with an AxioScan. Z1 digital slide scanner (Zeiss) equipped with a ×20 magnification objective. We scanned three sections per animal. Images were evaluated using the commercially available image analysis software Definiens Developer XD 2 (Definiens). First, regions of interest were annotated manually to select islets of Langerhans for analysis. A specific rule set was then defined to detect and quantify the cells within each defined region on the basis of the fluorescence intensity of DAPI, morphology, size and neighbourhood. The Ins-, Gcg-or Sst-expressing cells were classified automatically using the fluorescence intensity of each hormone.
EdU detection protocol. EdU staining was carried out according to the EdU imaging kit manual (Life Technologies) after staining with the secondary antibody. DAPI staining and mounting was performed as described above.
Pancreatic insulin content. Pancreatic insulin content was determined by an acid ethanol extraction. The pancreas was dissected, washed in 1× PBS and homogenised in an acid-ethanol solution (5 ml 1.5% HCl in 70% ethanol) followed by incubation at −20 °C for 24 h. After 2 rounds of acid-ethanol precipitation, the tissue was centrifuged (2000 r.p.m., 15 min, 4 °C) and the supernatant was neutralised with 1 M Tris pH 7.5. Insulin was measured using a mouse insulin ELISA (Crystal Chem) and was normalised over the protein concentration that was determined by BCA protein assay.
Islet isolation. Islet isolation was performed by collagenase P (Roche) digestion of the adult pancreas. In brief, 3 ml of collagenase P (1 mg ml −1 ) was injected into the bile duct and the perfused pancreas was consequently dissected and placed into another 3 ml of collagenase P for 15 min at 37 °C. Then, 10 ml of G-solution (HBSS (Lonza) + 1% BSA (Sigma-Aldrich)) was added to the samples, and this was followed by centrifugation at 1,600 r.p.m. at 4 °C. After another washing step with G-solution, the pellets were re-suspended in 5.5 ml of gradient preparation (5 ml 10% RPMI (Lonza) + 3 ml of 40% Optiprep (Sigma-Aldrich) per sample), and placed on top of 2.5 ml of the same solution. To form a 3-layer gradient, 6 ml of G-solution was added on the top. Samples were then incubated for 10 min at room temperature before centrifugation at 1,700 r.p.m. Finally, the interphase between the upper and the middle layers of the gradient was harvested and was filtered through a 70 μm Nylon filter then washed with G-solution. Islets were hand picked under the microscope.
Single-cell suspension. To achieve a single-cell suspension of islets, islets were hand picked into a 1.5 ml Eppendorf tube, pelleted (800 r.p.m., 1 min) washed with PBS (minus Mg or Ca, Gibco) and digested with 0.25% trypsin with EDTA (Gibco) at 37 °C for 8 min. Mechanical disaggregation every 2-3 min was required. The digestive reaction was then stopped and cells were pelleted (1200 r.p.m., 5 min).
Single-cell sequencing. Single-cell libraries were generated using the Chromium Single-cell 3′ library and gel bead kit v2 (PN 120237) from 10x Genomics. In brief, to reach a target cell number of 10,000 cells per sample 16,000 cells per sample were loaded onto a channel of the 10x chip to produce Gel Bead-in-Emulsions (GEMs). This underwent reverse transcription to barcode RNA before cleanup and cDNA amplification followed by enzymatic fragmentation and 5′ adaptor and sample index attachment. Libraries were sequenced on the HiSeq4000 (Illumina) with 150 bp paired-end sequencing of read2.
FACS sorting. FACS sorting of endocrine cells was performed using the FACS-Aria III (BD Biosciences). Single cells were gated according to their FSC-A (front scatter area) and SSC-A (side scatter area). Singlets were gated dependent on the FSC-W (front scatter width) and FSC-H (front scatter height) and dead cells were excluded using the marker 7AAD (eBioscience). The FVF endocrine populations were discriminated according to their Venus fluorescence emission at 488 nm and the βand α-lineages were discriminated according to their BFP emission at 405 nm (positive and negative respectively). To enrich for β cells the distinct SSC-A high populations were gated. FACS sorted cells were sorted directly into Qiazol (Qiagen) for RNA isolation. Compound, cytokine treatments and GSIS with human micro-islets. Dilution series of compounds were performed in 3D InSight Human Islet Maintenance Medium. Each of the assessed compounds was added to the culture medium at the concentration indicated below, one day prior to the start of the cytokine treatment. The cytokine cocktail contained: tumour necrosis factor alpha (TNFα, 10 ng ml −1 , Thermo Fisher Scientific PHC3016), interferon gamma (IFNγ, 10 ng ml −1 , Sigma-Aldrich I3265) and interleukin-1beta (IL-1β, 2 ng ml −1 , Sigma-Aldrich I17001), was prepared in PBS containing 0.1% BSA (Sigma-Aldrich A7888). The same concentration of PBS-BSA solution was maintained in each experimental condition. Regular redosings with cytokines and compounds were performed every 2-3 d. Prior to GSIS, culture medium was removed and islet microtissues were washed twice with Krebs Ringer Hepes Buffer (KRHB; 131 mM NaCl, 4.8 mM KCl, 1.3 mM CaCl 2 , 25 mM Hepes, 1.2 mMm KH 2 PO 4 , 1.2 mM MgSO 4 , 0.5% BSA) containing 2.8 mM glucose and equilibrated for 1 h in the same solution. GSIS was performed in KRHB containing indicated glucose concentrations for 2 h. Following GSIS, the tissues were lysed using the CellTiter-Glo Luminescent Cell Viability Assay (Promega G9241) with protease inhibitor cocktail (Promega G6521) and the luminescence was recorded with a microplate reader (Infinite M1000, TECAN) for the analysis of total ATP content. The lysates were then used for assessment of total insulin content. After the correct dilutions in KRHB were performed, total and secreted insulin was quantified using the Stellux Chemi Human Insulin ELISA (Alpco, 80-INSHU-CH10). The Caspase-Glo 3/7 Assay (Promega, G8090) was used to assess caspase-3 and caspase-7 activity in the islet microtissues following compound treatment.
Statistical analysis not including scRNA-seq data. Preliminary data processing and calculations during ongoing studies were carried out using Microsoft Excel 2016. All further statistical analyses were performed using GraphPad Prism 8. We used the one-way analysis of variances (ANOVA) followed by Tukey's post hoc analysis to determine significance among the different treatment groups. In the case of only two groups, the unpaired Student two-tailed t-test was used to detect significant differences. The human micro-islets were derived from three different donors, each of whom naturally varied in their GSIS. To compare the treatment effects among all donors, we used one-way ANOVA with the different donors as random effect followed by Tukey's post hoc analysis. This analysis was performed in R. Grubbs' test (α < 0.05) was used to detect significant outliers, which were then excluded from subsequent statistical analysis and figure drawing. P < 0.05 was considered to be statistically significant. All results are mean ± s.e.m. unless otherwise indicated.
Preprocessing of droplet-based scRNA-seq data. Demultiplexing of binary base call (BCL) files, alignment, read filtering, barcode and unique molecular identifier (UMI) counting were performed using the CellRanger analysis pipeline (v2.0.0) provided by 10x Genomics. High quality barcodes were selected on the basis of the overall distribution of total UMI counts per cell using the standard CellRanger cell detection algorithm. All further analyses were run with python3 using the scanpy package 62 (v1.0.4 + 92.g9a754bb, https://github.com/theislab/scanpy) except stated otherwise. Genes with expression in less than ten cells were excluded. Furthermore, as is also applied as standard preprocessing steps in scanpy tutorials, low quality or outlier cells were removed if they: (1) had a high fraction of counts from mitochondrial genes (40% or more); (2) expressed more than 7,000 genes; or (3) had more than 100,000 UMI counts. Cell-by-gene count matrices of all samples were then concatenated to a single matrix. To account for differences in sequencing depth, UMI counts of each cell were normalised by total counts of that cell (pp. normalize_per_cell with mean = TRUE) and values were log-transformed. Highly variable genes (n = 1,625) were selected on the basis of normalised dispersion using the setting the lower cutoffs for the mean to 0.0125 and for the dispersion to 0.5. This matrix was used as input for all further analyses unless otherwise indicated.
Embedding, clustering and cell type annotation. Clustering was performed on the full data set to reduce systematic biases such as batch effects as was recommended in ref. 63 . A single-cell neighbourhood graph (kNN-graph) was computed on the 50 first principal components using 15 neighbours. To minimise condition effects and to facilitate clustering we recomputed the kNN-graph using the first 15 diffusion components of the PCA-based graph as suggested in ref. 64 . For clustering and cell type annotation, Louvain-based clustering 65 was used as implemented in louvain-igraph (v0.6.1 https://github.com/vtraag/louvain-igraph) and adopted by scanpy (tl.louvain). The resolution parameter was varied in different parts of the data manifold to account for strong changes in resolution (for details, see Data availability). Clusters were annotated on the basis of the mRNA expression of the four main hormone genes Ins1 and Ins2, Gcg, Sst and Ppy (endocrine cells) and other known marker genes (non-endocrine cells) and were merged if they reflected heterogeneity only in a cell type outside the focus of this study.
Ductal cells (that express Krt19), acinar cells (that express Prss2), endothelial cells (that express Plvap), stellate cells (that express Col1a2) and small clusters of potential doublet-like cells that co-express endocrine and non-endocrine markers were removed from further analysis. Immune cells (that express Cd74) that infiltrate the islets were finely subclustered into macrophages (that express Cd86, Adgre1 and Cd14), dendritic cells (that express Cd86, Itgax and Iftim3), B cells (that express Cd79a and Cd79b) and T cells (that express Cd8a and Cd3d).
The hormone genes Ins1 and Ins2, Gcg, Sst and Ppy were expressed at very high levels and showed background level expression in all other endocrine subtypes and non-endocrine cell types. Such background expression is a common phenomenon in droplet-based scRNA-seq data. It is commonly said to be due to free-floating mRNA in the single-cell solution that comes from lysed cells and that is incorporated into all droplets. For annotation, only hormone expression that was well above the background level in non-endocrine cells, such as ductal, immune and endothelial cells, was considered.
For the identification of β-cell substates a new kNN-graph on the first 50 principal components was calculated and put into Louvain-based clustering of both Ins monohormonal and the connected Ins-PP cells. Similarly, Ins-Sst cells were subclustered from Ins-Sst-PP cells after recalculating the kNN-graph on the first 50 principal components. Ins-Gcg-Sst cells were assigned using a manual threshold for all three hormones that was well above ambient levels (threshold = 6 for normalised data).
For visualisation, UMAP was run as recommended in ref. 66 . For each UMAPplot the UMAP was newly calculated by recomputing the kNN-graph on the represented cell subset using the first 50 principal components.

Identification of polyhormonal singlets and doublet-like endocrine cell clusters.
Polyhormonal cells have previously been reported to exist in pancreatic islets [67][68][69] . However, the expression of multiple hormones in the same droplet can also be an indication of a doublet. It can therefore be difficult to distinguish polyhormonal singlets from doublets. A doublet rate of ~8-10% was measured in experiments with the same concentration of cells using the 10x technology 70 . This rate includes doublets with contributions from two different cell types (here, polyhormonal doublets) and from the same cell type (here, monohormonal doublets). The monohormonal doublets resemble monohormonal singlets and do not affect subsequent analyses 71 . We calculated the expected doublet frequency of polyhormonal doublets for a doublet rate of 10% using the frequency of monohormonal cell types that contribute to the doublet (doublet contributors) and assuming that doublets are generated by sampling singlet cells uniformly at random 71 . In every sample the proportion of observed polyhormonal cells clearly exceeded the expected polyhormonal doublet frequency. Therefore, in our data set it is unlikely that all detected polyhormonal cells are doublets. Application of doublet cell detection tools Scrublet 71 (v0.1, https://github.com/AllonKleinLab/scrublet) and DoubletDetection (https://github.com/JonathanShor/DoubletDetection) failed to resolve which clusters represent doublets and which clusters represent polyhormonal singlets. Predictions of the tools disagreed with each other and the doublet rate was consistently overestimated. We therefore used the following criteria to evaluate polyhormonal cell clusters and to distinguish between singlets and doublets: (1) Doublets should not express unique genes. All genes should also be expressed in at least one doublet contributor.
(2) Doublets should express marker genes or lineage-determining transcription factors of the doublet contributors. Downregulation of these genes indicates singlet populations.
(3) Previous reporting of polyhormonal singlet cells in the literature.
(4) Clusters of polyhormonal cells with a higher frequency than expected by our doublet simulation indicate polyhormonal singlet clusters (Extended Data Fig. 6g).
On the basis of these criteria we found sufficient evidence for Ins-PP, Ins-Sst-PP, Gcg-PP (low) and Gcg-PP (high) cells to be polyhormonal singlets, but we excluded Ins-PP-Gcg, Ins-Gcg, Ins-Gcg-Sst, Gcg-Sst-PP and Sst-PP(high) cells as probable doublets.
Cell-cycle classification. To classify cells into cycling and non-cycling cells, first, a score was assigned to each cell for a set of S and G2/M phase genes 72 as proposed 73 , and second, all cells with an S-score or a G2/M-score > 0.25 were classified as cycling. The threshold was chosen on the basis of the score distribution. The score for a given gene set was computed as described 74 and implemented in scanpy (tl. score_genes_cell_cycle).

Marker genes of the main endocrine cell types.
For the characterisation of the four endocrine cell types, specific marker genes were identified by comparing the gene expression profile of each cell type against all cells of the other three cell types using a test with overestimated variance as implemented in scanpy (tl.rank_genes_ groups). All genes that ranked within the top 300 genes, had a test score of >8 and were unique markers for one cell type were considered as marker genes.
Differential expression testing to describe subpopulations and treatment responses. Differential expression testing between treatments and for the characterisation of immature β cells and polyhormonal subpopulations was performed on quantile-normalised (quantile threshold = 0.95) and log-transformed data to account for extremely highly expressed genes (for example, the main hormones in endocrine cells) that may incorrectly alleviate the expression of other genes in a cell when total count normalisation is applied. By quantile normalisation, each cell is normalised by the total UMI count in the cell, of genes that account for less than 5% of the total UMI counts across all cells. Thus, very highly expressed genes are not considered for normalisation. For differential expression testing we used limma-trend 75,76 , as implemented in the Bioconductor package limma (v3.28.10), through an rpy2 (v2.9.1) interface as recommended in ref. 77 . In each test only genes expressed in >1% of cells in any of the two subsets tested were considered. Gene set enrichment was performed using EnrichR 78 through its web interface. Genes with a false discovery rate (FDR) of <0.01 and an estimated log(fold change) (output from limma model, not the actual log(fold change), as log-transformed data were the input) of >0.25 were used as input. Notably, the hormone genes Ins1 and Ins2, Gcg, Sst and Ppy, as well as other known cell type marker genes Pyy, Iapp, Ttr, Gpx3, Ctrb1 and Try5 were also differentially expressed in other cell-types in which they are expressed only at background levels (free-floating mRNA, see Embedding, clustering and annotation). These genes are indicated in the Supplementary Tables and were excluded from plotting.
Identification of specific β-cell dedifferentiation markers. Genes specific for dedifferentiated β cells (β-mSTZ) were extracted from the list of all significantly upregulated genes (FDR < 0.05, estimated log(fold change) > 0.25) in β cells from mSTZ-treated mice compared to β cells of healthy mice by two filtering steps. First, non-specific genes that were also differentially expressed in any of the other monohormonal endocrine cell types (α, δ and PP) were excluded. Second, only genes that were expressed in at least 25% of β cells from mSTZ-treated mice and in less than 5% of β cells from healthy mice were considered. Location was extracted from the GeneCards database (https://www.genecards.org).

Inference of β-cell maturation, dedifferentiation and regeneration trajectories.
Pseudotime of β-cell maturation in healthy islets and dedifferentiation upon mSTZ treatment was calculated using diffusion pseudotime (dpt) 79,80 as implemented in scanpy (tl.dpt), selecting a random root cell within the starting population. The choice of root cell did not affect the inferred pseudotemporal ordering strongly. Similarly, the dpt approach was used to model β-cell regeneration and to estimate the location of treated β cells along the path from dedifferentiated to mature β cells.
Here, dpt was used as a cell-to-cell distance metric across samples. Cycling cells as well as a small subpopulation of β cells were excluded from visualisation as they were not part of the linear trajectory and they showed very high pseudotime values.
Comparison of β-cell dedifferentiation trajectory to embryonic and postnatal maturation. To compare the dedifferentiation trajectory to embryonic and postnatal maturation we used a publicly available single-cell RNA-seq data set as a reference that contained cells that were sorted from Gcg-Venus and Ins-GFP reporter mice at 6 different time points (E17.5, P0, P3, P9, P18 and P60) 37 . The filtered and annotated raw count matrix was downloaded from the Gene Expression Omnibus (GEO) (GEO accession number: GSE87375). The analysis was run using the updated scanpy package v1.4.4, as only this version includes the necessary data integration methods. ERCC RNA spike-in and genes with expression in fewer than three cells were excluded. The data were normalised to total counts per cell using the pp.normalize_total function in scanpy with default parameters and excluding highly expressed genes, and were log-transformed (pp. log1p). For the scope of this manuscript we used a subset of the data that contained only β cells (Ins1-positive cells). Therefore, we computed a kNN-graph on the 50 first principal components using 15 neighbours and performed a first round of Louvain-based clustering (tl.louvain). As input, data were subset to the 3,000 top ranking highly variable genes (pp.highly_variable_genes). We excluded all clusters that showed high expression of Gcg, which is indicative of α cells, or that showed high expression of Mki67, which is indicative of proliferative cells. In addition, we filtered cells that showed high expression of the δ cell markers Sst or Hhex.
We integrated this reference data set with the data from our study, which was subset to β cells from Ctrl and mSTZ-Vehicle treatment, using BBKNN 81 from the scanpy external package (sce), which allowed us to compute an embedding and trajectory that included both data sets. To reduce noise, we excluded genes that were expressed in fewer than 15 cells in each data set. In addition, cycling cells (Mki67 > 1, 33 reference cells and 30 of the cells from this study) were excluded, as their expression profile is dominated by the expression of cell-cycle genes and these cells therefore formed a separate cluster which was not part of the linear maturation trajectory. To ensure that β-cell maturation dominates the gene expression variation we first considered only the 2,000 top-ranked highly variable genes (pp.highly_variable_genes) of the reference data set that were also expressed in our data, and, second, we reduced the contribution of heterogeneity within the β1 and β-mSTZ cluster to gene expression variation by randomly subsetting both clusters to 500 cells. We then scaled and zero centred each data set separately (pp. scale) and concatenated the two data sets, which resulted in a 1,788 cells by 1,654 genes count matrix. We computed a common kNN-graph on the first ten principal components using the sce.pp.bbknn function with default parameters and k = 5 within batch neighbours. As the data sets did not show a strong batch effect even without integration, assessed by visual inspection of the first principal components and diffusion components, and thus transitions between cells from the two data sets also showed a high probability, we were able to use diffusion pseudotime (tl.dpt) to infer the maturation trajectory with the common kNN-graph as input. The trajectory was calculated by selecting a random root cell from the embryonic cells that were sampled at E17.5. The choice of root cell did not affect the inferred ordering strongly. The ordering of the cells from each data set was largely consistent with the ordering obtained prior to integration, assessed by the distribution along the trajectory of the time points or β-cell subgroups. To quantify the cluster similarity, PAGA was applied to the common kNN-graph (tl.paga).
To compute an embryonic or immaturity cell score and a maturity cell score we extracted the gene signatures from the reference data set. The reference data was subset to the 2,000 top ranked highly variable genes and Louvain-based clustering was performed on the kNN-graph that was computed on the first 50 principal components with 15 neighbours. The cluster that consisted of cells that were sampled at E17.5 and P0 was annotated as 'embryonic or neonatal' , whereas the cluster consisting mainly of cells that were sampled at P60 was annotated as 'mature' . Differentially expressed genes between these two clusters were used for scoring. Genes that were upregulated in the mature cluster or in the embryonic or neonatal cluster were used as a gene set for maturity or for embryonic or immaturity score, respectively. For differential expression testing the t-test with overestimated variance implemented in the tl.rank_genes_groups function of scanpy was used. The top 500 ranked genes with a log(foldchange) > 0.25 and an adjusted P value < 0.01 were considered. Cell scores were computed using the tl.score_genes function in scanpy.
Inference of cluster-to-cluster distances, lineage relations and cell movement. PAGA 64 was performed to infer cluster and lineage relations using the tl.paga function of scanpy with a threshold on edge significance of 0.05. In a PAGA graph, paths represent cluster connections or relations that indicate potential routes of cell transitions. Edge weights represent the confidence of a connection calculated on the basis of a measure of cluster connectivity.
To infer the direction of possible transitions 82 and cell movements we estimated RNA velocity using a stochastic version implemented in the scVelo python package (v0.1.16.dev13 + c1a6dad, https://github.com/theislab/scvelo with scanpy v1.3.2). Splicing information of reads was extracted using the velocyto pipeline (v0.17.7, http://velocyto.org). We then followed the recommended steps described in scVelo to estimate RNA velocities and RNA velocity force fields. First, data were preprocessed by filtering genes with less than 30 spliced or 30 unspliced counts and both unspliced and spliced counts were normalised by total counts. Then the first-and secondorder moments for each cell were computed across its 15 nearest neighbours of the kNN graph in PC space (50 PCs). Next, RNA velocities were estimated using a stochastic model of transcriptional dynamics. To obtain a more conservative estimate a 95% quantile fit was used. Finally, to project the velocity vector of each cell into the low-dimensional UMAP embedding for visualisation and interpretation, the expected mean direction given all potential cell transitions on the kNN graph was computed. Each potential cell transition is assigned a probability corresponding to the correlation to the predicted transition by the velocity vector (velocity graph). For example, a high probability corresponds to a high correlation with the velocity vector. The projection results in a low dimensional map of RNA velocity which indicates the predicted cell state transitions. For computation of the velocity graph and embedding only genes with an r 2 > 0.1 of the velocity fit were considered.
The velocities of each gene were calculated over all treatments except for healthy β cells, for which only healthy cells were used. A treatment can here be considered as a process by which cells move from the diseased cells potentially towards healthy cells, as for the pseudotime inference described above. During this process genes are induced and/or repressed, as approximated by RNA velocity. Therefore, to also take into account these intermediate gene states, all treatments were included for model fitting and velocity estimation. Both PAGA and the RNA velocity graph and projection were instead only computed on the represented cell subset. For this, the kNN-graph was recalculated for the cell subset using the first 50 principal components and the highly variable genes as initially defined.
Reporting Summary. Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Data availability
Custom python scripts written for performing scRNA-seq data analysis are available in a github repository (https://github.com/theislab/pancreas-targeted_ pharmacology). Versions of packages that might influence numerical results are indicated in the scripts. Raw data and gene expression matrices of scRNA-seq are deposited in GEO under the accession number GSE128565. Source data for