Beta cell-derived cholecystokinin drives obesity-associated pancreatic adenocarcinoma development

Domain: chemistry
Published: 2026-02-27
Editorial outcome: major_revision

Abstract
Pancreatic endocrine-exocrine crosstalk plays a key role in normal physiology and disease and can be altered by host metabolic states, such as obesity. Classically, endocrine islet beta (β) cell secretion of insulin is thought to promote the development of obesity-associated pancreatic adenocarcinoma (PDAC), an exocrine cell-derived tumor. Here, we show that β cell expression of the peptide hormone cholecystokinin (CCK) is necessary and sufficient for obesity-associated PDAC progression in mice and that CCK expression – rather than insulin – correlates strongly with enhanced tumorigenesis. Single-cell RNA-sequencing, in silico latent-space archetypal and trajectory analysis, and experimental lineage tracing in vivo reveal that obesity induces the expansion of postnatal immature β cells, which adapt to express CCK via stress-responsive JNK/cJun signaling. Finally, obesity perturbs CCK-dependent peri-islet exocrine cell transcriptional states and enhances islet-proximal tumor formation. These results define endocrine-exocrine CCK signaling as a bona fide driver of obesity-associated PDAC development and uncover avenues to target the endocrine pancreas to subvert exocrine tumorigenesis.

Introduction
The pancreas is comprised of two functionally distinct cellular compartments: 1) the endocrine pancreas, composed of islet cells that secrete hormones to maintain glucose homeostasis, and 2) the exocrine pancreas, encompassing acinar cells that produce digestive enzymes and the ducts through which these enzymes traverse. Recent studies have challenged the longstanding dogma that the endocrine and exocrine compartments function independently, demonstrating endocrine-exocrine crosstalk in normal physiology and disease1,2. For example, maturity-onset diabetes of the young type 8 (MODY 8) is associated with mutations of CEL, a gene that encodes for carboxyl ester lipase produced and secreted by acinar cells3. Mutant CEL is acquired by islet beta (β) cells and aggregates intracellularly to induce endoplasmic reticulum (ER) stress, resulting in β cell dysfunction and diabetes, a disease of the endocrine pancreas4. Other exocrine diseases have similarly been associated with the development of diabetes, including chronic pancreatitis, Wolcott-Rallison syndrome, cystic fibrosis, and pancreatic ductal adenocarcinoma (PDAC)2. Conversely, obesity and diabetes – host metabolic states linked to β cell dysfunction – are associated with an increased risk of developing and dying of PDAC5-10, a highly lethal tumor thought to primarily arise from acinar cells of the exocrine pancreas11,12. However, the cellular and molecular mechanisms that govern endocrine-exocrine interactions in tumorigenesis are not completely understood.

To date, few studies have directly examined the importance of β cells and β cell-secreted hormones in PDAC pathogenesis. Recent experiments have demonstrated a functional role of insulin signaling in pancreatic tumorigenesis. Partial knockout of the insulin genes (Ins1 and Ins2) or acinar cell-specific knockout of the insulin receptor (InsR) reduced oncogenic Kras-driven tumor development in mice fed a high-fed diet (HFD)13-15. Mechanistically, insulin signaling in acinar cells increased digestive enzyme production and acinar-to-ductal metaplasia (ADM)15, an early prerequisite step in PDAC development12,16. Our lab previously showed that obesity (due to loss of the appetite suppression hormone leptin (Lepob/ob)17) induces aberrant β cell expression of the peptide hormone cholecystokinin (CCK) in mice and that islet CCK expression is associated with increased body-mass index (BMI) in humans18. CCK stimulates acinar cell proliferation, digestive enzyme production, and ADM18-20, arguing that CCK – like insulin – may drive PDAC development. Indeed, exogenous administration of the CCK analogue cerulein promotes oncogenic Kras-driven tumorigenesis in mice21. We crossed Lepob/ob mice with an autochthonous mouse model that mimics the genetic and histologic features of human PDAC progression (KC: Pdx1-Cre;KrasLSL-G12D)22. Resultant obese KCO (KC;Lepob/ob) mice exhibited enhanced preinvasive pancreatic intraepithelial neoplasia (PanIN) formation and progression to invasive PDAC, which was associated with increased CCK expression18. Furthermore, both weight loss and the antidiabetic dapagliflozin reduced β cell CCK expression and neoplasia18, arguing that β cell-derived CCK is a marker of both a dysregulated metabolic state and enhanced exocrine tumorigenesis.

Here, we demonstrate that β cell CCK expression is necessary and sufficient for obesity-associated PDAC development in vivo, validating CCK as a bona fide driver of exocrine tumor formation. Leveraging single-cell RNA sequencing (scRNA-seq) of multiple congenic obesity models and a suite of machine learning-based computational tools for in silico lineage tracing (TrajectoryNet)23, archetypal analysis (AAnet)24,25, batch integration across datasets (scMMGAN)26, optimal-transport distance-preserving embedding of patient data (DiffusionEMD)27, and a framework for uncovering cellular transitions and dynamic regulatory interactions (Cflows)28,we gain deep insights into β cell heterogeneity and the dynamic transcriptional changes that lead to β cell CCK expression in obesity (Supplementary Fig. 1). Through both in silico and experimental lineage tracing approaches and gene regulatory analysis, we further discover that CCK+ β cells emerge from stress-induced expansion and adaptation of a postnatal immature β cell population and that CCK expression is mediated by stress-responsive JNK/cJun signaling. Finally, we show that obesity is associated with CCK-dependent transcriptional alterations in peri-islet acinar cells and increased peri-islet tumorigenesis. Together, our results establish the critical importance of endocrine-exocrine signaling in PDAC pathogenesis and implicate β cell stress pathways and hormones that could be targeted to intercept exocrine tumorigenesis.

Results
β cell CCK expression promotes PDAC progression
While obesity is associated with increased β cell CCK expression18,29, the extent to which β cell CCK can independently drive tumorigenesis remains unclear. To determine this, we crossed a transgenic model of β cell-specific CCK (mouse gene Cck) overexpression (Ins1-Cck)30 to our previously generated KCO model of obesity-associated PDAC (Fig. 1a). As anticipated, 3-month-old lean KC;Lepob/+;Ins1-Cck/+ and obese KCO mice exhibited increased pancreatic CCK expression relative to lean KC;Lepob/+ littermates (Fig. 1b,c). Strikingly, KC;Lepob/+;Ins1-Cck/+ mice demonstrated a significant increase in disease burden (combined ADM, PanIN, and PDAC) and propensity to progress to PDAC compared to age-matched KC;Lepob/+ littermates (Fig 1d,e). These differences in tumorigenesis were independent of changes in weight, glucose, or endogenous insulin production (C-peptide) (Fig. 1f-h), and C-peptide levels were not associated with Cck expression nor disease burden in KC;Lepob/+;Ins1-CCK/+ mice (Fig 1i,j). Critically, there was no significant difference in overall disease burden or progression to PDAC between KC;Lepob/+;Ins1-Cck/+ and KCO littermates (Fig. 1d,e), suggesting that β cell CCK overexpression is sufficient to phenocopy the pro-tumorigenic effects of obesity. Together, these experiments argue that β cell CCK is a potential insulin-independent driver of pancreatic tumor development in obesity.

Pancreatic CCK drives obesity-associated PDAC development
To ascertain if CCK is necessary for obesity-associated PDAC development, we obtained a conditional model for Cck knockout (Cckflox/flox) generated by the International Knockout Mouse Project31. We first validated pancreas-specific Cck knockout in obese non-tumor-bearing 4-month-old Pdx1-Cre;Lepob/ob;Cckflox/flox mice (Supplementary Fig. 2a,b), as we observed pancreatic Cck expression was abolished while duodenal expression was preserved (Supplementary Fig. 2c-e). Since CCK has been reported to be a β cell survival factor29,30,32, we next evaluated the role of pancreatic Cck knockout on islet health. Pdx1-Cre;Lepob/ob;Cckflox/flox mice displayed comparable glucose levels and islet mass to Cckflox/+ and Cck+/+ controls (Supplementary Fig. 2f,g), corresponding to similar C-peptide levels under both fasted and fed conditions (Supplementary Fig. 2h). These data argue that β cell CCK expression is dispensable for maintaining islet homeostasis in obesity.

We next tested whether β cell Cck knockout altered tumor development by crossing Cckflox/flox mice into the KCO model (Fig. 2a). As expected, 3-month-old KCO;Cckflox/flox littermates demonstrated significantly decreased pancreatic CCK expression compared to age-matched littermate controls (Fig. 2b,c). Remarkably, KCO;Cckflox/flox mice exhibited significantly reduced disease burden and progression to adenocarcinoma relative to KCO;Cck+/+ littermates (Fig. 2d,e) despite increased weight and C-peptide (Supplementary Fig. 2i-k). Partial reductions in disease burden and progression were also observed in KCO; Cckflox/+ mice (Fig. 2d,e), consistent with a gene dosage effect. Strikingly, we observed that pancreatic Cck expression was inversely correlated with insulin secretion (C-peptide) in multiple cohorts of KCO mice, including those subject to Cck knockout (Fig. 1k and 2f). In fact, Cck gene expression strongly correlated with disease burden, whereas C-peptide was negatively correlated (Fig. 1l,m and 2g,h). Combined, these data indicate that β cell CCK expression is necessary for obesity-associated PDAC progression and that CCK – rather than insulin – is the primary driver of pancreatic tumorigenesis in these obesity-mediated tumor models.

CCK and insulin expression are inversely correlated in obesity
We next investigated the molecular mechanisms by which obesity reprograms β cells to express CCK and drive tumor progression. We performed scRNA-seq of isolated islet cells from aged-matched male congenic (C57/B6) wild-type (WT), moderately obese HFD-fed, and severely obese Lepob/ob mice that displayed increasing weight and glucose levels (Supplementary Fig. 3a-c). Consistent with these phenotypic observations, primary islets isolated from HFD-fed and Lepob/ob mice exhibited decreased glucose-stimulated insulin secretion (GSIS), a marker of β cell dysfunction (Supplementary Fig. 3d,e). After pre-processing of scRNA-seq data (see Methods), we visualized and clustered 23,469 single cells from all samples (Supplementary Fig. 3f). We annotated these clusters based on published endocrine and exocrine markers from human and mouse islet maps33 (Supplementary Fig. 3g) and subsetted the data to 20,294 endocrine cells for downstream analysis. We re-embedded the endocrine cells and annotated clusters corresponding to alpha (ɑ), β, delta (δ), and pancreatic polypeptide (PP) cells based on high expression of glucagon (Gcg), insulin (Ins1, Ins2), somatostatin (Sst), and polypeptide (Ppy), respectively (Supplementary Fig. 4a). Similar to prior studies34,35, a small population of polyhormonal cells (expressing insulin and a second hormone) was also observed. Cck was primarily expressed in β cells, which showed the greatest variation in transcriptional states across conditions (Supplementary Fig. 4a) and were the most abundant islet cell type (Supplementary Fig. 4b).

We subsequently re-embedded β cells from all samples (Supplementary Fig. 4c,d) and continuously inferred individual and heterogeneous cellular trajectories over the progression of obesity using an ODE-based dynamic optimal transport neural network termed TrajectoryNet and the Cflows framework23,28. TrajectoryNet is designed to interpolate continuous dynamics for every single cell across distinct timepoints, where the dynamics learned are biologically plausible through ensuring energy efficiency and modeling cellular proliferation. By learning individual cellular trajectories and leveraging real timepoint dynamics to guide trajectories, this approach represents an advancement over traditional trajectory-based inference methods, which infer a single pseudotemporal ordering over the entire population of cells based on expression similarity, which may not correspond to true latent dynamics36. Here, we learned individual cellular trajectories from the combined WT samples (timepoint 1) to HFD (timepoint 2) to Lepob/ob (timepoint 3) conditions, and transcriptional changes along the learned trajectories were analyzed by decoding them back into the gene expression space. This analysis revealed a clear trajectory corresponding to increasing obesity (Supplementary Fig. 5a,b). Along this trajectory, insulin (Ins1, Ins2) showed a reciprocal relationship with CCK, wherein the insulin genes decreased while Cck increased in expression (Supplementary Fig. 5c). These findings were concordant with the inverse relationship between pancreatic Cck expression and C-peptide in tumor-bearing KCO mice (Fig. 1k and 2f). Gene trends learned by TrajectoryNet showed a decrease in β cell identity, maturation, and insulin secretion markers and increased protein processing, dedifferentiation, and ER stress markers as obesity progresses (Supplementary Fig. 5d). These data argue that β cells retain secretory capacity but upregulate Cck at the expense of insulin. We took advantage of the Min6 murine insulinoma cell line to directly detect whether CCK is sorted into granules for secretion. This frequently employed β cell model expresses insulin and CCK at comparable levels to β cells of KCO mice (Supplementary Fig. 5e). Co-immunoelectron microscopy revealed that CCK was present in the same secretory granules as insulin, including those at the plasma membrane (Supplementary Fig. 5f). These findings indicate that β cell CCK may be secreted to drive tumor progression.

Archetypal analysis reveals obesity-dependent alterations in β cell heterogeneity
Recent studies have suggested that adult β cells are heterogeneous and that host metabolic conditions can modulate this heterogeneity34,37,38. To characterize obesity-associated β cell heterogeneity, we separately re-embedded the cells for each condition (WT, HFD, and Lepob/ob) and applied the archetypal analysis approach AAnet. AAnet is an autoencoder that embeds cells into a simplicial latent representation and defines them with respect to cellular “archetypes,” or extreme states24,25, termed “latent-space” archetypal analysis. Unlike traditional archetypal analysis, AAnet does not assume that cells are archetypal in gene expression space but instead finds a lower dimensional latent space that is amenable for such factor analysis. This enables the characterization of β cell states while preserving continuous variation between highly plastic ones. AAnet identified seven archetypes (AT) within cells from WT mice, four ATs from HFD-fed mice, and three ATs from Lepob/ob mice, where the number of ATs was determined based on the elbow point of reconstruction error20 (Fig. 3a and Supplementary Fig. 6a). ATs showed varying abundance across the three conditions (Supplementary Fig. 6b). Similarity scores calculated between all ATs identified three groups with shared signatures (Groups 1, 2, and 3) and differing proportions across conditions (Supplementary Fig. 6c,d).

We characterized these β cell ATs using marker genes defined by previous research on β cell heterogeneity37,39,40 (Fig. 3b), visualizing each marker gene and its variability (Supplementary Fig. 6e-g and Supplementary Fig. 7-9) and determining differential expression across all genes per condition for each AT (Supplementary Data 1). Group 1 ATs (WT AT 5 and 7; HFD AT 1; Lepob/ob AT 1) were significantly (q < 0.05) enriched for expression of insulin, insulin secretion markers, and β cell maturation markers, consistent with a mature β cell phenotype (Fig. 3b and Supplementary Data 1). Group 3 ATs (WT AT 1, 2, 3, and 6; HFD AT 2 and 3) were characterized by significantly (q < 0.05) higher mitochondrial expression and – compared to Group 1 – reduced insulin expression (Fig. 3b and Supplementary Data 1). For WT AT 1 and HFD AT 3, maturation and insulin secretion markers Ucn3, Slc2a2, Slc30a8 were significantly enriched. Group 2 ATs (WT AT 4; HFD AT 4; Lepob/ob AT 2 and 3) displayed significantly lower expression of insulin (Ins1) and maturation markers and higher expression of dedifferentiation/immaturity33,41, hormone processing, and secretion markers (Fig. 3b and Supplementary Data 1). Together, these findings capture β cell diversity and complexity and reveal archetypal clusters of immature (Group 2), intermediate (Group 3), and mature (Group 1) β cells across WT, HFD, and Lepob/ob mice.

In silico identification of the cell-of-origin for CCK-expressing β cells in obesity
We next sought to decipher how Cck+ β cells emerge in response to obesity. Obesity increases insulin demand, augmenting β cell mass in mice and humans42,43. Lineage tracing studies in mice have shown that new β cells primarily arise from self-duplication rather than transdifferentiation from other endocrine or exocrine cell types under a variety of physiologic scenarios, including homeostasis, pregnancy, and pancreatic injury44. In contrast, upon complete β cell ablation, α (Gcg+) and δ (Sst+) cells may acquire insulin (Ins1+/Ins2+) expression45,46, supporting islet cell plasticity in β cell regeneration. Additional studies have demonstrated that exocrine duct cells can transdifferentiate into β cells upon pancreatic injury47,48. To distinguish whether Cck+ cells arise from β cell duplication versus transdifferentiation in obesity, we leveraged TrajectoryNet for in silico lineage tracing. Since TrajectoryNet outputs a trajectory for each cellular state, we can trace trajectories of pathogenic states backwards to find their cell-of-origin. We traced back the most likely paths of Cck+ Lepob/ob cells, allowing us to identify “cell-of-origin” cells, defined as the cells at the first timepoint of the learned trajectories (Fig. 3c). We then compared the proportion of each cell type within the cells-of-origin versus the distribution in measured WT cells, terming the ratio the cell type “enrichment score” (see Methods). This analysis revealed that Cck+ Lepob/ob β cells derive from β cells in the WT condition, and Cck+ Lepob/ob polyhormonal cells derive from polyhormonal cells in the WT condition (Fig. 3d). These data raise the possibility that Cck+ cells arise from pre-existing β cells rather than by transdifferentiation of other cell types.

We next used learned cellular dynamics from TrajectoryNet to decipher which β cell ATs are most likely to be the cell-of-origin for the highest Cck-expressing (Cck-hi) Lepob/ob cells (Lepob/ob AT 2; Fig. 3b). TrajectoryNet utilizes an auxiliary network to learn the relative proliferation rates of cells from the WT condition, where WT cells predicted to be proliferating have a relative growth rate of >1. Across all WT ATs, cells committed to WT AT 4 had a higher proliferation rate than all other ATs (Fig. 3e,f). This argues that immature β cells represent a regenerative subpopulation in adult mice, consistent with prior data49,50. We retrieved the cells-of-origin for the Cck-hi Lepob/ob AT 2 cells and determined the enrichment for each WT AT within this population. Overwhelmingly, WT AT 4 was the most likely source for Cck-hi Lepob/ob AT 2 cells over all other ATs (Fig. 3g,h). We next leveraged TrajectoryNet to compute the enrichment of potential intermediate state ATs along the trajectory of cells originating from WT AT 4 and terminating at Lepob/ob AT 2 (see Methods). Interestingly, immature HFD AT 4 was enriched in the middle of the trajectory, followed by Lepob/ob AT 3 (Fig. 3i). Strikingly, the four ATs within this trajectory displayed increasing Cck expression: WT AT 4 was Cck-neg, HFD AT 4 and Lepob/ob AT 3 were Cck-mid, and Lepob/ob AT 2 was Cck-hi (Fig. 3j), consistent with Cck being a key marker of β cell adaptation to obesity. Collectively, these findings suggest that a subpopulation of immature adult β cells with proliferative potential (WT AT 4) may give rise to Cck-expressing β cells in obesity.

Experimental lineage tracing confirms that CCK+ cells originate from immature β cells
To validate our in silico findings, we performed in vivo lineage tracing by crossing tamoxifen-inducible Gcg-CreERT2, Ins1-CreERT, and Ngn3-CreSERT2 mouse lines51-53 with Lepob/ob mice harboring a Cre-inducible fluorescent reporter (Rosa26LSL-TdTomato)54 to label α, mature β, and immature β cells, respectively (Fig. 4a). Neurogenin 3 (Ngn3) drives β cell maturation by binding to the promoter of β cell-specific transcription factors (Pax4 and NeuroD), and its expression is lost once β cells differentiate to a mature phenotype55. Notably, Ngn3 was enriched in WT AT 4 (Fig. 3a). We first confirmed faithful labelling of α and β cell populations by treating mice with tamoxifen at ~4 weeks of age (before significant obesity onset), isolating pancreata one week later (5 weeks of age), and performing immunofluorescence (IF) for glucagon or insulin, respectively, on tissue sections to determine overlap with TdTomato+ cells (Supplementary Fig. 10a,b). As expected, Gcg-CreERT2 did not label β cells but demonstrated high efficiency labeling of α cells (Supplementary Fig. 10a). Ins1-CreERT and Ngn3-CreERT2 both labeled β cells, with Ins1-CreERT resulting in much greater average labeling (Fig. 4b and Supplementary Fig. 10b), consistent with Ngn3+ cells representing a subpopulation of β cells at this time point.

We next evaluated the effects of obesity by treating Lepob/ob; Rosa26LSL-TdTomato;CreERT mice with tamoxifen at 4 weeks of age and performing co-IF for CCK and insulin to determine overlap with TdTomato 12 weeks post-tamoxifen (16 weeks of age). Importantly, all end glucose levels and weights were comparable across groups (Fig. 4c,d). Gcg-CreERT2 mice showed CCK expression in most of the islet. However, only ~3% of CCK+ cells expressed TdTomato (Fig. 4e,f), indicating that CCK+ cells were not primarily arising from ɑ to β cell transdifferentiation. In contrast, most CCK+ cells were labeled with TdTomato in Ins1-CreERT (Fig. 4e,f). The incomplete labeling of CCK+ cells in Ins1-CreERT mice was likely due to inefficient labeling of β cells with this CreERT line in obese mice (Supplementary Fig. 10b), which has been previously observed in lean mice52. In support of this hypothesis, we discovered a positive linear correlation between TdTomato+/insulin+ cells (a marker of labeling efficiency) and TdTomato+/CCK+ cells (Supplementary Fig. 10c). Like Ins1-CreERT, the majority of CCK+ cells were TdTomato+ in obese Ngn3-CreERT2 mice (Fig. 4c,d). Strikingly, this frequency was comparable to Ins1-CreERT despite reduced initial β cell labeling (Fig. 4b and Supplementary Fig. 10b). This indicates that Ngn3-labeled β cells expanded, consistent with the proliferative potential of immature Ngn3+ adult β cells, as we observed in silico (WT AT 4; Fig. 3e,f). Indeed, the proportion of insulin+/TdTomato+ comparing 12 weeks vs. 1-week post-tamoxifen increased significantly in obese (Lepob/ob) compared to lean (Lepob/+) Ngn3-CreERT2 mice (Fig. 4b and Supplementary Fig. 10d). Conversely, the fraction of insulin+/TdTomato+ significantly decreased in obese Ins1-CreERT mice (Fig. 4b), showing that mature β cells present prior to obesity onset were likely replaced. Neither TdTomato+ nor CCK+ β cells in the Ngn3-CreERT2 or Ins1-CreERT lines at 16 weeks of age were Ki67+ (Supplementary Fig. 10e), suggesting that islet cells were largely post-mitotic at this time point. Together, these experimental data confirm the in silico analyses demonstrating that an immature subpopulation of postnatal β cells expand and adapt to adopt a CCK+ state.

Obesity is associated with transcriptional signatures of β cell stress
To characterize transcriptional dynamics along the axis of Cck-hi Lepob/ob AT 2 development (Fig. 3h), we clustered all highly variable gene trends to Lepob/ob AT 2 into genes that are decreasing (blue) or increasing (red) over the obesity progression (Fig. 5a). We then performed gene set enrichment analysis (GSEA) for each gene cluster56,57. As expected, genes decreasing with obesity were significantly (q < 0.05) associated with gene expression regulation in β cells and regulation of insulin secretion (Fig. 5b and Supplementary Data 2a). Upregulated genes over the obesity trajectory were significantly (q < 0.05) associated with diabetes and cellular stress responses including oxidative (electron transport chain, oxidative phosphorylation) and ER stress (Fig. 5b and Supplementary Data 2b). Several additional gene sets related to the unfolded protein response and protein processing in the ER were also significantly enriched for genes increasing over the obesity trajectory (Supplementary Data 2b). Together, these results argue that obesity promotes transcriptional signatures of β cell stress and that oxidative and ER stress are hallmarks of the obesity-induced stress response.

To confirm these findings in orthologous systems, we compared our data to recently published transcriptional signatures derived from primary human islets treated with various β cell stress agents58 including drugs that induce ER stress (brefeldin A (BFA) and thapsigargin (TG)) and inflammation stress-inducing cytokines (IFNɣ, TNFɑ, and/or IL1β). Enriched gene sets in the BFA- and TG-treated (vs. control) islets – including those associated with ER stress – highly and significantly overlapped with gene sets enriched over the obesity progression (p = 1.21 x 10-48, hypergeometric test) (Fig. 5c and Supplementary Data 3). While human islets treated with cytokines showed significant overlap with the obesity signatures (p = 7.88 x 10-5, hypergeometric test), the overlapping gene sets largely represented a general stress response (Fig. 5d and Supplementary Data 3). These data suggest that obesity induces transcriptional signatures of ER stress in β cells in mice.

We further analyzed the effects of pharmacologic β cell stress induction in vivo in mice subject to multiple low doses of streptozotocin (STZ)37. STZ is taken up specifically in β cells via the Glut2 transporter and induces metabolic (NAD+ depletion) and ER stress59-62. Like obesity, STZ treatment led to hyperglycemia, decreased pancreatic Ins1 expression, and increased β cell Cck expression in mice (Fig. 5e-g). To determine the molecular mechanisms of STZ-induced Cck expression, we reanalyzed a previously published scRNA-seq dataset of β cells from mice treated with STZ vs. vehicle control37. Visualization of these data showed distinct β cell clusters corresponding to treatment groups (β-mSTZ for STZ-treated mice) and two refined subclusters (β1 and β2) within the vehicle-treated cluster, as previously described37. We characterized β cell marker gene expression and variability of these subclusters (Supplementary Fig. 11). Concordant with obesity, Ins1 and Ins2 expression decreased and Cck expression increased with STZ treatment (Fig. 5h). β1 cells displayed high expression of maturation markers and β cell identity transcription factors (TFs) (Fig. 5i). In contrast, β2 cells reflected an immature state in vehicle-treated mice, with decreased expression of insulin secretion and maturation genes. β-mSTZ cells displayed features of Lepob/ob AT 2 cells, exhibiting loss of β cell maturation markers but increased expression of Cck and dedifferentiation/immaturity, ER stress, hormone processing, and secretion markers (Fig. 5i).

Given this similarity, we sought to map β1, β2, and β-mSTZ cells directly to β cells in our obesity dataset for comparison across all transcriptional measurements. However, due to systematic technical variation from differences in cell handling in distinct batches (batch effects), the datasets were not directly comparable. To overcome this, we employed scMMGAN26, a data integration approach that uses adversarial learning to align cells from different batches based on their data geometry. scMMGAN enables batch integration and preserves biological variation for direct comparison across diverse experimental setups and stressors. We used scMMGAN to learn the best mapping of vehicle-treated β cells (β1, β2) to WT β cell ATs and β-mSTZ cells to Lepob/ob β cell ATs using our trained AAnet models (Fig. 5j and Supplementary Data 4). The archetypal mapping of the vehicle-treated β cells revealed representation of all WT ATs, where the immature β2 state was enriched for WT AT 4 (Fig. 5j,k). β-mSTZ cells also displayed representation of all Lepob/ob ATs, where Cck-hi (>95th percentile expression) STZ β cells were enriched for Lepob/ob AT 2 (Fig. 5j,k). Based on results in our obesity models (Fig. 3g,h), we hypothesized that Cck-hi β-mSTZ cells originate from immature β2 cells. To test this, we ran TrajectoryNet to trace back the most likely paths from β-mSTZ cells assigned to Lepob/ob AT 2 to the vehicle-treated β cells. As expected, TrajectoryNet-learned relative proliferation rates determined that β2 cells have a higher growth rate than β1 cells (Fig. 5l). TrajectoryNet further demonstrated that β2 cells – rather than β1 – were far more enriched within the cells-of-origin for Cck-hi β-mSTZ cells (Fig. 5m). These data argue that both physiologic (obesity) and pharmacologic (STZ) stressors induce immature β cell expansion and adaptation to a CCK-expressing state in vivo.

The murine obesity trajectory correlates with type 2 diabetes in humans and mice.
We next evaluated whether obesity mirrored other physiologic perturbations of β cell function, including age and type 2 diabetes (T2D). Specifically, we used scMMGAN to map β cells from previously published non-diabetic (ND) and T2D mouse scRNA-seq datasets33 totaling 85,129 cells from 31 datasets integrated across conditions and developmental stages (Supplementary Fig. 12a and Supplementary Data 4). Along the obesity progression axis, Ins1 and Ins2 expression from the mouse reference cells decreased, and Cck expression increased (Supplementary Fig. 12b). Predominantly, cells from diabetes models (db/db and mSTZ) mapped to the latter end of the obesity progression, and cells from normal mice – including aged ~2-year-old) mice – mapped to the WT region of the obesity progression (Supplementary Fig. 12c). These data suggest that diabetes induction, and not age, aligns to β cell adaptation in obesity.

To further refine this β cell mapping using our trained AAnet models, we annotated the archetypal assignment of each cell and calculated the proportion of each AT for all datasets (Supplementary Fig. 12d). All ATs were represented in at least one ND dataset and at least one T2D dataset, and most ATs were represented in multiple datasets. Additionally, datasets from the same condition showed highly similar archetypal proportions. T2D datasets were predominantly enriched for Lepob/ob ATs (Supplementary Fig. 12d). ND mice from different developmental stages showed differences in archetypal commitment. Embryonic ND samples had a significantly higher proportion of Lepob/ob AT 2 versus other samples (Supplementary Fig. 12e-g), suggesting Cck-hi cells result from dedifferentiation towards an embryonic state. In contrast, postnatal P16 mice (when β cells begin to mature) had a significantly higher proportion of WT AT 4 (Supplementary Fig. 12e-g), an immature state. 2-month-old mice were enriched for archetypal states representing an intermediate maturation phenotype (WT AT 3 and HFD AT 2); 4-6-month-old mice were enriched for a mature phenotype (WT AT 1,5,7 and HFD AT 1); and 20-month-old mice were enriched for intermediate maturation ATs (WT AT 2 and 6; Supplementary Fig. 12g). These data argue that aging alone does not mimic the effects of obesity.

To determine the relevance of our findings to human biology, we used scMMGAN to map five published human scRNA-seq datasets consisting of 15 ND and 9 T2D donors to all cells from WT, HFD-fed, and Lepob/ob mice (Supplementary Data 4). The mapped insulin (INS) expression for human donors decreased along the obesity trajectory, recapitulating β cell changes along the mouse obesity progression (Fig. 6a). We annotated the archetypal assignment of each cell for each donor and calculated the proportion of each AT (Fig. 6b). All ATs were represented in at least one ND sample, and all ATs (except WT AT 2) were represented in at least one T2D sample. All ATs were found in both male and female donors, with most present in multiple ND and T2D donors (Supplementary Fig. 13). When comparing archetypal proportions between ND and T2D donors, mature β cells (WT AT 1) were notably more abundant in ND donors. In contrast, Lepob/ob AT 2 and Lepob/ob AT 3 cells were present in significantly greater proportions among T2D donors (Fig. 6c). To further assess the relationship between measured clinical variables in human donors (age, sex, body-mass index (BMI), and T2D status) and the obesity trajectory in our dataset, we leveraged the patient embedding approach DiffusionEMD27. As clinical variables are determined on the level of each donor, rather than each individual cell, it is not sufficient to compare the relationship between each variable and the cellular progression axis. DiffusionEMD overcomes this problem by using optimal transport to map an embedding of donors based on the distribution of each donor’s cells on the obesity trajectory. That is, donors with similar cell distributions are embedded close together, and donors with very different cell distributions (e.g., opposite sides of the obesity progression axis) are embedded far apart. The resulting DiffusionEMD embedding revealed a linear trajectory across donors, where donors on the left are mapped early on the obesity trajectory (high INS), and donors on the right are mapped later along the trajectory (low INS) (R=-0.93) (Fig. 6d). Visualizing clinical variables showed a low positive to negligible correlation between age (R=0.28), BMI (R=-0.18), and sex (R=-0.17) (Fig. 6e) with the obesity progression axis. Instead, we observed much stronger correlations between obesity progression and HbA1c (R=0.50) and T2D status (R=0.49) (Fig. 6e). These results suggest that the β cell adaptations associated with obesity progression in mice are concordant with T2D across species, donors, and scRNA-seq datasets.

JNK/cJun signaling modulates CCK expression in β cells
We next sought to determine the regulatory mechanisms that govern β cell Cck expression in obese mice. Transcription factors (TFs) previously associated with Cck regulation predominantly decreased in expression along the obesity trajectory except for Jun, Fos, and Crem, which increased (Fig. 7a). To predict which TFs regulate Cck expression, we leveraged Cflows, a framework that combines neural ODE networks with Granger causality to infer gene regulatory interactions from inferred continuous trajectories28. We first computed Granger causality scores between highly variable TFs and genes and subsequently pruned to TF-target interactions with high Granger causality and prior evidence for transcriptional regulatory interactions in a reference database63 (see Methods). We then subsetted to the top 100 regulatory TFs (based on Granger causality scores over all targets) and their targets with high Granger scores (Supplementary Data 5). For visualization, we further subsetted to only TF-target interactions where the TF and the target both increased in expression along the obesity trajectory. This resulted in a subnetwork of 39 interactions across 36 TFs and targets, including Cck (Fig. 7b). The subnetwork was enriched (q < 0.05) for the AP-1 transcription factor network, including Jun and Fos directly connected to Cck, as well as Egr1, Act1, Ccnd1, and Cdk1 within the subnetwork (BioPlanet). CCKR signaling map (Panther P06959, e.g., Cck, Jun, Fos, Clu, Itgb1, and Ccnd1) was also significantly enriched (Supplementary Data 6).

Jun encodes the TF cJun, which is activated and stabilized by phosphorylation by cJun N-terminal kinases (JNKs). JNK/cJun signaling is a stress-responsive mitogen-activated protein kinase (MAPK) pathway that can be induced by oxidative and ER stress64,65, raising the possibility that β cell Cck transcription may be regulated by JNK/cJun signaling in vivo. Consistent with this hypothesis, JNK inhibition (JNKi) decreased pancreatic Cck expression in Lepob/ob mice proportional to the level of total and phosphorylated cJun (Fig. 7c,d and Supplementary Fig. 14a). Similarly, we observed that JNKi can reduce Cck expression in Min6 insulinoma cells in vitro (Supplementary Fig. 14b,c). To determine whether cJun directly regulates Cck expression, we performed cJun CUT&PRUN analysis in Min6 cells, which revealed a reduction in overall peak number and intensity with JNKi (Supplementary Fig. 14d), consistent with reduced cJun levels (Supplementary Fig. 14b). As expected, cJun principally bound promoters or distal intergenic regions (Supplementary Fig. 14e), and the most frequent binding motif (27.2%, p = 5.4 x 10-18, E = 3 .8 x 10-16, binomial test) matched the canonical AP-1 binding site (5’-TCA[GC]TCA-3’). cJun was enriched at a region ~3 kb downstream of the 3’ end of the Cck gene, the peak signal intensity of which was selectively reduced with JNKi (Fig. 7e). Data from ENCODE66,67 derived from Cck-expressing tissues (mouse cortex and hippocampus) suggested that this region may be an enhancer marked by H3K27Ac and DNA hypersensitivity (Fig. 7e). Furthermore, analyses of human islet ATAC-seq68, H3K27Ac ChIP-seq68, and Hi-C69 data confirmed that this putative enhancer region is conserved in humans and resides at the boundary of a chromatin loop (Supplementary Fig. 14f). This higher-order chromatin structure may provide an additional level of transcriptional regulation of Cck. These data indicate that Cck expression may be mediated by stress-responsive JNK-cJun signaling via a conserved 3’ distal enhancer.

Functional β cell depletion suppresses pancreatic exocrine tumorigenesis in lean mice
Given the importance of β cells in obesity-associated PDAC development, we further tested whether functional β cell depletion could alter exocrine tumorigenesis even in lean mice. We crossed the rapidly progressive KPC (Pdx1-Cre; KrasLSL-G12D; Trp53LSL-R172H/+) PDAC model with Akita mice (Fig. 8a), which harbor a mutation in the Ins2 gene that causes protein misfolding, ER stress, and, consequently, non-inflammatory β cell loss70. We confirmed decreased β cell mass and induction of ER stress in islets of KPC-Akita mice without changes in islet immune cell infiltration (Fig. 8b,c). KPC-Akita mice exhibited significantly reduced circulating insulin and C-peptide levels and early-onset hyperglycemia compared to KPC controls (Fig. 8d,e), consistent with functional β cell depletion. All KPC and KPC-Akita mice developed PDAC without evidence of liver metastatic disease at 12 weeks of age. Importantly, there were no stage-specific differences in proliferation, immune cell infiltration, or fibrosis comparing PanIN and PDAC lesions between these models (Fig. 8f). Nonetheless, KPC-Akita mice exhibited significantly reduced overall disease burden compared to age-matched KPC littermates (Fig. 8g), arguing that tumor initiation rather than progression was partially compromised. These data showcase a basal role for functional β cells in PDAC development, even in lean mice.

Peri-islet acinar cell adaptation and enhanced tumorigenesis in obese mice
Finally, we sought to determine how β cell adaptations to obesity influence exocrine cells to drive tumor development by analyzing transcriptional changes in acinar cells, the putative cell-of-origin for PDAC11,12,16. We re-embedded scRNA-seq data for Cpa1+/Prss2+ cells and recovered 260 acinar cells that passed quality control metrics from all samples (Supplementary Fig. 15a). We suspected that these cells were enriched in those found in the peri-islet area, as they were obtained from highly pure islet isolations and represented ~1% of sequenced cells. Consistent with this hypothesis, acinar cells from obese (HFD-fed and Lepob/ob) versus lean (WT) mice were enriched in gene signatures previously observed in peri-islet acinar cells of the leptin receptor-deficient (Leprdb/db) obesity model71,72 (Supplementary Fig. 15b). Differential expression analysis revealed increased expression of multiple proteases in acinar cells of HFD-fed and Lepob/ob mice (Supplementary Fig. 15c), including those (Try4, Try5, Try10) confirmed to be upregulated in peri-islet acinar cells in obese Leprdb/db mice by single molecule fluorescence in situ hybridization (smFISH)72. Furthermore, a subset of these proteases (Cela3b, Cpb1, Ctrb1) were previously shown to be induced by CCK signaling in acinar cell models72, raising the possibility that this may be a consequence of β cell-derived CCK.

We next clustered acinar cells into three subclusters, the proportions of which were altered by obesity (Fig. 9a,b). Obesity was associated with a decrease in Tff2+ acinar cells (cluster 1) (Fig. 9c,d and Supplementary Fig. 15d), which have been shown to be resistant to oncogenic Kras-induced transformation73. Conversely, obesity was associated with an increase in the proportion of acinar cells expressing higher levels of proteases (clusters 2 and 3) (Fig. 9c and Supplementary Fig. 15d). Strikingly, cluster 3 was uniquely present in obesity models (Fig. 9b), and cells within this cluster exhibited high expression levels of the Regenerating (Reg) gene family (Fig. 9c,e and Supplementary Fig. 15d). Reg expression was also observed in a published analysis of peri-islet acinar cells in human pancreas samples74, the gene signature of which was enriched in cluster 3 versus cluster 1 (Fig. 9f), arguing that acinar cells from cluster 3 harbor features of human peri-islet acinar cells. Reg expression is associated with acinar cell stress and inflammation75, and Reg proteins can promote ADM76,77, making them both markers and drivers of a tumor-permissive acinar cell state. In concordance with our scRNA-seq data, bulk pancreata of obese Lepob/ob mice exhibited upregulation of Reg genes compared to lean congenic WT mice (Fig. 9g). Importantly, Reg expression was abolished by chronic CCK receptor antagonist treatment of Lepob/ob mice (Fig. 9g), indicating that cluster 3 cells may be induced by paracrine signaling from β cell-derived CCK in mice. Consistent with this hypothesis, we observed increased peri-islet acinar cell Reg2 protein expression in KC;Lepob/+;Ins1-Cck/+ and KCO mice compared to KC;Lepob/+ controls (Fig. 9h). Together, these data argue for the potential of β cell-derived CCK to promote a tumor-permissive peri-islet acinar cell state.

We reasoned that these CCK-dependent alterations in peri-islet transcriptional states may enhance Kras-driven tumorigenesis in proximity to islets in obesity. In two-dimensional analyses of tissue sections, we found that greater proportions of PanINs in obese KCO mice emerged in proximity to islets than in lean KC;Lepob/+ mice (Fig. 10a). To confirm whether the spatial relationships between islets and neoplasia were preserved in three dimensions (3D), we employed tissue clearing and immunofluorescence, coupled with 3D light-sheet microscopy imaging of intact pancreata of obese KCO mice and lean KC;Lepob/+controls at timepoints with comparable disease burden (8- and 12-weeks-old, respectively; Fig. 10b and Supplementary Movies 1-2). 3D reconstructions of the tissue (Fig. 10c) identified 329 tumors in lean KC;Lepob/+ mice and 281 in obese KCO mice with no significant difference in relative disease volume (Fig. 10d). Although islet radius was significantly larger in KCO mice (Fig. 10e), as previously reported for obese mice78, fewer total islets were identified (271 vs. 444 in KC;Lepob/+ mice). In support of islet proximity dependent tumorigenesis, both the neoplasia-to-islet and islet-to-neoplasia distances were significantly shorter in obese KCO mice (Fig. 10f). Cumulative frequency analysis confirmed a higher proportion of islets located closer to neoplastic lesions in obese mice (Fig. 10g), supporting endocrine-exocrine CCK signaling in acinar cell adaptation and obesity-associated pancreatic tumorigenesis.

Methods
All research performed in this study complies with relevant ethical regulations including animal studies approved under the Yale University Institutional Animal Care and Use Committee (IACUC) protocol #20170 and Yale Environmental Health and Safety recombinant DNA protocol #17-35.

Animal studies
Mice were housed in a specific-pathogen free (SPF) facility at room temperature and humidity with standard 12-hour day/night cycles. Akita (Stock #003548), KrasLSL-G12D (Stock #008179), Lepob/ob (Stock #000632), Ins1-CreERT (Stock #024709), Gcg-CreERT2 (Stock #042277), Ngn3-CreERT2 (Stock #028365), Pdx1-Cre (Stock #014647), Rosa26LSL-TdTomato (Stock #007914), Pgk1-FlpO (Stock #011065), and Trp53LSL-R172H (Stock #008652) mice were obtained from the Jackson Laboratory (JAX). Ins1-Cck mice were a gift from D.B. Davis (University of Wisconsin). Cryorecovery of Gcg-CreERT2 and Ngn3-CreERT2 embryos was performed by the Yale Genome Editing Center (YGEC). CCK tm2a mice (EuCOMM) sperm was obtained from Riken, and mice were generated by in vitro fertilization (IVF) by the YGREC. Cckflox mice were generated by crossing progeny with Pgk1-FlpO strain, which was later crossed out. KPC (+/- Akita), KCO (+/- Ins1-Cck or Cckflox), and CreERT;Lepob/ob;Rosa26LSL-TdTomato mice were generated by breeding. Genotyping of the mice was done via PCR using genomic DNA from the tail or ear as a template. Tail DNA was isolated using Hotshot extraction, and PCR was run using GO Taq Green Mastermix (Promega). Primers and protocols used were obtained from JAX or previously described18 with the exception of Cckflox mice for which the following primers were used: forward primer: 5’-ATGGGCCCACACTTAAAACG-3’, reverse primer: 5’-AGCCAAGCGAAGTTCCTAGA-3’. Primers were obtained from the W.M. Keck Biotechnology Resource Laboratory at Yale. For Lepob/ob breedings, Lepob/ob mice were administered AAV2/1-CAG-mLeptin via intramuscular injection once, which restores fertility18. For bulk RNA-seq experiments, Lepob/ob and congenic age-matched C57/B6 WT mice (Stock #000664) were purchased from JAX. Lepob/ob mice were treated with CCK receptor antagonists (proglumide (0.1 mg/mL; Sigma) or lorglumide sodium salt (0.05 mg/mL; Sigma)) or control drinking water continuously for 6 weeks starting at 6 weeks of age. Mice were euthanized at 12 weeks of age and isolated pancreata were snap frozen for RNA-seq analysis, as described below. CreERT mice were induced at ~4 weeks of age with 8mg/40g (obese Lepob/ob) or 6 mg/40g (non-obese Lepob/+) of tamoxifen (Sigma) in corn oil via oral gavage for 4 consecutive days. Mice were euthanized 1- or 12-weeks following start of tamoxifen administration for analysis. For STZ experiments, male mice were treated with 50 mg/kg per day of STZ or PBS (vehicle) by intra-peritoneal injection for 5 consecutive days and analyzed 7 weeks afterwards. For JNKi experiments, 17-week-old Lepob/ob mice were treated with either 20mg/kg SP-600125 (MedChemExpress) or vehicle control (5% DMSO, 40% PEG-300, 5% Tween 80, 50% ddH2O) by oral gavage once daily for five days prior to dissection. All mice were fed standard chow (Inotiv, 2018S) ad libitum except for mice high-fat diet (HFD), which received a 60% kcal lard-based HFD (Research Diets, D12492) for 10 weeks. Mouse weights were measured using a small animal scale, and blood glucose was measured using a OneTouch Ultra2 glucometer. Serum was collected at time of dissection under fed conditions unless otherwise noted. For fasted conditions, food was removed from the cage for 6 hours prior to sample acquisition. Serum C-peptide was measured by ELISA (Crystal Chem 90050) according to manufacturer’s instructions. Serum insulin was measured by radioimmunoassay (RIA, Linco RI-I3K) by the Yale Diabetes Research Center. Mice for KCO tumor studies were kept in a hybrid C57/B6:129 background. Mice for KPC tumor studies were kept in a C57/B6 background. Specific endpoints are listed in the figure legends and were selected based on prior studies on tumor progression in these models18,103. Sibling littermates were used as controls in all tumor studies. The maximal tumor size permitted was 1 cm3, and this was not exceeded in the study. Mice of both sexes were used for all experiments (and denoted in the figures with symbols) except for the STZ and scRNA-seq experiments for which only male mice were used due to reduced β cell effects of STZ and obesogenic effects of HFD, respectively, in female mice.

Cell culture
Min6 cells (at a passage of 31-50, obtained from Dr. G. Cline) were grown under high glucose DMEM (Corning) supplemented with 15% fetal bovine serum (FBS; ThermoFisher Scientific), penicillin-streptomycin (1%; ThermoFisher Scientific), L-glutamine (2mM; ThermoFisher Scientific), sodium pyruvate (1mM; ThermoFisher Scientific), and 2-mercaptoethanol (0.05mM; Sigma). Cells tested negative for mycoplasma by PCR performed by the Yale Molecular and Serological Diagnostics Core.

Immunoelectron microscopy
Cultured Min6 cells were fixed in 0.1% glutaraldehyde and 4% paraformaldehyde (PFA) in 0.1M sodium cacodylate buffer (pH 7.4) for 30 minutes at room temperature, followed by 4% PFA in the same buffer for one hour. After buffer rinse, cells were scraped off the culture dish and gently centrifuged. Pelleted cells were resuspended in 10% gelatin and submerged in 2.3M sucrose in PBS at 4C overnight prior to freezing in liquid nitrogen. 60-nm thick sections were cut using a Leica EM UC6/FC6 cryo-ultramicrotome and collected on nickel grids. After warming to room temperature, sections were immunolabeled with primary antibodies for 30 minutes, followed by 30-minute incubation with immunoglobulin G antibodies coupled to 5 nm and 10 nm protein A-gold particles (Utrecht University Medical Center, Netherlands). Specific antibodies used are listed in Supplementary Table 1. Images were acquired using a FEI Tecnai Biotwin transmission electron microscope (80 kV) and an AMT NanoSprint15 cMOS camera.

Western Blot
Mouse pancreata were snap-frozen in liquid nitrogen and manually pulverized using the BioPulverizer (BioSpec Products 59012MS). Pancreatic tissue (~40 mg) and Min6 whole cell lysates were harvested in RIPA buffer (Pierce), EDTA (1:100) (ThermoFisher Scientific), and Halt protease and phosphatase inhibitor cocktails (1:100) (ThermoFisher Scientific). Tissue samples were vortexed and homogenized using a pellet mixer (VWR) and rotated at 4C for one hour prior to centrifugation at maximum speed for 20 minutes. Cells were scraped and placed on a rotator at 4C for 15 min and then spun as above. Protein concentrations were measured using the BCA Protein Assay Kit (Pierce) following the manufacturer’s protocol. Equal amounts of protein (30 µg) were loaded to each well of mini-PROTEAN 4-20% TGX stain-free precast gels (Bio-Rad). Protein was transferred to nitrocellulose membranes using the TransBlot Turbo Transfer System (Bio-Rad). For pancreatic tissue lysates, membranes were stained with Ponceau S (0.5% (w/v) in 1% glacial acetic acid) to verify equal protein loading and then washed with PBS-T (PBS, Tween-20 (0.1%)). Membranes were blocked using Odyssey Blocking Buffer (Licor) in PBS (1:2) for one hour, then incubated with primary antibodies diluted in Odyssey Blocking Buffer and PBS-T (1:2) overnight at 4C. Membranes were washed for 10 minutes with PBS-T three times and then incubated in DyLight secondary antibodies (Cell Signaling Technologies) for one hour at room temperature. Membranes were washed with PBS-T for 10 minutes, three times each, prior to imaging. Specific antibodies used are listed in Supplementary Table 1. Membranes were imaged using the ChemiDoc Touch Imaging System (Bio-Rad). Band intensity was quantified using ImageLab software.

RNA Extraction, qRT-PCR, and bulk RNA-sequencing
Mouse pancreata were snap frozen in liquid nitrogen and manually pulverized using the BioPulverizer (BioSpec Products 59012MS). RNA was extracted from pulverized tissue (20-40 mg) and Min6 cells, homogenized via QIAshredder column (Qiagen), and purified using the Qiagen Mini RNeasy kit following the manufacturer’s protocol with DNAse I treatment to remove any genomic DNA contamination. RNA quality and quantity was assessed using an Agilent Bioanalyzer, and only high-quality RNA samples (RIN>7) were used in downstream analysis. 1 g of RNA was used to synthesize cDNA using the Applied Biosystems™ High-Capacity cDNA Reverse Transcription Kit. cDNA was diluted in nuclease-free water (1:5), and qPCR was run using TaqMan™ Universal PCR Master Mix and FAM probes (ThermoFisher Scientific) for Gapdh (Mm99999915_g1), Ins1 (Mm01950294_s1), Syn (Mm00436850_m1), and Cck (Mm00446170_m1). qPCR assays were run on a CFX Opus 384 Real-Time PCR System (Bio-Rad). RNA-seq libraries were prepared with polyA selection (Illumina) and sequenced (>25M 100-bp paired-end (2x100) reads on NovaSeq (Illumina)) with the Yale Center for Genome Analysis (YCGA). Low-quality bases and adapter sequences were trimmed using Trim Galore (v0.6.7) and CutAdapt (v3.5). Cleaned reads were aligned to the UCSC mouse genome (mm10) using the STAR aligner (v2.7.9). Differential gene expression analysis was performed using DESeq2 in R, and genes with an adjusted p-value < 0.05 were considered differentially expressed.

CUT&RUN analysis
Min6 cells were plated (2.5 million) and grown for a total of three days. The day after plating, cells were treated with the JNK inhibitor SP600125 (20 µM, MedChemExpress) or DMSO (control) for 48 hours. 550,000 cells were harvested per reaction and condition. Cells were processed using the CUTANA ChIC/CUT&RUN Kit (14-1048), JUN/cJun CUTANA CUT&RUN antibody (SKU:13-2019), and positive (anti-H3K4Me3) and negative (IgG) control antibodies from the kit following the manufacturer’s protocol. Library preparation was performed by the YCGA followed by sequencing to a depth of 30 million paired-end 150 bp reads (2x150) per sample on NovaSeq (Illumina). Adapter sequences from paired-end sequencing reads were trimmed with Trimmomatic and quality control was performed using FastQC. Bowtie2 (v2.4.2) was used to align reads to the GRCm38 genome using the –very-sensitive flag and allowing for one mismatch. SAMtools version 1.16 was used to remove unmapped reads, secondary alignment reads, and reads with a Phred score of <30. Macs2 (v2.2.7.1) was used to call narrow peaks with default peak-calling settings. H3K4me3 and cJun CUT&RUN peaks were called against IgG negative controls. IDR was used to determine significant peaks between the two biological replicates. Deeptools (v3.5.1) was used for data visualization. To visualize average CUT&RUN signal per sample, Deeptools bamCoverage was used to normalize bigWig signal based on E. coli spike-in control DNA. A scaling factor was calculated based on the percentage of reads mapped to the E. coli K12, MG1655 reference genome. The following additional bamCoverage flags were used: --extendReads --ignoreDuplicates --samFlagInclude 64. Differences in read coverage between samples were normalized based on CPM. In all plots, the mean normalized bigWig signal between two biological replicates is plotted. Peak intersections were performed using BEDTools (v2.30). Enriched motifs were identified using MEME v5.4.1 meme-chip. For all published and publicly available datasets, BED and/or bigWig files were converted to the mm10 genome or hg19 genome for mouse and human sequencing datasets, respectively. Peaks and bigWig signal were visualized with the Integrated Genomics Viewer (IGV). Hi-C data were visualized with the WashU Epigenome Browser104-107.

Tissue preparation for histology, immunohistochemistry, and immunofluorescence
Mice were euthanized by CO2 asphyxiation. Pancreata was isolated, fixed in 10% neutral-buffered formalin overnight, dehydrated in 70% ethanol, and embedded in paraffin by Yale Tissue Pathology Services. Adjacent 5-μm sections were cut and stained with hematoxylin and eosin (H&E) or subject to immunohistochemistry (IHC) following citrate-buffer antigen retrieval using a ThermoFisher Scientific Autostainer 360 or manual staining (for CCK). Mach 2 HRP-labeled micro-polymers (Biocare Medical) or HRP-conjugated secondary antibodies (Jackson ImmunoResearch or Vector Labs) were used for primary antibody detection. Tissues were scanned using an Aperio whole slide scanner at 20X magnification. For immunofluorescence (IF) analyses, mice were fixed with 4% PFA (Electron Microscopy Sciences) via cardiac perfusion. The pancreas, small intestines, and liver were harvested and fixed in 4% PFA overnight at 4C, washed in PBS three times for 5 min, and cryoprotected in 30% sucrose in PBS at 4C. Tissue was embedded in O.C.T. (Tissue-TEK) and stored at -80C prior to cryosectioning. 10-μm cryosections were obtained using a Leica Cryostat (CM1860). Slides were thawed at room temperature for one hour and washed three times for 10 min with PBS. Tissue was permeabilized and blocked for one hour with PBT (PBS, TritonX-100 (0.3%)) with 10% donkey serum (Jackson ImmunoResearch). Slides were incubated at 4C with primary (Day 1) and secondary (Day 2) antibodies (Jackson ImmunoResearch) overnight in PBT and 5% donkey serum in a humidified chamber and washed three times for 10 min after each antibody incubation with PBT. On Day 3, slides were stained with DAPI (1:500; ThermoFisher Scientific) in PBS and washed once with PBS. Slides were mounted using Vectashield Hardset Antifade Mounting Medium (VectorLabs). Slides were imaged at the Yale West Campus Imaging Core using a Leica SP8 Spinning Disk Confocal Microscope and an Andor Benchtop Confocal Microscope BC43. Monochromatic red, green, and blue images were merged using ImageJ software (NIH). Specific antibodies used for IHC and IF analyses are listed in Supplementary Table 1.

Whole-mount immunostaining, tissue clearing, and light sheet microscopy
Mouse pancreatic tissue was cleared using a modified iDISCO protocol108,109. Mice were transcardially perfused with 1× PBS followed by 4% PFA containing 10 μg/mL heparin. Pancreata were harvested and post-fixed overnight at 4 °C in 4% PFA. Samples were dehydrated through a graded methanol series (20%, 40%, 60%, 80%, and 100%) in B1N buffer (300 mM glycine, 0.1% Triton X-100, 0.01 N NaOH, 0.02% sodium azide in distilled water; 1 hr per step). Delipidation was performed in 100% dichloromethane (DCM), followed by three 30 min washes in 100% methanol. Tissues were bleached overnight at 4 °C in 5% hydrogen peroxide (H₂O₂) in methanol, then rehydrated through a reverse methanol gradient in B1N buffer. Permeabilization was achieved by overnight incubation in modified PTxwH buffer (0.5% Triton X-100, 0.1% Tween-20, 2 mg/mL heparin, 0.01% sodium azide in 1× PBS) supplemented with 5% dimethyl sulfoxide (DMSO) and 300 mM glycine, followed by three 30 min washes in PTxwH. Samples were blocked overnight at room temperature in blocking buffer (3% donkey serum in PTxwH), then incubated with anti-Glucagon antibody (Sigma-Aldrich, G2654; 1:2000 in blocking buffer) for 7 days at 37 °C. After washing in PTxwH, tissues were incubated with Alexa Fluor 568–conjugated anti-mouse secondary antibody under identical conditions. Post-immunostaining, samples were embedded in 0.75% agarose, dehydrated through a methanol gradient, incubated overnight in 100% methanol, delipidated in DCM, and refractive index–matched in dibenzyl ether. Cleared, agarose-embedded samples were imaged on an Ultramicroscope Blaze light-sheet microscope (Miltenyi Biotec) equipped with a scientific CMOS camera and an MI PLAN 4× objective lens (zoom 0.60). Samples were immersed in ethyl cinnamate (refractive index = 1.55) to minimize optical aberrations. Images were acquired in ImSpector Pro v7.7.2 (Miltenyi Biotec) in multicolor mosaic mode using 568 nm excitation (Glucagon) and 488 nm excitation (autofluorescence). The light sheet was generated by dual-sided illumination with a width of 50% and thickness of 5.00 µm; z-stacks were collected at 2.5 µm steps. Tiles were stitched using Imaris Stitcher (Oxford Instruments).

Immunocytochemistry
Min6 cells were grown in a 6-well plate with a 20 mm round cover glass (Cell Treat, Part#: 229173) for 3 days, washed with PBS (2x, 5 min each), fixed with 4% PFA (15 min), and washed with PBS (2x, 5 min each) at room temperature. Cells were blocked using PBT with 10% donkey serum and incubated overnight at 4° C in primary antibody in PBS-T with 2.5% donkey serum. On Day 2, cells were washed with PBS (3x, 5 min each) and incubated overnight at 4°C in secondary antibody in PBS-T with 5% donkey serum. On Day 3, cells were washed with PBS-T (3x, 5 min each), incubated with DAPI (1:500) for 5 min, washed in PBS for 5 min, and mounted on slides using Vectashield Hardset Antifade Mounting Medium. Slides were imaged with a Leica SP8 Spinning Disk Confocal Microscope. Specific antibodies used for immunocytochemistry (ICC) are listed in Supplementary Table 1.

Pancreatic islet isolation, perifusion, and single-cell RNA sequencing (scRNA-seq)
Male 16-week-old wild-type (WT) C57/B6 mice (Stock #000664) and male C57/B6 mice fed a 60% kcal fat diet (Research Diets 12492) for 10 weeks (Stock #380050) starting at 6 weeks of age were obtained from JAX. Islets were isolated, pooled for each group (n=4 mice per group), and processed for perifusion for GSIS measurements or dispersed into single cells for sequencing, as previously described18. For scRNA-seq, dead cells were removed using a Dead Cell Removal Kit (Miltenyi Biotec, 130-090-101). Cell concentration and viability were determined using a Countess II Automated Cell Counter (ThermoFisher Scientific). 10X Genomics Chromium Single Cell 3’ Reagent Kit v3 library preparation and sequencing on Illumina HiSeq were performed by the YCGA.

Quantification and statistical analysis of tumor development
All disease analyses were performed in a blinded manner on coded scanned H&E slides in QuPath (v0.6.0). Disease burden was quantified by measuring cross-sectional area of disease (ADM, PanIN, and PDAC including stroma) relative to total pancreas area. Histological analysis of tumor progression (presence of PanIN and/or PDAC) was performed by two independent blinded reviewers (C.F.R and M.D.M.) using previously established morphologic criteria110.

Quantification and statistical analysis of immunohistochemistry
Quantitative IHC analysis for synaptophysin, insulin, Ki67, CD45, and SMA was performed in QuPath (v0.6.0) using optimized marker-specific detection settings. Islet mass was quantified by measuring cross-sectional area of synaptophysin-positive cells relative to total pancreas area on tissue sections. Proportional β cell area was calculated by cross-sectional area of insulin-positive cells relative to synaptophysin-positive cells on tissue sections. For Ki67, CD45, and SMA quantification, representative annotations for each mouse were drawn for region of interests (ROIs) containing PDAC or PanIN lesions (PanIN epithelial fields for Ki67 and SMA; PanIN epithelial and surrounding stromal fields for CD45). Annotation measurements were exported from QuPath and processed in R (v4.3) using the tidyverse package. For Ki67, positive cell density was calculated as the number of Ki67+ nuclei per mm² of annotated lesion area. For CD45, stromal infiltration in PanINs was estimated as CD45+ cells per mm² of stroma, determined by subtracting epithelial counts and areas from total field measurements; PDAC CD45 density was calculated directly from lesion ROIs. For SMA, fractional positive area was computed as SMA+ DAB area divided by total annotation area. For each marker and classification, multiple ROIs per mouse were averaged to yield one value per mouse, and data were aggregated by genotype (e.g., KPC vs. KPC-Akita) to calculate mean ± SEM. Normality was confirmed by the Shapiro–Wilk test, and genotype comparisons were performed using unpaired two-sample t-tests (α = 0.05).

Quantification and statistical analysis of tumor proximity to islets
Scanned H&E sections of pancreata from lean 13-week-old KC;Lepob/+ and obese 7-week-old KC;Lepob/ob mice (n=4-5 mice per group) were digitally analyzed using QuPath (v0.6.0) in blinded fashion. Neoplastic lesions were largely low-grade PanINs in these models at the designated timepoints. All low-grade PanINs in the splenic lobe of the pancreas (where most islets are located in mice) of each mouse were quantified and classified as islet proximal or distant. Islets (diameter >100 µm) were defined as being proximal if they resided within a distance of 400 µm from the tumor epithelial cell edge to the center of an islet. This represents the approximate length of ten acinar clusters within a lobule, which is twice the expected diffusion limit but would allow hormone delivery via the interconnected islet-acinar capillary network111,112. Fisher’s exact test was performed to compare lean and obese conditions. Despite harboring fewer overall PanINs (n=58 (14.5 per mouse) vs. 146 (29.2 per mouse)) at this earlier timepoint, KC;Lepob/ob mice exhibited a great proportion of PanINs close to islets. There was no significant difference (p=0.23, Welch’s t-test) in the average number of islets per mouse (14 for KC;Lepob/+ and 19.75 for KC;Lepob/ob). For 3D analyses, islets and neoplasia were reconstructed in Imaris v10.2.0 (Oxford Instruments). Islets were segmented based on glucagon immunofluorescence, and neoplastic lesions were identified by differences in autofluorescence intensity113. The “surface” tool in Imaris was used to quantify the shortest distances from islet to neoplasia and neoplasia to islet, as well as islet radii. Statistical comparisons of two groups were performed using the Mann–Whitney test. Cumulative distributions were compared with the Kolmogorov-Smirnov test.

In vivo lineage tracing analysis
To analyze sufficient islets at single-cell resolution for each mouse, we obtained widefield fluorescence images using a Nikon Spinning Disk Confocal (CSU-W1) under 40x or 60x oil immersion with laser wavelengths of 405-nm, 488-nm, 561-nm, and 647-nm. Fluorescence channels were merged using ImageJ, and positive cells were counted manually. Briefly, 15-20 islets were analyzed for each 16-week-old experimental mouse per line, scored for presence or absence of TdTomato, insulin, or CCK expression, and each counted cell was marked in ImageJ to avoid double counting. The percentage of CCK+ or insulin+ cells also expressing TdTomato were calculated for each mouse for statistical comparisons between mice. For mouse line validation (5-week analysis timepoint), 10 islets were analyzed per mouse. For Ins1-CreERT and Ngn3-CreERT2 mouse lines, the proportion of insulin+ cells harboring TdTomato fluorescence was quantified. For the Gcg-CreERT2 line, we calculated the percentage of Gcg+/TdTomato+ cells and divided the result by the total Gcg+ cell population. Due to unequal variance, the Welch’s t-test was used for statistical analysis when comparing groups.

Preprocessing of scRNA-seq data
Murine obesity models
scRNA-seq data from 12 week-old C57/B6 Lepob/ob and WT mice were previously published by our group18, and raw data deposited in GEO (GSE137236) were reanalyzed in comparison to new HFD and WT data from 16-week-old mice, as described above. All sample data were processed with Cell Ranger (v3.0.2) using the mm10 mouse genome indices from 10X Genomics. We combined batches from all four scRNA-seq datasets and removed rare genes (expressed in fewer than 15 cells). We then filtered and L1 normalized cells by library size (maximum library size=15000) and square-root transformed the data. We finally filtered cells with high mitochondrial expression (maximum expression=12.5). All processed cells were visualized with PHATE (v1.0.11)114, with Ins1, Ins2, Sst, Ppy, and Gcg removed from PHATE calculation due to very high expression affecting cell-cell distances. All cells were then clustered with k-means clustering (n_clusters=20) and clusters were annotated based on published mouse and human islet marker genes33 to distinguish endocrine cells and β cells for downstream analysis. Clusters 13 and 14 were annotated as low-quality clusters due to lack of enrichment for key markers of pancreatic islet cell populations and distinct embedding visualization. Cluster 13 showed low overall counts, strong enrichment for hemoglobin genes, and low ribosomal counts, indicative of contaminating erythroid cells, and cluster 14 showed low ribosomal counts and enrichment primarily for marks such as Malat1, Neat1, and Gm42418, lncRNAs whose high expression reflects high nuclear fraction and cytoplasmic RNA loss due to membrane damage115. Cluster 11 was further split into clusters 11 and 20 due to enrichment of Schwann cell markers S100b, Ngfr, and Pmp22 in cluster 20. This clustering was chosen due to stability over five runs (average pairwise ARI=0.86) and balance between overclustering and merging different cell types. β cells were then subsetted and further filtered to remove outlier populations with large differences in the number of genes expressed per cell, indicating low quality of the cell populations. Endocrine cell clusters were retained for TrajectoryNet, re-embedded, and clustered with k-means clustering (n_clusters=10, grouped for annotation).

To verify transcriptional differences between models were driven by condition-specific effects (WT vs. HFD-fed vs. Lepob/ob) rather than batch effect or age-based differences, we calculated pairwise Earth Mover’s Distances (EMD) between the four samples, where a smaller distance reflects a higher similarity between cells in each sample. In endocrine cells, measured from the principal component (PC) space used for trajectory inference and analysis (see below), we confirmed the two WT samples (12- and 16-week-old) were more similar to each other (EMD=0.45) than to their respective age-matched samples from obese mice (EMD=0.76 and 0.80). Furthermore, both samples were closer to the HFD-fed sample than the Lepob/ob sample (EMD=0.63 vs. 0.76 for 12-week-old WT and 0.80 vs. 1.26 for 16-week-old WT). This trend held in β cells analyzed in the PC space, where the two WT samples (12- and 16-week-old) were more similar to each other (EMD=0.59) than their respective age-matched samples (EMD=0.94 and 0.66). Both samples were closer to the HFD-fed sample than the Lepob/ob sample (EMD=0.59 vs. 0.94 for the 12-week-old and 0.66 vs. 0.84 for 16-week-old). Finally, we computed the modified average silhouette width (ASW) of batch, as previously implemented116. This metric measures the silhouette of a given batch for each cell type. For the scaled metric we used here, 0 indicates separation of batches, and 1 indicates perfect overlap of batches. In the combined β cell PC space, the two age-matched batches had a modified ASW of 0.71, and the two WT samples (subsetting to WT only) had a modified ASW of 0.84. Together, these results suggest higher concordance between the WT samples than their respective age-matched samples.

Given the similarity between WT samples, we merged cells from both samples for all downstream analyses. A total of 17,336 β (Ins1+/Ins2+) cells and 260 acinar (Cpa1+/Prss2+) cells were retained. For gene cluster analysis, GSEA of gene clusters, and Granger causality analysis, we used highly variable genes defined using default parameters in Scanpy (v1.9.3)117. Endocrine and acinar cell embeddings were generated with PHATE. β cell embeddings were generated by first denoising with MAGIC (v3.0.0)118, then reducing dimensionality with PCA. This enabled visualization of heterogeneity similar to PHATE while maintaining invertibility to the gene space for trajectory and latent-space archetypal analysis.

Published datasets
Upregulated (log FC expression>1.5, FDR<0.05) acinar-specific genes from RNA-seq analysis of islets isolated from Leprdb/db versus WT mice71 were extracted from “Table S3” of Egozi et al.72. The gene signature associated with the acinar-REG+ population were derived from differentially expressed genes identified in Tosti et al72. Differentially expressed genes for patients exposed to pharmacologic β cell stressors were extracted from the Gene Expression Omnibus (GEO, GSE237448)58. We subsetted to DEGs per stressor that were enriched versus control (CTRL) in the majority of patients (at least 3 of 5). We then computed the gene sets for DEGs per stressor. For BFA and TG, we took the intersection of gene sets before comparing overlap with obesity. For CM, IFN, IL1B, IL1BIFNG, and TNFA, we took the intersection of gene sets before comparing overlap with obesity. All gene sets enriched per stressor and associated genes are provided in Supplementary Data 3.

Pre-processed and annotated scRNA-seq data from islets derived from mice treated with STZ (mSTZ) or control (vehicle) were obtained from GEO (GSE211799)37 and subsetted to β cells. β1, β2, and β-mSTZ annotations were used from original work, except for six control cells labeled as β-mSTZ, which were relabeled as β2 for cell-of-origin analysis to differentiate between an immature pre-existing population and a differentiated cell state.

Human islet datasets were harmonized across measured clinical variables, subsetted to β cells from ND and T2D donors, and pre-processed for library size normalization and square-root transformation. The mouse atlas33 was subsetted to only β cells from control/ND and type 2 diabetic (T2D) mouse models. Samples with fewer than 2 cells were excluded. Samples annotated “Embryonic” show embryo progression from E12.5 to E15.5. Samples annotated as P16, 4m, and aged were islet cells isolated from 16 day-old, 4-month-old, and 2-year-old mice and sorted according to the Fltp lineage-tracing model. mSTZ samples are derived from healthy adult control mice, an STZ-treated model, and mSTZ models subject to different antidiabetic treatments. db/db samples are from healthy adult control, Leprdb/db, and Leprdb/db subject to different antidiabetic treatments. Chem samples are healthy young adult control samples (only control retained for this analysis).

TrajectoryNet
Overview
Traditionally, characterizing progression axes across cells relies on “pseudotime” analysis, where each cell is assigned a given time value corresponding to its relative state on the main axis or axes of variation. Ordering of cells is determined by the choice of a starting cell or inferred by the method. However, due to the destructive nature of scRNA-seq, inferring the trajectory of each cell across timepoints remains experimentally challenging, and the inferred or chosen starting and ending timepoints may not reflect the true latent dynamics of the system. This motivated the development of TrajectoryNet, a manifold-aware neural ODE network that performs dynamic optimal transport between timepoints23,28. TrajectoryNet learns individual trajectories for each cell from the first timepoint (WT) to the last timepoint (Lepob/ob) by optimizing for distribution matching between the predicted and true timepoint, as well as energy efficiency, which enables dynamic optimal transport of cells. Additionally, TrajectoryNet is regularized for unbalanced dynamic optimal transport through modeling cellular proliferation and death with an auxiliary proliferation network. Together, this ensures TrajectoryNet can learn constrained continuous normalizing flows for each cell that are biologically plausible.

TrajectoryNet-based cell-of-origin definition
Each individual trajectory learned by TrajectoryNet defines the trajectory of a given cell in the last timepoint, where TrajectoryNet generates the most likely cells along the path to the first timepoint. Thus, the generated cells at the first timepoint represent the “cells-of-origin” for the learned trajectories. Notably, in addition to distribution matching, TrajectoryNet ensures cellular paths are energy efficient and models biological proliferation and death, meaning that the learned cells-of-origin may be enriched for a particular subset of cells within the first timepoint (in our dataset, cells in the WT population).

Cflows for TrajectoryNet-based granger causality analysis
Using the top 2,278 highly variable genes (defined in Preprocessing of scRNA-seq data above), we first computed the mean gene trends over time with TrajectoryNet for two runs averaged. For each transcription factor (TF) in the TRRUST v2 database63 and all highly variable genes, we computed the p-value of the Granger causality test, or the time-lagged correlation between the TF and the target with lag of 10 (10% of the trajectory), where the p-value is computed based on the chi square distribution. We then convert the p-value into a Granger causality score, as described in the Cflows framework28, by taking the log and adding a sign based on the correlation between two gene trends.

For building a regulatory network of genes strongly associated with the obesity progression trajectory, we first identified the top 100 transcription factors with the most regulatory effect (i.e., the highest Granger score across all highly variable genes) and subsetted all TF-target relationships from TRRUST to those TFs. Then, within this subset, we further filtered TF-target pairs to only those with high absolute Granger score (i.e., -log (p-value) > 25). We considered this network the full regulatory network. For visualization, we extracted those pairs where both genes are considered “increasing” over the trajectory based on the gene cluster analysis. We visualized the largest connected component, which includes 28 out of 36 genes in the “increasing” subnetwork.

Experimental details
We learned cellular trajectories to all cells in the Lepob/ob condition for all endocrine cells and to cells assigned to Lepob/ob AT 2 for the β cells in our obesity models and in the mSTZ model. TrajectoryNet (v0.2.4) was run twice in each setting, where enrichment scores (see below) were calculated separately and plotted together, and the mean expression from both runs was computed for visualizing gene trends, calculating gene clusters, and downstream gene set enrichment analysis. One run was used for visualizing individual cellular trajectories. We ran TrajectoryNet on embeddings built from running MAGIC118, then PCA (to ensure denoised trajectories as well as invertibility to the gene space), and with the proliferation regularization with alpha = 2 (other parameters default). TrajectoryNet proliferation rate was robust to choice of alpha within [0.01,10] (pairwise Spearman ρ between 0.81-1.0).

scMMGAN
Overview
Mapping one dataset to another or generating a representation of one dataset in another modality, species, condition, or strain, allows us to build a broader understanding of the broad relationships between cells from different settings toward identifying common gene signatures. scMMGAN is a generative adversarial network architecture for domain adaptation, where the model is trained to preserve relationships between data points and ensure the original data can be reconstructed from the generated output26. Unlike batch correction methods, scMMGAN allows us to map β cells from different experimental setups without the obesity progression manifold.

Experimental details
For all runs, we ran scMMGAN (v0.0) with 15,000 training steps, batch size of 128, learning rate of 0.0001, and λ of 1. The correspondence loss for scMMGAN can be flexibly defined. For all mappings, we enforce the preservation of the Pearson correlation of a cell prior to mapping and after mapping to the new domain. This ensures maintenance of the cell’s identity for accurate mappings. For mapping vehicle to WT cells and mSTZ to Lepob/ob cells, we ran scMMGAN after MAGIC, then PCA (to match the preprocessing of each condition for AAnet). We used λ of 10, and scMMGAN was run twice to show consistency of enrichment correspondence scores. For mapping the mouse atlas and the human samples, we ran scMMGAN after PCA (for the batch-integrated mouse atlas in 1 run, and for the human datasets individually due to batch effect). Due to the larger differences between these datasets and our obesity progression, versus the mSTZ dataset, we used λ of 15 to increase preservation of biological identity. The correspondence mapping was then visualized using PHATE114.

Evaluation of integration
To assess integration of batches, we assessed the tradeoff of integration and biological preservation. First, we computed the modified average silhouette width (ASW) of batch, as previously implemented116. This metric measures the silhouette of a given batch for each cell type (where we have only one cell type). For the scaled metric we used here, 0 indicates separation of batches, and 1 indicates perfect overlap of batches. To ensure we did not enforce mixing where cells are transcriptionally dissimilar, we also evaluated our integration based on preservation of biological signatures. We thus computed the cosine similarity of each cell before and after mapping (where high cosine similarity is 1 and low cosine similarity is -1). We termed this measure global distortion based on a related metric using cell correlation119. A good mapping has both batch integration and cosine similarity near 1, though generally, a tradeoff is necessary due to large technical variation between datasets. In particular, we expected some transformation of cellular gene expression to ensure mapping. Thus, we aimed for mild to moderate global distortion (cosine similarity = 0.2-0.6) and high batch integration (integration > 0.6) in order to ensure the nearest cells in the obesity progression were truly near in the embedding space.

For vehicle to WT mapping, scMMGAN balanced batch integration (mean Batch ASW=0.91) and biological preservation (mean cosine similarity before versus after correction=0.69; Supplementary Data 4a). For mSTZ to Lepob/ob, scMMGAN also observed high integration performance (mean Batch ASW=0.64, mean cosine similarity=0.71; Supplementary Data 4b). For the mouse integrated atlas, after scMMGAN alignment, the datasets showed high batch integration (Batch ASW=0.98) and mild global distortion across all cells (cosine similarity=0.58 before and after alignment), with no distortion specific to a particular archetype, developmental stage, disease state, or dataset (Supplementary Data 4c). For human samples, scMMGAN showed high batch integration (Batch ASW=0.97) and moderate global distortion across all cells (cosine similarity=0.35 before and after alignment), common when mapping between datasets with substantial differences, including species-specific differences, with no distortion specific to a particular archetype, sex, or disease state (Supplementary Data 4d).

AAnet
Overview
Defining cell states via cellular clustering remains a useful tool for single-cell analysis but proves insufficient for characterizing heterogeneity within cell types. Indeed, β cell heterogeneity is well appreciated in the field, but cell states are highly plastic and require an alternative approach to model this plasticity. AAnet24,25 is a latent-space archetypal analysis approach that learns cellular “archetypes”, or extreme states, and defines each cell with respect to the extreme states. This allows us to identify cells highly specific to a given state.

Archetypal commitment
Archetypal commitment is captured in the AAnet latent space, where each cell is defined by its affinity to each archetype (a vector of size number of archetypes, where the entries sum to 1 and higher value implies higher commitment). Thus, we determine cells that have commitment over 0.5 for a given archetype as “committed” to that archetype25.

Archetypal assignment for mapped cells
For the vehicle and mSTZ analysis37, we learned mappings directly onto the embedding space used for training the WT AAnet and Lepob/ob AAnet models, respectively. Thus, for these datasets, we used the trained AAnet models to encode archetypal assignments for cells. For other mappings, including with TrajectoryNet and with scMMGAN on the full trajectory of cells, we assigned each cell the archetype of its closest neighbor in the reference dataset.

Experimental details
We ran AAnet (v0.0) with default parameters on each condition embedding separately, defined by running MAGIC followed by PCA (again to enable denoising and invertibility to the gene space). To determine the number of archetypes for each condition24, we first generated a training and test set, where the training test was generated by density subsampling 80% of the dataset, and the test set was the held-out group. We then trained AAnet on the training set 3 times for each number of archetypes between 2 and 10. For each trained model, we calculated the mean squared error (MSE) of the input versus reconstructed cells for the held-out test set. We plotted the MSE and chose the elbow point as the number of archetypes (Supplementary Fig. 6a).

DiffusionEMD
Overview
Patient manifolds are multi-scale manifolds, where one level characterizes relationships between cells, and a higher level characterizes relationships between patients120. To compare single-cell datasets between patients, such as to compare cell type proportions, gene trends, or clinical variables across patient mappings, Diffusion Earth Mover’s Distance (DiffusionEMD)27 models each single-cell data set as a distribution supported on a common cell-cell graph enabling comparison via optimal transport.

Experimental details
We ran DiffusionEMD (v0.5.0) with default parameters to embed human patient samples after their mapping to the obesity progression with scMMGAN.

Quantitative and statistical analyses of scRNA-seq data
Enrichment score calculation
To compute the enrichment of a given label in a mapped population of cells, we annotated each mapped cell using the label of the nearest “real” cell in the reference dataset. We then calculated the proportion of each class in the mapped population. This alone does not sufficiently characterize the mapping of a population due to large class imbalances between labels. Therefore, we calculated the ratio between the mapped label proportions and the label proportions in the reference dataset. We term this the “enrichment score”. For each class, if the ratio is greater than 1, the class is “enriched” in the mapped population. If the ratio is less than 1, the class is enriched in the reference population. If the ratio is equal to 1, the class is equally represented in the reference and the mapped population. For all endocrine cells (Fig. 3), we calculated the enrichment of each cell type for cells-of-origin to Cck+ Lepob/ob β cells and for cells-of-origin to Cck+ Lepob/ob polyhormonal cells versus the distribution of cell types in the WT population. Within β cells (Fig. 3), we calculated the enrichment of each WT archetype for cells-of-origin versus the distribution of archetypes in the WT population. We also calculated the enrichment of each archetype for each timepoint versus the distribution of all archetypes to learn the trajectory of archetypes from WT AT 4 to Lepob/ob AT 2. For the mSTZ dataset (Fig. 5), we calculated the enrichment of each WT archetype for the β2 population versus the distribution of archetypes in the reference WT population. We also calculated the enrichment of each Lepob/ob archetype for the Cck-hi population versus the distribution of Lepob/ob archetypes in the reference. Finally, we calculated the enrichment of β1 versus β2 for mapped cells-of-origin from TrajectoryNet to Lepob/ob AT 2-annotated cells.

Differential expression
For differential gene expression of archetypes, we compared cells committed to each archetype versus all other cells and performed a two-sided Wilcoxon rank sum test, where differentially expressed genes have Benjamini-Hochberg (BH)-adjusted p-value (q-value) < 0.05 and log fold-change > 0.0.

Gene clusters
We clustered the 2,278 highly variable genes based on TrajectoryNet-inferred transcriptional changes from all WT to Lepob/ob AT 2. We computed the mean expression over all trajectories from two TrajectoryNet runs and clustered each gene’s association with continuous time. Genes were clustered using k means clustering with the number of clusters set to two.

Gene set enrichment analysis
Gene set enrichment analysis was performed using GSEApy (v1.1.2)56 and Enrichr57, using the BioPlanet 2019, Panther 2016, KEGG 2019, Reactome 2022, and GO Biological Process 2023 databases. Gene sets with BH-adjusted p-value (q-value) < 0.05 were retained. The hypergeometric test was used for determining statistical significance of overlapping gene sets.

Statistics and Reproducibility
Statistical tests used for each experiment including n, definitions of center, and precision measures are described in the figure legends. A p-value threshold 0.05 was used to denote statistical significance. No statistical method was used to predetermine sample size. In single-cell analyses, outlier populations of β cells with large differences in the number of genes expressed per cell, indicating low quality of the cell populations, were filtered. No other data were excluded. For allocation into experimental treatment groups, mice were randomized following pre-stratification for sex and weight. Investigators were not blinded to allocation but were blinded to outcome assessment (e.g., disease burden evaluation).

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "This paper explores the role of pancreatic endocrine-exocrine crosstalk in tumorigenesis using a combination of experimental and sequencing methods. Key methods include β cell ablation experiments in mice to assess its impact on PDAC development, and the induction of obesity to study its effects on β cell hormone dysregulation and tumor formation. Single-cell RNA-sequencing was used to analyze transcriptional changes in peri-islet exocrine cells, coupled with in silico latent-space archetypal and trajectory analysis to map cellular states and transitions. Additionally, genetic lineage tracing in vivo was used to investigate β cell adaptation and pro-tumorigenic transformation, identifying the JNK/cJun signaling pathway as a putative mediator of these effects.\nMajor weaknesses:\nOverall, the descriptions of the experiments in the paper are well written, but it is difficult to see the conceptual insights this paper contributes. Almost all panels are based on single-cell transcriptional data without functional or mechanistic follow-up. The conclusions are based on a few experiments measuring cellular expression without clear implications or next steps. The discussion does not address weaknesses or alternative interpretations, but rather lists the findings and references consistent with their own work.\nSpecific comments:\nFig 1d: The KPC-Akita group does not show any significant blood glucose levels changes relative the WT group which seems incorrect if they have no insulin-producing islets. The Akita group is hitting the roof at 600 mg/dl glucose and therefore this experiment is not reliable. This should be done under fasted conditions.\nAll heatmaps lack units in the scale bars. It is not possible for the reader to know the fold change based on the presented data in e.g. Fig 3f, 3d, 3h and many other figures. This is a major concern also because the variability is not shown either, and the values are based on a very low cell number extracted from the single-cell RNA seq analysis.\nThe impact of the proposed mechanism of cJun activation under beta cell stress leading to CCK expression is unclear without follow up experiments.\nWhile a lot of transcriptional associations are shown, there are often no or little explanations of how PDAC develops or progresses based on these data. Overall, these concerns dampens my enthusiasm for this study.", "reviewer_2": "Summary\nIn this manuscript, the authors use genetic and computational models to demonstrate that obesity drives a pro-tumorigenic beta cell state, and that endocrine-exocrine crosstalk is important for the development of pancreatic ductal adenocarcinoma. Using the Akita mouse line crossed into the KPC PDAC model, the authors establish a tumor-supporting role for endocrine beta cells at baseline. They then leverage deep bioinformatic analysis to evaluate beta cell heterogeneity in obesity and show that pro-tumorigenic CCK+ beta cells are enriched in obesity and arise from the proliferation of immature Ngn3+ beta cells. Importantly, the authors use lineage tracing in vivo to support these findings, demonstrating that Ngn3-CreERT is sufficient to label CCK+ cells that are enriched over time in the context of obesity. Mechanistically, obesity was associated with beta cell transcriptional signatures associated with cellular and ER stress, and pharmacological induction of beta cell stress phenocopied obesity-driven transcriptional changes and increased CCK expression. CCK expression was regulated by cJun N-terminal kinases (JNK) in obese mice in vivo through direct regulation of a CCK enhancer by cJun. Together, the central implications of this work are that (1) beta cells are important for pancreatic tumor development at baseline and (2) that obesity triggers a general beta cell stress response via cJun/JNK, leading to proliferation of immature pro-tumorigenic CCK+ beta cells. These findings are each important in their own right, and are supported by a wealth of computational analysis backed by biological validation. However, one weakness of the manuscript is the failure to unify these two findings. We can infer that the mechanisms of obesity-driven beta cell remodeling identified here may be important for obesity-driven PDAC development, but this is not directly shown. Instead, the deeper investigation of beta cell biology in obesity is not extended to the context of PDAC or obesity-associated PDAC. The authors have the opportunity to strengthen this manuscript by drawing some more direct throughlines between these otherwise compelling and interesting findings. Further, the finding that tumorigenesis is reduced in KPC-Akita mice is not supported by in-depth characterization of the model; building on this work could also greatly strengthen the manuscript and the authors’ claims.\nMajor comments\n1. The finding that beta cell loss in the KPC-Akita model impairs tumor development (Fig.1c) is central to the manuscript and this claim would be better supported by additional characterization of the model. Is tumor mass or volume reduced in KPC-Akita mice? How is survival impacted? Cellular markers such as proliferation, cell death? Metastasis? Assessment of stromal content would also be useful. As it stands, tumor area is only assessed at one time point and it is difficult to parse whether these effects are driven by impacts on tumor latency or proliferation. Although assessing each of these points may not be necessary, any additional details would help the reader better assess the magnitude of the effect of beta cell loss on tumor trajectory in the commonly-used KPC model, and help the reader assess whether beta cells are really “essential” for PDAC development.\n2. In Fig. 2j the authors quantify PanINs in proximity to islets in KCO mice. The text should be modified to reflect this; referring to PanIN as “early tumors” is somewhat misleading. This distinction matters in that low-grade PanIN will not all necessarily progress to PDAC and low-grade PanIN-like bystander lesions have even been detected in normal tissue adjacent to PDAC in mice. The increase in proximal PanIN in KCO mice is suggestive of beta cell promotion of early lesions but not necessarily the development of tumors. The authors could quantify tumors proximal to islets in KCO mice at later timepoints to support this point. Further, it seems possible that if obesity indeed drives up beta cell mass and PanIN content as would be expected, then the increase in peri-islet lesions could simply be driven by a greater density of islets or PanIN within the tissue in obese mice. It appears that the authors evaluated KC mice at 13 weeks and KCO at 7 weeks, possibly to account for this though it’s not clear if this is why these timepoints were chosen. The authors should include quantification of islets/beta cells or PanIN to illustrate the magnitude of the difference in these two conditions so the reader could assess the possible relevance of islet/PanIN density. As an additional control, the authors could consider quantifying PanIN or tumors proximal to islets in KPC vs KPC-Akita. This study would have the same issue with different density of islet/PanIN by condition but would provide further support that beta cell function is (or is not) required for peri-islet PanIN or tumor development. Assessment of islet-proximal lesions in KCO mice +/- a CCK inhibitor or other beta cell signaling modulator would probably be most appropriate to truly answer this question but is certainly not necessary within the scope of this manuscript.\n3. The authors set out to study both endocrine (beta cell) contributions to pancreatic tumorigenesis, and the underlying mechanisms driving pro-tumorigenic beta cell function in the context of obesity. However, the authors do not directly link the identified mechanisms of obesity-associated beta cell proliferation/phenotype to PDAC or obesity-associated PDAC. It would strengthen the manuscript to more directly demonstrate whether beta cell stress/proliferation and JNK-driven CCK expression contribute to baseline PDAC and obesity-associated PDAC. To extend these findings to the cancer context, the authors could:\n• Show whether PDAC induces a similar beta cell stress response in vivo at baseline by comparing and quantifying beta cell number or insulin/CCK expression in wild-type, wild-type Akita, KPC and KPC-Akita mice, and an obesity KPC or KCO model. How does baseline PDAC development impact beta cells? It’s unclear whether the obesity-driven changes to beta cell biology identified here are similarly present and important for baseline PDAC development. If changes to beta cell state/mass identified in obesity (for example increase in CCK+ beta cells) are not identified in KPC mice, what alternative mechanism explains the role for beta cells in baseline tumor development?\n• The authors demonstrated an important role for beta cell CCK on obesity-associated tumorigenesis in prior work. However, additional studies in KPC/KPC-Akita mice would validate the role of CCK+ beta cells in baseline PDAC development. The authors could test if CCKR inhibitor or JNKi reduce tumor burden in KPC mice to verify the functional relevance of the new mechanism in vivo and the role of CCK+ beta cells outside of the obesity setting. KPC-Akita mice could be used to control for beta cell-independent effects of inhibitors on PDAC tumors.\n• Evaluate the effect of obesity on tumorigenesis in the context of beta cell depletion (KPC-Akita). This experiment would delineate the extent to which obesity-driven cancer is dependent on beta cells. Is there a reason the authors did not pursue these studies, or why the KPC-Akita model might not be suitable to answer this question?\nMinor comments\n1. Could the authors please provide any quantification of beta cell depletion in the KPC-Akita model (Figure 1b)? It would be helpful to know the extent to which insulin+ beta cells are depleted.\n2. The authors should define “PP cells” in the section on “Dynamic changes in beta cell expression of insulin and CCK with obesity”.\n3. Reg expression was reduced by CCKRi treatment of Ob/ob mice; did this correspond to a loss of the Cluster 3 Reg+ cells? (Fig. 2i)\n4. “Thus, Cck+ cells arise from pre-existing beta cells….” In statements regarding the use of bioinformatic tools such as TrajectoryNet it might be more appropriate to soften the language to indicate that these tools really generate a hypothesis as to the lineage or cell of origin instead of truly indicating which population arose from another as illustrated instead by the authors’ later lineage tracing studies.", "reviewer_3": "The article by Garcia et al. crosses the KPC model of PDAC with the Akita model of beta-cell ablation and shows decreased tumour formation at 12 weeks in the KPC-AKITA model even in the presence of elevated blood glucose, suggested that more available free glucose in the bloodstream isn’t driving tumorigenesis but the beta-cells directly. Following this discovery, they switch focus to beta-cell transformation during obesity, which is a risk factor for PDAC development, and find that acinar cells in the peri-islet space have altered phenotype during obesity. They then characterise beta- heterogeneity during obesity, identifying an embryonic/immature beta-cell program in mice driven by JNK/cJun which leads to increased CCK expression, a known PDAC promoter.\nAlthough the topic is highly relevant to better understand exocrine and endocrine pancreas interaction/communication there is need for a reinterpretation of the results and in some cases refine some of the stated conclusions:\nMajor Comments\n• The author consistently jump between figures and the order of graphs within individual figures does not match the in-cite text. Re-ordering either the figures or the in-cite text to match each other would make it a lot easier for the reader to follow the data. This was especially frustrating with figures 3 and 4, and also the out-of-order graphs in extended Fig. 7. Furthermore, their reasoning and logic for using the various bioinformatic tools employed in the paper isn’t stated very clearly until the discussion.\n• Authors state, “Despite decreased tumorigenesis, KPC-Akita mice showed markedly elevated glucose levels (Fig. 1d), suggesting that hyperglycemia is not independently pro-tumorigenic in the absence of β cells”. Surprisingly KPC mice have a glucose basal level of 400 mg/dl, these glucose levels are extremely high, could the authors comment on that. 400 vs 600 mg/dl will not allow us to draw any conclusion about the roles of hyperglycemia in PDAC.\n• In extended figure 3. The authors show the single cell data of the 3 mouse models analyzed, but the number of cells per model is different. A proportion plot per model will be helpful in this regard as many conclusions are based on this single cell data.\n• As the Akita model causes inflammatory cell death of beta-cells, can the authors be confident that there isn’t an artificial recruitment of immune cells to the pancreas, which alters the tumour microenvironment of the pancreas, therefore priming it to remove early-stage tumours, and therefore the results they achieve is immune dependent and not beta-cell dependent?\n• What would further strengthen the authors argument that glucose isn’t driving PDAC tumour development is a time course of glucose levels over the disease progression of the two models KPC and KPC-AKITA. Moreover, further investigation of the role of insulin signalling pathway will also help elucidate the main driver of this conclusion: immature beta cell phenotype is activated by ER stress, glucose changes or insulin ones.\n• The study has focused on acinar derived PDAC, but they frequently state PDAC as a whole, ignoring ductal derived PDAC. A more nuance discussion specifying acinar vs ductal origin is needed throughout.\n• In the Peri-islet acinar cell adaptation section drawing any conclusion based on 260 acinar cells is an overstatement. Moreover, assuming that the acinar cells come from peri-islet region exclusively in not right as the islet preps contain also acinar contaminants not located in the peri-islet region. A detailed validation of this section is needed at list at the level of IF or similar. Tosto et al. already described 3 acinar populations in the human pancreas. The authors should look at Tosti or other scRNAseq data of acinar cells as well.\n• In the section of “identification of the cell-of-origin for pro-tumorigenic β cells in obesity”, the authors introduce the possibility of several cells of origin of new beta cells but then just analyze beta cell as possible pro-tumorigenic β of origin. Further elaborating this decision is important.\n• In extended Fig 4a, the authors have quantified the amount of TdTom+ cells that are glucagon positive, the same should be performed for the Ins1-CreERT and Ngn3-CreERT models in Supp Fig 4b, given the incomplete labelling hypothesis presented in the paper. I also think an explanation for the inefficient labelling should be stated in-text or hypothesised on. The interpretation of Ngn3-CreERT labelling results are unclear. Moreover, as far as it is known, the number of Neurog3 positive cells after weaning is extremally low, thus not sure about the labelling, are those cells in the islets or ducts? Several papers argue about it, a detailed characterization of the model and % of positive cells is needed.\n• The reasoning for choosing Lepob/ob AT2 to examine in TrajectoryNet out of the different ATs in group 2 is not clear. By choosing that as the endpoint, are you then forcing its evolutionary trajectory? Explanation is needed as to why the authors have chosen this comparison over other archetypes.\n• Extended Fig. 6 has not been expanded upon in the text. Please explain these results in more detail.\n• In Fig.6b, WT AT 4 is only represented in 3/15 ND patients. This is inconsistent with the mouse data presented earlier, and questions whether your proposed cell of origin is conserved in humans.\n• Does treatment with the JNKi inhibitor result in decreased tumour formation in obese PDAC models? While the embryonic program stimulated in B-cells during pharmacological (STZ) stress and obesity is interesting, the data isn’t related back to the PDAC model in the first figure.\nMinor Comments\n• Overall, the authors concluding statements at the end of every results section are very vague, and a more informative summary of each section's results would make it easier for the reader to interrupt and follow their conclusions.\n• In Fig.1b, the number of mice used to ascertain that insulin was decreased and ER stress induced in the KPC-Akita model relative to KPC is not stated in the figure legends or methods. Quantification of the insulin and BiP positive area within the islets would help support the conclusions the authors have made in the manuscript “We confirmed decreased β cell mass and induction of ER stress in islets of KPC-Akita mice “and therefore their conclusions that B-cell loss impedes PDAC tumorigenesis.\n• In Figure 1C, it is not stated in methods how the authors quantified tumour burden in the KPC and KPC-Akita models.\n• The author writes, “Having established that β cells are basally required for exocrine tumorigenesis, …”. This is a very strong conclusion, as the KPC-Akita mice still have 40% tumour burden. I think promote would be a more accurate interpretation of the author’s results.\n• In methods, it is stated that only male mice were used for STZ and scRNAseq experiments due to the reduced beta-cell effects of STD and a High Fed Diet. Please comment on whether these results are then relevant for human female patients.\n• In Fig.3D and the accompanying text, as you heavily discuss the biological functions of the groups of archetypes, I would rearrange the heatmap on Fig. 3D to follow the evolutionary tree in ext. Fig.5 and then colour code them by WT, HFD and Lepob/ob to make it easier for the reader to follow your conclusions.\n• In Fig.5J, you only show the mapping of WT AT 4 and OB/OB AT 2. What about the other archetypes, these should be included in the supplementary/extended figures.", "reviewer_4": "The paper presents an interesting study with extensive analyses. However, the current structure makes it somewhat confusing. The experiments are not clearly presented in a hypothesis-driven manner, which makes it difficult to follow the logical flow of the study. Additionally, the bioinformatics methodology needs more detail, particularly regarding the parameters used and how they were selected.\nMajor Concerns\nTumorigenesis in KPC-Akita Mice (Figure 1)\nSince tumors are still observed in KPC-Akita mice, can the authors definitively conclude that β-cell loss impedes pancreatic tumorigenesis? This claim should be carefully re-evaluated.\nSingle-Cell Data Generation and Experimental Design\nThe comparison between congenic obesity models (HFD-fed for 10 weeks, Lepob/ob) and age-matched WT controls raises concerns. The authors should clarify how they accounted for age-related differences to isolate the effects of the conditions studied.\nClustering and Data Quality (Figure 2 Extended Data)\nCluster 6 appears to express multiple markers from different cell types. How was cell-type annotation performed to ensure accuracy?\nWhat were the nFeature thresholds and percentage of mitochondrial genes for clusters 13 and 14, which were considered low-quality?\nMethods: Data Integration and Batch Effect Control\nHow was data integration performed, and how were batch effects controlled?\nWhich covariates were used to integrate the datasets and minimize technical noise?\nHow was n = 20 clusters determined?\nWhat methods were used to verify that the clusters retained biological identity and were not artifacts of overclustering?\nTrajectory Analysis Validity\nGiven the differences in age and conditions among the animals, could the trajectory analysis be artificially linking cell types that may not truly represent progenitor or intermediate states?\nThe authors hypothesize that acinar cells surround the islets (peri-islet localization). Can they compare their gene signature with spatial transcriptomics data from the literature to validate this claim?\nGene Expression Measurements (Figure 2)\nwhat metric is used for gene expression measurement?\nHow are normalized counts defined in this analysis?\nValidation of β-Cell Archetypes\nThe authors should validate β-cell archetypes using IHC or FACS.\nHow was the number of subpopulations determined?\nIn WT mice, how do the β-cell archetypes vary with age?\nBatch Effect Concerns and Best Practices\nThe authors state that datasets were not directly comparable due to batch effects. However, best practices in single-cell data integration should be followed to improve comparability: https://www.sc-best-practices.org/preamble.html.\nMethodological Reproducibility: Software, Databases, and Parameters\nThe methods section should explicitly state:\nSoftware versions used (e.g., Seurat, Scanpy, Cell Ranger).\nReference databases used for annotation.\nFlags and parameters for each function/tool and justification for their selection.", "reviewer_5": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.", "reviewer_6": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts."}, "author_response": "We thank the reviewer for the comments. In the revised manuscript, we have provided additional experiments to establish the functional importance of β cell CCK expression in obesity-driven cancer. Specifically, we have used a new mouse model that enables conditional knockout of CCK (CCKFlox/Flox) in the pancreas of KCO mice (Fig. 2). We find that pancreas-specific knockout of CCK reduces tumor formation (Fig. 2) without reducing body weight, decreasing insulin secretion, or compromising glucose homeostasis (Extended Data Fig. 2). Since our scRNA-seq and IF analyses do not identify significant CCK expression in any cells in the pancreas beyond β cells, our data argue that β cell-derived CCK is a bona fide driver of obesity-driven pancreatic cancer. We also confirmed that β cell CCK overexpression is sufficient to promote pancreatic tumorigenesis in lean KC;Lepob/+ mice to a degree that was non-significantly different to obese KCO mice (Fig. 1). Finally, our results establish that CCK – rather than insulin – more strongly correlates with disease burden across our models (Fig. 1l,m and Fig. 2g,h), challenging the conceptual dogma in the field regarding insulin as the primary β cell driver of obesity-associated PDAC. Together, these data implicate targeting β cell CCK as potential means to intercept obesity-driven PDAC development.\n\n[... Author response continues with detailed point-by-point responses to all reviewers ...]"}, {"round": 2, "reviewer_comments": {"reviewer_1": "The manuscript has been significantly improved by incorporating text changes and adding additional experiments. I have no further scientific comments.\nThe code has not been verified by any reviewer. This should be corrected.", "reviewer_2": "In this manuscript, the authors use genetic and computational models to demonstrate that obesity drives a pro-tumorigenic beta cell state, and that endocrine-exocrine crosstalk via CCK is important for obesity-associated pancreatic ductal adenocarcinoma. The authors use bioinformatic analysis of scRNA-seq followed with lineage tracing in vivo to show that obesity drives cellular stress via cJUN/JNK, leading to proliferation of immature (Ngn3+) CCK-expressing beta cells that could drive tumorigenesis. Through additional model characterization and modification/reorganization of the text the authors have sufficiently addressed our concerns. The authors have re-organized assessment of the KPC-Akita model which improves the clarity of the manuscript, highlighting the role of beta cells in tumor development in lean mice, and thus independent of obesity and CCK. Inclusion of a new CCKfl/fl allele crossed to the KCO mouse model provides compelling evidence that beta cell derived CCK is a driver of obesity-associated tumor development. The authors should be congratulated on this nice study that will expand our understanding regarding the relationship between obesity and PDAC.\nMy minor suggestions for improvements, are:\nThe authors should briefly note the timepoint at which each model was assessed for disease burden either in the text or legends. If CCKfl/fl KCO mice were assessed at the same timepoints as KC counterparts it would be helpful to present the disease burden in KCO +/- CCK side by side with KC counterparts to demonstrate the magnitude of the effect relative to lean mice.\nAdditionally, if the authors assessed ADM, PanIN, and PDAC across mouse models it would be more informative to report the burden of each rather than combined disease burden.\nIn Figure 8f, it is somewhat confusing that CD45 and SMA+ cells are quantified within PanIN, PDAC lesions; I’m not sure what that would mean. Were these cells quantified with FOV containing PanIN or PDAC lesions?", "reviewer_3": "The reviewed manuscript submitted by Garcia et al. has addressed most of our previous concerns and figures have improved. The results support much better the conclusions of the manuscript.\nMost strikingly, the addition of the CCK overexpression and CCK knockout models has really strengthened the conclusions of this manuscript and shows a much clearer role for CCK on PDAC progression then the first iteration of this manuscript. The serial glucose measurements between the two models strengthens the conclusions proposed in the original manuscript, that hyperglycaemia is not the main driver of PDAC progression in these models. And the correlation analyses showing CCK expression associates with disease burden instead of insulin strengthens these conclusions.\nThe revised manuscript is clearer, the figures are easier to interrupt, and the story has a much better flow. They have corrected many of the statements that we felt were overstated or not clear in the original manuscript. The enrichment of the Acinar-REG+ signature identified by Tosti et al from the human pancreas in cluster 3 better supports that the acinar cells are of peri-islet origin. They have addressed the question of sex bias (where they previously only used male mice for the STZ and scRNAseq experiment) by performing all their new in-vivo experiments with both sexes and showing no significant differences between the sexes, suggesting the vast majority of their results will still be relevant for female patients.\nOverall, we think the authors took great efforts to address many of the reviewer’s previous concerns, and the manuscript is a lot stronger for it. I endorse the acceptance of this paper.", "reviewer_4": "Overall, the revisions are appreciated. The revised organization and the new data linking β-cell CCK expression to obesity-associated PDAC progression improve the readability of the manuscript.\nPlease add minimal clarifications and diagnostics to strengthen rigor and reproducibility: specify the feature space and cell subset used for EMD. The added EMD analysis is a useful global similarity check, but it is not sufficient to support a claim of “lack of batch effect”: sample-level EMD can be driven by differences in cell-state composition and QC/capture differences, may miss subtle batch/age shifts within β cells (or β substates), and does not directly validate the absence of batch effects in the same latent space used for downstream clustering/trajectory unless that is explicitly shown. Therefore, please temper the statement “verifying lack of batch effect” and include one within–β cell batch/age mixing diagnostic on the trajectory input/latent space (e.g., ASW_batch, iLISI, or kBET) plus a β-cell embedding colored by sample, and show the distribution of trajectory position/pseudotime by sample to demonstrate the inferred trajectory is not simply tracking sample identity. Finally, it may be more appropriate to soften trajectory/archetype language to reflect that these tools generate testable hypotheses regarding lineage relationships or candidate cells of origin, whereas the definitive evidence for lineage directionality is provided by the subsequent in vivo lineage tracing experiments.", "reviewer_5": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.", "reviewer_6": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts."}, "author_response": "We thank the reviewer for their positive feedback on our work. [...] We thank the reviewer for their acknowledgment of the improvement of the manuscript and overall enthusiasm for our work. [...] We thank the reviewer for the acknowledgment of the improvement of the manuscript."}]

Decision letter