ROR2 Regulates Cellular Plasticity in Pancreatic Neoplasia and Adenocarcinoma

Pancreatic ductal adenocarcinoma (PDAC) is a deadly disease with few curative options for the vast majority of patients, in part because early detection methods are lacking. PDAC is marked by extensive tumor cell heterogeneity comprised of “classical” and “basal-like” molecular subtypes that, instead of acting as stable cellular communities, exist on a spectrum, with plasticity between subtypes likely contributing to therapeutic resistance (1–3). Cellular plasticity is a hallmark of pancreatic cancer throughout progression, including transdifferentiation of normal cells into precancerous metaplastic and neoplastic lesions accompanied by the acquisition of distinct cellular identities. Whether cellular plasticity in precancerous disease is reflected in late-stage cancer progression has not been thoroughly considered, though some studies suggest a connection (4, 5).

Both pancreatic duct and acinar cells can give rise to PDAC in mouse models (6–8). Acinar cells are more amenable to transformation in mice (8) and epigenetic analysis in patient samples suggests an acinar cell of origin (9). Inflammation or expression of oncogenic Kras^G12D induces transdifferentiation of acinar cells to duct-like progenitor cells, a process called acinar-to-ductal metaplasia (ADM; refs. 6, 7, 10). Metaplastic duct cells are characterized by an upregulation of genes that define normal duct cells, such as Krt19 and Sox9, and by loss of acinar differentiation markers, such as those encoding digestive enzymes (8, 11, 12).

The great majority of adult acinar cells are resistant to KRAS^G12D-induced transformation, but in combination with the ablation of transcription factors that otherwise maintain acinar cell identity, such as Hnf1a, Nr5a2, Mist1, Ptf1a, and Pdx1, transformation is significantly accelerated (13–17). Of these, genetic variants in HNF1A, NR5A2, and PDX1 are associated with susceptibility to PDAC development in patients (14, 18). These findings support the hypothesis that maintenance of acinar cell identity is a powerful suppressor of pancreatic cancer initiation.

In both injury and neoplasia models, metaplastic ducts show a high degree of cellular heterogeneity with subpopulations marked by senescent and proliferative characteristics and a variety of gastric cell identities including gastric chief-, pit-, tuft-, and enteroendocrine-like cells (19). In the stomach, chief cells secrete digestive enzymes and pit cells produce mucins that protect the epithelium from gastric juice (20); the functions of the chief and pit cell-like populations in pancreatic neoplasia are unclear. Notably, in experimental pancreatitis, a subset of metaplastic duct cells also assume a gastric SPEM (Spasmolytic Polypeptide-Expressing Metaplasia) identity (21). In the stomach, SPEM arises after glandular injury from corpus chief cells, and SPEM cells characteristically co-express markers of gastric chief and mucous neck cells, such as Pepsinogen C (Pgc) and Mucin6 (Muc6), respectively (22).

Although gastric cell identity is evident in pancreatic precancerous lesions, it has not been noted in carcinomas, which are instead marked by the classical and basal-like subtypes. The classical subtype exhibits higher expression of pancreatic lineage- and epithelial-defining transcription factors, whereas the basal-like subtype is characterized by an erosion of epithelial cell identity and acquisition of mesenchymal and basal cell features (23–27). Importantly, subtype identity of pancreatic cancer cells is associated with prognosis, with the basal-like subtype correlating to poorer overall survival (1).

The complex heterogeneity in both precancerous lesions and carcinoma has been of intense interest. Unfortunately, the study of the earliest molecular changes associated with the reprogramming of normal pancreatic acinar cells to neoplasia has been challenging using standard single-cell RNA-sequencing methods, largely due to the fragility of normal acinar cells during tissue disruption, together with high levels of digestive enzymes, especially nucleases. To overcome this, we performed single-nucleus ATAC-sequencing (snATAC-seq) of undigested, rapidly pulverized, frozen bulk mouse pancreas tissue to identify possible mechanisms of cellular reprogramming initiated by the inducible ablation of the acinar cell identity factor Pdx1, the expression of oncogenic KRAS^G12D, or both together. Strikingly, chromatin accessibility of marker genes was sufficient to determine cell identity of multiple cell types within the tissue, including epithelial, fibroblast, and immune cell populations. In transformation-sensitive Pdx1-null acinar cells, we identified increased accessibility and expression of several genes associated with a gastric SPEM signature, which was intensified in acinar cells also expressing Kras^G12D. Among these genes, we focused on Ror2, which encodes a receptor tyrosine kinase in the noncanonical WNT signal transduction pathway. Genetic ablation of Ror2 in a mouse model of pancreatic neoplasia shifted cell identity, giving rise to precancerous lesions dominated by mucinous gastric pit-like cells. In PDAC, ROR2 expression strongly correlated with a basal-like/mesenchymal phenotype in which it controls cell proliferation. Forcing ROR2 expression in classical subtype PDAC cell lines was sufficient to induce a complete conversion to rapidly proliferating, poorly differentiated mesenchymal cells. The concomitant shift in signal transduction induced by ROR2 expression conferred resistance to a KRAS^G12D inhibitor while introducing a profound sensitivity to AKT inhibition.

To initiate pancreatic tumorigenesis, oncogenic KRAS must overcome resistance to transformation posed by acinar cell identity (13–17, 28). In order to study the earliest steps in this process, we compared pancreata with conditional adult acinar cell specific expression of Kras^G12D alone, to those also lacking Pdx1, a developmental transcription factor that we showed previously to be important for maintaining acinar cell differentiation (17). We used Ptf1aCre^ERT-driven recombination to activate Kras^G12D expression and ablate Pdx1 by treating 8- to 9-week-old mice with Tamoxifen, followed by tissue collection 2 and 10 weeks later to capture early progressive cell identity changes. Sole expression of Cre^ERT (Ptf1a^ERT) and loss of Pdx1 (Ptf1a^ERT;Pdx1^f/f) did not cause any obvious morphologic changes in the epithelial compartment, with the majority of cells being Amylase-positive acinar cells (Fig. 1A and B). Oncogenic Kras expression in Ptf1aCre^ERT;Kras^LSL-G12D (Ptf1a^ERT;K*) mice led to sporadic ADM, characterized by loss of Amylase, expression of the ductal marker KRT19, and a surrounding stromal reaction indicated by Vimentin staining (Fig. 1A and B). In contrast to the low metaplastic transformation efficiency of KRAS^G12D alone, the concomitant loss of Pdx1 significantly accelerated cellular reprogramming in Ptf1aCre^ERT;Kras^LSL-G12D;Pdx1^flox/flox (Ptf1a^ERT;K*;Pdx1^f/f) animals, as previously reported (17). Only 10 weeks after Tamoxifen, acinar cells had completely transdifferentiated into metaplastic duct-like cells, accompanied by a fibroinflammatory stromal reaction (Fig. 1A and B).

To better define how Kras^G12D expression alone and in combination with Pdx1 ablation impacts acinar cell identity, we established a novel protocol allowing high-quality snATAC-seq of snap-frozen, rapidly pulverized bulk pancreatic tissue (Supplementary Fig. S1A). Overall, this protocol allows for very rapid processing of samples, strongly preserving cellular identity. After extensive quality control and removal of low-quality nuclei (Supplementary Fig. S1B and S1C), we performed downstream analysis with a total of 33,505 nuclei and compared chromatin accessibility states of Ptf1a^ERT and Ptf1a^ERT;Pdx1^f/f control animals to mice with Kras^G12D expression, sacrificed 2 and 10 weeks post Tamoxifen.

Based on chromatin accessibility at the gene body and promoter, we calculated predicted gene expression activity scores by using the Signac pipeline (29) and used top predicted differentially accessible genes as well as previously established marker genes to assign cell cluster identity (Supplementary Fig. S1D and S1E). We identified acinar, ductal, and meta/neoplastic cells as well as fibroblasts, mesothelial, endothelial, immune, and islet cells (Fig. 1C; Supplementary Fig. S1E). Surprisingly, although seemingly morphologically normal, acinar cells of different mouse genotypes grouped into very distinct cell clusters, whereas other cell types remained unperturbed across genotypes at this resolution. Notably, acinar cells of Ptf1a^ERT control (Acinar, yellow cluster) and Ptf1a^ERT;K* (K* Acinar, brown cluster) mice showed slight differences, whereas acinar cells of Pdx1 knockout mice (Pdx1^f/f Acinar, orange cluster) exhibited an obvious switch of cell identity, independent of KRAS status (Fig. 1C; Supplementary Fig. S1E). Acinar cell clusters were defined by predicted expression of digestive enzymes, such as Cela1 and Cpa1. Cela1 showed lower expression levels in acinar cells of Ptf1a^ERT;K* and Ptf1a^ERT;K*;Pdx1^f/f animals, an indicator of distinct differentiation states (Fig. 1D). Consistent with the histology, expression of Kras^G12D led to time- and genotype-dependent meta/neoplastic duct cell formation (teal cluster), characterized by expression of metaplastic and ductal marker genes, such as Onecut2 and Hnf1b, respectively (Fig. 1C–E). These data were consistent with the rampant metaplastic and neoplastic transformation in Ptf1a^ERT;K*;Pdx1^f/f animals 10 weeks post Tamoxifen, with acinar cells being almost completely absent, replaced by an immense increase in meta/neoplastic duct cells (Fig. 1C and E). These massive changes in the epithelial compartment were accompanied by a substantial increase of fibroblasts (light-green cluster), with little change in immune cell numbers (violet cluster; Fig. 1E). Of note, the desmoplastic reaction in these mice made tissue digestion and isolation of high-quality nuclei more challenging, resulting in a lower overall number (Fig. 1C).

Our snATAC-seq approach corroborated histologic observations and changes in cell populations that occur during the early stages of tumorigenesis of the pancreas. Importantly, we uncovered that acinar cells from different mouse genotypes segregate into distinct clusters with an obvious shift in identity in metaplastic transformation-sensitive Pdx1-null acinar cells, revealing previously elusive details of early cellular reprogramming. Next, we sought to investigate gene accessibility changes among our identified acinar cell clusters in more detail.

To unmask cell identity differences across our acinar cell clusters, we subsetted acinar cells from the 2-week time point based on their genotype (Fig. 2A). To validate the identification of differentially regulated genes, we used an integrative approach combining our snATAC-seq with bulk RNA-seq data of isolated acinar cells from matching mouse genotypes and sex. Genes that showed the highest correlation across the two techniques are depicted in Fig. 2B (Supplementary Fig. S2A). One gene confirming the validity of this approach was Kdm6a, which is located on the X chromosome and was both more accessible and more highly expressed in females (Fig. 2B). Acinar cells with Pdx1 loss exhibited an obvious decrease of accessibility and expression of acinar differentiation genes, such as Cela1 or Ptf1a (17), but elevated levels of the AP-1 transcription factor Junb, which is associated with transcriptional reprogramming in ADM (Fig. 2B; Supplementary Fig. S2A; ref. 30).

Among the top genes that showed the most significant enrichment in accessibility and expression in metaplastic transformation-sensitive Pdx1-null acinar cells was the receptor tyrosine kinase gene Ror2 (Receptor Tyrosine Kinase Like Orphan Receptor 2; Fig. 2B; Supplementary Fig. S2A). ROR2 regulates noncanonical WNT signaling and defines cancer cell identity in a variety of tumor types (31–33). Besides Ror2, unsupervised clustering analysis of the RNA-seq data revealed that gastric lineage-specific genes, including Cobalamin binding intrinsic factor (Cblif), Pepsinogen5 (Pga5), Gastrokine3 (Gkn3), and Aquaporin5 (Aqp5; refs. 34, 35), are upregulated in Pdx1-knockout acinar cells, whereas acinar differentiation genes are down-regulated (Fig. 2C; Supplementary Fig. S2B). Considering that the acquisition of gastric signatures has recently been identified as a hallmark of pancreatic tumorigenesis (19, 36), we asked if Ror2 itself is associated with gastric cell identity.

Analysis of published mouse scRNA-seq data derived from whole stomach tissue (19) indeed revealed Ror2 expression in three distinct gastric cell populations (Fig. 2D). For gastric cell type annotation, we used a previously annotated mouse scRNA-seq dataset generated from the combined gastric corpus and antral regions (34) and found that the Ror2-expressing population 1 is associated with an endocrine, enterochomaffin-like cell identity, whereas populations 2 and 3 show an enrichment of gastric neck cell markers (Fig. 2D; Supplementary Fig. S2C). In accordance with our RNA-seq data from acinar cells, we detected some cells that co-express Ror2 with Cblif and Aqp5 in the stomach (Fig. 2C and D). IHC analysis of mouse stomach demonstrated ROR2 protein expression at the base, including AQP5-positive deep antral neck cells in the antrum (Supplementary Fig. S2D). Of note, gene expression patterns of deep antral cells share characteristics of gastric metaplasia (SPEM), which forms in the corpus region upon gastric injury as part of a wound-healing program. SPEM cells arise from corpus chief cells, co-express chief, and neck cell markers and resemble normal deep antral cells (22). Researchers recently identified the SPEM cell identity in pancreatitis-induced metaplasia (21), a putative precursor to neoplasia.

In addition to Ror2, we found that SPEM-associated genes, such as Gkn3, Aqp5 (gastric neck cell markers) and Pgc (chief cell marker), were upregulated in acinar cells with Pdx1 ablation, especially in combination with oncogenic Kras expression (Fig. 2C). Gene set enrichment analysis (GSEA) for previously published SPEM markers (21, 37–39) confirm that the SPEM signature is significantly upregulated in Ptf1a^ERT;K*;Pdx1^f/f acinar cells (Fig. 2E). Thus, PDX1 actively maintains acinar cell identity, while it suppresses a gastric/SPEM identity and susceptibility to neoplastic transformation by KRAS^G12D. Notably, we identified Ror2 as a previously-unrecognized marker of the gastric cell identity-associated gene signature.

Given that oncogenic KRAS-induced meta/neoplastic cells in the pancreas acquire features reminiscent of gastric cell lineages (19, 36), we examined if ROR2 expression was associated with this identity shift. By performing RNA-ISH and IHC, we found that Ror2 is indeed expressed in a portion of meta/neoplastic cells in Ptf1a^ERT;K* animals, comprising ∼25% of KRT19-positive cells, especially in later time points (Fig. 3A). To explore whether Ror2 expression levels define meta/neoplastic duct cell identity, we used pancreas scRNA-seq data from Schlesinger and colleagues (19), which contained 3,309 lineage-traced meta/neoplastic cells isolated from Ptf1a^ERT;K*;ROSA^LSL-TdTomato mice (Fig. 3B). In concordance with our IHC data, Ror2 is expressed in subpopulations of tdTomato-positive meta/neoplastic cells (Fig. 3B). To determine the identity of ROR2-positive cells, we examined the expression of gastric cell type signatures from normal mouse gastric epithelia (34) as well as typical SPEM genes. We found that the ROR2^High cell population (Fig. 3B, circled) is indeed associated with higher expression of gastric neck cell and SPEM signatures. Consequently, we named this cluster neck/SPEM cell-like (Fig. 3C).

In accordance with Schlesinger and colleagues (19), we identified a pit cell-like population, confirmed by Tff1 expression, a dividing population marked by expression of Mki67, and a senescent population with high expression of tumor suppressor genes, such as Cdkn2a (Fig. 3C; Supplementary Fig. S3A–S3C). Three clusters did not exhibit clear segregation of marker gene expression, which we designated as “hybrid.” Although the “acinar-like” cluster retained expression of acinar genes, expression of ductal genes, such as Dcdc2a, was enriched in the “duct-like” cell cluster (Fig. 3C; Supplementary Fig. S3D and S3E). Of note, we found the highest expression of Ror2 in the neck/SPEM cell-like cluster, followed by dividing and acinar-like meta/neoplastic cell populations (Fig. 3D). Among the top 10 genes that showed the highest correlation to Ror2 expression, we found the gastric neck and SPEM marker Muc6, as well as the SPEM genes Aqp5, Gkn3, or Pgc (Fig. 3D and E), reminiscent of the gene signature found in Pdx1-null acinar cells (Fig. 2C). Validation studies using multiplex immunofluorescence (IF)/ISH staining confirmed co-localization of Ror2 with Muc6 expression in ductal lesions of Ptf1a^ERT;K* mice, 20 weeks post Tamoxifen (Fig. 3F). Moreover, we detected that a proportion of ROR2-positive meta/neoplastic cells indeed co-expressed Ki67, whereas co-localization with CDKN2A or TFF1 was rare (Fig. 3G–I).

Interestingly, ROR2-positive lesions showed lower abundance of Krt19 expression, consistent with the SPEM population found in injury-induced ADM (Fig. 3F; Supplementary Fig. S3F; ref. 21). In a model of chronic pancreatitis (21), a substantial proportion of “Mucin/Ductal” cells also expressed Ror2, again showing correlation with neck cell and SPEM marker gene expression, such as Muc6 or Aqp5 (Supplementary Fig. S3G–S3J). In summary, our data suggested that ROR2-positive meta/neoplastic cells are of a proliferative neck/SPEM cell identity, distinct from cells with senescent and pit cell like features. Given that the relevance of these cellular identities on pancreatic cancer development is unknown, we next sought to determine if ablation of the neck/SPEM marker ROR2 impacts carcinogenesis in mice.

To assess the functional relevance of Ror2 in early metaplasia and neoplasia, we generated Ptf1a^ERT;K* mice with a conditional knockout allele of Ror2 and a ROSA^LSL-YFP reporter. Ptf1aCre^ERT-mediated recombination led to the deletion of Exon3/4 of the Ror2 gene in the majority of acinar cells (Supplementary Fig. S4A). IHC for the YFP reporter expression suggested a recombination efficiency of ∼90%. We noted that some lesions retained ROR2 expression, suggesting either less efficient recombination of this locus or possible selection for cells that escaped recombination (Supplementary Fig. S4B). Although we observed slight differences in tissue remodeling, quantified as loss of acinar cells and formation of meta/neoplastic duct-like structures between Ptf1a^ERT;K* and Ptf1a^ERT;K*;Ror2^f/f animals at 10 weeks post Tamoxifen (Supplementary Fig. S4C), Ptf1a^ERT;K*;Ror2^f/f animals presented with significantly more transformation at 20 and 40 weeks post Tamoxifen (Fig. 4A). Notably, 40 weeks post Tamoxifen, 3 Ptf1a^ERT;K*;Ror2^f/f mice developed regions of carcinoma, a rapidity of progression not observed in our control cohort. Also, they showed a more pronounced stromal reaction with significantly elevated collagen deposition as represented by Picrosirius Red stain, which was also observed when Ror2 was ablated in a mammary cancer model (Supplementary Fig. S4D; ref. 40). When comparing lesion identity in mice that presented with neoplasia, we found that Ptf1a^ERT;K*;Ror2^f/f animals exhibit an increase in the number of mucinous meta/neoplastic lesions, the first indication of a possible shift in epithelial cell identity (Fig. 4A).

Based on our previous observation that ROR2 was mainly expressed in proliferative neck/SPEM subpopulations of meta/neoplastic cells, but negatively associated with senescent (CDKN2A-positive) and pit-like (TFF1-positive) cells, we analyzed Ror2 knockout mice for a possible shift in cellular identity. Characterization of KRT19-positive ductal structures at 20 weeks post Tamoxifen revealed that lesions in Ror2 knockout mice were less proliferative and had a significantly higher number of CDKN2A-positive senescent cells, with the senescent phenotype difference having abated by 40 weeks (Fig. 4B and C). The high prevalence of mucinous lesions at 40 weeks led us to stain for MUC6, a SPEM marker, as well as TFF1 and MUC5AC, both markers for gastric pit cell identity. We observed a slight reduction of MUC6-expressing KRT19-positive cells (Fig. 4D), accompanied by a higher proportion of ductal cells expressing MUC5AC and TFF1 when Ror2 was ablated (Fig. 4E and F). Also, this shift of cell identity was evident in immunoblot analyses from bulk pancreatic tissue lysates confirming increased TFF1 levels in Ror2-null, Kras^G12D-expressing mice (Fig. 4G).

Altogether, these data show that ablation of Ror2 promotes tissue transformation and fibrosis. Also, Ror2 ablation shifted cellular identity of precancerous lesions, with an increased number exhibiting a mucinous TFF1-positive gastric pit cell-like identity and an overall higher likelihood of progression to PDAC. These two concurrent phenotypes suggest a possible connection between the pit-like phenotype and progression to carcinoma, a possibility bolstered by TFF1 also being a well-established marker for the classical PDAC molecular subtype (24, 26).

To assess if Ror2 and gastric gene signatures define subtype identity in PDAC, we turned to the well-established KPC (Ptf1aCre;Kras^LSL-G12D;Trp53^R172H) mouse model, which recapitulates human PDAC subtype heterogeneity, including the epithelial, classical, and the more mesenchymal, basal-like subtypes (41). Dual IHC for ROR2 (brown) and TFF1 (teal) revealed that expression is very heterogeneous in neoplastic and carcinoma areas in KPC mice and confirmed that ROR2 and TFF1 define distinct, nonoverlapping subpopulations (Fig. 5A). Accordingly, analysis of previously published RNA-seq data generated from mouse pancreatic cancer cell lines (42) also showed a significant negative correlation between the two genes. Moreover, although Ror2^Low/Tff1^High cells are characterized by a differentiated epithelial cell identity (clusters C2a-C2c), Ror2^High/Tff1^Low expressors are associated with a more mesenchymal cell differentiation state (C1; Fig. 5B).

Using Pearson correlation analysis, we further uncovered that Ror2 does not only negatively correlate to Tff1, but to myriad genes that define meta/neoplastic pit-like cells. Importantly, the great majority of most enriched pit cell-like marker genes (n = 50) show a strong association with the classical subtype, including nine genes (highlighted in red) that are part of a recently defined 30 gene human classical PDAC signature (Fig. 5C; ref. 27). To validate these findings in human PDAC tissue, we surveyed publicly available scRNA-seq data from 174 human PDAC samples, to our knowledge, the largest aggregated collection of tissues from patients with PDAC (2, 43–51). Ductal cells (KRT18 and KRT19 enriched, negative for acinar and other cell lineage markers) were subsetted from total cells and re-clustered in Seurat. GSEA within the ductal population revealed a tight overlap between the pit cell-like and the classical PDAC gene signature (Fig. 5D; Supplementary Fig. S5A). These data suggest that cell fate traits of meta/neoplastic pit-like cells are maintained in classical PDAC in mouse models and importantly, in patient tumors.

To identify cellular identity features that show high correlation to Ror2, we re-visited the RNA-seq data of mouse pancreatic cancer cell lines (42) given its large number of samples that represent classical or basal-like PDAC subtype identity. Strikingly, genes that show highest correlation to Ror2 expression exhibit markedly elevated expression in the mesenchymal cell lines (C1) associated with the basal-like subtype, as defined by gene signatures from Moffitt and colleagues (Fig. 5E; ref. 23). To elucidate the biological function of genes that positively and negatively correlated with Ror2, we conducted GO pathway analysis. Ror2 positively-correlated genes were associated with cancer cell aggressiveness and epithelial-mesenchymal transition (EMT), whereas negatively correlated genes were regulators of apoptosis or epithelial differentiation (Supplementary Fig. S5B). Interestingly, the most significantly enriched GO term for positively-correlated genes “multicellular organism development” contained a myriad of well-known EMT-inducing transcription factors, such as PRRX1, SNAIL1, TWIST, or ZEB family members. Correlation plots present significant positive correlation of Ror2 with Snai1, Twist1, Zeb1, or Vimentin (Vim) expression, whereas expression of epithelial genes E-Cadherin (Cdh1) and Krt19 is negatively correlated with Ror2 (Fig. 5F). Interestingly, Ror2-positive pancreatic cancer cells were associated with lower Krt19 expression, a finding that we observed previously in ROR2-positive precancerous lesions.

Based on our finding that Ror2 correlated to the basal-like/mesenchymal PDAC subtype in mice, we asked if ROR2 showed a similar association in human PDAC. For this, we examined ROR2 expression in resected tumors, patient-derived organoids and cancer cell lines. Within ROR2-expressing tumor samples, we observed significant heterogeneity, with highest expression correlating to moderate to poorly differentiated carcinoma (Fig. 6A). Given that ROR2 is highly expressed in fibroblasts, which can confound analysis of bulk RNA-seq data from whole tumor tissue, we analyzed data from human micro-dissected PDAC epithelium (52) and found that ROR2 expression correlates with the expression of EMT genes, such as ZEB1 and VIM, whereas CDH1 and EPCAM inversely correlate (Fig. 6B). Of note, basal-like tumor samples are highly underrepresented in human PDAC, likely due to these aggressive tumors not usually qualifying for resection and downstream analysis. Accordingly, when we analyzed human organoid (Supplementary Fig. S6A) and pancreatic cancer cell lines, including primary patient lines, we found that the majority exhibited more epithelial characteristics and were negative for ROR2 expression. However, all lines that did express ROR2 also showed mesenchymal features, such as Vimentin and ZEB1 expression, and extremely low levels of the E-Cadherin (Fig. 6C and D). We conclude that ROR2 is strongly associated with a mesenchymal differentiation in human PDAC.

Next, we assessed if expression of ROR2 in the epithelial cell lines HPAF-II and the primary patient cell line UM32 is sufficient to cause cellular reprogramming and a switch of cell identity. Notably, overexpression of ROR2 clearly induced EMT, characterized by loss of the epithelial and gain of the spindle-shaped, mesenchymal morphology (Fig. 6E). IF and immunoblot analyses demonstrated concomitant expression of the EMT markers ZEB1 and Vimentin, whereas E-Cadherin was completely lost (Fig. 6E and F). Consequently, we compared RNA-seq data of HPAF-II and UM32 control and ROR2-overexpressing cells to the EMT signature established by Groger and colleagues (53). ROR2-overexpressing HPAF-II and UM32 cells displayed a significantly enriched EMT score compared with control cells (Fig. 6G). In addition, expression of previously identified classical PDAC subtype and pit cell marker genes were completely abrogated upon ROR2 expression, overall indicating that ROR2 erodes pancreatic cancer cell identity and promotes acquisition of a mesenchymal phenotype (Supplementary Fig. S6B).

Previous studies revealed that activation of EMT programs induce stem cell features, with ZEB1 being a key factor for stemness in various cancer entities (4, 54, 55). Our RNA-seq data unveiled that the gastric stem cell genes LGR5 or CXCR4 are greatly upregulated upon expression of ROR2 in HPAF-II and UM32 cells (Supplementary Fig. S6C). To test if this stem cell-like signature predicts function, we examined the ability of ROR2-expressing cells to form tumor spheres in low adhesion plates. Indeed, ROR2-expressing cells formed a significantly greater number and larger tumor spheres in vitro compared to control cells (Fig. 6H). Overall, our data identify ROR2 as a driver of human PDAC progression, inducing EMT and stemness.

Next, we studied the effects of ROR2 loss in the ROR2-expressing, mesenchymal PANC1 and UM5 cell lines. CRISPR/Cas9-mediated knockout using two different guide RNAs (ROR2 KO #1, ROR2 KO #2) had no consistent effect on EMT marker gene expression (Supplementary Fig. S7A), demonstrating that ROR2 signaling is not necessary for maintaining EMT once initiated. However, ablation of ROR2 in PANC1 and UM5 led to a significant decrease in cell proliferation (Fig. 7A), an effect that was particularly dramatic in the patient-derived primary cell line UM5, so much so that recovery of RNA/protein or execution of functional assays were extremely challenging (Supplementary Fig. S7B). To study immediate ROR2 target genes, we performed siRNA-mediated acute knockdown of ROR2 and performed RNA-seq. Top enriched GSEA gene sets in control versus ROR2 knockdown cells were associated with fundamental regulators of proliferation, such as cell cycle-regulating E2F targets and G2/M transition genes (Supplementary Fig. S7C), an effect that has been observed in other tumor types (56, 57). In contrast, ectopic overexpression of ROR2 in the ROR2-negative cell lines HPAF-II and UM32 induced an extensive increase of cell proliferation (Fig. 7B), correlating with an upregulation of E2F and G2/M target genes (Supplementary Fig. S7D).

Next, we used exploratory phospho-Arrays to identify downstream signaling of ROR2 and confirmed these findings by immunoblot analysis. Although ROR2 knockout cells, especially UM5, showed a reduction in phospho-AKT levels, overexpression of ROR2 led to a marked increase, accompanied by a virtual ablation of activated phospho-ERK (Fig. 7C and D), a surprising effect considering both cell lines harbor mutant KRAS^G12D. To test if this profound shift away from ERK and toward AKT signaling would alter response to pharmacologic inhibitors, we tested HPAF-II and UM32 ROR2-overexpressors’ response to the KRAS^G12D inhibitor MRTX1133 (Fig. 7E) and the AKT inhibitor MK-2206 (Fig. 7F). Indeed, we found that ROR2 expressing cells were more resistant to KRAS^G12D inhibition compared with controls, but newly vulnerable to AKT inhibition.

Pancreatic cancer is an almost universally fatal disease marked by immense tumor cell heterogeneity that contributes to resistance to therapy (1). This tumor cell heterogeneity, in turn, is driven by the profound plasticity of the exocrine pancreas (58), a quality evident in acinar cell transdifferentiation to metaplastic duct-like cells in chronic pancreatitis, and in their susceptibility to form precancerous lesions induced by the expression of oncogenic KRAS (19, 21). In this study, we applied single-nucleus ATAC-seq to capture early reprogramming of acinar cells induced by the expression of Kras^G12D, as well as in acinar cells with a compromised identity due to Pdx1 ablation. We discovered that ablation of Pdx1 alone initiates reprogramming toward a gastric cell identity, an effect also seen upon ablation of Ptf1a, another transcription factor critical for the maintenance of acinar cell identity (59). Gastric cell gene expression was intensified by the expression of Kras^G12D. Among the upregulated genes we identified was the noncanonical WNT receptor ROR2, which we found is expressed in the deep antral cells of the stomach, together with many markers of gastric metaplasia (SPEM) (​22).

In progressive neoplasia in the mouse, we found that ROR2 remains associated with and promotes a proliferative neck cell/SPEM-like identity, whereas its ablation resulted in more extensive neoplastic transformation. This is consistent with SPEM playing a wound-healing function in neoplasia as it does in normal stomach and pancreas, possibly suppressing the progressive loss of normal tissue (60). Early neoplastic lesions lacking Ror2 were highly senescent, whereas in late neoplasia, the tissue became dominated by a gastric pit cell-like cell identity, suggesting that Ror2 is actively involved in maintaining cell identity. Interestingly, a subset (3/11) of Ptf1a^ERT;K*;Ror2^f/f mice developed carcinoma at a very early 40-week time point, a progression speed not observed in control Ptf1a^ERT;K* animals, either in this study or historically. The observation that the classical subtype of PDAC maintains a mostly pit cell-like identity signature leads us to hypothesize that the pit cell-like neoplastic cells are the immediate precursors to carcinoma. If this is the case, either their sheer numbers or their selective ability to overcome senescence and low proliferation earlier in progression would make PDAC more likely to form.

In PDAC, ROR2 expression is associated with an aggressive basal-like molecular subtype and a mesenchymal phenotype. Its ectopic expression in human PDAC cell lines with a classical subtype identity was sufficient to promote extensive reprogramming toward a mesenchymal and highly proliferative phenotype, a function observed in other tumor types (31, 32). Though the possibility exists that the ROR2-expressing neck/SPEM neoplasia is a precursor to basal PDAC, other than ROR2, we did not observe substantial expression of other SPEM markers in this subtype. A higher likelihood exists that ROR2 actively drives mesenchymal differentiation in a manner that reflects its role in fibroblasts, in which it is consistently expressed in normal, neoplastic, and cancerous pancreas.

Overall, our data point to a complex function of ROR2, resisting progression in precancerous stages and promoting progression in PDAC. Linking these two functions is ROR2’s regulation of cellular plasticity. Cellular plasticity is increasingly recognized as a mechanism of therapeutic resistance in cancer (27, 61). ROR2 has been shown to induce EMT in several tumor types (31–33, 62), though its relationship to ERK (63–65) and AKT signaling pathways (65–67) seems to be highly context-dependent. We found that ROR2 expression in two KRAS^G12D-mutated PDAC cell lines exhibits a dramatic effect on the classic effector pathways downstream of KRAS, virtually eliminating ERK activity and drastically upregulating AKT activity. ROR2’s hyperactivation of AKT was consistent with observations in several tumor types (66, 67). This remarkable shift away from a dependency on ERK signaling presaged resistance to KRAS inhibition, whereas the hyperactivation of AKT created an increased vulnerability to AKT inhibition by MK-2206. To further explore its role in resistance to KRAS inhibition, we examined ROR2 expression in cell lines that were selected for resistance to MRTX1133 in vitro and in vivo and found its expression inconsistent, suggesting that ROR2 upregulation is not a common mechanism of acquired resistance. However, our data suggest that its expression within a subset of aggressive, basal-like PDACs is likely to be predictive of therapeutic efficacy of KRAS and AKT inhibitors.

In addition to targeting its downstream effectors, the recent advent of inhibitory small molecules and monoclonal antibodies (68, 69) will also allow for direct targeting of ROR2 itself in the future. Another approach to targeting ROR2 activity is indirectly, by targeting its ligands. We consistently observed expression of the ROR2 ligand WNT5A in PDAC epithelia and stroma (data now shown), in accordance with other studies demonstrating its prevalent expression in primary tumors where it promotes EMT (70), as well as in subsets of circulating tumor cells from patients with PDAC (71). However, treatment of our ROR2-expressing human PDAC cell lines with the Porcupine O-acyltransferase (PORCN) inhibitor LGK-974, an agent that inhibits posttranslational modification and secretion of WNT ligands, did not reduce cell proliferation as we observed in our ROR2 knockout cells. This suggests that ROR2 may demonstrate WNT ligand-independent activities or that inhibition of both canonical and noncanonical WNTs together demonstrates more complicated effects than targeting the noncanonical arm alone.

In summary, we found that erosion of pancreatic acinar cell identity induced by Pdx1 ablation was sufficient to confer a gastric identity gene signature that was maintained up to and including progression to the classical subtype of PDAC. Among the neck cell/SPEM subpopulation of metaplastic epithelium, we identified ROR2 to be a powerful regulator of cellular identity throughout progression. In early-stage disease, ROR2 suppresses neoplastic transformation and progression to PDAC. In carcinoma, ROR2 promotes progression to an aggressive basal-like subtype in which it modulates therapeutic responses to KRAS and AKT inhibition.

The human pancreatic cancer cell lines PANC1 (RRID:CVCL_0480) and HPAF-II (RRID:CVCL_0313) were purchased from ATCC and UM5 and UM32 were kindly provided by Dr. Costas Lyssiotis (University of Michigan). Identity of cell lines was confirmed by STR analysis (University of Arizona). PANC1 were cultivated in DMEM (10-013-CV, Corning), supplemented with 10% FBS (F2242, Thermo Fisher Scientific) and 0.5% Gentamicin (15-710-072, Fisher Scientific). HPAF-II cells were cultured in Minimum Essential Medium (MEM, 11-095-080, Fisher Scientific) containing 1% sodium pyruvate (11360070, Thermo Fisher Scientific), 1% nonessential amino acids (11140050, Thermo Fisher Scientific), 10% FBS, and 0.5% Gentamicin. UM5 and UM32 were cultivated in RPMI 1640 (10-040-CV, Corning) with 10% FBS and 0.5% Gentamicin. 293FT cells were grown in DMEM supplemented with 1% sodium pyruvate, 1% nonessential amino acids, 10% FBS, and 1% Penicillin/Streptomycin. Generally, cells were passaged at a confluency of ∼80% by using Trypsin/EDTA (25200114, Thermo Fisher Scientific). All cell lines were regularly tested for mycoplasma using the MycoStrip-Mycoplasma Detection Kit (Rep-mys-50, Invivogen) and maintained at 37°C in a humidified chamber with saturated atmosphere containing 5% CO₂.

Organoid cultures 163, 265, 281, and 286 were previously established from PDX tumors (72). Organoids were maintained in Cultrex UltiMatrix Reduced Growth Factor Basement Membrane Extract (BME001-10, R&D Systems) and incubated in Advanced DMEM/F12 (12-634-010, Thermo Fisher Scientific) culture medium, supplemented with human recombinant fibroblast growth factor-basic (FGF2; 100-18C, Peprotech), B27 (17504044, Thermo Fisher Scientific), and other factors (73). Culture medium was replaced every 2 to 3 days, and organoids were passaged every 12 to 14 days. For passaging, organoids were digested in 1 mg/mL collagenase-dispase solution (LS004106, Cedarlane) for 1.5 hours followed by 15 minutes TrypLE (12604021, Thermo Fisher Scientific) treatment. Dissociated organoid cells were washed, resuspended in Cultrex, and seeded at a density of 20,000 cells per well of a 24-well plate. After solidification, 1 mL of culture medium was added to each well.

ROR2 knockout in PANC1 and UM5 cells was induced by using the eSpCas9-LentiCRISPRv2 vector system (74) containing the following guide RNA sequences: ROR2 guide RNA (gRNA) 1 (AGCCGCGGCGGATCATCATC), ROR2 gRNA 2 (ATGAAGACCATTACCGCCAC) or scrambled control gRNA (ACGGAGGCTAAGCGTCGCAA). For the generation of lentiviruses, 293FT cells were transfected with the target vector and the packaging plasmids pLP1, pLP2, and pLP/VSVG (Thermo Fisher Scientific) by using Lipofectamine 2000 (11668019, Thermo Fisher Scientific). Lentiviral particles were collected 72 hours after transfection and concentrated by using Lenti-X Concentrator (631231, Takara). After transduction, PANC1 and UM5 cells were constantly kept in cell culture medium containing 5 and 2.5 µg/mL Puromycin (P8833, Millipore Sigma), respectively. ROR2 knockout was confirmed by Western Blot.

For overexpression of ROR2, pLV[Exp]-EGFP:T2A:Puro-EF1A>hROR2 (VB900004-2980bun, VectorBuilder) or control vector pLV[Exp]-EGFP/Puro-EF1A>ORF_stuffer, (VB900021-9574dpn, VectorBuilder) were packaged into lentiviral particles, concentrated supernatant was added to HPAF-II and UM32 cells, and cells were selected with Puromycin.

One day before siRNA transfection, cells were seeded in antibiotics-free cell culture medium. Cells were transfected with either 10 nmol/L Silencer Select Negative Control #1 (4390843, Thermo Fisher Scientific) or siROR2 #1 (s9760, Thermo Fisher Scientific) resuspended in Opti-MEM I Reduced Serum Medium (31985062, Thermo Fisher Scientific) containing Lipofectamine 2000. After 6 hours of incubation, cells were washed with PBS and cultivated in normal cell culture medium. Cells were harvested 3 and 4 days after transfection for RNA and protein isolation, respectively. Knockdown efficiency was ∼70%.

To estimate cell proliferation, 2,000 cells were seeded on a well of a 96-well plate, and brightfield images were taken every 2 hours by using the BioTek BioSpa 8 Automated Incubator and Citation 5 reader (Agilent). BioTek Gen5 imaging software (Agilent) was used to determine cell confluency, which is generally represented as the mean of four wells.

Spheroid formation was studied according to the clonal assay principle. For each experiment, single cells were seeded in wells of three, U-bottom shaped 96-well ultra-low attachment plates (CLS7007-24EA, Millipore Sigma). Cells were cultivated in serum-free RPMI 1640 medium containing B-27 supplement, 10 ng/mL FGF2, 20 ng/mL human recombinant epidermal growth factor (EGF; AF-100-15, Peprotech), 10 mmol/L HEPES (15-630-106, Thermo Fisher Scientific) and 0.5% Gentamicin. Spheroid cultures were kept for 14 days, number of spheres in relation to total number of seeded wells was assessed, and sphere forming efficiency was calculated. Only spheres with a diameter >75 µm were considered. Brightfield images were taken with the THUNDER Microscope Imager (Leica).

For drug response assays, 4,000 cells/well were seeded in 100 µL of appropriate media in a 96-well plate. After overnight incubation, cells received treatment with increasing concentrations of the KRAS^G12D inhibitor MRTX1133 (E1051, Selleckchem) and AKT inhibitor MK-2206 (S1078, Selleckchem). All vehicles were treated with 1% DMSO. Cell viability was assessed by using CellTiter-Glo 2.0 reagent (G9242, Promega). After 48 to 72 hours, 100 µL CellTiter-Glo 2.0 reagent was added to each well, followed by incubation of the plate on a shaker for 2 minutes using 120 rpm and 10 minutes incubation at room temperature. Luminescence was measured by using the SpectraMax iD3 (Molecular Devices). The average of each duplicate was calculated and normalized to the vehicle control. IC₅₀ values were calculated using GraphPad Prism version 9.

The strains Ptf1aCre^ERT, Pdx1^f/f, Kras^LSL-G12D, and R26R^LSL-YFP were previously described (17). Ror2^f/f animals were purchased from the Jackson Laboratory (strain #018354; ref. 75) and bred into Ptf1aCre^ERT/+;Kras^LSL-G12D animals. To induce Ptf1aCre^ERT-dependent recombination, animals were treated with Tamoxifen (T5648, Sigma). Tamoxifen was dissolved in corn oil and administered by oral gavage, 5 mg daily for a total of 5 days. All animal procedures were approved by the University of Michigan and Henry Ford Health Institutional Care and Use of Animal Committees.

Acinar cells were isolated from pancreas tissue as previously described (17). In brief, pancreas tissue was washed three times in Hank’s Balanced Salt Solution (HBSS; SH3026801, Thermo Fisher Scientific), minced and digested in HBSS containing 0.5 mg/mL Collagenase P (Roche, 11249002001) for 15 minutes at 37°C. Digestion reaction was stopped by washing the tissue three times in HBSS with 5% FBS. Acinar cells were isolated by straining the tissue through 500- and 100-µm filter meshes and centrifugating the cells through a gradient of 30% FBS in HBSS. Acinar cells designated for RNA-sequencing experiments were washed in PBS, resuspended in RLT buffer, and RNA was immediately isolated.

PCR-based genotyping of the Ror2^NULL allele was performed according to Ho and colleagues (75). In brief, acinar cells of Ptf1a^ERT;K* and Ptf1a^ERT;K*;Ror2^f/f mice were isolated 1 week after Tamoxifen administration, and genomic DNA was isolated using the DNeasy Blood and Tissue Kit (69504, Qiagen). PCR was performed using the Go Taq Master Mix (PRM7123, Promega).

For Western Blot analysis, cells were seeded at defined densities and lysed in RIPA buffer (10 mmol/L Tris-HCl (pH 7.4), 150 mmol/L NaCl, 0.1% SDS, 1% Na-deoxycholate and 1% Triton X-100) supplemented with protease (A32965, Thermo Fisher Scientific) and phosphatase inhibitors (4906837001, Millipore Sigma). Next, cell extracts were sonicated for 10 seconds, and residual debris was pelleted by centrifugation (10 minutes, 18,400 g, 4°C). For protein extraction from pancreata, ∼ 30 mg tissue was added to tissue lysis buffer (10 mmol/L Tris-HCl (pH 7.4), 150 mmol/L NaCl, 1 mmol/L EDTA, 1% NP-40, 0.1% SDS, 0.5% Na-deoxycholate supplemented with protease and phosphatase inhibitors). Tissue was homogenized by using the Qiagen tissue lyser (3 minutes, 50 oscillations/s). Samples were sonicated, and supernatant was extracted after centrifugation (20 minutes, 18,400 g, 4°C).

Protein concentration was quantified by using the Pierce BCA Protein Assay Kit (23227, Thermo Fisher Scientific). Protein lysates were denatured at 95°C for 5 minutes, 20 to 30 µg of protein were generally loaded on polyacrylamide gels, and electrophoresis was performed. Separated proteins were transferred onto a PVDF membrane (IPVH00010, Millipore Sigma) by wet-protein transfer. Membranes were blocked in Tris-buffered saline, 0.1% Tween 20 (TBST) containing 5% milk, and primary antibodies as listed in Supplementary Table S1 were applied overnight at 4°C. The secondary, anti-rabbit HRP-conjugated antibody (45-000-682, Thermo Fisher Scientific) was incubated for one hour at RT, and bands were visualized using Clarity or Clarity Max Western ECL substrate (1705061 or 1705062, Bio-Rad) with the ChemiDoc Imaging System (Bio-Rad).

Tissue designated for IHC analysis was fixed overnight in Z-Fix (NC9050753, Thermo Fisher Scientific), washed with PBS for 30 minutes at RT, processed, and embedded in paraffin. The use of human samples was approved by the Henry Ford Health Institutional Review Board. IHC single- and multiplex-staining were performed on 4 µm sections using the Discovery ULTRA stainer (Roche). Antibodies and dilutions are provided in Supplementary Table S1. Hematoxylin and Eosin (H&E) and Picrosirius Red staining were performed according to standard protocols. Stained sections were scanned with the PANNORAMIC SCAN II (3DHistech Ltd.), and the number of positively stained cells/area was quantified with HALO software v3.5.3577.214 (Indica Labs). Representative images were taken from scanned images or with the Stellaris 5 confocal microscope (Leica).

For IF analysis of adherent pancreatic cancer cells, cells were seeded on chamber slides. Cells were fixed with 4% formaldehyde for 10 minutes at RT, washed with PBS, and permeabilized using PBS containing 0.25% Triton X-100 (T8787, Millipore Sigma). After blocking in PBS supplemented with 1% BSA (BP9703, Thermo Fisher Scientific), 5% donkey serum (D9663, Millipore Sigma), and 0.02% Triton X-100, primary antibodies were incubated overnight at 4°C. Antibodies and dilutions are provided in Supplementary Table S1. Then, cells were stained with the secondary antibodies, anti-rabbit 488 (A21206, Thermo Fisher Scientific), and anti-mouse 568 (A10037, Thermo Fisher Scientific) as well as Hoechst 33342 (62249, Thermo Fisher Scientific).

Multiplexed RNA-ISH and IF staining were performed as previously described (76). In brief, FFPE sections were baked in a dry air oven at 65°C for 1 hour, deparaffinized, and hydrated in xylene and 100% ethanol. After antigen retrieval was performed in a steamer for 15 minutes using the Co-Detection Target Retrieval (323165, Advanced Cell Diagnostics) solution and blocking, the primary antibody was applied. Slides were postfixed in 10% Neutral Buffered Formalin for 30 minutes at RT and treated with RNAscopeProtease Plus (322381, Advanced Cell Diagnostics) at 40°C for 11 minutes. For application of RNA-ISH probes, the RNA-ISH Multiplex Fluorescent Assay v2 Kit (323110, Advanced Cell Diagnostics) was used according to manufacturer’s instructions. Following RNAscope staining, sections were stained with DAPI (4′,6-diamidino-2-phenylindole) at RT for 15 minutes, washed in PBS, secondary antibody was applied, and slides were embedded with ProLong Diamond (P36961, Thermo Fisher Scientific) mounting medium.

Pancreatic tissue was snap-frozen in liquid nitrogen and stored at −80°C. For isolation of nuclei, frozen tissue was pulverized by using the CP02 pulverizer (Covaris). Pulver was washed with PBS containing 0.04% BSA and 1× protease inhibitor and centrifuged (5 minutes, 500 g, 4°C). Supernatant was discarded, and pellet was resuspended in 0.1 U/μL DNase solution (EN0521, Thermo Fisher Scientific) and incubated on ice for 5 minutes. Tissue was washed twice with PBS containing 0.04% BSA and 1× protease inhibitor. After the last centrifugation step (5 minutes, 500 g, 4°C), supernatant was discarded, and pellet was resuspended in freshly prepared lysis buffer (10 mmol/L Tris-HCl, 10 mmol/L NaCl, 3 mmol/L MgCl₂, 0.1% Tween-20, 0.1% Nonident P40, 0.01% Digitonin (BN2006, Thermo Fisher Scientific), 1% BSA and 1× protease inhibitor) and incubated on ice for 10 minutes. Then, nuclei were isolated by using a 2 mL Dounce homogenizer, 25 strokes with pestle A and 25 strokes with pestle B. To remove connective tissue and residual debris, the sample was filtrated through a 70 μm cell strainer (352350, Thermo Fisher Scientific), followed by centrifugation (1 minute, 100 g, 4°C). Supernatant was transferred to a new tube and washed with wash buffer (10 mmol/L Tris-HCl, 10 mmol/L NaCl, 3 mmol/L MgCl₂, 0.1% Tween-20, 1% BSA and 1× protease inhibitor) followed by centrifugation (5 minutes, 500 g, 4°C) for three times. Before the final washing step, suspension was filtered through a 40 µm Flowmi cell strainer (BAH136800040-50EA, Millipore Sigma). Finally, nuclei were resuspended in Nuclei Buffer (10× Genomics). Library preparation and sequencing was performed at the Advanced Genomics Core of the University of Michigan. First, single-nucleus suspensions were counted on the LUNA Fx7 Automated Cell Counter (Logos Biosystems). Then, single-nucleus ATAC libraries were generated using the 10× Genomics Chromium instrument following the manufacturer’s protocol. Final library quality was assessed using the LabChip GX (PerkinElmer). Pooled libraries were subjected to paired-end sequencing according to the manufacturer’s protocol (Illumina NovaSeq 6000). 10,000 cells were profiled with a sequencing depth of 50,000 reads per cell. Bcl2fastq2 Conversion Software (Illumina) was used to generate demultiplexed fastq files. Raw sequencing data were processed using cellranger-atac-2.0.0 (10× Genomics), and reads were aligned to mm 10 reference genome.

Doublet removal was performed using scDblFinder (excluding cells with a −log₁₀(q-value); ref. 77) and AMULET (doublet score > 0.975; ref. 78). In addition, nuclei were excluded if the number of reads with low mapping quality was >40,000, percent of mitochondrial reads was >10%, percentage of reads in peaks was <20%, nucleosome signal was >0.63, activity was detectable in <3,200 features or >32,000 features, and TSS enrichment score <2.5. Next, features with <500 counts across all samples or with reads in <1% of cells were excluded. Gene activity scores and features were created using the GeneActivity function from the Signac R package (29). Genes with an activity score of <200 across all cells (per sample) or with activity scores in <1% of cells (per sample) as well as nuclei with activity scores in <1,000 genes were removed. Then, for each nucleus, gene activity scores were normalized, and natural log transformed. Uniform Manifold Approximation and Projection (UMAP) and shared nearest-neighbor graph were constructed using the Seurat R packages RunUMAP and FindNeighbors functions, respectively (79). Dimensional reduction was performed via latent semantic indexing (lsi) using dimensions 2:30. Cell clusters were identified using the FindClusters function with smart local moving algorithm (resolution = 0.4; ref. 80). UMAP projection created a central cluster of low-quality nuclei that demonstrated comparatively poorer QC metrics, such as significantly lower percentage of reads in peaks, lower TSS enrichment, and mapping quality scores. This cluster as well as clusters with <20 cells were removed. Due to a lymph node contamination in one of the 2-week snATAC-seq samples, sequencing of one sample was repeated. The minor batch effect was remedied via ComBat (81) using adjusted cluster labels in the model matrix. UMAP projects were recreated with batch-corrected data as described above. The Seurat R package FindAllMarkers function was used to identify differentially accessible peaks among cell clusters (79). Potential marker genes were required to be present in at least 25% of cells per cluster and show on average, at least a 0.25-fold difference (log-scale) among cell populations. Based on identified cluster-specific marker genes, clusters were manually annotated.

For total RNA extraction, the RNeasy Plus Mini Kit (74134, Qiagen) was used following manufacturer’s instructions. Genomic DNA was digested on column by using the RNase-Free DNase Set (79254, Qiagen). mRNA was reverse transcribed to cDNA with the iScript(tm) cDNA Synthesis Kit (1708891BUN, Bio-Rad). Quantitative real-time PCR (qRT-PCR) was performed on the QuantStudio 6 Pro qRT-PCR cycler (Thermo Fisher Scientific) using 20 ng cDNA, 500 nmol/L forward and reverse primers as well as 1× Fast SYBR Green Master Mix per reaction (4385612, Thermo Fisher Scientific). Primer sequences are provided in Supplementary Table S2. Relative gene expression was determined in relation to expression of the housekeeping genes GAPDH and HPRT1 and calculated according to the ΔΔCt method (82).

Total RNA was at least analyzed in triplicates. Library preparation and sequencing of RNA samples from isolated mouse acinar cells were performed at the University of Michigan Advanced Genomics Core. Libraries underwent paired-end sequencing with 37.5 million reads per sample. Human RNA samples were processed at MedGenome Inc.. RNA-seq library preparation was performed by using the Illumina TruSeq stranded mRNA kit. RNA-seq libraries were paired-end sequenced (PE100/150) with 40 million reads per sample.

Sequencing reads of mouse RNA-seq experiments from isolated acinar cells were assessed for quality and quantity metrics using FastQC and multiQC (83). Reads were aligned to the mouse reference genome (GRCm38/mm 10 Ensembl release 98) using Salmon (84). Next, sample-specific transcripts per kilobase million (TPM) estimates were manually converted to gene-specific TPM estimates by summing transcript TPMs per gene. Genes and their respective transcripts were collected using the tr2g_ensembl function from the BUSpaRse R package. The original sample by gene count matrix was first adjusted for per sample sequencing library inspired by DESeq2 (85). The geometric mean was calculated for each gene across all samples. Gene counts were divided by this geometric mean, resulting in gene count ratios. A sample-specific-median or “sample size factor” of these ratios was finally divided from each gene count per sample. Features with estimates <1 TPM were removed, gene counts were log₂ transformed, and the top 2,000 most variable genes were selected. The gene expression matrix was further normalized by centering and scaling for each gene using [E-mean(E)/sd(E)]. For unsupervised clustering of these genes, the geometric means of the normalized counts per treatment group were used. Cluster analysis was performed using the Self-Organizing Map Program Package (SOM PAK; version 3.1; ref. 86).

For RNA-seq experiments of human pancreatic cancer cell lines, adapter trimming and removal of ribosomal and mitochondrial genome was carried out with fastq-mcf, cutadapt, and Bowtie2 (87), respectively. Reads were mapped to the hg19 genome using STAR (88), and raw read counts were estimated with HTSeq (v0.11.2; ref. 89) and normalized with DESeq2 (85).

In order to compare the snATAC-seq data to bulk RNA-seq data of isolated acinar cells, a sample by gene pseudo-bulk count matrix was created. First, UMAP clusters identified as acinar cells were subset from the snATAC-seq Seurat object. Then, acinar cells were grouped by sample, and the trimmed arithmetic mean of raw counts was calculated to represent the pseudo-bulk counts; before the mean was computed, a fraction of 0.05 of observations was trimmed from each end of per gene counts. The pseudo-bulk sample by gene read count matrix (8 × 20,227) was adjusted for per sample sequencing library inspired by DESeq2 (85) and calculated as described above. Genes with <1 count across all samples were removed, gene counts were log₂ transformed, and the top 2,000 most variable genes were selected. The pseudo-bulk ATAC-seq sample by gene (8 × 2,000) transformed count matrix was compared with the RNA-seq sample by gene (12 × 2,000) transformed count matrix. Both matrices were adjusted for per sample sequencing library and log₂ transformed. Any missing values were imputed with a value of zero. After performing a sample-by-sample spearman correlation and a principal component analysis on the 374 shared genes across the 12 samples of the RNA-seq matrix (biological replicates per genotype n = 3), the most distant samples per experimental condition were removed. Then, both matrices including a total of eight samples were standardized by centering and scaling for each gene using [E-mean(E)]/sd(E). ATAC-seq genes were measured for their correlation (Spearman) to their respective RNA-seq genes. Genes that demonstrated a correlation ≥0.5 were subset (n = 117).

Single-cell gene expression data from mouse pancreas cells were downloaded from Gene Expression Omnibus (GEO; GSE141017, Samples GSM4293545-GSM4293554; ref. 19). Scanpy version 1.9.3 (90) was used for re-processing of the digital count matrices. Counts were normalized to a target sum of 10,000 after masking highly expressed genes, log-transformed, scaled, and the first 50 principal components were used for Leiden clustering and UMAP projection. The PCA space was batch corrected using Harmony (91) with the sample batch as a covariate and the diversity clustering penalty parameter (theta) = 4, and otherwise, default settings were used. Coarse clustering was achieved with a low Leiden resolution score (0.1) and epithelial cell types were selected after screening clusters for expression of the following marker genes: epithelial (Krt8, Krt18, Krt19), endothelial (Cdh5, Pecam1, Plvap), mesenchymal (Col3a1, Dcn, Sparc), and immune (Tyrobp, Cd52, Ptprc). Then, samples were re-processed with the same settings as above. Epithelial cells were again coarsely clustered, and the primary acinar cell cluster was distinguished from ductal-type cells by the presence of the exocrine marker genes (Ctrb1, Prss2, Try5). Meta/neoplastic cells were identified by tdTomato expression. In brief, K-means (K = 2, scikit-learn v1.2.2) clustering was carried out on the log-normalized tdTomato expression levels across all epithelial cells. Cells in the higher expression group that were not present in the predominant acinar cell cluster were assigned as meta/neoplastic cells. To calculate cell type scores, gene lists were compiled from Busslinger and colleagues (34) and scored using the Scanpy function scanpy.tl.score_genes().

Single-cell gene expression data from mouse stomach cells were downloaded from GEO (GSE141017, Sample GSM4293556; ref. 19). Samples were processed as described above, except that coarse clustering was achieved with a resolution score of 0.04, and cells with <3,000 UMIs were removed.

For the reanalysis of the injury-induced ADM dataset (Ma and colleagues: GSE172380; ref. 21), all YFP + cells (GSE172380_Cluster + CelltypeLabel_YFP+_all-samples_QCed.csv.gz) were processed with Scanpy using the same pipeline, except cells with <1,000 UMIs or >20% mitochondrial content were discarded.

For analyzing RNA-seq data from mouse PDAC cell lines, data (GSE107458; ref. 42) were downloaded from GEO and aligned to mm 10. Quality control and normalization was performed, and to assess correlation of Ror2 to the expression of other genes, normalized matrix was subjected to Pearson’s correlation. Top 50 positively and negatively correlated genes were selected for further analysis, and a heatmap depicting mean normalized expression values was generated using pheatmap. Distribution of normalized expression for candidate markers genes with respect to Ror2 were plotted. Moreover, expression of top 50 most significantly enriched genes in pit cell-like ADM from Schlesinger and colleagues (19) data was visualized in a heatmap.

To analyze previously published RNA-seq data of human PDAC samples (GSE93326; ref. 52), processed data were retrieved from GEO, and counts were normalized using the normTransform function from DESeq2. Only micro-dissected tumor samples were selected for further analysis. For PDAC tumor subtype annotation, Moffitt’s (23) top-25 gene signature was employed, and clustering was performed using ConsensusClusterPlus. Correlation analysis was performed as described above. Analyses were performed using R (version 4.3.1).

For the analysis of human PDAC scRNA-seq data, samples from 174 primary tumors were included (EGAS00001002543, GSE154778, GSE155698, GSE156405, GSE194247, GSE202051, GSE205013, GSE211644, phs001840.v1.p1, PRJCA001063; refs. 2, 43–51). Post quality control and data harmonization was performed using Harmony (91), and ductal cells based on expression levels of KRT18 and KRT19 were selected for further analysis. Using single sample GSEA (ref. 92), each cell was scored for the expression of previously defined classical, basal-like, intermediate (27), and the established pit cell-like ADM signature (50 genes). Cell specific scores were converted to z-scores for comparison between signatures, and cell clusters were classified according to their signatures. Using Pearson’s correlation, the correlations between the classical and pit cell-like ADM signature as well as the basal-like cell and pit cell-like ADM signature were assessed.

GSEA was performed using the GSEA software 4.3.2 (93). For analysis of mouse RNA-seq data, normalized read counts that were centered and scaled by feature were used as input and enrichment analysis was run for generated SPEM signature gene set. GSEA for human cell line RNA-seq data was performed with normalized read counts, and data were aligned to the H: hallmark gene sets. Pathway analysis was performed using the Database for Annotation, Visualization, and Integrated Discovery (DAVID; refs. 94, 95). For the generation of GO terms, genes that exhibited a Pearson correlation coefficient >0.4 and <−0.4, and a P value <0.01 were used. EMT score was calculated based on log₂ normalized read counts by using a custom-built function (https://github.com/umahajanatlmu/useful_commands/blob/main/EMT_signature_score.R) incorporating gene signatures from a meta-analysis of 18 independent gene expression studies (53).

Statistical significance was calculated by two-tailed, unpaired Student t test and Mann–Whitney U test as indicated; *, P < 0.05; **, P < 0.01; ***, P < 0.001; ****, P < 0.0001. Statistical analysis was performed with GraphPad Prism 9 software (GraphPad Software Inc.).

The sequencing data generated in this study were deposited in the GEO database and are accessible through GEO Series accession number GSE253600.

S. Benitz reports grants from the Sky Foundation, Inc. and the German Research Foundation (DFG) during the conduct of the study. L. Huang reports a patent for organoid generation (US20170267977A1) pending and licensed to StemCell Technologies. H.C. Crawford reports grants from the Sky Foundation, Inc. during the conduct of the study. No disclosures were reported by the other authors.

The authors would like to thank the University of Michigan Epigenomics Core and Advanced Genomics Core for tissue preparation and sequencing services, respectively. We would like to thank Dr. Ethan Abel for providing packaging plasmids. We would like to thank Drs. Linda C. Samuelson, Samantha Kemp and Stephen Parker for the helpful discussions. This work was supported in part by Michigan State University through computational resources provided by the Institute for Cyber-Enabled Research (ICER). This work was supported by R01CA247516, U01CA224145, and U01CA274154 (to H.C. Crawford). S. Benitz was supported by a postdoctoral fellowship from the German Research Foundation (BE 7224/1-1) and a research grant provided by the Sky Foundation. U.M. Mahajan was supported by the German Research Foundation (MA 4115/3-1) and the German Cancer Aid (70115409). D.J. Salas-Escabillas was supported by 1F31CA271576. N.G. Steele was supported by R00CA263154.