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

Abstract

Lung cancer, the leading cause of cancer mortality, exhibits heterogeneity that enables adaptability, limits therapeutic success, and remains incompletely understood. Single-cell RNA sequencing (scRNA-seq) of metastatic lung cancer was performed using 49 clinical biopsies obtained from 30 patients before and during targeted therapy. Over 20,000 cancer and tumor microenvironment (TME) single-cell profiles exposed a rich and dynamic tumor ecosystem. scRNA-seq of cancer cells illuminated targetable oncogenes beyond those detected clinically. Cancer cells surviving therapy as residual disease (RD) expressed an alveolar-regenerative cell signature suggesting a therapy-induced primitive cell-state transition, whereas those present at on-therapy progressive disease (PD) upregulated kynurenine, plasminogen, and gap-junction pathways. Active T-lymphocytes and decreased macrophages were present at RD and immunosuppressive cell states characterized PD. Biological features revealed by scRNA-seq were biomarkers of clinical outcomes in independent cohorts. This study highlights how therapy-induced adaptation of the multi-cellular ecosystem of metastatic cancer shapes clinical outcomes.

Keywords: ALK; EGFR; lung cancer; single-cell RNA sequencing; targeted therapy.

Graphical abstract

Figure 1

Patient Characteristics and Experimental Overview (A) Consort diagram. 56 biopsies were processed, 49 samples passed quality control. (B) Tissue processing pipeline for scRNA-seq. Patient samples were disaggregated into single cells and sorted into microtiter plates using FACS. cDNA synthesis was performed using the Smart-seq2 protocol, and libraries were sequenced on Illumina platforms. (C) Circle plot of the clinically identified oncogenic driver (outer circle) and treatment time point (inner circle) for each sample. (D) t-stochastic neighbor embedding (t-SNE) plot of all cells colored by their cellular identity (epithelial cells [n = 5,581], immune cells [n = 13,431], stromal cells [n = 4,249]). See also Figure S1 and Tables S1 and S2.

Figure S1

Related to Figure 1 (A) Bar plots of the time interval between treatment start and tissue acquisition for PD and RD tumor samples (B) t-SNE of all epithelial cells (n = 5,581), numbers correspond to individual clusters. (C) Inferred large-scale copy number variations (CNVs) help identify cancer (pink) and non-cancer cells (purple). Epithelial and spike in control cells are included in the x axis and chromosomal regions on the y axis. Amplifications (red) or deletions (blue) were inferred by averaging expression over 100-gene stretches on the respective chromosomes. (D) Bar plot of cell counts for annotated epithelial cells. (E) Bar plot of the number of unique genes across all annotated epithelial cell types. (F) Bar plot of unique gene count of cancer versus non-cancer epithelial cells.

Figure S2

Related to Figure 2 (A) t-SNE plot of 3,754 cancer cells from 44 samples, numbers indicate individual clusters. (B) Circle plot illustrating the clinically identified oncogenic driver (outer circle) and time points (inner circle) of each biopsy, only for cancer cells. C) Density distribution of cluster occupancy of cancer (red) and non-cancer (blue) epithelial cell clusters, calculated as the percentage of the highest contributing individual patient over the total number of cells for that cluster. (D) tSNE of cancer cells colored by patient. (E) tSNE of non-cancer epithelial cells colored by patient. (F) Illustration of heterogeneity of primary driver mutated cancer cells found in exemplary sample LTS47. (G) Clinical characteristics of the 44 NSCLC samples in which at least one cancer cell was identified. Columns indicate clinically identified mutated gene, treatment response time point (TN, RD, PD), biopsy site, and primary or metastatic sample origin, respectively. (H) Cancer cell mutational landscape for each patient sample as determined by scRNaseq represented as a heatmap. Color indicates the number of mutant reads for each genomic region and sample divided by the total number of reads for that region in that sample, NC:No Coverage over the specific genomic region. (I) Mutational landscape of COSMIC tier 1 genes. Color indicates the number of mutant reads for each genomic region and sample divided by the total number of reads for that region in that sample, NC:No Coverage over the specific genomic region. (J) Kaplan-Meier plot showing overall survival of 1269 NSCLC patients within the MSK-Impact dataset. Patients were stratified by high (> = 2) and low (< 2) mutations from the 141 mutations that are present in both the MSK-Impact dataset and panel (I).

Figure 2

scRNA-seq Infers Patient Mutational Status and Reveals a Complex Mutational Landscape in Cancer Cells (A) Clinical characteristics of the 44 NSCLC samples in which at least one cancer cell was identified. Columns indicate clinically identified mutated gene, treatment response time point (TN, RD, PD), biopsy site, and primary or metastatic sample origin, respectively. (B and C) Cancer cell mutational landscape for each patient sample as determined by scRNA-seq represented as a binarized heatmap across driver genes (B) and COMSIC tier 1 genes (C). Red indicates the presence of mutation while blue indicates that no mutation was identified. See also Figure S2 and Tables S2 and S3.

Figure S3

Related to Figure 3 (A) Dot plot of the relative expression of established cellular proliferation genes (x axis) across treatment time points (y axis). The color intensity scale reflects the average gene expression and the size scale indicates the number of cells expressing the gene within that treatment time point. Applying grouped, pairwise comparisons of treatment time points of the average scaled expression of all genes demonstrated significantly different expression (p < 0.0001) in all comparisons. (B) Heatmap showing the expression of genes in the alveolar signature. Cells are grouped by treatment time point. (C) Boxplot of Spearman correlations of cancer cells from all treatment time points and healthy AT2 cells to an external reference of healthy AT2 cells. Non-cancer AT2 cells from our dataset were more similar to the external, healthy AT2 cells than any of our cancer cells across all time points (mean ρ = 0.65, −0.10, 0.24, −0.19, for healthy AT2 cells, and TN, RD, PD cancer cells, respectively). ^∗∗∗ indicates a p value < 0.001 (D-F) Immunoreactivity score (IRS) for membrane AQP4 (D), membrane SUSD2 (E), and nuclear CTNNB1 (F) across all time points. (G) Pairwise comparison of nuclear CTNNB1 IRS for a subset of patients receiving neoadjuvant TKI treatment prior to surgical removal of tumors, allowing for controlled, matched sample pairs at TN and RD treatment time points. Samples with AZ identifiers refer to patients with EGFR mutant NSCLC receiving neo-adjuvant osimertinib treatment. Sample with NC identifier refer to patient with ROS1 fusion-positive NSCLC receiving neo-adjuvant crizotinib treatment. (H-O) High content microscopy screening of EGFR mutant PC9 cells and ALK fusion-positive H3122 cells showing treatment response to TKI in presence or absence of WNT/β-catenin inhibition. In comparison to DMSO control, upper panel (H-K) shows single agent treatment, lower panel (L-O) shows combinational treatment of TKI WNT/β-catenin inhibitors. Two WNT/β-catenin inhibitors have been tested, XAV-939 (H,I,L,M) and PRI-724 (J,K,N,O). Values are shown as percent confluency, with maximum cutoff for full well confluency (100%). p values are calculated for all end points (day 6) values compared to single agent TKI. Error bars represent mean ± standard error of the mean (SEM), n = 4 technical replicates. (P-R) Heatmaps showing the expression of genes within each signature (kynurenine, SERPINE1/plasminogen activation, and gap junction, respectively) grouped by treatment time point.

Figure 3

Differential Gene-Expression Analysis between Treatment Time Points Reveals Treatment Stage-Specific Transcriptional Signatures (A) Boxplots showing the expression level of the alveolar signature across treatment time points as well as non-cancerous AT2 cells from our cohort. ∗∗∗p < 0.001. (B) Fold-change expression of NKX2-1 as quantified by RT-PCR in EGFR mutant PC9 cells after specified treatment duration (see STAR Methods), ^∗∗∗p < 0.001. (C) Representative IHC images of TN, RD, and PD tumor tissue sections stained for AQP4 demonstrating increased expression at the RD time point. Red arrows indicate cancer cells of interest. Scale bars correspond to 50 μm. (D) Kaplan-Meier plot of the relationship between the alveolar signature and patient OS within the TGCA dataset. Patients were stratified by signature expression quartile (Q1 = 128, Q2 = 127, Q3 = 128, Q4 = 127), where Q1 is the lowest expression and Q4 is the highest expression. (E and F) Representative IHC images of TN, RD, and PD tumor tissue sections stained for SUSD2 (E) and CTNNB1 (F) demonstrating increased expression at the RD time point. Red arrows indicate example regions of interest. Scale bars correspond to 50 μm. (G–J) Treatment response upon inhibition of β-catenin activity in EGFR mutant PC9 cells and ALK fusion-positive H3122 cells. Relative viability is shown as percent confluency compared to DMSO control. PC9 cells were treated with XAV-939 (G) or PRI-724 (H) with or without the combination of 50nM osimertinib. H3122 were treated with XAV-939 (I) or PRI-724 (J) with or without the combination of 50 nM alectinib. p values were calculated for all end points (day 6) values compared to single agent TKI. Error bars represent mean ± standard error of the mean (SEM), n = 4 technical replicates. ^∗p < 0.05, ^∗∗p < 0.01. (K) Boxplots showing the expression levels of the kynurenine signature expression across different treatment time points. ∗∗∗p < 0.001. (L) Fold-change expression of QPRT as quantified by RT-PCR in PC9 cells after treatment with osimertinib as in (B) (see STAR Methods) (AR), ^∗p < 0.05. (M) Kaplan-Meier plot of the relationship between the kynurenine signature and patient OS within the TGCA dataset. Patients were stratified by signature expression quartile (Q1 = 128, Q2 = 127, Q3 = 128, Q4 = 127), where Q1 is the lowest expression and Q4 is the highest expression. (N) Boxplot showing the expression levels of the plasminogen activation pathway signature across different treatment time points. (O) Boxplot showing the expression levels of the SERPINE1 across different treatment time points. (P) Kaplan-Meier plot of the relationship between the plasminogen activating pathway signature and patient OS within the TGCA dataset, respectively. Patients were stratified by signature expression quartile (Q1 = 128, Q2 = 127, Q3 = 128, Q4 = 127), where Q1 is the lowest expression and Q4 is the highest expression. (Q) Kaplan-Meier plot of the relationship between SERPINE1 expression and patient OS within the TGCA dataset, respectively. Patients were stratified by signature expression quartile (Q1 = 128, Q2 = 127, Q3 = 128, Q4 = 127), where Q1 is the lowest expression and Q4 is the highest expression. (R) Boxplot showing the expression levels of the gap-junction signature across treatment time points. (S) Kaplan-Meier plots of relationships between the gap-junction signature and patient OS within the TGCA dataset. Patients were stratified by signature expression quartile (Q1 = 128, Q2 = 127, Q3 = 128, Q4 = 127), where Q1 is the lowest expression and Q4 is the highest expression. See also Figures S3, S4, and S5 and Tables S2, S4, S5, and S6.

Figure S4

Related to Figure 3 (A) Graphical summary of cancer cell expression changes across treatment time points. RD features include (1) Alveolar signature, and (2) various RD specific invasive signaling pathways. PD features include: (3) kynurenine signature, (4) plasminogen activation and SERPINE1 signatures, (5) gap junction proteins, (6) expression of pro-inflammatory chemokines, (7) loss of tumor suppressor genes, and (8) various PD specific invasive signaling pathways. (B-F) Boxplots of pathway signature changes (alveolar, kynurenine, plasminogen activating, SERPINE1, and gap junction, respectively) across treatment time points within only EGFR mutant cancer cells (^∗∗∗ indicates a p value < 0.0001). (G-K) Boxplots of pathway signature changes (alveolar, kynurenine, plasminogen activating, SERPINE1, and gap junction, respectively) across treatment time points within only ALK cancer cells (^∗∗∗ indicates a p value < 0.0001). (L-O) Heatmap of sample average expression with PD only cancer cells for each cancer derived signature gene (alveolar, kynurenine, plasminogen activating/SERPINE1, and gap junction, respectively).

Figure S5

Related to Figure 3 (A, B, C) Longitudinal timeline of patient treatment, (A) Chest CT scan at each clinical evaluation time point, (B) Biopsy time point with procedural CT scan, (C) Hematoxylin and eosin (H&E) staining from treatment naive and progression time points demonstrating adenocarcinoma and squamous cell carcinoma, respectively, scale bar indicates 50 μm. (D) Heatmap of mutation state in clinical driver and a subset of COSMIC tier 1 mutated genes (displayed COSMIC tier 1 mutations occur in at least two out of three samples). Color red indicates the presence of mutation whereas color blue indicates no presence of mutation. (E-I) Boxplots of pathway signature changes (alveolar, plasminogen activating, SERPINE1, gap junction and squamous histology, respectively) across treatment time points.

Figure 4

Changes in the Composition of the Tumor Microenvironment within Each Tumor (A) t-SNE plot of all immune cells colored by immune cell type. (B) Patient occupancy for each immune cell type. (C) Fractional changes for each immune cell type across the three treatment states. Error bars indicate the 95% confidence interval for the calculated relative frequencies. ^∗p < 0.01 using a chi-square test of independence. (D) Representative in situ immunofluorescence images of changes from TN to RD and TN to PD in tumor tissue sections from two separate samples at two separate time points; AZ003 (TN and RD), TH281 (TN and PD). Scale bars correspond to 50 μm. (E) Quantification of fractional changes of macrophages across treatment time points from the images in (D) and Figure S5F. (F) Quantification of fractional changes of T-cells across treatment time points from the images in (D) and Figure S5F. See also Figure S6.

Figure S6

Related to Figure 4 (Α) Pairwise Pearson correlations between each treatment group’s immune cell compositions which corresponds to the fraction of each immune cell type’s abundance in the total immune cell population. (B) Total immune cells for each biopsy of patient TH266. (C) Total immune cells for each biopsy of patient TH226. (D) Fraction of each immune cell subtype for the two biopsies of patient TH266. Error bars indicate the 95% confidence interval for the calculated relative frequencies. Asterisks next to the title of each cell type indicate significance (p < 0.01) when using a chi-square test of independence. Titles of non-significant cell types are colored red and lack an asterisk. (E) Fraction of each immune cell sub-type for the three biopsies of patient TH226. Error bars indicate the 95% confidence interval for the calculated relative frequencies. Asterisks next to the title of each cell type indicate significance (p < 0.01) when using a chi-square test of independence. Titles of non-significant cell types are colored red and lack an asterisk. (F) Representative in situ immunofluorescence images from two patients with matched samples at different treatment time points, demonstrating fractional changes in the immune populations of macrophages and T cells. Scale bars correspond to 50 microns. (G) Kaplan-Meier plot of deconvoluted TCGA lung adenocarcinoma data showing the relation between OS and the fraction of macrophages for each patient. Patients were stratified by high and low macrophage fraction.

Figure S7

Related to Figure 5 (Α) t-SNE plot of all lung-derived macrophage cells. (B) Heatmap showing the expression level of the top 10 differentially expressed genes for each macrophage cluster. (C) t-SNE plot of all lung-derived T cells. (D) Heatmap showing the expression level of the top 10 differentially expressed genes for each T cell cluster.

Figure 5

Immune Cell Subpopulations Demonstrate Unique Transcriptional Profiles within Each Treatment Time Point (A) Fraction of cells belonging to each treatment stage for each macrophage cluster in Figure S6. Error bars indicate the 95% confidence interval for the calculated relative frequencies. ^∗p < 0.01 using chi-square test of independence. (B) Violin plots showing the expression level distribution of notable individual genes. (C) Fraction of cells belonging to each treatment stage for each T cell cluster in Figure S6. Error bars indicate the 95% confidence interval for the calculated relative frequencies. ^∗p < 0.01 using chi-square test of independence. (D) Violin plots showing the expression level distribution of notable individual genes. (E) Graphical summary of immune microenvironment changes across treatment time points. See also Figure S7 and Tables S4 and S5.

Figure 6

scRNA-seq Profiles Reveal Clinical-State Specific Features of the Tumor Cellular Ecosystem Features common to all time points are shown in the top-left quadrant and include the presence of multiple oncogenic drivers (1). Features shared in RD and PD are shown in the top-right quadrant and include various invasive signaling pathways (2). Features unique at RD, shown in the bottom-right quadrant, include the Alveolar signature (3) and increased T cell fraction (4). Features unique to PD, shown in the bottom-left quadrant, include upregulation of the plasminogen activation pathway (5), expression of gap-junction proteins (6), loss of tumor suppressor genes (7), expression of pro-inflammatory chemokines (8), increased Treg population (9), and increased kynurenine signature expression (10).