Intratumoral heterogeneity and cellular plasticity have emerged as hallmarks of cancer, including pancreatic ductal adenocarcinoma (PDAC). As PDAC portends a dire prognosis, a better understanding of the mechanisms underpinning cellular diversity in PDAC is crucial. Here, we investigated the cellular heterogeneity of PDAC cancer cells across a range of in vitro and in vivo growth conditions using single-cell genomics. Heterogeneity contracted significantly in two-dimensional and three-dimensional cell culture models but was restored upon orthotopic transplantation. Orthotopic transplants reproducibly acquired cell states identified in autochthonous PDAC tumors, including a basal state exhibiting coexpression and coaccessibility of epithelial and mesenchymal genes. Lineage tracing combined with single-cell transcriptomics revealed that basal cells display high plasticity in situ. This work defines the impact of cellular growth conditions on phenotypic diversity and uncovers a highly plastic cell state with the capacity to facilitate state transitions and promote intratumoral heterogeneity in PDAC.

Significance:

This work provides important insights into how different model systems of pancreatic ductal adenocarcinoma mold the phenotypic space of cancer cells, highlighting the power of in vivo models.

Pancreatic ductal adenocarcinoma (PDAC) has a dismal prognosis, with a 5-year overall survival of 10% (1). Cancer cell plasticity—the capacity to undergo cell state transitions via nongenetic differentiation mechanisms—promotes tumor progression and is a major cause of treatment failure (2, 3). Plasticity is a requisite for the emergence and maintenance of intratumoral heterogeneity, the property of tumors to contain discrete cell states with distinct functional properties such as states with intrinsic or acquired resistance to therapy (3, 4). Thus, understanding the heterogeneity and plasticity of cancer cell states in PDAC and in other cancers is of great biological and clinical significance.

Extensive genomic characterization of PDAC has revealed a largely homogenous mutational landscape, with recurrent alterations in KRAS, TP53, CDKN2A, and SMAD4 as the dominant oncogenic and tumor suppressor alterations (5). In contrast to the relatively uniform genomic landscape, PDAC tumors exhibit considerable inter- and intratumoral transcriptional diversity. The intertumoral heterogeneity of PDAC may be in part due to differences in the cell of origin (6), the underlying tumor epigenetic landscape (7), or specific tumor suppressors (8). Multiple transcriptional subtypes derived from either bulk or microdissected tumor samples have been described (9–12). The most reproducible molecular subtypes identified across different classifications are a basal-like and a classical subtype; these expression signatures also correlate with prognosis and may serve as important clinical biomarkers (13, 14). Furthermore, a subset of PDAC cells undergoes epithelial–mesenchymal transition (EMT), which promotes tumor progression and treatment resistance (15–17). More recently, several studies have leveraged single-cell RNA sequencing (scRNA-seq) to demonstrate that the classical, basal, and mesenchymal cell states can coexist within individual tumors (18–21). Little is known about the lineage relationships and plasticity of these transcriptionally defined cell states in PDAC tumors in situ.

GEMMs provide a unique platform for functional interrogation of cancer cell states in situ in autochthonous tumors that, like human cancers, evolve in the context of the relevant tissue site and a functional immune system (22). In commonly used PDAC GEMMs, Cre or Flp recombinase is placed under the control of the pancreatic epithelial progenitor–specific Pdx1 promoter allowing somatic activation of oncogenic KrasG12D and directed perturbation of Trp53 (herein KPC and KPF mice, respectively) in pancreatic epithelial cells (23–25). These models mimic KRAS-driven human PDAC at molecular and histopathologic levels, including key hallmarks such as desmoplastic stroma, an immunosuppressive microenvironment, and cancer cell state heterogeneity (23, 24, 26, 27). Given their widespread use, understanding the extent to which these PDAC GEMMs recapitulate cell state heterogeneity of human tumors is of great importance.

In addition to autochthonous GEMMs, PDAC is commonly studied in two-dimensional (2D) and three-dimensional (3D) organoid cell culture conditions and by orthotopic transplantation of cell lines or tumor fragments (28, 29). The degree to which these models recapitulate the plasticity and heterogeneity of primary autochthonous tumors is incompletely understood. Elegant recent work by Raghavan and colleagues demonstrated that human PDAC cells isolated from primary tumors or metastases undergo dramatic phenotypic shifts in ex vivo culture conditions (20). Interestingly, the cell states acquired in culture shifted in response to changes in culture conditions and to the addition or withdrawal of specific growth factors, suggesting that PDAC cells are plastic in vitro (20). However, whether PDAC cells propagated in these in vitro culture systems are capable of reconstituting the heterogeneity of autochthonous tumors upon transplantation is unknown. This is of particular importance, given that 2D and 3D culture systems are often necessary intermediates for expanding the number of cells or incorporating candidate gene perturbation systems, reporters, or other genetic elements into cells. These questions have broad significance beyond PDAC; yet, to our knowledge, a systematic comparison of cell state heterogeneity across autochthonous tumors and cell culture conditions, as well as transplants derived from these growth conditions, has not been performed in any cancer context.

Animal studies

All animal studies were approved by the MSKCC Institutional Animal Care and Use Committee (protocol # 17-11-008). MSKCC guidelines for the proper and humane use of animals in biomedical research were followed. All genetically engineered mice were maintained on C57/black 6 and Sv129 mixed backgrounds. Autochthonous KPCT mice were generated by crossing previously published KrasLSL-G12D/+ (30), Trp53flox/flox (31), Pdx1-Cre (23), and Rosa26LSL-tdTomato mice (32). For Lgr5 coexpression studies, the KPCT mice above were crossed with Lgr5GFP-IRES-CreER/+ reporter mice (33). Basal-cell state lineage tracing was performed by generating autochthonous tumors by crossing previously published KrasFSF-G12D/+ (34), Trp53frt/frt (35), Pdx1FlpO/+ (36), Rosa26mTmG/+ (37), and Lgr5CreER (38) mice to generate KPF; RosamTmG/+; Lgr5CreER/CreER mice. Mice were genotyped at ∼2 weeks of age and tumor-bearing mice euthanized for tumor harvest at 7 to 9 weeks of age when demonstrating PDAC symptoms (distended abdomen, weight loss, hunching).

Orthotopic models

NOD–SCID gamma (NSG) or first-generation (F1) C57/black 6 × Sv129 hybrid mice were used as transplant recipients. The following types of transplants were generated: ∼3-mm thick PDAC fragments generated from KPCT primary tumors following previously described method (39); or by implantation of tumor cells derived from KPCT primary PDAC tumors following a previously described method (40). For the latter, either 10,000 unsorted tumor cells, 10,000 FACS isolated tdTomato+ tumor cells, or 200,000 cells from primary murine 2D or 3D (Matrigel) cell cultures collected after eight passages were orthotopically implanted into the tail of the pancreas. Presence of tumors was confirmed by small animal high-resolution ultrasound imaging (Vevo 2100, Visual Sonics). Animals were monitored weekly, and tumors were harvested when they reached 1 cm3 or if the animal fulfilled criteria for humane euthanasia.

Lineage tracing

Autochthonous KPF; Rosa26mTmG/+; Lgr5CreER/CreER tumor-bearing mice were genotyped and monitored as above. For autochthonous lineage-tracing experiments, tumor-bearing mice were treated with a single dose of tamoxifen (200 mg/kg via oral gavage) at 6.5 weeks of age and tumors were harvested either 3 days or 10 days after tamoxifen administration. For fragment-based lineage-tracing experiments, tumors were harvested from autochthonous KPF; Rosa26mTmG/+; Lgr5CreER/CreER tumor-bearing NSG mice. Tumor implantation and growth was confirmed by ultrasound 7 and 14 days after transplantation, at which point, fragment-bearing mice were administered a single dose of tamoxifen (200 mg/kg via oral gavage). Tumors were harvested at 3, 10, and 20 days following tamoxifen administration. Day 3 was chosen as the baseline time point to account for conversion of tamoxifen to its active metabolite, 4-hydroxytamoxifen (4-OHT), to allow time for successful recombination and expression of GFP in lineage-traced cells, and washout of 4-OHT. Tamoxifen was dissolved in corn oil slurry at 20 mg/mL and incubated at 55°C for 1 hour until fully dissolved.

Isolation of tumor cells

To dissociate tumors into single-cell suspensions, primary KPCT and KPF; RosamTmG/+; Lgr5CreER/CreER tumors were finely chopped with scissors and incubated with digestion buffer containing collagenase IV (17104019, ThermoFisher Scientific, 0.1 U/mL), Dispase (#354235, Corning, 0.6 U/mL), DNase I (#69182–3; Sigma Aldrich, 10 U/mL), and Soybean Trypsin Inhibitor (#T9003, Sigma Aldrich, 0.1 mg/mL) dissolved in HBSS with Mg2+ and Ca2+ (#14025076, Thermo Fisher Scientific) in gentleMACS C Tubes (Miltenyi Biotec) for 42 minutes at 37°C using the gentleMACS Octo Dissociator (#130–096–427, Miltenyi Biotech). Following enzymatic dissociation, samples were washed with HBSS + Mg2+ and Ca2+ and filtered through a 100 μmol/L cell strainer and spun at 300 g for 5 minutes at room temperature. Cells were then washed with media and pelleted at 300 g for 5 minutes at 4°C. The supernatant was removed, and the pellet resuspended in FACS buffer media (200 mmol/L EDTA with 2% of heat-inactivated FBS in PBS) before being passed through a 40-μm strainer. At this point cells were used for direct orthotopic transplantation, for generation of 2D and 3D cell lines, or were prepared for FACS. For FACS staining, cell suspensions were blocked for 5 minutes at room temperature with rat anti-mouse CD16/CD32 (Mouse BD Fc Block, #553142, BD Biosciences) in FACS buffer, and incubated for 30 minutes with a mix of four APC-conjugated antibodies binding CD45 (#559864, BD Biosciences, 1:500), CD31 (#561814, BD Biosciences, 1:500), CD11b (#561690, BD Biosciences, 1:500), and TER-119 (#561033, BD Biosciences, 1:500) as well as 300 nmol/L DAPI as a live-cell marker. Sorting was performed on a BD FACSAria sorter (BD Biosciences) for either (CD45/CD31/CD11b/TER119)/tdTomato+/DAPI live cells (KPCT mice) and/or for (CD45/CD31/CD11b/TER119)/tdTomato+/GFP+/DAPI live cells (KPF; RosamTmG/+; Lgr5CreER mice) in lineage-tracing experiments. Sorted cells were collected directly into S-MEM (1×, Gibco) with 2% of heat-inactivated FBS for scRNA-seq or for single-cell assay for transposase-accessible chromatin using sequencing (scATAC-seq). For scRNA-seq, individual tumor cell suspensions were incubated for 30 minutes with hashtag oligonucleotide-conjugated antibodies (see Supplementary Data) in addition to FACS antibodies.

Droplet-based scRNA-seq

PDAC tumors were harvested, dissociated, and tumor cells isolated by FACS as described above. The scRNA-seq of FACS-sorted cell suspensions was performed using the 10X Genomics Chromium platform according to user guide manual for 3′ v3. Briefly, FACS-sorted cells were washed once with PBS containing 1% BSA and resuspended in PBS containing 1% BSA to a final concentration of 700 to 1,300 cells per μL. The viability of cells in all experiments was above 80%, as confirmed with 0.2% (w/v) Trypan Blue staining (Countess II). Cells were captured in droplets. Following reverse transcription and cell barcoding in droplets, emulsions were broken and cDNA purified using Dynabeads MyOne SILANE, followed by PCR amplification as per manufacturer's instructions. Between 20,000 to 25,000 cells were targeted for each sample. Samples were multiplexed together in droplet formulation lanes following the TotalSeq B cell hashing protocol (41). Final libraries were sequenced on Illumina NovaSeq S4 platform (R1 – 28 cycles, i7 – 8 cycles, R2 – 90 cycles). Detailed computational analysis is described in the Supplementary Data.

ATAC-seq

Cancer cells were isolated from autochthonous KPCT PDAC tumors by flow cytometry as above, frozen in Bambanker Cell Freezing Medium (Lymphotec, catalog no. 302–14681), and stored in a liquid nitrogen freezer until use. Single-cell libraries were prepared using 10x Genomics scATAC-seq protocol from 10x Genomics Chromium [Single-Cell ATAC Reagent Kits User Guide (CG000168, Rev A - v1 chemistry)]. Briefly, cells were centrifuged (300 g, 5 minutes, 4°C) and permeabilized with 100 μL chilled lysis buffer (10 mmol/L Tris-HCl pH 7.4, 10 mmol/L NaCl, 3 mmol/L MgCl2, 0.1% Tween-20, 0.1% IGEPAL-CA630, 0.01% digitonin, and 1% BSA). The samples were incubated on ice for 3 to 5 minutes and resuspended in 1 mL chilled wash buffer (10 mmol/L Tris-HCl pH 7.4, 10 mmol/L NaCl, 3 mmol/L MgCl2, 0.1% Tween-20 and 1% BSA). After centrifugation (500 g, 5 minutes, 4°C), the pellets were resuspended in 100 μL chilled Nuclei buffer (2000153, 10× Genomics). Nuclei were counted using a hemocytometer, and the nucleus concentration was adjusted to 3,000 nuclei per μL. We used 15,360 nuclei as input for tagmentation. Nuclei were diluted to 5 μL with 1× nuclei buffer (10x Genomics) and mixed with ATAC buffer (10x Genomics) and ATAC enzyme (10x Genomics) for tagmentation (60 minutes, 37°C). Single-cell ATAC-seq libraries were generated using the Chromium Chip E Single Cell ATAC kit (10× Genomics, 1000086) and indexes (Chromium i7 Multiplex Kit N, Set A, 10× Genomics, 1000084) following the manufacturer's instructions. Final libraries were verified using a TapeStation (Agilent). Libraries were sequenced on a NovaqSeq6000 (Illumina) with the following read lengths: 50 + 8 + 16 + 50 (Read1 + Index1 + Index2 + Read2). Detailed computational analysis is described in the Supplementary Data.

Cell culture and reagents

Primary 2D mouse cell lines were generated by harvesting and digesting autochthonous KPCT PDAC tumors as described above, and plating 100,000 cells in advanced DMEM/F12 (Gibco) supplemented with 2 mmol/L glutamine, 500 ng/mL pen/strep, and 2% heat-inactivated FBS. After 24 hours, the media was changed to remove non-adherent cells and the media was changed to advanced DMEM/F12 with 2 mmol/L glutamine, 500 ng/mL pen/strep, and 10% FBS. Cells were passaged using 0.05% trypsin-EDTA (Gibco).

Primary 3D mouse cell lines were generated as previously described (28). In brief, dissociated cells were seeded in growth-factor-reduced Matrigel (BD) and grown in enriched medium [Advanced DMEM/F12 medium supplemented with HEPES (1×, Invitrogen), Glutamax (1×, Invitrogen), penicillin/streptomycin (1×, Invitrogen), B27 (1×, Invitrogen), N-acetyl-L-cysteine (1 mmol/L, Sigma), RSPO1-conditioned medium (10% v/v, a kind gift from Calvin Kuo), Noggin recombinant protein (0.1 μg/mL, Peprotech), epidermal growth factor (EGF, 50 ng/mL, Peprotech), Gastrin (10 nmol/L, Sigma), fibroblast growth factor 10 (FGF10, 100 ng/mL, Preprotech), and nicotinamide (10 mmol/L, Sigma)]. For a subset of experiments, organoids were grown in minimal medium [Advanced DMEM/F12 medium supplemented with HEPES (1×, Invitrogen), Glutamax (1×, Invitrogen), penicillin/streptomycin (1×, Invitrogen), and 2% heat-inactivated fetal bovine serum].

Tissue histology and immunofluorescence

Tumors were harvested at the experimental endpoints, fixed in 10% neutral buffered formalin (Richard-Allan Scientific), embedded in paraffin, and cut into 5-μm sections. Slides were heated for 30 minutes at 60°C, deparaffinized, rehydrated with an alcohol series and incubated in Tris-EDTA antigen retrieval buffer (#E1161, Sigma-Aldrich) for 20 minutes in a pressure cooker on 95°C. Sections were washed in PBS, permeabilized in 0.3% PBS-Triton X-100 for 40 minutes, then blocked in PBS/0.1% Triton X-100 containing 2% BSA and 5% donkey serum (#D9663, Sigma-Aldrich). Primary antibodies were incubated overnight at 4°C in blocking buffer. The following primary antibodies were used: RFP (#6g6, ChromTech, 1:500), GFP (#ab5450, Abcam, 1:500), vimentin (#ab92547, Abcam, 1:500), cytokeratin 17 (#4543, Cell Signaling Technology, 1:250), galectin-4 (#PA5-34913, Thermo Scientific, 1:500). After washing with PBS, tissues were incubated with secondary antibodies (#A-11055, #A-10042, #A-31571, ThermoFisher, 1:500) for 1 hour at room temperature. After staining slides were counterstained with DAPI (#D9542, Sigma-Aldrich, 5 μg/mL) for 10 minutes and coverslipped with Mowiol mounting reagent. Hematoxylin and eosin were performed using standard protocols. Images were acquired using 20× or 40× objective on a Zeiss Axio Imager Z2 and ZEN 2.3 software.

In situ hybridization

Single-molecule mRNA in situ hybridization was performed on formalin-fixed paraffin-embedded tissues using the Advanced Cell Diagnostics RNAscope 2.5 HD Detection Kit (#322360, ACD) and probe Hs-Lgr5 (#312178, ACD) according to the manufacturer's instructions.

Data availability

The data generated in this study are publicly available in Gene Expression Omnibus at GSE209599.

To systematically interrogate PDAC cell state heterogeneity across distinct growth conditions, we sought to employ a model that is genetically defined and incorporates a cell state–agnostic marker for isolation of cancer cells. Thus, we generated KrasLSL-G12D/+; Trp53flox/flox; Pdx1-Cre; Rosa26LSL-tdTomato/+ (KPCT) mice, which reproducibly develop PDAC by 6–8 weeks of age and where PDAC cells are marked by expression of tdTomato fluorescent protein (Fig. 1A; refs. 24, 42). This is in contrast to two other commonly used KPC models employing either Trp53flox/+ or Trp53R172H/+ allelic configurations, where time to adenocarcinoma transition is significantly more variable than in the Trp53flox/flox model (24). We isolated the tdTomato+ PDAC cells from 15 tumor-bearing mice (aged 7–8 weeks) and carried out an unsupervised survey of PDAC cell state heterogeneity using droplet-based scRNA-seq (Fig. 1A). A total of 14,392 cells passed stringent quality-control standards (Fig. 1B; Supplementary Fig. S1; and Methods). Across mice, we observed a reproducible continuum of transcriptional cell states along an epithelial–mesenchymal axis (Fig. 1B; Supplementary Fig. S2A). Unsupervised clustering identified three higher-order cell states: a classical epithelial cell state, a basal epithelial cell state, and a mesenchymal cell state (Fig. 1B). We found our classical classification was most similar to previously described human classical signatures identified by Moffit and colleagues, Collisson and colleagues, and Raghavan and colleagues (Fig. 1C; Supplementary Fig. S2B; refs. 10, 12, 20). Similarly, there was significant overlap with the murine basal cell state and previously defined human basal signatures reported by Moffit and colleagues and Raghavan and colleagues (12, 20). Finally, the mouse mesenchymal state showed concordance with human quasi-mesenchymal (10) and squamous (9) signatures (Supplementary Fig. S2B).

Figure 1.

Single-cell profiling of autochthonous murine PDAC reveals a continuum of transcriptional cell states along the epithelial–mesenchymal transitional axis. A, Experimental workflow. B, Uniform Manifold Approximation and Projection (UMAP) embedding of scRNA-seq profiles of 14,392 cells from 15 independent tumors, classified as classical (purple), basal (green), or mesenchymal (orange) cells. C, UMAP projection of previously described classical and basal expression profiles of scRNA-seq data. Scale represents gene signature score, with color scale from 1st to 99th percentile. D, Ternary plot depicting the cell state classification probability of each cell, calculated by a Markov absorption classifier trained on a small subset of classical, basal, and mesenchymal cells, colored by subtype (purple, classical; green, basal; orange, mesenchymal). E, Diffusion-based pseudotemporal ordering of cells along the classical–basal–mesenchymal axis. F, Expression of specific classical, basal, and mesenchymal gene signatures along the pseudotemporal axis. Columns represent the smoothed (sliding window smoother with n = 250 cells) expression of individual genes along the pseudotemporal axis. The panel on the leftmost side shows the fractions of classical (purple), basal (green), and mesenchymal (orange) cells along the same pseudotime axis using the same smoothing approach. G, Left, representative UMAP embedding displaying the expression of candidate marker genes, Lgals4 (classical), Krt17 (basal), and Vim (mesenchymal). Scale represents size factor normalized log2-transformed UMI counts, with color scale from 1st to 99th percentile. Right, ×400, representative immunofluorescence images from KPCT PDAC tumors showing galectin-4, CK17, and vimentin in PDAC cells (green in respective panels, from left to right); red, tdTomato+ cancer cells. H, Schematic depicting model of axis of heterogeneity in PDAC.

Figure 1.

Single-cell profiling of autochthonous murine PDAC reveals a continuum of transcriptional cell states along the epithelial–mesenchymal transitional axis. A, Experimental workflow. B, Uniform Manifold Approximation and Projection (UMAP) embedding of scRNA-seq profiles of 14,392 cells from 15 independent tumors, classified as classical (purple), basal (green), or mesenchymal (orange) cells. C, UMAP projection of previously described classical and basal expression profiles of scRNA-seq data. Scale represents gene signature score, with color scale from 1st to 99th percentile. D, Ternary plot depicting the cell state classification probability of each cell, calculated by a Markov absorption classifier trained on a small subset of classical, basal, and mesenchymal cells, colored by subtype (purple, classical; green, basal; orange, mesenchymal). E, Diffusion-based pseudotemporal ordering of cells along the classical–basal–mesenchymal axis. F, Expression of specific classical, basal, and mesenchymal gene signatures along the pseudotemporal axis. Columns represent the smoothed (sliding window smoother with n = 250 cells) expression of individual genes along the pseudotemporal axis. The panel on the leftmost side shows the fractions of classical (purple), basal (green), and mesenchymal (orange) cells along the same pseudotime axis using the same smoothing approach. G, Left, representative UMAP embedding displaying the expression of candidate marker genes, Lgals4 (classical), Krt17 (basal), and Vim (mesenchymal). Scale represents size factor normalized log2-transformed UMI counts, with color scale from 1st to 99th percentile. Right, ×400, representative immunofluorescence images from KPCT PDAC tumors showing galectin-4, CK17, and vimentin in PDAC cells (green in respective panels, from left to right); red, tdTomato+ cancer cells. H, Schematic depicting model of axis of heterogeneity in PDAC.

Close modal

The classical, basal, and mesenchymal cell states do not represent discrete populations. Instead, a Markov absorption-based classifier trained on a subset of cells (43, 44) revealed a predominance of intermediate or mixed cell state identities (Fig. 1D). Interestingly, while we observed individual cells with mixed classical and basal phenotypes, and cells with mixed basal and mesenchymal phenotypes, we did not observe cells with mixed classical and mesenchymal transcriptional states (Fig. 1D). To systematically analyze gene expression changes along the epithelial–mesenchymal axis, we created an unsupervised pseudotemporal ordering of PDAC cells (Fig. 1E). We identified 3969 genes differentially expressed along pseudotime, which could be divided into four main gene signatures broadly corresponding with the classical, basal, and mesenchymal cell states (Supplementary Fig. S2C and S2D; Supplementary Table S1). Gene signature analysis revealed that the classical epithelial cell state was associated with markers of differentiated pancreatic epithelium, such as mucins and glycosylation enzymes, whereas the basal cells were associated with increased keratin and laminin expression. Mesenchymal cells were characterized by the upregulation of genes associated with the extracellular matrix and EMT (Fig. 1F; Supplementary Fig. S2E).

To localize these cell states in their tissue context, we identified a set of marker genes selectively expressed in each state and colocalized them with the tdTomato+ cancer cells in KPCT PDAC tissues (Fig. 1G; Supplementary Fig. S2D and S2F). Immunofluorescence staining revealed the presence of the classical cells marked by galectin-4 (Lgals4) with a columnar epithelial morphology. In contrast, basal epithelial cells expressing cytokeratin 17 (Krt17) tended to show more pronounced basement membrane contact. As expected, most cells expressing the mesenchymal marker vimentin (Vim) showed a mesenchymal or sarcomatoid morphology (Fig. 1G).

Collectively, our data suggest that autochthonous mouse PDAC tumors are composed of three higher-order cell states: a highly differentiated cell state characterized by the expression of classical and ductal markers (classical state), an intermediate epithelial cell state expressing genes linked to basement membrane adhesion (basal state), and a mesenchymal cell state characterized by the expression of extracellular matrix genes (mesenchymal state). Notably, these cell states do not form discrete populations. Instead, the PDAC cells form a continuum of transcriptional states along an epithelial–mesenchymal axis, which defines cell state heterogeneity in KPCT PDAC tumors (Fig. 1H).

The autochthonous KPC model is powerful, but the interrogation of candidate mechanisms in sufficiently large cohorts of mice requires cost-, labor-, and time-intensive breeding. As such, multiple cell-based model systems have been derived from the KPC model, which enable tractable and rapid perturbation of candidates using genetic or pharmacologic approaches. To investigate how these models recapitulate heterogeneity observed in the autochthonous KPC model, we harvested autochthonous tumors and processed them in different ways for orthotopic transplantation into pancreases of immunocompromised NSG mice. These include (i) tumor fragments surgically attached onto the host pancreas, (ii) dissociated cancer + stromal cell mix obtained after enzymatic tumor digestion, (iii) tdTomato+ tumor cells isolated by FACS from the cancer + stromal cell mix, and primary, (iv) 2D, and (v) 3D organoid cells lines derived from the model. All transplant tumors were enzymatically digested and tdTomato+ PDAC cells were isolated by FACS, followed by scRNA-seq (Fig. 2A). To facilitate direct comparison with autochthonous tumors, we projected cancer cells isolated from the transplants into the autochthonous tumor gene expression space (Fig. 2B; Supplementary Fig. S1; and Materials and Methods). In parallel, we used the Markov absorption-based classifier (Fig. 1C) to measure the cell state composition of our transplants (Fig. 2B and C; Supplementary Fig. S1; Materials and Methods).

Figure 2.

The fidelity of cancer cell states is dependent on the model system. A, Orthotopic modeling workflow. Note that in addition to orthotopic tumors, 2D and 3D cell cultures were also processed independently for scRNA-seq directly from culture conditions (red dashed box). B, UMAP embedding of PDAC cells from autochthonous and orthotopic models after projection into the autochthonous gene expression space. Cells are colored according to subtype (purple, classical; green, basal; orange, mesenchymal) based on Markov adsorption classification. C, Representative cellular heterogeneity in each PDAC model system. First row, PDAC tissues stained with hematoxylin and eosin staining at ×200 magnification. Scale bar, 100 μm. Second row, UMAP embedding of each model system (colored dots) onto the full data set (gray dots) depicting the distribution of individual cells throughout phenotypic space. Third row, ternary plots of predicted subtype probability per cell, calculated by Markov absorption colored by subtype (purple, classical; green, basal; orange, mesenchymal). Note that all models recapitulate all three cell states to some degree. Fourth row, ternary plots depicting the overall subtype classification of each individual tumor sample (dots). Note some inherent sample-to-sample (intertumoral) heterogeneity in each model system. D, Box plot of phenotypic volume based on unbiased global transcriptional heterogeneity of each model system calculated. Vertical black line represents the mean phenotypic volume of the autochthonous samples as reference. Higher values are associated with increased spread of cells over phenotypic space. Each dot represents an individual sample. E, Box plot of the compositional heterogeneity based the relative frequency of classical, basal, and mesenchymal cells in each tumor. Values represent the Shannon entropy of the fractions of classical, basal, and mesenchymal cells. Higher values indicate a more balanced composition of cell states. Vertical black line represents the mean phenotypic volume of the autochthonous samples as reference. F, Paired phenotypic heterogeneity in 2D and 3D organoid cell culture and orthotopic transplant. Left, ternary plots of individual cells isolated from 2D (top) and 3D (bottom) organoid cell culture. Right, ternary plots of cells isolated from orthotopic tumors derived from the same cell culture conditions displayed on the left. Within each panel, the top plots display the cell state classification of all cells per model and their distribution in combined UMAP space. Note the in vitro enrichment of epithelial phenotypes and the robust reestablishment of all cell states after orthotopic transplantation. G, Box plots of the phenotypic volume (top) and compositional entropy (bottom) of cell in in vitro cell culture conditions and their matched orthotopic transplants.

Figure 2.

The fidelity of cancer cell states is dependent on the model system. A, Orthotopic modeling workflow. Note that in addition to orthotopic tumors, 2D and 3D cell cultures were also processed independently for scRNA-seq directly from culture conditions (red dashed box). B, UMAP embedding of PDAC cells from autochthonous and orthotopic models after projection into the autochthonous gene expression space. Cells are colored according to subtype (purple, classical; green, basal; orange, mesenchymal) based on Markov adsorption classification. C, Representative cellular heterogeneity in each PDAC model system. First row, PDAC tissues stained with hematoxylin and eosin staining at ×200 magnification. Scale bar, 100 μm. Second row, UMAP embedding of each model system (colored dots) onto the full data set (gray dots) depicting the distribution of individual cells throughout phenotypic space. Third row, ternary plots of predicted subtype probability per cell, calculated by Markov absorption colored by subtype (purple, classical; green, basal; orange, mesenchymal). Note that all models recapitulate all three cell states to some degree. Fourth row, ternary plots depicting the overall subtype classification of each individual tumor sample (dots). Note some inherent sample-to-sample (intertumoral) heterogeneity in each model system. D, Box plot of phenotypic volume based on unbiased global transcriptional heterogeneity of each model system calculated. Vertical black line represents the mean phenotypic volume of the autochthonous samples as reference. Higher values are associated with increased spread of cells over phenotypic space. Each dot represents an individual sample. E, Box plot of the compositional heterogeneity based the relative frequency of classical, basal, and mesenchymal cells in each tumor. Values represent the Shannon entropy of the fractions of classical, basal, and mesenchymal cells. Higher values indicate a more balanced composition of cell states. Vertical black line represents the mean phenotypic volume of the autochthonous samples as reference. F, Paired phenotypic heterogeneity in 2D and 3D organoid cell culture and orthotopic transplant. Left, ternary plots of individual cells isolated from 2D (top) and 3D (bottom) organoid cell culture. Right, ternary plots of cells isolated from orthotopic tumors derived from the same cell culture conditions displayed on the left. Within each panel, the top plots display the cell state classification of all cells per model and their distribution in combined UMAP space. Note the in vitro enrichment of epithelial phenotypes and the robust reestablishment of all cell states after orthotopic transplantation. G, Box plots of the phenotypic volume (top) and compositional entropy (bottom) of cell in in vitro cell culture conditions and their matched orthotopic transplants.

Close modal

A granular analysis of the orthotopic transplants using the cell state classification shows that all of the orthotopic models recapitulate the continuum of classical, basal, and mesenchymal cells states observed in the autochthonous model (Fig. 2C, 2nd and 3rd rows). Concordantly, we observed both epithelial and mesenchymal cell morphologies in all tumor transplant models (Fig. 2C, 1st row) as well as the expression of classical, basal, and mesenchymal marker genes (Fig. 2C, 3rd row; Supplementary Fig. S3A). Interestingly, cell states in the fragment-based and unsorted tumor digest transplant models skewed towards the mesenchymal phenotype (Fig. 2C, bottom; Supplementary Fig. S3B and S3C), suggesting that inclusion of stromal cell components influences cell state. To quantify intratumoral heterogeneity, we calculated the phenotypic volume and compositional entropy of tumors. Phenotypic volume serves as an unsupervised measure of the spread of cells over gene expression space (45). Meanwhile, compositional or cell state entropy shows the degree at which the classical, basal, and mesenchymal cells are represented within a tumor, using the mean cell state distribution in the autochthonous tumors as a reference. We observed a high degree of intratumoral heterogeneity in both autochthonous and orthotopic models (Fig. 2D and E). Globally, most models displayed similar phenotypic diversity as the autochthonous KPCT model with the exception of fragment transplants, which were not as phenotypically complex based on phenotypic volume comparison (Fig. 2D). Similar results were seen when examining the compositional entropy of each model system based on predefined classical, basal, and mesenchymal cell states according to the Markov classification (Fig. 2E). We compared the gene expression of each cell state in the autochthonous model to the various orthotopic models (Supplementary Fig. S4A and S4B; Supplementary Table S2). While we were able to identify differentially expressed genes in each cell state, the overall heterogeneity of the cell states was still predominantly explained by the continuum across the classical, basal, and mesenchymal cell states (Supplementary Fig. S4B). In sum, our data indicates that, similar to autochthonous tumors, orthotopic transplants are composed of a complex, heterogeneous mixture of cell states along a classical–basal–mesenchymal axis.

Establishing primary cell lines from dissociated tumors exerts selective pressure and a bottleneck effect, which may drastically influence phenotypic diversity. Yet given their low cost and tractability, 2D and 3D organoid culture models are widely used in the field to investigate PDAC biology. To investigate heterogeneity in these models, we performed scRNA-seq on tumor cell lines established in either 2D adherent or 3D organoid culture conditions. Consistent with the selective pressure associated with cell culture conditions, we observed markedly reduced heterogeneity in the in vitro conditions: KPCT cells grown in 2D culture were predominantly basal with rare classical cells (Fig. 2F, top left; Supplementary Fig. S3B and S3C), whereas 3D organoid cultures demonstrated the inverse (Fig. 2F, bottom left; Supplementary Fig. S3B and S3C). Notably, despite these classifications, cells in 2D and 3D culture did not exhibit the fully classical or basal extremes observed in vivo, but were rather concentrated at the center of the ternary plots (Fig. 2F, left). This was accompanied by reduced expression of the classical and basal markers in culture (Supplementary Fig. S3A). Of note, we did not observe significant differences in cell state heterogeneity when comparing media composition reported by Boj and colleagues (28) compared with low-serum medium not supplemented with specific growth factors (Supplementary Fig. S5A; ref. 20). In sum, our results indicate that in vitro culture conditions fail to recapitulate in vivo PDAC cell states.

Strikingly, the constrained phenotypic space of the in vitro culture conditions dramatically expanded upon orthotopic transplantation, with reemergence of the full spectrum of classical, basal, and mesenchymal phenotypes (Fig. 2F, right) as well as significant increases in both phenotypic volume and compositional entropy (Fig. 2G). This included the emergence of all three phenotypic extremes and robust re-expression of defined cell state markers (Fig. 2F, right; Supplementary Fig. S3A). Our data shows that the classical–basal–mesenchymal axis, the main axis of heterogeneity in autochthonous KPCT tumors, is established in orthotopic transplants. However, we observed that the gene expression distance between autochthonous and transplant cells consistently exceeds the internal variability of autochthonous cells (Supplementary Fig. S5B), indicating that other sources of heterogeneity distinguish autochthonous tumors from transplants. Specifically, the KPCT cells in autochthonous tumors display a higher number of differentially expressed genes across the three cell states when compared with tumor transplants. Many of these genes, especially in classical and basal cells, come from a branching acinar lineage in the autochthonous tumors or are related to ribosomal function (Supplementary Table S2 and not shown). Another class of genes are classical, basal, or mesenchymal markers that are expressed more strongly in autochthonous tumors when compared with the transplant models.

Immunocompromised NSG mice provide a convenient and rapid platform for orthotopic transplantation, however the resulting tumors evolve in a tumor microenvironment lacking an adaptive immune system. To investigate the degree to which the adaptive immune system might impact intratumor heterogeneity, we performed orthotopic transplants of tumor fragments, primary 2D, and primary 3D organoid cell lines into syngeneic immunocompetent (F1) hosts (Supplementary Fig. S6A). Similar to our results in NSG mice, tumors derived from orthotopic transplantation of tumor fragments, 2D and 3D murine cell lines into F1 mice exhibited the full spectrum of classical, basal, and mesenchymal cell states. We observed similar phenotypic volume and compositional entropy in F1 compared with NSG mice, with fragment-derived tumors again skewing towards a more mesenchymal cell state (Supplementary Fig. S6A–S6C). We were able to identify differentially expressed genes comparing the NSG with F1 models, most notably the enrichment of inflammatory pathways, in particular IFNα and IFNγ, in the F1 tumors (Supplementary Fig. S6D–S6F). However, the overall heterogeneity of the cell states continued to be predominantly defined by the continuum of classical, basal, and mesenchymal cell states.

We found that the cancer cells in the autochthonous KPC GEMM and in all KPC transplant models organized along an epithelial–mesenchymal continuum (Fig. 1 and Fig. 2). At the center of this principal axis of heterogeneity, we observed the basal cell state that was conserved across all in vivo model systems. In addition to a specific basal gene expression program, this state coexpresses low levels of both classical and mesenchymal signature genes and shows a high degree of gene expression entropy (Fig. 3A; Supplementary Fig. S2C and S2D). To develop insights into the PDAC cell state continuum at the chromatin level, we performed scATAC-seq on autochthonous KPCT tumors. Consistent with the scRNA-seq data, we observed a continuum of chromatin accessibility states spanning classical, basal, and mesenchymal phenotypes (Fig. 3B). This progression was associated with an increase in global chromatin accessibility, with classical cells having the least open loci and cells with the basal and mesenchymal phenotypes exhibiting substantially higher chromatin openness (Fig. 3C). As expected, classical cells exhibited accessibility at epithelial and classical marker genes, such as Epcam and Lgals4, respectively, whereas mesenchymal genes such as Vim were not accessible (Fig. 3D; Supplementary Fig. S7A). Conversely, cells in the mesenchymal state had decreased accessibility at epithelial genes. Interestingly, cells in the basal phenotype uniquely showed co-accessibility at epithelial/classical and mesenchymal genes (Fig. 3D; Supplementary Fig. S7A). This bivalent chromatin state and the co-expression of classical and mesenchymal genes suggested that the basal phenotype may be functionally highly plastic, poised to engage either the classical or mesenchymal differentiation programs. Further supporting high plasticity of the basal state, it demonstrated concordance with a high-plasticity cell state that we recently discovered in lung adenocarcinoma (Supplementary Fig. S7B; ref. 46).

Figure 3.

Lineage tracing of a subset of basal state cells demonstrates functional phenotypic plasticity. A, Pseudotime analysis of cell states demonstrating mixed classical and mesenchymal expression profiles within the basal cell phenotype. Colored lines represent the average model-fitted expression of all pseudotime-dependent genes belonging to either the classical (purple), basal (green), and mesenchymal w(orange) gene expression clusters established in Supplementary Fig. S2C. Shading represents the area between the 5th and 95th percentile of genes. The panel on the top shows the fraction of classical (purple), basal (green), and mesenchymal (orange) along the same pseudotime axis divided into ten bins. B, scATAC-seq of autochthonous KPCT tumors at 7 to 8 weeks of age (n = 4). C, Total chromatin accessibility, calculated as number of fragments per cell. D, Accessibility of representative loci in each cell state. Note coaccessibility at classical (Lgals4) and mesenchymal (Vim) marker genes in the basal state. E,Lgr5 overlaps with basal cell phenotype in the autochthonous (left) and fragment-based (right) tumor models. UMAPs displaying Lgr5 expression. Scale represents size factor normalized log2-transformed UMI counts, with color scale from 1st to 99th percentile. F, Schematic of KPF; Rosa26mTmG/+; Lgr5CreER/CreER lineage-tracing system. G, Lineage tracing of fragment-derived KPF; Rosa26mTmG/+; Lgr5CreER/CreER tumors. Day 3 (top) and day 10 (bottom) post-tamoxifen labeling. Left, average expression of Lgr5 in tdTomato+ versus GFP+ cells at day 3 and day 10. Note the initial enrichment of Lgr5 gene expression at day 3, which is no longer present at day 10, reflecting transition of the basal cells to other cell states. Right, UMAP embedding displaying the position of GFP+ traced cells either as single cells. Note how at day 3 the cells occupy a relatively restricted phenotypic space and diffuse over time. H, Phenotypic volume of fragment-derived GFP+ traced cells over time. I, Immunofluorescence of the classical marker galectin-4 and mesenchymal marker vimentin in the Lgr5-traced cells. At initial labeling, Lgr5-traced cells express GFP but do not express galectin-4 or vimentin. Second inset shows traced cells at day 10, with white arrows indicating co-expression of GFP with galectin-4 (top) or vimentin (bottom), consistent with the acquisition of classical and mesenchymal phenotypes, respectively. J, Graphical summary of lineage-tracing system demonstrating plasticity of the basal cell population.

Figure 3.

Lineage tracing of a subset of basal state cells demonstrates functional phenotypic plasticity. A, Pseudotime analysis of cell states demonstrating mixed classical and mesenchymal expression profiles within the basal cell phenotype. Colored lines represent the average model-fitted expression of all pseudotime-dependent genes belonging to either the classical (purple), basal (green), and mesenchymal w(orange) gene expression clusters established in Supplementary Fig. S2C. Shading represents the area between the 5th and 95th percentile of genes. The panel on the top shows the fraction of classical (purple), basal (green), and mesenchymal (orange) along the same pseudotime axis divided into ten bins. B, scATAC-seq of autochthonous KPCT tumors at 7 to 8 weeks of age (n = 4). C, Total chromatin accessibility, calculated as number of fragments per cell. D, Accessibility of representative loci in each cell state. Note coaccessibility at classical (Lgals4) and mesenchymal (Vim) marker genes in the basal state. E,Lgr5 overlaps with basal cell phenotype in the autochthonous (left) and fragment-based (right) tumor models. UMAPs displaying Lgr5 expression. Scale represents size factor normalized log2-transformed UMI counts, with color scale from 1st to 99th percentile. F, Schematic of KPF; Rosa26mTmG/+; Lgr5CreER/CreER lineage-tracing system. G, Lineage tracing of fragment-derived KPF; Rosa26mTmG/+; Lgr5CreER/CreER tumors. Day 3 (top) and day 10 (bottom) post-tamoxifen labeling. Left, average expression of Lgr5 in tdTomato+ versus GFP+ cells at day 3 and day 10. Note the initial enrichment of Lgr5 gene expression at day 3, which is no longer present at day 10, reflecting transition of the basal cells to other cell states. Right, UMAP embedding displaying the position of GFP+ traced cells either as single cells. Note how at day 3 the cells occupy a relatively restricted phenotypic space and diffuse over time. H, Phenotypic volume of fragment-derived GFP+ traced cells over time. I, Immunofluorescence of the classical marker galectin-4 and mesenchymal marker vimentin in the Lgr5-traced cells. At initial labeling, Lgr5-traced cells express GFP but do not express galectin-4 or vimentin. Second inset shows traced cells at day 10, with white arrows indicating co-expression of GFP with galectin-4 (top) or vimentin (bottom), consistent with the acquisition of classical and mesenchymal phenotypes, respectively. J, Graphical summary of lineage-tracing system demonstrating plasticity of the basal cell population.

Close modal

To functionally investigate the plasticity of the basal cells, we sought to identify a marker gene enabling their lineage tracing in situ in established PDAC tumors. Interestingly, we observed strong enrichment of the adult epithelial stem cell marker Leucine-rich repeat-containing G protein-coupled receptor-5 (Lgr5) in the basal cells (Fig. 3E; Supplementary Fig. S7C). In mice, we observed co-localization of an Lgr5-GFP reporter and the basal marker CK17 in PDAC tumors in KPCT; Lgr5GFP-CreER/+ mice (Supplementary Fig. S8A). Conversely, the classical epithelial marker galectin-4 or the mesenchymal marker vimentin were not coexpressed with Lgr5 (Supplementary Fig. S8A). Combined analysis of normal pancreas from Genotype-Tissue Expression Program and tumors from The Cancer Genome Atlas revealed that LGR5 is robustly upregulated in human PDAC when compared with normal pancreas (Supplementary Fig. S8B). In situ hybridization of primary human PDAC tissues uncovered spatially restricted LGR5 expression in a subset of tumor cells (Supplementary Fig. S8C). We observed similar restricted expression of LGR5 in published human PDAC scRNA-seq data (Supplementary Fig. S8D; ref. 47). Interestingly, in contrast to our mouse data where the Lgr5 cells had a strong basal phenotype, in human data the LGR5+ cells overlapped more with the classical cell state (Supplementary Fig. S8D, right).

Having established that Lgr5 is a marker of the murine PDAC basal cell state, we crossed Lgr5CreERT2 and Rosa26mTmG alleles into the KrasFSF-G12D/+; Tp53 frt/frt; Pdx1-Flp; (KPF) PDAC model. In these KPF; Lgr5CreERT2/CreERT2; Rosa26mTmG/+ mice, PDAC tumorigenesis is initiated by Flp recombinase and the Lgr5+ lineage can be inducibly labeled with a heritable genetic switch from membrane-associated (m)tdTomato to mGFP allele using a single pulse of tamoxifen in established tumors (Fig. 3F). To measure plasticity of the basal cells, we performed scRNA-seq on mGFP+ descendants of the Lgr5+ cells at 3 days (baseline after tamoxifen washout) and at 10 days after the tamoxifen pulse in the orthotopic tumor fragment model (Fig. 3F and G). TdTomato-expressing descendants of Lgr5- cells were also subjected to scRNA-seq as an internal control. Notably, 3 days after the initial tamoxifen pulse PDAC cells expressing mGFP mRNA were significantly enriched for Lgr5 expression and had a strong overlap with the basal phenotype (Fig. 3G, top row). In contrast, 10 days after tamoxifen pulse mGFP+ cells displayed decreased Lgr5 expression and had diffused away from the basal state, showing increased overlap with the classical and mesenchymal cell states (Fig. 3G, bottom row). This corresponded with an increase in phenotypic diversity and compositional heterogeneity of mGFP+ cells (Fig. 3H; Supplementary Fig. S9A). We confirmed these results in situ by immunofluorescence: GFP+ cells did not express classical or mesenchymal markers galectin-4 (Fig. 3I, top) or vimentin (Fig. 3I, bottom), respectively, at baseline (3 days post-tamoxifen). However, consistent with the scRNA-seq data, we detected coexpression of mGFP and galectin-4 or vimentin at 10 days post labeling in subsets of PDAC cells (Fig. 3I), indicating that the LGR5+ cells can differentiate into either classical or mesenchymal states, respectively (Fig. 3J). We observed similar results when performing scRNA-seq and lineage tracing in autochthonous PDAC tumors (Supplementary Fig. S9B–S9D). Taken together, these results demonstrate that the basal PDAC cell state is endowed with high plasticity.

We found that the classical, basal, and mesenchymal cancer cell states described here overlap with previously reported transcriptional subtypes of human PDAC (9, 10, 12, 19, 20). Thus, cell state heterogeneity in the genetically engineered KPC mouse model is concordant with the spectrum of cell states observed in human PDAC. The mouse model largely mirrored the human continuum of cell states, centering on a basal cell state that was characterized by coexpression of and accessible chromatin at both epithelial and mesenchymal genes. Interestingly, previous studies have demonstrated that the PDAC cells with the most metastatic capacity reside in a hybrid (or partial) EMT state (16, 17), suggesting that retention of both mesenchymal and epithelial features may confer a high capacity for adaptation, such as is required of cancer cells for colonizing a distant tissue site. Our results suggest that the basal state harboring these hybrid features is endowed with increased adaptability and plasticity also at the primary site. Similar high-plasticity cell states have recently been reported in cutaneous and lung squamous cell carcinomas (48), lung adenocarcinoma (46), small cell lung cancer (43), and prostate cancer (49). Interestingly, the basal cells shared features of a high-plasticity cell state that we recently uncovered in lung adenocarcinoma (46). These findings raise the tantalizing possibility that the molecular mechanisms controlling plastic or hybrid cell states may be shared across cancer types; such mechanisms are an important area of further investigation.

Our studies leverage a murine KP(fl/fl)C model to investigate tumor heterogeneity and plasticity, which has a less variable tumor time course and more rapid progression than other commonly used KPC models employing either Trp53flox/+ or Trp53R172H/+ allelic configurations (24). TP53 deletion is uncommon in human PDAC, and mouse models of PDAC have demonstrated that Trp53 mutations have novel gain-of-function properties that promote invasion and metastasis (50, 51). Furthermore, Trp53 mutations have been implicated in remodeling the tumor microenvironment (52, 53). The degree to which TP53 mutations influence plasticity and tumor heterogeneity warrants further inquiry.

A predominance of the basal cell state signature in PDAC bulk gene expression analysis is associated with poor prognosis and aggressive clinical features in patients, such as advanced stage and resistance to chemotherapy (12, 54). Our results demonstrate the basal cell state is a transition state between the classical and mesenchymal phenotypes. Expanding on this observation, elimination of the basal cell state would be expected to blunt the ability of the classical cells to transition into mesenchymal cells and vice versa. Based on previous work, targeting the basal cells could also suppress metastasis (16, 17). Furthermore, plasticity contributes significantly to the failure of chemo-, targeted-, and immunotherapies across cancers, allowing cancer cells to acquire states that are adapted to therapy (4, 17, 55). Taken together, targeting the basal cell state either via cell ablation therapies or by disrupting its molecular drivers may have therapeutic potential for suppressing plasticity in PDAC.

Although still incomplete, a picture of transcriptional regulators (56) and microenvironmental inputs (20) driving the basal state has begun to emerge. Interestingly, WNT/β-catenin signaling is one of the defining attributes of the human basal PDAC cells (20). In agreement, we found that the murine basal cells express Lgr5, an established WNT target gene and epithelial stem cell marker (33). Interestingly, LGR5 is also expressed in human PDAC, although its expression is enriched in the classical cell state. Despite this difference, mouse Lgr5 provided a useful cell state marker for the purposes of our study. By combining fluorescence switch-based lineage tracing of the Lgr5+ murine PDAC cells and scRNA-seq we demonstrated that the basal cells are bi-potent in their differentiation potential in situ in established tumors. These experiments highlight the power and sensitivity of fluorescence-based tracing combined with scRNA-seq for investigating cell state heterogeneity. LGR5 is not expressed in the normal pancreas but is induced in plastic epithelial progenitors upon pancreas injury (57). We found that the transition from the classical state to the basal state in the PDAC cell state continuum is accompanied by a drastic increase in overall chromatin accessibility. Taken together, the basal state in PDAC may be associated with regenerative programs, which are known to involve the acquisition of plasticity that is defined by increased chromatin accessibility (58). Regenerative pathways and molecular mechanisms controlling the expansion of chromatin accessibility in the basal cell state may present additional entry points for therapeutic interventions in PDAC.

A major challenge in developing novel therapeutic strategies for PDAC has been the limited success of in vitro and in vivo preclinical models to accurately predict responses to investigational therapeutics in clinical trials. Recent work has demonstrated that PDAC cells are more sensitive to chemotherapy in 3D organoid culture, which enriches for the classical state, than in the more basal-like 2D culture conditions (13, 20). Forcing basal gene expression in 3D organoids by stimulation with TGFβ promoted chemoresistance (20). Thus, the transcriptional state and the diversity of such states are critical determinants of treatment response in PDAC and other cancers (4). We observed that transplanted tumors that included stroma from the primary tumor, such as unsorted tumor cells or primary tumor fragments, skewed towards a more mesenchymal phenotype, whereas transplants of purified tumor cells had a more balanced intratumoral diversity. The evolution of the tumor stroma over time, the differences between stroma in the primary autochthonous and transplanted settings, and how these differences influence cell state heterogeneity, remains an important open question.

Our results indicate that the 2D culture and 3D organoid culture conditions produce homogenous cell states where basal and classical gene expression signatures are mixed. However, upon orthotopic transplantation cells from either condition reconstitute the classical, basal, and mesenchymal phenotypic extremes observed in the autochthonous tumors. This finding suggests that cell culture intermediates can be safely used to experimentally manipulate or expand cells and still generate transplant tumors that are biologically highly similar to the primary tumors. Although organoids incorporate many features of PDAC tumors that are lost in 2D culture, such as 3D cell–cell interactions (13, 28), our results indicate that the full spectrum of intratumoral heterogeneity can only be achieved in in vivo model systems. This has broad biological and clinical significance across the constellation of human cancers. Future work will elucidate whether human patient–derived PDAC models and models of other cancers similarly reconstitute cell state heterogeneity upon transplantation.

K.L. Pitter reports grants from NCI and NIH during the study. C.A. Iacobuzio-Donahue reports other support from Bristol-Myers Squibb outside the submitted work. D. Pe'er reports personal fees from insitro outside the submitted work. T. Tammela reports grants from NCI, NIH, Starr Cancer Consortium, Josie Robertson Foundation, American Cancer Society, Rita Allen Fellowship; and grants from V Foundation during the conduct of the study. No disclosures were reported by the other authors.

K.L. Pitter: Conceptualization, data curation, formal analysis, supervision, funding acquisition, validation, investigation, visualization, methodology, writing–original draft, project administration, writing–review and editing. O. Grbovic-Huezo: Conceptualization, data curation, formal analysis, supervision, funding acquisition, validation, investigation, visualization, methodology, writing–original draft, project administration, writing–review and editing. S. Joost: Conceptualization, data curation, software, formal analysis, funding acquisition, validation, investigation, visualization, methodology, writing–original draft, writing–review and editing. A. Singhal: Data curation, formal analysis, supervision, investigation, visualization. M. Blum: Investigation. K. Wu: Data curation, formal analysis, investigation, visualization. M. Holm: Data curation, validation, investigation. A. Ferrena: Data curation, software, formal analysis, visualization. A. Bhutkar: Data curation, software, formal analysis, visualization. A. Hudson: Data curation, investigation, visualization. N. Lecomte: Methodology. E. de Stanchina: Supervision, investigation. R. Chaligne: Resources, supervision, methodology. C.A. Iacobuzio-Donahue: Resources, supervision, funding acquisition, methodology. D. Pe'er: Resources, data curation, supervision, funding acquisition, methodology, writing–original draft, project administration. T. Tammela: Conceptualization, resources, data curation, formal analysis, supervision, funding acquisition, methodology, writing–original draft, project administration, writing–review and editing.

The authors thank the members of the Tammela Laboratory and V. Balachandran (MSKCC) for helpful discussions and comments on the manuscript; R. Gardner for FACS support; K. Manova for histology support; J. Chan, G. Hartmann, S. Torborg, and X. Zhuang for help with experiments; and D. Alonso-Curbelo and J.P. Morris IV for helpful discussions and support with tumor dissociation methodology. This work was supported by Cycle for Survival, NIH/NCI grant R37-CA244911, the Starr Cancer Consortium, and, in part, by the NIH/NCI Cancer Center Support Grant P30-CA08748 (MSKCC). T. Tammela is supported by Josie Robertson, American Cancer Society, Rita Allen, and V Foundation Scholarships. K.L. Pitter is supported by NIH Director's Early Independence Award DP5-OD031864. S. Joost is supported by the Hope Funds for Cancer Research. A. Singhal is supported by the T32 Investigational Cancer Therapeutics Training Program Grant (NIH MSK ICTTP T32-CA009207) to the Sloan Kettering Institute. R. Chaligne is supported by The Alan and Sandra Gerry Foundation. The authors acknowledge the use of the Integrated Genomics Operation Core, funded by CCSG P30-CA08748, Cycle for Survival, the Marie-Josée and Henry R. Kravis Center for Molecular Oncology at MSKCC, and the Flow Cytometry and Histology Core Facilities at Sloan Kettering Institute.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.

Note: Supplementary data for this article are available at Cancer Research Online (http://cancerres.aacrjournals.org/).

1.
National Cancer Institute: Cancer Statistics (2021)
. https://www.cancer.gov/about-cancer/understanding/statistics.
2.
Gupta
PB
,
Pastushenko
I
,
Skibinski
A
,
Blanpain
C
,
Kuperwasser
C
.
Phenotypic plasticity: driver of cancer initiation, progression, and therapy resistance
.
Cell Stem Cell
2019
;
24
:
65
78
.
3.
Torborg
SR
,
Li
Z
,
Chan
JE
,
Tammela
T
.
Cellular and molecular mechanisms of plasticity in cancer
.
Trends Cancer
2022
;
S2405-8033:00094-2
.
4.
Marine
JC
,
Dawson
SJ
,
Dawson
MA
.
Nongenetic mechanisms of therapeutic resistance in cancer
.
Nat Rev Cancer
2020
;
20
:
743
56
.
5.
Cancer Genome Atlas Research Network
.
Electronic address aadhe, cancer genome atlas research N: integrated genomic characterization of pancreatic ductal adenocarcinoma
.
Cancer Cell
2017
;
32
:
185
203
.
6.
Flowers
BM
,
Xu
H
,
Mulligan
AS
,
Hanson
KJ
,
Seoane
JA
,
Vogel
H
, et al
.
Cell of origin influences pancreatic cancer subtype
.
Cancer Discov
2021
;
11
:
660
77
.
7.
Espinet
E
,
Gu
Z
,
Imbusch
CD
,
Giese
NA
,
Buscher
M
,
Safavi
M
, et al
.
Aggressive PDACs show hypomethylation of repetitive elements and the execution of an intrinsic IFN program linked to a ductal cell of origin
.
Cancer Discov
2021
;
11
:
638
59
.
8.
Bardeesy
N
,
Aguirre
AJ
,
Chu
GC
,
Cheng
KH
,
Lopez
LV
,
Hezel
AF
, et al
.
Both p16(Ink4a) and the p19(Arf)-p53 pathway constrain progression of pancreatic adenocarcinoma in the mouse
.
Proc Natl Acad Sci U S A
2006
;
103
:
5947
52
.
9.
Bailey
P
,
Chang
DK
,
Nones
K
,
Johns
AL
,
Patch
AM
,
Gingras
MC
, et al
.
Genomic analyses identify molecular subtypes of pancreatic cancer
.
Nature
2016
;
531
:
47
52
.
10.
Collisson
EA
,
Sadanandam
A
,
Olson
P
,
Gibb
WJ
,
Truitt
M
,
Gu
S
, et al
.
Subtypes of pancreatic ductal adenocarcinoma and their differing responses to therapy
.
Nat Med
2011
;
17
:
500
3
.
11.
Martens
S
,
Lefesvre
P
,
Nicolle
R
,
Biankin
AV
,
Puleo
F
,
Van Laethem
JL
, et al
.
Different shades of pancreatic ductal adenocarcinoma, different paths towards precision therapeutic applications
.
Ann Oncol
2019
;
30
:
1428
36
.
12.
Moffitt
RA
,
Marayati
R
,
Flate
EL
,
Volmar
KE
,
Loeza
SG
,
Hoadley
KA
, et al
.
Virtual microdissection identifies distinct tumor- and stroma-specific subtypes of pancreatic ductal adenocarcinoma
.
Nat Genet
2015
;
47
:
1168
78
.
13.
Tiriac
H
,
Belleau
P
,
Engle
DD
,
Plenker
D
,
Deschenes
A
,
Somerville
TDD
, et al
.
Organoid profiling identifies common responders to chemotherapy in pancreatic cancer
.
Cancer Discov
2018
;
8
:
1112
29
.
14.
Driehuis
E
,
van Hoeck
A
,
Moore
K
,
Kolders
S
,
Francies
HE
,
Gulersonmez
MC
, et al
.
Pancreatic cancer organoids recapitulate disease and allow personalized drug screening
.
Proc Natl Acad Sci U S A
2019
;
116
:
26580
90
.
15.
Rhim
AD
,
Mirek
ET
,
Aiello
NM
,
Maitra
A
,
Bailey
JM
,
McAllister
F
, et al
.
EMT and dissemination precede pancreatic tumor formation
.
Cell
2012
;
148
:
349
61
.
16.
Simeonov
KP
,
Byrns
CN
,
Clark
ML
,
Norgard
RJ
,
Martin
B
,
Stanger
BZ
, et al
.
Single-cell lineage tracing of metastatic cancer reveals selection of hybrid EMT states
.
Cancer Cell
2021
;
39
:
1150
62
.
17.
Aiello
NM
,
Maddipati
R
,
Norgard
RJ
,
Balli
D
,
Li
J
,
Yuan
S
, et al
.
EMT subtype influences epithelial plasticity and mode of cell migration
.
Dev Cell
2018
;
45
:
681
95
.
18.
Juiz
N
,
Elkaoutari
A
,
Bigonnet
M
,
Gayet
O
,
Roques
J
,
Nicolle
R
, et al
.
Basal-like and classical cells coexist in pancreatic cancer revealed by single-cell analysis on biopsy-derived pancreatic cancer organoids from the classical subtype
.
FASEB J
2020
;
34
:
12214
28
.
19.
Krieger
TG
,
Le Blanc
S
,
Jabs
J
,
Ten
FW
,
Ishaque
N
,
Jechow
K
, et al
.
Single-cell analysis of patient-derived PDAC organoids reveals cell state heterogeneity and a conserved developmental hierarchy
.
Nat Commun
2021
;
12
:
5826
.
20.
Raghavan
S
,
Winter
PS
,
Navia
AW
,
Williams
HL
,
DenAdel
A
,
Lowder
KE
, et al
.
Microenvironment drives cell state, plasticity, and drug response in pancreatic cancer
.
Cell
2021
;
184
:
6119
37
.
21.
Lee
JJ
,
Bernard
V
,
Semaan
A
,
Monberg
ME
,
Huang
J
,
Stephens
BM
, et al
.
Elucidation of tumor-stromal heterogeneity and the ligand-receptor interactome by single-cell transcriptomics in real-world pancreatic cancer biopsies
.
Clin Cancer Res
2021
;
27
:
5912
21
.
22.
Tammela
T
,
Sage
J
.
Investigating tumor heterogeneity in mouse models
.
Annu Rev Canc Biol
2020
;
4
:
99
119
.
23.
Hingorani
SR
,
Wang
L
,
Multani
AS
,
Combs
C
,
Deramaudt
TB
,
Hruban
RH
, et al
.
Trp53R172H and KrasG12D cooperate to promote chromosomal instability and widely metastatic pancreatic ductal adenocarcinoma in mice
.
Cancer Cell
2005
;
7
:
469
83
.
24.
Perez-Mancera
PA
,
Guerra
C
,
Barbacid
M
,
Tuveson
DA
.
What we have learned about pancreatic cancer from mouse models
.
Gastroenterology
2012
;
142
:
1079
92
.
25.
Schonhuber
N
,
Seidler
B
,
Schuck
K
,
Veltkamp
C
,
Schachtler
C
,
Zukowska
M
, et al
.
A next-generation dual-recombinase system for time- and host-specific targeting of pancreatic cancer
.
Nat Med
2014
;
20
:
1340
7
.
26.
Lee
JW
,
Komar
CA
,
Bengsch
F
,
Graham
K
,
Beatty
GL
.
Genetically engineered mouse models of pancreatic cancer: the kpc model (LSL-Kras(G12D/+);LSL-Trp53(R172H/+);Pdx-1-Cre), its variants, and their application in immuno-oncology drug discovery
.
Curr Protoc Pharmacol
2016
;
73
:
14.39.1-14.39.20. doi: 10.1002/cpph.2
.
27.
Hosein
AN
,
Huang
H
,
Wang
Z
,
Parmar
K
,
Du
W
,
Huang
J
, et al
.
Cellular heterogeneity during mouse pancreatic ductal adenocarcinoma progression at single-cell resolution
.
JCI Insight
2019
;
5
:
e129212
.
28.
Boj
SF
,
Hwang
CI
,
Baker
LA
,
Chio
II
,
Engle
DD
,
Corbo
V
, et al
.
Organoid models of human and mouse ductal pancreatic cancer
.
Cell
2015
;
160
:
324
38
.
29.
Erstad
DJ
,
Sojoodi
M
,
Taylor
MS
,
Ghoshal
S
,
Razavi
AA
,
Graham-O'Regan
KA
, et al
.
Orthotopic and heterotopic murine models of pancreatic cancer and their different responses to FOLFIRINOX chemotherapy
.
Dis Model Mech
2018
;
11
:
dmm034793
.
30.
Jackson
EL
,
Willis
N
,
Mercer
K
,
Bronson
RT
,
Crowley
D
,
Montoya
R
, et al
.
Analysis of lung tumor initiation and progression using conditional expression of oncogenic K-ras
.
Genes Dev
2001
;
15
:
3243
8
.
31.
Marino
S
,
Vooijs
M
,
van Der Gulden
H
,
Jonkers
J
,
Berns
A
.
Induction of medulloblastomas in p53-null mutant mice by somatic inactivation of Rb in the external granular layer cells of the cerebellum
.
Genes Dev
2000
;
14
:
994
1004
.
32.
Madisen
L
,
Zwingman
TA
,
Sunkin
SM
,
Oh
SW
,
Zariwala
HA
,
Gu
H
, et al
.
A robust and high-throughput Cre reporting and characterization system for the whole mouse brain
.
Nat Neurosci
2010
;
13
:
133
40
.
33.
Barker
N
,
van Es
JH
,
Kuipers
J
,
Kujala
P
,
van den Born
M
,
Cozijnsen
M
, et al
.
Identification of stem cells in small intestine and colon by marker gene Lgr5
.
Nature
2007
;
449
:
1003
7
.
34.
Young
NP
,
Crowley
D
,
Jacks
T
.
Uncoupling cancer mutations reveals critical timing of p53 loss in sarcomagenesis
.
Cancer Res
2011
;
71
:
4040
7
.
35.
Lee
CL
,
Moding
EJ
,
Huang
X
,
Li
Y
,
Woodlief
LZ
,
Rodrigues
RC
, et al
.
Generation of primary tumors with Flp recombinase in FRT-flanked p53 mice
.
Dis Model Mech
2012
;
5
:
397
402
.
36.
Wu
J
,
Liu
X
,
Nayak
SG
,
Pitarresi
JR
,
Cuitino
MC
,
Yu
L
, et al
.
Generation of a pancreatic cancer model using a Pdx1-Flp recombinase knock-in allele
.
PLoS One
2017
;
12
:
e0184984
.
37.
Muzumdar
MD
,
Tasic
B
,
Miyamichi
K
,
Li
L
,
Luo
L
.
A global double-fluorescent Cre reporter mouse
.
Genesis
2007
;
45
:
593
605
.
38.
Huch
M
,
Dorrell
C
,
Boj
SF
,
van Es
JH
,
Li
VS
,
van de Wetering
M
, et al
.
In vitro expansion of single Lgr5+ liver stem cells induced by Wnt-driven regeneration
.
Nature
2013
;
494
:
247
50
.
39.
Russell
J
,
Grkovski
M
,
O'Donoghue
IJ
,
Kalidindi
TM
,
Pillarsetty
N
,
Burnazi
EM
, et al
.
Predicting gemcitabine delivery by (18)F-FAC PET in murine models of pancreatic cancer
.
J Nucl Med
2021
;
62
:
195
200
.
40.
Pylayeva-Gupta
Y
,
Das
S
,
Handler
JS
,
Hajdu
CH
,
Coffre
M
,
Koralov
SB
, et al
.
IL35-producing B cells promote the development of pancreatic neoplasia
.
Cancer Discov
2016
;
6
:
247
55
.
41.
Stoeckius
M
,
Zheng
S
,
Houck-Loomis
B
,
Hao
S
,
Yeung
BZ
,
Mauck
WM
3rd
, et al
.
Cell Hashing with barcoded antibodies enables multiplexing and doublet detection for single cell genomics
.
Genome Biol
2018
;
19
:
224
.
42.
Tammela
T
,
Sanchez-Rivera
FJ
,
Cetinbas
NM
,
Wu
K
,
Joshi
NS
,
Helenius
K
, et al
.
A Wnt-producing niche drives proliferative potential and progression in lung adenocarcinoma
.
Nature
2017
;
545
:
355
9
.
43.
Chan
JM
,
Quintanal-Villalonga
A
,
Gao
VR
,
Xie
Y
,
Allaj
V
,
Chaudhary
O
, et al
.
Signatures of plasticity, metastasis, and immunosuppression in an atlas of human small cell lung cancer
.
Cancer Cell
2021
;
39
:
1479
96
.
44.
Levine
JH
,
Simonds
EF
,
Bendall
SC
,
Davis
KL
,
Amir el
AD
,
Tadmor
MD
, et al
.
Data-driven phenotypic dissection of AML reveals progenitor-like cells that correlate with prognosis
.
Cell
2015
;
162
:
184
97
.
45.
Azizi
E
,
Carr
AJ
,
Plitas
G
,
Cornish
AE
,
Konopacki
C
,
Prabhakaran
S
, et al
.
Single-cell map of diverse immune phenotypes in the breast tumor microenvironment
.
Cell
2018
;
174
:
1293
308
.
46.
Marjanovic
ND
,
Hofree
M
,
Chan
JE
,
Canner
D
,
Wu
K
,
Trakala
M
, et al
.
Emergence of a high-plasticity cell state during lung cancer evolution
.
Cancer Cell
2020
;
38
:
229
46
.
47.
Peng
J
,
Sun
BF
,
Chen
CY
,
Zhou
JY
,
Chen
YS
,
Chen
H
, et al
.
Single-cell RNA-seq highlights intra-tumoral heterogeneity and malignant progression in pancreatic ductal adenocarcinoma
.
Cell Res
2019
;
29
:
725
38
.
48.
Pastushenko
I
,
Mauri
F
,
Song
Y
,
de Cock
F
,
Meeusen
B
,
Swedlund
B
, et al
.
Fat1 deletion promotes hybrid EMT state, tumor stemness, and metastasis
.
Nature
2021
;
589
:
448
55
.
49.
Quintanal-Villalonga
A
,
Chan
JM
,
Yu
HA
,
Pe'er
D
,
Sawyers
CL
,
Triparna
S
, et al
.
Lineage plasticity in cancer: a shared pathway of therapeutic resistance
.
Nat Rev Clin Oncol
2020
;
17
:
360
71
.
50.
Morton
JP
,
Timpson
P
,
Karim
SA
,
Ridgway
RA
,
Athineos
D
,
Doyle
B
, et al
.
Mutant p53 drives metastasis and overcomes growth arrest/senescence in pancreatic cancer
.
Proc Natl Acad Sci U S A
2010
;
107
:
246
51
.
51.
Weissmueller
S
,
Manchado
E
,
Saborowski
M
,
Morris
JPt
,
Wagenblast
E
,
Davis
CA
, et al
.
Mutant p53 drives pancreatic cancer metastasis through cell-autonomous PDGF receptor beta signaling
.
Cell
2014
;
157
:
382
94
.
52.
Blagih
J
,
Zani
F
,
Chakravarty
P
,
Hennequart
M
,
Pilley
S
,
Hobor
S
, et al
.
Cancer-specific loss of p53 leads to a modulation of myeloid and T-cell responses
.
Cell Rep
2020
;
30
:
481
96
.
53.
Maddalena
M
,
Mallel
G
,
Nataraj
NB
,
Shreberk-Shaked
M
,
Hassin
O
,
Mukherjee
S
, et al
.
TP53 missense mutations in PDAC are associated with enhanced fibrosis and an immunosuppressive microenvironment
.
Proc Natl Acad Sci U S A
2021
;
118
:
e2025631118
.
54.
Aung
KL
,
Fischer
SE
,
Denroche
RE
,
Jang
G-H
,
Dodd
A
,
Creighton
S
, et al
.
Genomics-driven precision medicine for advanced pancreatic cancer: early results from the COMPASS trial
.
Clin Cancer Res
2018
;
24
:
1344
54
.
55.
Porter
RL
,
Magnus
NKC
,
Thapar
V
,
Morris
R
,
Szabolcs
A
,
Neyaz
A
, et al
.
Epithelial to mesenchymal plasticity and differential response to therapies in pancreatic ductal adenocarcinoma
.
Proc Natl Acad Sci U S A
2019
;
116
:
26835
45
.
56.
Camolotto
SA
,
Belova
VK
,
Torre-Healy
L
,
Vahrenkamp
JM
,
Berrett
KC
,
Conway
H
, et al
.
Reciprocal regulation of pancreatic ductal adenocarcinoma growth and molecular subtype by HNF4alpha and SIX1/4
.
Gut
2021
;
70
:
900
14
.
57.
Huch
M
,
Bonfanti
P
,
Boj
SF
,
Sato
T
,
Loomans
CJ
,
van de Wetering
M
, et al
.
Unlimited in vitro expansion of adult bi-potent pancreas progenitors through the Lgr5/R-spondin axis
.
EMBO J
2013
;
32
:
2708
21
.
58.
Gola
A
,
Fuchs
E
.
Environmental control of lineage plasticity and stem cell memory
.
Curr Opin Cell Biol
2021
;
69
:
88
95
.