Advertisement

Single-Cell Transcriptomics Reveals a Conserved Metaplasia Program in Pancreatic Injury

Published:October 22, 2021DOI:https://doi.org/10.1053/j.gastro.2021.10.027

      Background & Aims

      Acinar to ductal metaplasia (ADM) occurs in the pancreas in response to tissue injury and is a potential precursor for adenocarcinoma. The goal of these studies was to define the populations arising from ADM, the associated transcriptional changes, and markers of disease progression.

      Methods

      Acinar cells were lineage-traced with enhanced yellow fluorescent protein (EYFP) to follow their fate post-injury. Transcripts of more than 13,000 EYFP+ cells were determined using single-cell RNA sequencing (scRNA-seq). Developmental trajectories were generated. Data were compared with gastric metaplasia, KrasG12D-induced neoplasia, and human pancreatitis. Results were confirmed by immunostaining and electron microscopy. KrasG12D was expressed in injury-induced ADM using several inducible Cre drivers. Surgical specimens of chronic pancreatitis from 15 patients were evaluated by immunostaining.

      Results

      scRNA-seq of ADM revealed emergence of a mucin/ductal population resembling gastric pyloric metaplasia. Lineage trajectories suggest that some pyloric metaplasia cells can generate tuft and enteroendocrine cells (EECs). Comparison with KrasG12D-induced ADM identifies populations associated with disease progression. Activation of KrasG12D expression in HNF1B+ or POU2F3+ ADM populations leads to neoplastic transformation and formation of MUC5AC+ gastric-pit-like cells. Human pancreatitis samples also harbor pyloric metaplasia with a similar transcriptional phenotype.

      Conclusions

      Under conditions of chronic injury, acinar cells undergo a pyloric-type metaplasia to mucinous progenitor-like populations, which seed disparate tuft cell and EEC lineages. ADM-derived EEC subtypes are diverse. KrasG12D expression is sufficient to drive neoplasia when targeted to injury-induced ADM populations and offers an alternative origin for tumorigenesis. This program is conserved in human pancreatitis, providing insight into early events in pancreas diseases.

      Graphical abstract

      Keywords

      Abbreviations used in this paper:

      ADM (acinar to ductal metaplasia), Co-IF (co-immunofluorescence), DEG (differentially expressed gene), EEC (enteroendocrine cell), EYFP (enhanced yellow fluorescent protein), IHC (immunohistochemistry), PanIN (pancreatic intraepithelial neoplasia), PDAC (pancreatic ductal adenocarcinoma), scRNA-seq (single-cell RNA sequencing), SEM (scanning electron microscopy), sNuc-seq (single nucleus RNA sequencing), SPEM (spasmolytic polypeptide-expressing metaplasia), TF (transcription factor)

       Background and Context

      Acinar to ductal metaplasia is a reparative process in the pancreas in response to injury and is a precursor lesion for pancreatic ductal adenocarcinoma.

       New Findings

      In response to chronic injury, acinar cells undergo pyloric metaplasia, transitioning into progenitor-like populations, tuft cells, and enteroendocrine cells, generating a gastric pyloris-like phenotype. Oncogenic Kras targeted to acinar to ductal metaplasia populations results in neoplasia and MUC5AC expression.

       Limitations

      Functional studies are required to identify the role of acinar to ductal metaplasia populations in pancreatic injury.

       Impact

      This study identifies significant acinar cell plasticity, the formation of metaplastic epithelial populations with the potential to fuel disease progression, and acinar to ductal metaplasia populations as a source of neoplasia.
      Damage to the exocrine compartment of the pancreas results in infiltration of immune cells, stromal deposition, loss of acinar cells, and the emergence of a metaplastic cell lineage pathologically defined by the aberrant appearance of ducts termed “acinar to ductal metaplasia” (ADM).
      • Giroux V.
      • Rustgi A.K.
      Metaplasia: tissue injury adaptation and a precursor to the dysplasia-cancer sequence.
      ADM is a reparative program in which digestive enzyme-producing acinar cells transdifferentiate to cells with ductal organization as part of a critical acute phase program (paligenosis) that enables tissue reconstruction following injury.
      • Storz P.
      Acinar cell plasticity and development of pancreatic ductal adenocarcinoma.
      ,
      • Willet S.G.
      • Lewis M.A.
      • Miao Z.F.
      • et al.
      Regenerative proliferation of differentiated cells by mTORC1-dependent paligenosis.
      Under conditions of sustained injury, inherited mutation, environmental factors, or idiopathic stimuli, patients can develop chronic pancreatitis, a known risk factor for pancreatic ductal adenocarcinoma (PDAC). ADM is a pathognomonic epithelial change during chronic pancreatitis and occurs in response to oncogenic mutation in the exocrine tissue, identifying it as a pre-neoplastic event thought to be critical for eventual malignant transformation.
      • Giroux V.
      • Rustgi A.K.
      Metaplasia: tissue injury adaptation and a precursor to the dysplasia-cancer sequence.
      ,
      • Murtaugh L.C.
      • Keefe M.D.
      Regeneration and repair of the exocrine pancreas.
      The induction of injury and inflammation-induced ADM, therefore, may be one factor by which chronic pancreatitis increases the risk of PDAC.
      We have recently shown that ADM results in the transdifferentiation of acinar cells to tuft cells, which are rare chemosensory cells typically absent from the normal mouse pancreas.
      • DelGiorno K.E.
      • Naeem R.F.
      • Fang L.
      • et al.
      Tuft cell formation reflects epithelial plasticity in pancreatic injury: implications for modeling human pancreatitis.
      Using a combination of low-input RNA-sequencing strategies, super-resolution microscopy, and genetically engineered mouse models, we demonstrated that tuft cells inhibit pancreatic tumorigenesis by modulating tissue stroma.
      • DelGiorno K.E.
      • Naeem R.F.
      • Fang L.
      • et al.
      Tuft cell formation reflects epithelial plasticity in pancreatic injury: implications for modeling human pancreatitis.
      • Delgiorno K.E.
      • Hall J.C.
      • Takeuchi K.K.
      • et al.
      Identification and manipulation of biliary metaplasia in pancreatic tumors.
      • DelGiorno K.E.
      • Chung C.Y.
      • Vavinskaya V.
      • et al.
      Tuft cells inhibit pancreatic tumorigenesis in mice by producing prostaglandin D2.
      These data demonstrate that rare ADM-derived cell types can significantly impact disease progression. Given the role of ADM in disease initiation and the potentially critical link between injury and tumorigenesis, we combined lineage tracing and single-cell RNA sequencing (scRNA-seq) to identify populations arising from ADM.
      Using this strategy, we identified the emergence of tuft cells, enteroendocrine cells, and mucin/ductal cell populations bearing gene signatures similar to spasmolytic polypeptide-expressing metaplasia (SPEM) in the stomach (Tff2+Muc6+Gif+). Gene expression analysis of these populations, in combination with trajectory inference and RNA velocity modeling, suggests varying proclivity for seeding disparate differentiated tuft and enteroendocrine cell lineages. Formation of distinct cell types was confirmed by marker immunostaining and electron microscopy for cell type–specific structures. Our data were compared with KrasG12D-induced ADM, and KrasG12D expression was activated in injury-induced ADM populations using several inducible Cre drivers. Finally, we compared our dataset with published single nucleus RNA sequencing (sNuc-seq) data from human pancreatitis and identified conserved processes.
      • Tosti L.
      • Hang Y.
      • Debnath O.
      • et al.
      Single nucleus and in situ RNA sequencing reveals cell topographies in the human pancreas.
      These studies identify the emergence of previously undescribed acinar-derived cell types in pancreatic injury, metaplasia transcriptional programs shared between organ systems, population changes between injury and tumorigenesis, and a conserved response to exocrine injury in human disease.

      Results

       Lineage Tracing and scRNA-seq Identifies Epithelial Heterogeneity in Injury-induced ADM

      The goal of these studies was to evaluate cellular heterogeneity in ADM resulting from chronic injury. To this end, we conducted lineage tracing by labeling adult pancreatic acinar cells with enhanced yellow fluorescent protein (EYFP) in Ptf1aCreERTM;RosaLSL-EYFP/+ (CERTY) mice. We then induced pancreatic injury with repeated injection of the cholecystokinin ortholog caerulein; resulting in metaplasia, cell death, and inflammation.
      • DelGiorno K.E.
      • Naeem R.F.
      • Fang L.
      • et al.
      Tuft cell formation reflects epithelial plasticity in pancreatic injury: implications for modeling human pancreatitis.
      ,
      • Sah R.P.
      • Garg S.K.
      • Dixit A.K.
      • et al.
      Endoplasmic reticulum stress is chronically activated in chronic pancreatitis.
      Mice were sacrificed after either 2 or 4 weeks of treatment and pancreata were collected (Figure 1A, Supplementary Figure 1A). We isolated EYFP+ cells by fluorescence-activated cell sorting and profiled the transcriptome with 10x Genomics 3ʹ scRNA-seq (Figure 1A, Supplementary Figure 1B). We used lineage tracing to identify metaplastic populations derived from acinar cells (ie, ADM populations), but did not perform additional selections in order to agnostically profile the EYFP+ population. The combined libraries from the 2- and 4-week timepoints (2 mice per timepoint) contained a total of 21,140 cells with an average of 2782 genes and 17,253 unique molecular identifiers detected per cell (Supplementary Table 1). Clusters were annotated by examining both classic cell type markers and previously published single-cell gene signatures in pancreas tissue.
      • Elyada E.
      • Bolisetty M.
      • Laise P.
      • et al.
      Cross-species single-cell analysis of pancreatic ductal adenocarcinoma reveals antigen-presenting cancer-associated fibroblasts.
      Cells clustered mainly by cell type and data from the 2 timepoints were found to largely overlap (Supplementary Figure 1C). We identified both epithelial and stromal populations in our dataset; however, EYFP transcripts were found to be specific to Krt8+ epithelial populations (Supplementary Figure 2A and B).
      Figure thumbnail gr1
      Figure 1Identification of ADM lineages using scRNA-seq and unbiased clustering. (A) Schematic of lineage tracing and scRNA-seq strategy to identify ADM populations. (B) Uniform Manifold Approximation and Projection (UMAP) showing annotated clusters in the EYFP+ population. (C) Immunostaining for EYFP (yellow), acinar marker amylase (cyan), endocrine marker chromogranin A (CHGA, magenta), and DAPI (blue), in normal and injured CERTY pancreata. White arrowheads, YFP+ EECs. I, islets. Scale bars, 100 μm. (D) SEM of the injured pancreas highlighting an acinar cell, (E) a tuft cell, (F) a mucinous cell, and (G and H) 2 EECs. Scale bars, 500 nm; tuft cell inset, 166.7 nm.
      We removed EYFP-negative stromal cells, potential doublets, and further filtered resulting in 13,362 EYFP+ epithelial cells. Cells from the 4 libraries were integrated using the FastMNN approach (see methods). We identified 16 clusters falling into 4 broad categories as well as intermediate populations. Acinar cells (Prss1+Ptf1a+Try4+Cpa1+; 4592 cells) and tuft cells (Pou2f3+Siglecf+Trpm5+Gnat3+; 1280 cells) were identified, consistent with previous studies, validating our methodology (Figure 1B, Supplementary Figure 2D, Supplementary Table 2).
      • DelGiorno K.E.
      • Naeem R.F.
      • Fang L.
      • et al.
      Tuft cell formation reflects epithelial plasticity in pancreatic injury: implications for modeling human pancreatitis.
      ,
      • Schlesinger Y.
      • Yosefov-Levi O.
      • Kolodkin-Gal D.
      • et al.
      Single-cell transcriptomes of pancreatic preinvasive lesions and cancer reveal acinar metaplastic cells' heterogeneity.
      To examine changes in YFP+ acinar cells, we compared our dataset with several published scRNA-seq datasets of normal murine acinar cells and found that injured acinar cells express higher levels of markers related to stress, proliferation, pyloric metaplasia, and ADM (Supplementary Figures 3 and 4, Supplementary Table 3). We also observed a decrease in the percentage of transcripts from digestive enzymes and an increase in the number of active genes detected (Supplementary Figure 4C and D).
      • Schlesinger Y.
      • Yosefov-Levi O.
      • Kolodkin-Gal D.
      • et al.
      Single-cell transcriptomes of pancreatic preinvasive lesions and cancer reveal acinar metaplastic cells' heterogeneity.
      ,
      • Hosein A.N.
      • Huang H.
      • Wang Z.
      • et al.
      Cellular heterogeneity during mouse pancreatic ductal adenocarcinoma progression at single-cell resolution.
      Reads were down-sampled to 3k per cell in all samples to account for the influence of read-depth differences between datasets (Supplementary Figure 4D, methods). In addition to these populations, we identified mucin/ductal cells (Muc6+Pgc+Krt19+Car2+; 3983 cells) and enteroendocrine cells (EECs, Neurog3+Pax6+Chga+ Neurod1+; 1208 cells). Intermediate ADM-derived populations include acinar/mucin-ductal cells (1955 cells) characterized by decreasing levels of acinar marker genes and increasing expression of mucin/ductal genes, as well as a Sox4+ EEC/tuft putative progenitor cell population (205 cells) (Figure 1B, Supplementary Figure 2C and D, Supplementary Table 2).
      We used co-immunofluorescence (Co-IF) for EYFP and cell type markers in normal and injured CERTY pancreata to validate transcriptomically identified populations. In uninjured, tamoxifen-treated CERTY pancreata, EYFP expression is restricted to acinar cells (amylase+) and is absent from ductal cells and islets (CHGA+), consistent with prior reports.
      • Schlesinger Y.
      • Yosefov-Levi O.
      • Kolodkin-Gal D.
      • et al.
      Single-cell transcriptomes of pancreatic preinvasive lesions and cancer reveal acinar metaplastic cells' heterogeneity.
      ,
      • Pan F.C.
      • Bankaitis E.D.
      • Boyer D.
      • et al.
      Spatiotemporal patterns of multipotentiality in Ptf1a-expressing cells during pancreas organogenesis and injury-induced facultative restoration.
      Under conditions of chronic injury, however, EYFP is expressed in both acinar and ductal cells (including tuft cells), as well as a large population of CHGA+ cells consistent with EEC formation (Figure 1C, Supplementary Figure 5).
      • DelGiorno K.E.
      • Naeem R.F.
      • Fang L.
      • et al.
      Tuft cell formation reflects epithelial plasticity in pancreatic injury: implications for modeling human pancreatitis.
      As an orthogonal approach to confirm these cell types, we used scanning electron microscopy (SEM) on ultra-thin tissue sections from the pancreata of caerulein-treated mice. SEM readily identified acinar cells by the presence of large zymogen granules and abundant endoplasmic reticulum. Tuft cells were identified by their distinct microvilli and deep actin rootlets.
      • Sato A.
      Tuft cells.
      We identified mucinous cells by mucin granules and several EEC subtypes, characterized by electron-dense granules (Figure 1DH). Collectively, these analyses reveal previously unrecognized epithelial heterogeneity in injury-induced ADM.

       Developmental Trajectory Analyses Identify Distinct Secretory Cell Lineages

      The cellular heterogeneity identified by scRNA-seq could arise by multiple mechanisms including, but not limited to, (1) direct transdifferentiation from acinar cells, (2) sequential generation from acinar cells to intermediate progenitors, or (3) generation of a common intermediate for tuft and EEC populations, which can then interconvert. We used several computational biology approaches to infer probable lineage trajectories (Figure 2A–D). Trajectory inference algorithms, like Monocle, traverse the collective heterogeneity arising from the asynchronous developmental progression of individual cells.
      • Trapnell C.
      • Cacchiarelli D.
      • Grimsby J.
      • et al.
      The dynamics and regulators of cell fate decisions are revealed by pseudotemporal ordering of single cells.
      Interestingly, the trajectory inferred by Monocle3 indicates that acinar cells first pass through a mucin/ductal state before forming tuft or EEC lineages (Figure 2A and B). A similar trajectory map is captured through the pCreode algorithm, which describes transitional routes in transcriptional variation in principal component space (Figure 2C, Supplementary Figure 6A).
      • Herring C.A.
      • Banerjee A.
      • McKinley E.T.
      • et al.
      Unsupervised trajectory analysis of single-cell RNA-seq and imaging data reveals alternative tuft cell origins in the gut.
      Figure thumbnail gr2
      Figure 2Lineage trajectory analysis of ADM populations. (A) Monocle3 trajectory analysis overlaid on the UMAP of EYFP+ cells from injured CERTY pancreata. (B) Pseudotime projection analysis and (C) pCreode approximation of lineage progression of EYFP+ cells. (D) RNA velocity analysis of all 4 datasets. Arrow direction and length indicate the probable lineage trajectory and velocity, respectively. (E) Brainbow reporter expression in normal and injured pancreata from Ptf1aCreERTM/+;RosaBrainbow.2/+ mice. DAPI, blue. Scale bar, 100 μm. (F) instant Structured Illumination Microscopy (iSIM) analysis of a whole RFP+ clone containing a tuft cell (phalloidin) and EECs (chromogranin A) within the same lesion. Scale bar, 10 μm. (G) Three-dimensional EM reconstruction of ADM showing a heterogeneous consortium of cell types encapsulated by basement membrane with closely apposed vasculature. From left to right, 2 partially annotated cross-sections from the 3-dimensional EM volume (volume: 58 × 29 × 17.5 μm; voxel size: 8–8–70 nm, x-y-z). Volumetric reconstructions of cell plasma membranes and nuclei, as well as a portion of the basement membrane rendered in gray, demonstrating the morphological variability between cell types. Scale bars, 10 μm.
      To independently validate these trajectories, we used RNA velocity, which annotates transitional cell states by leveraging splicing kinetics.
      • La Manno G.
      • Soldatov R.
      • Zeisel A.
      • et al.
      RNA velocity of single cells.
      We calculated RNA velocity individually in all 4 samples; all samples show a consistent vector field from acinar to mucin/ductal, then branching into EEC and tuft cells (Figure 2D). This trajectory is further supported by the activity of key transcription factors (TFs) (Supplementary Figure 6B and C). We applied gene regulatory network, or regulon, analysis using the SCENIC package (Supplementary Table 4). Regulon activity scores reflect the enrichment of a TF together with its correlated direct-target genes.
      • Aibar S.
      • Gonzalez-Blas C.B.
      • Moerman T.
      • et al.
      SCENIC: single-cell regulatory network inference and clustering.
      We observed a switch from acinar-associated TFs (Ptf1a+Rbpjl+Bhlha15/Mist1+) to Mucin (Spdef+Creb3l1+Creb3l4+) and Ductal (Sox9+Hnf1b+Onecut1+) associated TF activation in the Acinar/MucinDuctal intermediate population (Supplementary Figure 6B and C).
      • Gregorieff A.
      • Stange D.E.
      • Kujala P.
      • et al.
      The ets-domain transcription factor Spdef promotes maturation of goblet and paneth cells in the intestinal epithelium.
      • Seymour P.A.
      Sox9: a master regulator of the pancreatic program.
      • Pin C.L.
      • Rukstalis J.M.
      • Johnson C.
      • et al.
      The bHLH transcription factor Mist1 is required to maintain exocrine pancreas cell organization and acinar cell identity.
      Collectively, these data are consistent with a model in which injured acinar cells form a mucin/ductal progenitor population that seeds tuft cells and EECs as distinct lineages.
      We tested this hypothesis by conducting lineage tracing with a Brainbow reporter. Cre recombination in Ptf1aCreERTM;RosaBrainbow.2/+ mice results in the random expression of 1 of 4 reporters (YFP, GFP, RFP, CFP) in acinar cells. Mice were first treated with a very short course of tamoxifen followed by caerulein (as shown in Figure 1A). Analysis of uninjured pancreata shows solitary cells expressing Brainbow reporters. Injured pancreata, in contrast, are characterized by reporter expression in partial or entire lesions, demonstrating clonal expansion from single acinar cells (Figure 2E, Supplementary Figure 6D and E). To examine cellular heterogeneity, we conducted Co-IF on whole clones stained for markers of tuft cells (phalloidin, labels the microvilli and actin rootlets) and EECs (CHGA) and identified both cell types residing within the same clones (Figure 2F). To further evaluate cellular heterogeneity within individual ADM lesions, we performed SEM on 250 serial sections taken of a single metaplastic duct. Segmentation and 3-dimensional reconstruction of the processed image stack revealed an acinar cell adjacent to mucinous cells, a tuft cell, and an EEC with delta cell morphology (Video, Figure 2G).
      • Arrojo E.D.R.
      • Jacob S.
      • Garcia-Prieto C.F.
      • et al.
      Structural basis for delta cell paracrine regulation in pancreatic islets.
      Interestingly, the identified acinar cell is binuclear; how this factors into acinar cell plasticity is currently unknown, but it has been shown that binucleate cells are resistant to cell cycle entry.
      • Wollny D.
      • Zhao S.
      • Everlien I.
      • et al.
      Single-cell analysis uncovers clonal acinar cell heterogeneity in the adult pancreas.

       Injury-induced ADM as a Pyloric-type Metaplasia

      Top differentially expressed genes (DEGs) in the mucin/ductal population include Muc6, Tff2, and Gif, which are characteristic of SPEM in the stomach (Figure 3A, Supplementary Figure 2D, Supplementary Table 2).
      • Schmidt P.H.
      • Lee J.R.
      • Joshi V.
      • et al.
      Identification of a metaplastic cell lineage associated with human gastric adenocarcinoma.
      SPEM is a process similar to ADM in which digestive-enzyme–secreting chief cells undergo metaplasia in response to injury or oncogenic mutation.
      • Nam K.T.
      • Lee H.J.
      • Sousa J.F.
      • et al.
      Mature chief cells are cryptic progenitors for metaplasia in the stomach.
      ,
      • Choi E.
      • Hendley A.M.
      • Bailey J.M.
      • et al.
      Expression of activated Ras in gastric chief cells of mice leads to the full spectrum of metaplastic lineage transitions.
      SPEM occurs within the overall reorganization, known as pyloric metaplasia, characterized by gastric units in the body (main portion) of the stomach reprogramming to resemble the distal (pyloric) region.
      • Nam K.T.
      • Lee H.J.
      • Sousa J.F.
      • et al.
      Mature chief cells are cryptic progenitors for metaplasia in the stomach.
      ,
      • Saenz J.B.
      • Mills J.C.
      Acid and the basis for cellular plasticity and reprogramming in gastric repair and cancer.
      ,
      • Goldenring J.R.
      Pyloric metaplasia, pseudopyloric metaplasia, ulcer-associated cell lineage and spasmolytic polypeptide-expressing metaplasia: reparative lineages in the gastrointestinal mucosa.
      We identified a number of canonical SPEM markers (Cd44, Aqp5, Gkn3, Dmbt1, Wfdc2) enriched in our mucin/ductal population, as well as Gata5, a transcription factor that is widely expressed in gastric/intestinal epithelium but not in normal pancreas epithelium (Figure 3A, Supplementary Figure 7A).
      • Lee S.H.
      • Jang B.
      • Min J.
      • et al.
      Up-regulation of Aquaporin 5 defines spasmolytic polypeptide-expressing metaplasia and progression to incomplete intestinal metaplasia.
      We also identified markers of chief cells (Gif or Cblif, Pga5, Pgc) (Figure 3A, Supplementary Figure 7B). To confirm these findings, we conducted immunohistochemistry (IHC) and/or Co-IF for CD44v9, GKN3, AQP5, WFDC2, and GIF. Consistent with our sequencing data, we identified the expression of SPEM and chief cell markers in injured pancreata, whereas expression is absent from the normal pancreas (Figure 3B, Supplementary Figure 7C). To confirm that this phenotype is ADM-specific, we compared our mucin/ductal population to the EYFP-negative ductal cells in our dataset as well as published normal pancreatic ductal scRNA-seq datasets and found that the injury-induced ADM signature is distinct (Supplementary Figure 3, Supplementary Figure 4E and F, Supplementary Table 3).
      • Schlesinger Y.
      • Yosefov-Levi O.
      • Kolodkin-Gal D.
      • et al.
      Single-cell transcriptomes of pancreatic preinvasive lesions and cancer reveal acinar metaplastic cells' heterogeneity.
      ,
      • Hendley A.M.
      • Rao A.A.
      • Leonhardt L.
      • et al.
      Single-cell transcriptome analysis defines heterogeneity of the murine pancreatic ductal tree.
      Finally, we conducted SEM on ADM and identified cells with admixed mucin and digestive enzyme–containing granules, similar to the phenotype of SPEM associated with gastric injury (Figure 3C).
      Figure thumbnail gr3
      Figure 3ADM as a pyloric-type metaplasia. (A) Expression of SPEM marker genes Tff2, Cd44, and Aqp5, chief cell marker gene Gif, and TF Gata5 overlaid on the UMAP of EYFP+ cells. (B) Co-IF for SPEM markers TFF2 (yellow), CD44v9 (cyan), and AQP5 or GIF (magenta) and DAPI (blue). Scale bars, 25 μm. (C) EM of mature mucus neck cells (black arrowheads) in the normal stomach, SPEM in an acute model of gastric injury (white arrowheads), and mucin/ductal cells from pancreatic ADM. Analogous granules in ADM are marked by arrowheads. Scale bars, 2 μm. (D) UMAP of scRNA-seq data collected from normal murine stomach and a SPEM model of acute gastric injury, from Bockerstett et al.
      • Bockerstett K.A.
      • Lewis S.A.
      • Wolf K.J.
      • et al.
      Single-cell transcriptional analyses of spasmolytic polypeptide-expressing metaplasia arising from acute drug injury and chronic inflammation in the stomach.
      (E) Expression of chief cell and acute SPEM signatures overlaid on the UMAP of EYFP+ cells. Gene signatures generated from (D). EC, enterochromaffin; Prolif, proliferating.
      Based on these findings, and prior studies noting several of these markers in pancreatic intraepithelial neoplasia (PanIN), we hypothesized that acinar cells undergo a metaplasia within the general rubric of pyloric metaplasia with specific parallels to SPEM.
      • Schlesinger Y.
      • Yosefov-Levi O.
      • Kolodkin-Gal D.
      • et al.
      Single-cell transcriptomes of pancreatic preinvasive lesions and cancer reveal acinar metaplastic cells' heterogeneity.
      ,
      • Prasad N.B.
      • Biankin A.V.
      • Fukushima N.
      • et al.
      Gene expression profiles in pancreatic intraepithelial neoplasia reflect the effects of Hedgehog signaling on pancreatic ductal epithelial cells.
      To test this hypothesis, we compared our ADM dataset to a recently published scRNA-seq study by Bockerstett et al,
      • Bockerstett K.A.
      • Lewis S.A.
      • Wolf K.J.
      • et al.
      Single-cell transcriptional analyses of spasmolytic polypeptide-expressing metaplasia arising from acute drug injury and chronic inflammation in the stomach.
      comparing normal gastric epithelium with that in several models of SPEM (Figure 3D, Supplementary Figure 8A–C). When we overlaid DEGs from normal gastric cell types onto our EYFP+ dataset, we found the chief cell transcriptome to be most closely related to pancreatic acinar cells (Figure 3E, Supplementary Figure 8A, Supplementary Table 5). To generate SPEM signatures, we identified DEGs between SPEM populations and normal chief cells in either a mouse model of acute gastric injury (73 genes) or of chronic autoimmune injury (54 genes), and overlaid these signatures on our ADM dataset (Supplementary Table 6). We found both SPEM signatures to be enriched in our mucin/ductal population (Figure 3E, Supplementary Figure 8C).
      We also identified TFs common to ADM and SPEM. Gata5 and several mucin-related TFs (Spdef, Creb3l1, and Creb3l4) are enriched in the mucin/ductal cluster (Supplementary Figure 6B and 8D).
      • Gregorieff A.
      • Stange D.E.
      • Kujala P.
      • et al.
      The ets-domain transcription factor Spdef promotes maturation of goblet and paneth cells in the intestinal epithelium.
      ,
      • Sakamoto N.
      • Fukuda K.
      • Watanuki K.
      • et al.
      Role for cGATA-5 in transcriptional regulation of the embryonic chicken pepsinogen gene by epithelial-mesenchymal interactions in the developing chicken stomach.
      ,
      • Asada R.
      • Saito A.
      • Kawasaki N.
      • et al.
      The endoplasmic reticulum stress transducer OASIS is involved in the terminal differentiation of goblet cells in the large intestine.
      To determine if these transcriptional programs are also activated in SPEM, we examined regulon module activity in both the acute and chronic SPEM scRNA-seq datasets.
      • Aibar S.
      • Gonzalez-Blas C.B.
      • Moerman T.
      • et al.
      SCENIC: single-cell regulatory network inference and clustering.
      We found that these regulons, identified in ADM, are also elevated in SPEM, suggesting similar regulatory programs in both disease states (Supplementary Figure 8E and F).

       ADM Results in EEC Heterogeneity

      We next explored heterogeneity within our lineage-traced Tuft/EEC population. Reclustering of the Tuft+EEC subset revealed 9 cell clusters separated into Tuft/EEC progenitors, EECs, or tuft lineages (Figure 4A, Supplementary Table 7). Trajectory inference using Monocle3 and RNA velocity suggests that EEC and tuft cell lineages originate from a Sox4+ common progenitor population (Figure 4B, C, and F), which has previously been reported in the intestines.
      • Gracz A.D.
      • Samsa L.A.
      • Fordham M.J.
      • et al.
      Sox4 Promotes Atoh1-independent intestinal secretory differentiation toward tuft and enteroendocrine fates.
      Within the EEC population, we identified 5 distinct clusters, including a progenitor-like cluster (Neurog3+, Figure 4D), and putative enterochromaffin (Ddc+Tac1+Chgb+), delta (Sst+), epsilon (Ghrl+), and PP/gamma cells (Ppy+) (Figure 4E, Supplementary Figure 9A).
      • Pan F.C.
      • Bankaitis E.D.
      • Boyer D.
      • et al.
      Spatiotemporal patterns of multipotentiality in Ptf1a-expressing cells during pancreas organogenesis and injury-induced facultative restoration.
      Gast (gastrin) expression is detected in a subset of Ddc+Tac1+ cells; Cck (cholecystokinin) in a subset of Ghrl+ cells. Insulin is not detected; minor glucagon expression is present within the Ppy+ cluster (Supplementary Figure 9A and B)
      Figure thumbnail gr4
      Figure 4ADM results in substantial enteroendocrine cell heterogeneity. (A) UMAP of cell clusters identified in the combined tuft and EEC dataset identified in B and labeled by cell type. (B) Monocle3 trajectory analysis overlaid on the UMAP in (A). (C) Expression of TFs Sox4 and (D) Neurog3 overlaid on the UMAP from (A). (E) Expression of epsilon cell marker Ghrl, enterochromaffin cell marker Ddc, delta cell marker, Sst, and PP/gamma cell marker Ppy overlaid on the UMAP from (A). (F) RNA velocity analysis of all 4 datasets. Arrow direction and length indicate probable lineage trajectory and velocity, respectively. (G) IHC for ghrelin (Ghrl), serotonin (generated by Ddc), somatostatin (Sst), and pancreatic polypeptide (Ppy) in normal and injured pancreata. Black arrows, positive cells. I, islets. Scale bar, 25 μm. (H) Transmission electron microscopy of an enterochromaffin cell or a delta cell in the injured pancreas. Scale bars, 5 and 2 μm, respectively.
      To confirm expression of these hormones at the protein level, we performed IHC and detected ghrelin, serotonin (Ddc and Tph1), somatostatin, pancreatic polypeptide, and gastrin in injury-induced ADM (Figure 4G, Supplementary Figure 9F). To determine if these hormone-producing EECs are truly associated with ADM and not islets, we conducted lineage tracing with CERTY mice and identified coexpression of EYFP in all 4 EEC subtypes (Supplementary Figure 9G). Analysis of EM conducted on injured pancreata identified cells with enterochromaffin cell (serotonin-producing, pleomorphic granules) and delta cell (somatostatin-producing, round granules) morphology (Figure 4H).
      • Gunawardene A.R.
      • Corfe B.M.
      • Staton C.A.
      Classification and functions of enteroendocrine cells of the lower gastrointestinal tract.
      To identify putative regulatory processes involved in EEC subtype specification, we examined the activity of key TFs and signaling molecules. We found that the Neurog3+ EEC progenitor state is associated with expression of Notch pathway inhibitors and subsequent branching into enterochromaffin and epsilon cell lineages is associated with changes in key TFs.
      • Bastidas-Ponce A.
      • Tritschler S.
      • Dony L.
      • et al.
      Comprehensive single cell mRNA profiling reveals a detailed roadmap for pancreatic endocrinogenesis.
      Enterochromaffin lineage formation correlates with Nkx6–1 activation and a switch from Pax4 to Pax6 expression.
      • Oster A.
      • Jensen J.
      • Edlund H.
      • et al.
      Homeobox gene product Nkx 6.1 immunoreactivity in nuclei of endocrine cells of rat and mouse stomach.
      The epsilon lineage is associated with activation of Arx and Isl1. Subsequent branching into delta and PP/gamma cells is accompanied by reactivation of Pax6. Delta cells are characterized by a lack of Arx, reactivation of Sox9, and expression of Pou3f1 (Supplementary Figure 9C–E).
      • Baron M.
      • Veres A.
      • Wolock S.L.
      • et al.
      A single-cell transcriptomic map of the human and mouse pancreas reveals inter- and intra-cell population structure.
      • Heller R.S.
      • Jenny M.
      • Collombat P.
      • et al.
      Genetic determinants of pancreatic epsilon-cell development.
      • Andralojc K.M.
      • Mercalli A.
      • Nowak K.W.
      • et al.
      Ghrelin-producing epsilon cells in the developing and adult human pancreas.

       Transcriptomic Changes Associated With Tuft Cell Development

      A scRNA-seq study of intestinal epithelium previously demonstrated that tuft cells can be categorized as either neuronal or inflammatory-like based on marker expression.
      • Haber A.L.
      • Biton M.
      • Rogel N.
      • et al.
      A single-cell survey of the small intestinal epithelium.
      To identify tuft subpopulations in our dataset, we conducted DEG and GO Term analysis of tuft cell clusters (6, 5, 1, 0), including a total of 1280 cells (Figure 5A, Supplementary Figure 10A). Consistent with RNA velocity results (Figure 4F), we observed a transition from a Sox4+ progenitor population (cluster 6) to a Pou2f3+Spib+ population (clusters 5 and 1); Pou2f3 is the master regulator for tuft cell formation (Figure 5A–C).
      • Yamashita J.
      • Ohmoto M.
      • Yamaguchi T.
      • et al.
      Skn-1a/Pou2f3 functions as a master regulator to generate Trpm5-expressing chemosensory cells in mice.
      Subsequent clusters 1 and 0 are enriched for known tuft cell markers and signaling molecules (ie, Trpm5, Il25) (Figure 5A and J), suggesting that our dataset reflects the transcriptomic changes associated with tuft cell development, rather than separate, functional subtypes.
      Figure thumbnail gr5
      Figure 5Transcriptomic changes associated with tuft cell formation. (A) Heatmap of scaled DEGs for the 4 clusters of the tuft cell lineage (from A). A selection of known marker genes is labeled on the left. Cells were ordered by increasing pseudotime within each cluster group. (B) Smoothed regulon activity score trend of the tuft cell lineage in Monocle3 pseudotime. (C and D) Regulon activity scores of Pou2f3, Spib, and Ascl1 overlaid on the Tuft+EEC UMAP. (E) Heatmap of scaled RNA expression of tuft-related TFs. Cells are ordered by increasing pseudotime within each cluster group. (F) The aggregated fate probability toward terminal tuft lineage formation. (G–K) Expression of select genes overlaid on the Tuft+EEC UMAP, labeled by Seurat cluster.
      To test this hypothesis, we conducted regulon analysis and identified candidate TFs and downstream regulatory networks associated with tuft cell development.
      • Aibar S.
      • Gonzalez-Blas C.B.
      • Moerman T.
      • et al.
      SCENIC: single-cell regulatory network inference and clustering.
      We examined regulon activities as a function of Monocle3 Pseudotime, moving from cluster 6 → 5 → 1 → 0. Consistent with RNA velocity and DEG analysis, we identified a decrease in Sox4 activity, and an increase in Pou2f3, Spib, and Ascl1 activity (Figure 5B–D). Additional regulators include neuronal development-associated TFs (Tead2, Hmgb3, Mef2b), TFs associated with hematopoietic differentiation and inflammatory cell function (Runx1, Myb, Gfi1b, Nfatc1, and Foxe1), and TFs associated with epithelial differentiation (Ehf, Klf4) (Figure 5E).
      • Ramsay R.G.
      • Gonda T.J.
      MYB function in normal and cancer cells.
      • Kas K.
      • Finger E.
      • Grall F.
      • et al.
      ESE-3, a novel member of an epithelium-specific ets transcription factor subfamily, demonstrates different target gene specificity from ESE-1.
      • Ghaleb A.M.
      • Yang V.W.
      Kruppel-like factor 4 (KLF4): What we currently know.
      We next examined cell fate probability using a directed single-cell fate mapping approach. Cells were decomposed into macrostates matching our seurat clusters (Supplementary Figure 10B). Directed fate probabilities reveal a continuous process of tuft cell formation, which is shown by a steady gain in tuft cell lineage fate probabilities moving from cluster 6 → 0, consistent with tuft cell development and maturation (Figure 5F, Supplementary Figure 10C). This is consistent with transcriptional activity changes in many tuft cell marker genes. The first group we identified includes genes enriched in cluster 5, in which expression fades throughout development (Gnat3, Mef2b, Nrep) (Figure 5G, Supplementary Figure 10E). The second group represents a steady level of gene expression throughout the pseudotemporal ordering from cluster 5 to cluster 0 (Pou2f3, Gng13, Ptgs1, and so on) (Figure 5H, Supplementary Figure 10F). The last group includes genes associated with the late stage of tuft cell formation where expression increases along tuft cell development (Siglecf, Ptprc, Fcna) (Figure 5I, Supplementary Figure 10G).
      Finally, we explored the tuft cell classifications previously reported in small intestines in our dataset by overlaying consensus gene signatures.
      • Haber A.L.
      • Biton M.
      • Rogel N.
      • et al.
      A single-cell survey of the small intestinal epithelium.
      We found that the inflammatory signature is enriched late in tuft cell development, whereas the neuronal signature is enriched at a relatively early stage and remains expressed (Supplementary Figure 10D). Recently, Manco et al
      • Manco R.
      • Averbukh I.
      • Porat Z.
      • et al.
      Clump sequencing exposes the spatial expression programs of intestinal secretory cells.
      used “clump sequencing” to identify spatial expression programs of intestinal secretory cells. The authors found that tuft cells associated with intestinal crypts, where cells are more stem-like, express Sox4 and are enriched for the neuronal tuft gene signature. Tuft cells associated with the intestinal villus, where cells are most differentiated, are enriched for functional genes like Il17rb, Fapb1, and the inflammatory tuft signature. Consistent with this, we identified expression of Il17rb, Fabp1, and a number of other secretory signaling molecules and hormones enriched in differentiated tuft cell cluster 0 (Figure 5J–L, Supplementary Figure 10G). Collectively, these data identify regulatory programs and gene expression changes associated with tuft cell maturation.

       KrasG12D Expression Drives PanIN-specific Changes in ADM

      ADM occurs in response to oncogenic KrasG12D and may serve as a precursor lesion for PDAC. Recently, Schlesinger et al
      • Schlesinger Y.
      • Yosefov-Levi O.
      • Kolodkin-Gal D.
      • et al.
      Single-cell transcriptomes of pancreatic preinvasive lesions and cancer reveal acinar metaplastic cells' heterogeneity.
      published an scRNA-seq study of pancreatic tumorigenesis including lineage tracing using LSL-KrasG12D;Ptf1aCreERTM;LSL-tdTomato (KCT) mice. To identify similarities and differences between injury and oncogene-induced ADM, we compared EYFP+ cells (13,362 cells) with tdTomato+ cells that were fluorescence activated cell sorting–collected from a 6-month-old KCT mouse (2719 cells, Supplementary Figure 11A). We integrated the 2 datasets using fastMNN and identified 18 distinct clusters (Supplementary Figure 11B and C). Using the aforementioned gene signatures, as well as those in Schlesinger et al,
      • Schlesinger Y.
      • Yosefov-Levi O.
      • Kolodkin-Gal D.
      • et al.
      Single-cell transcriptomes of pancreatic preinvasive lesions and cancer reveal acinar metaplastic cells' heterogeneity.
      we found a number of populations common to both disease states (eg, tuft and EECs; Figure 6A). Metaplasia transcription factors Foxq1 and Onecut2 (described in Schlesinger et al
      • Schlesinger Y.
      • Yosefov-Levi O.
      • Kolodkin-Gal D.
      • et al.
      Single-cell transcriptomes of pancreatic preinvasive lesions and cancer reveal acinar metaplastic cells' heterogeneity.
      ) are also expressed in injury-induced ADM (Supplementary Figure 11D and E). Acinar populations in KrasG12D ADM are enriched with regenerating (REG) protein family members (Reg2+Reg3a+Reg3g+ “Acinar-REG+”), whereas those associated with injury-induced ADM express higher levels of ribosome transcripts; consistent with the recently identified “secretory acinar cell” phenotype (“Acinar-S”) thought to function primarily in production of digestive enzymes (Supplementary Figure 11F).
      • Tosti L.
      • Hang Y.
      • Debnath O.
      • et al.
      Single nucleus and in situ RNA sequencing reveals cell topographies in the human pancreas.
      Figure thumbnail gr6
      Figure 6KrasG12D expression drives PanIN-specific changes in ADM. (A) UMAP of FastMNN integrated datasets from the injury-induced ADM scRNA-seq dataset and an oncogene-induced ADM dataset from a 6-month-old KCT mouse derived from Schlesigner et al.
      • Schlesinger Y.
      • Yosefov-Levi O.
      • Kolodkin-Gal D.
      • et al.
      Single-cell transcriptomes of pancreatic preinvasive lesions and cancer reveal acinar metaplastic cells' heterogeneity.
      (B) Expression of pit cell markers Gkn1 and Muc5ac and (C) pre-tumor marker Msln. (D) IHC for MUC5AC or MSLN in the normal stomach or pancreas as well as in injury or oncogene (KrasG12D;Ptf1aCre/+)-induced ADM and neoplasia. (E and F) Hematoxylin-eosin (H&E), Periodic acid-Schiff (PAS)/AB (mucins), or MUC5AC expression in KHERT or KPERT mice treated with tamoxifen, HERT or PERT mice treated with caerulein or caerulein followed by tamoxifen, and KHERT and KPERT mice first treated with caerulein followed by tamoxifen. Scale bars, 50 μm.
      Differences between the ADM populations in the 2 datasets include a lack of senescence markers, and “metaplastic pre-tumor” populations in the injury model (Figure 6A).
      • Schlesinger Y.
      • Yosefov-Levi O.
      • Kolodkin-Gal D.
      • et al.
      Single-cell transcriptomes of pancreatic preinvasive lesions and cancer reveal acinar metaplastic cells' heterogeneity.
      Further, we identified distinct differences in the mucin/ductal population between the 2 disease states including a shift from Muc6-high,Tff2-high,Krt19-low cells to Muc6-low,Tff2-low,Krt19-high cells in the KCT pancreas, possibly reflecting enhanced disease progression (Figure 6A, Supplementary Figure 11B). DEG analysis of the mucin/ductal cell clusters between the 2 disease states identified significantly higher expression of SPEM markers in injury-induced ADM, consistent with the presence of a larger Muc6-high,Tff2-high,Krt19-low population. Oncogene-induced ADM, in contrast, is characterized by significantly higher expression of proto-oncogenes such as Areg, Fos, and Junb, known downstream targets of KrasG12D (Supplementary Figure 11G).
      • Deng T.
      • Karin M.
      c-Fos transcriptional activity stimulated by H-Ras-activated protein kinase distinct from JNK and ERK.
      Schlesinger et al
      • Schlesinger Y.
      • Yosefov-Levi O.
      • Kolodkin-Gal D.
      • et al.
      Single-cell transcriptomes of pancreatic preinvasive lesions and cancer reveal acinar metaplastic cells' heterogeneity.
      described the formation of a gastric pit cell-like population (Tff1+Gkn1+Muc5ac+) in oncogene-induced ADM (Figure 6A and B). To assay injury-induced ADM for pit cells, we conducted IHC for MUC5AC and found expression to be specific to KrasG12D pancreata. Interestingly, a number of pre-tumor markers described by Schlesinger et al
      • Schlesinger Y.
      • Yosefov-Levi O.
      • Kolodkin-Gal D.
      • et al.
      Single-cell transcriptomes of pancreatic preinvasive lesions and cancer reveal acinar metaplastic cells' heterogeneity.
      (Msln, Mmp7, Marckls1, and Igfbp7) are, in fact, expressed in injury-induced ADM (Figure 6C and D, Supplementary Figure 11H and I).
      • Crawford H.C.
      • Scoggins C.R.
      • Washington M.K.
      • et al.
      Matrix metalloproteinase-7 is expressed by pancreatic cancer precursors and regulates acinar-to-ductal metaplasia in exocrine pancreas.
      To confirm expression of these pre-tumor markers, we conducted IHC for MSLN and detected positive cells within the epithelium of injured pancreata (Figure 6D). These data suggest that “pre-tumor” cells represent a transitional state co-opted by oncogenic Kras during tumorigenesis.
      To determine if KrasG12D expression is sufficient to drive the aforementioned differences in injury- vs oncogene-induced ADM, we used several inducible Cre drivers to target KrasG12D expression to ADM populations. Hnf1b:CreERT2 mice have been described previously.
      • Solar M.
      • Cardalda C.
      • Houbracken I.
      • et al.
      Pancreatic exocrine duct cells give rise to insulin-producing beta cells during embryogenesis but not after birth.
      In the adult pancreas, Hnf1b expression is restricted to ductal cells, and Bailey et al
      • Bailey J.M.
      • Hendley A.M.
      • Lafaro K.J.
      • et al.
      p53 mutations cooperate with oncogenic Kras to promote adenocarcinoma from pancreatic ductal cells.
      have shown that 2 alleles of Trp53R172H in addition to KrasG12D are required to induce transformation. KrasG12D expression alone results in no discernable phenotype in the pancreas. To determine if KrasG12D expression targeted to injury-induced Hnf1b+ ductal cells is sufficient to drive PanIN-associated changes, we bred KrasG12D;Hnf1b:CreERT2;RosaLSL-EYFP/+ (KHERTY) mice, induced injury and ADM with caerulein first, and then activated KrasG12D expression after ADM had formed (Supplementary Figure 12A–C). We found that, although the Hnf1b:CreERT2 pancreas heals after injury and tamoxifen treatment, targeting KrasG12D expression to ADM results in 100% (n = 6) of KHERTY pancreata forming MUC5AC+ PanIN (Figure 6E).
      As the Hnf1b:CreERT2 mouse model drives KrasG12D expression from both normal and ADM-induced ductal cells, we next conducted this experiment using ADM-specific Pou2f3-creErt2;RosamT/mG mice.
      • McGinty J.W.
      • Ting H.A.
      • Billipp T.E.
      • et al.
      Tuft-cell-derived leukotrienes drive rapid anti-helminth immunity in the small intestine but are dispensable for anti-protist immunity.
      We have previously shown that POU2F3 expression is absent from the normal pancreas, that 99% of POU2F3+ tuft cells are derived from ADM, and that tuft cells are rare in the injured pancreas (making up ∼0.6% of the epithelium).
      • DelGiorno K.E.
      • Naeem R.F.
      • Fang L.
      • et al.
      Tuft cell formation reflects epithelial plasticity in pancreatic injury: implications for modeling human pancreatitis.
      ,
      • DelGiorno K.E.
      • Chung C.Y.
      • Vavinskaya V.
      • et al.
      Tuft cells inhibit pancreatic tumorigenesis in mice by producing prostaglandin D2.
      As shown in Figure 5H, Pou2f3 is expressed early in injury-induced EEC/tuft specification. To determine if ADM-induced Pou2f3+ cells are susceptible to oncogenic transformation and if KrasG12D expression is sufficient to drive the changes identified in our scRNA-seq analyses, we bred KrasG12D; Pou2f3-creErt2;RosamT/mG (KPERT) mice. We found that the Pou2f3-creErt2 pancreas heals after recovery from injury and tamoxifen treatment; however, 40% (4 of 10 mice) of KPERT pancreata contained MUC5AC+ PanIN (Figure 6F, Supplementary Figure 12D–F). PanIN are sparse in this model, consistent with the rarity of POU2F3-positive cells in injury-induced ADM. Collectively, these data identify changes in the epithelium between injury and oncogene-induced ADM and demonstrate that expression of KrasG12D targeted to injury-induced ADM populations is sufficient to drive these changes.

       Human Pancreatic Metaplasia as a Pyloric-type Transition

      To determine if the epithelial composition of human chronic pancreatitis reflects the diversity identified in mouse models, we evaluated marker expression in sNuc-seq data collected from human pancreata and conducted immunofluorescence on tissue samples. Recently, Tosti et al
      • Tosti L.
      • Hang Y.
      • Debnath O.
      • et al.
      Single nucleus and in situ RNA sequencing reveals cell topographies in the human pancreas.
      generated a comprehensive human pancreas cell atlas consisting of approximately 113,000 nuclei from healthy adult donors and 2726 nuclei from 2 patients with chronic pancreatitis. Consistent with mouse data, the pancreatitis dataset reflects an increase in a mucin/ductal population (MUC5B+/KRT19+; “MUC5B+Ductal”) and the appearance of tuft cells, as compared with the normal pancreas (Figure 7A, Supplementary Figure 13A). To directly compare the human and murine mucin/ductal populations, we overlaid their respective gene signatures on their counterpart species dataset and found enrichment in the Mucin/Ductal clusters (Supplementary Figure 13B).
      Figure thumbnail gr7
      Figure 7Human chronic pancreatitis as a pyloric-type transition. (A) UMAPs showing annotated clusters from sNuc-seq data collected from normal human pancreata (∼113k nuclei) or the injured pancreata of 2 patients with chronic pancreatitis (2726 nuclei) derived from Tosti et al.
      • Tosti L.
      • Hang Y.
      • Debnath O.
      • et al.
      Single nucleus and in situ RNA sequencing reveals cell topographies in the human pancreas.
      (B). Expression of the humanized acute SPEM signature (D) overlaid on the UMAPs from (A). (C) Co-IF for SPEM markers (TFF2, magenta; AQP5, yellow; CD44v9, cyan), (D) epsilon and delta cell markers (Ghrl, magenta; SST, yellow; γActin, cyan), and (E) gamma and enterochromaffin cell markers (PPY, magenta; 5-HT, yellow; γActin, cyan). DAPI, blue. Scale bars, 25 μm.
      To determine if human ADM represents a pyloric-type metaplasia transition, we humanized and overlaid the murine acute and chronic SPEM signatures on the human pancreatitis dataset. In both analyses, we found enrichment of the SPEM signature in the MUC5B+Ductal population (Figure 7B, Supplementary Figure 13C). Analysis of individual markers identifies TFF2, WFDC2, and CD44 (Supplementary Figure 13A). To validate these findings, we conducted immunostaining on pancreata from 15 patients with chronic pancreatitis, prescreened by a pathologist for the presence of metaplasia, but varying in terms of etiology and severity. First, we evaluated the expression of SPEM markers TFF2, AQP5, and CD44v9 (generated by CD44). As compared with normal pancreas, we identified up-regulation of CD44v9 and AQP5 in acinar cells in areas of tissue damage, in all 15 samples (Supplementary Figure 13D). Occasionally, TFF2 was coexpressed in these regions (triple-positive clusters present in 13 of 15 samples), but was more often found in metaplastic structures lacking acinar morphology, in PanIN-like lesions, and in duct glands associated with large ducts (Figure 7C).
      We next evaluated the presence of tuft cells and EECs in human chronic pancreatitis. Consistent with prior reports, we identified tuft cells by expression of marker phospho-Girdin in areas of metaplasia associated with injured acini, in normal ducts, and in incidental PanIN (10 of 15 patient samples; Supplementary Figure 13E).
      • Delgiorno K.E.
      • Hall J.C.
      • Takeuchi K.K.
      • et al.
      Identification and manipulation of biliary metaplasia in pancreatic tumors.
      ,
      • Kuga D.
      • Ushida K.
      • Mii S.
      • et al.
      Tyrosine phosphorylation of an actin-binding protein girdin specifically marks tuft cells in human and mouse gut.
      When we interrogated the sNuc-seq dataset for tuft cell TFs identified in our EYFP+ dataset, we found expression of SOX4POU2F3, SPIB, and RUNX1 in human tuft cells (Supplementary Figure 13F). To assess EEC heterogeneity, we conducted IF for ghrelin (Ghrl), somatostatin (SST), pancreatic polypeptide (PPY), and serotonin (5-HT). All 4 hormones could be detected in islets, although the proportions of each varied by patient. All hormones were associated with aberrant ductal structures, rather than injured acini, consistent with formation later in ADM. Ghrelin and PPY expression is rare in metaplasia, but both can be found in large ducts. 5-HT and SST are expressed in metaplasia, PanIN-like lesions, and in intralobular ducts (Figure 7D and E). Finally, we assessed the expression of pre-tumor markers in the pancreatitis sNuc-seq dataset. We found expression of MMP7 and MARCKSL1, but not IGFBP7, in MUC5B+Ductal cells, although it is not specific. MSLN, though, could be detected only in MUC5B+Ductal cells (Supplementary Figure 13G).
      Overall, the staining patterns in humans are largely reminiscent of the murine model, including the formation of distinct cell types (mucin/ductal, tuft, EECs) and expression of SPEM signatures. Differences include the lack of TFF2 in injured acini and expression of individual genes (ie, Muc6 vs MUC5B). Further studies are required to determine if these differences reflect the heterogeneity of etiology, chronicity, and severity of chronic pancreatitis in patients selected to undergo surgical resection.

      Discussion

      Using lineage tracing, scRNA-seq, multiple computational biology approaches (that rely on nonoverlapping assumptions), immunostaining, and high-resolution ultramicroscopy, we have identified significantly more epithelial heterogeneity in injury-induced pancreatic ADM than previously appreciated. Brainbow lineage tracing demonstrates that acinar cells can clonally expand, forming the secretory cell types identified in these analyses. Although scRNA-seq provides information regarding cell-type markers, high-resolution ultramicroscopy is critical to the identification of these populations, as it provides structural information, corroborating the formation of defined cell types. Collectively, these data are consistent with a model by which acinar cells undergo a pyloric-type metaplasia resulting in mucinous ductal cells; it is from this state that differentiated tuft cells or EECs emerge.
      We have discovered that chronic injury-induced ADM results in substantial EEC heterogeneity. ADM-derived EEC formation was previously described by Farrell et al
      • Farrell A.S.
      • Joly M.M.
      • Allen-Petersen B.L.
      • et al.
      MYC regulates ductal-neuroendocrine lineage plasticity in pancreatic ductal adenocarcinoma associated with poor outcome and chemoresistance.
      in PanIN resulting from expression of oncogenic KrasG12D by lineage tracing and expression of marker synaptophysin. The term “enteroendocrine cell,” however, represents a lineage populated by multiple subtypes characterized by the expression, or coexpression, of various hormones. Although best known for regulating intestinal function, insulin secretion, nutrient assimilation, and food intake in normal tissues, evidence points to poorly understood roles for several of these hormones in inflammation and tumorigenesis.
      • Gunawardene A.R.
      • Corfe B.M.
      • Staton C.A.
      Classification and functions of enteroendocrine cells of the lower gastrointestinal tract.
      ,
      • Gribble F.M.
      • Reimann F.
      Function and mechanisms of enteroendocrine cells and gut hormones in metabolism.
      In particular, somatostatin is known to have diabetogenic properties, which, when combined with the observation that adult-onset diabetes is a known association of pancreatic adenocarcinoma, suggests the intriguing idea that tumor-derived EECs have an important role in the endocrine/metabolic symptoms associated with this malignancy.
      • Molina-Montes E.
      • Coscia C.
      • Gomez-Rubio P.
      • et al.
      Deciphering the complex interplay between pancreatic cancer, diabetes mellitus subtypes and obesity/BMI through causal inference and mediation analyses.
      These data underscore the importance of understanding the contribution of EECs, including their regulatory processes and secretory products, in the contexts of pancreatic injury and tumorigenesis.
      Altogether, we propose that this carefully orchestrated process of acinar cell de- and trans-differentiation exists to mitigate injury and may invoke a subset of a generalized pyloric metaplasia response that can occur in many organs in the gastrointestinal tract.
      • Goldenring J.R.
      Pyloric metaplasia, pseudopyloric metaplasia, ulcer-associated cell lineage and spasmolytic polypeptide-expressing metaplasia: reparative lineages in the gastrointestinal mucosa.
      Ever more evidence is accumulating that there are shared patterns of regenerative changes that depend on cellular plasticity. Here, we observe striking transcriptional overlap between metaplastic cells in the pancreas and those of the stomach. Although care should be taken when comparing datasets derived from different strains of mice, our results are in line with recent studies indicating that mature cells like gastric chief cells and pancreatic acinar cells have a shared injury-induced plasticity program. This “paligenosis” program governs their conversion from differentiated to regenerative states.
      • Willet S.G.
      • Lewis M.A.
      • Miao Z.F.
      • et al.
      Regenerative proliferation of differentiated cells by mTORC1-dependent paligenosis.
      ,
      • Nam K.T.
      • Lee H.J.
      • Sousa J.F.
      • et al.
      Mature chief cells are cryptic progenitors for metaplasia in the stomach.
      ,
      • Miao Z.F.
      • Lewis M.A.
      • Cho C.J.
      • et al.
      A dedicated evolutionarily conserved molecular network licenses differentiated cells to return to the cell cycle.
      Cycles of paligenosis that fuel pyloric metaplasia may allow cells to accumulate mutations and increase the risk for progression to cancer.
      Interestingly, many of the epithelial populations identified in injury-induced ADM, such as tuft and EECs, are shared by oncogenic KrasG12D-induced ADM, suggesting that they represent an early event in tumorigenesis. It is the differences between the 2 disease states, however, that may provide important clues as to which population(s) signal transformation or progression to cancer. For example, the formation of pit-like cells (Muc5ac+Gkn1+Tff1+) is specific to KrasG12D-induced ADM.
      • Schlesinger Y.
      • Yosefov-Levi O.
      • Kolodkin-Gal D.
      • et al.
      Single-cell transcriptomes of pancreatic preinvasive lesions and cancer reveal acinar metaplastic cells' heterogeneity.
      Here, we demonstrate that activation of KrasG12D in injury-induced ADM is sufficient to drive the phenotypic changes associated with PanIN formation, as well as the generation of MUC5AC+ pit-like cells. These data have important implications with regard to tumor cell of origin. Previous studies have focused on the role of acinar and ductal cells in tumorigenesis; however, our data suggest that the plethora of cell types that form in response to injury are susceptible to oncogenic transformation. The particular cell state or cell type in which oncogenic mutation occurs may factor into the speed, aggressiveness, or sub-type of PDAC that develops. Further, these data shed new light into mechanisms by which pancreatitis predisposes patients to PDAC.

      Acknowledgments

      We thank Richard DiPaolo, Christopher Wright, Eunyoung Choi, Ela Contreras Panta, and members of the Goldenring, Choi, and DelGiorno laboratories for helpful discussions. We thank Dr Tom Bartol for his help with 3DEM stack alignment using SWiFT-IR (funded by National Institutes of Health P41GM103712 and National Science Foundation DBI-1707356 ). Finally, we thank the scientists whose work this study builds on, but cannot be discussed here due to space constraints. Sequencing data that support this study are available in the GEO archive (GSE172380). Biorxiv doi: https://doi.org/10.1101/2021.04.09.439243

      CRediT Authorship Contributions

      Zhibo Ma, PhD (Conceptualization: Supporting; Data curation: Lead; Formal analysis: Lead; Investigation: Lead; Visualization: Equal; Writing – original draft: Equal; Writing – review & editing: Equal)
      Nikki K. Lytle, PhD (Conceptualization: Supporting; Data curation: Lead; Investigation: Lead; Methodology: Lead; Writing – original draft: Supporting; Writing – review & editing: Supporting)
      Bob Chen, BS (Data curation: Equal; Formal analysis: Equal; Investigation: Equal; Methodology: Equal; Writing – original draft: Supporting; Writing – review & editing: Supporting)
      Nidhi Jyotsana, PhD (Data curation: Equal; Formal analysis: Supporting; Writing – original draft: Supporting; Writing – review & editing: Supporting)
      Sammy Weiser Novak, MS (Conceptualization: Equal; Data curation: Lead; Formal analysis: Lead; Investigation: Equal; Methodology: Equal; Resources: Equal; Visualization: Equal; Writing – original draft: Supporting; Writing – review & editing: Supporting)
      Charles Cho, MD, PhD (Formal analysis: Equal; Investigation: Equal; Writing – original draft: Supporting; Writing – review & editing: Supporting)
      Leah Caplan, BS (Data curation: Equal; Formal analysis: Equal; Writing – original draft: Supporting; Writing – review & editing: Supporting)
      Olivia Ben-Levy, BS (Data curation: Equal; Formal analysis: Supporting; Writing – original draft: Supporting; Writing – review & editing: Supporting)
      Abigail C. Neininger, BS (Data curation: Equal; Formal analysis: Equal; Investigation: Equal; Methodology: Supporting; Writing – original draft: Supporting; Writing – review & editing: Supporting)
      Dylan T. Burnette, PhD (Data curation: Equal; Formal analysis: Equal; Supervision: Equal; Visualization: Equal; Writing – original draft: Supporting; Writing – review & editing: Supporting)
      Vincent Q. Trinh, MD (Formal analysis: Equal; Resources: Equal; Writing – original draft: Supporting; Writing – review & editing: Supporting)
      Marcus C.B. Tan, MD (Data curation: Equal; Formal analysis: Equal; Supervision: Equal; Writing – original draft: Supporting; Writing – review & editing: Supporting)
      Emilee Patterson, BS (Data curation: Equal; Writing – original draft: Supporting; Writing – review & editing: Supporting)
      Rafael Arrojo e Drigo, PhD (Data curation: Supporting; Formal analysis: Equal; Supervision: Equal; Visualization: Supporting; Writing – original draft: Supporting; Writing – review & editing: Supporting)
      Rajshekhar Giraddi, PhD (Conceptualization: Equal; Writing – original draft: Supporting; Writing – review & editing: Supporting)
      Cynthia Ramos, MS (Data curation: Equal)
      Anna L. Means, PhD (Conceptualization: Equal; Writing – original draft: Supporting; Writing – review & editing: Supporting)
      Ichiro Matsumoto, PhD (Resources: Equal; Writing – review & editing: Supporting)
      Uri Manor, PhD (Funding acquisition: Supporting; Supervision: Equal; Visualization: Supporting; Writing – original draft: Supporting; Writing – review & editing: Supporting)
      Jason C. Mills, MD, PhD (Data curation: Equal; Supervision: Equal; Writing – original draft: Equal; Writing – review & editing: Equal)
      James R. Goldenring, MD, PhD (Conceptualization: Equal; Funding acquisition: Supporting; Resources: Equal; Writing – original draft: Equal; Writing – review & editing: Equal)
      Ken S. Lau, PhD (Data curation: Equal; Formal analysis: Equal; Funding acquisition: Supporting; Investigation: Equal; Methodology: Equal; Software: Equal; Supervision: Equal; Writing – original draft: Supporting; Writing – review & editing: Supporting)
      Geoffrey M. Wahl, PhD (Funding acquisition: Supporting; Resources: Equal; Supervision: Equal; Writing – original draft: Supporting; Writing – review & editing: Supporting)
      Kathleen E. DelGiorno, PhD (Conceptualization: Lead; Data curation: Equal; Formal analysis: Equal; Funding acquisition: Lead; Investigation: Equal; Project administration: Lead; Supervision: Lead; Validation: Equal; Visualization: Lead; Writing – original draft: Lead; Writing – review & editing: Lead)

      Supplementary Material

      • Loading ...

      Supplementary Materials and Methods

       Mice

      Mice were housed in accordance with National Institutes of Health guidelines in Association for Assessment and Accreditation of Laboratory Animal Care–accredited facilities at Vanderbilt University, the Salk Institute, and at Washington University. The Institutional Animal Care and Use Committees at Vanderbilt University, the Salk Institute, and Washington University approved all mouse procedures. Ptf1aCre-ERTM, RosaYFP/+, RosaBrainbow.2/+, RosamT/mG, and Hnf1b:CreERT2 mice were purchased from Jackson Laboratories (Bar Harbor, ME).
      • Pan F.C.
      • Bankaitis E.D.
      • Boyer D.
      • et al.
      Spatiotemporal patterns of multipotentiality in Ptf1a-expressing cells during pancreas organogenesis and injury-induced facultative restoration.
      ,
      • Solar M.
      • Cardalda C.
      • Houbracken I.
      • et al.
      Pancreatic exocrine duct cells give rise to insulin-producing beta cells during embryogenesis but not after birth.
      LSL-KrasG12D/+; Ptf1aCre/+ mice have been previously described in detail.
      • Hingorani S.R.
      • Petricoin E.F.
      • Maitra A.
      • et al.
      Preinvasive and invasive ductal pancreatic cancer and its early detection in the mouse.
      Mice were backcrossed into the CD-1 strain (Charles River, Charleston, SC). Pou2f3-creErt2 mice have been reported previously, were provided by Dr Ichiro Matsumoto (Monell Chemical Senses Center, Philadelphia, PA), and were backcrossed into the CD-1 strain.
      • McGinty J.W.
      • Ting H.A.
      • Billipp T.E.
      • et al.
      Tuft-cell-derived leukotrienes drive rapid anti-helminth immunity in the small intestine but are dispensable for anti-protist immunity.
      For stomach studies, wild-type C57BL/6 mice were purchased from The Jackson Laboratory and maintained in a specified-pathogen-free barrier facility under a 12-hour light cycle.

       Human Samples

      The Institutional Review Board of Vanderbilt University approved the distribution and use of all human samples.

       Lineage Tracing and Pancreatitis Induction

      Lineage tracing was conducted using the mouse model Ptf1aCre-ERTM;ROSAYFP/+, as previously described.
      • Pan F.C.
      • Bankaitis E.D.
      • Boyer D.
      • et al.
      Spatiotemporal patterns of multipotentiality in Ptf1a-expressing cells during pancreas organogenesis and injury-induced facultative restoration.
      ,
      • Delgiorno K.E.
      • Hall J.C.
      • Takeuchi K.K.
      • et al.
      Identification and manipulation of biliary metaplasia in pancreatic tumors.
      ,
      • DelGiorno K.E.
      • Naeem R.F.
      • Fang L.
      • et al.
      Tuft cell formation reflects epithelial plasticity in pancreatic injury: implications for modeling human pancreatitis.
      In this model, tamoxifen treatment induces Cre activity, which then initiates expression of enhanced yellow fluorescent protein (EYFP) specifically in Ptf1a+ acinar cells. Acinar cells were labeled with daily tamoxifen treatment (Sigma, St Louis, MO; 5 mg/d, 5 d/wk for 2 weeks) delivered in corn oil (Sigma) by oral gavage. Pancreatitis was induced with caerulein (Bachem, Bubendorf, Switzerland) administered intraperitoneally (IP) at 250 μg/kg, twice a day, for 5 days a week for either 2 or 4 weeks and mice were then allowed to recover for 2 days.
      • DelGiorno K.E.
      • Naeem R.F.
      • Fang L.
      • et al.
      Tuft cell formation reflects epithelial plasticity in pancreatic injury: implications for modeling human pancreatitis.
      Lineage tracing was conducted in Ptf1aCre-ERTM;ROSABrainbow.2/+mice with 2 daily treatments of 5 mg tamoxifen followed by the caerulein treatment protocol described previously. Pancreatitis was induced in KrasG12D;Hnf1b:CreERT2 mice (with or without ROSAYFP/+), Hnf1b:CreERT2;RosaYFP/+ mice, KrasG12D;Pou2f3-creErt2 mice, or Pou2f3-creErt2;RosamTmG mice with caerulein administered at 250 μg/kg, twice a day, for 5 days a week for 4 weeks followed by tamoxifen treatment (5 mg/d, 5 d/wk or tamoxifen chow for 5 days). Mice were allowed to recover for a minimum of 1 month or until they had to be killed due to lung tumors (KrasG12D;Hnf1b:CreERT2 mice) or large oral papillomas (KrasG12D;Pou2f3-creErt2 mice).

       Isolation of Pancreatic Single Cell Suspensions

      EYFP+ cells were isolated from the pancreata of caerulein-treated mice. The pancreas was quickly dissected and minced in Hank’s Balanced Salt Solution. Supernatant and fat were removed and pancreatic tissue was then incubated in 10 mL Dulbecco’s modified Eagles’ medium supplemented with 1 mg/mL collagenase I (Sigma), 1 mg/mL soybean trypsin inhibitor (Gibco, Waltham, MA), 0.5 to 1 mg/mL hyaluronidase (Sigma), and 250 μL of DNAse I, shaking gently at 37°C for a maximum of 30 minutes. Digestion was monitored and tissue was further digested mechanically by pipetting. Digested tissue was passed through a 100-μm filter, washed with phospate-buffered saline (PBS) containing 1 mM EDTA and 0.5% bovine serum albumin (BSA), and incubated with ACK lysing buffer (Gibco) to remove red blood cells. Single cell suspensions were incubated on ice with 4′,6-diamidino-2-phenylindole (DAPI) (molecular probes, Eugene, OR; 1:1000), used to exclude dead cells. EYFP+ cells were fluorescence-activated cell sorting (FACS) purified at the Salk Institute’s Flow Cytometry core facility on a BD Biosciences (San Jose, CA) Influx and an Aria Fusion cell sorter (100-μm size nozzle, 1 × PBS sheath buffer with sheath pressure set to 20 PSI). Cells were sorted in 1-drop Single Cell sort mode for counting accuracy.

       Single Cell Preparation and 10 × 3ʹ Library Construction

      FACS-isolated EYFP+ cells were loaded onto a microfluidic chip with barcoded beads according to 10X Chromium Next GEM Single Cell 3ʹ (v3.1, Catalog # 1000128; 10x Genomics Inc., Pleasanton, CA). Subsequent cell lysis, first-strand complementary DNA (cDNA) synthesis, and amplification were carried out according to the 10X v3.1 protocol, with cDNA amplification set for 11 cycles. Following sample indexing and bead-based library purification, the final library size distribution and concentration were measured using TapeStation (Agilent Biosystems, Santa Clara, CA) and Qubit (ThermoFisher, Waltham, MA). The libraries were pooled in equal molar ratio and sequenced on Illumina (San Diego, CA) NovaSeq 6000.

       10X Single Cell Data Processing

      The raw sequencing data from each mouse was first processed separately in Cell Ranger 4.0.0 using a customized reference based on refdata-gex-mm10-2020-A-R26 to allow quantification of EYFP expression. The EYFP reporter transgene was added to the refdata-gex-mm10-2020-A-R26 reference and rebuilt by running cellranger mkref with default parameters (10x Genomics). Each dataset captured similar numbers of cells (∼9k cells/mouse) with comparable sequencing depth and quality (see library QC summary in Supplementary Table 1). The filtered gene-by-cell count matrices from 10x Cellranger were further QCed and analyzed in R package Seurat (v3.2.3).
      • Stuart T.
      • Butler A.
      • Hoffman P.
      • et al.
      Comprehensive integration of single-cell data.
      Each library was processed separately in Seurat with NormalizeData(normalization.method = "LogNormalize") function. The normalized data were further linear transformed by ScaleData() function before dimension reduction. Principal components analysis (PCA) was performed on the scaled data by only using the most variable 2000 genes (identified using the default “vst” method). Cells were examined across all clusters to determine the low-quality cell threshold that accommodates the variation between cell types. Low-quality cells were removed with the same filtering parameters across libraries (percent.mt <=6 & nCount_RNA >= 1000). The QCed datasets were reprocessed as described previously, and clusters were labeled with cell types based on marker gene expression and previously published scRNA-seq derived gene signatures generated in the pancreas. Doublets were marked using the doubetFinder R package (largely consistent with the results from Scrublet package) and removed. In addition to the doublets labeled by douletFinder, 2 small clusters of cells with a high percentage of inferred doublets label, high number of UMI counts, and shared expression of marker genes were also removed as doublets. QCed datasets were integrated by using the fastMNN algorithm included in the SeuratWrapper package to minimize the technical batch variations.
      • Haghverdi L.
      • Lun A.T.L.
      • Morgan M.D.
      • et al.
      Batch effects in single-cell RNA-sequencing data are corrected by matching mutual nearest neighbors.
      A union of the most variable 1500 genes from each sample was used as integration features for fastMNN. The first 45 MNN principal components were used in UMAP embedding and cluster identification. New clusters were re-labeled by checking both known individual marker gene expression and the module scores of previously generated cell-type-specific gene signatures from a pancreas scRNA-seq study.
      • Elyada E.
      • Bolisetty M.
      • Laise P.
      • et al.
      Cross-species single-cell analysis of pancreatic ductal adenocarcinoma reveals antigen-presenting cancer-associated fibroblasts.
      EYFP-negative clusters of stromal contamination were removed and the epithelial population was further filtered for EYFP transcripts resulting in 13,362 EYFP+ cells (detectable EYFP transcripts) for downstream analysis.
      Tuft and enteroendocrine cells were subsequently subset, re-scaled, and re-clustered to better reveal cellular heterogeneity. A small cluster of 107 cells became evident as doublets after re-clustering and was removed from subsequent analysis. This cluster expresses tuft, acinar, and mucin/Ductal cells and has a higher number of genes and UMI counts detected per “cell.” Differentially expressed genes were identified by using FindAllMarkers() function with logfc.threshold=0.5 and min.pct=0.2.
      The processed human pancreas sNuc-seq dataset in seurat format from Tosti et al
      • Tosti L.
      • Hang Y.
      • Debnath O.
      • et al.
      Single nucleus and in situ RNA sequencing reveals cell topographies in the human pancreas.
      was obtained from http://singlecell.charite.de/pancreas/. Gene signature enrichment was assessed by seurat geneset module score.

       GO Enrichment Analysis

      GO analysis was done with R package topGO v2.42.0 (algorithm = "elim", statistic = "fisher") ( https://bioconductor.org/packages/release/bioc/html/topGO.html). Genes meeting the filtering threshold in DEGs analysis (pct.min=0.2, min.diff.pct=-inf, min.cell.feature =3, logfc.threshold > 0) were used as allGenes to create the topGOdata object. The top 20 enriched GO terms were plotted in ggplot2. Enrichment scores were calculated by taking -log10(P-value of Fisher’s exact test). In generating heatmaps of selected genes, cells were grouped by clusters and cells within a cluster group were re-ordered by increasing pseudotime within each cluster to better reveal expression pattern changes.

       RNA Velocity Estimation

      RNA velocity estimation was performed separately in each sample using the Velocyto.R package (v0.6).
      • La Manno G.
      • Soldatov R.
      • Zeisel A.
      • et al.
      RNA velocity of single cells.
      The spliced/unspliced/ambiguous reads classification was done by running the velocyto run10x pipeline on cellranger output. The resulting loom file was converted to Seurat object and subset to keep the 13362 EYFP+ cells. RNA velocity estimation was performed by using RunVelocity() function in Velocyto.R package with the following parameters: deltaT=1, ambiguous="ambiguous", kCells=25, fit.quantile=0.02,reduction = "mnn", ncores=30. Cell-cell distance used in the Runvelocity() step was calculated from the fastMNN corrected principal components space. Calculated velocity was visualized on the UMAP embedding through the show.velocity.on.embedding.cor() function with the parameters setting: min.grid.cell.mass=0.5, show.grid.flow=TRUE, grid.n=40, do.par=T, n.cores=30,cell.border.alpha=0.1, arrow.scale=1.

       CellRank Analysis

      The aggregate fate probabilities of tuft+eec lineage were calculated using the CellRank package (v.1.2.0) (doi: https://doi.org/10.1101/2020.10.19.345983). Briefly, the count matrix from the RNA assay in the Seurat object of Tuft+EEC data subset was converted into h5ad format using the SeuratDisk package and subsequently loaded using scranpy for CellRank analysis. New neighborhood graph was recomputed in scranpy using the first 30 fastMNN principal components. To get cell-cell transition probabilities, we used a PalantirKenel on monocle3 pseudotime. Nine macrostates corresponding to the 9 seurat clusters were identified using the compute_macrostate() function in the GPCCA estimator of the CellRank package. Macrostates corresponding to the Seurat clusters 0, 4, 7, and 8 were labeled as terminal states and the fate probabilities toward the tuft terminal state (seurat cluster 0) was calculated using the pl.cluster_fates()function. Similarly, the regulon activity score matrix was imported into scranpy for plotting regulon activity score over monole3 pseudotime. The plot was generated using the pl.gene.trends(n_test_point=500) function by fitting a GAM model between regulon activity score and monocle3 pseudotime in the cellrank package.

       Regulon Analysis

      Gene regulatory network (Regulon) analysis in this study was performed using the SCENIC R package (v1.2.0).
      • Aibar S.
      • Gonzalez-Blas C.B.
      • Moerman T.
      • et al.
      SCENIC: single-cell regulatory network inference and clustering.
      The mouse mm10 500bp_up_100bp_down and 10kb_up_and_down databases were downloaded from the cistarget database; the TSS database was downloaded from the Aerts lab Web site (https://resources.aertslab.org/cistarget/) and used in this analysis. The filtered count matrix of the 13,362 EYFP+ cells were further filtered to exclude non-expressed genes using the function geneFiltering(exprMat, scenicOptions=scenicOptions, minCountsPerGene=100, minSamples=25) in the SCENIC package. Spearman correlation between each TF and other genes was run by runCorrelation() function with the default settings. TF-target gene co-expression module was identified in pySCENIC by running the Run_arboreto_with_multiprocessing.py pipeline with default setting and random seed 915. The adjacency file was imported into SCENIC R package and the normalized weight was used as weight subsequent steps. Regulon and regulon activity was calculated by running the runSCENIC_1_coexNetwork2modules(), runSCENIC_2_createRegulons(),runSCENIC_3_scoreCells() functions with default settings. Filtered Count matrix was used in step 3 for regulon activity scoring using the AUCell algorism. The resulting regulon activity score matrix was subsequently imported into Seurat through CreateAssayObject() function. Regulon activity score was plotted in Seurat without further normalization or scaling.

       Monocole3 Trajectory Analysis

      Seurat object of the 13362 EYFP+ cells was converted into Monocle3 CellDataSet object by using the as.cell_data_set()function in the SeuratWrapper Package. Subsequent trajectory analysis is done in the monocle3 package.
      • Cao J.
      • Spielmann M.
      • Qiu X.
      • et al.
      The single-cell transcriptional landscape of mammalian organogenesis.
      Cells were re-clustered in monocle by running the cluster_cells() function with setting resolution= 2.5e-4. Then cells were ordered in pseudotime by running the order_cells() function with default setting. Root state was manually set to acinar cells. Trajectory plot was generated using the plot_cells() function. Tuft+EEC monocle3 trajectories were done similarly as described previously; 2586 Tuft+EEC cells subset was used and the cluster_cells() function was run with resolution = 1e-3.The root state was manually set to Cluster 6 based on the RNA velocity information.

       p-Creode Trajectory Analysis

      Count matrices were first filtered using the dropkick algorithm to identify cells vs ambient barcodes.
      • Heiser C.N.
      • Wang V.M.
      • Chen B.
      • et al.
      Automated quality control and cell identification of droplet-based single-cell data using dropkick.
      Counts data were then normalized to the median total read count across all retained barcodes, transformed using an inverse hyperbolic sine function, and Z-score standardized through functions available in Scanpy (v1.6.1).
      • Wolf F.A.
      • Angerer P.
      • Theis F.J.
      SCANPY: large-scale single-cell gene expression data analysis.
      p-Creode (v2.2.0) was performed on normalized data in supervised mode with a max of 10 principal components across 1000 scored bootstrapped runs. The highest scoring run, producing the most representative graph from the ensemble, was chosen to visualize overlays and pseudotime gene expression dynamics. These gene expression dynamics were smoothened and plotted using a generalized additive model, implemented through pyGAM (v0.8.0, zenodo.org). p-Creode density parameters were chosen according to automatic thresholding, or the best_guess function. Noise and target density parameters varied between the different datasets analyzed, resulting in the use of approximately 20% to 25% of cells in each bootstrapped replicate.

       Data Integration of Injury 10x scRNA-seq Data and Published Datasets

      To avoid artifacts in differential gene analysis between our dataset and previous published datasets, fastq files were downloaded from the SRA using accession number GSE141007 (Schlesinger et al
      • Schlesinger Y.
      • Yosefov-Levi O.
      • Kolodkin-Gal D.
      • et al.
      Single-cell transcriptomes of pancreatic preinvasive lesions and cancer reveal acinar metaplastic cells' heterogeneity.
      ), GSE159343 (Hendley et al
      • Hendley A.M.
      • Rao A.A.
      • Leonhardt L.
      • et al.
      Single-cell transcriptome analysis defines heterogeneity of the murine pancreatic ductal tree.
      ), and subsequently reprocessed using the same version of mm10 reference genome annotation (refdata-gex-2020-A reference genome from 10x Genomics) as our injury dataset. The tdTomato transgene was added to the reference genome to allow for quantification. The KrasG12D datasets were QCed using similar thresholds as were used in the original study (minimum of 200 genes and 1000 UMIs per cell) except with a slightly higher cutoff threshold for percentage of mitochondrial reads (<=7.5%). The KrasG12D dataset and our Injury dataset were merged in seurat v3 and further integrated using the fastMNN method. For integration of our injury dataset and the sorted 6-month KrasG12D dataset, top 2000 variable genes identified using the “vst” method on the merged Injury+KrasG12D dataset were used as the integration feature for fastMNN. The first 30 PCs returned by the fastMNN algorism were used for 2-dimensional UMAP embedding and Cluster identification. The normal pancreas scRNA-seq dataset from Hendley et al
      • Hendley A.M.
      • Rao A.A.
      • Leonhardt L.
      • et al.
      Single-cell transcriptome analysis defines heterogeneity of the murine pancreatic ductal tree.
      was QCed with percent.mt <=10 & nCount_RNA>=2000 & nCount_RNA < 30000. Cells were clustered and annotated using classic markers.
      For comparison between injury and normal acinar and ductal cells, populations of interest were extracted, log normalized in Seurat, and then subsequently integrated with our injury-induced ADM dataset. In addition, one cluster of EYFP-negative ductal cells was identified. This population was filtered out from our analysis of EYFP+ injury-induced ADM populations. In this integrated injury vs normal comparative analysis, we also included these 305 EYFP-negative ductal cells. Merged data from 4 different sources were integrated using the fast Mutual Nearest Neighbor (fastMNN) batch correction method and the first 40 “MNN” components were further embedded in 2 dimensions using the Uniform Manifold Approximation and Projection (UMAP) algorithm.
      • Tosti L.
      • Hang Y.
      • Debnath O.
      • et al.
      Single nucleus and in situ RNA sequencing reveals cell topographies in the human pancreas.
      Cells were clustered using the Seurat FindNeighbors(dims=1:40, reduction=”mnn”) and FindClusters(resolution=0.1). Differentially expressed gene analysis was performed on the integrated dataset between the injury group and normal group using the FindMarkers() function in Seurat with logfc.threshold=0.4, min.pct=0.3 used for the ADM-MucinDuctal vs EYFP-negative ductal cell comparison and logfc.threshold=0.4, min.pct=0.2 used for the injury acinar vs normal cells from (Hosein et al
      • Hosein A.N.
      • Huang H.
      • Wang Z.
      • et al.
      Cellular heterogeneity during mouse pancreatic ductal adenocarcinoma progression at single-cell resolution.
      and Schlesingler et al
      • Schlesinger Y.
      • Yosefov-Levi O.
      • Kolodkin-Gal D.
      • et al.
      Single-cell transcriptomes of pancreatic preinvasive lesions and cancer reveal acinar metaplastic cells' heterogeneity.
      ) comparison. DEGs were visualized in volcano plot using the EnhancedVocano package (https://github.com/kevinblighe/EnhancedVolcano). Top 100 up-regulated and top 100 down-regulated genes (only 36 injury acinar vs normal acinar down-regulated genes) were used to generate heatmaps. To get visually comparable patterns in the heatmap from the KRasG12D sample, the top DEGs were re-ordered based on its scaled average expression in the KRasG12D-6M sample prior to running the DoHeatmap()function in Seurat.
      To calculate the percentage of digestive enzyme reads, gene symbols of the human pancreatic digestive enzyme genes were obtained from Tosti et al
      • Tosti L.
      • Hang Y.
      • Debnath O.
      • et al.
      Single nucleus and in situ RNA sequencing reveals cell topographies in the human pancreas.
      and converted to mouse gene symbols using the biomaRt package. The total number of UMIs that belong to the digestive enzyme genes were divided by the total number of reads in each cell to obtain the percent of digestive enzyme reads. To compare the number of expressed genes per cell, we downsampled the number of UMIs per cell to 3k to avoid artifacts introduced by different sequencing depths between different datasets. The number of genes detected in the down-sampled acinar cells were examined and compared visually using a violin plot. We also want to point out that the number of expressed genes detected may be influenced by other technical factors, such as potential saturation effects from a shift in zymogens between the different conditions, which was not considered in this comparison.
      The datasets used for these comparisons were generated from different strains of mice (with different technical variations) and care should be taken when assigning a nonvalidated marker to a disease state rather expression being due to mouse strain or batch effects.

       Histological Staining

      Tissues were fixed overnight in zinc-containing, neutral-buffered formalin (Fisher Scientific, Waltham, MA), embedded in paraffin, cut in 5-μm sections, mounted, and stained. Sections were deparaffinized in xylene, rehydrated in ethanol, and then washed in PBST and PBS. Endogenous peroxidase activity was blocked with a 1:50 solution of 30% H2O2:PBS followed by microwave antigen retrieval in 100 mM sodium citrate, pH 6.0. Sections were blocked with 1% BSA and 5% normal goat or rabbit serum in 10 mM Tris (pH 7.4), 100 mM MgCl2, and 0.5% Tween-20 for 1 hour at room temperature. Primary antibodies were diluted in blocking solution and incubated overnight. Anti-Gkn3 was generated by the Mills Laboratory (Washington University, St. Louis, MO), anti-Tff2 was a generous gift from Sir Nicholas Wright (Barts Centre, London, UK), and anti-Gif was a generous gift from the Alberts Laboratory (Washington University, St. Louis, MO). Information on additional primary antibodies is provided in Supplementary Table 8. Slides were then washed, incubated in streptavidin-conjugated secondaries (for rabbit or mouse antibodies, Abcam [Cambridge, UK], for rat or goat antibodies, Vector [Olean, NY]) and developed with DAB substrate (Vector). Hematoxylin and eosin (H&E) staining was performed to assess tissue morphology. All slides were scanned and imaged on an Olympus VS-120 or Olympus VS-200 Virtual Slide Scanning microscope.

       Fluorescence Microscopy

      Tissues were fixed overnight in 4% paraformaldehyde, washed 3 times with PBS, and floated overnight in 30% sucrose. Tissues were then incubated in a 1:1 mixture of 30% sucrose and Tissue-Tek optimal cutting temperature compound (OCT, VWR) for 30 minutes, embedded in OCT and frozen at −80°C; 7-μm tissue sections were cut, permeabilized with 0.1% Triton X-100 in 10 mM PBS, and blocked with 5% normal donkey serum and 1% BSA in 10 mM PBS for 1 hour at room temperature. Tissue sections were stained with primary antibodies in 10 mM PBS supplemented with 1% BSA and 0.1% Triton X-100 overnight (Supplementary Table 8). Sections were then washed 3 x 15 minutes in PBS with 1% Triton X-100, incubated in Alexa Fluor secondary antibodies and/or phalloidin (Invitrogen, Carlsbad, CA), washed again for 3 x 5 minutes, rinsed with distilled water, and mounted with Prolong Gold containing DAPI (Invitrogen). Immunofluorescence on paraffin-embedded tissues followed the IHC protocol until the blocking step. Instead, tissues were blocked in the donkey serum block described previously and then followed the protocol for fluorescence microscopy described here. All slides were imaged on an Olympus VS-120, an Olympus VS-200 Virtual Slide Scanning microscope, or Zeiss 880 Airyscan Super-Resolution microscope.
      For multiplex staining featured in Figure 2B, primary antibodies (Supplementary Table 8) were used at the listed concentrations and incubated in tris buffered saline (TBS) overnight. The next day, slides were washed 4 x 5 minutes in TBS and incubated for 1 hour at room temperature with DAPI (H-1000; Invitrogen) and secondary antibodies raised in donkey or goat and conjugated to AlexaFluor 488, 555, 594 or 647 (Invitrogen or JacksonImmunoresearch) were used at 1:400 dilution. After incubation, samples were washed 4 x 5 minutes in TBS, covered with Vectashield (Invitrogen), and cover-slipped. Multi-plant confocal images were acquired using an upright Leica Stellaris X5 microscope with a ×20/0.75 N.A. multi-immersion objective. After acquisition, images were processed using ImageJ/Fiji (National Institutes of Health, Bethesda, MD) and denoised using a median filter of 1px. Images shown are maximum projections of confocal stacks.

       instant Structured Illumination Microscopy

      Image-stacks of Brainbow lineage tracing were captured with a system assembled by BioVision Technologies (Exton, PA) consisting of a VisiTech International (Sunderland, United Kingdom) instant Structured Illumination Microscopy (iSIM) module attached to a Nikon (Melville, NY) Ti2 stand equipped with a Nikon 100x 1.49 Plan Apo objective. Raw image stacks were then deconvolved using a Microvolution (Cupertino, CA) plugin in Fiji (ImageJ).

       3DEM Methods

      Samples of injured pancreas tissue were processed and imaged as described previously with some modifications.
      • Horstmann H.
      • Korber C.
      • Satzler K.
      • et al.
      Serial section scanning electron microscopy (S3EM) on silicon wafers for ultra-structural volume imaging of cells and tissues.
      Materials were sourced from Electron Microscopy Sciences (Hatfield, PA) unless noted otherwise. Briefly, fresh biopsies of injured tissues were dissected into small pieces (<1 mm) and quickly immersed in 37°C buffered fixative (2.5% glutaraldehyde, 2% paraformaldehyde, 3 mM CaCl2, 0.1M sodium cacodylate) before overnight storage at 4°C in the same solution.
      • Seligman A.M.
      • Wasserkrug H.L.
      • Hanker J.S.
      A new staining method (OTO) for enhancing contrast of lipid--containing membranes and droplets in osmium tetroxide--fixed tissue with osmiophilic thiocarbohydrazide (TCH).
      Samples were serially rinsed with ice-cold cacodylate/CaCl2 buffer and stained with buffered reduced osmium (1.5% Osmium, 1.5% potassium ferrocyanide, 3 mM CaCl2, 0.1M sodium cacodylate) for 45 minutes in the dark at room temperature. Samples were serially rinsed with ice-cold distilled water before 30 minutes of treatment with filtered 1% aqueous thiocarbohydrazide (Ted Pella, Redding, CA) that was prepared at 60°C for an hour before use. Samples were serially rinsed with ice-cold distilled water before staining with 1.5% aqueous osmium tetroxide for 45 minutes in the dark at room temperature. Samples were rinsed serially with ice-cold water before staining overnight with 1% uranyl acetate at 4°C. Samples were rinsed with distilled water and stained with Walton’s lead aspartate for 30 minutes at 60°C.
      • Walton J.
      Lead aspartate, an en bloc contrast stain particularly useful for ultrastructural enzymology.
      Samples were rinsed again with distilled water and serially dehydrated in ice cold ethanol. Samples were then infiltrated with Durcupan resin (hard formulation) and polymerized at 60°C for 3 days.
      A sample block was prepared for serial sectioning as previously described and a ribbon of 248 serial sections, each of dimension ∼150 μm × 400 μm × 70 nm (x-y-z), were collected onto a silicon wafer diced to a dimension of 35 × 6 mm using Diatome diamond knives mounted in a Leica UC8 ultramicrotome.
      • Harris K.M.
      • Perry E.
      • Bourne J.
      • et al.
      Uniform serial sectioning for transmission electron microscopy.
      ,
      • Burel A.
      • Lavault M.T.
      • Chevalier C.
      • et al.
      A targeted 3D EM and correlative microscopy method using SEM array tomography.
      The chip was imaged using a Zeiss Sigma VP scanning electron microscope equipped with a Gatan backscattered electron detector and Atlas5 (Fibics) control software. Multi-resolution maps were generated for a subset of sections throughout the ribbon, of which images were used for Figure 2, and a region of interest (ROI) for 3DEM was identified. Images were acquired with an accelerating voltage of 3 kV, a 20-μm aperture, at a working distance of 9 mm, and with pixel size of 8 nm. Contiguous images of the ROI were assembled into a series of 8-bit .tiff and aligned.
      A delta cell, a tuft cell, 2 mucinous cells, an acinar cell, and their associated basement membrane were identified by their salient ultrastructural features. Their boundaries, along with their respective nuclei, were traced using VAST Lite annotation software.
      • Berger D.R.
      • Seung H.S.
      • Lichtman J.W.
      VAST (Volume Annotation and Segmentation Tool): efficient manual and semi-automatic labeling of large 3D image stacks.
      Cells and nuclei were exported using the Matlab tools bundled with VAST Lite as .obj files, which were imported into Blender 2.79 (Blender Foundation, Amsterdam, The Netherlands) with add-ons BlendGAMer and Neuromorph both installed.
      • Lee C.T.
      • Laughlin J.G.
      • Angliviel de La Beaumelle N.
      • et al.
      3D mesh processing using GAMer 2 to enable reaction-diffusion simulations in realistic cellular geometries.
      ,
      • Jorstad A.
      • Blanc J.
      • Knott G.
      NeuroMorph: a software toolset for 3D analysis of neurite morphology and connectivity.
      Cells and nuclei were processed with GAMer to optimize geometry for visualization, before superimposition with 3DEM data using Neuromorph and animation and rendering using Blender. For video figures, rendered scenes were assembled into .avi files using ImageJ, post-processed using After Effects (Adobe, San Jose, CA), and compressed using Handbrake (Handbrake.fr).
      • Schneider C.A.
      • Rasband W.S.
      • Eliceiri K.W.
      NIH Image to ImageJ: 25 years of image analysis.

       Transmission Electron Microscopy of Murine Pancreas

      Tissues were fixed in a solution of 2% paraformaldehyde, 2.5% glutaraldehyde, and 2 mM CaCl2 in 0.15 M sodium cacodylate buffer (pH 7.4) for 2 hours at room temperature. They were then post-fixed in 1% osmium tetroxide for 40 minutes and 1.5% potassium ferricyanide in sodium cacodylate buffer for 1 hour at 4°C in the dark. Tissues were stained en bloc in 1% aqueous uranyl acetate (4C in the dark), dehydrated in a series of graded ethanols, and embedded in Eponate12 resin (Ted Pella). Ultra-thin sections (70 nm) were obtained using a diamond knife (Diatome) in an ultramicrotome (Leica EM UC7) and placed on copper grids (300 mesh). Sections were imaged on a Zeiss Libra 120 transmission electron microscope (TEM) operated at 120 kV.

       Generation of SPEM and TEM of Murine Stomach

      Gastric SPEM was induced with tamoxifen (5 mg/20 g mouse body weight; Toronto Research Chemicals Inc., Toronto, Canada) administered by intraperitoneal injection for three consecutive days.
      • Saenz J.B.
      • Burclaff J.
      • Mills J.C.
      Modeling murine gastric metaplasia through tamoxifen-induced acute parietal cell loss.
      Tamoxifen was prepared by sonication in a 90% sunflower seed oil (Sigma) and 10% ethanol solution. For TEM, stomachs were prepared and imaged as previously described.
      • Ramsey V.G.
      • Doherty J.M.
      • Chen C.C.
      • et al.
      The maturation of mucus-secreting gastric epithelial progenitors into digestive-enzyme secreting zymogenic cells requires Mist1.
      Briefly, stomach corpus was collected as described, fixed overnight at 4˚C in fixative (modified Karnovsky’s), and sectioned. Tissue was processed by the Washington University in St. Louis Department of Pathology and Immunology Electron Microscopy Facility.
      All images (IHC, IF, EM) were digitally enhanced to edit the color, brightness, and contrast levels using Zen (Carl Zeiss), ImageJ/Fiji (NIH), and/or Photoshop (Adobe) software.

       Figures

      Models shown in Figures 1A and in the graphical abstract were generated using BioRender.com.
      Figure thumbnail fx2
      Supplementary Figure 1Isolation and analysis of lineage-traced EYFP+ cell from injured pancreata. (A) H&E analysis of the normal mouse pancreas or pancreata treated with caerulein for either 2 or 4 weeks. Scale bar, 50 μm. (B) FACS strategy for isolation of lineage traced EYFP+ cells from injured CERTY pancreata. EYFP-high cells included in the red circle were collected for sequencing. (C) UMAPs showing annotated clusters from FACS collected EYFP+ cells from the injured pancreata of 4 CERTY mice treated with caerulein for either 2 or 4 weeks.
      Figure thumbnail fx3
      Supplementary Figure 2Annotation of lineage-traced EYFP+ cells. (A) Expression of EYFP and (B) epithelial cell marker Krt8 overlaid on the UMAP of all FACS collected EYFP+ cells. (C) UMAP of EYFP+ cells from B annotated with Seurat-identified clusters. (D) Expression of acinar cell (Prss1, Try4, and Cpa1), mucin/ductal cell (Muc6, Spp1, Car2, and Krt19), tuft/EEC progenitor cell (Sox4), tuft cell (Pou2f3, Trpm5, and Siglecf), EEC (Neurog3, Neurod1, Pax6, and Chga), and proliferation marker Mki67 on the UMAP of EYFP+ cells shown in (C). Color intensity in (A–C) indicates the normalized gene expression level for a given gene in each cell.
      Figure thumbnail fx4
      Supplementary Figure 3Comparison of normal and ADM-associated acinar and ductal cells. UMAP plot of integrated datasets colored by Seurat clusters in (A) and major cell type and data source in (B). Margins of clusters that correspond to the acinar cells from injury, normal acinar cells, and normal ductal cells, are marked with dashed curved lines. (C) A separate view of the major cell types in (B) split by data source.
      Figure thumbnail fx5
      Supplementary Figure 4Comparison of normal and ADM-associated acinar and ductal cells. (A) Volcano plot of the differentially expressed genes between injury acinar cells and normal acinar cells (Schlesinger et al,
      • Schlesinger Y.
      • Yosefov-Levi O.
      • Kolodkin-Gal D.
      • et al.
      Single-cell transcriptomes of pancreatic preinvasive lesions and cancer reveal acinar metaplastic cells' heterogeneity.
      NormalAcinar.YS and Hosein et al,
      • Hosein A.N.
      • Huang H.
      • Wang Z.
      • et al.
      Cellular heterogeneity during mouse pancreatic ductal adenocarcinoma progression at single-cell resolution.
      NormalAcinar.ANH). (B) Heat-map of the top 100 up-regulated and all 34 down-regulated genes in (A). The y-axis represents DEGs that were reordered based on increasing average level of scaled expression in KrasG12D_6M. The x-axis represents acinar cells grouped by sample. (C) Barplot of the average percentage of transcripts from digestive enzymes in individual cells. Error bars represent the standard error of the mean. (D) Violin plots of the number of UMI per cell and the number of genes detected per cell across down-sampled (max nUM I= 3k) samples. (E) Volcano plot of the differentially expressed genes between EYFP-positive ADM-induced mucin/ductal cells and EYFP-negative ductal cells. (F) Heatmap of the top 200 differentially expressed genes in (E). The y-axis represents DEGs that were reordered based on increasing average level of scaled expression in KrasG12D_6M. The x-axis represents ADM-induced mucin/ductal cells, EYFP-negative ductal cells, ADM cells from KrasG12D_6M sample (Schlesinger et al
      • Schlesinger Y.
      • Yosefov-Levi O.
      • Kolodkin-Gal D.
      • et al.
      Single-cell transcriptomes of pancreatic preinvasive lesions and cancer reveal acinar metaplastic cells' heterogeneity.
      ) and ductal cells grouped by sample (Hendley et al,
      • Hendley A.M.
      • Rao A.A.
      • Leonhardt L.
      • et al.
      Single-cell transcriptome analysis defines heterogeneity of the murine pancreatic ductal tree.
      NormalDuctal.AMH and Schlesinger et al,
      • Schlesinger Y.
      • Yosefov-Levi O.
      • Kolodkin-Gal D.
      • et al.
      Single-cell transcriptomes of pancreatic preinvasive lesions and cancer reveal acinar metaplastic cells' heterogeneity.
      NormalDuctal.YS).
      Figure thumbnail fx6
      Supplementary Figure 5Acinar cell lineage tracing in normal and injured CERTY pancreata. Co-IF for EYFP (yellow), ductal cell marker PanCK (cyan), and membrane and tuft cell marker phalloidin (magenta) in normal or injured CERTY pancreata. Tuft cells highlighted in insert. DAPI, blue. Scale bars, 100 μm.
      Figure thumbnail fx7
      Supplementary Figure 6Lineage progression of ADM-derived populations. (A) Expression of acinar cell marker Cpa1, mucin/ductal marker Muc6, tuft cell marker Pou2f3, and EEC marker Chga overlaid on the pCreode plot shown in C. (B) Heat map with hierarchical clustering showing regulons enriched in EYFP+ populations. (C) Expression of key transcription factors overlaid on the UMAP from B. Color intensity indicates the relative enrichment of signature gene sets in each cell. (D) Expression of Brainbow reporters (YFP, green, cytoplasmic; GFP, green, nuclear; RFP, red, cytoplasmic) in injured pancreata from Ptf1aCreERTM;RosaBrainbow.2/+ mice treated with 10 or (E) 5 doses of tamoxifen. Scale bars, 50 μm and 20 μm, respectively. Phalloidin staining in (E) highlights a tuft cell. DAPI, blue.
      Figure thumbnail fx8
      Supplementary Figure 7ADM as a pyloric-type metaplasia. (A) Expression of spasmolytic polypeptide-expressing metaplasia (SPEM) markers Gkn3, Dmbt1, and Wfdc2 and (B) chief cell markers Pga5 and Pgc overlaid on the UMAP of EYFP+ cells from Figure 1B. Color intensity indicates the normalized gene expression level for a given gene in each cell. (C) IHC for SPEM markers CD44v9, GKN3, AQP5, and WFDC2, and chief cell marker GIF in the normal stomach or pancreas in injury- or oncogene- (KrasG12D;Ptf1aCre/+) induced ADM and PanIN. I, islet. Scale bar, 50 μm.
      Figure thumbnail fx9
      Supplementary Figure 8Comparison of ADM to normal and injured gastric epithelium. (A) Expression of normal gastric epithelial cell type signatures derived from Bockerstett et al,
      • Bockerstett K.A.
      • Lewis S.A.
      • Wolf K.J.
      • et al.
      Single-cell transcriptional analyses of spasmolytic polypeptide-expressing metaplasia arising from acute drug injury and chronic inflammation in the stomach.
      and shown in D overlaid on the UMAP of EYFP+ cells from B. (B) UMAP of scRNA-seq data collected from normal stomach and a SPEM forming model of chronic gastric injury, derived from Bockerstett et al.
      • Bockerstett K.A.
      • Lewis S.A.
      • Wolf K.J.
      • et al.
      Single-cell transcriptional analyses of spasmolytic polypeptide-expressing metaplasia arising from acute drug injury and chronic inflammation in the stomach.
      (C) Expression of the chronic injury-induced gastric SPEM signature overlaid on the UMAP of EYFP+ cells from B. (D) Regulon scores (transcription factor activity) overlaid on the UMAP from Figure 1B. (E) Gene module scores for TFs shown in (D) for normal gastric chief cells, acute SPEM cells, or (F) SPEM populations associated with chronic gastric injury. Color intensity in (A), (C), and (D) indicates the normalized gene expression level for a given gene in each cell. EC, enterochromaffin cell.
      Figure thumbnail fx10
      Supplementary Figure 9ADM results in substantial enteroendocrine cell heterogeneity. (A) Expression of EEC markers overlaid on the UMAP of tuft cells and EECs from A. (B) Expression of islet-derived hormones insulin (Ins1 and Ins2) and glucagon (Gcg) overlaid on the UMAP from A. (C) Heat map with hierarchical clustering showing expression of Notch pathway inhibitor genes enriched in EEC progenitor cells. (D) Heat map with hierarchical clustering showing key TFs expressed in Tuft/EEC progenitors, tuft cells, and EECs. (E) Expression of TFs enriched in different EEC populations, overlaid on the UMAP from A. Color intensity indicates the normalized gene expression level for a given gene in each cell. (F) IHC for gastrin in the antrum of the normal stomach, duodenum, or pancreas and in the injured pancreas. Scale bar, 25 mm. (G) Lineage tracing and Co-IF for EYFP, green, hormones (Ghrelin, 5-HT, SST, PPY, magenta), pan-CK or phalloidin (white), or DAPI, blue. Arrows indicate positive cells. Scale bars, 10 μm.
      Figure thumbnail fx11
      Supplementary Figure 10Tuft cell heterogeneity as a reflection of differentiation status. (A) Barplots of the enrichment scores of the top 20 GO terms associated with the DEGs in the 4 tuft cell clusters. The x-axis represents the enrichment score calculated by -log10 (P value of Fisher's exact test). (B) Nine macrostates computed from Schur decomposition. The Seurat cluster labels were added to the associated macrostates and colored according to the Seurat cluster colors in A. (C) Fate probability scores towards the terminal tuft cell lineage overlaid on the Tuft+EEC UMAP. (D) Seurat geneset module scores of the neuronal and inflammatory signatures in murine small intestinal tuft cells overlaid on the Tuft+EEC UMAP. (E–G) Expression of select genes overlaid on the Tuft+EEC UMAP, labeled by Seurat cluster.
      Figure thumbnail fx12
      Supplementary Figure 11Comparison of injury and oncogene-induced ADM. (A) UMAP showing annotated clusters of scRNA-seq data generated from the pancreas of a 6-month-old KCT mouse. (B) FastMNN integrated EYFP+ injury-induced ADM and tdTomato+ cells from a 6-month-old KCT mouse, annotated by sample or (C) Seurat cluster. (D) Expression of metaplasia-specific transcription factors Foxq1 and (E) Onecut2 overlaid on the UMAP from A. (F) Volcano plots showing the top 20 differentially expressed genes between the injury- and oncogene-induced acinar and (G) mucin/ductal clusters. Blue, KrasG12D; red, injury. (H) Expression of pre-tumor markers Marcksl1 and (I) Igfbp7 overlaid on the UMAP from A. Color intensity indicates the normalized gene expression level for a given gene in each cell.
      Figure thumbnail fx13
      Supplementary Figure 12KrasG12D induction in ADM populations using inducible Cre mice. (A) Expression of Hnf1b overlaid on the UMAP from A. (B) Lineage tracing and EYFP expression in HERT or KHERT mice treated with tamoxifen or KHERT mice first treated with caerulein to induce injury, followed by tamoxifen to induce KrasG12D expression in ADM. Scale bar, 50 μm. (C) Co-IF for EYFP (green) and ductal cell marker cytokeratin (pan-CK, red) in the pancreata of KHERT mice treated with caerulein then tamoxifen. DAPI, blue. Scale bar, 50 μm. (D) Tomato (red) and GFP (green) expression in the intestines and (E) pancreas of PCERT;MTMG mouse treated with tamoxifen. DAPI, blue. Scale bars, 50 μm. (F) Co-IF for GFP (green) and tuft cell markers POU2F3 or DCLK1 (both yellow) in the pancreas of a PCERT;MTMG mouse treated with caerulein to induce ADM and tuft cell formation, followed by tamoxifen to induce GFP expression in Pou2f3+ cells. DAPI, blue. Scale bars, 25 μm.
      Figure thumbnail fx14
      Supplementary Figure 13ADM cell type markers in human pancreatitis. (A) Violin plots of MUC5B+Ductal cell markers MUC5B and KRT19, SPEM markers MUC6, TFF2, WFDC2, and CD44, and mucin/ductal transcription factors FOXQ1 and CREB3L1 in sNuc-seq data from human pancreatitis. (B) Expression of the human MUC5B+Ductal signature overlaid on the mouse ADM dataset (left) and the murine mucin/ductal gene signature overlaid on the human chronic pancreatitis dataset (right). (C) Expression of the chronic SPEM signature shown in humanized and overlaid on the sNuc-seq data from normal pancreas (left) and chronic pancreatitis (right). Color intensity indicates the relative gene expression for a given cell type in each cluster. (D) Co-IF for SPEM markers AQP5 (magenta), and CD44v9 (cyan), in an area of relatively normal pancreas (left) and injured pancreas (right) from the same patient sample. Scale bar, 25 μm. (E) Co-IF for SPEM marker TFF2 (magenta), tuft cell marker phospho-Girdin (yellow), and g-actin (cyan) in tissue specimens from patients with chronic pancreatitis. Scale bars, 25 μm. (F) Violin plots of tuft cell markers and TFs (SOX4, POU2F3, SPIB, RUNX1) shown in , or (G) pre-tumor markers in sNuc-seq data from human pancreatitis.
      Supplementary Table 8Primary Antibodies Used in Immunofluorescence (IF) and IHC Studies
      AntibodyCompanyCatalog #Dilution, IFDilution, IHC
      AmylaseAbcamab21156N/A1:5000
      Aqp5SigmaHPA-0650081:5001:500
      Cd44v9 (anti-mouse)Cosmo Bio Co.LKG-M0021:10,0001:10,000
      Cd44v9 (anti-human)Cosmo Bio Co.LKG0M0031:10,000N/A
      ChgASanta Cruzsc-14881:100N/A
      Gamma actinSanta Cruzsc-656381:100N/A
      GastrinBiogenexPU019-UPN/A1:100
      GFPAbcamab6673N/A1:1000
      GhrelinCell Signaling318651:10001:2000
      GIFN/AN/A1:10001:1000
      Gkn3N/AN/AN/A1:150
      MSLN (anti-mouse)LS BioLS-C407883N/A1:2000
      Muc5acThermofisherMA5-12178N/A1:500
      Pan-CKAbcamab93771:100N/A
      Phospho-Girdin (pY1798)IBL

      America
      281431:100N/A
      PPYAbcamab2558271:80,0001:80,000
      SerotoninSanta Cruzsc-580311:2001:500
      SomatostatinMilliporeMAB3541:10001:2000
      SynaptophysinCell Marque336R-941:500
      TFF2N/AN/A1:500N/A
      Wfdc2InvitrogenPA5-80227N/A1:2000
      N/A, not applicable.

      References

        • Giroux V.
        • Rustgi A.K.
        Metaplasia: tissue injury adaptation and a precursor to the dysplasia-cancer sequence.
        Nat Rev Cancer. 2017; 17: 594-604
        • Storz P.
        Acinar cell plasticity and development of pancreatic ductal adenocarcinoma.
        Nat Rev Gastroenterol Hepatol. 2017; 14: 296-304
        • Willet S.G.
        • Lewis M.A.
        • Miao Z.F.
        • et al.
        Regenerative proliferation of differentiated cells by mTORC1-dependent paligenosis.
        EMBO J. 2018; 37e98311
        • Murtaugh L.C.
        • Keefe M.D.
        Regeneration and repair of the exocrine pancreas.
        Annu Rev Physiol. 2015; 77: 229-249
        • DelGiorno K.E.
        • Naeem R.F.
        • Fang L.
        • et al.
        Tuft cell formation reflects epithelial plasticity in pancreatic injury: implications for modeling human pancreatitis.
        Front Physiol. 2020; 11: 88
        • Delgiorno K.E.
        • Hall J.C.
        • Takeuchi K.K.
        • et al.
        Identification and manipulation of biliary metaplasia in pancreatic tumors.
        Gastroenterology. 2014; 146: 233-244.e5
        • DelGiorno K.E.
        • Chung C.Y.
        • Vavinskaya V.
        • et al.
        Tuft cells inhibit pancreatic tumorigenesis in mice by producing prostaglandin D2.
        Gastroenterology. 2020; 159: 1866-1881.e8
        • Tosti L.
        • Hang Y.
        • Debnath O.
        • et al.
        Single nucleus and in situ RNA sequencing reveals cell topographies in the human pancreas.
        Gastroenterology. 2020; 160: 1330-1344.e11
        • Sah R.P.
        • Garg S.K.
        • Dixit A.K.
        • et al.
        Endoplasmic reticulum stress is chronically activated in chronic pancreatitis.
        J Biol Chem. 2014; 289: 27551-27561
        • Elyada E.
        • Bolisetty M.
        • Laise P.
        • et al.
        Cross-species single-cell analysis of pancreatic ductal adenocarcinoma reveals antigen-presenting cancer-associated fibroblasts.
        Cancer Discov. 2019; 9: 1102-1123
        • Schlesinger Y.
        • Yosefov-Levi O.
        • Kolodkin-Gal D.
        • et al.
        Single-cell transcriptomes of pancreatic preinvasive lesions and cancer reveal acinar metaplastic cells' heterogeneity.
        Nat Commun. 2020; 11: 4516
        • Hosein A.N.
        • Huang H.
        • Wang Z.
        • et al.
        Cellular heterogeneity during mouse pancreatic ductal adenocarcinoma progression at single-cell resolution.
        JCI Insight. 2019; 5e129212
        • Pan F.C.
        • Bankaitis E.D.
        • Boyer D.
        • et al.
        Spatiotemporal patterns of multipotentiality in Ptf1a-expressing cells during pancreas organogenesis and injury-induced facultative restoration.
        Development. 2013; 140: 751-764
        • Sato A.
        Tuft cells.
        Anat Sci Int. 2007; 82: 187-199
        • Trapnell C.
        • Cacchiarelli D.
        • Grimsby J.
        • et al.
        The dynamics and regulators of cell fate decisions are revealed by pseudotemporal ordering of single cells.
        Nat Biotechnol. 2014; 32: 381-386
        • Herring C.A.
        • Banerjee A.
        • McKinley E.T.
        • et al.
        Unsupervised trajectory analysis of single-cell RNA-seq and imaging data reveals alternative tuft cell origins in the gut.
        Cell Syst. 2018; 6: 37-51.e9
        • La Manno G.
        • Soldatov R.
        • Zeisel A.
        • et al.
        RNA velocity of single cells.
        Nature. 2018; 560: 494-498
        • Aibar S.
        • Gonzalez-Blas C.B.
        • Moerman T.
        • et al.
        SCENIC: single-cell regulatory network inference and clustering.
        Nat Methods. 2017; 14: 1083-1086
        • Gregorieff A.
        • Stange D.E.
        • Kujala P.
        • et al.
        The ets-domain transcription factor Spdef promotes maturation of goblet and paneth cells in the intestinal epithelium.
        Gastroenterology. 2009; 137: 1333-1345.e1–3
        • Seymour P.A.
        Sox9: a master regulator of the pancreatic program.
        Rev Diabet Stud. 2014; 11: 51-83
        • Pin C.L.
        • Rukstalis J.M.
        • Johnson C.
        • et al.
        The bHLH transcription factor Mist1 is required to maintain exocrine pancreas cell organization and acinar cell identity.
        J Cell Biol. 2001; 155: 519-530
        • Arrojo E.D.R.
        • Jacob S.
        • Garcia-Prieto C.F.
        • et al.
        Structural basis for delta cell paracrine regulation in pancreatic islets.
        Nat Commun. 2019; 10: 3700
        • Wollny D.
        • Zhao S.
        • Everlien I.
        • et al.
        Single-cell analysis uncovers clonal acinar cell heterogeneity in the adult pancreas.
        Dev Cell. 2016; 39: 289-301
        • Schmidt P.H.
        • Lee J.R.
        • Joshi V.
        • et al.
        Identification of a metaplastic cell lineage associated with human gastric adenocarcinoma.
        Lab Invest. 1999; 79: 639-646
        • Nam K.T.
        • Lee H.J.
        • Sousa J.F.
        • et al.
        Mature chief cells are cryptic progenitors for metaplasia in the stomach.
        Gastroenterology. 2010; 139: 2028-2037.e9
        • Choi E.
        • Hendley A.M.
        • Bailey J.M.
        • et al.
        Expression of activated Ras in gastric chief cells of mice leads to the full spectrum of metaplastic lineage transitions.
        Gastroenterology. 2016; 150: 918-930.e13
        • Saenz J.B.
        • Mills J.C.
        Acid and the basis for cellular plasticity and reprogramming in gastric repair and cancer.
        Nat Rev Gastroenterol Hepatol. 2018; 15: 257-273
        • Goldenring J.R.
        Pyloric metaplasia, pseudopyloric metaplasia, ulcer-associated cell lineage and spasmolytic polypeptide-expressing metaplasia: reparative lineages in the gastrointestinal mucosa.
        J Pathol. 2018; 245: 132-137
        • Lee S.H.
        • Jang B.
        • Min J.
        • et al.
        Up-regulation of Aquaporin 5 defines spasmolytic polypeptide-expressing metaplasia and progression to incomplete intestinal metaplasia.
        Cell Mol Gastroenterol Hepatol. 2021; 13: 199-217
        • Hendley A.M.
        • Rao A.A.
        • Leonhardt L.
        • et al.
        Single-cell transcriptome analysis defines heterogeneity of the murine pancreatic ductal tree.
        Elife. 2021; 10e67776
        • Prasad N.B.
        • Biankin A.V.
        • Fukushima N.
        • et al.
        Gene expression profiles in pancreatic intraepithelial neoplasia reflect the effects of Hedgehog signaling on pancreatic ductal epithelial cells.
        Cancer Res. 2005; 65: 1619-1626
        • Bockerstett K.A.
        • Lewis S.A.
        • Wolf K.J.
        • et al.
        Single-cell transcriptional analyses of spasmolytic polypeptide-expressing metaplasia arising from acute drug injury and chronic inflammation in the stomach.
        Gut. 2020; 69: 1027-1038
        • Sakamoto N.
        • Fukuda K.
        • Watanuki K.
        • et al.
        Role for cGATA-5 in transcriptional regulation of the embryonic chicken pepsinogen gene by epithelial-mesenchymal interactions in the developing chicken stomach.
        Dev Biol. 2000; 223: 103-113
        • Asada R.
        • Saito A.
        • Kawasaki N.
        • et al.
        The endoplasmic reticulum stress transducer OASIS is involved in the terminal differentiation of goblet cells in the large intestine.
        J Biol Chem. 2012; 287: 8144-8153
        • Gracz A.D.
        • Samsa L.A.
        • Fordham M.J.
        • et al.
        Sox4 Promotes Atoh1-independent intestinal secretory differentiation toward tuft and enteroendocrine fates.
        Gastroenterology. 2018; 155: 1508-1523.e10
        • Gunawardene A.R.
        • Corfe B.M.
        • Staton C.A.
        Classification and functions of enteroendocrine cells of the lower gastrointestinal tract.
        Int J Exp Pathol. 2011; 92: 219-231
        • Bastidas-Ponce A.
        • Tritschler S.
        • Dony L.
        • et al.
        Comprehensive single cell mRNA profiling reveals a detailed roadmap for pancreatic endocrinogenesis.
        Development. 2019; 146: dev173849
        • Oster A.
        • Jensen J.
        • Edlund H.
        • et al.
        Homeobox gene product Nkx 6.1 immunoreactivity in nuclei of endocrine cells of rat and mouse stomach.
        J Histochem Cytochem. 1998; 46: 717-721
        • Baron M.
        • Veres A.
        • Wolock S.L.
        • et al.
        A single-cell transcriptomic map of the human and mouse pancreas reveals inter- and intra-cell population structure.
        Cell Syst. 2016; 3: 346-360.e4
        • Heller R.S.
        • Jenny M.
        • Collombat P.
        • et al.
        Genetic determinants of pancreatic epsilon-cell development.
        Dev Biol. 2005; 286: 217-224
        • Andralojc K.M.
        • Mercalli A.
        • Nowak K.W.
        • et al.
        Ghrelin-producing epsilon cells in the developing and adult human pancreas.
        Diabetologia. 2009; 52: 486-493
        • Haber A.L.
        • Biton M.
        • Rogel N.
        • et al.
        A single-cell survey of the small intestinal epithelium.
        Nature. 2017; 551: 333-339
        • Yamashita J.
        • Ohmoto M.
        • Yamaguchi T.
        • et al.
        Skn-1a/Pou2f3 functions as a master regulator to generate Trpm5-expressing chemosensory cells in mice.
        PLoS One. 2017; 12e0189340
        • Ramsay R.G.
        • Gonda T.J.
        MYB function in normal and cancer cells.
        Nat Rev Cancer. 2008; 8: 523-534
        • Kas K.
        • Finger E.
        • Grall F.
        • et al.
        ESE-3, a novel member of an epithelium-specific ets transcription factor subfamily, demonstrates different target gene specificity from ESE-1.
        J Biol Chem. 2000; 275: 2986-2998
        • Ghaleb A.M.
        • Yang V.W.
        Kruppel-like factor 4 (KLF4): What we currently know.
        Gene. 2017; 611: 27-37
        • Manco R.
        • Averbukh I.
        • Porat Z.
        • et al.
        Clump sequencing exposes the spatial expression programs of intestinal secretory cells.
        Nat Commun. 2021; 12: 3074
        • Deng T.
        • Karin M.
        c-Fos transcriptional activity stimulated by H-Ras-activated protein kinase distinct from JNK and ERK.
        Nature. 1994; 371: 171-175
        • Crawford H.C.
        • Scoggins C.R.
        • Washington M.K.
        • et al.
        Matrix metalloproteinase-7 is expressed by pancreatic cancer precursors and regulates acinar-to-ductal metaplasia in exocrine pancreas.
        J Clin Invest. 2002; 109: 1437-1444
        • Solar M.
        • Cardalda C.
        • Houbracken I.
        • et al.
        Pancreatic exocrine duct cells give rise to insulin-producing beta cells during embryogenesis but not after birth.
        Dev Cell. 2009; 17: 849-860
        • Bailey J.M.
        • Hendley A.M.
        • Lafaro K.J.
        • et al.
        p53 mutations cooperate with oncogenic Kras to promote adenocarcinoma from pancreatic ductal cells.
        Oncogene. 2016; 35: 4282-4288
        • McGinty J.W.
        • Ting H.A.
        • Billipp T.E.
        • et al.
        Tuft-cell-derived leukotrienes drive rapid anti-helminth immunity in the small intestine but are dispensable for anti-protist immunity.
        Immunity. 2020; 52 (e7): 528-541
        • Kuga D.
        • Ushida K.
        • Mii S.
        • et al.
        Tyrosine phosphorylation of an actin-binding protein girdin specifically marks tuft cells in human and mouse gut.
        J Histochem Cytochem. 2017; 65: 347-366
        • Farrell A.S.
        • Joly M.M.
        • Allen-Petersen B.L.
        • et al.
        MYC regulates ductal-neuroendocrine lineage plasticity in pancreatic ductal adenocarcinoma associated with poor outcome and chemoresistance.
        Nat Commun. 2017; 8: 1728
        • Gribble F.M.
        • Reimann F.
        Function and mechanisms of enteroendocrine cells and gut hormones in metabolism.
        Nat Rev Endocrinol. 2019; 15: 226-237
        • Molina-Montes E.
        • Coscia C.
        • Gomez-Rubio P.
        • et al.
        Deciphering the complex interplay between pancreatic cancer, diabetes mellitus subtypes and obesity/BMI through causal inference and mediation analyses.
        Gut. 2021; 70: 319-329
        • Miao Z.F.
        • Lewis M.A.
        • Cho C.J.
        • et al.
        A dedicated evolutionarily conserved molecular network licenses differentiated cells to return to the cell cycle.
        Dev Cell. 2020; 55: 178-194.e7