Multifaceted roles of tumor-associated neutrophils in NSCLC

Introduction

Neutrophils are the most abundant population of immune cells in human peripheral circulation, and the defense of the innate immune system at the frontline against invading pathogens.^1–4 Notwithstanding, the roles played by neutrophils in regulating tumor immunity by interacting with other immune cells in the tumor microenvironment (TME) had been ignored for a long time. Recently, neutrophil-based indexes [eg, neutrophil-to-lymphocyte ratio (NLR) and derived NLR (dNLR)] are receiving increasing attention for their potential predictive and prognostic value in anti-cancer treatment. For example, a retrospective cohort study performed by Valero et al⁵ illustrated that an elevated pretreatment blood NLR was associated with inferior progression-free survival, overall survival (OS), response rate, and clinical benefit across multiple cancer types administered with immune checkpoint inhibitors including non-small cell lung cancer (NSCLC), implying its potential predictive and prognostic value in this scenario. In addition, the function of neutrophils in tumor immunity is gradually being disclosed, albeit to some extent contradictory. A number of studies have demonstrated that tumor-associated neutrophils (TANs) perform a completely different role in different periods of tumorigenesis and development. TANs are anti-tumorigenic in the early stages of malignant tumors, but pro-tumorigenic in the late stages. Eruslanov et al⁶ found that TANs exhibited CD62L^loCD54^hi activated phenotype at early stages of lung cancer with a broad variety of chemokine receptors (eg, CCR5, CCR7, CXCR3, and CXCR4) and produced large amounts of proinflammatory factors (eg, MCP-1, IL-8, and MIP-1α). Consequently, TANs stimulated T cell responses including T cell proliferation and IFN-γ release to inhibit tumor growth.⁶ Albanesi et al⁷ reported that eliminating TANs from TME at early stages of tumors could promote tumor growth, which also confirmed that TANs inhibited the growth of tumor in early periods of tumor initiation and progression. Interestingly, Mishalian et al⁸ observed in LLC (lung cancer) and AB12 (mesothelioma) tumor-bearing murine models that TANs acquire pro-tumorigenic properties during tumor progression in a time-dependent manner, denoting the functional phenotype of TANs might be educated by the ongoing evolvement of TME.

Nowadays, we describe this phenomenon as “the plasticity of neutrophils” that neutrophils being divided into two functional phenotypes, an anti-tumorigenic one (ie, N1-TAN) or a pro-tumorigenic one (ie, N2-TAN), which could be educated and modulated by the constantly changing TME.^9–15 For example, Fridlender et al¹⁰ proved that in tumor-bearing murine models, tumor-infiltrating TAN could polarize to N2-TAN under the induction of transforming growth factor (TGF)-β signaling, and once TGF-β inhibitors were administered, TAN could polarize to N1-TAN conversely. Similarly, Andzinski et al⁹ revealed that type I interferon (IFN) could induce the transformation from N2 to N1 in both mouse and human TANs. Meanwhile, peripheral circulating neutrophils possess identical plasticity to tumor-infiltrating neutrophils. Sagiv et al¹⁵ identified heterogeneous circulating neutrophils in mouse and human based on density gradient [ie, low/high-density neutrophils (LDN/HDN)], and they observed that the immunosuppressive LDN consisted of both mature and immature neutrophils. Subsequently, a BrdU pulse-chase experiment was conducted manifesting that at least part of the LDNs were transitioned from HDNs, and ultimately they further verified the plasticity of peripheral circulating neutrophils with tumor progression by a series of in vivo experiments.¹⁵

In general, the two-faced roles of TANs in tumor development and progression should be attributed to their plasticity and heterogeneity.¹³ Nonetheless, neutrophil subtypes might be far more intricate than dichotomous functional phenotypes (eg, N1/N2 TAN or LDN/HDN) due to the lack of in-depth exploration. For example, Xie et al¹⁶ clustered eight distinct neutrophil subtypes in infected mice via single-cell RNA-sequencing (scRNA-seq), and distinct subtypes presented different transcriptional profiles and functions, respectively. Given this, we presume that tumor-infiltrating TANs might be likewise complex and comprise multiple distinct subtypes similarly. However, there have been insufficient studies on the transcriptional landscape of TANs that have been reported so far. With the growing interest in the functions of TANs in regulating tumor immunity, it is imperative to comprehensively profile the transcriptome signature and discover diverse TAN subsets at single-cell resolution.

Our previously published paper¹⁷ mapped the cell type-specific transcriptome landscape of cancer cells and their TME in advanced NSCLC, and the correlation of tumor heterogeneity with TANs was revealed, which we believe is worth further investigation. However, as far as we know, neutrophils are fragile, especially sensitive to handling procedures, and additionally express low amount of mRNA molecules, which collectively made it difficult to characterize TANs at single-cell resolution. Accordingly, not all samples possess enough TANs for further scRNA-seq analysis. Ultimately, nine samples with a sufficient number of TANs were enrolled in this study, and we identified subgroups of TANs, clustering and annotating them individually to further clarify the heterogeneous transcriptional landscape of TANs. In this manuscript, we would like to report our findings and hope to provide a theoretical basis for subsequent functional research in the future.

Materials and Methods

Patients

All samples were obtained from patients with pathologically confirmed advanced and unresectable NSCLC from November 2018 through August 2019. As this is an extension of our previous research,¹⁷ nine samples (P19, P22, P26, P27, P34, P36, P37, P40, and P42 were re-named as P1-9 here, respectively) with sufficient TANs were selected for inclusion in this study. Specimens from primary lung tumors were collected using diagnostic procedures including CT-guided transcutaneous needle biopsy or transbronchial lung biopsy. Actionable gene alterations (eg, EGFR, KRAS, ALK, ROS1, RET, HER2, BRAF, and MET exon 14 skipping) were detected. Clinicopathological characteristics including sex, age at sampling, smoking history, Eastern Cooperative Oncology Group performance status (ECOG PS), TNM stage, and comorbidities were extracted from electronic medical records. Smoking history was categorized into smokers (including both current and former smokers) and non-smokers (individuals who had smoked <100 cigarettes in their lifetime). Patients with concomitant autoimmune diseases, interstitial lung disease, chronic obstructive pulmonary disease, or active infectious diseases (such as, tuberculosis, AIDS, hepatitis B, etc.) were excluded from this study. All patients enrolled provided written informed consent. The Ethical Committee of Shanghai Pulmonary Hospital Affiliated to Tongji University approved this study (No. K18-089-1).

Tissue Dissociation and the Preparation of Single-Cell Suspensions

The fresh tissues collected were immediately stored in the GEXSCOPE^® Tissue Preservation Solution (Singleron Biotechnologies, Nanjing, China) at 2°C to 8°C. Hanks Balanced Salt Solution-washed tissue samples were then cut into smaller pieces and digested in 2 mL GEXSCOPE^® Tissue Dissociation Solution (Singleron Biotechnologies, Nanjing, China) per the manufacturer’s instructions. More detailed procedures were provided in our previous paper.¹⁷ TC20™ Automated Cell Counter (Bio-Rad Laboratories, Inc.) was utilized to measure the cell concentration and viability of prepared single-cell suspensions.

Steps of Single-Cell RNA Sequencing and Bioinformatics Analysis

In our study, scRNA-seq was implemented via the GEXSCOPE^® platform (Singleron Biotechnologies, Nanjing, China) as recorded in our previous research.¹⁷ Briefly, single-cell suspensions were added onto a microfluidic chip (SCOPE-chip™, Singleron Biotechnologies) and scRNA-seq libraries were set up referring to the manufacturer’s instructions (GEXSCOPE^® Single-Cell RNAseq Library kits, Singleron Biotechnologies). Subsequently, the sequencing of resulting libraries was performed on Illumina HiSeq ×10 instrument with 150 bp paired end reads. Using scopetools, gene expression matrices were generated from raw reads.

Following strict quality control, putative cell doublets, stressed cells, and dead cells were filtered out, and neutrophils were identified and clustered by their canonical marker genes (CSF3R,^18,19 S100A8/9,^20–22 FCGR3B^18,23). Uniform Manifold Approximation and Projection (UMAP) was applied for two-dimensional visualization. Eventually, a total of seven neutrophil subtypes were defined by the top differentially expressed genes (DEGs) detected in each cluster.

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) Enrichment Analysis

Based on the DEGs mentioned above, GO and KEGG enrichment analysis was conducted. The GO enrichment analysis includes biological process, cellular component, and molecular function. KEGG analysis can profile pathways related to gene functions. R package clusterProfiler²⁴ was adopted for enrichment analysis, and R package ggplot2 was used for visualization.

Developmental Trajectory Analysis

The pseudo-temporal trajectory of neutrophil differentiation and maturation was reconstructed by Monocle2.²⁵ Marker gene expression along the developmental trajectory of neutrophil subtypes was depicted by a heatmap.

Composition Analysis of Neutrophil Subtypes Between Patient Groups

R package ggpubr was applied for statistical determination and visualization to evaluate whether the composition of neutrophil subtypes was significantly different between diverse patient groups. The t-test was adopted to test the statistical significance of the comparison between two groups. P < 0.05 (two-sided) was considered statistically significant.

Intercellular Interaction Analysis

Detailed methodology was reported previously.¹⁷ Briefly, CellphoneDB²⁶ was used to reveal the interaction between cell types. Permutation number for calculating the null distribution of average ligand-receptor pair expression in randomized cell identities was set to 1000. Individual ligand or receptor expression was thresholded by a cutoff based on the average log gene expression distribution for all genes across each cell type. Predicted interaction pairs with P 0.1 were considered to be significant and visualized by heatmap plot and dot_plot in CellphoneDB. In the diagram of cell interaction network, nodes represented cell types, and edge weights were calculated based on the number of interactions between two cell types. Visualization of the network was done with Cytoscape.²⁷

Overall Survival Analysis

The clinical profile and RNA sequencing profile of lung cancer, including lung adenocarcinoma (LUAD) and lung squamous cell carcinoma (LUSC), were downloaded from The Cancer Genome Atlas data portal (https://tcga-data.nci.nih.gov/tcga/). After matching clinical information and RNA expression, 1026 patients were included for further survival analysis. Cox regression was applied to discover genes with statistically significant difference in survival analysis. The target genes were illustrated by Kaplan–Meier survival plots and P.value was displayed. All these analyses were performed with R (4.0.2) and RStudio (2022.07.2).

Results

The Clinicopathological Characteristics of Enrolled Patients

A diagram delineating the study procedure is shown in Figure 1A. In total, nine patients were included for this study, and most were male (n = 7, 78%), LUSC (n = 7, 78%), stage IV (n = 5, 56%), wild type (n = 7, 78%), and systemic anti-cancer treatment naïve (n = 8, 89%). The proportion of smokers and non-smokers was roughly the same. The detailed clinicopathological characteristics of enrolled patients are provided in Figure 1B and Supplementary Table 1. A total of 8659 cells from 9 patients were annotated as 10 major cell types, including cancer cells, epithelial cells other than cancer cells, immune cells (mononuclear phagocytes, T cells, B cells, neutrophils, mast cells, and follicular dendritic cells), and stromal cells (fibroblasts and endothelial cells) (Figure 1C). Compared with immune and stromal cells, cancer cells presented higher heterogeneity and sample-specific transcriptomic signatures (Figure 1D). Moreover, the proportions of distinct cell types varied greatly among samples (Figure 1E).

Figure 1 Study design and patients enrolled. (A) A diagram delineating the study workflow for single-cell profiling of tumor-associated neutrophils in advanced non-small cell lung cancer. (B) The clinicopathological characteristics of enrolled patients (n = 9), including sex, smoking history, TNM stage, histology, and actionable gene alteration status. (C) UMAP plot of 8659 cells from 9 patients, colored by 10 major cell types. (D) UMAP plot of all cells, colored by patients. (E) Major cell-type composition of each patient. Abbreviations: LUAD, lung adenocarcinoma; LUSC, lung squamous cell carcinoma; UMAP, Uniform Manifold Approximation and Projection; EGFR, epidermal growth factor receptor; MP, mononuclear phagocyte; fDC, follicular dendritic cells.

Significant Heterogeneity of TANs in Human NSCLC

Eventually, 1820 TANs in total were analyzed. Neutrophils were clustered according to their canonical cell markers, CSF3R,^18,19 S100A8/9,^20–22 and FCGR3B^18,23 as previously published. Next, data were subclustered and visualized using UMAP plots. As shown in Figure 2A, seven distinct neutrophil subtypes were identified, denoted as N_1-7. Interestingly, we observed that neutrophils formed a discontinuum of states with N₇ being separate from N_1-6, and all seven subtypes were represented by merely three of nine patients in different proportions (Figures 2B and C, Supplementary Table 2). In addition, the proportions of seven subtypes were dramatically disparate among the nine patients (Figure 2C), which might be due to different histologies, various disease progressions, and the heterogeneity of specimens. The overwhelming majority (87%) of N₇ were detected in P7, a stage IV LUSC patient harboring L858R point mutation (Figures 2B and C, Supplementary Table 2).

Figure 2 Transcriptional landscape of TANs in advanced NSCLC via scRNA-seq analysis. (A) UMAP plots of 1820 neutrophils analyzed from 9 biopsy specimens, clustered into 7 distinct neutrophil subtypes by different colors. (B) UMAP plots of all TANs, colored according to specimen sources. (C) Various composition of distinct neutrophil subtypes in different patients. (D) A heatmap exhibiting top DEGs of each neutrophil cluster. Abbreviations: TANs, tumor-associated neutrophils; NSCLC, non-small cell lung cancer; scRNA-seq, single-cell RNA-sequencing; UMAP, Uniform Manifold Approximation and Projection; DEGs, differentially expressed genes.

To clarify the possible roles played within TME by different neutrophil subtypes, we inspected the DEGs of each neutrophil cluster. A heatmap exhibiting top DEGs of different neutrophil clusters is presented in Figure 2D. N₁ expressed masses of canonical neutrophil marker genes (S100A8 and S100A9), suggesting that N₁ was most likely to be a cluster of mature neutrophils. N₂ seemed to be involved in lipid metabolism (PLIN2 and LRPAP1) and could also function to promote angiogenesis and endothelial cell growth (VEGFA). A series of cytokines-related genes (CCL3, CCL4, CCL20, CCL3L1, CCL4L2, CXCL8, and CXCL2) were expressed by N₃, indicating that N₃ was able to regulate the status of TANs and further promote tumor growth. Considering that they both expressed genes encoding the Heat Shock Protein (HSP) family members (HSP90AB1, HSP90AA1, HSPA1A, and HSPH1) and certain proteases (CTSB, MMP12), we believed that N₄ and N₆ were the most closely related clusters. Additionally, N₅ might act as a regulator of innate immune responses to viral infections, since multiple genes encoding IFN-inducible antiviral proteins were expressed (ISG15, IFIT3, IFIT2, IFI6, and RSAD2). Last but not least, N₇ formed a cluster apart from the continuum of the other TAN states. Indeed, N₇ cells represent a small portion of TANs but show a unique transcriptional signature, expressing genes related to multiple biological processes such as DNA repair, RNA synthesis, and cell growth. A detailed list with associated functions of top DEGs detected is summarized in Figure 3A. GO and KEGG enrichment analysis further verified the possible functions of DEGs and related pathways (Figures 3B–H).

Figure 3 Top DEGs examined among seven distinct TAN subtypes with their possible functions, and GO and KEGG enrichment analysis of these DEGs for signaling pathways prediction. (A) A list of top DEGs detected among diverse TAN subtypes and their associated functions. GO and KEGG enrichment analysis of top DEGs from (B) N1, (C) N2, (D) N3, (E) N4, (F) N5, (G) N6, and (H) N7 cluster to predict possible signaling pathways, respectively. The vertical axis represents the possible signaling pathway, the bubble size indicates the number of DEGs in the pathway, and the bubble color corresponds to different p.adjust values. Abbreviations: DEGs, differentially expressed genes; TAN, tumor-associated neutrophil; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; IFN, interferon; DNA, deoxyribonucleic acid; RNA, ribonucleic acid.

TAN Subtypes Distribution Between LUAD and LUSC

As we previously published, neutrophils were significantly depleted in patients with LUAD regardless of driver gene alteration status.¹⁷ Accordingly, we aim to elaborate on the neutrophil-subtype composition in different NSCLC histologies. We found that the number of TANs, especially the proportion of N₃, was significantly higher in LUSC than that in LUAD. Conversely, the proportion of N₅ in LUSC was significantly lower than that in LUAD. No statistically significant differences were detected in the distribution between NSCLC histologies regarding other neutrophil clusters, possibly due to the small sample size in the LUAD group (Supplementary Figure 1).

Developmental Pseudo-Time Analysis of Marker Gene Expression Among Diverse TAN Subtypes

Arranging single cells in a pseudo-temporal way (pseudo-time analysis) based on our scRNA-seq data was conducted with Monocle2,²⁵ in order to reconstruct the developmental trajectory of neutrophil differentiation in NSCLC. Overall, a total of three branches were observed in the pseudo-time path (Figures 4A and B), and cells from various clusters were widely distributed along the pseudo-time path (Figure 4B). We believed that the 1-end was the start of neutrophil differentiation, and the 3-end was the end of neutrophil differentiation in terms of the expression levels of mature neutrophil marker genes (eg, CXCR2, S100A8, DUSP1, and FPR1)¹⁶ (Figure 4A). Thus, N₁ and N₅ might be neutrophil subtypes in a relatively mature status (Figure 4B).

Figure 4 Pseudo-time analysis and the developmental trajectories of TANs in human advanced NSCLC. (A) The developmental trajectories of the pooled TANs, and the inferred direction to differentiation and maturation was from the 1-end to the 3-end. (B) The developmental trajectories of TANs colored by cluster identities. (C) Distribution of each patient on the developmental trajectories, with one color indicating a TAN subtype. (D) Distribution of different TAN subtypes on the developmental trajectories. (E) The dynamics of HVG expression along the pseudo-time with each dot representing the average expression of the specified gene at the specified pseudo-time. (F) The dynamic expressions of genes related to neutrophil differentiation and maturation along pseudo-time axis. Abbreviations: TANs, tumor-associated neutrophils; NSCLC, non-small cell lung cancer; HVG, highly variable genes.

Next, we delineate the developmental trajectories of neutrophil differentiation in each patient (Figure 4C), and the distribution of each neutrophil cluster in the pseudo-time path (Figure 4D), illustrating the temporal and spatial heterogeneity of TANs in NSCLC intuitively. Interestingly, N₇ was mainly concentrated close to the start of differentiation, and almost none was distributed at the end of the pseudo-time path (Figure 4D), which might explain the great independence of N₇ compared with N_1-6 clusters. Moreover, N₄ and N₆ share similar developmental trajectories (Figure 4D), which was consistent with the fact that they share similar DEGs and possible functions as regulators to HSPs in TME.

The top marker genes in each TAN subtype were selected to determine their pseudo-time patterns. As shown in Figure 4E, the pseudo-time patterns of selected marker genes could be classified into two styles. The first style was that the expression of all marker genes increased over pseudo-time, and conversely, the other was that the expression of all marker genes decreased along the pseudo-time axis. The marker genes from predominantly N₁, N₃, and N₅ clusters conformed to the pseudo-time pattern of the first style, suggesting that they were at a later stage of neutrophil differentiation. The second style could be subdivided into three sub-styles: first, the expression of marker genes from the N₆ cluster elevated slightly at the beginning and then slightly declined; second, the expression of marker genes mainly from the N₄ cluster was lower during the entire differentiation process and further decreased to a minimum level at the last phase; third, the expression of marker genes from N₂ and N₇ clusters declined sharply from a maximum level along the pseudo-time axis.

Furthermore, pseudo-time analyses of marker genes closely related to the process of neutrophil differentiation in NSCLC were performed (Figure 4F). Specifically, genes such as DUSP1, FOS, MNDA, RGS2, and S100A8/9 were upregulated as neutrophils matured, which resembled the pseudo-time pattern of a widely recognized mature neutrophil marker gene, CXCR2.^28,29 However, MMP12 and PLN2 show a tendency to be downregulated along the pseudo-time axis.

Cell–Cell Interaction Analysis Within TME

For exploring the interplay between neutrophils and other cells within TME, we performed cell–cell interaction analysis via CellphoneDB, an online platform, which is based on ligand-receptor signaling database.²⁶ We display prominent interactions between TAN subtypes and macrophages, monocytes, tumor cells and T cells (Figure 5A, and Supplementary Figure 2). Notably, TAN subtypes interacted most with multiple macrophage subtypes, because macrophages expressed significantly larger number of receptors corresponding to ligands from neutrophils. The results indicated that all TAN subtypes, except N₂, could be recruited by macrophages, such as CCL13_Macrophage, CCL3_Macrophage, CXCL9_Macrophage, MARCO_CXCL5_Macrophage, and MMP12_Macrophage (Figure 5B), through CXCL8-CXCR2, CXCL3-CXCR2, CXCL2-CXCR2, CXCL1-CXCR2, CXCL8-CXCR1, CXCL1-CXCR1, and CCL3-CCR1 axes (Figure 5C).

Figure 5 Cell and gene interaction networks. (A) Cell interaction networks among all cell types. Nodes represent cell types, and edge weights were calculated based on the number of ligand-receptor pairs between two cell types. The colors of edges correspond to ligand cells. (B) A heatmap displaying the pair numbers of interactions among the predicted cell types. (C) A bubble plot showing selected chemokine-related interactions among the critical immune cell subtypes. Dot size represents values of –log10 (P), and dot color corresponds to log2 mean (Molecule 1, Molecule 2).

N₃ Subtype Was Considered as Predictors for Inferior OS in Advanced NSCLC

As mentioned above, N₃ subtype was determined to express a series of cytokines-related genes, and thus we reasoned that N₃ subtype was likely related to dismal prognosis in advanced NSCLC. Due to the limited sample size, we sought to verify our hypothesis via publicly available datasets as mentioned above. In total, 1026 patients were included for further survival analysis. Eight genes, ABHD5, ARHGDIA, CCL20, IL1A, IRAK2, PDLIM7, PTTG1IP, and SERPINB9, were N₃-specific DEGs compared with other TAN subtypes per our scRNA-seq data. Thereafter, we analyzed the correlation between expression level of those genes and OS among this online dataset cohort. As it showed in Figures 6A–H, high expression of N₃-specific DEGs were all associated with statistically significant shorter OS (ARHGDIA, p = 0.0014; PDLIM7, p = 0.0043; PTTG1IP, p = 0.014; IL1A, p = 0.018; CCL20, p = 0.025; ABHD5, p = 0.031; SERPINB9, p = 0.037; IRAK2, p = 0.038), indicating that N₃ subtype might be pro-tumorigenic predicting inferior OS in advanced NSCLC.

Figure 6 N3 subtype was considered as predictors for inferior OS in advanced NSCLC. High expression of eight N3-specific DEGs, (A) ARHGDIA, (B) PDLIM7, (C) PTTG1IP, (D) IL1A, (E) CCL20, (F) ABHD5, (G) SERPINB9, and (H) IRAK2, correlated with inferior OS among advanced NSCLC. Abbreviations: OS, overall survival; NSCLC, non-small cell lung cancer; DEG, differentially expressed genes.

Discussion

It has been realized that neutrophils can not only play a key role in the inflammatory response and tissue damage repair but also infiltrate into tumors and participate in the adaptive immune response.^{2–4,11–14,30–33} Collectively, we believe that it is necessary to profile neutrophils at single-cell resolution for further understanding of this specific immune cell population. Therefore, in this study, we analyzed the TAN population from nine NSCLC biopsy specimens using unbiased scRNA-seq to depict the transcriptional landscape of TANs in human advanced NSCLC.

In our study, seven diverse TAN subtypes were defined. In addition, the intratumoral and intertumoral heterogeneity was again verified in terms of the composition of TAN subtypes across specimens. Based on the expression of top DEGs, the N₃ cluster was considered inflammatory phenotype expressing genes encoding multiple chemotactic cytokines, and significantly accumulated in LUSCs. Meanwhile, high expression of N₃-specific DEGs was correlated with inferior OS, indicating that N₃ might be pro-tumorigenic. In addition, N₁ and N₅ clusters were deemed to be well differentiated and mature neutrophils on the basis of CXCR2 expression and pseudo-time analyses, and both accounted for relatively high proportions in LUAD specimens.

Certain genes connected with neutrophil differentiation were discovered. Thereinto, DUSP1, FOS, MNDA, RGS2, and S100A8/9 showed an upward trend with TAN differentiation and maturation like CXCR2 did; yet MMP12 and PLIN2 presented reversely a trend of descending in this scenario. Moncho-Amor et al³⁴ revealed that DUSP1 overexpression was involved in angiogenesis, invasion, and metastasis in NSCLC and downregulating DUSP1 in tumor-bearing mice led to attenuation of tumor growth. FOS was proved by Shibahara et al³⁵ to induce PD-L2 expression contributing to tumor immune escape in NSCLC. MNDA and S100A8/9 shared correlative expressions in MAPK and PI3K/AKT signaling pathways and together participated in immune modulation in Type 1 Diabetes³⁶ and myeloproliferative neoplasms.³⁷ Kim et al³⁸ observed that loss of Med1/TRAP220 promoted the invasion and metastasis abilities of human NSCLC cells via RGS2 upregulation. Hofmann et al³⁹ reported that MMP12 protein was only detected in tumor samples and correlated with recurrence and metastasis in NSCLC; moreover, MMP12 knockdown could suppress LUAD growth and invasion.⁴⁰ PLIN2 was examined as an independent predictive factor for OS in NSCLC.⁴¹ In addition, Pang et al⁴² established a six-neutrophil-differentiation-related-gene-based prognostic risk model consisting of MS4A7, CXCR2, CSRNP1, RETN, CD177, and LUCAT1 by integrating scRNA-seq data, and the model was verified to be suitable for predicting the prognosis and immunotherapy response of NSCLC patients. In other words, neutrophil differentiation-related genes might serve as both prognostic/predictive markers and potential therapeutic targets.

Previous studies have likewise defined several distinct neutrophil subtypes. However, most studies were conducted using mouse models.^{16,28,43–45} Salcher et al⁴⁶ proposed subpopulations of tissue-resident neutrophils including TANs and normal adjacent tissue-associated neutrophils (NANs) in NSCLC, which were further divided into four TAN subsets (TAN-1 to TAN-4) and three NAN subsets (NAN-1 to NAN-3). Zilionis et al¹⁹ revealed five subtypes of TANs in human lung cancer samples, which were further classified into three categories based on DEGs: 1) neutrophils expressing canonical neutrophil markers (hN₁), like N₁ in our study; 2) tumor-specific neutrophils that promote tumor growth (hN₅), resembling our N₃ cluster; and 3) a neutrophil subset, similar to our N₅ subtype, with a gene signature of type I IFN response (hN₂).

It was revealed that TANs and multiple macrophage subtypes interplay most closely through chemokine signaling pathways, particularly the axis of CXCL1/2/3/8-CXCR2. Specifically, Glu-Leu-Arg-positive CXCL1/2/3/8 recruit CXCR2⁺ myeloid-derived suppressor cells within TME facilitating tumor growth and metastasis, which is correlated with poor outcomes in cancers.^47,48 Intriguingly, macrophages behave bifacially within TME,⁴⁹ which is similar to neutrophils, functioning as anti-tumorigenic or pro-tumorigenic factors. This inspires us that the regulation of tumor immunity by neutrophils is a highly complex process, in which various unsolved mechanisms might be involved. In light of published literature, we posit that it might be attributed to the diversity and plasticity of neutrophils in TME during anti-cancer therapy and tumor progression.

In 2009, Fridlender et al¹⁰ first proposed two phenotypes of TANs, N1 TAN (antitumor effect) and N2 TAN (protumor effect) and proved that the transformation of TANs to the N2 phenotype was modulated by the overexpression of TGF-β on tumor cells. In addition, the conversion of TANs to the N1 phenotype under induction of type I IFN was later verified.⁹ Thus, the phenotypic switches between N1 and N2 have indicated an antagonistic signaling pathway between TGF-β and type I IFN, and meanwhile revealed the plasticity of neutrophils in TME. Sagiv et al¹⁵ identified different circulating neutrophils, LDN (protumor effect) and HDN (antitumor effect), on a density gradient. Although the switch from HDN to LDN was similarly dominated by TGF-β and type I IFN along with tumor progression, LDN did consist of both immature and mature neutrophils through morphological observation and functional verification.

The possible molecular mechanisms of neutrophils modulating tumor immunity, especially the immunosuppressive effects, could be divided into direct suppression and indirect suppression. The direct ways to promote tumor progression, metastasis, and angiogenesis are as follows: 1) Reactive oxygen or nitrogen species released by neutrophils can lead to DNA oxidative damage and genetic instability;⁵⁰ 2) The released enzymes produced by neutrophils,^51–53 such as neutrophil collagenase (ie, MMP-8) and gelatinase B (ie, MMP-9), help rebuild the extracellular matrix to facilitate angiogenesis and stabilize integrins to enhance migration and invasion;⁵¹ 3) Modifying certain molecules, such as a downstream molecule of CSF3R signaling pathway, nicotinamide phosphoribosyltransferase, can stimulate pro-angiogenic activity.⁵⁴ Moreover, the indirect ways include 1) Neutrophils secrete a variety of cytokines in TME to regulate other immune cells, and ultimately promote tumor progression. For example, CCL17 and Arginase 1 released by neutrophils can recruit regulatory T cells and compromise T cell responses, respectively.^55,56 2) Cytokines or growth factors, such as IL-8, secreted by neutrophils may boost exudation and metastasis of tumor cells.^57,58 3) Neutrophils can also secrete a substance called oncostatin M, which has a variety of chemotactic and pro-angiogenic effects.^59,60

There are some limitations in our study. First is the limitation of our small sample size, which leads to the infeasibility of analyzing the correlation between TAN subtypes and clinical outcomes. In addition, provided immune cells or neutrophils were enriched during sample preparation, many more neutrophils could be collected for further refining TAN heterogeneity. Neutrophils are a group of immune cells with high heterogeneity and plasticity within TME regulating tumor immunity. Therefore, a more refined classification of neutrophils is an urgent need, which is conducive to the exploitation of therapeutic strategies targeting immunosuppressive neutrophils.

Conclusion

In this study, the functional, spatial, and temporal heterogeneity of TANs was verified in NSCLC at single-cell transcriptome resolution. Seven distinct TAN subtypes were defined, and each had a distinct DEG profile. Among these, the N₃ subtype was associated with an inferior response to anti-tumor treatment. Neutrophil differentiation/maturation-related genes were also discovered, which may further provide novel targets for immunotherapy. Limitations exist due to the small sample size, and large-scale analysis with more enriched neutrophils is warranted in the future to fully understand the multifaceted roles of TANs in NSCLC.

Abbreviations

DEG, differentially expressed gene; DNA, deoxyribonucleic acid; dNLR, derived neutrophil-to-lymphocyte ratio; ECOG PS, Eastern Cooperative Oncology Group performance status; EGFR, epidermal growth factor receptor; fDC, follicular dendritic cells; GO, Gene Ontology; HDN, high-density neutrophil; HSP, Heat Shock Protein; HVG, highly variable gene; IFN, interferon; KEGG, Kyoto Encyclopedia of Genes and Genome; LDN, low-density neutrophil; LUAD, lung adenocarcinoma; LUSC, lung squamous cell carcinoma; MP, mononuclear phagocyte; NAN, normal adjacent tissue-associated neutrophil; NLR, neutrophil-to-lymphocyte ratio; NSCLC, non-small cell lung cancer; OS, overall survival; RNA, ribonucleic acid; scRNA-seq, single-cell RNA-sequencing; TAN, tumor-associated neutrophils; TGF, transforming growth factor; TME, tumor microenvironment; UMAP, Uniform Manifold Approximation and Projection.

Data Sharing Statement

The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request.

Acknowledgments

Abstract of this paper was presented at the ESMO Immuno-Oncology Congress 2022 as a poster presentation with interim findings. The poster’s abstract was published in “Poster Abstracts” in Annals of Oncology: https://www.esmoiotech.org/article/S2590-0188(22)00251-9/fulltext.

Ethics Approval and Informed Consent

This study was performed in line with the principles of the Declaration of Helsinki. The Ethical Committee of Shanghai Pulmonary Hospital Affiliated to Tongji University approved this study (No. K18-089-1). Informed consent was obtained from all individual participants included in the study.

Funding

This work was supported by grants from National Natural Science Foundation of China (No. 81902314), Natural Science Foundation of Shanghai (No. 20ZR1447100), Clinical Research Plan of Shanghai Hospital Development Center (No. SHDC2020CR4001), Shanghai Municipal Key Clinical Discipline Construction Project (Respiratory Medicine), and Shanghai Municipal Multidisciplinary Collaboration Capacity Construction Project for Major Diseases.

Disclosure

Fengying Wu serves as the associate editor-in-chief of Lung Cancer: Targets and Therapy. The other authors have no conflicts of interest to declare for this work.

References

1. Adrover JM, Nicolas-Avila JA, Hidalgo A. Aging: a temporal dimension for neutrophils. Trends Immunol. 2016;37(5):334–345. doi:10.1016/j.it.2016.03.005

2. Galdiero MR, Varricchi G, Loffredo S, Mantovani A, Marone G. Roles of neutrophils in cancer growth and progression. J Leukoc Biol. 2018;103(3):457–464. doi:10.1002/JLB.3MR0717-292R

3. Nemeth T, Sperandio M, Mocsai A. Neutrophils as emerging therapeutic targets. Nat Rev Drug Discov. 2020;19(4):253–275. doi:10.1038/s41573-019-0054-z

4. Zhou J, Nefedova Y, Lei A, Gabrilovich D. Neutrophils and PMN-MDSC: their biological role and interaction with stromal cells. Semin Immunol. 2018;35:19–28. doi:10.1016/j.smim.2017.12.004

5. Valero C, Lee M, Hoen D, et al. Pretreatment neutrophil-to-lymphocyte ratio and mutational burden as biomarkers of tumor response to immune checkpoint inhibitors. Nat Commun. 2021;12(1):729. doi:10.1038/s41467-021-20935-9

6. Eruslanov EB, Bhojnagarwala PS, Quatromoni JG, et al. Tumor-associated neutrophils stimulate T cell responses in early-stage human lung cancer. J Clin Invest. 2014;124(12):5466–5480. doi:10.1172/JCI77053

7. Albanesi M, Mancardi DA, Jonsson F, et al. Neutrophils mediate antibody-induced antitumor effects in mice. Blood. 2013;122(18):3160–3164. doi:10.1182/blood-2013-04-497446

8. Mishalian I, Bayuh R, Levy L, Zolotarov L, Michaeli J, Fridlender ZG. Tumor-associated neutrophils (TAN) develop pro-tumorigenic properties during tumor progression. Cancer Immunol Immunother. 2013;62(11):1745–1756. doi:10.1007/s00262-013-1476-9

9. Andzinski L, Kasnitz N, Stahnke S, et al. Type I IFNs induce anti-tumor polarization of tumor associated neutrophils in mice and human. Int J Cancer. 2016;138(8):1982–1993. doi:10.1002/ijc.29945

10. Fridlender ZG, Sun J, Kim S, et al. Polarization of tumor-associated neutrophil phenotype by TGF-beta: “N1” versus “N2” TAN. Cancer Cell. 2009;16(3):183–194. doi:10.1016/j.ccr.2009.06.017

11. Giese MA, Hind LE, Huttenlocher A. Neutrophil plasticity in the tumor microenvironment. Blood. 2019;133(20):2159–2167. doi:10.1182/blood-2018-11-844548

12. Mantovani A. The yin-yang of tumor-associated neutrophils. Cancer Cell. 2009;16(3):173–174. doi:10.1016/j.ccr.2009.08.014

13. Mukaida N, Sasaki SI, Baba T. Two-faced roles of tumor-associated neutrophils in cancer development and progression. Int J Mol Sci. 2020;21(10):3457. doi:10.3390/ijms21103457

14. Ng LG, Ostuni R, Hidalgo A. Heterogeneity of neutrophils. Nat Rev Immunol. 2019;19(4):255–265. doi:10.1038/s41577-019-0141-8

15. Sagiv JY, Michaeli J, Assi S, et al. Phenotypic diversity and plasticity in circulating neutrophil subpopulations in cancer. Cell Rep. 2015;10(4):562–573. doi:10.1016/j.celrep.2014.12.039

16. Xie X, Shi Q, Wu P, et al. Single-cell transcriptome profiling reveals neutrophil heterogeneity in homeostasis and infection. Nat Immunol. 2020;21(9):1119–1133. doi:10.1038/s41590-020-0736-z

17. Wu F, Fan J, He Y, et al. Single-cell profiling of tumor heterogeneity and the microenvironment in advanced non-small cell lung cancer. Nat Commun. 2021;12(1):2540. doi:10.1038/s41467-021-22801-0

18. Mistry P, Nakabo S, O’Neil L, et al. Transcriptomic, epigenetic, and functional analyses implicate neutrophil diversity in the pathogenesis of systemic lupus erythematosus. Proc Natl Acad Sci U S A. 2019;116(50):25222–25228. doi:10.1073/pnas.1908576116

19. Zilionis R, Engblom C, Pfirschke C, et al. Single-cell transcriptomics of human and mouse lung cancers reveals conserved myeloid populations across individuals and species. Immunity. 2019;50(5):1317–1334 e10. doi:10.1016/j.immuni.2019.03.009

20. Gate D, Saligrama N, Leventhal O, et al. Clonally expanded CD8 T cells patrol the cerebrospinal fluid in Alzheimer’s disease. Nature. 2020;577(7790):399–404. doi:10.1038/s41586-019-1895-7

21. Madissoon E, Wilbrey-Clark A, Miragaia RJ, et al. scRNA-seq assessment of the human lung, spleen, and esophagus tissue stability after cold preservation. Genome Biol. 2019;21(1):1. doi:10.1186/s13059-019-1906-x

22. Sun Z, Chen L, Xin H, et al. A Bayesian mixture model for clustering droplet-based single-cell transcriptomic data from population studies. Nat Commun. 2019;10(1):1649. doi:10.1038/s41467-019-09639-3

23. Stewart BJ, Ferdinand JR, Young MD, et al. Spatiotemporal immune zonation of the human kidney. Science. 2019;365(6460):1461–1466. doi:10.1126/science.aat5031

24. Yu G, Wang LG, Han Y, He QY. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS. 2012;16(5):284–287. doi:10.1089/omi.2011.0118

25. Trapnell C, Cacchiarelli D, Grimsby J, et al. The dynamics and regulators of cell fate decisions are revealed by pseudotemporal ordering of single cells. Nat Biotechnol. 2014;32(4):381–386. doi:10.1038/nbt.2859

26. Vento-Tormo R, Efremova M, Botting RA, et al. Single-cell reconstruction of the early maternal-fetal interface in humans. Nature. 2018;563(7731):347–353. doi:10.1038/s41586-018-0698-6

27. Shannon P, Markiel A, Ozier O, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003;13(11):2498–2504. doi:10.1101/gr.1239303

28. Evrard M, Kwok IWH, Chong SZ, et al. Developmental analysis of bone marrow neutrophils reveals populations specialized in expansion, trafficking, and effector functions. Immunity. 2018;48(2):364–379 e8. doi:10.1016/j.immuni.2018.02.002

29. Zhu YP, Padgett L, Dinh HQ, et al. Identification of an early unipotent neutrophil progenitor with pro-tumoral activity in mouse and human bone marrow. Cell Rep. 2018;24(9):2329–2341 e8. doi:10.1016/j.celrep.2018.07.097

30. Dumitru CA, Moses K, Trellakis S, Lang S, Brandau S. Neutrophils and granulocytic myeloid-derived suppressor cells: immunophenotyping, cell biology and clinical relevance in human oncology. Cancer Immunol Immunother. 2012;61(8):1155–1167. doi:10.1007/s00262-012-1294-5

31. Hassani M, Hellebrekers P, Chen N, et al. On the origin of low-density neutrophils. J Leukoc Biol. 2020;107(5):809–818. doi:10.1002/JLB.5HR0120-459R

32. Tecchio C, Cassatella MA. Neutrophil-derived chemokines on the road to immunity. Semin Immunol. 2016;28(2):119–128. doi:10.1016/j.smim.2016.04.003

33. Tsai CY, Hsieh SC, Liu CW, et al. Cross-Talk among polymorphonuclear neutrophils, immune, and non-immune cells via released cytokines, granule proteins, microvesicles, and neutrophil extracellular trap formation: a novel concept of biology and pathobiology for neutrophils. Int J Mol Sci. 2021;22(6):3119. doi:10.3390/ijms22063119

34. Moncho-Amor V, Ibanez de Caceres I, Bandres E, et al. DUSP1/MKP1 promotes angiogenesis, invasion and metastasis in non-small-cell lung cancer. Oncogene. 2011;30(6):668–678. doi:10.1038/onc.2010.449

35. Shibahara D, Tanaka K, Iwama E, et al. Intrinsic and extrinsic regulation of PD-L2 expression in oncogene-driven non-small cell lung cancer. J Thorac Oncol. 2018;13(7):926–937. doi:10.1016/j.jtho.2018.03.012

36. Jin Y, Sharma A, Carey C, et al. The expression of inflammatory genes is upregulated in peripheral blood of patients with type 1 diabetes. Diabetes Care. 2013;36(9):2794–2802. doi:10.2337/dc12-1986

37. Cokic VP, Mossuz P, Han J, et al. Microarray and proteomic analyses of myeloproliferative neoplasms with a highlight on the mTOR signaling pathway. PLoS One. 2015;10(8):e0135463. doi:10.1371/journal.pone.0135463

38. Kim HJ, Roh MS, Son CH, et al. Loss of Med1/TRAP220 promotes the invasion and metastasis of human non-small-cell lung cancer cells by modulating the expression of metastasis-related genes. Cancer Lett. 2012;321(2):195–202. doi:10.1016/j.canlet.2012.02.009

39. Hofmann HS, Hansen G, Richter G, et al. Matrix metalloproteinase-12 expression correlates with local recurrence and metastatic disease in non-small cell lung cancer patients. Clin Cancer Res. 2005;11(3):1086–1092. doi:10.1158/1078-0432.1086.11.3

40. Lv FZ, Wang JL, Wu Y, Chen HF, Shen XY. Knockdown of MMP12 inhibits the growth and invasion of lung adenocarcinoma cells. Int J Immunopathol Pharmacol. 2015;28(1):77–84. doi:10.1177/0394632015572557

41. Tang D, Zhao YC, Liu H, et al. Potentially functional genetic variants in PLIN2, SULT2A1 and UGT1A9 genes of the ketone pathway and survival of nonsmall cell lung cancer. Int J Cancer. 2020;147(6):1559–1570. doi:10.1002/ijc.32932

42. Pang J, Yu Q, Chen Y, Yuan H, Sheng M, Tang W. Integrating Single-cell RNA-seq to construct a Neutrophil prognostic model for predicting immune responses in non-small cell lung cancer. J Transl Med. 2022;20(1):531. doi:10.1186/s12967-022-03723-x

43. Giladi A, Paul F, Herzog Y, et al. Single-cell characterization of haematopoietic progenitors and their trajectories in homeostasis and perturbed haematopoiesis. Nat Cell Biol. 2018;20(7):836–846. doi:10.1038/s41556-018-0121-4

44. Muench DE, Olsson A, Ferchen K, et al. Mouse models of neutropenia reveal progenitor-stage-specific defects. Nature. 2020;582(7810):109–114. doi:10.1038/s41586-020-2227-7

45. Olsson A, Venkatasubramanian M, Chaudhri VK, et al. Single-cell analysis of mixed-lineage states leading to a binary cell fate choice. Nature. 2016;537(7622):698–702. doi:10.1038/nature19348

46. Salcher S, Sturm G, Horvath L, et al. High-resolution single-cell atlas reveals diversity and plasticity of tissue-resident neutrophils in non-small cell lung cancer. Cancer Cell. 2022;40(12):1503–1520. doi:10.1016/j.ccell.2022.10.008

47. Wang D, Sun H, Wei J, Cen B, DuBois RN. CXCL1 is critical for premetastatic niche formation and metastasis in colorectal cancer. Cancer Res. 2017;77(13):3655–3665. doi:10.1158/0008-5472.CAN-16-3199

48. Zhang H, Ye YL, Li MX, et al. CXCL2/MIF-CXCR2 signaling promotes the recruitment of myeloid-derived suppressor cells and is correlated with prognosis in bladder cancer. Oncogene. 2017;36(15):2095–2104. doi:10.1038/onc.2016.367

49. Christofides A, Strauss L, Yeo A, Cao C, Charest A, Boussiotis VA. The complex role of tumor-infiltrating macrophages. Nat Immunol. 2022;23:1148–1156. doi:10.1038/s41590-022-01267-2

50. Gungor N, Knaapen AM, Munnia A, et al. Genotoxic effects of neutrophils and hypochlorous acid. Mutagenesis. 2010;25(2):149–154. doi:10.1093/mutage/gep053

51. Ardi VC, Kupriyanova TA, Deryugina EI, Quigley JP. Human neutrophils uniquely release TIMP-free MMP-9 to provide a potent catalytic stimulator of angiogenesis. Proc Natl Acad Sci U S A. 2007;104(51):20262–20267. doi:10.1073/pnas.0706438104

52. Houghton AM, Rzymkiewicz DM, Ji H, et al. Neutrophil elastase-mediated degradation of IRS-1 accelerates lung tumor growth. Nat Med. 2010;16(2):219–223. doi:10.1038/nm.2084

53. Park J, Wysocki RW, Amoozgar Z, et al. Cancer cells induce metastasis-supporting neutrophil extracellular DNA traps. Sci Transl Med. 2016;8(361):361ra138. doi:10.1126/scitranslmed.aag1711

54. Pylaeva E, Harati MD, Spyra I, et al. NAMPT signaling is critical for the proangiogenic activity of tumor-associated neutrophils. Int J Cancer. 2019;144(1):136–149. doi:10.1002/ijc.31808

55. Mishalian I, Bayuh R, Eruslanov E, et al. Neutrophils recruit regulatory T-cells into tumors via secretion of CCL17–a new mechanism of impaired antitumor immunity. Int J Cancer. 2014;135(5):1178–1186. doi:10.1002/ijc.28770

56. Yachimovich-Cohen N, Even-Ram S, Shufaro Y, Rachmilewitz J, Reubinoff B. Human embryonic stem cells suppress T cell responses via arginase I-dependent mechanism. J Immunol. 2010;184(3):1300–1308. doi:10.4049/jimmunol.0804261

57. Chen MB, Hajal C, Benjamin DC, et al. Inflamed neutrophils sequestered at entrapped tumor cells via chemotactic confinement promote tumor cell extravasation. Proc Natl Acad Sci U S A. 2018;115(27):7022–7027. doi:10.1073/pnas.1715932115

58. Huh SJ, Liang S, Sharma A, Dong C, Robertson GP. Transiently entrapped circulating tumor cells interact with neutrophils to facilitate lung metastasis development. Cancer Res. 2010;70(14):6071–6082. doi:10.1158/0008-5472.CAN-09-4442

59. Kerfoot SM, Raharjo E, Ho M, et al. Exclusive neutrophil recruitment with oncostatin M in a human system. Am J Pathol. 2001;159(4):1531–1539. doi:10.1016/S0002-9440(10)62538-2

60. Queen MM, Ryan RE, Holzer RG, Keller-Peck CR, Jorcyk CL. Breast cancer cells stimulate neutrophils to produce oncostatin M: potential implications for tumor progression. Cancer Res. 2005;65(19):8896–8904. doi:10.1158/0008-5472.CAN-05-1734