Advertisement

Single-Cell Analysis Reveals Transcriptomic Features of Drug-Tolerant Persisters and Stromal Adaptation in a Patient-Derived EGFR-Mutated Lung Adenocarcinoma Xenograft Model

Open AccessPublished:December 16, 2022DOI:https://doi.org/10.1016/j.jtho.2022.12.003

      Abstract

      Introduction

      Targeted therapies require life-long treatment, as drug discontinuation invariably leads to tumor recurrence. Recurrence is mainly driven by minor subpopulations of drug-tolerant persister (DTP) cells that survive the cytotoxic drug effect. In lung cancer, DTP studies have mainly been conducted with cell line models.

      Methods

      We conducted an in vivo DTP study using a lung adenocarcinoma patient-derived xenograft tumor driven by an EGFR mutation. Daily treatment of tumor-bearing mice for 5 to 6 weeks with the EGFR inhibitor erlotinib markedly shrunk tumors and generated DTPs, which were analyzed by whole exome, bulk population transcriptome, and single-cell RNA sequencing.

      Results

      The DTP tumors maintained the genomic clonal architecture of untreated baseline (BL) tumors but had reduced proliferation. Single-cell RNA sequencing identified a rare (approximately 4%) subpopulation of BL cells (DTP-like) with transcriptomic similarity to DTP cells and intermediate activity of pathways that are up-regulated in DTPs. Furthermore, the predominant transforming growth factor-β activated cancer-associated fibroblast (CAF) population in BL tumors was replaced by a CAF population enriched for IL6 production. In vitro experiments indicate that these populations interconvert depending on the levels of transforming growth factor-β versus NF-κB signaling, which is modulated by tyrosine kinase inhibitor presence. The DTPs had signs of increased NF-κB and STAT3 signaling, which may promote their survival.

      Conclusions

      The DTPs may arise from a specific preexisting subpopulation of cancer cells with partial activation of specific drug resistance pathways. Tyrosine kinase inhibitor treatment induces DTPs revealing greater activation of these pathways while converting the major preexisting CAF population into a new state that may further promote DTP survival.

      Keywords

      Introduction

      Although EGFR tyrosine kinase inhibitors (TKIs) improve the survival of patients with lung cancer carrying EGFR-sensitizing mutations,
      • Hsu W.H.
      • Yang J.C.
      • Mok T.S.
      • Loong H.H.
      Overview of current systemic management of EGFR-mutant NSCLC.
      ,
      • Greenhalgh J.
      • Boland A.
      • Bates V.
      • et al.
      First-line treatment of advanced epidermal growth factor receptor (EGFR) mutation positive non-squamous non-small cell lung cancer.
      most patients relapse after 1 and 2 years. Approximately 50% of therapy resistance to first- and second-generation EGFR TKIs is due to an emergence of a secondary T790M EGFR mutation.
      • Chong C.R.
      • Jänne P.A.
      The quest to overcome resistance to EGFR-targeted therapies in cancer.
      • Lim S.M.
      • Syn N.L.
      • Cho B.C.
      • Soo R.A.
      Acquired resistance to EGFR targeted therapy in non-small cell lung cancer: mechanisms and therapeutic strategies.
      • Pallis A.G.
      • Serfass L.
      • Dziadziusko R.
      • et al.
      Targeted therapies in the treatment of advanced/metastatic NSCLC.
      Third-generation EGFR TKIs overcome T790M resistance, with disease control rates greater than 90%.
      • Ahn M.J.
      • Tsai C.M.
      • Shepherd F.A.
      • et al.
      Osimertinib in patients with T790M mutation-positive, advanced non-small cell lung cancer: long-term follow-up from a pooled analysis of 2 phase 2 studies.
      Nevertheless, patients eventually become resistant to these drugs, as well.
      • He J.
      • Huang Z.
      • Han L.
      • Gong Y.
      • Xie C.
      Mechanisms and management of 3rd-generation EGFR-TKI resistance in advanced nonsmall cell lung cancer (review).
      ,
      • Passaro A.
      • Jänne P.A.
      • Mok T.
      • Peters S.
      Overcoming therapy resistance in EGFR-mutant lung cancer.
      Some EGFR-mutant tumors have a minor population of T790M co-mutated cancer cells, which become dominant during first- and second-generation EGFR TKI therapy.
      • Su K.Y.
      • Chen H.Y.
      • Li K.C.
      • et al.
      Pretreatment epidermal growth factor receptor (EGFR) T790M mutation predicts shorter EGFR tyrosine kinase inhibitor response duration in patients with non-small-cell lung cancer.
      Rare cancer cells can also survive drug treatment without a genetic alteration, by entering a reversible state of drug tolerance.
      • Echeverria G.V.
      • Ge Z.
      • Seth S.
      • et al.
      Resistance to neoadjuvant chemotherapy in triple-negative breast cancer mediated by a reversible drug-tolerant state.
      ,
      • Rehman S.K.
      • Haynes J.
      • Collignon E.
      • et al.
      Colorectal cancer cells enter a diapause-like DTP state to survive chemotherapy.
      These drug-tolerant persisters (DTPs) may preexist or may stochastically be treatment induced. The DTPs typically divide slowly and eventually acquire genetic changes that enable better growth leading to tumor relapse. Understanding the DTP state and its origin is critical to inform strategies that reduce risk of tumor relapse.
      • Mikubo M.
      • Inoue Y.
      • Liu G.
      • Tsao M.S.
      Mechanism of drug tolerant persister cancer cells: the landscape and clinical implication for therapy.
      Most lung cancer DTP modeling has been performed using the PC9 cell line, which contains a TKI-sensitive EGFR exon 19 deletion. In this model, TKI-resistant cells can emerge from rare preexisting T790M-comutated cancer cells. Nevertheless, when these cells arise late, they are more drug resistant than those arising earlier.
      • Hata A.N.
      • Niederst M.J.
      • Archibald H.L.
      • et al.
      Tumor cells can follow distinct evolutionary paths to become resistant to epidermal growth factor receptor inhibition.
      This additional resistance is now known to reflect TKI induction of a nongenetic DTP phenotype that is reversible,
      • Sharma S.V.
      • Lee D.Y.
      • Li B.
      • et al.
      A chromatin-mediated reversible drug-tolerant state in cancer cell subpopulations.
      • Oren Y.
      • Tsabar M.
      • Cuoco M.S.
      • et al.
      Cycling cancer persister cells arise from lineages with distinct programs.
      • Kurppa K.J.
      • Liu Y.
      • To C.
      • et al.
      Treatment-induced tumor dormancy through YAP-mediated transcriptional reprogramming of the apoptotic pathway.
      similar to patient resensitization to TKIs after a drug holiday.
      • Kurata T.
      • Tamura K.
      • Kaneda H.
      • et al.
      Effect of re-treatment with gefitinib (‘Iressa’, ZD1839) after acquisition of resistance.
      • Yano S.
      • Nakataki E.
      • Ohtsuka S.
      • et al.
      Retreatment of lung adenocarcinoma patients with gefitinib who had experienced favorable results from their initial treatment with this selective epidermal growth factor receptor inhibitor: a report of three cases.
      • Shaffer S.M.
      • Dunagin M.C.
      • Torborg S.R.
      • et al.
      Rare cell variability and drug-induced reprogramming as a mode of cancer drug resistance.
      The DTP phenotype involves epigenetic changes to histone methylation and acetylation that remodel the transcriptome.
      • Sharma S.V.
      • Lee D.Y.
      • Li B.
      • et al.
      A chromatin-mediated reversible drug-tolerant state in cancer cell subpopulations.
      ,
      • Kurppa K.J.
      • Liu Y.
      • To C.
      • et al.
      Treatment-induced tumor dormancy through YAP-mediated transcriptional reprogramming of the apoptotic pathway.
      ,
      • Guler G.D.
      • Tindell C.A.
      • Pitti R.
      • et al.
      Repression of stress-induced LINE-1 expression protects cancer cell subpopulations from lethal drug exposure.
      Transcription factors such as NF-κB, TEAD, YAP, NRF2, and aldehyde dehydrogenases help establish the resistant state, potentially through increases in mitochondrial respiration, antioxidant responses, fatty acid metabolism (FAM), epithelial-mesenchymal transition, suppression of apoptosis, and expression of repetitive DNA sequences.
      • Mikubo M.
      • Inoue Y.
      • Liu G.
      • Tsao M.S.
      Mechanism of drug tolerant persister cancer cells: the landscape and clinical implication for therapy.
      ,
      • Oren Y.
      • Tsabar M.
      • Cuoco M.S.
      • et al.
      Cycling cancer persister cells arise from lineages with distinct programs.
      ,
      • Kurppa K.J.
      • Liu Y.
      • To C.
      • et al.
      Treatment-induced tumor dormancy through YAP-mediated transcriptional reprogramming of the apoptotic pathway.
      ,
      • Guler G.D.
      • Tindell C.A.
      • Pitti R.
      • et al.
      Repression of stress-induced LINE-1 expression protects cancer cell subpopulations from lethal drug exposure.
      ,
      • Raha D.
      • Wilson T.R.
      • Peng J.
      • et al.
      The cancer stem cell marker aldehyde dehydrogenase is required to maintain a drug-tolerant tumor cell subpopulation.
      One challenge that has limited translation of DTP knowledge into curative precision oncology has been a poor understanding of DTP origins. Different groups have reached differing conclusions about whether DTPs arise stochastically or from specific populations of preexisting cells.
      • Sharma S.V.
      • Lee D.Y.
      • Li B.
      • et al.
      A chromatin-mediated reversible drug-tolerant state in cancer cell subpopulations.
      • Oren Y.
      • Tsabar M.
      • Cuoco M.S.
      • et al.
      Cycling cancer persister cells arise from lineages with distinct programs.
      • Kurppa K.J.
      • Liu Y.
      • To C.
      • et al.
      Treatment-induced tumor dormancy through YAP-mediated transcriptional reprogramming of the apoptotic pathway.
      An additional challenge lies in dissecting potential heterogeneity among DTPs from different tumors or even within the same tumor. For example, osimertinib-induced PC9 DTPs can be separated into two subsets with distinct cycling properties and resistance mechanisms.
      • Oren Y.
      • Tsabar M.
      • Cuoco M.S.
      • et al.
      Cycling cancer persister cells arise from lineages with distinct programs.
      More work with patient samples and in vivo settings is needed to validate cell line findings and to discover potential patient-relevant mechanisms absent from in vitro models. This can be challenging, as exemplified by a recent study that had to pool lung cancer samples with different drivers and targeted therapies to create pan-treatment–naive and pan-residual disease (pan-RD) cohorts, which yielded a phenotypic rather than mechanistic description of potential DTPs.
      • Maynard A.
      • McCoach C.E.
      • Rotow J.K.
      • et al.
      Therapy-induced evolution of human lung cancer revealed by single-cell RNA sequencing.
      Here, we characterize the in vivo DTP state generated by the EGFR TKI erlotinib, using an EGFR-mutant patient-derived xenograft (PDX) that recapitulates patient sensitivity to first-generation EGFR TKIs.
      • Stewart E.L.
      • Mascaux C.
      • Pham N.A.
      • et al.
      Clinical utility of patient-derived xenografts to determine biomarkers of prognosis and map resistance pathways in EGFR-mutant lung adenocarcinoma.
      Through transcriptome profiling, we construct a new model for EGFR TKI DTP generation and survival. Our data support DTPs arising from a specific preexisting DTP-like (DTPL) population with partial elevation in pathways associated with drug resistance, with the TKI completing induction into a DTP state with higher expression/activation of these pathways. These findings were partly reproduced in another PDX model bearing EGFR exon 19 deletion and T790M mutation and treated by osimertinib. Concurrently, TKI treatment indirectly caused phenotypic changes in cancer-associated fibroblasts (CAFs), which could contribute to DTP survival through STAT3 activation.

      Materials and Methods

      DTP Modeling

      Animal studies were conducted using approved protocols (AUP.743). Lung cancer PDX models included PHLC137 with EGFR exon 19 delE746_A750 mutation and PHLC164 with EGFR exon 19 deletion plus T790M mutation.
      • Stewart E.L.
      • Mascaux C.
      • Pham N.A.
      • et al.
      Clinical utility of patient-derived xenografts to determine biomarkers of prognosis and map resistance pathways in EGFR-mutant lung adenocarcinoma.
      PHLC137 and PHLC164 were used at passages 12 to 15 and 10, respectively. Replicate NOD SCID mice (NOD.Cg-Prkdcscid/J) seeded with donor tumor fragments were randomized into baseline (BL) or DTP groups at average tumor volumes of 200 to 300 mm3. DTP PDXs were treated daily with 50 mg/kg erlotinib or 25 mg/kg osimertinib (UHN Shanghai Research & Development Inc., Toronto, Ontario, Canada) by oral gavage (excluding weekends) for 30 days. Formalin-fixed, paraffin-embedded tumors were stained using the BenchMark XT autostainer (Ventana Medical System, Tucson, AZ) (Supplementary Table 1). Stained slides were scanned with the Leica Aperio Scanscope AT2 (20×, 0.5 μm/pixel) and images processed with QuPath
      • Bankhead P.
      • Loughrey M.B.
      • Fernández J.A.
      • et al.
      QuPath: open source software for digital pathology image analysis.
      or HALLO software (Indica Labs, Albuquerque, NM) for automated counting and H-scoring.

      Genomics and RNA Expression Analyses

      The PHLC137 BL and DTP tumors were snap frozen and stored at −80°C. DNA and RNA were co-extracted using the Qiagen AllPrep DNA/RNA Micro kit. For whole-exome sequencing, 100 ng of genomic tumor DNA was used along with matching patient normal lung. Somatic mutation and copy number variation detection are detailed in the Supplementary Methods where PhyloWGS was used to reconstruct subclonal tumor composition on the basis of mutation calling.
      • Deshwar A.G.
      • Vembu S.
      • Yung C.K.
      • Jang G.H.
      • Stein L.
      • Morris Q.
      PhyloWGS: reconstructing subclonal composition and evolution from whole-genome sequencing of tumors.
      Microarray transcriptome profiling was performed on four BL and four DTP tumors using Illumina Human HT-12-V4 BeadChips and the cDNA-mediated annealing, selection, extension, and ligation assay (DASL) kit (Illumina) (Supplementary Methods, data have been uploaded on GEO Omnibus site with accession number GSE198672).

      Single-Cell RNA Sequencing, Pseudotime, and Trajectory Analysis

      The BL and DTP tumors were digested with liberase (Sigma Aldrich, St. Louis, MO). Single-cell libraries were prepared using 10× Genomics Single-Cell 3′ Reagent Kits version 2. Volumes for each sample were calculated for target capture of 1000 (DTP) or 6000 (BL) cells. The cDNA was amplified for 12 or 14 cycles (6000 or 1000 cells, respectively) and sequenced (HiSeq 2500). Cell Ranger (version 2.0, 10× Genomics) was used to process barcodes, remove empty “cells,” align reads to GRCh38-mm10 genomes, and generate gene-cell matrices after their standard pipeline.
      • Macosko E.Z.
      • Basu A.
      • Satija R.
      • et al.
      Highly parallel genome-wide expression profiling of individual cells using nanoliter droplets.
      Seurat (version 2.3.4) was used for further quality control, cell filtering, normalization, batch correction, and clustering. Clusters were visualized using t-distributed Stochastic Neighbor Embedding plots (Supplementary Methods, data have been uploaded on GEO Omnibus site with accession number GSE199716). Pseudotime and trajectory analysis were performed using Monocle (version 2.14.0)
      • Trapnell C.
      • Cacchiarelli D.
      • Grimsby J.
      • et al.
      The dynamics and regulators of cell fate decisions are revealed by pseudotemporal ordering of single cells.
      (Supplementary Methods).

      Identification of Differentially Expressed Genes Between DTP and BL Tumors, Gene Set Enrichment Analysis

      Differentially expressed genes (DEGs) were identified by the “limma” package (version 3.34.4) (microarray) or MAST in Seurat (version 2.3.4)
      • Finak G.
      • McDavid A.
      • Yajima M.
      • et al.
      MAST: a flexible statistical framework for assessing transcriptional changes and characterizing heterogeneity in single-cell RNA sequencing data.
      and used for Gene Set Enrichment Analysis (GSEA).
      • Subramanian A.
      • Tamayo P.
      • Mootha V.K.
      • et al.
      Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles.
      The GSEA was performed with the GSEA (Broad Institute, version: 4.0.3) or Enrichr tools,
      • Kuleshov M.V.
      • Jones M.R.
      • Rouillard A.D.
      • et al.
      Enrichr: a comprehensive gene set enrichment analysis web server 2016 update.
      using MSigDB version 7 (Supplementary Methods).
      • Reich M.
      • Liefeld T.
      • Gould J.
      • Lerner J.
      • Tamayo P.
      • Mesirov J.P.
      GenePattern 2.0.

      Signature Scoring in Cancer Cells and Mouse CAFs

      Cellular functions and biological processes were more robustly characterized by cell type and gene/pathway-specific transcriptional gene sets (signatures) (Supplementary Methods).

      Establishment of CAF Cultures

      Mouse and human CAFs were established from PHLC137 or human lung cancer patient tumors, respectively (Supplementary Methods).

      Results

      DTP Generation in an EGFR-Mutant PDX Model

      Erlotinib shrinks PHLC137 PDX tumors, which after 1 month of treatment enter a minimal RD (MRD) state (Fig. 1A).
      • Stewart E.L.
      • Mascaux C.
      • Pham N.A.
      • et al.
      Clinical utility of patient-derived xenografts to determine biomarkers of prognosis and map resistance pathways in EGFR-mutant lung adenocarcinoma.
      As confirmation that MRD contains DTPs that have not yet progressed to acquired resistance, tumors regrew when TKI treatment ceased and re-responded to drug rechallenge.
      • Stewart E.L.
      • Mascaux C.
      • Pham N.A.
      • et al.
      Clinical utility of patient-derived xenografts to determine biomarkers of prognosis and map resistance pathways in EGFR-mutant lung adenocarcinoma.
      BL and erlotinib-treated DTP tumors had mucinous differentiation, with greater acinar formation in DTP tumors (Supplementary Fig. 1), indicating DTPs did not dedifferentiate. In addition, the DTP tumors exhibited increased stromal fibrosis (Fig. 1B), resembling post-neoadjuvant EGFR TKI-treated lung cancers with EGFR mutations,
      • Lara-Guerra H.
      • Chung C.T.
      • Schwock J.
      • et al.
      Histopathological and immunohistochemical features associated with clinical response to neoadjuvant gefitinib therapy in early stage non-small cell lung cancer.
      and had reduced Ki67 expression (Fig. 1C and D).
      Figure thumbnail gr1
      Figure 1Generation of erlotinib-induced DTP cells in an EGFR-mutant primary PDX model. (A) Sensitivity of the PHLC137 PDX model to erlotinib. Red circle illustrates when BL samples were collected (before treatment). Blue circle illustrates when DTP samples were collected. Arrows indicate when drug treatment was first started in the erlotinib-treated arms and when it was stopped to reveal that DTPs are still viable. (B) H&E histologic views of BL and DTP tumors. (C) Immunohistologic comparison of expression of the Ki67 proliferation marker between BL and DTP tumors. (D) Quantification of Ki67 immunohistologic staining of BL and DTP tumors (n = 5) using the Cytonuclear Algorithm in HALLO. p value was calculated using an unpaired t test. Scale bars in (B) and (C) denote 100 μm. BL, baseline; DTP, drug-tolerant persister; H&E, hematoxylin and eosin; PDX, patient-derived NSCLC xenograft.

      DTPs Are Not Genetically Distinct From BL Tumor Cells, But Are Transcriptomically Different

      The DTP and BL tumors did not differ in exomic mutations (Supplementary Table 2), copy number alterations (Supplementary Table 3), or subclonal architecture (Fig. 2A and B). Thus, after erlotinib treatment, a genetically distinct drug-tolerant subclone did not emerge. To explore transcriptomic differences, we used microarrays to identify 1051 DEGs that differed by at least twofold (539 up-regulated genes, 512 down-regulated genes) in DTPs relative to BL cells. (Fig. 2C, Supplementary Table 4). The GSEA using the Hallmark gene sets found that down-regulated DTP DEGs were enriched for cell cycle, E2F, and MYC activities (Fig. 2D and Supplementary Table 5), consistent with reduced DTP Ki67 expression (Fig. 1C and D). Up-regulated DTP DEGs were most enriched for the NF-κB pathway activation (Fig. 2D and Supplementary Table 5). At the gene level (Supplementary Table 4), the most up-regulated DTP gene was APOE (Fig. 2C). The APOE is up-regulated during brain and lung injury to protect against oxidative damage,
      • Laskowitz D.T.
      • Horsburgh K.
      • Roses A.D.
      Apolipoprotein E and the CNS response to injury.
      ,
      • Yao X.
      • Gordon E.M.
      • Figueroa D.M.
      • Barochia A.V.
      • Levine S.J.
      Emerging roles of apolipoprotein E and apolipoprotein A-I in the pathogenesis and treatment of lung disease.
      which is triggered by EGFR TKIs.
      • Oren Y.
      • Tsabar M.
      • Cuoco M.S.
      • et al.
      Cycling cancer persister cells arise from lineages with distinct programs.
      ,
      • Raha D.
      • Wilson T.R.
      • Peng J.
      • et al.
      The cancer stem cell marker aldehyde dehydrogenase is required to maintain a drug-tolerant tumor cell subpopulation.
      MUC5AC also was up-regulated, along with several other mucins (MUC4, MUC20, MUC1), supporting DTPs maintaining or increasing mucinous lineage differentiation. In addition, several alcohol dehydrogenases (ALDHs), ALDH1A1 and ALDH1A3, were expressed at higher levels, which may help DTPs tolerate TKI-generated reactive oxygen species (ROS).
      • Raha D.
      • Wilson T.R.
      • Peng J.
      • et al.
      The cancer stem cell marker aldehyde dehydrogenase is required to maintain a drug-tolerant tumor cell subpopulation.
      Figure thumbnail gr2
      Figure 2BL and DTP tumors from PHLC137 are genetically similar but transcriptomically different. (A) Totality of clones across all BL and DTP samples, as predicted by PhyloWGS. Clones are indicated by numbered nodes, with branch points starting at node 1 illustrating progressive restriction of genetic alterations to fewer subsets of cells. (B) Relative prevalence of each clone node from (A) in each tumor sample. (C) Volcano plot of DEGs between DTP and BL tumor cells identified from microarray data. (D) GSEA of up- and down-regulated microarray DEGs (DTP relative to BL) using the Hallmark gene sets. Top significant results plotted by their NES and ranked by q-value are illustrated. (E) Immunohistologic comparison of expression of the candidate DTP markers, ECM1 and AKAP12, between BL and DTP tumors. Scale bars are 100 μm. (F) Quantification of ECM1 and AKAP12 immunohistologic staining of BL and DTP tumors. H-score (staining intensity × percentage of cells) and percentage of positive cells were quantified using QuPath software. p values were calculated using unpaired t tests. BL, baseline; DEG, differentially expressed gene; DTP, drug-tolerant persister; FC, fold change; FDR, false discovery rate; GSEA, Gene Set Enrichment Analysis; NES, normalized enrichment scores; NS, not significant.
      To determine whether DEGs translated into protein expression differences that could mark the DTP state, we performed IHC on two up-regulated DTP genes, ECM1 and AKAP12. Both proteins were more highly expressed in DTP versus BL tumor cells (Fig. 2E and F).

      A DTPL Subpopulation Preexists in Treatment-Naive Tumors

      Although we profiled approximately 3000 high-quality BL cancer cells by single-cell RNA sequencing (scRNA-seq), only 94 high-quality cells were obtained from the microscopic DTP tumors, consistent with their MRD nature. Combining BL and DTP data revealed seven clusters (Fig. 3A and B. Cluster 6 included 77% (72 of 94 cells) of DTPs and 4% (122 of 2989 cells) of BL cells. Because this cluster contained most DTPs, we tentatively denoted the 122 BL cells in cluster 6 as DTPL. In the pseudotime trajectory analysis, tumor cells organized along a tripartite trajectory (Fig. 3C, top panel). All DTPs were found exclusively at the terminal end of a single state, suggesting a stable state. Most DTPs and DTPL cells clustered closely together at the terminal end of the same state, whereas some DTPL cells were located on the backbone path, suggesting a transition state. To determine the direction of this transition, we plotted the pseudotime along the trajectory (Fig. 3C, lower panel). DTP or DTPL cells had the highest pseudotime value, whereas BL cells had the lowest values, suggesting the transition progressed from BL to DTPL to DTP. This analysis supports DTPL cells not only being similar to DTPs but also potentially being precursors to DTPs (Fig. 3C).
      Figure thumbnail gr3
      Figure 3scRNA-seq analysis of BL and DTP tumors reveals a preexisting DTPL subpopulation in treatment-naive PHLC137 tumors. (A) tSNE plots revealing clustering of individual BL (coral) and DTP (teal) cells. (B) Seurat assignment of BL and DTP cells to distinct clusters. Cells within cluster 6 are outlined in plots found in both (A) and (B). (C) Trajectory-pseudotime analysis of the potential relatedness between BL, DTP, and DTPL cells. Top panel, 2D plane projection depicting divergent trajectories among BL (coral), DTP (green), and DTPL (blue) cells using Monocle. Lower panel, pseudotime analysis overlaid on 2D plane projection depicting divergent trajectories among BL, DTP, and DTPL cells. The darkness of the blue color depicts the pseudotime (dark is early and light is late). In both panels, the solid black line indicates the main diameter path. Pseudotime trajectories indicate DTP (light blue) as a defined final stable state that DTPL cells can transition to at a low frequency. (D) Fraction of BL, DTP, and DTPL cells in different cell cycle phases, separated by cluster 6 versus noncluster 6. (E) Venn diagrams of overlaps between DTP and DTPL DEGs (relative to BL) on the basis of scRNA-seq data. (F) Quantification of expression of signatures associated with DTPs or RD from other studies in PHLC137 BL, DTP, and DTPL cells. See for signature origins and descriptions. 2D, two-dimensional; BL, baseline; DTP, drug-tolerant persister; DTPL, DTP-like; FAM, fatty acid metabolism; RD, residual disease; ROS, reactive oxygen species; scRNA-seq, single-cell RNA sequencing; tSNE, t-distributed Stochastic Neighbor Embedding.
      The cycling properties of BL and DTP cells were explored by applying cell cycle signatures of G1, S, and G2/M phases to the scRNA-seq data.
      • Tirosh I.
      • Izar B.
      • Prakadan S.M.
      • et al.
      Dissecting the multicellular ecosystem of metastatic melanoma by single-cell RNA-seq.
      Within cluster 6, 99% of the DTPs and 97% of the DTPL cells were in G1, whereas outside of cluster 6, 95% of DTP and only 48% of BL cells were in G1 (Fig. 3D). Thus, DTPL cells share the reduced cycling profile of DTPs. To further investigate similarity between DTPs and DTPL cells, DEGs relative to BL cells were identified from scRNA-seq data; 1024 DEGs were identified between DTP cells in cluster 6 and BL cells outside of cluster 6 (310 down-regulated, 714 up-regulated in DTPs) (Supplementary Fig. 2, Supplementary Table 6). For the DTPL cells, 524 DEGs (238 down-regulated, 286 up-regulated) were identified relative to BL cells outside of cluster 6 (Supplementary Table 7); 70% of DTPL DEGs (both directions) were also DTP DEGs (up-regulated: [203 of 286 genes], down-regulated [167 of 238 genes]) (Fig. 3E), further supporting their similarity.
      We next explored gene signatures elevated in DTPs from other lung cancer models and patients. The FAM and ROS signatures from the Hallmark gene sets were reported to be higher in DTPs generated after osimertinib treatment of PC9 cells
      • Oren Y.
      • Tsabar M.
      • Cuoco M.S.
      • et al.
      Cycling cancer persister cells arise from lineages with distinct programs.
      ; both signature scores were lowest in our BL cells outside cluster 6, intermediary for DTPL cells, and highest in DTPs (Fig. 3F). Using pooled tumor cell cohorts bearing EGFR and other driver mutations that were treated with their respective kinase inhibitors, an alveolar signature and a set of 629 genes that are differentially overexpressed in RD compared with treatment-naive tumors (RD signature) have previously been reported.
      • Maynard A.
      • McCoach C.E.
      • Rotow J.K.
      • et al.
      Therapy-induced evolution of human lung cancer revealed by single-cell RNA sequencing.
      These genes were also expressed lowest in our BL cells, intermediary in DTPL cells, and highest in DTPs (Fig. 3F). Overall, DTPs and DTPL cells are similarly enriched for the G1 cell cycle phase, share multiple gene and signature expression profiles reported as being overexpressed in other DTPs or post-treatment residual tumors, and are suggested to be developmentally related by trajectory and pseudotime analysis.
      By scRNA-seq, the top DEG enriched in DTPs was ECM1, which was also found in 13% of the BL tumor cells (Supplementary Table 6) and confirmed to be expressed at the protein level in a small subset of BL tumor cells (Fig. 2E and F). To investigate the possibility that ECM1 could also mark preexisting DTPL cells in human EGFR-mutant lung tumors, we evaluated ECM1 expression by IHC in treatment-naive resected primary lung adenocarcinomas (LUADs) with EGFR mutations. Of all cases, 71% had some ECM1 expression, with most of these tumors (79%) having rare or focal expression (Supplementary Fig. 3), consistent with ECM1 marking a minor population of cancer cells in most EGFR-mutant lung cancers.

      Top Pathway Changes That Characterize DTPs and DTPL Cells

      The GSEA was used to identify biological pathways that characterize DTP and DTPL populations. For DTPs, 24.4% (256 of 1051) of the microarray and 25% (256 of 1024) of the scRNA-seq DEGs overlapped with the other platform (Supplementary Fig. 4). The 256 common DEGs included ECM1 and AKAP12, cell cycle, E2F, MYC, and mucin (MUC5AC, MUC4) genes (Supplementary Fig. 2 and Supplementary Table 8). Using these common “high-confidence” DTP DEGs, GSEA with the Hallmark gene sets still identified NF-κB as the top up-regulated pathway, whereas E2F, cell cycle, and MYC remained the top down-regulated pathways in DTPs (Fig. 4A and Supplementary Table 9). The same pathways were the top hits for DTPL DEGs (Fig. 4B and Supplementary Table 10). Because NF-κB is activated by EGFR TKI treatment and is a survival factor,
      • Blakely C.M.
      • Pazarentzos E.
      • Olivas V.
      • et al.
      NF-κB-activating complex engaged in response to EGFR oncogene inhibition drives tumor cell survival and residual disease in lung cancer.
      we corroborated the NF-κB finding by quantifying expression of another NF-κB signature consisting of direct targets in EGFR-mutant lung cancer cells.
      • Blakely C.M.
      • Pazarentzos E.
      • Olivas V.
      • et al.
      NF-κB-activating complex engaged in response to EGFR oncogene inhibition drives tumor cell survival and residual disease in lung cancer.
      This signature was highest in DTPs, intermediary in DTPL cells, and lowest in BL cells (Fig. 4C). We also found that this NF-κB signature and the Hallmark NF-κB gene set were more highly expressed in PC9 osimertinib-induced DTPs relative to treatment-naive cells (Supplementary Fig. 5), suggesting NF-κB activation is common to EGFR TKI-induced DTPs.
      Figure thumbnail gr4
      Figure 4Analysis of biological pathway changes in DTP and DTPL cells from PHLC137. (A, B, D, E) GSEA of up- and down-regulated DEGs in DTP and DTPL cells (relative to BL). GSEA was performed with either the overlapped DEGs from microarray and scRNA-seq data (DTP) or scRNA-seq data alone (DTPL). The top hits based on adjusted p values are found with the mined database for each plot indicated along the y axis. (C, F) Quantification of expression of signatures associated with NF-κB (C) or NRF2 (F) activity in BL, DTP, and DTPL cells from PHLC137. See for signature details and . Because of the low number of up-regulated DEGs that overlap with the TGF-β gene set, and the inclusion of both induced and repressed targets, including positive and negative regulators of TGF-β signaling within this overlapping group, we have low confidence that TGF-β signaling is increased in DTP or DTPL cells. BL, baseline; DEG, differentially expressed gene; DTP, drug-tolerant persister; DTPL, DTP-like; GSEA, Gene Set Enrichment Analysis; scRNA-seq, single-cell RNA sequencing; TGF-β, transforming growth factor-β.
      We next performed GSEA with the ENCODE and ChEA transcription factor target gene sets (Supplementary Tables 11 and 12). Of the factors whose activity changed in only one direction, GATA2 was the top up-regulated transcription factor in both DTP and DTPL cells (Fig. 4D and E). GATA2 has been reported to promote lung cancer survival at least partly through NF-κB activation,
      • Kumar M.S.
      • Hancock D.C.
      • Molina-Arcas M.
      • et al.
      The GATA2 transcriptional network is requisite for RAS oncogene-driven non-small cell lung cancer.
      although GATA2’s role in survival has been challenged.
      • Tessema M.
      • Yingling C.M.
      • Snider A.M.
      • et al.
      GATA2 is epigenetically repressed in human and mouse lung tumors and is not requisite for survival of KRAS mutant lung cancer.
      MYC was the top down-regulated transcription factor in DTP and DTPL cells, with E2F also significantly down-regulated in both populations (Fig. 4D and E). We also found that NFE2L2/NRF2 activity was predicted to be up-regulated in DTP/DTPL cells, which promotes generation of cycling PC9 DTPs and EGFR TKI resistance.
      • Oren Y.
      • Tsabar M.
      • Cuoco M.S.
      • et al.
      Cycling cancer persister cells arise from lineages with distinct programs.
      ,
      • Foggetti G.
      • Li C.
      • Cai H.
      • et al.
      Genetic determinants of EGFR-driven lung cancer growth and therapeutic response in vivo.
      The NRF2 result was corroborated by our finding that a signature consisting of direct NRF2 targets in A549 cells
      • Namani A.
      • Liu K.
      • Wang S.
      • et al.
      Genome-wide global identification of NRF2 binding sites in A549 non-small cell lung cancer cells by ChIP-Seq reveals NRF2 regulation of genes involved in focal adhesion pathways.
      was highest in DTPs, intermediary in DTPL cells, and lowest in BL cells (Fig. 4F).
      FOXA1, IGF1R, β-catenin, ERK, AXL, and YAP have also been implicated in EGFR TKI resistance.
      • Kurppa K.J.
      • Liu Y.
      • To C.
      • et al.
      Treatment-induced tumor dormancy through YAP-mediated transcriptional reprogramming of the apoptotic pathway.
      ,
      • Maynard A.
      • McCoach C.E.
      • Rotow J.K.
      • et al.
      Therapy-induced evolution of human lung cancer revealed by single-cell RNA sequencing.
      ,
      • Del Bubba M.
      • Di Serio C.
      • Renai L.
      • et al.
      Vaccinium myrtillus L. extract and its native polyphenol-recombined mixture have anti-proliferative and pro-apoptotic effects on human prostate cancer cell lines.
      • Phuchareon J.
      • McCormick F.
      • Eisele D.W.
      • Tetsu O.
      EGFR inhibition evokes innate drug resistance in lung cancer cells by preventing Akt activity and thus inactivating Ets-1 function.
      • Taniguchi H.
      • Yamada T.
      • Wang R.
      • et al.
      AXL confers intrinsic resistance to osimertinib and advances the emergence of tolerant cells.
      • Wang R.
      • Yamada T.
      • Kita K.
      • et al.
      Transient IGF-1R inhibition combined with osimertinib eradicates AXL-low expressing EGFR mutated lung cancer.
      Nevertheless, we did not detect up-regulation of FOXA1 or IGF1R transcripts in DTPs (Supplementary Tables 4 and 6), which is the reported mode of activation.
      • Wang R.
      • Yamada T.
      • Kita K.
      • et al.
      Transient IGF-1R inhibition combined with osimertinib eradicates AXL-low expressing EGFR mutated lung cancer.
      Furthermore, we did not detect increases in β-catenin nuclear localization, ERK phosphorylation, or AXL or YAP protein levels (Supplementary Figs. 6–8). In fact, AXL expression was reduced in DTP tumors. We were unable to reliably detect phosphorylation of AURKA or FRS2, a downstream signaling effector of FGFRs, by IHC in our PHLC137 BL and DTP tumor sections. Therefore, we do not have a conclusive answer as to whether these pathways are activated. Thus, there are commonalities and differences between potential drug tolerance mechanisms that are activated in the PDX versus cell line models.
      • Suda K.
      • Mitsudomi T.
      Drug tolerance to EGFR tyrosine kinase inhibitors in lung cancers with EGFR mutations.

      Similarities in DTP Generation Between Different EGFR-Mutant Lung Cancer PDX Models

      To investigate potential similarities in the heterogeneity of treatment-naive tumor cell populations and DTPs between different EGFR-mutant lung cancer PDXs, we also conducted scRNA-seq analysis on another TKI-treated PDX model. PHLC164 has EGFR exon 19 deletion plus T790M mutations
      • Stewart E.L.
      • Mascaux C.
      • Pham N.A.
      • et al.
      Clinical utility of patient-derived xenografts to determine biomarkers of prognosis and map resistance pathways in EGFR-mutant lung adenocarcinoma.
      and was treated for 50 days with osimertinib to generate DTPs (Fig. 5A). The combined scRNA-seq analysis of BL (n = 3608) and DTP (n = 172) cells revealed that 49% (84 cells) of DTPs are found in cluster 3, which is the cluster that harbors the most DTPs (Fig. 5B). Cluster 3 also contained a minor subpopulation (4.9%, 177 cells) of BL cells (Fig. 5B), putatively representing DTPL cells. Although the DTPL DEG analysis was underpowered, we nevertheless determined that at least 24% (245 of 1034) and 18% (34 of 188) of the DEGs identified in PHLC164 DTPs and DTPL cells, respectively, were shared with the corresponding populations in PHLC137 (Fig. 5C). Overall, the PHLC164 DTP DEGs were also enriched for similar up- and down-regulated pathways as the PHLC137 DTP DEGs (Fig. 5D). Finally, as with PHLC DTPs, PHLC164 DTPs had elevated expression of the RD and alveolar signatures identified in patient samples, and the PHLC164 DTPL cells could be confirmed to have intermediate expression of the alveolar signature (Fig. 5E). Thus, different EGFR-mutant lung cancer PDX models treated with distinct TKIs have some similarities in preexisting DTPL cells and DTPs.
      Figure thumbnail gr5
      Figure 5Analysis of BL and DTP tumor cells in PHLC164 PDXs treated with osimertinib. (A) Plots of tumor growth without drug treatment (control) and with osimertinib treatment initiated at “day 0.” At day 50, DTP tumors were collected from some animals in the treatment arms whereas other mice in that arm were taken off drug to reveal reversibility of the DTP state. (B) scRNA-seq analysis of BL and osimertinib-induced DTPs. Left panel, tSNE plot revealing clusters in the combined BL and DTP cell populations. Right panel, annotations of BL (red) and DTP (blue) cell populations, focusing on cluster 3, which is the dominant DTP cluster. (C) Comparisons of DTP and DTPL DEGs (both relative to BL cells) between PHLC137 and PHLC164. Upper Venn diagrams reveal 1034 differentially expressed DTP genes (DEGs) in PHLC164, with 24% (245 of 1034) of total DTP DEGs overlapping with those in PHLC137 DTPs. Lower Venn diagram illustrates 188 DEGs for DTPL cells in PHLC164, with 18% (34 of 188) overlapping with those of PHLC137 DTPL cells. (D) GSEA analysis of PHLC164 DTP DEGs revealing up-regulated (estrogen response, UV response up, NF-κB, and TP53) and down-regulated (E2F, MYC) pathways. The small number of DEGs identified in PHLC164 DTPL cells limited the effectiveness of GSEA pathway analysis for these genes. (E) Analysis of expression of two patient RD signatures (RD and alveolar)
      • Maynard A.
      • McCoach C.E.
      • Rotow J.K.
      • et al.
      Therapy-induced evolution of human lung cancer revealed by single-cell RNA sequencing.
      in PHLC164 DTP and DTPL cells. BL, baseline; DEG, differentially expressed gene; DTP, drug-tolerant persister; DTPL, DTP-like; FDR, false discovery rate; GSEA, Gene Set Enrichment Analysis; PDX, patient-derived xenograft; RD, residual disease; scRNA-seq, single-cell RNA sequencing; UV, ultraviolet.

      Unique Carcinoma-Associated Fibroblast Phenotype in DTP Tumor Stroma

      Carcinoma-associated fibroblasts (CAFs) are heterogeneous stromal cells that can promote tumor progression across cancers
      • Kalluri R.
      • Zeisberg M.
      Fibroblasts in cancer.
      and have been implicated in NSCLC recurrence post-TKI treatment.
      • Hu H.
      • Piotrowska Z.
      • Hare P.J.
      • et al.
      Three subtypes of lung cancer fibroblasts define distinct therapeutic paradigms.
      How TKI treatment changes CAF identity in treatment naive versus DTP cancer cells is poorly understood. We used our scRNA-seq data (8023 mouse cells: 4478 in BLs, 3545 in DTPs) to compare CAFs in BL versus DTP tumors of PHLC137 model (Fig. 6A). Approximately 3900 mouse cells across both tumor states were annotated as fibroblasts (Fig. 6B). Collagen gene expression was highest in these annotated cells,
      • Tsukui T.
      • Sun K.H.
      • Wetter J.B.
      • et al.
      Collagen-producing lung cell atlas identifies multiple subsets with distinct localization and relevance to fibrosis.
      confirming fibroblast identity (Supplementary Fig. 9). CAFs were clustered into four subpopulations (Fig. 6C), with cluster 4 and cluster 5 being shared at similar percentages between states (approximately 20% each). Cluster 2 was the most abundant CAF subpopulation in BL (53%) but disappeared almost completely in DTP tumors (0.6%). Conversely, cluster 0 was a minor subpopulation in BL (7%) and the dominant subpopulation in the DTP state (59%).
      Figure thumbnail gr6
      Figure 6Erlotinib treatment changes CAF phenotypes in PHLC137 and triggers emergence of a novel CAF population. CAFs were analyzed from the murine scRNA visualization of all murine stromal cells. (B) tSNE plot revealing cell type annotation of the stromal cell fibroblasts (outlined in [A] and [B]). (C) Seurat assignment of all murine cells to distinct clusters. Clusters 0, 2, 4, and 5 consist of CAFs, with their proportions in BL and DTP tumors indicated. (D) Evaluation of the similarity of each of the four CAF clusters to 57 published CAF transcriptomic signatures (identified by last name of first author of the study and the type of CAF or marker described in the study). Enrichment of each signature in each cluster was quantified by scoring expression of all signature genes in all cells for the cluster and by quantifying the separation between top and second-highest scoring cluster. See for details. (E) Trajectory-pseudotime analysis of the potential relatedness between the four CAF clusters. Left and right panels, 2D plane projections depicting divergent trajectories among the four CAF subsets (right panel) in BL and DTP tumors (left panel) using Monocle. Middle panel, pseudotime analysis overlaid on 2D plane projection depicting divergent trajectories among the four CAF subsets. Darkness of the blue color depicts the pseudotime (dark is early and light is late). In all panels, the solid black line indicates the main diameter path. (F) Secreted factors that are candidates for paracrine signals (categorized by function) that were found among cluster-enriched up-regulated DEGs. 2D, two-dimensional; BL, baseline; CAF, cancer-associated fibroblast; DEG, differentially expressed gene; DTP, drug-tolerant persister; scRNA, single-cell RNA; tSNE, t-distributed Stochastic Neighbor Embedding.
      We next evaluated the similarity of the CAF subpopulations to CAF subtypes previously described with transcriptomic signatures in NSCLC and pancreatic and breast cancers (Supplementary Tables 13 and 14). By comparing enrichment scores (Supplementary Table 15) of each signature to each CAF cluster, we could estimate how well any particular signature matched a particular cluster (Fig. 5D and Supplementary Methods). Cluster 2 CAFs were overwhelmingly connected to myofibroblasts “myCAFs” (8/10 myCAF signatures), transforming growth factor (TGF)-β CAFs, or CAFs with increased ECM gene expression. The top signature for cluster 2 CAFs describes a subset of mouse CAFs with a high capacity for ECM remodeling, collagen secretion, and ACTA2 (α-smooth muscle actin or SMA) expression.
      • Dominguez C.X.
      • Müller S.
      • Keerthivasan S.
      • et al.
      Single-cell RNA sequencing reveals stromal evolution into LRRC15+ myofibroblasts as a determinant of patient response to cancer immunotherapy.
      These features are common to myofibroblast-like CAFs, also known as “myCAFs,” which are the most consistently described CAF population in treatment-naive cancers and are induced by TGF-β signaling.
      • Biffi G.
      • Oni T.E.
      • Spielman B.
      • et al.
      IL1-induced JAK/STAT signaling is antagonized by TGFβ to shape CAF heterogeneity in pancreatic ductal adenocarcinoma.
      ,
      • Sperb N.
      • Tsesmelis M.
      • Wirth T.
      Crosstalk between tumor and stromal cells in pancreatic ductal adenocarcinoma.
      The top signature for cluster 4 CAFs describes myCAFs as active in wound healing and T-cell infiltration.
      • Kieffer Y.
      • Hocine H.R.
      • Gentric G.
      • et al.
      Single-cell analysis reveals fibroblast clusters linked to immunotherapy resistance in cancer.
      For cluster 5, several enriched signatures included the term “iCAF,” which was initially coined to represent “non-myCAFs” with a more proinflammatory phenotype owing to cytokine secretion.
      • Öhlund D.
      • Handly-Santana A.
      • Biffi G.
      • et al.
      Distinct populations of inflammatory fibroblasts and myofibroblasts in pancreatic cancer.
      Another signature enriched in this cluster was one for fibroblasts found in multiple tissues across different steady-state and disease conditions.
      • Buechler M.B.
      • Pradhan R.N.
      • Krishnamurty A.T.
      • et al.
      Cross-tissue organization of the fibroblast lineage.
      Only one signature identified cluster 0 as the top hit, but this was a low-confidence descriptor.
      Through pseudotime analysis (Fig. 6E), cluster 2, cluster 5, and cluster 0 CAFs were found as distinct loci at the end of trajectories, suggesting they represent well-defined solid or stable states. The states in cluster 2 and cluster 5 could relate to the CAF identities described by the top signatures for these clusters. Cluster 4 CAFs were dispersed along all trajectories, suggesting they may comprise multiple distinct intermediate activation states.

      CAFs in Treatment-Naive Versus Post-Treatment Stroma Have Enrichment of Distinct Biological Functions

      To further characterize CAF subtypes, we identified DEGs between DTP and BL mouse fibroblast cells (Supplementary Table 16) and DEGs enriched in each cluster (Supplementary Table 17) by scRNA-seq. Using IHC for ACTA2/α-SMA and ICAM-1, DEGs for cluster 2 and cluster 0 CAFs, respectively, we confirmed that some DEGs were differentially expressed at the protein level in BL and DTP tumors (Supplementary Fig. 10). ACTA2 also marks the dominant CAF population in treatment-naive human EGFR-mutant LUADs, which had widespread and intense ACTA2 expression in all samples and complete absence of ICAM-1, the post-treatment cluster 0 marker (Supplementary Fig. 11).
      The GSEA on the top 300 cluster-enriched DEGs identified several biological features in each cluster (Supplementary Table 18). For cluster 2 CAFs, TGF-β signaling was the third most significantly enriched term from the Hallmark collection, whereas extracellular matrix organization, collagen fibril organization, and supramolecular fiber organization were the top three terms from Gene Ontology Biological Processes, all well-known features of myCAFs.
      • Elyada E.
      • Bolisetty M.
      • Laise P.
      • et al.
      Cross-species single-cell analysis of pancreatic ductal adenocarcinoma reveals antigen-presenting cancer-associated fibroblasts.
      Additional myCAF markers, including Acta2, Tagln, Ctgf, and Lrrc15, and various collagen genes and Tnc, were among the top cluster 2 DEGs. For cluster 4, cluster 5, and cluster 0 CAFs, NF-κB activation was found in the Hallmark gene sets, which was supported for cluster 5 and cluster 0 by identification of the NF-κB family member RELA with the ENCODE/ChEA chromatin immunoprecipitation data sets. The Hallmark and ENCODE/ChEA data sets collectively identified IL-6/JAK/STAT3 signaling in cluster 0 CAFs, and Il-6, a target of NF-κB,
      • Liu T.
      • Zhang L.
      • Joo D.
      • Sun S.C.
      NF-κB signaling in inflammation.
      and Stat3, an effector of IL-6 signaling,
      • Johnson D.E.
      • O’Keefe R.A.
      • Grandis J.R.
      Targeting the IL-6/JAK/STAT3 signalling axis in cancer.
      were both up-regulated DEGs in cluster 0 CAFs.
      Analysis of CAF DEGs for secreted signaling factors (Supplementary Table 19) revealed that every CAF cluster expressed growth-promoting and immune-modulatory factors, with cluster 4, cluster 5, and cluster 0 also secreting proangiogenic signals and cluster 2 producing TGF-β ligands (Fig. 6F). Although Fgf7 was up-regulated in cluster 0 CAFs, ERK activation was not increased in DTP cancer cells (Supplementary Fig. 7), making it unclear whether Fgf7 was expressed at sufficient levels to promote EGFR TKI resistance, as has been suggested in studies of human CAFs.
      • Hu H.
      • Piotrowska Z.
      • Hare P.J.
      • et al.
      Three subtypes of lung cancer fibroblasts define distinct therapeutic paradigms.
      Immune-modulatory factors were cluster specific, highlighting that each CAF subtype may exert its own distinct form of immunomodulation.
      In support of IL6-JAK-STAT3 signaling being active in cluster 0 CAFs, we found phospho-STAT3 in approximately 40% of the DTP stromal cells, a twofold increase over BL stroma (Fig. 7A and B). Similarly, STAT3 phosphorylation in DTP cancer cells was increased relative to BL cancer cells (Fig. 7A and B). IL6 has been found to promote survival of EGFR-mutant lung cancer cell lines to TKI treatment.
      • Kim S.M.
      • Kwon O.J.
      • Hong Y.K.
      • et al.
      Activation of IL-6R/JAK1/STAT3 signaling induces de novo resistance to irreversible EGFR inhibitors in non-small cell lung cancer with T790M resistance mutation.
      To verify that PHLC137 cancer cells could be responding to murine Il6 produced by cluster 0 CAFs, we established an in vitro culture of these cancer cells and tested the ability of murine Il6 to activate STAT3 phosphorylation in these cells. Murine Il6 robustly induced STAT3 phosphorylation (phospho-STAT3), which was retained for at least 24 hours (Supplementary Fig. 12). Altogether, these findings suggest that Il6 secreted by cluster 0 CAFs could activate STAT3 in PHLC137 DTPs to promote their survival in vivo.
      Figure thumbnail gr7
      Figure 7Cluster 0 CAFs may promote DTP survival through Il6 and their state is regulated by the balance between TGF-β and NF-κB activating signals. (A) STAT3 is activated in DTP tumor and stromal cells. Immunohistochemical staining of pSTAT3 in PHLC137 BL and DTP tumors. Scale bars are 100 μm. H-score (staining intensity × percentage of cells) and percentage of positive cells were quantified using the Multiplex Algorithm in HALLO. p values were calculated using unpaired t tests. (B) CAF states are plastic and regulated by TGF-β and NF-κb activating signals. Mouse CAFs isolated from PHLC137 BL tumors were treated with TGF-β1 or IL-1α, and expression of CAF cluster-enriched markers was analyzed by RT-qPCR. Marker expression was normalized to β-actin and is relative to untreated control cultures. Data are from two separate experiments, each with three biological replicates. p values were calculated using an unpaired t test. (C) Erlotinib treatment of PHLC137-derived human EGFR-mutant lung cancer cell/mouse CAF cocultures changes the CAF phenotype. Cocultures were treated with 1 μM erlotinib for 12 days, after which expression of the CAF cluster 0 and cluster 2-enriched markers, Il6 and Lrrc15, respectively, was quantified in the murine CAFs by RT-qPCR. Data are the mean ± SD from three biological replicates. p values were calculated using an unpaired t test. (D) Schema of origin of DTPs, changes in CAF phenotypes, and crosstalk between cancer cells and CAFs during EGFR TKI treatment. DTPs are induced by the TKI from a specific subset of DTPL baseline cells. BL cancer cells promote the TGF-β–driven cluster 2 CAF phenotype, which in turn, promotes BL proliferation. As BL cancer cells are ablated (or changed) by the TKI, cluster 2 CAFs convert to cluster 0 CAFs, which promote STAT3 activation in DTP cancer cells by means of IL6. Each CAF cluster secretes its own unique complement of growth and immunomodulatory factors. BL, baseline; CAF, cancer-associated fibroblast; DTP, drug-tolerant persister; DTPL, DTP-like; IL, interleukin; pSTAT3, phosphorylated-STAT3; RT-qPCR, reverse transcriptase-quantitative polymerase chain reaction; S, stromal cells; T, tumor cells; TGF-β, transforming growth factor-β; TKI, tyrosine kinase inhibitor.

      CAF Phenotypes Are Plastic and Regulated by Secreted Factors, Interactions With Cancer Cells, and EGFR TKI Treatment

      The striking transitions between cluster 2 and cluster 0 CAFs during TKI treatment suggest that these populations interconvert. To test this possibility, we first derived CAF cultures from PHLC137 BL and DTP tumors. Under BL conditions, BL CAFs expressed higher levels of cluster 2 CAF markers, whereas DTP CAFs had higher expression of cluster 0 CAF markers (Supplementary Fig. 13), consistent with the in vivo transcriptomic profiling of these CAFS. Because of difficulties in growing DTP CAFs, we could only test the response of BL CAFs to secreted factors. The TGF-β stimulation of BL CAFs further increased expression of cluster 2 CAF markers (Fig. 7B, top panel), consistent with the cluster 2 phenotype having similarity to TGF-β–regulated myCAFs. Conversely, stimulation of BL CAFs with IL1α, an NF-κB activator, induced expression of cluster 0 markers (Fig. 7B, bottom panel). These findings were paralleled by similar experiments with human CAFs that were also immortalized in a cluster 2 myCAF state
      • Navab R.
      • Strumpf D.
      • Bandarchi B.
      • et al.
      Prognostic gene-expression signature of carcinoma-associated fibroblasts in non-small cell lung cancer.
      and therefore less susceptible to increases in the cluster 2 phenotype (Supplementary Fig. 14). Nevertheless, overall, the data support a general CAF plasticity and interconversion mechanism that is regulated by TGF-β and NF-κB activators. Finally, we tested whether in co-cultures of EGFR-mutant lung cancer cells and CAFs, erlotinib treatment could convert cluster 2 CAFs to a cluster 0 phenotype. Indeed, after 12 days of erlotinib treatment, expression of one of the more robust markers for the cluster 0 phenotype, Il6, was increased 12-fold, whereas expression of the cluster 2 marker Lrrc15 was down-regulated (Fig. 7C). Thus, the major CAF populations pre- and post-TKI treatment likely reflect a single type of CAF population whose phenotype is regulated by TGF-β and NF-κB activators and their modulation by the absence or presence of a TKI (Fig. 7D).

      Discussion

      Persistence of drug-tolerant cancer cells after TKI treatment is a major obstacle to achieving better outcomes for patients with EGFR-mutant lung cancer. Here, we investigated DTP origins, their potential heterogeneity and biological underpinnings, and stromal changes that may contribute to DTP survival. To address these questions, we performed a comprehensive in vivo study with erlotinib in an EGFR-mutant lung cancer PDX and used scRNA-seq to characterize individual cancer and stromal cells.
      Most DTP work using PC9 cells collectively supports DTPs being stochastically induced by TKI treatment.
      • Sharma S.V.
      • Lee D.Y.
      • Li B.
      • et al.
      A chromatin-mediated reversible drug-tolerant state in cancer cell subpopulations.
      ,
      • Oren Y.
      • Tsabar M.
      • Cuoco M.S.
      • et al.
      Cycling cancer persister cells arise from lineages with distinct programs.
      ,
      • Guler G.D.
      • Tindell C.A.
      • Pitti R.
      • et al.
      Repression of stress-induced LINE-1 expression protects cancer cell subpopulations from lethal drug exposure.
      Our data from PHLC137 support a DTP origin model that is a hybrid of the preexisting and induced concepts. We found evidence for specific preexisting cancer cells being in an intermediate DTPL state, which is further induced by the TKI to achieve the full DTP state. This intermediate state is identifiable by FAM, ROS, NRF2, and NF-κB signatures, which represent pathways that promote DTP generation in other lung cancer models,
      • Oren Y.
      • Tsabar M.
      • Cuoco M.S.
      • et al.
      Cycling cancer persister cells arise from lineages with distinct programs.
      ,
      • Blakely C.M.
      • Pazarentzos E.
      • Olivas V.
      • et al.
      NF-κB-activating complex engaged in response to EGFR oncogene inhibition drives tumor cell survival and residual disease in lung cancer.
      and by signatures that characterize targeted-treatment RD in patients with NSCLC,
      • Maynard A.
      • McCoach C.E.
      • Rotow J.K.
      • et al.
      Therapy-induced evolution of human lung cancer revealed by single-cell RNA sequencing.
      also supporting relevance of the PDX DTPs to clinical MRD.
      In some cancer models, DTPs have also been found to arise from specific preexisting treatment-naive cancer cells. DTPs generated from MET-amplified gastric cancer cell lines after crizotinib treatment preferentially arise from rare preexisting cells with high aldehyde dehydrogenase activity.
      • Raha D.
      • Wilson T.R.
      • Peng J.
      • et al.
      The cancer stem cell marker aldehyde dehydrogenase is required to maintain a drug-tolerant tumor cell subpopulation.
      Furthermore, a study of vemurafenib resistance in melanoma cell lines proposed that DTPs arise from specific preexisting DTPL populations with an intermediate DTP phenotype, with drug treatment completing induction of the DTP state,
      • Shaffer S.M.
      • Dunagin M.C.
      • Torborg S.R.
      • et al.
      Rare cell variability and drug-induced reprogramming as a mode of cancer drug resistance.
      similar to our model.
      DTPs generated from PHLC137 have biological similarities and differences with DTPs described in other lung cancer studies. PC9 cell DTPs can be separated into two classes depending on their propensity to cycle, with each class having up-regulation of distinct biological pathways.
      • Oren Y.
      • Tsabar M.
      • Cuoco M.S.
      • et al.
      Cycling cancer persister cells arise from lineages with distinct programs.
      Cycling DTPs uniquely have induction of FAM and antioxidant responses by means of NRF2. Nevertheless, PHLC137 DTP and DTPL cells, which are mostly in the G1 phase of the cell cycle, also had signs of up-regulation of these activities. Thus, these pathways may operate more broadly in noncycling DTPs depending on the context. Because the preexisting DTPL cells in PHLC137 had intermediate elevation of FAM, ROS, and NRF2 signatures, these pathways may be partially activated, independent of TKI exposure. PHLC137 DTP and DTPL cells also had signs of increased NF-κB signaling. NF-κB activation by TKIs has been documented in EGFR-mutant lung cancer models, including a PDX established from oligometastatic disease, genetically engineered mice, and cell lines, where it promotes TKI resistance.
      • Blakely C.M.
      • Pazarentzos E.
      • Olivas V.
      • et al.
      NF-κB-activating complex engaged in response to EGFR oncogene inhibition drives tumor cell survival and residual disease in lung cancer.
      The broadness of this activation may stem from direct activation of NF-κB through crosstalk with TKI-engaged EGFRs.
      • Blakely C.M.
      • Pazarentzos E.
      • Olivas V.
      • et al.
      NF-κB-activating complex engaged in response to EGFR oncogene inhibition drives tumor cell survival and residual disease in lung cancer.
      Nevertheless, we did not detect transcriptional changes in IGF-1R expression or AXL protein levels, nor did we find evidence for β-catenin, ERK, or YAP activation, which have been variously reported to be important to EGFR TKI tolerance in other in vitro settings.
      • Suda K.
      • Mitsudomi T.
      Drug tolerance to EGFR tyrosine kinase inhibitors in lung cancers with EGFR mutations.
      Thus, either model or patient-specific differences may also affect mechanisms governing DTP generation. In further support of this possibility, we could confirm some similarities between PHLC137 DTP and DTPL cells with DTP and DTPL cells found in a second EGFR-mutant PDX model treated with a different TKI. These similarities included a preexisting DTPL population in treatment-naive tumors and evidence of increased NF-κB activity and expression of signatures associated with RD in patients.
      Our study is the first to use a NSCLC PDX model to investigate the potential role of CAFs in promoting survival of EGFR TKI-induced DTPs. Clustering of scRNA-seq data from PHLC137 identified four CAF populations, most notably myCAFs (cluster 2) as distinct to treatment-naive tumors, and a novel population as distinct to post-treatment tumors (cluster 0). myCAFs were defined by TGF-β signaling and secretion of growth factors and ECM components, which likely promote BL tumor growth. Post-treatment CAFs were the least myCAF like and were associated with NF-κB and STAT3 activation. Although NF-κB activation was found in cluster 4 and cluster 5 CAFs, the number of DEGs connected to NF-κB signaling in post-treatment CAFs was markedly higher. Post-treatment CAFs also had increased IL1R expression, which contributes to NF-κB activation in other CAF populations,
      • Biffi G.
      • Oni T.E.
      • Spielman B.
      • et al.
      IL1-induced JAK/STAT signaling is antagonized by TGFβ to shape CAF heterogeneity in pancreatic ductal adenocarcinoma.
      Tnfsf9, an NF-κB–activating ligand,
      • Wang C.
      • Lin G.H.
      • McPherson A.J.
      • Watts T.H.
      Immune regulation by 4-1BB and 4-1BBL: complexities and challenges.
      and IL-6, an NF-κB target. Post-treatment CAFs also expressed ICAM-1. Although both IL-6 and ICAM-1 have been described as iCAF markers, no iCAF signature and no signature overall had good similarity to cluster 0 CAFs, indicating the iCAF designation, as described across multiple studies, does not apply to post-treatment CAFs. Rather, these CAFs seem to be specific to the post-drug treatment state, which may represent a distinct type of CAF biology.
      Cluster 4 and cluster 5 CAF proportions did not change during TKI treatment. Cluster 4 CAFs may represent a collection of intermediately activated states, given their positioning along all pseudotime trajectories. Signature analysis identified cluster 5 CAFs as the most iCAF-like by some signatures, although other iCAF signatures were most enriched in other clusters. Furthermore, all four CAF clusters expressed their own specific complement of immunomodulatory molecules among their top DEGs. These observations further highlight the nonspecificity of the iCAF designation. The emergent number of post-treatment CAFs quantitatively mirrored the change in myCAFs, suggesting the former are derived from myCAFs by interconversion. Indeed, treatment of mouse or human CAFs with TGF-β induced the myCAF phenotype, whereas treatment of the same CAFs with IL1α induced cluster 0 marker expression.
      Adaptive changes to TKI treatment in both cancer cells and CAFs likely contribute to the emergence and survival of DTPs (Fig. 6E). In BL, TGF-β–producing tumor cells may condition fibroblasts to adopt a myCAF phenotype, which may be strengthened by TGF-β production from myCAFs. The TKI kills most cancer cells, except for a small number of DTPL cells that are enriched for the G1 phase of the cell cycle and have partial activation of survival pathways. The TKIs further activate pathways such as NF-κB (and FAM and ROS, in some contexts) to complete induction of the DTP state. In DTP tumors, TGF-β signaling in CAFs is reduced, either by elimination of cancer cells or a change in their pro–TGF-β state that may be caused by TKI exposure, allowing other signals to induce new CAF states such as the cluster 0 phenotype. This aspect of the model is supported by our in vitro coculture experiments and by data with pancreatic ductal adenocarcinoma models where TGF-β signaling overrides NF-κB–driven non-myCAF fates.
      • Biffi G.
      • Oni T.E.
      • Spielman B.
      • et al.
      IL1-induced JAK/STAT signaling is antagonized by TGFβ to shape CAF heterogeneity in pancreatic ductal adenocarcinoma.
      Thus, BL cancer cell depletion or changes induced by TKIs could reduce TGF-β signaling enough to allow NF-κB to drive a post-treatment CAF fate that includes increased IL-6 expression, which protects EGFR-mutant lung cancer cells from TKI-induced death by means of STAT3.
      • Kim S.M.
      • Kwon O.J.
      • Hong Y.K.
      • et al.
      Activation of IL-6R/JAK1/STAT3 signaling induces de novo resistance to irreversible EGFR inhibitors in non-small cell lung cancer with T790M resistance mutation.
      ,
      • Nan J.
      • Du Y.
      • Chen X.
      • et al.
      TPCA-1 is a direct dual inhibitor of STAT3 and NF-κB and regresses mutant EGFR-associated human non-small cell lung cancers.
      While recognizing the limitation that our results are derived from one to two PDX models, to our knowledge, this is the first report of an EGFR-mutant DTP study conducted in PDXs. We are encouraged by our findings that many of the phenotypic changes found in these models are consistent with those observed in cell line and post–TKI-treated patient tumors. Nevertheless, we believe that additional studies using more PDX models are warranted to investigate potential mechanistic heterogeneity in DTP generation across EGFR-mutant tumors. This may be challenging given the difficulty of establishing EGFR-mutated PDX models, especially from early stage-resected patient tumors.
      • Stewart E.L.
      • Mascaux C.
      • Pham N.A.
      • et al.
      Clinical utility of patient-derived xenografts to determine biomarkers of prognosis and map resistance pathways in EGFR-mutant lung adenocarcinoma.
      ,
      • Johnson D.E.
      • O’Keefe R.A.
      • Grandis J.R.
      Targeting the IL-6/JAK/STAT3 signalling axis in cancer.
      Similar to cell line models, a further limitation of PDX studies is the incomplete immune microenvironment, whose role in DTP generation cannot be comprehensively assessed.
      In conclusion, we suggest that the DTP and post-treatment CAF pathways identified in this work, especially NF-κB and IL-6/JAK/STAT3, be further investigated as targets for treatments that may reduce recurrence in patients. In this context, Blakely et al.
      • Blakely C.M.
      • Pazarentzos E.
      • Olivas V.
      • et al.
      NF-κB-activating complex engaged in response to EGFR oncogene inhibition drives tumor cell survival and residual disease in lung cancer.
      have reported that PBS-1086, a selective inhibitor targeting both the canonical and noncanonical NF-κB pathways increased the response to erlotinib of tumors formed by EGFR-mutant lung cancer cell lines and an erlotinib-resistant PDX. Furthermore, in the 11-18 EGFRL858R-mutant lung cancer cell line transfected with an IL-6 expression vector to mimic exposure to IL-6 from the microenvironment, co-treatment with the JAK2 inhibitor ruxolitinib lowered the IC50 to erlotinib and increased its effectiveness in killing cancer cells.
      • Blakely C.M.
      • Pazarentzos E.
      • Olivas V.
      • et al.
      NF-κB-activating complex engaged in response to EGFR oncogene inhibition drives tumor cell survival and residual disease in lung cancer.
      Further studies in additional PDX models are warranted to assess the broader importance of these pathways to DTP generation.

      CRediT Authorship Contribution Statement

      Nadeem Moghal: Conceptualization, Data curation and analysis, Funding and supervision, Manuscript writing and editing, Final manuscript approval.
      Quan Li: Conceptualization, Data curation and analysis, Manuscript writing and editing, Final manuscript approval.
      Erin L. Stewart: Conceptualization; Data generation, Curation and analysis, Manuscript writing and editing, Final manuscript approval.
      Roya Navab: Conceptualization, Data generation, Curation and analysis, Manuscript writing and editing, Final manuscript approval.
      Masashi Mikubo: Data generation, Curation and analysis, Final manuscript approval.
      Elisa D’Arcangelo: Data curation and analysis, Manuscript writing and editing, Final manuscript approval.
      Sebastiao N. Martins-Filho: Data generation, Curation and analysis, Manuscript writing and editing, Final manuscript approval.
      Vibha Raghavan: Data curation and analysis, Final manuscript approval.
      Nhu-An Pham: Data curation and analysis, Final manuscript approval.
      Ming Li: Data generation, Curation and analysis, Final manuscript approval.
      Frances A. Shepherd: Manuscript writing and editing, Final manuscript approval.
      Geoffrey Liu: Conceptualization, Data analysis, Funding and supervision, Manuscript writing and editing, Final manuscript approval.
      Ming-Sound Tsao: Conceptualization, Data analysis, Funding and supervision, Manuscript writing and editing, Final manuscript approval.

      Acknowledgments

      This work was supported by the Canadian Institute of Health Research Foundation grant FDN-148395 (MST) and the Canadian Cancer Society Research Institute IMPACT grants 701595 (MST) and 703206 (GL). EL Stewart was supported by the Excellence in Radiation Research for the 21st Century (EIRR21) scholarship and the Frederick Banting and Charles Best Canada Graduate Scholarship Doctoral Award. Drs. S. Martins-Filho and M. Mikubo were supported by the Terry Fox Foundation Training Program in Molecular Pathology of Cancer. Dr. Shepherd is the Scott Taylor Chair in Lung Cancer Research. Dr. Liu is the Alan B. Brown Chair in Molecular Genomics. Dr. Tsao is the M. Qasim Choksi Chair in Lung Cancer Translational Research.

      Supplementary Data

      References

        • Hsu W.H.
        • Yang J.C.
        • Mok T.S.
        • Loong H.H.
        Overview of current systemic management of EGFR-mutant NSCLC.
        Ann Oncol. 2018; 29: i3-i9
        • Greenhalgh J.
        • Boland A.
        • Bates V.
        • et al.
        First-line treatment of advanced epidermal growth factor receptor (EGFR) mutation positive non-squamous non-small cell lung cancer.
        Cochrane Database Syst Rev. 2021; 3: CD010383
        • Chong C.R.
        • Jänne P.A.
        The quest to overcome resistance to EGFR-targeted therapies in cancer.
        Nat Med. 2013; 19: 1389-1400
        • Lim S.M.
        • Syn N.L.
        • Cho B.C.
        • Soo R.A.
        Acquired resistance to EGFR targeted therapy in non-small cell lung cancer: mechanisms and therapeutic strategies.
        Cancer Treat Rev. 2018; 65: 1-10
        • Pallis A.G.
        • Serfass L.
        • Dziadziusko R.
        • et al.
        Targeted therapies in the treatment of advanced/metastatic NSCLC.
        Eur J Cancer. 2009; 45: 2473-2487
        • Ahn M.J.
        • Tsai C.M.
        • Shepherd F.A.
        • et al.
        Osimertinib in patients with T790M mutation-positive, advanced non-small cell lung cancer: long-term follow-up from a pooled analysis of 2 phase 2 studies.
        Cancer. 2019; 125: 892-901
        • He J.
        • Huang Z.
        • Han L.
        • Gong Y.
        • Xie C.
        Mechanisms and management of 3rd-generation EGFR-TKI resistance in advanced nonsmall cell lung cancer (review).
        Int J Oncol. 2021; 59: 90
        • Passaro A.
        • Jänne P.A.
        • Mok T.
        • Peters S.
        Overcoming therapy resistance in EGFR-mutant lung cancer.
        Nat Cancer. 2021; 2: 377-391
        • Su K.Y.
        • Chen H.Y.
        • Li K.C.
        • et al.
        Pretreatment epidermal growth factor receptor (EGFR) T790M mutation predicts shorter EGFR tyrosine kinase inhibitor response duration in patients with non-small-cell lung cancer.
        J Clin Oncol. 2012; 30: 433-440
        • Echeverria G.V.
        • Ge Z.
        • Seth S.
        • et al.
        Resistance to neoadjuvant chemotherapy in triple-negative breast cancer mediated by a reversible drug-tolerant state.
        Sci Transl Med. 2019; 11eaav0936
        • Rehman S.K.
        • Haynes J.
        • Collignon E.
        • et al.
        Colorectal cancer cells enter a diapause-like DTP state to survive chemotherapy.
        Cell. 2021; 184: 226-242.e21
        • Mikubo M.
        • Inoue Y.
        • Liu G.
        • Tsao M.S.
        Mechanism of drug tolerant persister cancer cells: the landscape and clinical implication for therapy.
        J Thorac Oncol. 2021; 16: 1798-1809
        • Hata A.N.
        • Niederst M.J.
        • Archibald H.L.
        • et al.
        Tumor cells can follow distinct evolutionary paths to become resistant to epidermal growth factor receptor inhibition.
        Nat Med. 2016; 22: 262-269
        • Sharma S.V.
        • Lee D.Y.
        • Li B.
        • et al.
        A chromatin-mediated reversible drug-tolerant state in cancer cell subpopulations.
        Cell. 2010; 141: 69-80
        • Oren Y.
        • Tsabar M.
        • Cuoco M.S.
        • et al.
        Cycling cancer persister cells arise from lineages with distinct programs.
        Nature. 2021; 596: 576-582
        • Kurppa K.J.
        • Liu Y.
        • To C.
        • et al.
        Treatment-induced tumor dormancy through YAP-mediated transcriptional reprogramming of the apoptotic pathway.
        Cancer Cell. 2020; 37: 104-122.e12
        • Kurata T.
        • Tamura K.
        • Kaneda H.
        • et al.
        Effect of re-treatment with gefitinib (‘Iressa’, ZD1839) after acquisition of resistance.
        Ann Oncol. 2004; 15: 173-174
        • Yano S.
        • Nakataki E.
        • Ohtsuka S.
        • et al.
        Retreatment of lung adenocarcinoma patients with gefitinib who had experienced favorable results from their initial treatment with this selective epidermal growth factor receptor inhibitor: a report of three cases.
        Oncol Res. 2005; 15: 107-111
        • Shaffer S.M.
        • Dunagin M.C.
        • Torborg S.R.
        • et al.
        Rare cell variability and drug-induced reprogramming as a mode of cancer drug resistance.
        Nature. 2017; 546: 431-435
        • Guler G.D.
        • Tindell C.A.
        • Pitti R.
        • et al.
        Repression of stress-induced LINE-1 expression protects cancer cell subpopulations from lethal drug exposure.
        Cancer Cell. 2017; 32: 221-237.e13
        • Raha D.
        • Wilson T.R.
        • Peng J.
        • et al.
        The cancer stem cell marker aldehyde dehydrogenase is required to maintain a drug-tolerant tumor cell subpopulation.
        Cancer Res. 2014; 74: 3579-3590
        • Maynard A.
        • McCoach C.E.
        • Rotow J.K.
        • et al.
        Therapy-induced evolution of human lung cancer revealed by single-cell RNA sequencing.
        Cell. 2020; 182: 1232-1251.e22
        • Stewart E.L.
        • Mascaux C.
        • Pham N.A.
        • et al.
        Clinical utility of patient-derived xenografts to determine biomarkers of prognosis and map resistance pathways in EGFR-mutant lung adenocarcinoma.
        J Clin Oncol. 2015; 33: 2472-2480
        • Bankhead P.
        • Loughrey M.B.
        • Fernández J.A.
        • et al.
        QuPath: open source software for digital pathology image analysis.
        Sci Rep. 2017; 716878
        • Deshwar A.G.
        • Vembu S.
        • Yung C.K.
        • Jang G.H.
        • Stein L.
        • Morris Q.
        PhyloWGS: reconstructing subclonal composition and evolution from whole-genome sequencing of tumors.
        Genome Biol. 2015; 16: 35
        • Macosko E.Z.
        • Basu A.
        • Satija R.
        • et al.
        Highly parallel genome-wide expression profiling of individual cells using nanoliter droplets.
        Cell. 2015; 161: 1202-1214
        • Trapnell C.
        • Cacchiarelli D.
        • Grimsby J.
        • et al.
        The dynamics and regulators of cell fate decisions are revealed by pseudotemporal ordering of single cells.
        Nat Biotechnol. 2014; 32: 381-386
        • Finak G.
        • McDavid A.
        • Yajima M.
        • et al.
        MAST: a flexible statistical framework for assessing transcriptional changes and characterizing heterogeneity in single-cell RNA sequencing data.
        Genome Biol. 2015; 16: 278
        • Subramanian A.
        • Tamayo P.
        • Mootha V.K.
        • et al.
        Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles.
        Proc Natl Acad Sci U S A. 2005; 102: 15545-15550
        • Kuleshov M.V.
        • Jones M.R.
        • Rouillard A.D.
        • et al.
        Enrichr: a comprehensive gene set enrichment analysis web server 2016 update.
        Nucleic Acids Res. 2016; 44: W90-W97
        • Reich M.
        • Liefeld T.
        • Gould J.
        • Lerner J.
        • Tamayo P.
        • Mesirov J.P.
        GenePattern 2.0.
        Nat Genet. 2006; 38: 500-501
        • Lara-Guerra H.
        • Chung C.T.
        • Schwock J.
        • et al.
        Histopathological and immunohistochemical features associated with clinical response to neoadjuvant gefitinib therapy in early stage non-small cell lung cancer.
        Lung Cancer. 2012; 76: 235-241
        • Laskowitz D.T.
        • Horsburgh K.
        • Roses A.D.
        Apolipoprotein E and the CNS response to injury.
        J Cereb Blood Flow Metab. 1998; 18: 465-471
        • Yao X.
        • Gordon E.M.
        • Figueroa D.M.
        • Barochia A.V.
        • Levine S.J.
        Emerging roles of apolipoprotein E and apolipoprotein A-I in the pathogenesis and treatment of lung disease.
        Am J Respir Cell Mol Biol. 2016; 55: 159-169
        • Tirosh I.
        • Izar B.
        • Prakadan S.M.
        • et al.
        Dissecting the multicellular ecosystem of metastatic melanoma by single-cell RNA-seq.
        Science. 2016; 352: 189-196
        • Blakely C.M.
        • Pazarentzos E.
        • Olivas V.
        • et al.
        NF-κB-activating complex engaged in response to EGFR oncogene inhibition drives tumor cell survival and residual disease in lung cancer.
        Cell Rep. 2015; 11: 98-110
        • Kumar M.S.
        • Hancock D.C.
        • Molina-Arcas M.
        • et al.
        The GATA2 transcriptional network is requisite for RAS oncogene-driven non-small cell lung cancer.
        Cell. 2012; 149: 642-655
        • Tessema M.
        • Yingling C.M.
        • Snider A.M.
        • et al.
        GATA2 is epigenetically repressed in human and mouse lung tumors and is not requisite for survival of KRAS mutant lung cancer.
        J Thorac Oncol. 2014; 9: 784-793
        • Foggetti G.
        • Li C.
        • Cai H.
        • et al.
        Genetic determinants of EGFR-driven lung cancer growth and therapeutic response in vivo.
        Cancer Discov. 2021; 11: 1736-1753
        • Namani A.
        • Liu K.
        • Wang S.
        • et al.
        Genome-wide global identification of NRF2 binding sites in A549 non-small cell lung cancer cells by ChIP-Seq reveals NRF2 regulation of genes involved in focal adhesion pathways.
        Aging (Albany NY). 2019; 11: 12600-12623
        • Del Bubba M.
        • Di Serio C.
        • Renai L.
        • et al.
        Vaccinium myrtillus L. extract and its native polyphenol-recombined mixture have anti-proliferative and pro-apoptotic effects on human prostate cancer cell lines.
        Phytother Res. 2021; 35: 1089-1098
        • Phuchareon J.
        • McCormick F.
        • Eisele D.W.
        • Tetsu O.
        EGFR inhibition evokes innate drug resistance in lung cancer cells by preventing Akt activity and thus inactivating Ets-1 function.
        Proc Natl Acad Sci U S A. 2015; 112: E3855-E3863
        • Taniguchi H.
        • Yamada T.
        • Wang R.
        • et al.
        AXL confers intrinsic resistance to osimertinib and advances the emergence of tolerant cells.
        Nat Commun. 2019; 10: 259
        • Wang R.
        • Yamada T.
        • Kita K.
        • et al.
        Transient IGF-1R inhibition combined with osimertinib eradicates AXL-low expressing EGFR mutated lung cancer.
        Nat Commun. 2020; 11: 4607
        • Suda K.
        • Mitsudomi T.
        Drug tolerance to EGFR tyrosine kinase inhibitors in lung cancers with EGFR mutations.
        Cells. 2021; 10: 1590
        • Kalluri R.
        • Zeisberg M.
        Fibroblasts in cancer.
        Nat Rev Cancer. 2006; 6: 392-401
        • Hu H.
        • Piotrowska Z.
        • Hare P.J.
        • et al.
        Three subtypes of lung cancer fibroblasts define distinct therapeutic paradigms.
        Cancer Cell. 2021; 39: 1531-1547.e10
        • Tsukui T.
        • Sun K.H.
        • Wetter J.B.
        • et al.
        Collagen-producing lung cell atlas identifies multiple subsets with distinct localization and relevance to fibrosis.
        Nat Commun. 2020; 11: 1920
        • Dominguez C.X.
        • Müller S.
        • Keerthivasan S.
        • et al.
        Single-cell RNA sequencing reveals stromal evolution into LRRC15+ myofibroblasts as a determinant of patient response to cancer immunotherapy.
        Cancer Discov. 2020; 10: 232-253
        • Biffi G.
        • Oni T.E.
        • Spielman B.
        • et al.
        IL1-induced JAK/STAT signaling is antagonized by TGFβ to shape CAF heterogeneity in pancreatic ductal adenocarcinoma.
        Cancer Discov. 2019; 9: 282-301
        • Sperb N.
        • Tsesmelis M.
        • Wirth T.
        Crosstalk between tumor and stromal cells in pancreatic ductal adenocarcinoma.
        Int J Mol Sci. 2020; 21: 5486
        • Kieffer Y.
        • Hocine H.R.
        • Gentric G.
        • et al.
        Single-cell analysis reveals fibroblast clusters linked to immunotherapy resistance in cancer.
        Cancer Discov. 2020; 10: 1330-1351
        • Öhlund D.
        • Handly-Santana A.
        • Biffi G.
        • et al.
        Distinct populations of inflammatory fibroblasts and myofibroblasts in pancreatic cancer.
        J Exp Med. 2017; 214: 579-596
        • Buechler M.B.
        • Pradhan R.N.
        • Krishnamurty A.T.
        • et al.
        Cross-tissue organization of the fibroblast lineage.
        Nature. 2021; 593: 575-579
        • Elyada E.
        • Bolisetty M.
        • Laise P.
        • et al.
        Cross-species single-cell analysis of pancreatic ductal adenocarcinoma reveals antigen-presenting cancer-associated fibroblasts.
        Cancer Discov. 2019; 9: 1102-1123
        • Liu T.
        • Zhang L.
        • Joo D.
        • Sun S.C.
        NF-κB signaling in inflammation.
        Signal Transduct Target Ther. 2017; 217023
        • Johnson D.E.
        • O’Keefe R.A.
        • Grandis J.R.
        Targeting the IL-6/JAK/STAT3 signalling axis in cancer.
        Nat Rev Clin Oncol. 2018; 15: 234-248
        • Kim S.M.
        • Kwon O.J.
        • Hong Y.K.
        • et al.
        Activation of IL-6R/JAK1/STAT3 signaling induces de novo resistance to irreversible EGFR inhibitors in non-small cell lung cancer with T790M resistance mutation.
        Mol Cancer Ther. 2012; 11: 2254-2264
        • Navab R.
        • Strumpf D.
        • Bandarchi B.
        • et al.
        Prognostic gene-expression signature of carcinoma-associated fibroblasts in non-small cell lung cancer.
        Proc Natl Acad Sci U S A. 2011; 108: 7160-7165
        • Wang C.
        • Lin G.H.
        • McPherson A.J.
        • Watts T.H.
        Immune regulation by 4-1BB and 4-1BBL: complexities and challenges.
        Immunol Rev. 2009; 229: 192-215
        • Nan J.
        • Du Y.
        • Chen X.
        • et al.
        TPCA-1 is a direct dual inhibitor of STAT3 and NF-κB and regresses mutant EGFR-associated human non-small cell lung cancers.
        Mol Cancer Ther. 2014; 13: 617-629