High-resolution single-cell atlas reveals diversity and plasticity of tissue-resident neutrophils in non-small cell lung cancer

Abstract

Non-small cell lung cancer (NSCLC) is characterized by molecular heterogeneity with diverse immune cell infiltration patterns, which has been linked to therapy sensitivity and resistance. However, full understanding of how immune cell phenotypes vary across different patient subgroups is lacking. Here, we dissect the NSCLC tumor microenvironment at high resolution by integrating 1,283,972 single cells from 556 samples and 318 patients across 29 datasets, including our dataset capturing cells with low mRNA content. We stratify patients into immune-deserted, B cell, T cell, and myeloid cell subtypes. Using bulk samples with genomic and clinical information, we identify cellular components associated with tumor histology and genotypes. We then focus on the analysis of tissue-resident neutrophils (TRNs) and uncover distinct subpopulations that acquire new functional properties in the tissue microenvironment, providing evidence for the plasticity of TRNs. Finally, we show that a TRN-derived gene signature is associated with anti-programmed cell death ligand 1 (PD-L1) treatment failure.

Keywords: cell-cell communication; patient stratification; single-cell sequencing; therapy response; tissue-resident neutrophils.

Graphical abstract

Figure 1

Schematic outline of the overall concept used in this study (A) Summary of the data integration and analysis workflow. (B) Overview of the core NSCLC atlas and the epithelial, immune, and stromal/endothelial components depicted as uniform manifold approximation and projection (UMAP) plots. (C) Fractions of depicted cell types per scRNA-seq platform. (D) UMAP of the UKIM-V dataset (n = 17) colored by cell type. (E) Core atlas extended by Leader and UKIM-V-2 datasets. (F) Cell-type composition by histopathological tumor type (LUAD, LUSC). FDR = 0.1. (G) Immunohistochemistry staining of neutrophils (ASD⁺ cells), macrophages (CD68⁺ cells), and T cells CD4 (CD4⁺ cells) per high-power field (HPF) in LUAD (n = 55) versus LUSC (n = 55). Evaluation was performed by two separate expert lung pathologists (C.K. and S.P.). The horizontal line represents the median, and whiskers extend to the interquartile range; Wilcoxon test, ^∗p < 0.05, ^∗∗∗∗p < 0.0001. See also Figures S1 and Tables S1, S2, and S3.

Figure 2

Tumor immune phenotypes in NSCLC (A) Patient characteristics and stratification of the tumor immune phenotypes. Tumor type (histopathological) refers to the histological subtypes as provided by the original datasets based on pathological assessment; tumor type (transcriptomic) is based on the most abundant transcriptomically annotated cancer-cell subtype in the scRNA-seq atlas. (B and C) Differential activation of (B) PROGENy cancer pathways and (C) CytoSig cytokine signaling signatures in cancer cells between the four tumor immune phenotypes. Heatmap colors indicate the deviation from the overall mean, independent of tumor histology and stage. White dots indicate significant interactions at different FDR thresholds. Only cytokine signatures with an FDR <0.1 in at least one patient group are shown. See also Figure S2.

Figure 3

Cellular crosstalk analysis (A) Circos plot of the cellular crosstalk of cancer cells toward the major immune cells in LUAD versus LUSC. Shown are the top 10 differentially expressed cancer cell ligands. Red interactions are upregulated in LUAD, and blue interactions are upregulated in LUSC. (B) Cancer-immune cell crosstalk in each patient subtype. Top panel: differentially expressed ligands of cancer cells in each subtype (B, M, T, ID) (DESeq2 on pseudo-bulk, FDR < 0.1). Bottom panel: respective receptors and the expression by cell type. Dot sizes and colors refer to the fraction of cells expressing the receptor and gene expression, respectively, averaged over all patients. Dots are only shown for receptors that are expressed in at least 10% of the cells of the respective cell types. See also Figure S3.

Figure 4

Association of cellular composition and distinct genotypes and survival in the TCGA data (A–E and G) SCISSOR analysis relating phenotypic information from bulk RNA-seq data from TCGA with single cells. UMAP plots indicate the position of cells positively (blue) or negatively (red) associated with mutation or better survival. A log2 ratio >0 indicates a positive association with mutation or better survival, respectively. Shown are cell types with a log2 ratio significantly different from 0 at an FDR <0.01 (paired Wilcoxon test). (A) Association of cellular composition with KRAS mutation in patients with LUAD (n = 156). (B) Association of cellular composition with EGFR mutation in patients with LUAD (n = 98). (C) Association of cellular composition with STK11 mutation in patients with LUAD (n = 141). (D) Association of cellular composition with STK11 mutation in patients with LUSC (n = 83). (E) Association of cellular composition with overall survival. (F) Kaplan-Meyer plot of patients with high (top 25%) and low (bottom 25%) B cell fractions of TCGA patients with lung cancer as determined by deconvolution with EPIC. p value has been determined using CoxPH regression using tumor stage and age as covariates. (G) Association of cellular composition with overall focusing on CD8⁺ T cell subclusters. CD8⁺ T cell subclusters were annotated based on gene sets from Oliveira et al. See also Figure S4.

Figure 5

Characterization of tissue-resident neutrophils using scRNA-seq (A and B) UMAP of tissue-resident neutrophils (TRNs) from the extended atlas, (A) classified into tumor-associated neutrophils (TANs) and normal-adjacent associated neutrophils (NANs) and (B) colored by histology (as defined by histopathological assessment). (C) Neutrophil fractions (as percentage of leucocytes) by flow cytometry of LUAD and LUSC tumor tissue (LUAD n = 47, LUSC n = 16; Wilcoxon test, ^∗∗p < 0.01). The horizontal line represents the median, and whiskers extend to the interquartile range. (D) Candidate TAN genes. Each dot refers to a patient with at least 10 neutrophils in both NAN and TAN groups. Lines indicate the mean per study. p values are derived from a paired t test and adjusted for FDR. (E) Expression levels of VEGFA in various cell types in primary tumor samples. Each dot represents a patient with at least 10 cells (median values, boxes represent the interquartile range [IQR], whisker data points within 1.5 times the IQR). (F) Transcription factor analysis of TAN versus NAN using DoRothEA. Each dot represents a single patient, and bars are the mean of all patients. p values are derived using a paired t test and are FDR adjusted. Shown are transcription factors with a mean score difference >0.2 and an FDR <0.1. (G) Comparison between tumor and normal-adjacent samples for selected candidate genes using flow cytometry. Each dot represents a patient that was not part of the scRNA-seq dataset. Paired Wilcoxon test, ^∗p < 0.05, ^∗∗∗∗p < 0.0001. CD16: the horizontal line represents the median, and whiskers extend to the IQR. (H) Selected multiplex immunofluorescence (M-IF) staining of LOX-1 (red) and pancytokeratin (green) in tumor tissue and matched normal-adjacent lung tissue of a patient with LUSC. Scale bar: 100 μm. See also Figure S5.

Figure 6

Tissue-resident neutrophil subtypes in NSCLC (A) UMAP of all TRNs colored by TAN and NAN subclusters. The neutrophil clusters derive from 85 patients, 42 of whom have >10 neutrophils. (B) Top 5 markers for each TAN and NAN cluster. The marker gene quality is reflected by the area under the receiver operator characteristics curve (AUROC; 1 = marker gene perfectly distinguished the respective cluster from other clusters in all patients; AUROC 0.5 = no better than random). (C) Quantification of HLA-DR expression by flow cytometry of tumor and normal-adjacent tissue. Each dot represents the mean of each patient. Paired Wilcoxon test, ^∗∗p < 0.01. (D) UMAP of TRNs from the UKIM-V dataset with RNA velocity vectors projected on top. (E) Partition-based graph abstraction (PAGA) based on RNA velocities, projected on the UMAP plot. Edges represent cell-type transitions called by PAGA. (F) Outgoing interactions of TRN subclusters with cancer cells and CD8⁺ T cells. Top panel: differentially expressed ligands in each subcluster (FDR <0.01, abs. log2 fold change >1). Heatmap colors clipped at ±3. Bottom panel: respective receptors and the expression by cell type. Dots are only shown for receptors that are expressed in at least 10% of the respective cell types. (G) UMAP of the extended atlas colored by the score of the TRN gene signature (38 genes with high specificity for TRNs) (H) Heatmap of the TAN and NAN gene signatures across the TRN subclusters. Colors indicate the mean gene expression across patients in the respective clusters. (I–K) Predictive value of the TRN signature in bulk RNA-seq data from the OAK and POPLAR cohorts of patients with NSCLC treated with atezolizumab (anti-PD-L1) or docetaxel (chemotherapy). (I) Comparison of non-responders (progressive disease) with responders (complete response, partial response) treated with atezolizumab, shown for each histotype. (J) Kaplan-Meyer plot comparing patients treated with atezolizumab with high (top 25%) and low (bottom 25%) TRN signature scores. p value has been determined using CoxPH regression including cohort and histology as covariates. (K) Kaplan-Meyer plot comparing patients treated with docetaxel with high (top 25%) and low (bottom 25%) TRN signature scores. See also Figure S6 and Tables S4 and S5.