Endocrine islet β-cell subtypes with differential function are derived from biochemically distinct embryonic endocrine islet progenitors that are regulated by maternal nutrients

Endocrine islet b cells comprise heterogenous cell subsets. Yet when/how these subsets are produced and how stable they are remain unknown. Addressing these questions is important for preventing/curing diabetes, because lower numbers of b cells with better secretory function is a high risk of this disease. Using combinatorial cell lineage tracing, scRNA-seq, and DNA methylation analysis, we show here that embryonic islet progenitors with distinct gene expression and DNA methylation produce b-cell subtypes of different function and viability in adult mice. The subtype with better function is enriched for genes involved in vesicular production/trafficking, stress response, and Ca2+-secretion coupling, which further correspond to differential DNA methylation in putative enhancers of these genes. Maternal overnutrition, a major diabetes risk factor, reduces the proportion of endocrine progenitors of the b-cell subtype with better-function via deregulating DNA methyl transferase 3a. Intriguingly, the gene signature that defines mouse b-cell subtypes can reliably divide human cells into two sub-populations while the proportion of b cells with better-function is reduced in diabetic donors. The implication of these results is that modulating DNA methylation in islet progenitors using maternal food supplements can be explored to improve b-cell function in the prevention and therapy of diabetes.


Introduction
Endocrine islet b cells regulate whole-body glucose homeostasis via glucose-stimulated-insulin secretion (GSIS).Insu cient GSIS results in diabetes, a metabolic disease that a icts over 6% of the world population.Thus, understanding how each individual makes a su cient number of b cells with speci c levels of secretory function, known as functional b-cell mass, holds the key to prevent and to cure this disease.
The majority of b cells in adults are made through two pathways.During embryogenesis, islet progenitors that transiently express transcription factor (TF) Neurogenin3 (Neurog3 or Ngn3) differentiate into immature b cells (Gu et al., 2002).At late gestation and early postnatal stages, these cells proliferate and mature to make up the functional b-cell mass needed for physiology.Intriguingly, previous work has clearly established that adult b cells comprise cell subsets with differential levels of GSIS (Bader et al., 2016;Dorrell et al., 2016;Salomon and Meda, 1986), proliferation (Bader et al., 2016), metabolism (Pipeleers et al., 2017), and Ca 2+ in ux (Johnston et al., 2016).Yet when/how each b-cell subset is determined, how stable they are, and how altering the proportions of different b-cell subsets affect the risk of diabetes have not been directly tested.
In contrast, some reported b cell subsets are relatively stable.For example, a few "virgin b cells" that arise via postnatal a-ell transdifferentiation can retain their immaturity for months (Lee et al., 2021).
Additionally, two recent studies suggest that differential DNA methylation in several enhancers and/or global histone methylation (H3K27me3) can mark relatively stable b-cell subtypes, whose proportions can be modi ed by the dosages of genes that render these modi cations (Dror et al., 2023;Parveen et al., 2023).However, virgin b cells represent only a few percentage of all b cells (van der Meulen et al., 2017).Metabolic signaling have been shown to alter the DNA and histone methylation states in b cells (Avrahami and Kaestner, 2019;Avrahami et al., 2015;Brasacchio et al., 2009;Hall et al., 2018;Hall et al., 2019;Wortham et al., 2023).Thus, questions regarding the long-term stability of large b-cell subsets, in addition to stage when these heterogeneities are induced, remain to be answered.
Besides the issues on b-cell subset stability and origin, the causal relations between b-cell heterogeneity and diabetes are not known.Several studies have shown that the proportions of b-cell subsets with better-functionality are usually reduced in islets of diabetic donors (Dorrell et Weng et al., 2023).Yet it is not clear if b-cell subsets interconvert under diabetes-related metabolic stress due to b-cell plasticity (Zhang et al., 2022), or that b cells with better function were made at reduced proportions that are promoted by diabetes-predisposing factors.Addressing these possibilities is important, because environmental factors contribute to ~half of all human diabetes with unknown mechanisms (Laakso and Fernandes Silva, 2022).
Our previous work revealed that Neurog3 + islet progenitors are heterogenous in respect to their expression levels of TF Myt1, a Neurog3 target that is detected in all differentiated islet cells (Liu et al., 2019).The result is the co-presence of transient Myt1 + Ngn3 + (M + N + ) and Myt1 -Ngn3 + (M -N + ) cells.The M + N + progenitors have higher levels of DNA methylation in a DNA enhancer of Arx, a key fate-specifying gene for a cells (Dhawan et al., 2011).They favor b-cell fate while the M -N + progenitors a fate, although they both can give rise to b cells (Liu et al., 2019).Because enhancer DNA methylation is also important in b-cell maturation and identity (Dhawan et al., 2011;Dhawan et al., 2015;Liu et al., 2019;Papizan et al., 2011), we postulated that M + N + and M -N + progenitors give rise to b-cell subtypes with differential enhancer methylation associated with their heterogeneity in gene expression and function
Pseudo-islets were made from FACS-puri ed tdT + and tdT -islet cells that include all endocrine cell types and assayed for GSIS (Fig. 1D, E).Two-months old tdT + pseudo-islets secrete signi cantly higher levels of insulin than tdT -samples in response to 20 mM glucose but not to 30 mM KCl, from both male (Fig. 1F) and female (Fig. 1G) samples.Similar results were observed in 8-months-old tdT + and tdT -cells (Fig. 1H).In addition, pseudo-islets made from pure b cell subsets of 2-months-old Myt1 cCre ; Neurog3 nCre ; Rosa26 eYFP ; Ins2 Apple mice (noted as MNYI) showed similar results, with descendant b cells of M + N + progenitors (eYFP + ) having higher GISS compared to eYFP -cells (Fig. S1A -C).In MNYI mice, the Rosa26 eYFP allele labels the descendants of M + N + progenitors with eYFP production, whereas Ins2 Apple labels b cells with uorescent protein Apple (Stancill et al., 2019), enabling isolation of pure b-cell subsets.
When producing pseudo-islets using MNYI b cells, we noticed that those from eYFP -cells are usually smaller than from eYFP + cells, although the same number of starting cells were used (Fig. 1I, J).We therefore tested if these b-cell subtypes have different viability in vitro.After 4 -7 days in culture under 11 mM glucose, signi cantly more eYFP + pseudo-islets/islet cells survive compared to eYFP -samples (Fig. 1K-M).Similarly, the co-presence of high glucose and palmitate, a condition that was widely used to mimic metabolic stress in T2D, leads to more pseudo-islet cell death in eYFP -than in eYFP + b cells (Fig. 1N, O).
We last compared the mitochondrial activity in these b-cell subsets.The mitochondria of adult eYFP + b cells have a statistically signi cant but only slightly higher (~6%) level of transmembrane potential (p=0.011) compared to eYFP -cells (Fig. S1D, E).But they have similar levels of mitochondrial mass and ADP/ATP ratios (Fig. S1E, F).The signi cance of these ndings are not known and not further pursued in this study.
These above results suggest that different islet progenitors give rise to stable b-cell subtypes with different ability to proliferate, function, and survive, leading us to determine the genetic/epigenetic basis for the production and activity of the two b-cell subtypes.
Neonatal b-cell subtypes have different transcriptional pro les that cannot be attributed to asynchronous maturation scRNA-seq was coupled with lineage tracing to compare the transcriptomes of b-cell subtypes derived from the M + N + and M -N + progenitors.Islet cells from P2 MNA mice were used for inDrops sequencing.Male and female samples were combined because no sex-based difference in b cell function has been reported at this stage.After ltering out low quality droplets and cells (Fig. S2A, B), we identi ed 1,275 tdT + and 2,183 tdT -single b-cell transcriptomes from four batches of biological replicates (Fig. 2A, Table S1).There are a total of 919 differentially expressed genes (DEGs) between these cell subtypes, including 362 down-and 557 up-regulated in the tdT + cells (Table S2).Gene ontology (GO) analysis showed that the upregulated genes comprise ER, Golgi function, Response to ER stress, Calcium ion transport, and Apoptosis (Fig. 2B).Those downregulated associate with Electron transport, Oxidative phosphorylation, DNA replication, and Cell cycle (Fig. 2C).All these processes are known to regulate b-cell proliferation, function, and death, corresponding to thei detected functional differences between these b-cell subtypes.
Intriguingly, supervised analyses showed that the DEGs include neither key maturation-de ning gene such as Cfap126, HK1, Mafa, NeuroD, Nkx6.1, Pdx1, Slc2a2, Syt4, or Ucn3, nor important aging markers such as Igf1r and Trp53BP1 (Fig. 2D).Instead, the DEGs include cell cycle regulators (Anapc7, Nek4, and Ccng1), cell death genes (Bax and Bcl2l1),K + or Ca 2+ channel genes (Cacna1a, Cacna2d1, and Kcnk3), and vesicular Ca 2+ sensors (Syt7 and Syt16) (Fig. 2E).In addition, we found that DNMT3a, the de novo DNA methylation enzyme, has higher expression level in the tdT + b cells, consistent with our report that a key difference between the M + N + and M -N + progenitors is their differing DNA methylation patterns (Liu et al., 2019).These ndings suggest that the DEGs cannot be attributed to asynchronous maturation or aging of b-cell subsets.Rather, they likely arise from differential genetic/epigenetic program(s) in their respective progenitors, which propagate in postnatal stages to de ne stable b-cell subtypes.We therefore compared the expression pro les of adult tdT + and tdT -b cells.
The adult b cell subtypes from M + N + or M -N + progenitors maintain differential gene expressions scRNA-seq was used to compare the transcriptomes of P60 tdT + and tdT -b cells in MNA mice.The male and female samples were initially processed separately, because b-cell dimorphism has been established at this stage.We obtained 1,319 tdT + and 1,164 tdT -high quality single b-cell transcriptomes in three male islet samples, in addition to 1,639 tdT + and 1,230 tdT -b cells in three female replicas (Fig. 3A.Fig. S2C, D. Table S1).When male and female samples were analyzed separately, both showed similar pathways that are differentially expressed between the tdT + and tdT -b cells (Table S3).We therefore analyzed and presented results after combining the male and female samples.Note that the portion of tdT + b cells is higher in P60 than P2 (Fig. 3B), consistent with the higher proliferation rate of tdT + b cells in early postnatal stages (Fig. 1A-C).
The P60 tdT + and tdT -b cells have 1,867 DEGs, with the majority (1,614) upregulated in the tdT + cells (Table S4).Consistent with the higher GSIS of tdT + b cells, the down-regulated genes in these cells are enriched in several processes that inhibit GSIS, including Actin-binding (Kalwat and Thurmond, 2013) and Cellular senescence (Aguayo-Mazzucato et al., 2017; Thompson et al., 2019) (Fig. 3C); while several upregulated processes are known to promote GSIS, including Protein translation/transport, Insulin secretion, ER chaperone complex, Mitochondrion, Calcium ion-regulated exocytosis, and cAMP-dependent protein kinase complex (Fig. 3D).Notably, the DEGs include down-regulation of CD9 and up-regulation of CD24a, Gpx3, and Rbp4 (Fig. 3E-I We next examined if any of the DEGs between P2 b-cell subtypes are preferentially retained at P60.Among the 919 P2 DEGs, 193 were retained at P60 (Table S4), a 2.5-fold enrichment (p=2.37E-33,hypergeometric analysis).Most of these 193 genes are upregulated in P60 tdT + b cells and they are enriched for ER/Golgi function, vesicle production/secretion, and quality controls (Fig. 3I).These ndings suggest that despite the dynamic gene expression changes during postnatal b-cell maturation, a substantial number of genes maintain their differential expression between the b-cell subtypes from newly-born to adult stages.We therefore tested the functional contribution of some DEGs in b-cell subtypes.
Several DEGs may contribute to differential GSIS in tdT + and tdT -b cells For functional tests, we chose the families of Myt TF (Myt1, Myt1L, and St18), synaptotagmins, and voltage-gated calcium channels (Cacna1a, Cacna1c, and Cacna2d1).Multiple members in each of these families display upregulation in the more functional P60 b-cell subtype (Fig. 4A, Table S4 Myt1 F/+ ; Myt1l F/+ ; St18 F/+ ; Pdx1 Cre mice were derived.This reduces the levels of Myt1, Myt1L, and ST18 proteins by 34-63% in islets compared with controls (Fig. S3B, C).These levels of reduction did not impact the overall morphology of islets (Fig. 4B), but compromised glucose clearance in both male and female mice (Fig. 4C, D), accompanied by reduced islet GSIS (Fig. 4E).Note that mutant males displayed glucose intolerance earlier than females, results that usually occur when islet function was partially compromised (Gannon et al., 2018).
To test the effect of reduced Syts expression in b-cells, CRISPRi was used to reduce the level of Syt7 by ~50% in mouse islets (Fig. 4F).This was achieved using two Syt7 guide RNAs to recruit a fusion protein of dCas9 with transcriptional repressor domain KRAB to the transcription starting sites of Syt7 (Fig. S3A) (Jensen et al., 2021).Reduced Syt7 expression did not alter the islet morphology (Fig. 4G) but compromised GSIS in male islets (Fig. 4H).Note that this level of Syt7 reduction did not compromise GSIS in female islets (Fig. 4I), consistent with the ndings that male b cells are more liable for dysfunction (Gannon et al., 2018).
For testing the contribution of Ca 2+ channels, we directly compared glucose-induced Ca 2+ in ux in tdT + and tdT -b cells in real time.The cytoplasmic Ca 2+ in ux in response to 11 mM glucose, assayed with bcell speci c GCaMP6 as a surrogate (Fig. S3B), is more robust in the tdT + b cells than in tdT -cells (Fig. 4J, K, S3C, Movie S1).These results, combined with known roles of voltage-gated channels in Ca 2+ in ux/secretion coupling and the identical ADP/ATP ratio in the two b-cell subtypes (Fig. S1F), suggest that higher Ca 2+ channel expression in tdT + cells contributes to their higher GSIS as well.We conclude that, multiple factors, by virtue of their differential expression, contribute to the different GSIS of b-cell subtypes.We next explored the mechanisms that establish this gene expression heterogeneity in b cells.
The neonatal b-cell subtypes from M + N + or M -N + progenitors display methylation differences at a subset of putative enhancers Work from our lab and others has shown that cell-type speci c transcriptional enhancers are selectively hypomethylated during progenitor differentiation (Barnett et al., 2020).Many hypomethylated regions (HMRs) are maintained in subsequent developmental stages, resulting in patterns that record the histories of cell fate decisions (Scott et al., 2023).Corresponding to these ndings, we now detected higher levels of DNMT3a in tdT + b cells at both P2 and P60 (Tables S2 and S4).In addition, acute inhibition of DNMTs in embryonic pancreatic buds signi cantly reduced the proportions of M + N + progenitors (Fig. 5A, B).Thus, we postulated that differential enhancer methylation in islet progenitors drive cell-subtype production and maintenance.
Whole genome bisul te sequencing (WGBS) was done on sorted eYFP + and eYFP -b cells from P2 MNYI males and female islets.This identi ed 78,634 HMRs, falling into 6 clusters based on their methylation level changes between the cell subtypes (Fig. S4A).Further analysis de ned 4,091 differentially methylated regions (DMRs) by gating those with a minimum length of 50 bp and 2 differentially methylated CpG sites (p<0.05).Among these DMRs, 1,905 have lower levels of methylation in eYFP -cells and 2,186 lower in eYFP + cells (Fig. 5C.Table S5).Each DMR contains 2-15 differentially methylated CpGs, spreading up to 3kb (Fig. S4B, C).Roughly 60% of the DMRs localize within the putative enhancer regions (from 2 kb upstream of the transcription starting site to the 3' UTR, Fig. 5D), while 80% within 100 kb, of their putative target genes (Fig. S4D).Note that the mean-difference in methylation of these DMRs is ~1.6-folds, lower than the usually larger (>2-folds) differences observed between diverse cell types and developmental timepoints (Fig. 5C) (Barnett et al., 2020;Liu et al., 2019;Scott et al., 2023).This is expected because of the higher levels of similarity between b subtypes.Also note that we detected DMRs in the Arx locus in these subtypes (Fig. 5E), which is differentially methylated in the M + N + and M -N + islet progenitors (Liu et al., 2019).This result is consistent with our hypothesis that the differential DNA methylation patterns are passed down do differentiated islet cells.).The implication is that b-cell subtypes may express same levels of these TFs (Table S2), but they can transcribe different levels of their target genes due to the presence of these DMRs, leading to differential function.Supporting this possibility, genes associated with P2 DMRs (3,143 and 2,871 corresponding to lower and higher levels of methylation in eYFP + cells, respectively, Table S5) are enriched for pathways that regulate Golgi, Wnt Signaling, Ion transport, Circadian entrainment, Ca 2+ transport, and Actin binding (Fig. 5G, H).Note that 902 genes are associated with DMRs with lower methylation while also with DMRs higher methylation in eYFP + cells (Table S5).This may suggest some coordinated regulation of these genes by multiple gene regulatory elements.

The DMRs between P2 b-cell subtypes associate with DEGs of both neonatal and adult b-cell subtypes
We next asked if any P2 DMR-associated genes are differentially expressed between the two b-cell subtypes at P2 and P60.Among the 919 DEGs at P2, 150 and 152 have DMRs with lower or higher levels of methylation in eYFP + cells, respectively (Table S2), signi cant enrichments over non-differentially expressed genes (p=0.0004,Fisher Exact test).The ratios of down-and up-regulated genes are similar in those associated with lower-or higher-methylated DMRs (Fig. 5I, p=0.4.Fisher Exact test), consistent with the idea that changes in enhancer DNA methylation can be both activating and repressing (Dhar et al., 2021).Corroborating this conclusion, ectopic methylation of a DMR in the rst intron of Syt7 (Fig. 5J) using a proven dCas9-DNMT3a construct (Liu et al., 2019) signi cantly reduces Syt7 expression in a MIN6 b-cell line (Fig. 5K, M).Note that this DMR was chosen because Syt7 expression levels can impact b-cell function and the methylation level in this DMR changes during b-cell aging, implying its functional signi cance (Avrahami et al., 2015).
Similar to P2 gene expression, 566 of the 1,867 P60 DEGs are associated with the DMRs at P2 (Table S4.P<0.0001, Fisher Exact test).These data support an idea that DMRs at P2 can impact the gene expression at later stages, likely by maintaining their DNA methylation landscape and/or promoting stable gene networks.We therefore tested how DNA methylomes evolve in the two b-cell subtypes from newly-born to adult stages.DNA methylomes in postnatal b cells are highly dynamic but adult b-cell subtypes maintain DMRs that associate with their DEGs WGBS in P60 b-cell subtypes identi ed 68,485 HMRs (Fig. S5A), 59,327 of which overlap with that at P2.They constitute 6 clusters, with three (a total of 31,698 HMRs) showing >1.8 fold methylation changes between P2 and P60 total-b cells (Fig. 6A), underscoring the dynamic DNA methylation changes from neonatal to adult stages.In contrast, the P60 b-cell subtypes have relatively small DNA methylation differences, with 1,113 and 1,253 DMRs that have higher or lower methylation in eYFP + b cells compared to eYFP -cells, respectively (Fig. 6B).These DMRs have 2-14 CpG dinucleotides differentially methylated, spanning 50 bp to 3.5 kb (Fig. S5-B, C), with the majority (>70%) locate in the putative regulatory regions of their predicted target genes (Fig. S5D).
There are 206 and 173 DEGs between P60 b-cell subtypes that are associated with P60 DMRs with lower or higher levels of methylation in eYFP + cells, respectively (Table S4 and Fig. 6D), a signi cant enrichment (p<0.0001,Fisher Exact test).At this stage, signi cantly more up-regulated DEGs are associated with DMRs with lower levels of methylation (p=0.036,Fisher exact test.Fig. 6D).The DEGs that are associated with lower levels of DNA methylation are enriched for Organism development, Organ morphogenesis, Growth hormone synthesis, secretion, and action, Anchoring junction, Vesicle-mediated transport, Microtubule associated complex, and etc. (Fig. 6E); while those with higher levels of methylation are enriched for Metal ion binding, FoxO signaling, Insulin signaling, G-protein, Positive regulation of insulin secretion, Cellular senescence, and etc. (Fig. 6F).These ndings suggest that differential methylation at a subset of putative enhancers distinguish b-cell subtypes at P60.We therefore examined whether any DMRs are maintained between b-cell subtypes from P2 to P60 to account for their stability.
A signi cant portion of genes that are associated with P2 DMRs between b-cell subtypes also associate with those at P60 Hierarchical analysis did not show any cluster that are differentially methylated at both stages (Fig. 6G).Only 102 DMRs remain unchanged (Fig. 6H), representing no enrichment (p=0.22,hypergeometric test).However, about one third (1,648) of the DMR-associated genes overlaps between P2 and P60 (Fig. 6I.Table S6), a 2.2-folds enrichment (p=1.11E-295,hypergeometric test).In other words, putative enhancers of 1,648 genes, e.g., DNMT3a, are differentially methylated between the b-cell subtypes at both P2 and P60, although the exact locations of the DMRs have shifted (Fig. 6J).In addition, 193 of the 1,648 common DMR-associated genes display differential expression at P60, a signi cant overlap (p=0.0001,Fisher exact test.Table S6).They are enriched for GSIS-related processes such as Heterotrimeric Gprotein complex, Positive regulation of insulin secretion, Calcium signaling, Circadian entrainment, Actin cytoskeleton organization, andetc.(Fig.6K).Thus, our combined ndings so far, together with the known roles of DNA methylation in b cells, suggest that differential enhancer methylation in early progenitor subsets can propagate after birth, producing stable b-cell subtypes at adult stages.This conclusion sets the stage to test if any factor modulates the risk of diabetes by changing the proportions of b-cell subtypes.

Maternal overnutrition leads to reduced DNMT3a expression in and reduced proportion of M + N + progenitor cells
We tested whether maternal overnutrition can alter the proportions of M + N + and M -N + islet progenitor cells.Maternal nutrition is a prominent risk factor for diabetes in offspring (Friedman, 2018), supporting a theory known as Developmental Origin of Health and Disease (DOHaD) in human diabetes.This theory posits that maternal factors program b cells to certain levels of activities that maybe un t for postnatal GSIS when nutrient levels change (Wadhwa et al., 2009).Yet the exact cellular and molecular mechanism remains unclear (Elsakr et al., 2019).
We compared the gene expression in E15.5 Neurog3 + islet progenitors that are maternally exposed to control and high-fat-diet (CD and HFD) (Fig. S6).HFD treatment results in downregulation of 2,155 genes and upregulation of 1,414 in Neurog3 + cells (Table S7).The downregulated genes are enriched for processes including Ubl conjugation pathways, microtubules, Ca 2+ transport, ER, Golgi, chromatin regulators, and Ca 2+ -regulated exocytosis (Fig. S6B), similar to that of P2 tdT + cell-enriched terms (Fig. 3C).The upregulated genes are enriched for protein translation, mitochondrion, unfolded protein response, and several biosynthetic/metabolism pathways (Fig. S6C).Intriguingly, maternal HFD reduces the expression levels of Myt1, several channels, Ca 2+ sensors, and DNMT3a in Neurog3 + progenitors (Fig. 7A), genes that are enriched in P2 and P60 b-cell subtypes with better function (Fig. 2E).Corresponding to these results, there is a trend of reduced M + N + progenitor-derived b-cell subtype in perinatal MNYI mice exposed to maternal HFD (Fig. 7B, C).This weak phenotype likely re ects the reduced b-cell subtype with better function (corresponding to the reduced Myt1 expression in Neurog3 + cells), which is mitigated by compensatory proliferation of these cells in response to higher nutrient levels (Mosser et al., 2015).Unfortunately, because b cells re-activate Neurog3 under metabolic stress during diabetes development (Talchai et al., 2012), which we expect to label M -N + progenitor-derived b cells after birth, we could not compare the function of the adult b-cell subtypes from M -N + and M -N + progenitors in the maternal HFD model.Therefore, we evaluated if the equivalent b-cell subtypes exist in human islets and if their proportion changes under diabetic conditions.
DEGs from mouse b cell subtypes subdivide human b cells into distinguishable populations.
Currently, there is a lack of consensus markers to reproducibly subdivide human b cells across all studies, likely due to variable genetic/epigenetic predisposition, asynchronous maturation, proliferation, or physiological stimuli (Mawla and Huising, 2019).We reasoned that a combination of multiple markers can overcome these sources of variation and increase the reliability with which we can identify human bcell subpopulations.We therefore selected 20 mouse DEGs as a gene signature for this goal (Fig. 7D), based on their known function and on their relatively high levels of expression in human b cells (Table S4) (Blodgett et al., 2015).
The gene signature were applied to two available single-cell transcriptome datasets, one includes all control b cells from the Human Pancreas Analysis Program (HPAP) (Kaestner et al., 2019) and the second from a cohort of 12 non-diabetic patients (Xin et al., 2018).Our gene signature consistently classify both b-cell datasets into visually notable subpopulations -population 1 (Pop.1) and Pop. 2 (Fig. 7E, F).Gini index measurement veri ed this observation, showing that the signature genes were expressed in a clustered manner (Fig. 7G) (O'Hagan et al., 2018).Pre-ranked GSEA revealed strikingly consistent pathways that are enriched in Pop. 1 of both datasets as those activated in the adult mouse tdT + b cell subtype, which include protein/peptide translation, tra cking, and mitochondria-associated pathways (Fig. 7H, I).
To explore if the proportions of Pop. 1 and Pop. 2 change in dysfunctional T2D human islets, we examined their relative representation in the T2D data set of (Kaestner et al., 2019).Our gene signature was highly effective in classifying disease state b cells into Pop. 1 and Pop. 2 b cells as in control samples using Gini index analysis (Fig. 7J).More importantly, we found a decline in Pop. 1 b cells (corresponding to the better function mouse b-cell subtype) in the T2D samples as compared to the controls (Fig. 7K).Given the known reduction of functional secretion and insulin-producing b cells in T2D, a decrease in this population associated with these same features is expected.Pre-ranked GSEA further corroborates this with Pop. 1 b cells in the T2D population having functional changes, with the upregulation of catabolic processes, a noted side-effect of T2D disease (Zanetti et al., 2020) (Fig. 7L).

Discussion
We show here that embryonic islet progenitors with differential gene expression and DNA methylation give rise to b-cell subtypes of different function in postnatal mice.Maternal overnutrition, a diabetes risk factor, alters the proportion of these progenitors.Human donor islets contain equivalent b-cell subsets, with the one of higher functionality reduced in diabetic donors.Thus, deregulating the proportions of bcell subtypes during early embryogenesis could predispose the risk of diabetes.The implication is that manipulating DNA methylation in islet progenitors, e.g. with food supplements during gestation (Li et al., 2014), may reduce the risk of diabetes.In addition, controlling nutrient levels in culture could aid the production of functional b-like cells for transplantation-based diabetes therapy.Broadly speaking, we suggest that progenitor heterogeneity may underly the presence of cell subtypes in many other organs, which can be explored with our re ned lineage tracing.
By far, most studies on b-cell heterogeneity rely on sampling single b cells at particular stages.These led to a popular conclusion that most observed heterogeneity results from unstable cell states (Bonner-Weir Our combinatorial lineage tracing enables longitudinal study of roughly half of all b cells that are derived from an islet progenitor subset.We clearly demonstrated the stability of functionally/genetically/epigenetically distinct b-cell subtypes.We further veri ed several cell markers that were previously shown to divide rodent and human b cells into functionally distinct subtypes (Dorrell et al., 2016;Dror et al., 2023), including cell surface markers CD9 and CD24 that can be used for live cell subtype puri cation.We also identi ed several DEGs (Myt TFs, Syt genes, and Ca 2+ channels) that may contribute to the functional differences of these cell subtype, underscoring the multifactorial nature of bcell heterogeneity.
For mechanistic understanding, we suggest that differential DNA methylation in islet progenitors is an important regulator of b-cell subtype production and risk of diabetes.First, inhibiting DNMT activity reduces the portion of the M + N + progenitors in embryonic pancreas.Second, DNMT3a is detected at higher levels in the b-cell descendants of M + N + cells (with higher levels of GSIS) compared with that of M - N + cells.Third, the thousands DMRs between b-cell subtypes are enriched for motifs recognized by TFs of the GATA, FoxA, FoxM1, and the HNF TF families, which facilitate b-cell production, function, and proliferation.In addition, these DMRs preferentially associate with DEGs between the two b-cell subtypes, implying their functiaonla signi cance.Fourth, the Neurog3 + progenitors that are exposed to maternal overnutrition, an established risk factor for diabetes, express lower levels of DNMT3a and Myt1.In other word, maternal overnutrition reduces the proportion of M + N + progenitors, likely resulting in a reduction of better-function b-cell subtype and a higher risk of diabetes in offspring.Consistent with this latter possibility, we show that donor islets from diabetic patients have lower proportion of the b-cell subset that is equivalent to the M + N + progenitor-derived cells.
Our studies left several pressing questions unanswered.First, it is unclear if each of our identi ed b-cell subtypes contain sub-subtypes.For example, we show that the tdT + b cells (with higher secretory function) express higher levels of Rbp4.But its expression was reported to anti-correlate with exocytotic activities (Camunas-Soler et al., 2020).In addition, we did not detect differential expression of CD63, which marks b cells of higher GSIS (Rubio-Navarro et al., 2023).Neither did we detect higher levels of mKi67 transcripts in the tdT + b cells.One explanation for these discrepancies is that the tdT + and tdT - b cells can be further subdivided, leading to varied expression levels in different cell sub-subtypes.
Second, why M + N + and M -N + b-cell progenitors have different DNMT3a/DNA methylation and how postnatal b cells regulate their methylome are not known.Future studies examining how methylating (DNMTs) and demethylating (TETs) enzymes are targeted to speci c loci will be informative.Lastly, our combinatorial lineage marking is not inducible.Reactivation of Neurog3, which happens when b cells are under stress (Talchai et al., 2012), will label M -N + progenitor-derived b cells.Thus, although Myt1 coexpression with Neurog3 during embryogenesis represents the progenitors of better-function cells, this model is not t for studying b-cell subtype stability in diabetes-prone models at postnatal stages.

Materials and Methods
Animal models and use Mouse production and usage followed protocol approved by the Vanderbilt IACUC for Dr. Gu, in compliance with the policies of AALAC.Wild type CD1 mice were purchased from Charles River.Routine crosses were used derive mice of desired genotypes, determined by gene speci c PCR reactions.

Islet β cell preparation
Islets were hand-picked after the pancreata were digested with 0.5 mg/ml type IV collagenase (dissolved in HBSS with Ca 2+ /Mg 2+ ) as previously described (Huang et al., 2018).Brie y, neonatal pancreata (younger than 1 week old) were dissected and cut into 4-6 pieces and digested with collagenase Type IV for ~ 10 minutes at 37 °C.The pancreata were broken into pieces with pipetting and washed 4 times with complete RPMI1640 media (with 11 mM glucose, 1X Penicilin/Streptomycin, and 10% FBS).Islets were then picked under a dissecting scope.Islets isolation from older mice used similar method except pancreata were perfused with the collagenase via pancreatic duct.
For single cell-preparation (used for scRNA-seq and FACS), islets were washed 3-times (3-5 minutes each) with Ca 2+ /Mg 2+ -free HBSS.They were then dissociated to single cells (at 37 °C, with 0.25% trypsin, usually 4-6 minutes), washed 2X with complete RPMI1066 media, and resuspended for encapsulationbased scRNA-seq or for FACS-based cell isolation using BD FACSAriaTMIII.Note that for FACS, DAPI or propidium iodide was used for counter-staining to exclude dead cells.Sorted live cells were immediately used for down-stream experiments.Also note that when making MNA pseudo-islets, the collected cell population was gated to minimize including non-endocrine cells in islets, which are usually smaller than endocrine cells.

Islet β cell proliferation, pseudo-islet production, and secretion assays
To assay for β-cell proliferation, islets from each MNA animal were processed as one sample, cyto-spun onto 2-4 slides for staining.After Immuno urescence (IF) staining for Ki67 and insulin (usually with 2 slides), confocal images were taken to count the portion of Ki67 + β-cell subtypes.A Leica TCS-SP5 scanning or an Olympus FV-1000 confocal microscope were used for image capture.Pseudo-islets preparation followed published protocols, with some modi cation (Friedlander et al., 2021).
In this case, sorted cells were resuspended in enriched complete RPMI media (15% FBS) at 30,000 cells per ml.Then 30 microliters were dropped on the cover of 10 cm plates, which was inverted to cover plates with ~ 5 ml 1X HBSS.The plates were left at 37 °C in an incubator for 4 days un-disturbed.Pseudo-islets were then collected, left in fresh complete RPMI 1066 media for another day for recovery.Afterwards, pseudo-islets were used for GSIS assays followed established protocol for while islet [see below and (Hu et al., 2020)].For viability assays, they were cultured with different levels of glucose and/or palmitate (see Fig. 1-related text).They were then assayed for viability via DAPI or trypan-blue intact or size.

Mitochondrial activity assays
For comparing mitochondrial mass and activity, islets from 2-3 month-old MNYI mice (both males and females) were isolated and dissociated into single cells.They were then stained with MitoView ™ 633 or MitoView ™ 650 (Biotium) at 37 °C for 15 minutes in the presence of 20 mM glucose, followed by Flowcytometry analysis (Ernst et al., 2023).These dyes assay the mitochondrial transmembrane potential and mitochondrial mass, respectively.The β-cells were gated into eYFP + and eYFP -β-cell subtypes in order to compare their MitoView ™ signals.For ADP/ATP ratio in puri ed β-cell subtypes, an EnzyLight ™ ADP/ATP ratio Assay Kit (BioAssay Syetems) was used, following the manufacturer's protocols.

Insulin secretion assays
For insulin secretion using pseudo-islets or intact islets, the % of total insulin secreted within a 45-minute window was measured unless noted.(Pseudo-)islets were allowed to recover in RPMI-FBS for 2 hours or overnight.Islets were washed twice with pre-warmed KRB solution (2.8 mM glucose, 102 mM NaCl, 5 mM KCl, 1.2 mM MgCl 2 , 2.7 mM CaCl 2 , 20 mM HEPES, 5 mM NaHCO 3 , and 10 mg/ml BSA, pH 7.4) and then incubated in KRB (37 °C) for one hour, washed with pre-warmed KRB once more.10-15 (pseudo-)islets were then transferred into each of the wells of 12-well plates with 0.2 (pseudo-iselts) 1 ml (islets) prewarmed KRB to start the secretion assays.The basal solution contains 2.8mM glucose in KRB.The stimulatory conditions use 20 mM glucose or (20 mM glucose + 30 mM KCl).For all assays, four or more mice of each genotype were used for islet isolation, with islets from two or more mice mixed and examined as 2-3 technical replicas.If Pseudo-islets were used, a starting volume of 200 microlitter with 20-30 pseudo-iselts was used due to their small size.Total insulin was assayed after ethanol-acid extraction as in (Huang et al., 2018).Insulin was measured with an Elisa kit from ALPCO.Assays from human islets use similar method, except the basal glucose was 3.3 mM.

Intraperitoneal glucose tolerance test (IPGTT) in mice and IF assays in tissue sections
IPGTT followed previously described method (Hu et al., 2020).Brie y, mice were fasted overnight.Glucose was injected at 2mg/kg.Blood glucose was then measured via tail nip.IF assays follow routine method using both frozen and para n sections, using antibodies listed in the reagent table.

Ca 2+ recording with GCaM6
Ca 2+ in ux in individual β cells were monitored using GCaMP6 as a surrogate using a similar approach as in (Dickerson et al., 2019).Brie y, islets were isolated from 2-months old MNA mice.They were allowed to attach to bronectin-coated plates with glass-bottoms overnight in RPMI1640 media.The islets were then infected with adenoviruses that express GCaMP6 speci cally in β cells.Three days later, the islets were incubated (20 minutes, 37°C) in KRB supplemented with 1 mM glucose for acclimation.
Real-time uorescence imaging then followed to record was then performed using a Nikon Eclipse TE2000-U microscope equipped with an epi uorescence illuminator (Sutter, Inc), a CCD camera (HQ2; Photometrics Inc), and Nikon Elements software.Islets were perifused at 37°C at a ow of 2 mL/min with appropriate KRB-based solutions that contained glucose concentrations and compounds speci ed in the gures.The ratios of emitted uorescence intensities at excitation wavelengths of 340 and 380 nm (F340/F380) were determined every 5 seconds with aid of agarose.For insulin secretion, the % of total insulin secreted within a 45-minute window was measured unless noted.(Pseudo-)islets were allowed to recover in RPMI-FBS for 2 hours or overnight.Islets were washed twice with pre-warmed KRB solution (2.8 mM glucose, 102 mM NaCl, 5 mM KCl, 1.2 mM MgCl 2 , 2.7 mM CaCl 2 , 20 mM HEPES, 5 mM NaHCO 3 , and 10 mg/ml BSA, pH 7.4) and then incubated in KRB (37 °C) for one hour, washed with pre-warmed KRB once more.10-15 islets were then transferred into each of the wells of 12-well plates with 1 ml prewarmed KRB to start the secretion assays.For all assays, four or more mice of each genotype were used for islet isolation, with islets from two or more mice mixed and examined as 2-3 technical replicas.Total insulin was assayed after ethanol-acid extraction as in (Huang et al., 2018).Insulin was measured with an Elisa kit from ALPCO.Assays from human islets use similar method, except the basal glucose was 3.3 mM.If Pseudo-islets were used, a starting volume of 200 microlitter with 20-30 psuedo-iselts was used due to their small size.
Then a short DNA fragment that drives Syt7 guide RNA (5'-GAACATGTGGGGCTCCAGCG-3') under the control of a human U6 promoter was inserted upstream of the CMV enhancer.The resulting plasmid was electroporated into MIN6 cells, which was characterized 5 days after culture.For mRNA assays, cells from three independent plates were FACS-ed based on eGFP expression and studied for mRNA expression using routine reverse transcription coupled with real-time-PCR (see below).

Real-time RT-PCR-based gene expression assays
Total RNA samples were extracted from isolated islets (from at least three animals) of puri ed cells (3 batches from independent FACS collection).Pools of cDNA were then prepared using a High Capacity cDNA Reverse Transcription Kit (Thermal Fisher), followed by CYBR-green-based RT-PCR assays (BioRad).The input for each sample was normalized using GAPDH levels.independent plates were sorted based on eGFP expression and studied for mRNA expression using routine reverse transcription coupled with realtime RT-PCR, using GAPDH as a normalizer (oligos used: 5'-AACTTTGGCATTGTGGAAGG-3' and 5'-GGATGCAGGGATGATGTTCT-3').Oligos for Syt7 are (5'-CTTTCGCCTTCGACATACCC-3' and 5'-GGCTGAGCTTGTCTTTGTCC-3').
In vitro pancreas bud culture assay As shown in (Liu et al., 2019), E14.5 pancreas buds were dissected and cut into 3-4 pieces.They were then cultured in RPMI1066 supplemented with 10% FBS on lters in an air-liquid interface format, with DMSO or 2 µM AzaC.About 48 hours later, the tissue were prepared for frozen section and used for immuno uorescence assays using antibodies listed in the reagent table.Immunostaining follows standard protocols as in (Liu et al., 2019).

Maternal diet treatment for bulk RNAseq and labeled cell counting
Six-seven weeks-old female CD1 mice were fed with CD or HFD (Bio-ServTM F3282, with 60% calories from fat compared to 6% in CD) for 2 weeks.They were then crossed with Ngn3 eGFP or MNYI mice under CD or HFD feeding until embryonic collection or birth.
For isolating islet progenitors, pregnant females were dissected at E15.5.Pancreata with eGFP expression were pooled, dissociated, and used for cell isolation via FACS.Embryos from 2-3 litters of mice were pooled and processed as one sample for RNA preparation and sequencing as previously published (Huang et al., 2018).For counting β cells that were labeled with eYFP expression in MNYI mice, newly born mice and DAM were switched to CD feeding.Pups were then euthanized at P4 for pancreas collection.Islet from those with eYFP and Apple will be used for islet isolation, which were imaged using LSM for cell counting.At least three litters of mice were used for counting.

Single cell RNA-sequencing and analysis
Murine scRNA-seq, single-cell encapsulation and library generation Dissociated cells from islets were washed and checked for viability with Trypan Blue.Samples with > 95% viability were counted and mixed with 10% K562 cells to evaluate doublet formation rate.Cells were washed using 0.02%BSA in DPBS before resuspension in a solution of DPBS with 15% Optiprep at 100,000 cells/mL before proceeding to inDrops (1CellBio) encapsulation as described in (Klein et al., 2015;Liu et al., 2019;Simmons and Lau, 2022).In brief, cells were owed into a micro uidic chip where cells were mixed with an RT/Lysis buffer and barcoded capture beads immediately before being portioned into droplets.Cell concentrations and ow rates were carefully controlled to minimize doublets and ensure droplet uniformity.The droplet emulsion was collected into a 1.5 ml tube in fractions of ~ 3,500 estimated captured cells.Droplets then proceeded to UV exposure using a UV Cleaver to release captured oligos from the beads, a 1-hour incubation at 50 o C to facilitate reverse transcription, and nally a 5-minute incubation at 85 o C to stop the reaction.Emulsions were then cooled on ice and demulsi ed to isolate the aqueous component containing barcoded cellular cDNA, which is stored at -80 o C for library preparation.Libraries were prepared for sequencing utilizing the TruDrop protocol (Southard-Smith et al., 2020).Samples were then sequenced using the Novaseq 6000 (Illumina) in a customized format.For each sample both tdT + and tdT − β cells were processed and sequenced together and contained islets from 3-4 mice.During sequencing we targeted 120 million reads.
scRNA-seq, alignment, droplet matrix generation, and droplet matrix quality control Sequences and droplets were aligned and ltered as described in (Osumi-Sutherland et al., 2021).Brie y, the libraries were processed utilizing the DropEst pipeline (Petukhov et al., 2018) (Fig. S2A, B) and the STAR aligner with Ensemble reference genome, CRCh38 release 25 (Dobin et al., 2013).High-quality single-cell-containing droplets were identi ed employing the dropkick algorithm (Heiser et al., 2021) in conjunction with prior-knowledge in gene expression pro ling (Fig. S2B).Normalization, inverse hyperbolic sine transformation, and Z-score scaling was performed on cells as an AnnData object utilizing the Scanpy toolkit (Wolf et al., 2018).These normalized matrices were then projected into two dimensions using a UMAP initialization based on their 50 principal component decomposition.Geneexpression consistency was con rmed and the nal selection of cell-containing droplets was made based on quality metrics such as total counts, mitochondrial count percentages, and gene-expression pro les.

scRNA-seq, count matrix normalization and UMAP visualization
The Scanpy toolkit and numpy functions were used on the raw count data to normalize by median library size, perform log-like transformation with Arcsinh, and standardize each gene using Z-score (Chen et al., 2021;Wolf et al., 2018).Two models of UMAP visualization were employed after normalization: one based on the standard 50 principal component decomposition (as described in scRNA-seq, alignment, droplet matrix generation, and droplet matrix quality control section), and another using Harmonycorrected principal components to correct any sample speci c batch effects.UMAP visualizations for murine data were generated using the 'scanpy.tl.umap' function, while human UMAP visualizations utilized this function or their own embedded x_UMAP pro le.Sample numbers were visualized utilizing the same x_UMAP pro le to ensure even sample distribution before further analysis (donor_ID Fig. S7).Both human and murine inputs utilized their 50 principal component decompositions.Murine data were integrated with the Harmony algorithm (Korsunsky et al., 2019) which adjusts these principal components.
scRNA-seq, cell type labeling and unsupervised clustering Utilizing the scMRMA (Li et al., 2022) algorithm in R, the PanglaoDB was utilized to label major pancreatic cell types.These algorithm-de ned labels were then cross-checked with prior-knowledge geneexpression pro les and corrected to label major pancreatic cell types including: α cell, β cells, δ cells, endothelial cells, and enteric neurons.Employing the Scanpy toolkit, the Leiden algorithm single-cell subpopulations were labeled as described in Chen et al (Fig. 2A) (Chen et al., 2021).
Murine scRNA-seq, pseudobulk analysis Cells were subset based on the 'β cell' identity (as de ned by the scRNA-seq, cell type labeling and unsupervised clustering section).Utilizing the Scanpy toolkit, cells were normalized, Arcsinh transformed, then categorized into tdTomato-positive (tdT + ) and tdTomato-negative (tdT − ) subsets based on histogram plots of expression with the cutoff for positivity being 0.9.tdT + and tdT − cells were separately averaged by sample before being averaged by sex or in combination then the difference between tdT + and tdT − was determined and numerically sorted, where tdT was the most highly expressed gene across all β cell samples.To identity differentially expressed genes, low-expressing genes were rst ltered out (arti cially set at Arcsinh-transformed reading below 0.002 per sample, roughly one read in each cell).Paired t-test was then used to identify differentially expressed genes [p < 0.05, with a difference of at least 5% (log-transformed) between the two samples].

Bioinformatic gene set enrichment analysis
After identi cation of DEGs, the Database for Annotation, Visualization and Integrated Discovery (DAVID) was primarily utilized for comprehensive functional clustering analysis (Sherman et al., 2022), which generates clusters of gene ontogeny/processes/molecular signature for functional interpretation.In all cases, similar/overlapping pathways or processes were consolidated and summarized for presentation.
Only top 15 or those having [-Log10(p value)] > 1.5 (if less than 15 processes) were presented.In addition, we did not include several general terms such as "transcription" and "cancer" in our lists, because their functional implication in β cells is not informative.In addition, we also utilized the pre-ranked DEGs in GO term enrichment and GSEA utilizing the GSEApy package (Fang et al., 2023).Pathway analyses included GO_Biological_Process_2023, WikiPathways_2019_Mouse, KEGG_2019_Mouse, Reactome_2022, and MSIGDB_Hallmark_2020.
Human scRNA-seq, processing and gene signature scoring Two publicly available datasets with focuses on adult T2D disease and control β cell populations were selected for analysis (Fig. S7).The data were normalized and transformed as described using leidon and donor ID-based mapping (in the section on scRNA-seq, cell type labeling and unsupervised clustering) (Fig. S7A, B, D, E, G, H).We selected a set of 20 functionally known genes upregulated in the murine pseudobulk analysis.These 20 genes were utilized to create a gene signature score employing the 'scanpy.tl.score_genes' function as described in (Osumi-Sutherland et al., 2021; Wolf et al., 2018).Each signature was calculated and standardized using a reference sample size of 2000 genes.Scores were subset based on positivity, with 'population 1' (Pop. 1) β cells containing a score higher than 0.2 and 'Pop.2' β cells containing a score lower than 0.2 (gene score human Fig. S7C, F, I).To determine if the 'Pop.1' β cell genes were signi cantly clustered together, rather than randomly, we utilized the package GiniClust3 on the pre-processed data with our selected gene list (Dong and Yuan, 2020).This identi ed Gini index values and a corresponding p-value, where a p-value < 0.05 indicates the expression of our gene signature is signi cantly clustered rather than randomly.
Human scRNA-seq, differential gene-expression and gene set enrichment analysis Utilizing the Scanpy toolkit, differential testing of gene expression was performed based on population labels (as de ned by the Human scRNA-seq, processing and gene signature scoring section).The 'scanpy.tl.rank_genes_groups' function was utilized using Wilcoxon rank-sum and only those statistically signi cant were placed through GSEA.These pre-ranked genes were placed through GSEA using the GSEApy function (as described in the scRNA-seq, gene set enrichment analysis and differential geneexpression testing section).Pathways analyzed included GO_Biological_Process_2023, WikiPathways_2019_Human, KEGG_2021_Human, Reactome_2022, and MSigDB_Hallmark_2020.

scRNA-seq, population quanti cation and gene-expression analysis
Percent change in human 'population 1' β cells was quanti ed per patient in both control and T2D β cells and plotted using the Matplotlib and Seaborn toolkits.Utilizing the Scipy toolkit, signi cance was tested using the 'scipy.stats.ttest_ind'assuming equal variance.To quantify gene-expression between tdTomato populations, a gene signature score for each gene of interest was calculated utilizing the methods described above (in Human scRNA-seq, processing and gene signature scoring).In order to analyze the scores, a scatterplot was utilized to visualize each score and statistical tests were conducted on their distributions.An initial Kruskal-Wallis test was performed and if the null hypothesis was rejected, a post hoc Mann-Whitney U test was conducted.P-values of less than 0.05 were considered signi cant.
DNA methylome assays and analysis DNA methylome assays largely followed published protocols for Tagmentation-based whole genome bisul te sequencing (T-WGBS) (Adey and Shendure, 2012; Barnett et al., 2020;Wang et al., 2013), which relies on a hyperactive Tn5 transposase to fragment DNA while adding adaptor oligo sequences for sequencing (Wang et al., 2013).The genomic DNA prepared from FACS puri ed β-cell subtypes were used.The key steps for these assays are:

Transposome Preparation
T-WGBS adapters were prepared by annealing the oligonucleotides in PCR tubes (10µl 100µM Tn5mC-Apt1 oligo, 10µl 100µM Tn5mC1.1-A1blockoligo, 80µl nuclease free water).Oligos were incubated in a PCR thermocycler as follows: 95°C for 3 min, 65°C for 3 min, ramp down to 24°C at a rate of − 1°C/second, hold at 24°C.Annealed oligos were then combined with 100 µl glycerol to create a 5µM, 50% glycerol adapter mixture.Transposomes were assembled by combining equal parts puri ed Tn5 transposase enzyme and annealed oligos, followed by a 25°C incubation for 60 min at room temperature.

T-WGBS
T-WGBS libraries were prepared as previously reported (Adey and Shendure 2012, Barnentt et al. 2020, Wang et al., 2013).Genomic DNA from puri ed β-cell subtypes were extracted using a kit from Zymo (D3020).For each stage/sex, at least three batches of cells from FACS (4-7 mice for each batch, at least 30,000 cells) were combined, yielding 100-150 ng DNA.Genomic DNAwas diluted in a 50µl tagmentation reaction (40 ng genomic DNA, 2.5µl transposome assembled with T-WGBS adapters, 10µl 5X Tris-DMF, XuL Nuclease-free (NF) water to obtain a total reaction volume of 50µl) and incubated at 55°C for 8 min in a PCR thermocycler.Tagmentation reactions were then stopped by adding 250µl Zymo DNA binding buffer from the DNA Clean & Concentrator-5 kit (Zymo).Reactions were puri ed per the manufacturer's protocol instructions in the DNA Clean & Concentrator-5 kit (Zymo) and eluted in 15µl NF-water.DNA eluate from the tagmentation reaction was used as input into the subsequent oligo replacement and gap repair reaction (11µl DNA eluate, 2µl 10µM Tn5mC-Repl01 oligo, 2µl 10x ampligase buffer, 2µl dNTPs 2.5mM each).This gap repair reaction was assembled in a PCR tube and incubated as follows in a PCR thermocycler: 50°C for 1 minute, 45°C for 10 minutes, ramp down to 37°C at a rate of − 0.1°C/second, hold at 37°C.Upon reaching 37°C, 1µl T4 DNA polymerase and 2.5µl ampligase were both added separately without removing the tube from the thermocycler.The reaction was mixed using a pipette without removal of the tube from thermocycler, minimizing bubbles.The gap repair reaction was subsequently incubated as follows: 37°C for 30 minutes, hold at 4°C.Optionally, 2µl of the gap repair reaction can be reserved for a test PCR ampli cation to determine library distribution pre-bisul te conversion.The gap repair reaction was stopped with the addition of 102.5 µl of Zymo DNA binding buffer from the DNA Clean & Concentrator-5 kit (Zymo).Reactions were puri ed per the manufacturer's protocol instructions in the DNA Clean & Concentrator-5 kit (Zymo) and eluted in 20.5µl NF-water.Gap repaired, tagmented DNA was subsequently bisul te converted according to the manufacturer's protocol instructions using the EZ DNA Methylation-Gold kit (Zymo).For the bisul te conversion reaction, 20.5µl gap repaired DNA was mixed with 130µl CT conversion reagent in a PCR tube and incubated in a PCR thermocycler as follows: 98°C for 10 min, 64°C for 2.5 hrs, hold at 4°C.Final puri cation/desulfonation was performed as directed by the kit manufacturer's protocol instructions.Final elution was in 25µl NFwater.Eluted bisul te-converted DNA was ampli ed and barcoded in 50µl PCR reactions (25µl 2x KAPA HiFi HotStart Uracil + ReadyMix, 20µl eluted bisul te-converted DNA, 1.5µl 10µM i5 index primer, 1.5µl 10µM i7 index primer) as follows: 98°C 45 sec; 10 cycles of 98°C 15 sec, 62°C 30 sec, 72°C 30 sec; nal extension 72°C 2 min; hold at 12°C.Post-ampli cation PCR reactions were puri ed using the DNA Clean and Concentrator-5 kit (Zymo) and eluted in 22µl NF-water.Library concentrations and size distributions were evaluated using an Agilent 2200 TapeStation with a D5000 screentape.T-WGBS DNA libraries were sequenced at Vanderbilt University's genomics research core using 2×150bp paired-end reads on the Illumina NovSeq 6000, obtaining approximately 300-500 million reads per library.(D-H) Insulin secretion assays in pseudo-islets from b-cell subsets.tdT + and tdT -islet cells were separated by FACS (D), made into pseudo-islets (E), and assayed for insulin secretion in response to G2.8 (2.8 mM glucose), G20, or (G20 + 30 mM KCl).The % of total insulin secreted during a 45 minutes window was presented.(F, G) Secretion from pseudo-islets of 2-months old MNA male or female mice, respectively.(H) Pseudo-islets secretion from 8-months old MNA mice, with males and females mixed.Presented are (mean + SEM).p values are calculated with unpaired t-test.
The P2 DMRs are enriched for binding motifs of several TFs belonging to the FOXA, GATA, and HNF families (Fig.5F), all of which have established roles in b-cell production and/or function(Golson et  al., 2015; Gu et al., 2010; Gupta et al., 2007; Lorberbaum and Sussel, 2017; Reizel et al., 2021; Weng et al., 2023 and Aguayo-Mazzucato, 2016;Xin et al., 2018;Zeng et al., 2017).Some lineage-marking studies support this conclusion.For example, Bader et al. detected two b-cell subsets that express different levels of Fltp/CFAP126 in mice, but one of them convert to the other over time(Bader et al., 2016).In contrast, another lineage study showed the presence of virgin b cells that stay immature for nearly half the lifespan in mice(Lee et al., 2021).Along a similar line, two recent studies, using DNA/histone methylation as lineage-surrogates, concluded that stable b-cell subtypes exist(Dror et  al., 2023; Parveen et al., 2023).However, the virgin b cells represent only a few percent of all b cells, while DNA/histone methylation patterns change according to physiological conditions [our results here and (Avrahami and Kaestner, 2019; Avrahami et al., 2015; Hall et al., 2018; Hall et al., 2019; Wortham et al., 2023)].Thus, the stability of b-cell subtypes needs direct substantiation, a prerequisite for determining the stage when subtypes arise and whether/how deregulating proportions of cell subtypes causes diabetes.

Figures
Figures

(
I-O) b-cell death assays in pseudo-islets.(I) Examples of newly-formed pseudo-islets (5 days in hanging drops) from eYFP + and eYFP -b cells of MNYI mice.(J) Relative pseudo-islet size (represented by pseudoislet area measured from microscopic images).(K) Images of typical live pseudo-islet (white arrow points at one example) 7-days after culture at 11 mM glucose, stained with trypan blue to visualize dead cells (red arrows).The numbers (L) and areas (M) of the surviving pseudo-islets (white arrow) were used to compare the viability of two types of pseudo-islets.Three batches of pseudo-islets were used for this assay.(M) The proportions of islets (using area as a parameter) that survived.(N, O) Pseudo-islet cell death induced by 24-hours treatment by 20 mM glucose and 0.6 mM palmitate in eYFP + and eYFP - pseudo-islets.Presented in (M) are (mean + SEM), with p values calculated with unpaired t-test.Scale bars, 20 mm.For all assays, at least two batches of pseudo-islets were used, derived from different islet isolation, dissociation, cell sorting, and reaggregation plates.

Figure 4 Several
Figure 4

Figure 5 M
Figure 5

Figure 6 DNA
Figure 6