Single-cell profiling of tumor heterogeneity and the microenvironment in advanced non-small cell lun

Lung cancer is a highly heterogeneous disease. Cancer cells and cells within the tumor microenvironment together determine disease progression, as well as response to or escape from treatment. To map the cell type-specific transcriptome landscape of cancer cells and their tumor microenvironment in advanced non-small cell lung cancer (NSCLC), we analyze 42 tissue biopsy samples from stage III/IV NSCLC patients by single cell RNA sequencing and present the large scale, single cell resolution profiles of advanced NSCLCs. In addition to cell types described in previous single cell studies of early stage lung cancer, we are able to identify rare cell types in tumors such as follicular dendritic cells and T helper 17 cells. Tumors from different patients display large heterogeneity in cellular composition, chromosomal structure, developmental trajectory, intercellular signaling network and phenotype dominance. Our study also reveals a correlation of tumor heterogeneity with tumor associated neutrophils, which might help to shed light on their function in NSCLC.

Tumor ecosystems are comprised of cancer cells, infiltrating immune cells, stromal cells, and other cell types together with noncellular tissue components, which interact and collectively determine disease progression as well as response to therapy1,2. It is well known that cancer patients elicit very individualized responses to different treatments, demanding better characterization of the whole tumor ecosystem beyond currently applied clinical typing of somatic mutations in cancer cells. Furthermore, precisely targeted therapies against well-defined oncogenic drivers reveal a wide spectrum of responses in different settings. For example, the KRAS G12C inhibitors seemed to induce tumor responses in the majority of lung cancers but much less in pancreatic cancers, which differ in their tumor micromilieu dominated by cancer-associated fibroblasts3. Antibodies targeting PD1 or PD-L1 have achieved substantial overall survival improvements in advanced non-small cell lung cancer (NSCLC). The 5-year survival rate could be prolonged from less than 5–29.6% in PD-L1-positive patients4,5. However, major challenges still remain, including low response rate in unselected patients, lack of reliable predictive biomarkers, and identification of more immunotherapeutic targets. Thus, comprehensive understanding of NSCLC ecosystems holds the promise to improve personalized treatment strategies6.

Conventional ‘bulk’ RNA-sequencing methods process a mixture of all cells, averaging out underlying differences in cell-type-specific transcriptomes. In contrast, single-cell RNA-sequencing (scRNA-seq) profiles the gene expression pattern of each individual cell and decodes its intercellular signaling networks. This unbiased characterization provides clear insights into the entire tumor ecosystem, such as mechanisms of intratumoral and intertumoral heterogeneity, as well as cell–cell interactions through ligand-receptor signaling7. Thus, several studies deeply characterized the lung tumor microenvironment (TME) at single-cell resolution. An extensive taxonomy of stromal cells with different pathway activities in NSCLC patients presented a first-ever lung cancer TME cell atlas8. Isolated infiltrating T cells in NSCLC were classified according to their functional states and dynamics and a subset of regulatory T cells (Tregs) was found to correlate with the poor prognosis in lung adenocarcinoma9. Tumor-infiltrating myeloid cells (TIMs), including monocyte, macrophage, dendritic, and granulocyte cell lineages, were categorized into at least 25 different states by scRNA-seq10. Subsets of TIMs defined by unique markers have been associated with patient prognosis. Heterogeneity of tumor endothelial cells was also studied for both human and mouse11. All reports mentioned above focused on early stage, resectable lung cancers, which may not reflect the cellular profiles of tumors at advanced stages that have undergone intense and exhaustive interactions with stromal and immune cells. Focusing on the evolutional dynamics of lung adenocarcinoma, Kim’s study was performed on the lung adenocarcinoma samples from early-stage tissues to advanced stage biopsies including both primary and metastatic sites12. Another recent study uncovered transcriptional signatures specific to various targeted therapies and clinical states on primary and metastatic lung biopsies by low throughput Smart-seq2 technology, which only included one squamous carcinoma patient13. Until now, the late-stage landscape of lung squamous carcinoma was mostly absent.

In this study, we apply scRNA-seq to analyze the cancer and TME landscape of advanced NSCLC for both lung adenocarcinoma and squamous carcinoma. We identify distinct cell populations and cellular signals that are differentially enriched in tumors depending on the pathological types, presence or absence of driver mutations, and degree of tumor heterogeneity. Our data provide a comprehensive scRNA-seq profiling on a large number of small biopsies and may be used to improve diagnostics and prognosis in clinical settings.

Establishment of advanced NSCLC cell atlas

We applied scRNA-seq analyses to biopsy samples from 42 advanced NSCLC patients with diverse histological and molecular phenotypes and treatment history (Fig. 1a, Supplementary Table 1). Following multiple quality control and filtering steps, a total of 90,406 cells were analyzed with respect to their transcriptomes. By characteristic canonical cell markers, eleven major cell types were detected, classified as carcinoma cell types, epithelial cells others than carcinoma cells, immune cell types (T cells, B lymphocytes, myeloid cells, neutrophils, mast cells, and follicular dendritic cells) and stromal cell types (fibroblasts and endothelial cells) (Fig. 1b, c and Fig. S1). Similar to the observations in previous studies, stromal and immune cells of different patients clustered together by cell types, while cancer cells showed higher heterogeneity and patient-specific expression signatures (Fig. 1d)14,15. Similar to the observations from previous studies8,12, the portions of cancer, stromal, and immune cells varied greatly among samples, which could be intrinsic to different tumor phenotypes or related to locations within the tumor where biopsies were taken (Fig. 1e and Supplementary Data 1). For example, tumor specimen P42 (lung adenocarcinoma mixed with sarcomatoid carcinoma) and P7 revealed a strongly inflammatory micromilieu with almost 50% T cells in contrast to specimen P2, P3, and P17, which were practically T cell depleted.

Lung Squamous Carcinoma has higher inter- and intratumor heterogeneity than lung adenocarcinoma

Based on single-cell expression levels of genes commonly used as markers for immunohistochemistry-based NSCLC classification, namely NAPSA, TTF-1 (NKX2-1) for lung adenocarcinoma (LUAD), and TP63, CK5 (KRT5) for lung squamous carcinoma (LUSC), subtype classifications aligned well with the histopathological classifications (Fig. S2). Next, we used the scRNA-seq data to infer copy number alterations (CNAs) in cancer cell populations. The inferred CNA profiles of 42 patients showed both interpatient and intrapatient heterogeneity (Fig. 2a). For LUAD patients, prominent arm-level insertions were found in chromosome 7 and 8q, with deletions in chromosome 10. Noteworthy, LUAD with known driver mutations have additional amplifications in the 1q and 5p arms. In contrast, LUSC patients mostly have 3q insertions and 5q deletions. Interestingly, some of the LUAD patients without driver mutations have similar CNA profiles to LUSC. Although expression profiles and composition of the cancer cell transcriptomes were largely patient-specific, carcinoma cells from some patients were more similar than others (Fig. 2band Fig. S3A, B). In most cases, cancer cells from LUAD and LUSC patients partitioned into separate clusters. More than half of the LUAD patients clustered into one group, while most LUSC tumors formed patient-specific clusters, indicating higher intertumor differences in LUSC than in LUAD. Most patients, especially patients with LUAD e.g., P16, P20, and P32, had dominant clones, while in a few LUSC such as P27 and P37 the malignant cells spread across multiple clusters (Fig. 2b and Fig. S3C). LUSC patients showed significantly higher clonality than LUAD patients (Fig. S3D).

To quantify the intratumoral heterogeneity, we defined both a CNA-based and an expression-based intratumor heterogeneity score, denoted as ITH_CNA and ITH_GEX (see Methods for their definitions). We observed various degrees of heterogeneity within the tumor (Fig. S4A, B). ITH_CNA and ITH_GEX showed a moderate correlation (Fig. 2d), potentially due to the nondriver genomic alternations or the microenvironment shaped tumor phenotypes. We further divided patients into three groups according to both the cancer type and mutation: LUAD patients with driver mutation (n = 12), denoted as LUADm, LUAD patients without driver mutation (n = 6), denoted as LUADn, and LUSC patients without driver mutation (n = 16), denoted as LUSCn. Interestingly, LUSCn patients have significant higher ITH_CNA comparing to patients of LUADm, while there was no statistical significance in terms of ITH_GEX (Fig. 2e). This finding also suggested patients with driver mutations may be phenotypically influenced beyond genomic alternations. The comparison between this cohort and a cohort from public data revealed increased ITH_GEX scores of late-stage patients12(Fig. S4C).

Plasticity of lung epithelial cells and their developmental trajectories into malignant tumor cells

All identified alveolar cells express the canonical markers (CLDN18, SFTPA1, SFTPC) of Alveolar Type 2 cells (AT2) without expressing Alveolar Type 1 cell markers (CAV1, AGER). Further clustering analysis unveiled two distinct cluster of AT2 cells, denoted as AT2-1 and AT2-2 (Fig. 3a). AT2-2 resembled a normal AT2 phenotype with common AT2 markers SFTPA and transporter ABCA3 upregulated (Fig. 3b). In contrast, AT2-1 expressed cell proliferation and cell migration related genes, such as CEACAM6, KITLG, and FOXC1, implying a phenotypic change towards malignancy. Epithelial cells could be further separated into ciliated epithelial cells, club cells, and basal cells (Fig. 3c, d).

Previous studies showed that AT2 cells and club cells could both develop into LUAD cells, while basal cells and club cells are potential progenitors of LUSC16,17. Therefore, we organized AT2 cells, club cells and LUAD cancer cells according to their developmental trajectories (Fig. 3e). The inferred pseudotime paths showed AT2 cells and club cells transited into LUAD tumors independently. In contrast, basal cells seemed to act as a transitional state between club cells and LUSC tumor cells (Fig. 3f). Besides such distinct signatures, we found tumor cells of some patients clustered closely at the end of the branches, implying a homogeneous and terminal phenotype, while others have more diverse and heterogeneous profiles spreading along cancer developmental trajectories.

Advanced NSCLC TME revealed a rich program of stromal and immune components

To further identify subgroups of each stromal and immune cell types, we clustered and annotated them individually. We identified five subtypes of endothelial cells (EC) including lymphatic, venous, and arterial endothelial cells (LEC, VECs, and AECs), tip cells, and an EC cluster enriched with interferon induced genes (Fig. S5and Supplementary Data 2). Furthermore, we divided fibroblasts into pericytes and fibroblasts, including six subclusters of fibroblasts (Fig. S6 and Supplementary Data 3). For immune cells, our data revealed two B cell subgroups and seven different plasma cells (Fig. S7). Myeloid cells, especially macrophages, have a broad range of phenotypes and could be divided into 10 different groups (Fig. S8). Dendritic cells (DCs), including plasmacytoid dendritic cells (pDCs), conventional type 1 and 2 DCs (cDC1 and cDC2), and mature DCs were also discovered. Neutrophils have two distinct clusters, expressing potential polymorphonuclear myeloid-derived suppressor cells (PMN-MDSCs) related genes such as LOX-1 (OLR1) to different extents (Fig. S9).

Detailed analysis of T cells uncovered Th17-like cells and their potential interconversion with Tregs

Within tumor-infiltrating T cells, we identified CD4+ naïve T cells, CD4+ Tregs, CD4+ T helper 17-like T cells (Th17-like), CD8+ effector T cells, CD8+ exhausted T cells, and Natural Killer (NK) cells according to expression of their respective markers (Fig. 4a). T cell subtypes were confirmed by supervised cell-type annotation based on previously studied T subtype expression profiles9 (Fig. 4a). To further characterize two NK clusters (CD3D−, KLRD1+, NKG7+), we referred to the CD16+ (FCGR3A) cluster as NK-1 and CD16− cluster as NK-2 (Fig. 4b). NK-1 contains upregulated transcripts encoding fractalkine receptor (CX3CR1) and fibroblast growth factor binding protein 2 (FGFBP2), both involved in lymphocyte cytotoxic functions. NK-2 had higher expression of tissue-resident markers such as CD49a (ITGA1), CD103 (ITGAE), and ZNF683. Co-inhibitory immune checkpoints including CTLA4 and TIGIT were enriched in CD4+ Tregs and CD8+ exhausted T cells (Fig. 4c). However, LAG3 was mainly expressed in CD8+ exhausted T cells, which is consistent with previous findings9.

We then performed trajectory analysis on CD4+ T cells to determine their developmental pathways within TME using both Slingshot18 and Monocle19. Slingshot revealed a transitional relationship between Tregs and Th17-like cells, originated from naïve cells (Fig. 4d). Uncovered by Monocle, naïve cells differentiated into two major branches, Tregs and a proliferating population (Fig. 4e, f). Interestingly, Th17-like cells, confirmed by expression of their master transcription factor RORC, showed a transitional phenotype spreading along the developmental pathway from naïve cells to Tregs (Fig. 4f). The CD4+ Th17-like population marked by high expression of gene KLRB120 is, to our knowledge, the first report of Th17-like cells identified in NSCLC tumor environments by scRNA-seq. As supported by literature21, natural Tregs (nTregs), a subset of Tregs, are believed to interconvert with Th17-like cells. This result revealed a complex and delicate interplay between Tregs and Th17-like cells and highlighted the importance of their balance in adaptive immune responses to tumor antigens22.

Both NSCLC subtypes and ITH shaped the immune landscape in TME

To investigate if tumor subtypes and their ITH levels affect their microenvironment, we compared the cell-type composition of NSCLC by their histology and their driver mutation status. We found that neutrophils were significantly depleted in all LUAD patients (Fig. 5a). While comparing LUAD patients with and without oncogenic driver mutations, a macrophage cluster with highly expressed CCL13 was enriched in the group of mutated tumors (Fig. 5b). The proportions of the tissue-resident macrophage cluster expressing the scavenger receptor MARCO and CXCL5, and cDCs also exhibited significant differences among the three groups (Fig. 5b). Interestingly, cDC2 displayed Langerin (CD207) expression, which was inferred to be dictated by the environment23. TCGA survival analysis revealed that CD207 is a prognostic marker for LUAD, but not for LUSC (Fig. 5c). However, MARCO is not associated with clinical outcomes of both LUSC and LUAD. Since we identified two subtypes of MARCO+ alveolar macrophages, these results combined implied the multifunctional roles of tissue-resident myeloid cells. We next investigated the correlation between ITH scores and the immune cell compositions. We found neutrophils and two subtypes of macrophages were positively correlated with ITH_GEX, while plasma cells were negatively correlated with ITH_GEX (Fig. 5d). This finding potentially suggested high immunosuppressive environment and low cancer killing ability for patients with high ITH_GEX. Overall, we found the myeloid compartment is the mostly affected by tumor subtypes and ITH levels instead of tumor-infiltrating lymphocytes.

Divergent intercellular networks observed among LUADn, LUADm, and LUSC

In order to explore the interplay among cell types within the tumor microenvironment, we performed a cell–cell interaction analysis and showed a prominent interaction between cancer cells and endothelial cells, fibroblasts and macrophages (Fig. 6a). Analysis of the interacting molecules across cells showed a complex network with the interplay of oncogenic pathways as EGFR, NOTCH, WNT, with PDGF and inflammatory signaling pathways, in particular affecting TNF-a and chemokine responsive pathways (Fig. 6b). Notably, VEGFA-mediated protein–protein interactions and two analogous immune checkpoint pathways CD226-TIGIT-CD96 and CD274 (PD-L1)-CTLA4-CD28 were also identified within the interaction network.

Cancer cells expressed high levels of ligands CXCL1, CXCL2, CXCL3, and CXCL8, signaling to the receptors CXCR1 and CXCR2 expressed by neutrophils (Fig. 6c). Some of the LUSCn and LUADn patients showed increased interactions between DCs and T cells, including CXCR3 and its partners, suggesting strong effector T cell trafficking and recruiting24. We also confirmed activation of the CXCL12-CXCR4 pathway between tumor and sprouting endothelial cells (AECs and tip cells) described in Fig. S5. Regarding growth factors, the majority of tumors, regardless of subtypes, had very strong signals of VEGF interactions between tumor and various types of endothelial cells (Fig. 6d). PDGF signaling, on the other hand, was activated between tumor and cancer-associated fibroblast (CAF) cells. A distinct pattern for LUSC patients is the activation of FGF pathways among stromal cells and tumor cells, also supported by previous studies25. Since only a portion of LUSC patients have FGF pathways activated, patient stratification may be important for the usage of drugs targeting FGF pathways.

For patients’ immune environment, we found macrophages, instead of cancer cells, played a major role in inhibiting T cell functions through checkpoint pathways (Fig. 6e). Different subgroups showed different dominant pathways. For example, LUADm have high levels of activation of the TIGIT pathway but low levels of activation of TIM3 (HAVCR2) pathway. We did not detect any significant activation of the PD1/PD-L1 axis except for a few LUSC patients, potentially due to the low expression of PD1/PD-L1 on the transcriptomic level. Interaction analysis performed on a public dataset confirmed the similar activation state of checkpoint pathways for late-stage LUAD (Fig. S10A). Interestingly, an early-stage LUAD patient in the same dataset showed opposite activation status of TIGIT and TIM3 with respect to the late-stage patient (Fig. S10B). Nevertheless, by cellular network analysis of the scRNA-seq data, we generated a comprehensive view of patients’ TME including angiogenesis, CAF activation, recruitment of immunosuppressive cells, T cell activation, and detailed activation profiles of checkpoint pathways. Therapeutics related interactions were heterogeneous even within the same subtype of lung cancer, highlighting the needs for more precise biomarkers to increase the drug efficacies.

In this study, we present the valuable comprehensive landscape of cancer cells, immune cells, and stromal cells in advanced NSCLC by scRNA-seq analysis. We were able to identify 11 major cell types from advanced NSCLC, including 48 subtypes besides cancer cells, the majority of which are consistent with previous studies. We focused on the cancer cells, which were not studied extensively at single-cell level in the previous literature. The shared arm-level CNAs were consistent with the observations from previous genomic sequencing data26, indicating a representative cohort of advanced NSCLC tumor types. Based on a quantitative approach to define inter- and intratumor heterogeneity, we unmasked a broad range of clonality, homogeneity, and the complexity beyond current classification systems of advanced NSCLC. In general, LUSC has higher ITH than LUAD. However, our data call for a more precise profiling of individual patients on the cellular levels beyond the traditional pathological definitions. For example, specimen P7 is a LUSC tumor with strong TP63/CK expressions and weak NAPSA expression (Fig. S2A). Interestingly, the major clone of this patient only represents less than 75% of its cancer cells (Fig. S3C). Further investigation showed one of its minor clones clustered together with many LUAD patients.

We also identified rare cell types such as FDC and Th17-like lymphocytes. The existence of FDC indicated the formation of lymphoid follicles, which usually correlates with favorable clinical outcomes27. In tumor-infiltrating CD8+ T cells, there are more exhausted T cells than cytotoxic T cells, which is opposite to what is observed in early stage, resectable NSCLC patients9. Notably, the myofibroblast to fibroblast ratio in our study was remarkably high compared to healthy or asthma lungs28. Thus, CAFs with myofibroblast characteristics may act as an important malignant signature for advanced stage lung cancer. Certain cell subtypes identified in this study were previously determined to be associated with drug responses. For example, CXCL9+ Mac was enriched in patients responding to immunotherapy29.

From the cellular composition analysis, we showed neutrophils to be enriched in LUSC. This phenomenon has been demonstrated by previous studies in NSCLC that neutrophils are more abundant in human LUSC compared to LUAD due to differences in TME30. In a subsequent study, the master transcription factor SOX2, a lineage-specific oncogene for squamous cell carcinomas, was found to be overexpressed and to promote tumor associated neutrophil (TANs)-accumulation by upregulating CXCL5(the mouse homolog of human CXCL6) expression31. Therefore, different neutrophil infiltration features were proposed to be regulated by tumor-intrinsic driver mechanisms. On the other hand, cancer and neutrophils have stronger interactions in LUAD patients. The combined observations suggested complex and diverse functions of neutrophils within TME. Our study also revealed a correlation of tumor heterogeneity and neutrophil contents. Given that separate studies showed high neutrophil contents and high tumor heterogeneity related to immunotherapy failure respectively32,33,34, our data further bridged the gap between high neutrophil contents and high tumor heterogeneity. The mechanism of interplay between tumor heterogeneity and tumor-infiltrating neutrophils might be a key element to explain the difference in immunotherapy efficacy.

By mapping cells in TME and their possible functions, by identifying more cell types and their marker genes, and by highlighting intratumor cell–cell interactions, we presented here a comprehensive collection of datasets, which provide deep insights on advanced NSCLC. Some of our findings are unreported and will need further functional validation. Despite this limitation, it can serve as valuable resources and a proof-of-concept study for future research to identify biomarkers and targets for treatment and enable personally tailored therapeutic decisions for patients with advanced NSCLC.