Therapy-Induced Evolution of Human Lung Cancer Revealed by Single-Cell RNA Sequencing

1Heterogeneity is a property of many biological systems and diseases such as cancer. Biological plasticity in cancer cells is one form of heterogeneity that allows for early adaptation to treatment and limits the success of precision approaches for cancer treatment. In addition to cancer-cell intrinsic heterogeneity, cells within the tumor microenvironment (TME) further contribute to tumor heterogeneity in a cancer cell extrinsic manner. While these tumor compartments and tumor heterogeneity have been characterized in many cancer subtypes, our understanding of how these properties evolve and interact longitudinally in response to systemic treatment remains incomplete, particularly in metastatic tumors.

2Many oncogene-driven cancers such as those with alterations in EGFR, ALK, ROS1, and BRAF are treated with targeted therapies against the cognate oncoprotein. This has led to improvements in the clinical outcomes of metastatic solid cancers such as lung cancer and melanoma as well as hematologic malignancies. However, tumors typically respond incompletely and then regrow after acquiring drug resistance. Bulk tumor sampling after progression on targeted therapy has identified resistance mechanisms and demonstrated that tumors become increasingly molecularly heterogeneous following treatment.

3Single-cell RNA sequencing (scRNA-seq) is one approach to dissect the heterogeneity of complex biological systems. There is currently a paucity of single-cell studies that sample metastatic malignancies and prior scRNA-seq studies of metastatic disease largely focused on single treatment time points. This is due, in part, to challenges associated with obtaining high-quality samples of metastatic human tumors, particularly at multiple treatment time points.

4By developing a custom pipeline, we performed scRNA-seq analyses on advanced-stage NSCLC samples that were obtained from patients before initiating systemic targeted therapy (TKI naive [TN]), at the residual disease (RD) state, which includes samples taken at any time during treatment with targeted therapy while the tumor was regressing or stable by clinical imaging (RD), and upon subsequent progressive disease as determined by clinical imaging, at which point the tumors showed acquired drug resistance (progression [PD]).

5We used scRNA-seq to profile 49 samples (45 lung adenocarcinomas, 1 squamous cell carcinoma, and 3 tumor adjacent tissues [TATs]) (Figure 1A), corresponding to 30 individual patients. We used a customized workflow to isolate viable single cells primarily from small tissue samples as well as surgical resections (Figure 1B). Samples were categorized into three separate time points (TN, RD, or PD) and further subcategorized by oncogenic driver (Figure 1C). Collection time for RD samples is illustrated in Figure S1A. Additional sample details and patient demographics are included in Table S1.

6Gene-expression profiles of 23,261 cells were retained after quality control filtering. Following gene-expression normalization, we performed principal-component analysis (PCA) and clustered cells using graph-based clustering on the informative PCA space (n = 20). The resulting cell clusters were annotated as immune, stromal (fibroblasts, endothelial cells, and melanocytes), or epithelial cells (Figure 1D) by established marker genes (Table S2). Epithelial cells (n = 5,581) were subsetted and re-clustered into 26 discrete epithelial clusters (Figure S1B). The number of cells for each cell type and the analyses that each sample was utilized for are detailed in Table S1.

7Given the association between cancer and large-scale chromosomal alterations, we utilized copy-number variation (CNV) from RNA expression to classify epithelial cells as either cancer or non-cancer; compared to fibroblasts and endothelial cells (controls), cancer cells displayed larger changes from relative expression intensities across the genome (Figure S1C). Three TAT samples were included in this analysis and the majority of cells originating from these samples were classified as non-cancerous (Table S1). We compared the average CNV score of samples among treatment time points (TN, RD, PD) and found it to be consistent. The non-cancer epithelial cell clusters (n = 16) were further annotated into cell subtypes (Figures S1D and S1E).

8As noted by others, we found that cancer cells expressed an elevated number of unique genes compared to non-cancer cells (Figure S1F). The difference in the number of uniquely expressed genes was not explained by sequencing depth (Pearson correlation = 0.19).

9Cancer cells were identified in 44 of the original 49 tumor biopsy samples including a small fraction of cells originating from each of the TAT samples (0.57%–1.8% of total TAT obtained cells). Given that TAT cells may represent an intermediate cellular state between normal and cancer cells, their presence at low frequency is unsurprising and has been described previously.

10All cancer cells (n = 3,754) were re-clustered, resulting in 25 unique clusters (Figures S2A and S2B). For each of the 25 clusters, we calculated the number of cells of the highest contributing individual patient over the total number of cells for that cluster for both non-cancer and cancer epithelial cells (patient occupancy) (Figures S2C–S2E). The majority of cancer cell clusters were patient specific, having high patient occupancy scores, similar to prior reports. Conversely, non-cancer cell types exhibited lower patient occupancy (Figure S2E). Thus, patient-specific malignant cell clustering reflected the unique molecular signatures of an individual patient’s tumor rather than technical artifact.

11Additional genetic alterations can co-exist with a primary targetable oncogenic “driver” alteration (e.g., oncogenic EGFR, ALK, KRAS) and may help promote tumor progression and limit therapy response. We queried scRNA-seq transcripts from each cancer cell to identify somatic alterations (Figures 2A–2C; Figures S2G–S2I). In the 44 of 49 biopsy samples that contained cancer cells, we identified 20 samples harboring the clinically known oncogenic driver (Figure 2B; Table S3). This percentage is consistent with the potential drop-out occurrence in scRNA-seq analyses. In 24 samples where we did not identify the clinically known oncogenic driver, no cells expressed the gene of interest, thus not allowing mutation detection for that gene (Figure S2H; Table S3). In 11 of the 20 samples (55%) where we identified the clinically actionable oncogene, we also identified an additional oncogenic alteration that was not detected in clinical-grade bulk nucleic acid testing of tumor from the same patient (i.e., occult genetic alterations) (Figure 2B; Figure S2H). An example is sample LTS47. This tumor was determined to harbor an EML4-ALK oncogenic gene rearrangement by clinical-grade bulk DNA analysis. scRNA-seq additionally revealed that this sample contained cancer cells harboring KRAS G13D and KRAS G12C occult mutations (Figure S2F). Given the potential of dropout in scRNA-seq data, we cannot conclude whether the ALK fusion and the KRAS mutations co-exist within the same cell. However, neither population of KRAS mutant cells showed evidence of the ALK gene rearrangement. This sample was obtained from the patient after multiple lines of therapy, which could have allowed for evolution of multiple mechanisms of resistance. Loss of an oncogenic driver is also a mechanism of resistance, although given the limitation of scRNA-seq we were not able to determine whether this mechanism of resistance applies to this case.

12We also queried scRNA-seq data for mutations from the COSMIC (Catalogue of Somatic Mutations In Cancer) lung adenocarcinoma tier 1 mutations (Table S2),. Many of the mutations we identified had not been previously reported by the clinical-grade assay conducted on a patient’s tumor despite having been included in the clinical panel (Figure 2C; Figure S2I; Tables S2 and S3). Though this may reflect differences in biopsy technique or tumor clonality at the time of clinical testing, these results also demonstrate that clinical-grade bulk DNA-based testing may underestimate tumor heterogeneity.

13To assess the clinical outcomes of patients harboring multiple oncogenic alterations and to determine the broader translational impact of our findings, we utilized the MSK-Impact NSCLC dataset. Those patients whose tumor showed greater than or equal to 2 mutations from the tier 1 COSMIC mutation set detected in the scRNA-seq profiling (mutation high) had significantly lower overall survival (OS) compared to those patients whose tumor had less than 2 COSMIC tier 1 mutations (mutation low) (p < 0.01; Figure S2J). Thus, scRNA-seq analysis can provide increased granularity into cancer cell genomic heterogeneity and provides insight into the transcriptionally represented mutational landscape.

14We hypothesized that defining the biological programs activated in cancer cells during therapy response may identify signaling pathways that promote adaptation and survival of cancer cells that comprise RD during initial treatment. We compared the transcriptional profiles of individual cancer cells obtained from tumor samples from TN to RD (Table S4) and focused on the 629 significantly (p < 0.05) upregulated genes in RD cancer cells as a proxy for evidence of pathway activation. We found numerous genes associated with cancer-associated pathways (Table S5). Importantly, we found that RD cancer cells expressed decreased proliferation marker genes compared to TN and PD, consistent with the expectation that during targeted treatment persisting cancer cells are generally less proliferative (Figure S3A).

15Interestingly, we identified an alveolar cell gene-expression signature composed of 17 established gene markers of alveolar cells that showed significantly increased expression in RD versus TN time points (p < 0.0001; Figure 3A; Figure S3B; Table S2). Alveolar cells are comprised of alveolar type 1 (AT1) and type 2 (AT2) subtypes and form the lining of the lung alveoli. AT2 cells produce surfactants and can act as stem-like progenitor cells, which become active and proliferate in the setting of diverse types of lung injury and are suspected to be the cell of origin in oncogene-driven lung cancers. AT1 cells are the dominant population in alveoli and mediate gas exchange and, when injured or dying, can release proliferation and regenerative signals. AT1 cells contain two population subtypes HOPX+/IGFBP2+ and HOPX+/IGFBP2–, the latter representing the cell population which maintains cellular plasticity and can proliferate as well as trans-differentiate into AT2 cells allowing for tissue regeneration after injury. The alveolar signature we detected in the cancer cells at RD includes both AT1- and AT2-associated genes (Table S2), including AQP4, SFTPB/C/D, CLDN18, FOXA2, NKX2-1, and PGC for AT2 cells and AGER, HOPX, and IGFBP2 for AT1 cells (;; Figure S3B). Additionally, the alveolar cell state we identified in cancer cells was not derived from mis-annotated non-cancer alveolar cells within our cancer cell populations (Figure S3C).

16We validated the activation of the alveolar cell signature at RD using orthogonal approaches. First, we used an established preclinical model consisting of patient-derived EGFR mutant NSCLCs (PC9) to develop analogs of the TN, RD, and PD clinical states. Using RT-PCR, we measured the expression of NKX2-1, a hallmark alveolar cell signature gene upregulated in RD clinical samples and found a significantly higher (p < 0.001) expression in the persister state cells compared to control and acquired resistance state cells (Figure 3B). This suggests that the alveolar signature identified from the clinical scRNA-seq analysis can be reproduced under controlled conditions in vitro. Furthermore, immunohistochemical (IHC) analysis showed induction of AQP4 protein expression, another marker of the alveolar cell signature, at the plasma membrane of RD clinical samples compared to both TN and PD clinical samples (Figure 3C; Figure S3D).

17We tested whether the alveolar signature is a clinically relevant biomarker of patient survival in the TCGA lung adenocarcinoma bulk RNA-seq dataset (https://www.cancer.gov/about-nci/organization/ccg/research/structural-genomics/tcga, ). We found a significant (p < 0.0001) association between high expression of our alveolar signature in these tumors and improved patient OS when compared to patients whose tumors showed a lower alveolar expression signature (Figure 3D; Table S6).

18These findings support the assertion that there is a distinct alveolar gene-expression signature characterizing RD cancer cells, associated with improved patient survival. A plausible model is that the identified alveolar signature that is activated in RD cancer cells reflects a cell injury and repair signature, reminiscent of non-cancerous AT1 and AT2 cells. Increased expression of this signature could lead to repair and escape of cell death during treatment to support cancer cell persistence, while at the same time constituting a less aggressive malignant state. This is consistent with the notion that RD represents a persister cell state observed in preclinical models of slow-cycling cancer cells that survive without rapid proliferation (as in Figure S3A), as a prelude to the onset of aggressive tumor progression upon absolute drug resistance.

19The molecular details of the alveolar and cell injury repair signature are notable. In our RD cohort, the WNT/β-catenin-associated pathway genes SUSD2 and CAV1 exhibited increased expression (Table S5). We used IHC analysis of both SUSD2 and CTNNB1 (β-catenin) protein expression (Figures S3E and S3F) to validate the observed transcriptional changes. In agreement with our scRNA-seq findings, we found significantly increased membrane SUSD2 and significantly increased nuclear CTNNB1 (β-catenin) in the RD state compared with both TN and PD. Additionally, the comparison of nuclear localization of CTNNB1 in a unique series of paired TN and RD samples obtained from EGFR(AZ) or ROS1(NC) patients treated with neoadjuvant TKI on one of two clinical trials (osimertinib: NCT03433469 or crizotinib: NCT03088930) is shown in Figure S3G. SUSD2 is an activated downstream target of the WNT pathway, while CAV1 can promote nuclear localization of β-catenin (CTNNB1) and transcriptional activation of the WNT/β-catenin pathway. In NSCLCs, the WNT/β-catenin signaling pathway contributes to tumorigenesis, repair, and regeneration after cell injury. The self-renewal and injury response in AT2 cells specifically can utilize the WNT/β-catenin signaling pathway. Additionally, in EGFR mutant NSCLC activation of the WNT/β-catenin pathway may limit EGFR inhibitor response and may promote survival of a persister cell population during EGFR inhibitor therapy in vitro. Overall, the RD state is characterized by signals of cellular injury and survival, which act, in part, through the WNT/β-catenin pathway, which may be therapeutically targetable.

20The clinical data suggest that WNT/β-catenin activation is engaged early during treatment to facilitate the development of RD and drug tolerant persister cells during primary EGFR or ALK targeted therapy. To further explore the therapeutic potential of the WNT pathway findings, we utilized patient-derived PC9 cells as an EGFR mutant NSCLC model and H3122 cells as a model for ALK fusion-driven NSCLCs. We tested the hypothesis that upfront blockade of WNT/β-catenin signaling together with oncogenic EGFR or ALK would decrease the number of cells surviving initial treatment and increase the depth of response from the outset of therapy. Parental cells were treated with an IC50 (inhibitor concentration yielding a 50% decrease in cell number) dose of the appropriate EGFR or ALK TKI (osimertinib or alectinib, respectively). Two different WNT/β-catenin pathway inhibitors XAV939 and PRI-724 in four previously reported concentrations or combination therapy thereof were tested. Our in vitro results support our hypothesis by demonstrating that the upfront inhibition of the WNT/β-catenin pathway in combination with the cognate TKI led to a significant and dose dependent decrease in cell confluency and increased depth of response (Figures 3G–3J; Figures S3H–S3O).

21When comparing cancer cells from TN and PD samples, we found 901 differentially upregulated genes in PD cancer cells (Table S4). Within those genes, we identified genes involved in the kynurenine pathway and multiple genes and pathways associated with oncogenesis and inflammation (Table S5).

22We observed a significant (p < 0.0001) increase in the expression of IDO1, KYNU, and QPRT genes involved in the kynurenine pathway, in PD and TN cancer cells (Figure 3K; Figure S3P; Table S2) Expression of these genes can result in immunosuppressive behavior indicating that cancer cells within PD tumors may directly inhibit the activity of the immune system. The identification of this pathway as a mediator of immune suppression within PD tumors has important potential therapeutic implications, as IDO1 is upregulated in many cancers. Multiple clinical trials have attempted to block this pathway using IDO1 inhibitors as a monotherapy as well as in combination with immune checkpoint inhibitors or hormone therapy, albeit with limited success. QPRT also exhibited increased expression (p < 0.05) specifically at acquired resistance of EGFR inhibitor treatment (i.e., the analog of clinical PD) in our in vitro model using PC9 cells (Figure 3L), reinforcing the assertion that this pathway is indicative of cancer progression under the selective pressure of treatment.

23To further demonstrate the clinical relevance of the kynurenine pathway, we again used the TCGA lung adenocarcinoma RNA-seq dataset. Higher tumor expression of this signature was a biomarker of worse OS in patients (p < 0.05; Figure 3M; Table S6). This is consistent with the notion that activation of this pathway leads to immunosuppression and an inability of the immune system to effectively surveil and eradicate cancer cells.

24We compared cancer cells from RD and PD patient samples to elucidate the differences that occur during the outgrowth of PD from RD and found a total of 2,182 genes which had significantly (p < 0.001) increased expression in either RD or PD (NRD = 1,121, NPD = 1,061) (Table S4). Among the differentially overexpressed genes at RD were genes associated with the alveolar cell signature, cell growth, differentiation, cell motility, and tumor suppression (Table S5). RD cancer cells overexpress surfactant genes (SFTPB/C/D and SFTA3), which are part of the alveolar cell signature (Figure 3A; Figure S3B). Furthermore, NKX2-1 and NFIX were overexpressed in RD cancer cells and are associated with decreased cell motility. Low expression of NKX2-1 leads to loss of differentiation and enhanced tumor seeding ability. The collective findings arising from this and the previous RD cancer cell analyses suggest that an injury-repair and regenerative cell state may promote cancer cell indolence, increased tumor control, and improved clinical outcomes.

25By contrast, PD cancer cells differentially overexpressed genes associated with invasion, cell-to-cell communication, differentiation, and immune modulation (Table S5). Several genes in the plasminogen activation pathway were significantly overexpressed (ANXA2, PLAT, PLAUR, PLAU) (Figure 3N) along with the plasminogen inhibitor SERPINE1 (PAI1) (p < 0.0001, Figure 3O; Figure S3Q). ANXA2 and PLAUR are the receptor proteins in the plasminogen activation cascade and involved in inflammation, angiogenesis, invasion, and metastasis, via degradation of the extracellular matrix. Signaling is initiated when ANXA2 or PLAU binds to PLAT (uTa) or PLAU (uPa), respectively. Plasminogen is then degraded to plasmin through the activity of PLAT and/or PLAU leading to activation of metalloproteinases and degradation of fibrin. SERPINE1 shows increased expression in a number of cancer subtypes and plays important roles in cell adhesion, invasion, tumor vascularization, radio-resistance, and immunosuppression. High expression of the plasminogen activation signature correlated with worse patient OS (p < 0.01) within the TCGA lung adenocarcinoma RNA-seq dataset and cohort (Figure 3P; Table S6). Similarly, in this independent dataset high expression of SERPINE1 was associated with worse OS (p < 0.05) (Figure 3Q; Table S6). EGFR inhibitor therapy can induce expression of SERPINE1 and EGFR mutant patients with greater than 2-fold induction of SERPINE1 (PAI1) plasma levels during EGFR inhibitor treatment demonstrated shorter progression-free survival. Collectively, our scRNA-seq findings shed light on the clinical relevance and potential role of the plasminogen activation cascade in inferior clinical outcomes and targeted therapy resistance.

26Additionally, we found several gap-junction proteins differentially overexpressed in PD cancer cells compared to RD cancer cells (p < 0.0001, Figure 3R; Table S2). Gap-junction proteins (e.g., connexins) are integral membrane proteins that allow for cytosolic exchange of ions, metabolites and secondary messengers between cells. While some have been identified as tumor suppressors, we found that high expression of GJB2/3/5 (Figure S3R; Table S6) was linked to worse survival in the TCGA lung adenocarcinoma RNA-seq dataset (p < 0.001) (Figure 3S). These collective findings suggest a pro-tumor effect not only in our cohort but also in NSCLCs more generally.

27Within cancer cells, we identified a rich complexity of clinically relevant, expressed mutations that may impact therapy response. Furthermore, evaluation of transcriptional profiles of individual cancer cells across different treatment time points identified several clinically relevant cell-state changes (Figure S4A). We found that the identified treatment time point signatures largely persist irrespective of the type of oncogenic driver mutation (EGFR or ALK) (Figures S4B–S4K) and of biopsy site (primary or metastatic) (data not shown). While we found these cancer cell signatures are robust, it is important to acknowledge that there is patient heterogeneity among samples (Figures S4L–S4O).

28Obtaining consecutive clinical tumor biopsies from individual advanced-stage lung cancer patients before and during treatment is challenging given that most tumors regress by 50% or greater, albeit incompletely, during TKI treatment. Nevertheless, we obtained samples from the same primary tumor site from 3 treatment time points from a patient (TH226) whose tumor contained a standard EGFR exon 19 deletion oncogenic mutation and was treated with the EGFR inhibitor osimertinib (Figures S5A–S5C). In all 3 biopsies, we identified by scRNA-seq RNA expression of the EGFR exon 19 driver mutation in the cancer cells and several other mutations of interest (Figure S5D).

29When comparing TH226 to the rest of the scRNA-seq dataset, we found overlapping differentially expressed genes and signatures (Table S4; Figures S5E–S5H). Intriguingly, we also found numerous genes associated with squamous cell differentiation (KRT16, KRT14, KRT6A, KRT5, CLCA2, PKP1, ANXA8, DSG3) overexpressed at PD compared to TN and RD time points (p < 0.0001, Figure S5I, Tables S4 and S5). This is particularly interesting given that the patient’s lung tumor biopsy at PD demonstrated a histologic shift to squamous cell carcinoma from that of prior biopsies that showed pure adenocarcinoma histology (Figure S5C). Histologic transformation to squamous cell carcinoma is a mechanism of EGFR inhibitor resistance in EGFR mutant NSCLCs. Thus, scRNA-seq has the power to provide a high resolution, gene and pathway level view of biological and histological plasticity that arises during cancer drug treatment.

30We next addressed the evolution of the TME during targeted treatment. Immune cells (n = 13,431) were clustered and annotated (Figure 4A; Tables S1 and S4). In contrast to clusters of cancer cells, which were predominantly patient specific (Figures S2C and S2D), immune cell-type clusters showed low patient occupancy (Figure 4B). This is consistent with the expectation of finding common immune cell phenotypes across patients and samples.

31We compared the immune cell composition across all 3 time points, expressed as the correlation between fractional immune cell abundance vectors. The immune composition within RD was the most dissimilar from the other two treatment states (r = 0.78 versus TN samples, r = 0.82 versus PD samples, Pearson’s correlation coefficient) (Figure S6A). Across all treatment time points, T cells and macrophages were the dominant cell populations and demonstrated an inversion in relative abundance during tumor response and resistance to treatment, a finding we examined further as described below (Figure 4C). T cells comprised a larger fraction of all immune cells within the TME at RD compared to TN or PD samples (27% T cells TN, 46% RD, 31% PD). Macrophage infiltration followed the inverse pattern, with a decrease in macrophages at RD compared to TN and PD (37% macrophages TN, 21% RD, 37% PD).

32In 2 patients, we examined immune cells from available matched tumor biopsies obtained at different treatment time points (TH226 and TH266, Figures S6B and S6C, respectively). In 2 tumor biopsies available for patient TH266, both macrophages and T cells showed reduction in the fraction of macrophages and an increase in the fraction of T cells from TN to RD, findings which match the entire cohort (Figure S6D). TH226 exhibited a similar pattern with the fraction of macrophages decreasing at RD after initiation of treatment and increasing again at PD (Figure S6E). We validated our findings on tissue samples using immunofluorescence staining (Figures 4D–4F; Figure S6F). Additionally, we deconvoluted TCGA bulk transcriptome data for NSCLCs into fractions of immune cells types (see STAR Methods) and found that TCGA samples with high fractions of macrophages had significantly worse OS (p < 0.01) (Figure S6G). This supports the clinical relevance of our observations and is consistent with prior reports associating macrophage infiltration with poor prognosis in patients who undergo surgical resection of early-stage NSCLCs, as well as with worse progression-free survival during EGFR TKI therapy.

33These findings are particularly intriguing given their similarity to melanoma tumors treated with PD-1 inhibitor, albeit here in the distinct context of oncoprotein-targeted therapy in lung cancer. Specifically, an increase in the number of CD8+ T cells and natural killer (NK) cells and a decrease in M1 macrophages were observed in melanoma during PD-1 inhibition. There may be common responses in NK/T cells and macrophages during treatment across different tumor histologies and treatments. Hence, conserved approaches to targeting RD across different cancer subtypes and therapeutic modalities may exist, an area for future investigation.

34Macrophages from lung tumor biopsies (n = 1,379) were clustered into 5 distinct groups (Figure S7A) followed by differential gene expression in each resulting cluster (Figure S7B; Table S4). In addition, we calculated the fraction of cells originating from each of the three treatment groups in each of the 5 macrophage clusters (Figure 5A).

35Cluster MF0, which was slightly enriched in TN cells, was characterized by expression of genes associated with an immunosuppressive phenotype (C1QC, GPNMB, APOE, TREM2, FOLR2) (Figure 5B; Figure S7B). Clusters MF1 and MF3 were enriched at RD. Macrophage cluster MF1 expressed features associated previously with tumor-infiltrating myeloid derived suppressor cells (FCN1, VCAN, S100A8, S100A9) and with THBS1 and PTX3, which are associated with resolution of inflammation, wound healing, and with inhibition of IL-1β (Figure 5B; Figure S7B). Cluster MF3 expressed genes associated with pro-inflammatory response to tissue damage (CLEC2D, IL7R, OGT) and with promoting inflammatory signaling (FYN, DUSP4, FOXO1). Macrophages in this cluster express CCL5, a cytokine that has previously been associated with promotion of residual HER2+ breast cancer survival following HER2-targeted therapy. Cluster MF4 consisted of proliferating myeloid cell populations (TOP2A, MKI67) and did not significantly differ between groups.

36Macrophages at PD were overrepresented in group MF2 (Figure 5A) and expressed pro-inflammatory cytokines CXCL9, CXCL10, and CXCL11 (Figure 5B; Figure S7B), which favor lymphocyte recruitment into the TME. Top differentially expressed genes in this population also included the guanylate-binding family proteins GBP1 and GBP5, which are induced in IFN-γ-activated macrophages and promote inflammatory signaling within the innate immune system via inflammasome assembly (Figure 5B; Figure S7B). Despite the expression of pro-inflammatory genes within the MF2 macrophages, the top differentially expressed gene within this group of PD-specific macrophages was IDO1 (Figure 5B). IDO1 is induced by inflammation within the TME and promotes a tolerogenic environment through immunosuppressive myeloid cell populations, regulatory T cell (Treg) differentiation, and an immunosuppressive cytokine milieu.

37T cells and NK cells (n = 2,226) were analyzed in the same manner as macrophages and resulted in 5 distinct T/NK cell populations (Figure 5C; Figure S7C). These included two populations (TC0, TC4) enriched in TN samples and 2 populations (TC1, TC2) enriched at PD (Figures 5C and 5D; Figure S7D). There was a high overall fraction of T cells in RD tumors (Figure 4C), and there was no single T cell cluster that demonstrated an excess of T cells in RD (Figure 5C).

38Both TN and PD T cells demonstrated a relative decrease in T cell infiltration (Figure 4C). The T cells which were present in the TN state were enriched for T cell populations TC4 and TC0. TC4 expressed markers consistent with a natural killer (NK) or natural killer T cell (NKT) phenotype, including NK cell markers (KIR2DL3, FCGR3A) as well as moderate expression of T cell markers (CD3, CD8). TC0 reflected a naive-like CD8+ phenotype with expression of CCR7, IL7R, and SELL (Figure 5D; Figure S7D). While overall T cell infiltrate remained limited at PD (Figure 4C), there was relative enrichment for T cell phenotypes with immunosuppressive features, including T cell clusters TC1 and TC2 (Figure 5C). TC1 was identified as a T cell cluster with a dysfunctional or exhausted phenotype, characterized by expression of the inhibitory receptors PDCD1 (which encodes for the PD-1 protein) and CTLA4 (Figure 5D). TC2 was composed of Treg cells (expressing FOXP3, IL2RA). Consistent with a relatively immunosuppressive environment, there was additionally reduction at PD in infiltration by the NK/NKT cell cluster TC4.

39Tumor biopsies obtained at RD revealed the presence of a more pro-inflammatory, “hot,” TME, which was absent in TN or PD biopsy samples as manifested by increased overall proportion of T cells and reduced infiltration by regulatory T cells (TC2) (Figure 5E). Compared to the PD state, at RD there were fewer dysfunctional T cells (TC1) and greater NK/NKT cell (TC4) infiltration (Figure 5D; Figure S7D). A population of proliferating tumor infiltrating T cells was shared across all treatment states (TC3) and was slightly enriched within PD samples. This T cell population was characterized by both cytotoxic phenotypes (CD8, GZYMB) and PDL1/CTLA4 expression and may reflect a pre-dysfunctional cytotoxic T cell population.

40In summary, both the TN and PD TME were characterized by the relative predominance of macrophage over T cell infiltration; however, the phenotypic characteristics of these infiltrating immune cells differ between the two groups. At PD, there was infiltration by an IDO1+ macrophage population, of proliferating regulatory T cells, and of dysfunctional T cells which were minimally present at TN and RD. In contrast, the TN state was characterized by a predominance of more classically immunosuppressive M2-like macrophages. By distinction, in RD there was increased infiltration of T cell populations without dysfunctional or immunosuppressive gene-expression patterns and decreased immunosuppressive macrophage infiltration (Figure 5E).

41There remains an incomplete catalog of single-cell transcriptional data that can be used to understand cell states and the therapy-induced evolution of biological heterogeneity of diseases such as cancer, particularly for advanced-stage solid malignancies. Patients with metastatic disease do not routinely receive surgical resection as part of their treatment. Thus, techniques for single-cell profiling that require larger amounts of tissue are not suitable for the interrogation of tissue samples from metastatic disease. Our scRNA-seq analyses of advanced-stage NSCLC biopsies obtained at different treatment time points from patients elucidate the rich mutational and transcriptional diversity within individual tumor samples and the dynamic changes in the transcriptional profiles of cancer cells and the TME composition during treatment. Our findings provide a roadmap that highlights the underlying cellular ecosystem and mechanisms that can inform efforts to better treat oncogene-driven cancers. Our study offers a rare view of the clinically relevant biological processes that characterize RD, which is a treatment phase that is infrequently captured in human solid malignancies.

42Our scRNA-seq data revealed widespread intra-tumoral heterogeneity in oncogenic alterations that are expressed in cancer cells (Figure 2) by demonstrating expression of not only the putative oncogenic driver but also additional oncogenic mutations (Figure 6, #1). This provides a potential explanation for why complete responses to treatment are rare. Tumors harbor the appropriate genetic framework and evolutionary playbook to evolve resistance. These “hard-wired” properties can go undetected by current bulk sampling analysis. Tumor resilience and evolution during therapy are bolstered by the therapy-induced transcriptional plasticity that we demonstrated by scRNA-seq profiling.

43We uncovered transcriptional signatures specific to different treatment time points and clinical states (Figures 3 and Figure 6, #3, #5, #6, and #7). The majority of these signatures were biomarkers of significantly worse OS and were most pronounced at PD. Conversely, we found the alveolar cell signature was enriched at RD and was associated with improved survival. This signature exhibited features consistent with cellular plasticity and injury response, perhaps indicating a treatment-inducedadaptive phenotype that permits the survival of cancer cells, albeit in a less aggressive state. Our data also highlight a connection from the alveolar cell signature to the WNT/β-catenin pathway as a mechanism of injury response and regeneration. Though the WNT/β-catenin pathway is potentially therapeutically targetable,, it will be critical to determine how to best modulate this pathway to impact residual cancer cell survival for clinical benefit.

44A general principle our data highlight is that by employing targeted treatments that take advantage of specific cell states we may be able to engineer cancer (or TME) cell fate(s) to improve therapeutic responses in metastatic solid malignancies. If deployed at the appropriate time, treatments that target liabilities of a specific cell state or prevent further adaptation may help improve patient survival by constraining continued tumor evolution toward complete drug resistance (Table 1).

45The recognition of targetable oncogene-driven NSCLCs as a subset of lung cancer that is poorly responsive to current checkpoint inhibitor immunotherapies necessitates an improved understanding of the immunologic milieu in this patient population. We found a relatively low T cell infiltration in the TME of TN and PD patients (Figures 5C and 5D), consistent with prior reports of low cytotoxic T cell infiltration in treatment-naive EGFR mutant NSCLCs and with an association between EGFR activation and an immunosuppressive phenotype in preclinical models. Our results uncovered an induction of a more inflammatory phenotype during RD on targeted therapy, hallmarked by infiltration of T cells (Figure 6, #4) and decreased infiltration of immunosuppressive macrophages (Figure 5E). This inflammatory state may represent a complement to the alveolar cell, injury repair, and regenerative state present in the cancer cell compartment (described above), with the potential for crosstalk between the cancer cells and TME. These TME changes were transient, as at PD there was enrichment for IDO1-expressing macrophages, regulatory T cells, and other immunosuppressive T cell populations. These are all features of an environment hostile to the establishment of an effective immunologic response (Figure 6, #9 and #10). These findings confirm those reported within preclinical models of EGFR mutant NSCLCs, which demonstrated the potential for a transient immunostimulatory effect after initial EGFR TKI exposure despite the immunosuppressive phenotype observed following longer-term TKI treatment. The induction of a more immunostimulatory phenotype during targeted therapy (i.e., in RD) may offer a window-of-opportunity to introduce novel TME target-based combination therapies earlier during treatment, perhaps around the time of RD in the context of a more favorable TME to increase initial response and consolidate the anti-tumor responsein a multi-modal therapeutic approach.

46Given that cancer cell signaling and the TME are linked, there may be treatment strategies that target both compartments concurrently. The kynurenine pathway is one example. We identified increased kynurenine pathway activation in cancer cells and myeloid cells at PD (Figure 6, #5). IDO1, as a rate-limiting enzyme in the kynurenine pathway, can influence diverse components of the TME including T cell and myeloid cell populations as well as angiogenesis in favor of immunosuppression. The use of IDO1 inhibitors as part of a combination immunotherapy strategy with PD1/PDL1 (programmed cell death protein 1/programmed death-ligand 1) checkpoint inhibitors showed promise in early-phase studies,yet failed to demonstrate improved outcomes in advanced-stage melanoma. We demonstrated distinct evolving TME states, suggesting that there may be a window-of-opportunity at which point kynurenine pathway inhibitors may be more effective (Figure 6; Table 1).

47The scRNA-seq dataset presented here demonstrates the feasibility of performing scRNA-seq on tumor biopsies obtained longitudinally at clinically relevant time points during active targeted treatment of advanced-stage solid malignancy patients. The data provide a foundation to develop strategies for the elimination or neutralization of RD to induce more durable responses for patients with advanced-stage NSCLCs and potentially other solid malignancies across different therapeutic modalities.

48Limitations of our study include the number and diversity of cells and genotypes of individual tumor biopsies due to the use of small needle biopsies or fluid collections versus larger surgical resections. Due to the real-world challenges of tissue acquisition, we acquired matched samples from a small number of individual patients. Because RD is not sampled during standard treatment, there were fewer samples at this disease state. Single-cell-derived transcriptomes are relatively sparse as a consequence of a combination of factors including transcriptional stochasticity, rarity of sampling mRNA molecules, uneven amplification of mRNA molecules during cDNA synthesis and library preparation, and sparse read coverage of library molecules. These challenges limit our ability to perform a saturating mutation analysis using single-cell data.

49Further information and requests for resources should be directed to and will be fulfilled by the Lead Contact, Trever Bivona (https://cran.r-project.org/web/packages/REdaS/index.html).

50This study did not generate new unique reagents.

51The dataset generated during this study is available as an NCBI BioProject #PRJNA591860. All code used to generate the results of this study can be found on github at czbiohub/scell_lung_adenocarcinoma and czbiohub/cerebra. The below methods reference specific code notebooks (script xx) available at czbiohub/scell_lung_adenocarcinoma to analyze data.

52All patients gave informed consent for collection of clinical correlates, tissue collection, research testing under Institutional Review Board (IRB)-approved protocols (CC13-6512 and CC17-658, NCT03433469). Patient demographics are listed in Table S1. Patient studies were conducted according to the Declaration of Helsinki, the Belmont Report, and the U.S. Common Rule.

53PC9 (EGFRexon19del) and H3122 (EML4-ALKv1) cells were purchased from ATCC and grown in a 5% CO2, humidified atmosphere at 37°C. Cultures were maintained using RPMI 1640 medium (GE Healthcare) supplemented with 10% (v/v %) fetal bovine serum (VWR), 100 IU/mL penicillin and 100 μg/mL streptomycin (GIBCO).

54Formalin-fixed paraffin embedded (FFPE), frozen, and fresh tissue samples were obtained according to the safety standards of the interventional radiologist, pulmonologist, or surgeon. Demographic and clinical history for each patient was obtained from chart review and is listed in Table S1. Days until progression were determined based on imaging studies which demonstrated definitive growth of a known tumor site or new extra-CNS metastatic deposits. Residual disease state was determined by serial imaging demonstrating continued reduction or stability tumor with no evidence of progression. Complete details of each patient sample acquisition are outlined in Table S1 and Figure S1A. Additionally, the timing of each sample acquisition is shown in Figure S1A.

55Tissue was first cut into small pieces and placed into a 1.5 mL tube (or multiple tubes if necessary). 1.5 mL of collagenase buffer (10mL DMEM (GE Life Sciences, SH30081.01), 0.20 g Collagenase Type 2 (Worthington Biochemical, LS004176)) was added to the tube and the sample was digested for 30 minutes at 37°C, shaking in a thermomixer @ 800-1000 rpm. The sample was manually agitated by pipetting up and down 15 times then returned to the thermomixer for 25 minutes. After incubation, the sample was removed from the thermomixer, agitated again by pipetting the sample up and down 15 times before passing the sample through a 100-micron filter (Fisherbrand, 22363548) into a new 15 mL falcon tube. The filter was washed with 1-2 mL of collagenase buffer before the sample was spun in the centrifuge at 500xg for 10 minutes. If the resulting cell pellet was red, 0.5 mL RBC lysis buffer (Thermo Fisher Scientific, A1049201) was added to sample tubes and allowed to sit at room temperature for 3 minutes before quenching with 1.0 mL DMEM (GE Life Sciences, SH30081.01) + 6% FBS (Omega Scientific, Inc, FB-11) and spun in the centrifuge at 500xg for 5 minutes. Remaining cells were stained with 10 μl CD45-FITC (Miltenyi Biotec, 130-080-202) and 1 μl of Hoechst stain (Thermo Fisher Scientific, H3570). Samples incubated on ice in the dark for 20 minutes. 1mL of FACS Buffer (1:20 dilution of BSA (Miltenyi Biotec, 130-091-221) in Running Buffer (Miltenyi Biotec, 130-091-221)) was then added to the stained cells and spun at 500xg for 10 minutes before aspirating off supernatant. Cells were resuspended with 0.5 mL of FACS Buffer. PI (Life Technologies, P3566) was added immediately prior to sorting.

56Cells were filtered through a 100 μm strainer (Fisherbrand 22363548), pelleted (500xg, 5 min, 4°C), and resuspended in FACS buffer. Cells were then stained with CD45-FITC (Miltenyi Biotec, 130-080-202) for 20 min at 4°C in the dark. Cells were then pelleted (500xg, 5 min, 4°C) and resuspended in FACS buffer before being transferred to a FACS tube (Falcon 14-956-3C). Sytox Blue dead cell stain (Thermo Fisher Scientific, S34867) was added immediately prior to sorting.

57Lysis plates were created by dispensing 0.4 μL lysis buffer (0.5U Recombinant RNase Inhibitor (Takara Bio, 2313B), 0.0625% Triton™ X-100 (Sigma, 93443-100ML), 3.125 mM dNTP mix (Thermo Fisher, R0193), 3.125 μM Oligo-dT30VN (IDT, 5′AAGCAGTGGTATCAACGCAGAGTACT30VN-3′) and 1:600,000 ERCC RNA spike-in mix (Thermo Fisher, 4456740)) into 384-well hard-shell PCR plates (Biorad HSP3901) using a Tempest liquid handler (Formulatrix). All plates were then spun down for 1 minute at 3220xg and snap frozen on dry ice. Plates were stored at −80°C until used for sorting.

58Cells were sorted into 384-well plates using SH800S (Sony) sorter. Cells were sorted using the “Ultra purity” setting on the sorter. For a typical sort, a FACs tube containing 0.3-1ml the pre-stained cell suspension was vortexed gently and loaded onto the FACS machine. A small number of cells were flowed at low pressure to check cell concentration and amount of debris. Then the pressure was adjusted, flow was paused, the first destination plate was unsealed and loaded. Single-cell sorting was done where half the plate was sorted for CD45+/PI-/Hoechst+ while the second half was sorted for CD45-/PI-/Hoechst+. Immediately after sorting, plates were sealed with a pre-labeled aluminum seal, centrifuged and flash frozen on dry ice.

59cDNA synthesis was performed using the Smart-seq2. Briefly, 384-well plates containing single-cell lysates were thawed on ice followed by first strand synthesis. 0.6 μL of reaction mix (16.7 U/μl SMARTScribe Reverse Transcriptase (Takara Bio, 639538), 1.67 U/μl Recombinant RNase Inhibitor (Takara Bio, 2313B), 1.67X First-Strand Buffer (Takara Bio, 639538), 1.67 μM TSO (Exiqon, 5′-AAGCAGTGGTATCAACGCAGACTACATrGrG+G-3′), 8.33 mM DTT (Bioworld, 40420001-1), 1.67 M Betaine (Sigma, B0300-5VL), and 10 mM MgCl2 (Sigma, M1028-10X1ML)) was added to each well using a Tempest liquid handler or Mosquito (TTP Labtech). Reverse transcription was carried out by incubating wells on a ProFlex 2x384 thermal-cycler (Thermo Fisher) at 42°C for 90 min and stopped by heating at 70°C for 5 min.

60Subsequently, 1.5 μL of PCR mix (1.67X KAPA HiFi HotStart ReadyMix (Kapa Biosystems, KK2602), 0.17 μM IS PCR primer (IDT, 5′-AAGCAGTGGTATCAACGCAGAGT-3′), and 0.038 U/μl Lambda Exonuclease (NEB, M0262L)) was added to each well with a Mantis liquid handler (Formulatrix) or Mosquito, and second strand synthesis was performed on a ProFlex 2x384 thermal-cycler by using the following program: 1. 37°C for 30 minutes, 2. 95°C for 3 minutes, 3. 23 cycles of 98°C for 20 s, 67°C for 15 s, and 72°C for 4 minutes, and 4. 72°C for 5 minutes.

61The amplified product was diluted with a ratio of 1-part cDNA to 10 parts 10mM Tris-HCl (Thermo Fisher, 15568025). 0.6 μL of diluted product was transferred to a new 384-well plate using the Viaflow 384 channel pipette (Integra). Illumina sequencing libraries were prepared as described in. Briefly, tagmentation was carried out on double-stranded cDNA using the Nextera XT Library Sample Preparation kit (Illumina, FC-131-1096). Each well was mixed with 0.8 μL Nextera tagmentation DNA buffer (Illumina) and 0.4 μL Tn5 enzyme (Illumina), then incubated at 55°C for 10 min. The reaction was stopped by adding 0.4 μL “Neutralize Tagment Buffer” (Illumina) and spinning at room temperature in a centrifuge at 3220xg for 5 min. Indexing PCR reactions were performed by adding 0.4 μL of 5 μM i5 indexing primer, 0.4 μL of 5 μM i7 indexing primer, and 1.2 μL of Nextera NPM mix (Illumina). All reagents were dispensed with the Mantis or Mosquito liquid handlers. PCR amplification was carried out on a ProFlex 2x384 thermal cycler using the following program: 1. 72°C for 3 minutes, 2. 95°C for 30 s, 3. 12 cycles of 95°C for 10 s, 55°C for 30 s, and 72°C for 1 minute, and 4. 72°C for 5 minutes.

62Following library preparation, wells of each library plate were pooled using a Mosquito liquid handler. Pooling was followed by two purifications using 0.7x AMPure beads (Fisher, A63881). Library quality was assessed using high sensitivity capillary electrophoresis on a Fragment Analyzer (Agilent) or Tapestation (Agilent), and libraries were quantified by qPCR (Kapa Biosystems, KK4923) on a CFX96 Touch Real-Time PCR Detection System (Biorad). Plate pools were normalized to 2 nM and equal volumes from library plates were mixed together to make the sequencing sample pool.

63Libraries were sequenced on the NextSeq or NovaSeq 6000 Sequencing System (Illumina) using 2 × 100bp paired-end reads and 2 × 8bp or 2 × 12bp index reads. NextSeq runs used high output kits, whereas NovaSeq runs used either a 200 or 300-cycle kit (Illumina, 20012860). PhiX control library was spiked in at ∼1%.

64All specimens were acquired from individuals with NSCLC as noted above. 4-micron thick formalin-fixed paraffin embedded (FFPE) human tissue sections were processed using previously published method and Agilent-Dako manufacturer recommendations were followed for antigen retrieval. All wash steps were performed at room temperature for three minutes each, unless otherwise noted. Briefly, slides were deparaffinized in xylenes (2 washes, 5 min each), and rehydrated in graded dilutions of aqueous ethanol (2 washes in 100% EtOH; 2 washes in 95% EtOH; 1 wash in 70% EtOH). Slides were washed once in ddH2O before being placed in an antigen target retrieval solution, 1x pH 6.1 Citrate retrieval solution (Dako) and pressure cooked using one cycle (2 hours) for antigen retrieval. Slides were allowed to cool to room temperature, washed three times with 1x PBS, then the tissue was blocked for endogenous peroxidase activity for 10 minutes using 0.3% H2O2. Slides were washed three times with 1x PBS, then incubated for 1 hour in a prepared protein blocking buffer solution (1X PBS containing 1% (w/v) BSA, 15 mM sodium azide, 0.05% (w/v) Tween-20). Slides were incubated overnight at 4C with either β-catenin (CTNNB1) rabbit monoclonal antibody (#8480S, Cell Signaling Technology, 1:100 dilution), AQP4 rabbit monoclonal antibody (#59678, Cell Signaling Technology, 1:100 dilution), or SUSD2 rabbit polyclonal antibody (HPA004117, Sigma Aldrich, 1:400 dilution). The following morning, the slides were washed three times with 1x PBS, incubated using commercial anti-rabbit and anti-mouse labeled polymer-HRP solution (Agilent Dako) for 30 minutes. Slides were washed three times in 1x PBS before incubation with freshly prepared 3,3-diaminobenzidine chromogen solution (Agilent Dako) for < 1 minute. Slides were washed twice in ddH2O and were counterstained using a commercial hematoxylin solution (VWR Biosciences). Excess dye was removed using three washes in ddH2O, and the hematoxylin was developed by incubating for 1 minute in 0.1% (w/v) sodium bicarbonate solution, and washed once in ddH2O. Tissues were dehydrated in aqueous ethanol (2 washes in 95% EtOH; 2 washes in 100% EtOH) and incubated in xylene for 5 minutes before being coverslipped. Stained slides were digitized using an Aperio ScanScope XT Slide Scanner (Leica Biosystems) using a 40X objective.

65Multiplex immunofluorescence staining was performed on the Opal IHC Multiplex Assay (NEL821001KT, Akoya Biosciences). Sequential 4 micron sections mounted on glass slides were sequentially stained for panel 1: PD-L1 (clone E1L3N, dilution 1:50, Cell Signaling Technologies), CD68 (clone KP1, dilution 1:500, Dako), IDO (clone D5J4E, dilution 1:100, Cell Signaling Technologies), HLA-DR (clone CR3/43, dilution 1:250, Abcam), CD14 (clone SP192, dilution 1:100, Abcam), and cytokeratin (polyclonal Z0622, dilution 1:250, Dako); or panel 2: CD3 (clone LN10, Leica), PD-1 (clone NAT105, dilution 1:100, Abcam), CD14, CD8 (clone C8/144B, dilution 1:100, Dako), FoxP3 (clone 236A/E7, dilution 1:200, Abcam), and cytokeratin on a Bond RX autostainer (Leica Biosystems). Slides were dewaxed (Leica), heat treated in Epitope retrieval solution 1 or 2 (AR9961/AR9640, Lecia) buffer depending on the antibody for 20 min at 93C, blocked in Antibody (Ab) Diluent (ARD1001EA, Akoya Biosciences), incubated for 30 min with the primary Ab, 10 min with horseradish peroxidase (HRP)-conjugated secondary polymer (anti-rabbit and anti-mouse, ARH1001EA, Akoya Biosciences), and 10 min with HRP-reactive OPAL fluorescent reagents (NEL821001KT, Akoya Biosciences). Slides were washed between staining steps with Bond Wash (AR9590, Leica) and stripped between each round of staining with heat treatment in antigen retrieval buffer. After the final heat treatment in antigen retrieval buffer, the slides were stained with spectral DAPI (FP1490, Akoya Biosciences), and coverslipped with Prolong Diamond mounting media (P36961, Thermo Fisher). Whole slide scans were acquired using the 10x objective via the Vectra imaging system (Perkin Elmer, version 3.0).

66For validation of candidate gene expression via a RT2 Profiler PCR array (QIAGEN, CLAH34795), human lung cancer PC9 cells (5 × 105) were treated for 48 hours (day 2) with DMSO (TN) or for 7 and 19 days with 2μM Osimertinib (Selleck Chemicals, S7297) with replenishment of drug every 3-4 days (Persister cells that evade drug-induced apoptosis by being in a low- to no-proliferative state, in patients this corresponds to the RD state), respectively. PD samples were derived from an acquired resistant PC9 cell line (Osimertnib IC50 = 89μM), that was generated by continuous treatment with 2μM Osimertinib with replenishment of drug every 3-4 days and presented active proliferation under drug at which time they were considered to be resistant and in the PD state. RNA was extracted via RNeasy Mini Kit (QIAGEN, 74104). RNA quality was confirmed as RIN > 7.5 via Bioanalyzer RNA 6000 Pico kit (Agilent, 5067-1514) and RNA was quantified via Qubit RNA HS Assay kit (Thermo Fisher Scientific, Q32852). A total of 400ng of RNA was reverse transcribed using the First Strand Synthesis Kit (QIAGEN, 330401) and then loaded into a custom 384 well RT2 profiler array (QIAGEN, CLAH34795).

67Small molecule inhibitors were all purchased commercially from Selleck Chemicals, and included Osimertinib (S7297), Alectinib (S2762), XAV-939 (S1180), and PRI-724 (S8262). Dimethyl sulfoxide (DMSO) (Fisher Scientific) was used to dissolve small molecule inhibitors according to manufacturer’s recommendations for use in in vitro experiments. PC9 and H3122 cells (5 × 103) were seeded in 96-well plate format (μclear CellStar, Greiner) and rested for 24 hours before treatment. Treatment included: i) DMSO, ii) tyrosine kinase inhibitors (TKI) Osimertinib (PC9 cells) or Alectinib (H3122 cells), iii) Wnt/β-catenin inhibitors PRI-724 or XAV-939, and v) indicated combination therapies of TKI and Wnt/ β-catenin inhibitors. All conditions were plated in technical quadruplicate and cells were retreated every 3 days. At each imaging interval, cellular nuclei were stained with Hoechst 33342 (Thermo Fisher Scientific) and scanned using a CellInsight High-Content Microscope (Version 6.4.3 Build 7204, Thermo Fisher Scientific) with a 4X objective.

68Sequences from the Illumina sequencing were demultiplexed using bcl2fastq version 2.19.0.316 (Illumina). Reads were aligned using the hg38 genome using STAR version 2.5.2b with parameters TK. Gene counts were produced using HTSEQ version 0.6.1p1 with default parameters except stranded was set to false and mode was set to intersection-nonempty.

69Standard procedures for filtering, variable gene selection, dimensionality reduction, and clustering were performed using the Seurat v3 in RStudio using R, where cells with fewer than 500 genes and 50,000 reads were excluded. We used DoubletFinder to identify potentially sorted doublet cells. 218 doublets were excluded from further analysis. Samples with less than 10 total cells were filtered from the analysis. Counts were log-normalized, then scaled by linear regression against the number of reads. Variable genes (Ngenes = 2,000) were selected using a threshold for dispersion, with z-scores normalized by expression level. The variable genes were projected onto a low-dimensional subspace using principal component analysis. The number of principal components (Npcs) were selected based on inspection of the plot of variance explained (Npcs = 20). A shared-nearest-neighbors graph was constructed based with metric the Euclidean distance in the low-dimensional subspace. Cells were visualized using a 2-dimensional tSNE on the same distance metric (Res = 0.5, Kparam = 30, script 03). Cell types were assigned to each cluster of cells using the abundance of known marker genes (Table S2, script S01-03 and script NI01).

70Cells previously annotated as epithelial (n = 5,581) were subset and re-clustered using methods described above and the following parameters: Ngenes = (2,000), Npcs = 20, Res = 0.7, Kparam = 30 (script NI02). Malignant epithelial cells were identified using inferCNV. inferCNV which works by finding cells with large copy number variations as determined by sorting expressed genes by their chromosomal location and applying a moving average, a sliding window of 100 genes within each chromosome, to the relative expression values. All epithelial cells as well as 300 fibroblasts and 300 endothelial cells were used as input (script NI03). An additional 500 fibroblasts and 500 endothelial cells were used as reference controls. We scored each cell for the extent of CNV signal and plotted cells on a dendrogram which was then cut at the highest point in which all the spiked in endothelial and fibroblasts cells belonged to one cluster (k = 6, one fibroblast control was misassigned). All cells that clustered together with spiked in controls were labeled “nontumor,” whereas the remaining two clusters were labeled as “tumor.”

71Noncancerous epithelial cells (n = 1,827), as determined as those cells lacking large chromosomal aberrations from InferCNV analysis, were subset and re-clustered using the following parameters: Ngenes = (2,000), Npcs = 20, Res = 0.5, Kparam = 20 (script NI05). Cell types were assigned to each cluster of cells using the abundance of known marker genes (Table S2) and differentially expressed genes as found by using the Seurat function FindAllMarkers using the default Wilcoxon rank sum test.

72Cancerous epithelial cells (n = 3,754), as determined as those cells harboring large chromosomal aberrations from InferCNV analysis, were subset and re-clustered using the following parameters: Ngenes = 2,000, Npcs = 20, Res = 0.9, Kparam = 10 (script NI04). We found the differences in gene expression between the three treatment time points (TN, RD, and PD) using the Seurat function FindMarkers using the MAST test and sample_name as the latent variable. Three separate tests were used to ascertain the differences between: 1) TN and RD, 2) TN and PD and 3) RD and PD (Table S5). Resulting differential gene lists were then filtered to limit patient specific effects. This is achieved by setting a threshold for non-zero expressing cells per patient (RD = 3 of RD patients and PD = 6 of PD patients) and removing differentially expressed genes explained by less than the thresholds set. The top 100 genes from each comparison were manually curated to evaluate for pathway activation. Decreased expression could indicate lack of detection due to the stochasticity of scRNaseq and thus for analysis of activated pathways we focused on upregulated genes. Gene signatures (Table S2) were compiled using differential expressed as well as known cell marker genes. Specifically, the alveolar signature is made of differentially expressed AT1/AT2 genes among the cancer cell time point comparisons as well has additional known AT1/AT2 genes. The remaining signatures were identified directly from top differentially expressed genes.

73To ensure that we were not misclassifying healthy AT2 cells as cancer cells, we compared the expression levels of our combined alveolar gene signature between the three time points (TN, RD, PD) and non-cancer AT2 cells from our dataset as well as additional non-cancer AT2 cells from an external dataset. Non-cancer AT2 cells from our dataset were more similar to the external AT2 cells than any of our cancer cells across all time points (average spearman correlation coefficient = 0.65, −0.10, 0.24, −0.19, for non-cancer AT2 cells, and TN, RD, PD cancer cells respectively).

74Cancer cells from EGFR and ALK driven tumor samples were subset separately. We compared all three treatment time points (TN, RD, PD) for EGFR patients where we only compared two treatment time points for ALK (TN, PD) as only one ALK+ driven sample represented the RD time point. We then compared the five cancer cell signatures derived from the grouped analysis (alveolar, kynurenine, plasminogen activation, SERPINE1, and gap junction). Pairwise wilcoxon tests were calculated between each treatment time point (TN, RD, and PD).

75To understand the PD sample heterogeneity all cancer cells from PD samples were subset. Each sample’s average expression of genes included in gene signatures (alveolar, kynurenine, plasminogen activation, SERPINE1, and gap junction, Figures S3B and S3F–S3H) and overall signature score was calculated and plotted using the R pheatmap package

76Longitudinal analysis of a single patient was done by subsetting all cells originating from patient TH226. As above, the differences in gene expression between the three treatment time points (TN, RD, and PD) was found by applying the Seurat function FindMarkers using the MAST test with sample_name as the latent variable. Three separate tests were used to ascertain the differences between: 1) TN and RD, 2) TN and PD and 3) RD and PD (script NI07-08, Table S4)

77TCGA LUAD data were downloaded from https://xenabrowser.net/datapages/. Metadata was downloaded from An Integrated TCGA Pan-Cancer Clinical Data Resource. Mean expression of each cancer cell expression signature (alveolar, kynurenine, plasminogen activating, SERPINE1, and gap junction) was calculated per TCGA sample. TCGA samples were then split by quartile groups. Quartiles were plotted using library packages survival and survminer in R (script NI10). Log rank p values are reported for each signature across four expression quartiles. Cox hazard regression model was computed for comparison of quartile 1 (low expressors) versus quartile 4 (high expressors) for all signatures.

78Tumor populations were annotated, then scored in a blinded, randomized analysis by a clinical pathologist for percent tumor positivity and subcellular staining intensity at the membrane, cytosolic, and nuclear compartments. SUSD2 membrane staining was graded by two reviewers in a blinded, randomized fashion using the slides annotated for tumor presence. Staining intensity was graded as negative, weak, intermediate, or strong and received scores of 0, 1, 2, or 3 respectively. Percent tumor positivity coefficient was graded as 0, negative; 1, less than 10% immunopositive; 2, between 10%–50% immunopositive; 3, between 51%–80% immunopositive; 4, greater than 80% immunopositive. Calculation of immunoreactivity scores was performed by multiplying the staining intensity score (0-3) with the percent tumor positive coefficient (0-4) to yield a value between 0 and 12.

79Alignment bams for all non-immune cells (stroma and epithelial) were passed to GATK HaplotypeCaller which was run from the latest available Docker container (broadinstitute/gatk:4.0.11.0) using the following options:–disable-read-filter MappingQualityReadFilter–disable-read-filter GoodCigarReadFilter–disable-read-filter NotSecondaryAlignmentReadFilter–disable-read-filter MappedReadFilter–disable-read-filter MappingQualityAvailableReadFilter–disable-read-filter NonZeroReferenceLengthAlignmentReadFilter–disable-read-filter NotDuplicateReadFilter–disable-read-filter PassesVendorQualityCheckReadFilter–disable-read-filter WellformedReadFilter

80–disable-read-filter MappingQualityReadFilter

81–disable-read-filter GoodCigarReadFilter

82–disable-read-filter NotSecondaryAlignmentReadFilter

83–disable-read-filter MappedReadFilter

84–disable-read-filter MappingQualityAvailableReadFilter

85–disable-read-filter NonZeroReferenceLengthAlignmentReadFilter

86–disable-read-filter NotDuplicateReadFilter

87–disable-read-filter PassesVendorQualityCheckReadFilter

88–disable-read-filter WellformedReadFilter

89Disabling these specific read filters proved necessary for scRNaseq, as inherent low-coverage causes the vast majority of reads to be flagged for removal otherwise. The full human variant set (dbSNP) was downloaded from NCBI (https://www.ncbi.nlm.nih.gov/variation/docs/human_variation_vcf/), and every variant call was assessed for its presence/absence in the human variant database. dbSNP is a public, living catalog of ∼674 million human somatic SNPs and indels that have been reported by peer-reviewed publications.

90Cloud-based parallelization of HaplotypeCaller jobs was achieved with Reflow, a workflow engine for distributed, incremental data processing in the cloud (GRAIL, https://github.com/grailbio/reflow). HaplotypeCaller outputs a separate variant calling format file (VCF) for each cell, which were processed with the python package cerebra (https://github.com/czbiohub/cerebra). Variants found in dbSNP were removed, not to be included in further analysis. We reasoned that by removing ‘common’, population-level variants, we could better hone-in on disease specific variation.

91In addition to scRNaseq reads, we obtained bulk DNA reads from peripheral blood for the majority of our patients (with the exception of three). These PBMC reads were run through HaplotypeCaller to establish ‘germline’ mutation profiles for each of our patients. Germline mutations were then subtracted out from each of that patient’s single cell VCFs. This filtering step was omitted for the three patients for which we did not obtain peripheral blood, however, these single cell VCFs were still passed through our dbSNP filter.

92We also applied a fathmm filter to all cells. fathmm takes a machine learning approach to predict the likelihood of a given SNP to be pathogenic, integrating ENCODE annotations for things like transcription factor binding sites, histone modifications, cross-species sequence alignment and conservation scores, etc. Only variants computationally predicted to be pathogenic were included in our analysis, i.e., those variants with a fathmm score > 0.7.

93The remaining variants were then filtered through the COSMIC (Catalogue of Somatic Mutations in Cancer) complete mutation–genome screens database (https://cancer.sanger.ac.uk/cosmic/download). Only SNPs/indels associated with ‘Lung’ as per their COSMIC annotation were kept. Variant calls were mapped to their corresponding genes, and per-patient / per-sample mutational profiles were established. We used the ERCCs spiked into each cell sample as a negative control for false positive mutations, which can arise due to technical artifacts such as PCR errors. We found the median false positive mutation rate to be 0.000256% per.

94Fusion transcripts were detected with STAR-fusion (https://github.com/STAR-Fusion/STAR-Fusion/wiki) version 1.6.0, run from a Docker container (trinityctat/ctatfusion:1.5.0). The following options were used: –FusionInspector validate, –examine_coding_effect, –denovo_reconstruct. Distributed processing of STAR-fusion jobs was accomplished with Reflow. Output files were processed with cerebra, then combined with variant calls to create per-cell and per-sample summary tables.

95Mutation information from cerebra outputs were summarized by sample. Coverage information was provided by a secondary output from cerebra summarized by sample and gene. Where all cells are summarized by sample and all fathmm filtered ROIs are summarized by corresponding gene (script NI06). Plots were generated using the R pheatmap package. Two comprehensive tables, Tables S3 and S7, detail mutations and fusions per cell.

96MSK-Impact data was downloaded from cBioPortal (and subset to only NSCLC samples MSK-Impact data was subset to only those mutations that were also found in the scRNaseq dataset of mutations (n = 141 unique mutations)). We stratified MSK-Impact samples by those with greater than or equal to 2 mutations from the tier one COSMIC mutations found in the scRNaseq dataset (mutation high), and those less than 2 mutations (mutation low) (Figure 2D). Kaplan-Meier plots were visualized with the lifelines package in python version 3.4 (script NI12).

97All cells annotated as immune (n = 13,431) were subset and clustered as described above (script IM01) using the following parameters (Ngenes = 2000, Npc = 20, Res = 0.7). The resulting 18 clusters were assigned to different major immune cells types using a list of curated gene markers (Table S2) and by manual curation of differentially expressed genes for each cluster (Table S4). The different cell types and number of cells belonging to each type are described in the main text.

98To assess changes in fractional abundance of different immune cell populations we used all cells though excluded thoracentesis and brain samples due to difference in the immune makeup of these tumor environments which would skew the data. The function freqCI from the R package REdaS (script IM02) was used to calculate confidence intervals for relative frequencies.

99Macrophages (n = 1,379) and T cells (n = 2,226) from lung biopsies were subset and clustered as described above (script IM03 and IM04 respectively) using the following parameters for MFs (Ngenes = 2000, Npc = 10, res = 0.3) and T cells (Ngenes = 2000, Npc = 10, res = 0.3). The resulting clusters are discussed in the main text and the lists of differentially expressed genes are provided (Table S4). We repeated this analysis where we subset the data to only patients with multiple biopsies and sufficient cells (TH226 and TH266) (script IM05).

100Three to six regions from each slide containing tumor and stroma were selected utilizing Phenochart (v1.0.8, Perkin Elmer) for high resolution multispectral acquisition on the Vectra system at 20X magnification. The images were analyzed with inForm software (v2.4.8, Akoya) to unmix adjacent fluorochromes, subtract autofluorescence, segment the tissue into tumor and stroma regions, segment the cells into nuclear, cytoplasmic, and membrane compartments, and to phenotype the cells according to morphology and cell marker expression. Fractions of macrophage and T cell populations were calculated as: (population of interest) / (macrophage + T cell populations) and plotted using ‘ggplot2′ in R.

101As with the survival analysis using cancer cell gene signatures, we used the downloaded TCGA LUAD dataset and metadata to access patient survival outcomes as they pertain to the fractional changes of immune populations within a given tumor. We used CIBERSORT to deconvolute the bulk TCGA samples into relative fractions of immune cell populations as determined by using the LM22 reference. The total macrophage population was found by combining fractions for Monocytes, Macrophages.M0, Macrophages.M1, and Macrophages.M2. The total T cell population was found by combining fractions of T.cells.CD8, T.cells.CD4.naive, T.cells.CD4.memory.resting, T.cells.CD4.memory.activated, T.cells.follicular.helper, T.cells.regulatory.Tregs, T.cells.gamma.delta, NK.cells.resting, and NK.cells.activated. TCGA samples were then split by quartile groups. Quantiles were plotted using library packages survival and survminer in R (script NI10). Log rank p values are reported across four expression quartiles. Cox hazard regression model was computed for comparison of quartile 1 (low expressors) versus quartile 4 (high expressors).

102Fold Change was calculated by determining the ratio of mRNA levels to control (day 2) values using the delta threshold cycle (Ct) method (DCt). A t test was used to find the significance of change between baseline (day 2) and treated time points (days 7, 19 and 70) based on normalized Cts to baseline (script NI14). Plots were made using ‘ggplot2′ in R.

103Analysis was performed using CellInsight (Thermo Fisher Scientific) companion software across technical quadruplicates. The DMSO treated condition and single agent Wnt/β-catenin inhibitors reached confluency after 3 days in culture (100% maximum cut-off value). Significance values were calculated using a Student’s t test calculated at treatment endpoints (day 6).

104Detailed protocols for single cell dissociation of small tumor biopsies (https://doi.org/10.17504/protocols.io.65rhg56) and high throughput smartseq2 libraries (https://doi.org/10.17504/protocols.io.2uwgexe) are available at protocols.io.