Single-Cell Profiling of Tumor-Associated Neutrophils in Advanced Non-Small Cell Lung Cancer

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.05 and of average log expression >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.