Persister cell phenotypes contribute to poor patient outcomes after neoadjuvant chemotherapy in PDAC

Neoadjuvant chemotherapy can improve the survival of individuals with borderline and unresectable pancreatic ductal adenocarcinoma; however, heterogeneous responses to chemotherapy remain a significant clinical challenge. Here, we performed RNA sequencing (n = 97) and multiplexed immunofluorescence (n = 122) on chemo-naive and postchemotherapy (post-CTX) resected patient samples (chemoradiotherapy excluded) to define the impact of neoadjuvant chemotherapy. Transcriptome analysis combined with high-resolution mapping of whole-tissue sections identified GATA6 (classical), KRT17 (basal-like) and cytochrome P450 3A (CYP3A) coexpressing cells that were preferentially enriched in post-CTX resected samples. The persistence of GATA6hi and KRT17hi cells post-CTX was significantly associated with poor survival after mFOLFIRINOX (mFFX), but not gemcitabine (GEM), treatment. Analysis of organoid models derived from chemo-naive and post-CTX samples demonstrated that CYP3A expression is a predictor of chemotherapy response and that CYP3A-expressing drug detoxification pathways can metabolize the prodrug irinotecan, a constituent of mFFX. These findings identify CYP3A-expressing drug-tolerant cell phenotypes in residual disease that may ultimately inform adjuvant treatment selection.

Article https://doi.org/10.1038/s43018-023-00628-6 (FFPE) tissues from 122 individuals and on normal pancreas tissues from 9 organ donors. RNA-seq and multiplex IF were performed on identical samples from 48 individuals. Laser-capture microdissection (LCM) was performed on 32 chemo-naive samples, which were analyzed by RNA-seq. Assessment of clinical stage 19 demonstrated that chemo-naive and post-CTX individuals had similar numbers of late-stage tumors (stage III and IV), although a higher number of stage IIB tumors were present in the chemo-naive group (Extended Data Fig. 1a). Chemo-naive samples selected for RNA-seq and IF analysis received GEM monoadjuvant chemotherapy almost exclusively, with a small number of individuals receiving either mFFX (RNA-seq: 5 of 60; IF: 11 of 76) or GEM nab-paclitaxel (RNA-seq: 3 of 60; IF: 5 of 76). Post-CTX samples received either mFFX or GEM nab-paclitaxel. Most individuals who received adjuvant postoperative chemotherapy were placed on the original neoadjuvant chemotherapy regimen. Demographic characteristics are presented in Supplementary Tables 1 and 4.
Post-CTX samples were also retrospectively reviewed for College of American Pathologists (CAP) scoring for histological treatment effect in response to neoadjuvant chemotherapy as a prognostic factor for individuals undergoing surgical resection for localized PDAC 14 . CAP scoring was performed by a specialist pancreatic cancer pathologist (F.B.). Only one individual out of a total of 32 evaluable participants receiving neoadjuvant chemotherapy exhibited a near-complete response (score 1), 19 exhibited a partial response (score 2), and the 12 remaining individuals exhibited no response (score 3; Supplementary Table 4).

RNA-seq profiling of chemo-naive and post-CTX PDAC-HD samples
RNA-seq was initially performed to define the neoplastic and stromal changes associated with neoadjuvant chemotherapy (Supplementary Table 6).  Tables 7 and 8). Importantly, the removal of LCMderived samples from the analysis did not change the significance of the results, suggesting that sample preparation was not a confounding factor in this study (Supplementary Table 8).
WGCNA identified 15 coordinately expressed gene programs (GPs) representing distinct biological processes that could discriminate chemo-naive and post-CTX samples (  Table 8). Nine core GPs encompassing the regulation of β-cell development, extracellular matrix organization, collagen formation, immunoregulatory interactions and chemokine signaling were significantly enriched in the post-CTX samples. By contrast, six GPs encompassing cell cycle checkpoints, resolution of sister chromatid cohesion, xenobiotics and phase I functionalization of compounds, cell-cell communication and Eph-ephrin, RAF and MAPK signaling pathways were significantly repressed in response to neoadjuvant chemotherapy. to 30-50% in conjunction with adjuvant chemotherapy, most individuals relapse within a median of 12. 8-21.6 months (refs. 5,8,9).
Individuals with borderline resectable disease seem to have a better survival benefit from neoadjuvant gemcitabine (GEM) with capecitabine or mFOLFIRINOX (mFFX) rather than with chemoradiotherapy [10][11][12] . Induction therapy may also increase resectability and improve survival in individuals with otherwise unresectable local disease 9,13,14 . Although second-line cytotoxic therapies are also delivered after disease progression or relapse in metastatic, locally advanced and postresection settings, response to treatment and overall survival are disappointing compared to other tumor types 5,8,15 .
Dissociated responses are observed in individuals with metastatic disease, in which some metastases respond to treatment whereas others remain stable or progress. Chemoradiotherapy may also be implicated in the differential prognosis of pathological treatment effects in individuals who have received neoadjuvant chemotherapy 14 .
Single-nucleus and spatial transcriptomic profiling of chemo-naive and post-therapy (chemoradiotherapy and chemotherapy) samples has identified distinct neoplastic cell phenotypes that exist in untreated samples and persist after therapy 16 . Systematic profiling of metastatic biopsies and matched organoid models has identified a continuum of transcriptional states spanning classical and basal-like phenotypes that are joined by an intermediate coexpressor (IC) or 'hybrid' transcriptional state 17 . Ex vivo studies have further revealed that transitions between classical and basal-like phenotypes may significantly impact drug responses, with basal-like states exhibiting broadly decreased sensitivity to chemotherapy 17 .
Despite these advances, the contribution of distinct neoplastic cell phenotypes to outcomes following neoadjuvant chemotherapy remains largely underexplored. Here, we identify GATA6-, KRT17-and cytochrome P450 3A (CYP3A)-coexpressing cells that are enriched in resected PDAC samples following neoadjuvant chemotherapy. We also reveal a tumor cell-intrinsic role for CYP3A-expressing detoxification pathways in the persistence of drug-tolerant cells in minimal residual disease.

The PDAC Heidelberg (PDAC-HD) sample cohort
To define the impact of neoadjuvant chemotherapy, we established a PDAC-HD cohort with validated PDAC tissue and defined pathological clinical characteristics, including long-term follow-up in 171 individuals comprising (1) chemo-naive PDAC tissue after primary resection in individuals with radiologically resectable tumors with most receiving adjuvant chemotherapy (n = 115) after resection and (2) postchemotherapy (post-CTX) PDAC tissue after resection following neoadjuvant therapy in individuals with radiologically borderline resectable or unresectable tumors (n = 56 (refs. 13,18) ; Fig. 1a). Individuals who received chemoradiotherapy were excluded from the study.

Transcriptomic subtypes and outcomes post-CTX
Initial indications have suggested that transcriptomic subtypes can prognosticate in individuals with resectable and metastatic disease 24,25 .
Subtyping analysis of chemo-naive PDAC-HD samples demonstrated that Moffitt and Notta subtypes were prognostic for PDAC-HD chemo-naive samples ( Fig. 1e and Extended Data Fig. 1b,c). Consistent with earlier findings, Bailey, Collison and Notta subtypes exhibited considerable overlap with the Moffitt classical and basal-like subtypes (Extended Data Fig. 1b). Notably, 71% (n = 45) of the PDAC-HD chemo-naive samples were identified as belonging to the classical subtype, whereas the remaining 29% (n = 19) were identified as basal-like.
Subtyping analysis of post-CTX resected PDAC-HD samples demonstrated that both classical and basal-like subtypes persisted after therapy (Fig. 1c). In comparison to chemo-naive PDAC-HD samples, the basal-like subtype was found to be significantly overrepresented in post-CTX samples compared to the classical subtype (P < 0.001, chi-squared test, two sided; Fig. 1d). Enrichment of the basal-like subtype in post-CTX samples was independent of clinical stage at time of treatment (Extended Data Fig. 1d). Survival analysis revealed that although Moffitt subtypes were prognostic for chemo-naive PDAC-HD samples (log-rank P = 0.016), they did not discriminate between good and poor outcomes in PDAC-HD post-CTX samples (log-rank P = 0.38; Fig. 1e and Extended Data Fig. 1c). Similarly, the Bailey (log-rank P = 0.43), Collisson (log-rank P = 0.68) and Notta (log-rank P = 0.061) subtyping schemes did not predict outcome in post-CTX PDAC-HD groups. To corroborate these findings, purity independent subtyping of tumors (PurIST) was applied to the PDAC-HD cohort 25 . Although high PurIST scores exhibited significant overlap with the Moffitt basal subtype, high and low PurIST scores similarly did not prognosticate in PDAC-HD post-CTX samples (log-rank P = 0.65: high 13 versus low 15). These findings, although representative of the participants included in this study, will require further validation in suitably matched independent cohorts.

GATA6 + and KRT17 + cell phenotypes persist post-CTX
The persistence of subtype-specific programs in post-CTX samples suggested that neoplastic cell phenotypes may contribute to outcomes following neoadjuvant chemotherapy. To characterize neoplastic cell phenotypes in PDAC-HD samples, we performed immunofluorescence (IF) staining using validated antibodies to GATA6, HNF1A, KRT5, KRT17, KRT81 and S100A2 (Fig. 3a-c and Supplementary Table 13). Biomarker expression was determined alone or in the context of KRT19 coexpression, a ductal biomarker that is significantly upregulated in atypical ductal cells and PDAC 16 (Fig. 3a-c). chemotherapy impacts the TME. a, Top, significantly enriched immune cell types. Middle, significantly expressed immunosuppressive genes. Bottom, significantly expressed immunostimulatory genes. Significance was determined by two-sided Wilcoxon rank-sum test adjusted for multiple testing (P ≤ 0.05); T H 1, type 1 helper T cell. b, Immune cell types that exhibit significant enrichment in chemo-naive (n = 32) and post-CTX (n = 33) samples. The heat map represents median immune cell enrichment, and the bar chart represents significance of enrichment as -log 10 (Wilcoxon rank-sum test two-sided P value adjusted for multiple testing (P adj )); HSC, hematopoeitc stem cell; aDC, activated dendritic cell; T CM , central memory T cell; GMP, granulocyte-monocyte progenitor; T reg , regulatory T cell; CMP, common myeloid progenitor. c, Bar charts showing significant enrichment of specific stromal signatures in post-CTX samples; myCAFs, myofibroblast-like CAFs; iCAFs, inflammatory CAFs. The significance is provided as -log 10 (Wilcoxon rank-sum test two-sided P value adjusted for multiple testing). The dotted line represents -log 10 (P adj ≤ 0.05). d, Volcano plots showing the enrichment of immunomodulatory and myofibroblastic cell signature genes in samples treated preoperatively with GEM or mFFX. Genes significantly enriched (log 2 (fold change) of >1 and -log 10 (P adj ) of >2) in post-CTX mFFX samples are shown. P adj represent the significance of a two-sided Wald test adjusted for multiple testing; FC, fold change. e, Kaplan-Meier survival analysis for high (greater than median) and low (less than median) macrophage M2 enrichment values in post-CTX (GEM and mFFX) and mFFX PDAC-HD samples. Participant numbers for each group are provided under 'Numbers at risk'. A log-rank (two-sided) P value of ≤0.05 is considered significant. f, Immune cell types that exhibit significant enrichment in mFFX post-CTX samples (n = 23) relative to GEM post-treated samples (n = 10). The heat map represents median immune cell enrichment, and the bar chart represents the significance of enrichment as -log 10 (Wilcoxon rank-sum test twosided P value adjusted for multiple testing). g, Immunostimulatory genes that are significantly and differentially expressed between post-CTX GEM (n = 10) and mFFX (n = 23) samples. The bar chart provides the significance of enrichment as -log 10 (Wilcoxon rank-sum test two-sided P value adjusted for multiple testing). Correlation heat map showing correlations between immunomodulatory factors in post-CTX samples. Pearson's correlations are shown in the plot. Significance was determined by two-sided Pearson's correlation test. P values were not adjusted for multiple testing. All correlations shown are significant.    GATA6 is broadly considered a surrogate marker of the classicallike subtype 27,29-34 . GATA6 protein expression by immunohistochemistry is predictive of better survival with 5-fluorouracil/folinic acid in the adjuvant (chemo-naive) setting but not with adjuvant GEM 33   Time (months) Survival probability 10  (although these biomarkers may not represent the full spectrum of basal/squamous cell states). GATA6 protein expression was determined by IF in chemo-naive and post-CTX PDAC-HD samples (Fig. 3b). This analysis demonstrated that the percentage of nuclear GATA6 expression was significantly elevated post-CTX (considering GEM and mFFX samples together) compared to the chemo-naive group. This trend of increased GATA6 expression was also observed when the percent expression of nuclear GATA6 in KRT19 + cells was considered. Stratifying post-CTX samples by chemotherapy regimen demonstrated that GATA6/KRT19 expression was significantly higher in GEM samples than in the mFFX group, suggesting that chemotherapy regimen may influence the number of GATA6-expressing cells (Fig. 3c). Imaging of representative classical GATA6 hi KRT17 low and basal-like KRT17 hi GATA6 low samples clearly demonstrated dominant biomarker expression in post-CTX samples ( Fig. 3d and Extended Data Fig. 3a). Taken together, these data revealed that GATA6 hi and KRT17 hi cell phenotypes persisted after chemotherapy.
To determine whether persistent GATA6 hi classical and KRT17 hi basal-like phenotypes in resected post-CTX samples are associated with outcomes, we categorized IF levels according to robust tertiles of protein expression and performed a survival analysis (Fig. 3e,f). High GATA6 protein expression in post-CTX samples was associated with significantly worse outcomes (Fig. 3e). To further understand the relationship between GATA6 protein expression and neoadjuvant chemotherapy, we examined whether the association between GATA6 hi expression and poor outcome was common or unique to either GEM or mFFX neoadjuvant chemotherapy. Importantly, we found that GATA6 hi expression was associated with significantly worse outcomes after mFFX, but not GEM, neoadjuvant chemotherapy ( Fig. 3e). A similar assessment of KRT17 IF expression tertiles demonstrated that KRT17 hi samples were significantly associated with worse outcomes after mFFX, but not GEM, neoadjuvant chemotherapy (Fig. 3f). Exploratory survival analysis using basal/squamous biomarkers KRT5, KRT81 and S100A2 revealed additional significant associations, although validation in independent cohorts will be required (Extended Data Fig. 4; please note that exploratory univariate analyses were not corrected for multiple testing). Collectively, these data demonstrate that GATA6 hi and KRT17 hi phenotypes persist post-CTX and contribute to poor outcomes in individuals after neoadjuvant chemotherapy.

Complex neoplastic heterogeneity defines post-CTX samples
Classical and basal-like subtype cell populations may coexist intratumorally 17,28 . Single-cell RNA-seq of chemo-naive and biopsied liver metastases has identified hybrid or IC neoplastic cell states that share biomarkers common to both the classical-like and basal-like/squamous subtypes, namely GATA6 and KRT17 (ref. 17). Enrichment analysis using gene signatures for single-cell classical (scClassical), scBasal and scIC cell states demonstrated that the scIC cell state is preferentially enriched in post-CTX samples (Extended Data Fig. 5a,b), suggesting that neoadjuvant chemotherapy may promote 'hybrid' cell states.
To assess whether GATA6 hi and KRT17 hi persister cell phenotypes are mutually exclusive or coexist in PDAC-HD samples, we performed multiplexed co-staining (colocalization) for GATA6, KRT17 and KRT19 ( Fig. 4a and Supplementary Table 9). This analysis revealed complex patterns of intratumoral expression. These included individual samples exhibiting dominant GATA6 hi KRT17 low or dominant KRT17 hi GATA6 low expression and samples with interspersed mosaic GATA6 hi KRT17 low and KRT17 hi GATA6 low expression foci within the same tissue section. A subset of individual samples also exhibited hybrid staining with GATA6 and KRT17 expressed in the same cells. Most strikingly, a number of samples exhibited a gradient of biomarker expression wherein dominant GATA6 hi KRT17 low , GATA6 hi KRT17 hi hybrid and dominant KRT17 hi GATA6 low expression foci were arranged serially within the same tissue section (Fig. 4a).
Deconvolution of biomarker intensities in whole-tissue sections allowed us to count individual cells according to biomarker expression with n obs = 7,545,622 cells designated as GATA6 + KRT17 -, GATA6 -KRT17 + or GATA6 + KRT17 + (Fig. 4c,d). Comparison of cell-type abundance between chemo-naive and post-CTX samples revealed that GATA6 + KRT17cells are enriched in chemo-naive samples, while , chemo-naive (n = 77) and post-CTX (n = 45) samples stained with GATA6 (red), HNF1A (red), KRT5 (red), KRT17 (red), KRT81 (red), S100A2 (red) and KRT19 (green) antibodies. b, Box plots showing relative whole-section protein expression of the indicated biomarkers in chemo-naive (n = 77) and post-CTX (n = 45) samples. Biomarker protein expression is considered alone or in the context of KRT19 coexpression. Kruskal-Wallis rank-sum test (two-sided) P values are provided at the top of each plot. c, Box plot showing relative whole-section protein expression of nuclear GATA6 in KRT19 + cells between GEM (n = 19) and mFFX (n = 25) post-CTX samples. The Kruskal-Wallis rank-sum test (two-sided) P value is provided at the top. d, Multiplexed IF images of representative classical and basal post-CTX samples stained with GATA6 (red), KRT17 (red) and KRT19 (green) antibodies. Representative images are presented in rows, with the leftmost image showing the entirety of the imaged region. The top right image and bottom right image show selected regions (i and ii) at increased magnification. e,f, Kaplan-Meier survival analysis for high, medium and low GATA6 and KRT17 protein expression tertiles in post-CTX PDAC-HD samples. Kaplan-Meier survival analyses for post-CTX samples representing combined (GEM and mFFX) treatment, GEM alone or mFFX alone are shown. Participant numbers for each group are provided under 'Numbers at risk'. A log-rank P value of ≤0.05 is considered significant. All box plots show the median (line), the interquartile range (IQR) between the 25th and 75th percentiles (box) and 1.5× the IQR ± the upper and lower quartiles. P values were not adjusted for multiple testing.  GATA6 + KRT17 + and GATA6 -KRT17 + cell populations are preferentially enriched post-CTX (Fig. 4c,d). Samples treated preoperatively with mFFX exhibited significantly greater numbers of GATA6 + KRT17 + cells than those treated preoperatively with GEM (Fig. 4e,f).
Survival analysis demonstrated no significant association between GATA6 + KRT17 + cell enrichment and outcome in either the chemo-naive or post-CTX groups (Fig. 4b). To disentangle the relationship between participant outcomes and GATA6/KRT17 cell phenotypes, post-CTX samples were dichotomized using GATA6 expression tertiles, as described earlier, and GATA6/KRT17 cell populations were compared. In the context of mFFX neoadjuvant chemotherapy, samples with low GATA6 expression (associated with longer survival) had increased numbers of GATA6 + KRT17cells, whereas samples with high GATA6 expression (associated with shorter survival) had increased numbers of GATA6 + KRT17 + hybrid cells but significantly fewer GATA6 + KRT17cells (Fig. 4g). GATA6 protein expression was not prognostic for GEM neoadjuvant chemotherapy, and we failed to observe a similar shift toward higher GATA6 + KRT17 + hybrid cell numbers or near loss of GATA6 + KRT17cells in the GATA6 hi samples treated preoperatively with GEM. These findings suggest that poor outcomes associated with GATA6 hi and KRT17 hi expression in post-CTX samples involve a higher relative    enrichment of GATA6 + KRT17 + and GATA6 -KRT17 + cells and near loss of GATA6 + KRT17cells.

Drug-tolerant persister cell phenotypes
Network analyses of GPs enriched in chemo-naive samples identified a subnetwork of coexpressed genes that comprised key pancreatic transcription factors GATA6, HNF4A, HNF1A, FOXA2 and FOXA3 and genes involved in drug metabolism, that is, CYP450 enzymes and phase I functionalization of compounds (Fig. 5a,b). The coexpression of this subnetwork of genes, while highly expressed in chemo-naive samples, was retained in a subset of post-CTX samples (Fig. 5a). Of particular interest was the identification of CYP450 family genes (CYP3A4, CYP3A5 and CYP2C9), which were coexpressed with transcription factors GATA6, HNF4A and NR1I2 (also known as pregnane X receptor (PXR); Fig. 5a,b). Transcription factors HNF4A and NR1I2 (PXR) have been shown to regulate both steady-state and substrate-induced expression of CYP3A5 in both classical and basal-like pancreatic cancer cell lines 37 . A reanalysis of existing data 31 demonstrated that GATA6 and HNF4A are required for the expression of not only CYP3A5 but also several other genes involved in drug detoxification (Fig. 5c). The CYP3A subfamily (encoded by 4 genes, CYP3A4, CYP3A5, CYP3A7 and CYP3A43) has been shown to metabolize common cytotoxic chemotherapeutics, including irinotecan, paclitaxel, docetaxel, anthracyclines, vinca alkaloids and tyrosine kinase inhibitors, such as erlotinib and gefitinib, in the liver and intestinal epithelia 37 . The contribution of these genes to extrahepatic drug tolerance is poorly understood [38][39][40][41] . Recent evidence, however, suggests that tumor cell-intrinsic expression of CYP3A5 may underpin intrinsic drug resistance in colorectal cancer and PDAC 37,42 . CYP3A5 activity has been associated with resistance to paclitaxel in pancreatic cancer cell lines 37 and resistance to irinotecan in colorectal cancer 42 . Given that irinotecan, a constituent of mFFX, is a substrate of CYP3A family proteins, we surmised that CYP3A family proteins might underpin the persistence of drug-resistant GATA6 hi phenotypes in mFFX post-CTX samples.
Given the tight relationship between CYP3A4/CYP3A5 and CYP3A7 gene expression, we stained for CYP3A5 protein (referred to herein as CYP3A) as a biomarker reflecting expression of the broader CYP3A network (Fig. 5d). CYP3A protein expression was initially assessed in normal pancreatic donor tissue, chemo-naive and post-CTX PDAC-HD samples (Fig. 5d, left, and Supplementary Table 5). Importantly, we observed a significant increase in CYP3A protein expression in samples treated preoperatively with mFFX, but not GEM, further supporting the role of CYP3A as an important mediator of mFFX drug resistance (Fig. 5e).
To determine whether CYP3A protein expression was prognostic in PDAC-HD samples, we performed a survival analysis using dichotomized CYP3A values representing high (highest 25% of IF values; n = 11 participants) and low (remainder of IF values; n = 33 participants) expression. This expression cutoff was used, as CYP3A values were negatively skewed ( Fig. 5f and Extended Data Fig. 6). In the post-CTX group, high CYP3A protein expression was associated with significantly worse outcomes. Mirroring the results obtained for GATA6, high CYP3A protein expression in mFFX, but not GEM, was associated with significantly worse outcomes. Exploratory analysis using antibodies specific to hENT1, an established biomarker of GEM resistance in PDAC, found no significant association between hENT1 expression and participant outcome in either setting, although validation in independent cohorts will be required [43][44][45] (Extended Data Fig. 6; as before, exploratory univariate analyses were not corrected for multiple testing).
PDOs resistant to irinotecan exhibited predominant Moffitt classical transcriptomic profiles, whereas five of the six most susceptible PDOs exhibited Moffitt basal transcriptomic profiles (Fig. 7e). Importantly, CYP3A protein expression, as determined by western blotting and IF, was positively correlated with irinotecan resistance in PDOs (Fig. 7a-e and Extended Data Fig. 8). Similarly, genes involved in drug metabolism and phase I functionalization of compounds were enriched in irinotecan-resistant PDOs, suggesting that coordinated networks of drug-metabolizing genes mediate irinotecan resistance (Fig. 7e).
MS analysis demonstrated that resistant and susceptible PDOs take up irinotecan. Exposure of PDOs to 2 μM irinotecan demonstrated that less than 1% of the intracellular irinotecan was converted to intracellular SN-38 (Extended Data Fig. 9f). Importantly, irinotecan-resistant PDOs showed significantly lower metabolic ratios (intracellular SN-38:intracellular irinotecan) than irinotecan-sensitive organoids, suggesting that the resistance phenotype involves lower conversion of irinotecan to SN-38 (Fig. 8b,c), effects that were independent of proliferation (Extended Data Fig. 9a-c). In addition, resistant PDOs tended to show higher SN-38 concentrations in the supernatant, suggesting that the active transport of SN-38 may also   50 values are provided for the indicated treatments. e, Alternative models to explain drug tolerance and persistence in post-CTX samples. Intrinsic resistance: treatment-mediated selection of preexisting drug-tolerant phenotypes may shape residual disease. This process may involve non-genetic mechanisms, including intrinsic and/or drug-induced expression of drug-detoxifying genes and/or a transition toward basal-like or 'hybrid' states from predominant classical states due to intrinsic neoplastic plasticity. Acquired resistance: rare subclones acquire a drug-resistant driver alteration before or during therapy. These resistant clones expand and eventually drive relapse due the clonal acquisition of the preexisting drug-resistant mechanism.  contribute to irinotecan resistance (Fig. 8b,c). This was especially true for the irinotecan-resistant PDO h20, which exhibited the highest concentration of SN-38 in the supernatant (Fig. 8b,c and Extended Data Fig. 9g). The CYP3A-mediated metabolite APC was inconclusively quantifiable in the majority of PDOs; however, physiologically relevant amounts were recorded in the supernatant, suggesting that APC is actively transported from PDOs (Fig. 8c). PDO h20 exhibited quantifiable levels of APC in both cells and supernatant and had the highest supernatant concentration of APC for the PDOs tested. Notably, PDO h20 was derived from a participant who received neoadjuvant mFFX and was enriched for GATA6 -CYP3A + KRT17 + and GATA6 -CYP3A + KRT17cell phenotypes. PDO h3, which shows relative susceptibility to irinotecan and was derived from a metastatic liver recurrence after first-line GEM monotherapy, was the only other organoid to produce high supernatant concentrations of APC.
To establish a direct role for CYP3A in irinotecan drug resistance, we treated selected resistant and susceptible PDOs with optimal concentrations of irinotecan in combination with either ketoconazole or cobicistat, potent pan-CYP3A inhibitors (Fig. 8d and Extended Data Fig. 9d,e,h). Treatment of resistant PDOs with irinotecan in combination with ketoconazole significantly increased sensitivity to irinotecan. Importantly, the ability of ketoconazole to sensitize resistant cells to irinotecan was not observed when the active metabolite SN-38 was used as a substitute (Fig. 8d and Extended Data Fig. 9h). CYP3A activity has also been shown to mediate paclitaxel resistance in pancreatic cancer cell lines 37 . Consistent with these findings, ketaconzole-mediated inhibition of CYP3A increased sensitivity to paclitaxel in irinotecan-resistant PDOs. Together these findings implicate CYP3A as an important mediator of irinotecan drug resistance in PDAC.

Discussion
Drug resistance to standardized chemotherapy remains a significant challenge for PDAC, with relapse occurring in most individuals 5,8,9 . It is currently unclear how drug-resistant clones emerge and evolve during therapy. Two major mechanisms have been proposed to explain drug resistance (Fig. 8e). The first of these supposes that drug resistance arises from rare subclones that acquire a drug-resistant driver alteration (genetic) during therapy. In response to therapy, these resistant clones expand and eventually drive relapse due to the clonal acquisition of the preexisting drug-resistant mechanism. The second model proposes that drug-tolerant cells, or 'persisters', initially emerge from a preexisting subpopulation of cells that do not harbor classical drug-resistant driver alterations (non-genetic). These cells survive initial drug treatments by an epigenetic and/or transcriptional adaption that allows a drug-tolerant, slow-cycling 'persister' state to emerge. Clinically, this persister state resembles minimal residual disease from which relapse can occur if treatment is discontinued or if persister cells acquire a drug-resistant driver alteration due to continued drug therapy 48 .
The evidence presented herein strongly suggests that preexisting subpopulations of GATA6-, KRT17-and CYP3A-coexpressing cells survive the initial rounds of mFFX treatment to emerge as a persistent population of drug-resistant cells. Recent evidence has revealed that CYP3A expression in pancreatic cancer cells may mediate resistance to paclitaxel and tyrosine kinase inhibitors 37 . Here, we extend these findings to demonstrate that CYP3A activity is a mediator of irinotecan resistance. Pharmacological inhibition of CYP3A activity using the pan-CYP3A inhibitor ketoconazole sensitized resistant PDOs to irinotecan but not the active metabolite SN-38, directly implicating CYP3A activity in irinotecan resistance. Strikingly, GATA6-, KRT17-and CYP3A-coexpressing 'hybrid' cells were enriched in resistant PDOs, further implicating these cell phenotypes in chemotherapy resistance.
Recent studies have demonstrated that classical and basal-like phenotypes exist as a gene expression continuum with a 'hybrid' or IC state acting as a plastic intermediate 17 . Modulation of cell state by the addition of stromal cues or chemotherapy has been shown to shift this continuum from a classical state toward an induced 'hybrid' and/or basal-like (mesenchymal) state. Accordingly, the enrichment of 'hybrid' and basal-like states at the expense of classical-like states observed in post-CTX samples may reflect this underlying plasticity. Understanding how resistant neoplastic cell populations emerge following chemotherapy will be critical for the development of effective first-line and adjuvant therapies. Future studies comparing matched before and after treatment samples should provide critical new insights into the genetic and non-genetic mechanisms underpinning therapy resistance.
The number of individuals pretreated with chemotherapy is increasing. Neoadjuvant chemotherapy may improve outcomes for individuals with locally advanced and borderline cancer, and studies in resectable individuals are ongoing [10][11][12][13][14][15] . Resected tumor material represents a unique opportunity to tailor adjuvant treatment to residual cancer cells. Development of effective treatments for the persister cells identified in this study is therefore of the utmost importance. As increased CYP3A expression may increase resistance to both irinotecan and nab-paclitaxel, alternative chemotherapy combinations without these drugs, such as GEM-capecitabine, might be considered. Moreover, the enrichment of suppressive macrophage signatures and/or specific immunoregulators in neoadjuvant samples suggests that combination approaches targeting myeloid cells 49,50 and/or major inhibitory checkpoint molecules 51 may also enhance benefit.

Participant characteristics
The study was approved by the Ethics Committee of Heidelberg University for use of pancreatic cancer tissue (project numbers S-018/2020, S-708/2019 and S-083/2021). Participants were selected according to the following criteria: (1) individuals with locally advanced borderline unresectable PDAC (without metastases) who received surgical resection after mFFX or GEM-based therapy without any chemoradiation neoadjuvant therapy and (2) individuals with resectable PDAC at presentation who had upfront surgical resection without prior chemotherapy and/or chemoradiation (adjuvant chemo-naive). Participants after resection could receive adjuvant chemotherapy, but those who had chemoradiation were again excluded. Consecutive pseudoanonymized participants with prior consent and with usable samples of sufficient quality for investigation were identified based on the eligibility criteria. Additional pseudoanonymized IDs were then ascribed to each participant for experimental investigation, so blinded to all of the investigators and investigations. Clinical correlates were only ascribed once the experimental procedure results were obtained.
Cryopreserved and FFPE tissues were processed as described previously 52 . Hematoxylin and eosin (H&E) staining of cryopreserved and FFPE tissue was examined by a specialist pancreas pathologist (F.B.) for staging 53 and tumor content in the Department of Pathology. RNA-seq was performed on cryopreserved PDAC tissues, and multiplexed IF was performed on PDAC FFPE tissues.
Article https://doi.org/10.1038/s43018-023-00628-6 Medium was changed twice a week. For continued passaging, organoids were recovered from the Matrigel using cold Cell Recovery Solution (Corning), further dissociated into single cells using TrypLE (Gibco) and embedded with fresh Matrigel. Organoid cell lines were checked for KRAS mutations by DNA Sanger sequencing.
The following PCR conditions were used as previously described 54

Pharmacological assay of organoids
Organoids were dissociated into single cells before being plated in 10 μl of Matrigel with 1,000 cells per well into white 96-well plates (Greiner). Chemotherapeutics were tested in triplicate: 5-fluorouracil (Sigma), irinotecan (Sigma) and oxaliplatin (Selleckchem) with concentrations ranging from 1.0 × 10 −7 to 1.0 × 10 −3 mol liter -1 , GEM (Sigma), paclitaxel (Selleckchem) and SN-38 (Sigma) with concentrations ranging from 1.0 × 10 −10 to 1.0 × 10 −6 mol liter -1 and ketoconazole (Sigma) and cobicistat (MCE) ranging from 1.0 × 10 −6 to 2.0 × 10 −5 mol liter -1 . In combination treatment, 5.0 × 10 −6 mol liter -1 ketoconazole or cobicistat was combined with a corresponding dose of irinotecan, paclitaxel or SN-38. Reagents were dissolved in DMSO and added 72 h after plating. All treatment wells were normalized to 0.25% DMSO. After 96 h of treatment, cell viability was assessed using the CellTiter-Glo 3D cell viability assay (Promega), as per the manufacturer's instructions, on a FLUOstar plate reader (BMG Labtech). Therapeutic results (viability versus dose) were analyzed with GraphPad software. A four-parameter log-logistic function with an upper limit equal to the mean of the DMSO values was fit to the drug response curve and IC 50 value calculated.

CYP3A activity assay in organoids
After 96 h of culture with or without CYP3A inhibitor, the pan-CYP3A activity of organoids was measured using a P450-Glo kit (Promega), as per the manufacturer's instructions, on a FLUOstar plate reader (BMG Labtech). All activity measurements from each well were normalized to cell numbers using the CellTiter-Glo assay mentioned above.

Compound analysis by UPLC-MS/MS
Organoids were treated for 96 h with 0.1 μM or 2 μM irinotecan (<IC 50 ) to prevent selection bias from killing most of the organoid population. Non-lethal drug concentrations were chosen to allow metabolic phenotyping of intracellular and extracellular (supernatant) concentrations of irinotecan, SN-38, APC and NPC.
Intracellular and supernatant concentrations were quantified with a validated UPLC-MS/MS assay following the guidelines of the European Medicines Agency (EMA) and US Food and Drug Administration (FDA) and on bioanalytical method validation (http://www.ema.europa. eu/docs/en_GB/document_library/Scientific_guideline/2011/08/ WC500109686.pdf and http://www.fda.gov/downloads/Drugs/Guid-anceComplianceRegulatoryInformation/Guidances/ucm070107.pdf). Each of four performed validation runs included blank and internal standard controls, seven calibration samples (twofold) and four quality control (QC) concentrations (sixfold). The assays fully complied with the applicable parts of the recommendations of the US FDA and EMA.
Optimized MS/MS parameters for the detection of irinotecan and its metabolites can be found in Supplementary Table 3. The calibrated range was 10-10,000 pg ml -1 , showing linear regression coefficients of >0.99. Overall accuracies (interday and intraday) were between 90.0 and 114.0% with a corresponding precision of <15%. A Xevo-TQ-S tandem mass spectrometer (Waters) coupled to an Acquity classic UPLC (Waters) and equipped with a heated electrospray ionization source was used for quantification. Determinations were performed with selected reaction monitoring using collision-induced dissociation with argon in the positive ion mode. Chromatographic separation was performed on a BEH C18 column (50 × 2.1 mm and 1.7 μm, Waters) with a linear gradient from 5 to 75% acetonitrile (ACN) + 0.1% formic acid in 1.2 min (corresponding decrease of aqueous eluent: 19:1 water:ACN + 0.1% formic acid) at a flow rate of 0.5 ml min -1 .
Organoids were lysed using 300 μl of 5% aqueous NH 3 , from which 100 μl was withdrawn for analysis (study samples). Calibration and QC samples were prepared by spiking 25 μl of calibration or QC spike solution (corresponding sample concentration calibration: 10, 30, 100, 300, 1,000, 3,000 and 10,000 pg ml -1 ; QC: 10, 30, 3,750 and 7,500 pg ml -1 ) into 100 μl of cell lysate. All samples (lysed organoids or supernatants) were spiked with 25 μl of internal standard solution (irinotecan-D 6 and SN-38-D 6 ), and study samples were additionally spiked with 25 μl of blank solution for volume compensation. Irinotecan and metabolites were extracted with protein precipitation using 300 μl of acetonitrile containing 0.1% formic acid. After shaking and centrifugation at 16,100g for 5 min, extracts were transferred to 96-well collection plates and evaporated with a blowdown evaporator (Ultravap, Porvair Sciences). After addition of 100 μl of a mixture of water:ACN (3:1 (vol/vol)) containing 0.1% formic acid, 20 μl was injected into the UPLC-MS/MS system for analysis.
Measured drug concentrations in cellular lysates were normalized to protein content of the sample, which were evaluated using a commercial BCA assay kit (Pierce BCA Protein Assay kit, Thermo Scientific). Measurements were performed as previously described 56 . Specifically, Article https://doi.org/10.1038/s43018-023-00628-6 as noted in that paper, a calibration curve with nine standard sample concentrations of BSA (0-2,000 μg ml -1 ) was prepared. Wells of a 96-well plate were loaded with 8 μl of BSA standard or the respective sample and 64 μl of the working solution (a mixture of reagent A and reagent B contained in the kit). After 30 min of light-protected incubation at 37 °C, absorption was read at 562 nm using a Spectramax plate reader (Molecular Devices), and protein content was calculated.
Tissue processing and next-generation sequencing LCM was performed on cryopreserved tissue samples from individuals after resection and before any adjuvant therapy, according to previously described methods 29,57 . RNA from bulk and LCM tumor specimens was isolated using an AllPrep DNA/RNA/miRNA Universal kit (Qiagen). RNA samples with an RNA integrity number of ≥8, a 28S/18S ratio of ≥1.0 and a DV200 of >70% were considered suitable for high-throughput sequencing. High-throughput sequencing library preparation was performed with an Illumina TruSeq stranded mRNA kit (Illumina, 20020595) with IDT unique dual indices (Illumina, 20022371) following the manufacturer's recommendations. Five hundred nanograms of total RNA was used as input. For each library, at least 57 million mapped reads were produced for downstream analysis.

IF assays
IF staining was accomplished using 4-μm-thin FFPE tissue sections, as described earlier 58 . Whole sections were captured using a TissueGnostics Fluorescence Imaging System (TissueGnostics), with a fluorescence microscope unit (Observer. Z1, Zeiss) and a Lumencor Sola SE III 365-to 730-nm LED light source (AHF Analysentechnik). Captured images were analyzed using StrataQuest software version 7.0 (TissueGnostics), which calculated the intensity of the fluorescence signals in each single cell within individual tissue sections. Cells containing high-intensity target signals were selected by partitioning cell fluorescent scattergrams against DAPI or other referenced co-stained target signals into upper expression quantiles. The expression level of each marker protein was calculated as the percentage of stained cells relative to DAPI and/or stated reference protein. Biomarker protein expression in individual cells was gated into high, medium or low expression tertiles, and cell-specific expression was quantified for each sample. The antibodies used in this study are listed in Supplementary Tables 13 and 14.

RNA-seq analysis
RNA-seq data were aligned using STAR version 2.5.3a via a DKFZ internal next-generation sequencing data-processing pipeline 59 . Downstream transcriptomic analysis was performed as previously described 24,[26][27][28] . Briefly, merged count data obtained from the DKFZ processing pipeline and representing chemo-naive and post-CTX samples were batch corrected using the sva R package. Batch-corrected count data were subsequently logR transformed using the DESeq2 R package to generate normalized gene expression values. The logR-normalized data were used for all downstream analyses unless otherwise specified.

Subtyping analysis
Subtyping analysis was performed on samples using gene expression signatures representing Moffitt 24,46 , Collisson 26 , Bailey 27 and Notta 28 subtypes. To classify samples according to subtype, logR-normalized gene expression values were clustered using the ConsensusClusterPlus R package. Samples were subsequently assigned to a specific subtype based on the results of the clustering analysis. Heat maps representing subtype clusters and showing representative subtype-specific genes were generated using the ComplexHeatmap R package. The PurIST algorithm was used as an orthogonal measure of basal-like subtype status in samples and was performed using normalized gene expression data, as previously described 25 . PurIST scores approaching 1 indicate an increased likelihood that the sample is basal like.
Gene expression signatures representing Bailey GPs, single-cell states or malignant lineage or state programs were obtained from Bailey et al. 27 , Raghavan et al. 17 and Hwang et al. 16 , respectively. These gene expression signatures were used to cluster samples and/or generate gene set enrichment scores. Gene set enrichment scores were generated using the GSVA R package. t-SNE analysis of gene expression or gene set enrichment scores was performed using the Rtsne R package. Volcano plots were generated using the ggplot2 R package. These plots represent the set of differentially expressed genes between post-CTX samples treated with either GEM or mFFX. Hwang et al. 16 malignant lineage and state program genes significantly enriched in the indicted sample group (log 2 (fold change) > 1 and -log 10 (adjusted P value) > 2) were highlighted in relevant plots using the gghighlight R package. The ggstatsplot R package was used to assess the significance of single-cell subtype signature enrichment.
WGCNA WGCNA, as implemented by the WGCNA R package, was performed using logR-normalized gene expression data. WGCNA was performed as previously described 27 . Specifically, 'WGCNA clusters genes into network GPs using a topological overlap measure (TOM) that represents a highly robust measure of network interconnectedness and provides a measure of connection strength between two adjacent genes and all other genes in a network. Genes were clustered using 1 - TOM as the distance measure and GPs defined as branches of the resulting cluster tree using a dynamic branch-cutting algorithm. Module eigengene values were used as a measure of GP expression in each sample.' Module eigengene values were used to identify GPs significantly enriched in either chemo-naive or post-CTX samples. Pathway enrichment analysis (see Supplementary Table 7 for the list of genes representing each GP) was performed using the clusterProfiler and ReactomePA R packages. Networks representing GPs (Fig. 4b) were generated using the Reactome FI Cytoscape plugin 8.05 (https://apps.cytoscape.org/apps/ reactomefiplugin) in Cytoscape version 3.9.1 (https://cytoscape.org) 60 .

Stromal cell enrichment analysis
Stromal cell type and/or phenotype enrichment in samples was estimate using xCell 61 , as implemented by the immunedeconv R package (Fig. 2). The logR-normalized gene expression values were used to obtain xCell enrichment scores. Enrichment scores for gene signatures representing Grünwald subTMEs 20 , Öhlund CAF phenotypes 62 or Hwang et al. 16 fibroblast programs were generated using the GSVA R package. Volcano plots were generated as described above using the set of differentially expressed genes between post-CTX samples treated with either GEM or mFFX. Immunomodulatory gene expression values were compared between samples as indicated. Correlations between immunomodulatory factors were generated and visualized using the corrplot R package.

Multiplexed IF cell-type analysis
GATA6, KRT17 and CYP3A cell counts obtained from StrataQuest image processing, as described above, were used for enrichment analysis. Bar and pie statistical plots were generated from individual cell counts using the ggstatsplot R package. Ternary plots were generated using the ggtern R package 63 . The relative percent enrichment of GATA6, KRT17 and/or CYP3A protein expression was calculated by dividing individual cell-type counts by the sum of all cell-type counts in each sample and multiplying by 100. Bar plots representing percent tumor cell enrichment were generated using the ggpubr R package.

Survival analysis
Survival was estimated using the Kaplan-Meier method, as implemented by the survminer R package. Participants were stratified by transcriptomic subtype or protein expression values as indicated. Forest plots were generated using the ggforestplot R package.

Statistics and reproducibility
No statistical methods were used to predetermine sample sizes, but our sample sizes are similar to those published in previous publications 16,17,24,26,27 . Participant selection was performed blind to clinical variables. Data collection and analysis were not performed blind to the conditions of the experiment. Informed consent was obtained from the participants of this study. Research findings do not apply to one sex or gender only. The gender of each participant was collected by consent and was self-reported. Information on gender is provided in Supplementary Tables 1, 2, 4 and 5. For tissue-based findings, 171 unique participants were included in the study, with 100 males and 71 females taking part. For RNA-seq, 56 males and 41 females were included in the analysis ( Supplementary Tables 1 and 4). For multiplexed IF ( Supplementary Tables 1 and 5), 71 males and 51 females were included in the analysis. For PDOs, organoids derived from 11 males and 20 females were included in the analysis (Supplementary Table 10). No gender-based analyses are shown, and no significant associations with gender were observed in the data. All source data comprise a participant identifier that can be used to disaggregate the data based on gender. All experiments are representative of at least three independent biological experiments. H&E and IF images for samples or PDOs are representative of at least three independent IF experiments on the same region of interest or PDO. No data points were excluded from the analyses. Data distributions were assumed to be normal, but this was not formally tested. P values of ≤0.05 were considered significant. Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Reporting summary
Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.

Data availability
All data relevant to this study are available from the corresponding authors upon request. All processed data, including normalized expression data for participant samples and PDOs, multiplexed IF and experimental results, are provided in the Supplementary Tables. RNA-seq data are available at the European Genome-Phenome Archive under accession number EGAS00001007143. RNA-seq data published in Hwang et al. 16 and Brunton et al. 31 were obtained from Gene Expression Omnibus accession number GSE202051 and BioProject accession number PRJNA630992, respectively. Source data are provided with this paper.

Code availability
No custom algorithms were used in this study. Analyses used existing code made available through referenced software packages. Code used in this study is available upon request.  Fig. 2 | Enrichment of subtype-specific gene programs (GP) and neoplastic cell phenotypes in chemo-naïve and post-CTX samples. a) Heatmap and bar plot of previously defined Bailey GP signatures that are significantly enriched in either the Classical/Progenitor/Immunogenic, Basallike/Squamous or ADEX subtypes. Heatmap z scores represent changes in median GP scores between chemo-naïve (n = 64) and post-CTX (n = 33) samples. Bar plots represent -Log 10 (Kruskal Wallis rank sum test two-sided P-value adjusted for multiple testing (pAdj)). b) Single sample gene set enrichment analysis using GP signatures (as shown in a) followed by tSNE. Sample clusters are identical between tSNE plots with individual samples highlighted according to treatment type, chemotherapy regimen and indicated GP score. Analyses for chemo-naïve (n = 64) and post-CTX (n = 33) samples are shown. c) Bar charts showing the differential enrichment of specific neoplastic cell populations in chemo-naïve and post-CTX samples as indicated. The significance of the enrichment is provided as -Log10(Wilcoxon rank sum test two-sided pAdj-value adjusted for multiple testing). The dotted line represents -Log10(pAdj ≤ 0.05). d) Volcano plots showing the differential enrichment of genes associated with the indicated neoplastic cell populations in post-CTX samples treated with either GEM or mFFX. Genes significantly enriched (Log 2 Fold change > 1 and -Log 10 (pAdj) >2) in samples treated pre-operatively with GEM or mFFX are shown. pAdj represent the significance of a two-sided Wald test adjusted for multiple testing.