Skip to main content

Integrated analysis of single-cell RNA-seq and bulk RNA-seq unravels the molecular feature of tumor-associated neutrophils of head and neck squamous cell carcinoma

Abstract

Background

Head and neck squamous cell carcinoma (HNSCC) is a lethal malignancy with a high recurrence and distant metastasis rate, posing significant challenges to patient prognosis. Recent studies suggest that tumor-associated neutrophils (TANs) can modulate immune cell infiltration and influence tumor initiation and progression. However, the potential clinical significance of TANs in HNSCC remains insufficiently explored.

Methods

TANs-specific marker genes were identified via single-cell sequencing data from HNSCC. Based on data from The Cancer Genome Atlas (TCGA), a prognostic risk model was constructed using TANs cell marker genes, and the model was validated with data from the Gene Expression Omnibus (GEO) database. The associations between the TANs signature and clinical characteristics, functional pathways, immune cell infiltration, immune checkpoint expression, and responses to immunotherapy and chemotherapy, were then investigated. Cell counting kit‐8(CCK-8), Transwell, and wound healing assays were conducted to assess the functional role of TANs marker molecules.

Results

TANs characteristic genes were identified from single-cell sequencing data from HNSCC patients. On the basis of these characteristic genes, a tumor-associated neutrophils-associated signature (NRS) was developed and validated across internal and external cross-platform cohorts through comprehensive procedures. The NRS demonstrated robust and reliable performance in predicting overall survival. Additionally, patients with a low NRS showed enhanced immune cell infiltration, active lipid metabolism, and increased sensitivity to immunotherapy. In contrast, patients with a high NRS exhibited poor prognostic outcomes, advanced clinical stages, and significant associations with HNSCC metastasis and progression. Furthermore, we identified a TANs-associated biomarker, OLR1, and validated that OLR1 promotes HNSCC proliferation, invasion, and migration through CCK-8, Transwell invasion, and wound healing assays.

Conclusion

This study has developed a promising TANs-based tool that may aid in personalized treatment and prognostic management for patients with HNSCC.

Peer Review reports

Introduction

Head and neck cancers rank as the sixth most prevalent malignant tumors globally, with head and neck squamous cell carcinomas (HNSCC) originating from the nasal cavity, pharynx, larynx, and oral cavity, accounting for over 90% of head and neck malignancies [1]. These malignancies demonstrate considerable heterogeneity and aggressiveness [2]. Projections indicate that by 2030, the global burden of HNSCC will rise by 30% [3], and by 2040, there will be nearly 600,000 new cases of HNSCC diagnosed each year [4]. Traditional HNSCC treatments include a combination of surgical intervention, radiotherapy, and chemotherapy, with recent progress in tumor immunotherapy. For instance, pembrolizumab has been approved as a first-line treatment for advanced unresectable HNSCC [5]. Nevertheless, the therapeutic efficacy of immune checkpoint inhibitors in the majority of HNSCC patients remains suboptimal [6], with the current five-year survival rate for patients with middle-to-advanced stage HNSCC persisting at less than 50% [7]. Developing valuable biomarkers and therapeutic targets based on the current situation is of significant clinical importance for effective diagnosis, treatment, and improving prognosis.

The stromal cells and immune cells in the tumor microenvironment (TME) are the main cellular components mediating immune tolerance and immune escape [8]. The role of neutrophils within the TME has not been as extensively explored as that of other myeloid cell types, like dendritic cells (DCs) and macrophages [9, 10]. However, growing evidence supports the notion that tumor-associated neutrophils(TANs) are critical mediators of the immunosuppressive environment in cancer development and are driving factors in tumor progression. For instance, research conducted by Sebastien et al. illustrates that TANs inhibit T cell proliferation and natural killer (NK) cell cytotoxicity by releasing reactive oxygen species (ROS) and inducible nitric oxide synthase (iNOS), along with arginase 1 (ARG1) activity [11,12,13]. Neutrophil-released MMP- 9 degrades the extracellular matrix, releasing vascular endothelial growth factor (VEGF) and promoting angiogenesis to facilitate tumor growth [14]. Additionally, the damage to tumor-infiltrating neutrophils on the cytotoxic T lymphocyte (CTL) activity within the tumor bed is an important reason for the resistance to immunotherapy [13]. Therefore, targeting or manipulating neutrophils may render tumors more sensitive to conventional immunotherapy, with several neutrophil-targeting therapies for hepatocellular carcinoma (and other malignancies) currently entering clinical trials [15,16,17]. Considering the significant role of TANs in tumor progression and response to immunotherapy, developing neutrophil-related features to predict survival outcomes and immunological characteristics in HNSCC patients holds significant value.

The rapid advancement of high-throughput sequencing technology has significantly influenced research in the field of biology. Many studies have leveraged sequencing data to predict clinical outcomes across various cancer types [18]. However, TME is a complex environment with high heterogeneity, and traditional transcriptomics studies may overlook biological differences between different cell types. Compared to conventional RNA sequencing, single-cell RNA sequencing (scRNA-seq) technology facilitates the identification of heterogeneity in tumor cells and stromal cells from a cellular perspective, distinguishing gene expression characteristics of different cell types, thereby accurately identifying key regulatory genes in cell subpopulations [19]. Currently, the role of neutrophils in the TME of HNSCC and in immune therapy responses is under-researched. Therefore, this study aims to elucidate the impact of TANs on HNSCC prognosis by integrating scRNA-seq and transcriptomics, evaluate their regulation of the immune microenvironment in HNSCC, and validate the functional roles of their marker molecules in HNSCC through in vitro experiments.

Material and method

Data sources used for analysis

The study employed HNSCC scRNA-seq data obtained from the GSE181919 database, comprising 20 tumor samples and 9 control samples [20]. For transcriptome data, we initially acquired 44 control samples and 519 HNSCC samples from the TCGA-HNSCC cohort available in the TCGA database. Additionally, we utilized two cohorts of HNSCC patients (GSE65858 and GSE41613) from the GEO database as external validation cohorts [21, 22].

Single cell sequencing data processing

The Seurat software package (version 4.3.0) was employed for further data reduction and cell clustering analysis on single-cell sequencing data derived from GSE181919 [23]. Single-cell gene expression profiles were filtered to exclude cells exhibiting more than 10% mitochondrial gene content or fewer than 200 transcripts per cell, as these are considered low quality. Subsequently, the RunPCA function was applied to conduct linear dimensionality reduction through principal component analysis (PCA), and the RunTSNE function was executed to generate visual representations of the clusters. Afterward, the"SingleR"package was employed for cluster annotation, thus enabling the identification of distinct cell types [24]. The"FindAllMarkers"function was employed to discern genes associated with TANs [25], defined as those genes exhibiting a |log2(FC)| value greater than 0.585 and a P-value less than 0.05.

The development and validation of prognostic risk models utilizing TANs marker genes

Univariate Cox regression analysis was conducted on TANs marker genes to identify prognostic genes. Subsequently, LASSO regression was employed for further gene selection to prevent model overfitting. The findings from the multivariate Cox regression analysis were employed to develop a risk score model, which subsequently facilitated the stratification of patients into high-risk and low-risk categories according to the median risk score.

$$\text{Risk score}={\text{expr}}_{1}\text{ x }{\upbeta }_{1}+{\text{expr}}_{2}\text{x }{\upbeta }_{2}+\dots +{\text{expr}}_{\text{n}}\text{ x }{\upbeta }_{\text{n}}$$

The expression of key genes, represented by expr, and the corresponding regression coefficient β are utilized in the multivariate genetic Cox regression analysis.

The time-dependent area under the ROC curve was utilized to assess the model's prognostic capability. Kaplan–Meier survival analysis was conducted to estimate prognostic differentiation. Two GEO datasets were independently validated, and the"rms"R package was utilized to construct the nomogram. Calibration curves were constructed to assess the nomogram's predictive accuracy.

Function analysis

In order to elucidate the underlying biological functions associated with each subtype, we performed the gene set enrichment analysis (GSEA) [26], considering terms with a P-value of less than 0.05 as statistically significant. A total of 9,997 gene sets were obtained from the Molecular Signatures Database (MSigDB) and subsequently transformed into a pathway enrichment score matrix utilizing the gene set variation analysis (GSVA) algorithm [27]. The Edge package was employed to further identify characteristic pathways within the two NRS groups, applying thresholds of a P-value less than 0.05 and a log2 (FC) value exceeding 0.2, respectively.

Immune landscape

The infiltration abundance of 28 immune cells in bulk tumor samples was quantified using the single-sample gene set enrichment analysis (ssGSEA) method [28]. The evaluation of the correlation between risk signature and the TME involved the collection of 27 key immune checkpoints, including molecules from the TNF superfamily [29], B7-CD28 family [30], and other relevant factors [31]. The"estimate"R package was employed to compute stromal scores, immune scores, and estimate scores in order to assess disparities in the TME among HNSCC patients. The Cancer Immune Cycle (CIC) [32] was proposed by Karasaki et al. as a systematic approach to assessing the immune status of different subsets.

Immunotherapy and chemotherapy prediction

The efficacy of NRS in predicting the response to immunotherapy was assessed using the Tumor Immune Dysfunction and Exclusion (TIDE) score. Simultaneously, three immunotherapy cohorts (GSE140901, GSE126044, and the IMvigor210 cohort) were enlisted to further validate the assessment of immunotherapeutic responses using the Subclass Mapping (SubMap) algorithm [33]. To further investigate the relationship between NRS characteristics and chemotherapy, the Genomics of Drug Sensitivity in Cancer (GDSC) [34] database was employed to determine the half-maximal inhibitory concentrations (IC50) [35] of frequently utilized clinical chemotherapy medications for HNSCC. The pharmacogenomic datasets of the Profiling Relative Inhibition Simultaneously in Mixtures (PRISM) and Cancer Therapeutics Response Portal (CTRP) [36] encompass extensive drug screening and transcriptomic data, spanning numerous cancer cell lines, aiding in the development of predictive models for drug response. Drug sensitivity was assessed using the area under the dose–response curve (AUC), where lower AUC values signify increased sensitivity to treatment. Differential analyses of drug sensitivity were then conducted based on various NRS characteristics to identify potential therapeutic targets.

Random Forest assessment

In this study, the Random Forest (RF) algorithm was used for gene screening, aiming to further identify key genes closely related to prognosis from the model. A random forest model was constructed, in which each decision tree was trained based on a randomly selected subset of features. During the construction process, the feature selection of each tree was random, which enabled the random forest to identify genes that contributed greatly to drug response. The model ranked the genes according to their importance scores and selected genes with high importance scores as candidate genes for subsequent analysis.

BEST

BEST (https://rookieutopia.com/app_direct/BEST/) is an online analysis platform that includes data from the TCGA and GEO databases and provides visualized analysis results. Using BEST, we analyzed the protein expression of OLR1, survival analysis, immune infiltration, the correlation of immune cells, and survival analysis in immunotherapy cohorts. We utilized the BEST Web database to perform immune infiltration analysis using four algorithms: CIBERSORT, CIBERSORT_ABS, Quantiseq, and MCPcounter [37, 38]. The detailed methodology has been compiled in the “supplementary methods” file.

Nomogram method

The Nomogram method is used to predict survival. First, each clinical feature is assigned a corresponding score according to its weight. Then, the scores of each clinical feature are added together according to the patient's specific situation to obtain a total risk score. Finally, based on the total score, the patient's survival probability in 1-, 3 –, and 5-year is predicted.

Cell culture and siRNA transfection

Human cell lines HOK, CAL- 27, SCC- 4, SCC- 25 and HL- 60 were acquired from Wuhan Pricella Biotechnology Co., Ltd (Wuhan, China). The HOK, CAL- 27, SCC- 4, and SCC- 25 cell lines were cultured in DMEM medium (Gibco, USA). The HL- 60 cell line was cultured in RPMI 1640 medium (Gibco, USA). All the cells were cultured in media supplemented with 10% fetal bovine serum (FBS) and 1% penicillin–streptomycin solution at 37 °C with 5% CO2. siRNA and Lipo3000 were mixed in serum-free medium and left to stand for 15 min before being added to the well designated for transfection. Cells were harvested 48 h later for subsequent procedures.

Quantitative real‐time PCR (qRT-PCR)

Total RNA was isolated from logarithmically growing cells using TRIZOL (Takara, Japan). After measuring RNA concentration, reverse transcription was performed to synthesize complementary DNA (cDNA) utilizing a qRT-PCR kit (Novoprotein Scientific Inc, China). qRT-PCR was performed on a PCR amplification instrument (Bio-Rad) using SYBR qPCR SuperMix Plus (Novoprotein Scientific Inc, China). The expression levels of all target genes were normalized to GAPDH and analyzed using the 2 − ΔΔCq method.

CCK‐8 assay

After cell digestion, cells were seeded onto a 96-well plate at a density of 3 × 103 cells per well. The culture solution was exchanged for a solution containing CCK- 8 reagent. After incubating at 37 °C for 3 h, the optical density (OD) was detected at 450 nm.

Western Blot assay

Protein samples isolated from cells were separated by gel electrophoresis and transferred to a PVDF membrane. After blocking the membrane with 5% non-fat milk was incubated overnight at 4 °C with the primary antibody on a shaker. Finally, the membrane underwent incubation with the secondary antibody at ambient temperature for a duration of 1 h. The chemiluminescent signal was generated using a high-sensitivity ECL substrate (Abbkine, China).

Wound healing assay

Cells were cultured in six-well plates until they achieved approximately 95% confluence. Subsequently, a sterile 10 µL pipette tip was employed to create a linear scratch on the plates'substrate, simulating a wound. The detached cells were then removed by washing with phosphate-buffered saline (PBS), and serum-free medium was introduced into the wells. The area of the scratch was monitored and documented through photography at the initial time point (0 h) and after 24 h.

Transwell migration and invasion assays

The cell migration capacity was evaluated utilizing the Transwell assay. Specifically, 2 × 104 cells were seeded into the upper chamber containing 200 µL of serum-free medium, while 500 µL of complete medium was introduced into the corresponding wells of the plate. Following a 48-h incubation period, the cells that had migrated to the lower membrane were fixed, stained, and subsequently quantified. The same experimental procedure was followed for invasion assays after pre-coating the chambers with Matrigel (BD Biosciences).

Chemotaxis assay

1 × 10^5 dHL- 60 cells were seeded into the upper chamber of a Transwell device, and different treatment media were added to the lower chamber. After incubation for 4 h in a cell culture incubator, dHL- 60 cells in the lower chamber were observed and counted using an optical microscope.

Statistical analysis

Survival graphs were developed via the Kaplan–Meier technique and evaluated using log-rank analysis. The Wilcoxon test was employed to contrast the two groups. The R software (version 4.3.0) was utilized for statistical evaluation and data depiction, with a p-value less than 0.05 deemed to hold statistical significance.

Results

Identification of feature genes related to TANs in HNSCC

After quality control and standardized processing, 32,223 cells were obtained based on ScRNA-seq data from GSE181919. As shown in Figure S1, the details of data preprocessing are demonstrated. After logarithmic normalization and dimensionality reduction, we identified 16 distinct clusters (Fig. 1A). Figure 1B shows the top five marker genes for each cell cluster. Then, each cluster was annotated and grouped into 10 cell types (Fig. 1C). Among the 16 clusters, cluster 3 was identified as a neutrophil subpopulation. Figures 1D and 1E illustrate the differences in cell distribution across all samples, and we calculated the percentage composition of each sample and cell cluster. Notably, neutrophils were markedly enriched in tumor tissues relative to normal tissues. Furthermore, we extracted neutrophil-related genes (NRGs) from the tumor samples of single-cell data for further analysis using the FindAllMarkers function, with filtering criteria of p-value < 0.05 and |log2(FC)|> 0.585.

Fig. 1
figure 1

The identification of neutrophil cluster according to scRNA data of HNSCC patients. (A) uMap plot shows the regrouping of HNSCC single cells into 16 separate clusters. (B) The expression patterns of the cell type-specific marker genes across the cell clusters. (C) uMap plot showing cell types recognized by marker genes. (D-E) Subgroups in cancer and adjacent tissue and proportion as well as cell number calculation

Development and validation of tumor-associated neutrophil-related riskscore

Firstly, 41 NRGs significantly associated with the clinical outcomes of HNSCC patients were identified through univariate Cox regression analysis (Fig. 2A). To mitigate the risk of overfitting in the model, the outcomes of the univariate Cox regression analysis were further refined using LASSO regression analysis (Fig. 2B and 2 C). According to lasso regression analysis, we identified 11 TANs genes that are most predictive of prognosis, which help us develop biomarkers that affect HNSC prognostic. The LASSO regression outcomes were examined using multivariate Cox regression, while 11 genes contributed to the creation of the tumor-associated neutrophil-related risk score model (NRS): RiskScore = (− 0.212 × CCR7) + (− 0.192 × MMP19) + (− 0.059 × LGALS9) + (− 0.055 × RNF144B) + (− 0.009 × LYZ) + (0.068 × HBEGF) + (0.082 × FTH1) + (0.102 × VSIG4) + (0.117 × OLR1) + (0.150 × THBS1) + (0.152 × TPP1). In the prognostic model, CCR7 exhibited the strongest protective regulatory effect. In contrast, TPP1 emerged as the most significant risk factor (Fig. 2D). The processed TCGA-HNSCC dataset, incorporating survival data from the TCGA database, was randomly partitioned into a training set and a test set in a 3:1 ratio. HNSCC patients were stratified into high NRS (H-NRS) and low NRS (L-NRS) groups according to the median risk score derived from the TCGA training set. Survival distribution analysis demonstrated a significant difference in survival rates between the two groups, with the H-NRS associated with a worse prognosis and a 5-year survival rate below 50% (Fig. 2E-F). Furthermore, the ROC curve results for both the training and test sets indicate the model demonstrates good validation performance (Fig. 2I-J), highlighting the significant potential of NRS in predicting clinical outcomes in HNSCC. Consistent with prior research findings, the external independent validation datasets (GSE41613, GSE65858) demonstrated that the H-NRS group exhibited poorer prognostic outcomes (Fig. 2G-H and 2 K-L). Taken together, the aforementioned findings demonstrate that our NRS is both robust and reproducible across cross-platform cohorts.

Fig. 2
figure 2

Construction and validation of neutrophil-related risk score using the TCGA cohort. (A) Forest plot of survival related neutrophil-related genes based on univariate Cox regression analysis. (B) The trajectory of each independent variable with lambda. (C) Plots of the produced coefficient distributions for the logarithmic series for parameter selection (lambda). (D) Visualization of gene-specific weights in the model score. (EH) Kaplan–Meier survival curve of the risk signature in TCGA training dataset (E), TCGA testing dataset (F), GSE65858 (G), and GSE41613 (H). (I-L) Time-dependent ROC curve of the risk signature in TCGA training dataset (I), TCGA testing dataset (J), GSE65858 (K), and GSE41613 (L)

A nomogram integrating NRS and clinical parameters for precision prediction

Subsequently, we performed univariate and multivariate Cox analyses, and the results indicated that compared with other common clinicopathological parameters, NRS was significantly associated with the overall survival (OS) of HNSCC patients. It might serve as an independent prognostic factor (Fig. 3A-B). The results of the decision curve analysis (DCA) curve indicated that NRS outperformed common clinical parameters in the prognostic prediction of HNSCC patients (Fig. 3D). To enhance the clinical validity of the risk profile, a column line graph predicting OS was constructed in TCGA-HNSCC using NRS and clinical data (Fig. 3C). The 1-, 3-, and 5-year survival rates of each patient can be predicted based on the total score of each patient's NRS score and clinical feature score. Calibration graphs for the nomogram at 1-, 3-, and 5-year intervals demo pRRopheticnstrated a strong alignment with the ideal predictive curve, indicating a high level of congruence between the projected OS rates and the actual observed data (Fig. 3E). The ROC curves further demonstrated that the column-line graph was superior for predicting prognosis (Fig. 3F). Subsequently, we compared the distribution of clinical features between the two both groups. (Fig. 3G). H-NRS was associated with advanced stage, T-stage, and N-stage, consistent with the poor prognosis of H-NRS. In contrast, L-NRS was associated with high microsatellite instability (MSI-H) levels.

Fig. 3
figure 3

Clinical Applications of NRS characteristics. A Univariate Cox analysis of risk score and clinicopathological parameters. B Multivariate Cox analysis of risk score and clinicopathological parameters. C Nomogram model integrating the risk score and clinical parameters was constructed. D DCA plots were created for comparisons between the risk score and the clinical and pathological variables. E Calibration curves for 1-, 3-, and 5-year of nomogram. F ROC curves for 1-, 3- and 5-year survival based on the nomogram's predictions. G Correlations of two NRS groups with clinical characteristics in the TCGA-HNSCC dataset

Biological functions related to the NRS

We further calculated the variations in 9,997 pathways using the GSVA algorithm to demonstrate the biological functional differences in HNSCC patients with different NRS characteristics (Fig. 4A). H-NRS patients are primarily enriched in various biological processes related to tumor progression, such as invasion, metastasis, and proliferation (e.g., angiogenesis, KRAS signaling, and epithelial-mesenchymal transition), which aligns with their poor prognostic outcomes. In contrast, L-NRS patients are mainly associated with lipid metabolism processes and immune pathways. GSEA analysis indicated that the ECM-receptor interaction, TGF-β signaling, WNT signaling, and JAK-STAT signaling pathways are significantly activated in the H-NRS group (Fig. 4B).The L-NRS group is enriched in B cell receptor signaling, T cell receptor signaling, and fatty acid metabolism pathways (Fig. 4C). Overall, the H-NRS group shows a stronger correlation with tumor-promoting biological processes, whereas the L-NRS group is characterized by tumor suppression and a heightened immune response.

Fig. 4
figure 4

Functional enrichment analysis between H-NRS and L-NRS patients. A Heatmaps of biological processes for two NRS groups in TCGA-HNSCC datasets. High and low ssGESA scores are represented in red and blue, respectively. Positive prediction of signatures is indicated in red and absence in white. (B-C) Gene set enrichment analysis (GSEA) between H-NRS (B) and L-NRS (C) groups

Immune characteristics related to the NRS

Prior research has shown that the composition of immune cells in HNSCC is heterogeneous, with key components continually affected by changes in the TME [39]. Neutrophils are essential components of TME to control tumor progression and treatment resistance [40, 41]. Therefore, to explore the correlation between NRS characteristics and the TME, we measured the abundance of 28 immune cells infiltrated using ssGSEA (Fig. 5A). The findings demonstrated that the majority of immune cells were predominantly infiltrated in the L-NRS group, especially those activating in anti-tumor immunity, such as activated CD4 T cells, CD8 T cells, and NK cells. Based on the ESTIMATE algorithm, the L-NRS features also showed high immune scores (Fig. 5C). These analyses suggest that L-NRS profile status may be associated with immune activation status. Furthermore, immune checkpoints, including proteins from the B7-CD28 family and TNF superfamily, exhibited elevated expression levels in the L-NRS group (Fig. 5B). This observation implies that individuals within the L-NRS profile population may exhibit increased sensitivity to immune checkpoint inhibitor (ICI) therapy. The H-NRS group, on the other hand, had higher tumor purity (Fig. 5D) and a lower immune infiltration. Moreover, macrophages were highly infiltrated in the H-NRS group, suggesting that a H-NRS profile is associated with immunosuppressive status. Concurrently, the Cibersort algorithm was utilized to evaluate the prevalence of various immune cell infiltration (Figure S2). The findings revealed that the H-NRS group is enriched with M2 macrophages, indicating that H-NRS characteristics imply an immunosuppressive state.

Fig. 5
figure 5

Immune infiltration characteristics of the NRS. A The abundance of each TME-infiltrating cell type was quantified by the ssGSEA algorithm. B Assessment of infiltration abundance of 27 immune checkpoints in two NRS groups. C Immunity score, ESTIMATE score, and stroma score were used to quantify the different immune statues between the two NRS groups. D Comparison of tumorpurity between the two NRS groups. E Tide score prediction sample for treatment response to immunotherapy. F Radar plots showed the immunogram patterns of the two NRS groups. GThe SubMap algorithm evaluated the expression similarity between the two NRS groups and the patients with different immunotherapy responses. NR: no response to treatment. R: responds to treatment. ns > 0.05, *P ≤ 0.05, **P < 0.01, ***P < 0.001

Immunotherapy and potential drug targets

Subsequently, the evaluation through Tide revealed that the population with L-NRS characteristics responded better to immunotherapy than those with H-NRS characteristics (Fig. 5E). To systematically assess the relationship between NRS characteristics and the potential for immunotherapy, we established a comprehensive tumor immune evaluation, an immunogen of a quantitative cancer immune cycle proposed by Karasaki et al. (Fig. 5F). The results indicated that H-NRS characteristics exhibited lower immune cycle scores, consistent with previous findings suggesting insufficient immune infiltration. Conversely, the L-NRS group displayed higher anti-tumor-related and immune-activated immune cycle scores. Therefore, individuals with L-NRS characteristics are more likely to benefit from immunotherapy. Therefore, individuals with L-NRS characteristics are more likely to experience positive outcomes from immunotherapy. Additionally, results based on the Submap algorithm indicated that patients with L-NRS characteristics showed higher expression similarity with those responsive to PD- 1 and PD- 1/PD-L1 inhibitors (Fig. 5G), suggesting an increased sensitivity to immunotherapy. Overall, the research findings indicate that patients with L-NRS characteristics may derive clinical benefits from immunotherapy, while those with H-NRS characteristics may not be suitable for it due to high costs and the potential for immune-related adverse events associated with a low response rate. Consequently, it is imperative to develop and implement additional therapeutic strategies to improve clinical outcomes for the H-NRS group. Using the oncoPredict package, we assessed the sensitivity of four chemotherapy drugs for HNSCC based on NRS characteristics (Fig. 6A). A decrease in IC50 values indicates increased sensitivity to individual drugs, and the analysis showed that the H-NRS group was more sensitive to chemotherapy drugs, including Docetaxel, Gemcitabine, Cisplatin, and Paclitaxel. Additionally, drug evaluation models based on the CTRP and PRISM databases suggested that dasatinib and gefitinib could be viable treatment options for patients with H-NRS characteristics (Fig. 6B-C).

Fig. 6
figure 6

Prediction of drug sensitivity. (A) Variations in IC50 value of cisplatin, gemcitabine, docetaxel, and paclitaxel across high- and low-risk patients. (B-C) Identification of 2 potential therapeutic agents for H-NRS group. ns > 0.05, *P ≤ 0.05, **P < 0.01, ***P < 0.001

Identification of key regulatory genes in the NRS model

To identify the key regulatory factors in the NRS model, we selected genes with relative importance > 0.3 as final markers through random survival forest analysis, showing the importance order of seven genes (Fig. 7A Figure S3A). Next, we used ROC diagnostic curves to further filter key regulatory factors, identifying genes with an ROC > 0.60, which included CCR7, LYZ, OLR1, THBS1, and TPP1 (Fig. 7B, Figure S3B). Survival analysis revealed that TPP1, OLR1, and THBS1 were risk genes associated with poor prognosis (Figure S3C). We examined the expression distribution of key regulatory genes across multiple cell subpopulations (Fig. 7C-D), revealing that CCR7 is enriched in B cells, LYZ and OLR1 are highly enriched in neutrophils, and TPP1 is enriched in dendritic cells and neutrophils. Furthermore, enrichment analysis shows that the expression level of OLR1 in key regulatory genes is notably higher in pro-cancer pathways such as KRAS signaling up, epithelial-mesenchymal transition, and angiogenesis (Fig. 7E).

Fig. 7
figure 7

Identification of key regulatory genes in the NRS model. A Ranking the importance of genes in NRS characteristics. B The AUC value for the survival analysis of key regulatory genes. C the expression of key regulatory genes in different cells. D the distribution and amount of five key regulatory genes. E Enrichment analysis of key regulatory genes

OLR1: A predictive biomarker of NRS characteristics indicates tumor invasion

As a risk gene specifically enriched in neutrophils, we chose OLR1 for further analysis. Firstly, compared to the normal controls in the TCGA-HNSCC dataset, the expression level of OLR1 was significantly higher in tumor tissues. (Fig. 8A). Compared to the normal control cell line (HOK), OLR1 mRNA levels were significantly elevated in HNSCC cell lines (CAL- 27, SCC- 4, and SCC- 25) (Fig. 8B). Survival analysis revealed a significant association between high OLR1 expression and poor prognosis (Fig. 8C), and high OLR1 expression was also associated with higher T staging and pathological grade (Fig. 8D). Moreover, patients with high OLR1 expression showed better responses to radiotherapy (Fig. 8E). We further explored the association between TME and OLR1 by using several algorithms such as CIBERSORT, CIBERSORT_ABS, Quantiseq, and MCPcounter to evaluate the correlation between immune cell infiltration levels and OLR1 expression. (Fig. 8F). The results of these algorithms all indicate that OLR1 expression shows a positive association with multiple cell types in the TCGA-HNSCC data, especially macrophages, regulatory T cells, and neutrophils, while it is negatively correlated with specific T cell subtypes, such as CD8 + T cells (Fig. 8G). At the same time, in terms of treatment sensitivity assessment, HNSCC patients with high OLR1 who received immunotherapy had better prognoses (Fig. 8H). The above analysis indicates the potential of the OLR1 gene to serve as a transcriptional marker of neutrophil behavior in the TME and a prognostic and therapeutic marker for HNSCC patients.

Fig. 8
figure 8

OLR1 correlated with HNSCC prognosis, immune infiltration status, and treatment response. A Comparison of OLR1 expression between tumor and normal samples in HNSCC. B qRT-PCR analysis of OLR1 expression in HNSCC cell lines. C Differences in prognosis between high and low OLR1 patients. D OLR1 expression in HNSCC tissues with different T stages. E Differences in radiotherapy response between high OLR1 and low OLR1 patients. F Correlation analysis between OLR1 expression and immune cell infiltration. G Correlation Analysis between OLR1 expression and M2 macrophages, Tregs, and CD8 T cells in the TCGA-HNSCC data. H Differences in progression-free survival between high OLR1 and low OLR1 patients in the immunotherapy cohort. ns > 0.05, *P ≤ 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001

The impact of OLR1 gene on neutrophil infiltration, tumor proliferation, and Metastasis in HNSCC

We established an OLR1 knockdown HNSCC cell model using siRNA. Subsequently, after inducing HL- 60 cell differentiation, we used them as a substitute for neutrophils in chemotaxis assays (Fig. 9A), and flow cytometry along with Giemsa staining confirmed their successful differentiation (Figure S4A-B). Next, we analyzed the role of OLR1 in HNSCC progression. GSEA analysis indicated that OLR1-related molecules were predominantly involved in the epithelial-mesenchymal transition (EMT) pathway (Fig. 9B), and OLR1 expression was significantly correlated with EMT (R = 0.571, p < 0.001, Figure S4C). To further investigate the relationship between OLR1 and EMT, we performed Western blot analysis to detect the expression of EMT-related proteins in HNSCC cell lines (CAL- 27 and SCC- 4). The results showed that OLR1 knockout increased the expression of E-cadherin and decreased the expression of N-cadherin in CAL- 27 and SCC- 4 cells. These findings suggest that silencing OLR1 reduces the EMT ability of HNSCC cells (Fig. 9C). The results of CCK- 8, Transwell, and wound healing assays revealed that the loss of OLR1 significantly decreased the proliferation, migration, and invasion abilities of HNSCC cells (Figs. 9D-F). We propose that OLR1 may serve as a novel predictive biomarker for HNSCC progression and metastasis related to NRS characteristics.

Fig. 9
figure 9

A The molecules significantly associated with OLR1 were enriched in epithelial-mesenchymal transition pathway. B Correlation analysis of OLR1 with epithelial-mesenchymal transition mesenchymal transition pathways score. C Expression levels of OLR1, N-cadherin, E-cadherin, and vimentin in OLR1-knockdown and control HNSCC cells. D CCK- 8 analysis in HNSCC cell lines. E Transwell migration and invasion assay for HNSCC cells. F Wound healing assay for HNSCC cells. ns P > 0.05, *P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001

Discussion

The majority of inflammatory cells present in solid tumors are neutrophils, and a high intratumoral density of these cells is significantly associated with tumor grading, lymph node metastasis, and tumor staging [41]. It is posited that both the tumors and the TME regulate neutrophil recruitment, with TANs potentially contributing to tumor pathogenesis through multiple mechanisms. TANs have been documented to generate substantial quantities of reactive nitrogen species (RNS) and reactive oxygen species (ROS) within the TME [42,43,44]. These reactive species not only induce DNA damage but also promote tumor initiation and progression. In addition to their direct carcinogenic effects, neutrophils within the TME contribute to tumor immune evasion through the secretion of inhibitory cytokines and the induction of immunosuppressive cells [45]. Additionally, TANs contribute to tumor cell invasion and metastasis by releasing matrix metalloproteinases (MMPs), which degrade the extracellular matrix [46]. A comprehensive understanding of the precise mechanisms by which neutrophils influence tumor dynamics is essential for developing novel cancer treatment strategies. Presently, several therapies targeting neutrophils are advancing to the clinical trial phase in the context of tumor treatment [47,48,49]. Unfortunately, there is a paucity of systematic research concerning TANs in HNSCC. In response, we have developed a multi-biomarker model centered on TANs, designed to aid in assessing prognosis and the TME in HNSCC patients. This model provides a theoretical foundation for implementing personalized precision treatment strategies.

The inherent limitations in transcriptome sequencing methodologies impede the comprehensive investigation of neutrophil differential expression. Nevertheless, neutrophils show considerable heterogeneity within solid tumors [41, 50]. We identified tumor-associated neutrophil-related genes (NRGs) utilizing scRNA-seq data to specifically investigate the neutrophils from tumor samples. Leveraging these NRGs, we developed an optimal neutrophil risk scoring model, which delineated two distinct HNSCC subgroups with significant prognostic disparities. Notably, elevated NRS characteristics were correlated with a poorer prognosis. Validation sets derived from both internal and cross-platform data confirm the model's robustness, with outstanding AUC values underscoring its precision. The risk score exhibits strong concordance with malignant clinical pathology, as high NRS features are associated with advanced T staging and stage classification. Moreover, multifactorial analysis results indicate that NRS features can function as independent prognostic factors. These findings suggest that NRS features can reliably predict the prognosis of HNSCC patients. Further construction of a nomogram enhances the clinical utility of NRS features, as observed through calibration curves showing a high consistency between predicted and actual values, proving the nomogram's ability to provide accurate survival rate estimates for HNSCC patients.

Enrichment analysis shows that cell cycle regulatory pathways such as DNA repair, E2 F pathway, and G2/M checkpoint are significantly activated in the H-NRS group. Regarding the cell cycle process, researchers point out that activated neutrophils cause target cells to accumulate in G2/M, consistent with the installation of DNA damage checkpoints [51, 52]. These findings suggest that TANs drive tumorigenic transformation by increasing cell replication errors through DNA damage induction [53]. Additionally, the H-NRS group shows increased activity in pathways promoting tumor progression and metastasis, such as ECM, hypoxia, angiogenesis, and epithelial-mesenchymal transition. Hypoxia within tumors exerts significant selective pressure, thereby facilitating tumor cell angiogenesis, the establishment of metastatic niches, invasion, and disseminating metastases [54]. Research findings further indicate that TANs facilitate tumor cell migration and invasion by promoting epithelial-to-mesenchymal transition, corroborating our results [55,56,57]. In the L-NRS group, many immune pathways are predominantly enriched, including the T cell receptor signaling pathway, adaptive immune activation, and B cell receptor signaling pathway. They all play crucial roles in the transmission of immune signals, immune responses, and killing of tumor cells, indicating that L-NRS patients have higher immunogenicity.

Immune therapy provides new opportunities for the treatment of various cancers, including HNSCC. In recent years, immune therapy represented by nivolumab and pembrolizumab has become the first-line treatment choice for metastatic or recurrent HNSCC patients, and these interventions have to some extent improved the prognosis of some HNSCC patients, but the overall response rate is still below 20% [58, 59]. In addition, immune therapy is relatively expensive and may pose small but real risks, leading to harmful side effects. Therefore, personalized immune therapy is particularly important, as it can significantly reduce the incidence of adverse events [60]. The immune status of patients will affect their individual response to immune therapy, ultimately determining the prognosis. Our study shows that the H-NRS group is characterized by M2 macrophage infiltration. It has been shown that the HNSCC microenvironment is more prone to inducing M2 polarization, promoting immune suppression, angiogenesis, metastasis, and invasion, leading to poor prognosis [61]. TANs induce CD4 + T cells to express the immune-suppressive transforming growth factor-bata (TGF-β) and IL- 10 and promote the recruitment of M2 macrophages through the expression of CCL2 and CCL17, generating an immune-suppressive TME, promoting local disease progression [62, 63]. In contrast, the L-NRS group exhibited a marked increase in immune scores and high infiltration of NK cells, CD4 + T cells, and CD8 + T cells. Previous studies have shown that proliferating and activated T cells release numerous cytokines and chemokines, which can enhance the host's anti-tumor defense capabilities and improve the effectiveness of immunotherapy. This has been validated in our predictions of immunotherapy, as L-NRS has shown better treatment response in multiple immunotherapy cohorts. Currently, chemotherapy remains a crucial strategy for treating HNSCC patients. Therefore, we assessed the variance in chemotherapy outcomes between the two risk groups. Compared to immunotherapy, traditional chemotherapy regimens such as cisplatin, gemcitabine, docetaxel, and cisplatin are more effective for H-NRS patients, indicating that H-NRS patients may benefit more from monotherapy or combination chemotherapy. Overall, the results indicate that NRS characteristics can serve as predictive indicators for chemotherapy sensitivity in HNSCC patients. Furthermore, this study predicted the chemotherapy drugs that may exhibit high sensitivity in patients, thus providing potential treatment options for HNSCC patients.

We finally screened five key genes (CCR7, LYZ, OLR1, THBS1 and TPP1). CCR7 is a member of the chemokine receptor family and is closely associated with lymph node metastasis in various cancers. Studies have shown that CCR7 is highly expressed in HNSCC and promotes lymph node metastasis by regulating tumor cell migration and invasion [64]. Additionally, CCR7 plays a role in recruiting immune cells within the tumor microenvironment, potentially contributing to tumor immune evasion [65, 66]. LYZ is a lysozyme associated with innate immunity and serves as a crucial component of the innate immune system [67]. Studies have shown that LYZ can act as an immune-related biomarker for predicting both the response to immunotherapy and the prognosis of various cancers [68, 69]. OLR1 has been reported to be involved in malignant progression across various cancers [70]. It also influences multiple immune cells, such as tumor-associated macrophages and T cells, facilitating tumor immune suppression [71, 72]. Our study reveals the functional role of OLR1 in HNSCC, demonstrating that its high expression significantly enhances the malignant phenotype of tumor cells and is associated with tumor-associated neutrophils. Thrombospondin- 1 (THBS1) is a matricellular protein that plays a crucial role in angiogenesis and immune suppression [73]. Studies have shown that THBS1 contributes to immune suppression and metastasis in various types of cancer [74]. Our study further supports the pivotal role of THBS1 in HNSCC and suggests its potential as a therapeutic target. TPP1 is a lysosomal enzyme that has been found to be abnormally expressed in various cancers. Research suggests that TPP1 plays a crucial role in maintaining telomeres, thereby promoting tumor growth [75]. Additionally, TPP1 has been identified as a potential biomarker for gastric cancer [76].

Prior research findings have demonstrated that OLR1 is overexpressed across various cancer types [77]. We hypothesize that the upregulation of OLR1 could be instrumental in the pathogenesis of HNSCC. Our research indicates that OLR1 is overexpressed in HNSCC tissues, particularly in later clinical pathological stages, which is also a potential indicator of poor prognosis for HNSCC patients. TME plays a crucial role in both tumor immune evasion and immune clearance, significantly influencing the progression of HNSCC. Studies have shown that OLR1 expression is closely associated with the infiltration levels of various immune cells, particularly tumor-associated macrophages (TAMs) and regulatory T cells (Tregs) [70], both of which contribute to tumor immune suppression, and metastasis. OLR1 signaling may activate key pathways such as NF-κB and PI3 K/Akt [78,79,80], which are central to the induction of pro-inflammatory and immunosuppressive gene programs, thereby promoting tumor progression. The interaction between oxidized low-density lipoprotein (oxLDL) and OLR1 can trigger these pathways. In HNSCC, OLR1 expression is strongly correlated with tumor-associated macrophage infiltration, with high OLR1 levels predicting poorer prognosis. Neutrophil populations with high OLR1 expression exhibit increased levels of nitric oxide (NO), reactive oxygen species (ROS), and the immunosuppressive factor ARG1, which inhibit T cell proliferation and enhance tumor immune suppression [81]. Analysis of the single-cell RNA sequencing dataset for HNSCC reveals a specific enrichment of OLR1 expression in neutrophils. This observation suggests that the OLR1 gene may serve as a transcriptional marker indicative of neutrophil activity within the TME of HNSCC. In addition, our research revealed that OLR1 is involved in the regulation the malignant phenotype of HNSCC, and its expression is involved in epithelial-mesenchymal transition, tumor angiogenesis, and hypoxia signaling pathways. In vitro experiments have shown that silencing OLR1 expression significantly inhibits the malignant biological behavior of HNSCC, consistent with previous studies. Furthermore, the knockout of OLR1 significantly reduces tumor cell proliferation, migration, and invasion. In conclusion, we have established the tumor-associated neutrophil characteristics of HNSCC, which hold significant clinical value in prognosis assessment and immune therapy response prediction. This study expands our understanding of the role of TANs in HNSCC and identifies a potential therapeutic target for HNSCC. However, our research has some limitations. Firstly, the study mainly relied on publicly available datasets, and although we extensively validated the constructed scores, it is undeniable that predictive ability needs further validation in large-scale prospective clinical studies. Furthermore, this study only confirmed the mRNA and protein expression of characteristic genes, and conducted cell function experiments of OLR1 in HNSCC cells. Therefore, further experimental exploration is needed to elucidate the molecular mechanisms of how neutrophil characteristic genes regulate the prognosis of HNSCC.

Data availability

Data used in this work can be acquired from the TCGA (https://portal.gdc.cancer.gov/) and GEO (https://www.ncbi.nlm.nih.gov/geo/). Gene expression data from GDSC (https://www.cancerrxgene.org/), CTRP (http://portals.broadinstitute.org/ctrp/) and PRISM (http://cosbi.ku.edu.tr/prism) were also obtained. Other data supporting the findings of this study are available from the corresponding author upon reasonable request.

Abbreviations

HNSCC:

Head and neck squamous cell carcinoma

TANs:

Tumor-associated neutrophils

TCGA:

The Cancer Genome Atlas

GEO:

Gene Expression Omnibus

CCK- 8:

Cell counting kit‐8

NRS:

Neutrophils-associated signature

TME:

The tumor microenvironment

DCs:

Dendritic cells

NK:

Natural killer

ROS:

Reactive oxygen species

Inos:

Inducible nitric oxide synthase

ARG1:

Arginase 1

VEGF:

Vascular endothelial growth factor

CTL:

Cytotoxic T lymphocyte

scRNA-seq:

Single-cell RNA sequencing

PCA:

Principal component analysis

GSEA:

Gene set enrichment analysis

MSigDB:

Molecular Signatures Database

GSVA:

Gene set variation analysis

ssGSEA:

Single-sample gene set enrichment analysis

CIC:

Cancer Immune Cycle

TIDE:

Tumor Immune Dysfunction and Exclusion

SubMap:

Subclass Mapping

GDSC:

Genomics of Drug Sensitivity in Cancer

IC50:

Half-maximal inhibition concentrations

PRISM:

Profiling Relative Inhibition Simultaneously in Mixtures

CTRP:

Cancer Therapeutics Response Portal

AUC:

Area under the dose–response curve

FBS:

Fetal bovine serum

cDNA:

Complementary DNA

OD:

Optical density

PBS:

Phosphate-buffered saline

References

  1. Siegel RL, Miller KD, Jemal A. Cancer statistics, 2019. CA Cancer J Clin. 2019;69(1):7–34.

    Article  PubMed  Google Scholar 

  2. Zhou Q, Yuan O, Cui H, Hu T, Xiao GG, Wei J, Zhang H, Wu C. Bioinformatic analysis identifies HPV-related tumor microenvironment remodeling prognostic biomarkers in head and neck squamous cell carcinoma. Front Cell Infect Microbiol. 2022;12:1007950.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  3. Johnson DE, Burtness B, Leemans CR, Lui VWY, Bauman JE, Grandis JR. Head and neck squamous cell carcinoma. Nat Rev Dis Primers. 2020;6(1):92.

    Article  PubMed  PubMed Central  Google Scholar 

  4. Liu Y, Zhang N, Wen Y, Wen J. Head and neck cancerpathogenesis and targeted therapy. MedComm (2020). 2024;5(9):e702.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  5. Burtness B, Harrington KJ, Greil R, Soulières D, Tahara M, de Castro G, Jr., Psyrri A, Basté N, Neupane P, Bratland Å, et al. Pembrolizumab alone or with chemotherapy versus cetuximab with chemotherapy for recurrent or metastatic squamous cell carcinoma of the head and neck (KEYNOTE-048): a randomised, open-label, phase 3 study. Lancet. 2019;394(10212):1915–28.

    Article  CAS  PubMed  Google Scholar 

  6. Machiels JP, Tao Y, Licitra L, Burtness B, Tahara M, Rischin D, Alves G, Lima IPF, Hughes BGM, Pointreau Y, et al. Pembrolizumab plus concurrent chemoradiotherapy versus placebo plus concurrent chemoradiotherapy in patients with locally advanced squamous cell carcinoma of the head and neck (KEYNOTE-412): a randomised, double-blind, phase 3 trial. Lancet Oncol. 2024;25(5):572–87.

    Article  CAS  PubMed  Google Scholar 

  7. Leemans CR, Snijders PJF, Brakenhoff RH. The molecular landscape of head and neck cancer. Nat Rev Cancer. 2018;18(5):269–82.

    Article  CAS  PubMed  Google Scholar 

  8. Jin MZ, Jin WL. The updated landscape of tumor microenvironment and drug repurposing. Signal Transduct Target Ther. 2020;5(1):166.

    Article  PubMed  PubMed Central  Google Scholar 

  9. Fu LQ, Du WL, Cai MH, Yao JY, Zhao YY, Mou XZ. The roles of tumor-associated macrophages in tumor angiogenesis and metastasis. Cell Immunol. 2020;353:104119.

    Article  CAS  PubMed  Google Scholar 

  10. Wculek SK, Cueto FJ, Mujal AM, Melero I, Krummel MF, Sancho D. Dendritic cells in cancer immunology and immunotherapy. Nat Rev Immunol. 2020;20(1):7–24.

    Article  CAS  PubMed  Google Scholar 

  11. Jaillon S, Ponzetta A, Di Mitri D, Santoni A, Bonecchi R, Mantovani A. Neutrophil diversity and plasticity in tumour progression and therapy. Nat Rev Cancer. 2020;20(9):485–503.

    Article  CAS  PubMed  Google Scholar 

  12. 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–67.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  13. Zhang Y, Guoqiang L, Sun M, Lu X. Targeting and exploitation of tumor-associated neutrophils to enhance immunotherapy and drug delivery for cancer treatment. Cancer Biol Med. 2020;17(1):32–43.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  14. Christoffersson G, Vågesjö E, Vandooren J, Lidén M, Massena S, Reinert RB, Brissova M, Powers AC, Opdenakker G, Phillipson M. VEGF-A recruits a proangiogenic MMP-9-delivering neutrophil subset that induces angiogenesis in transplanted hypoxic tissue. Blood. 2012;120(23):4653–62.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  15. Ikeda M, Morimoto M, Tajimi M, Inoue K, Benhadji KA, Lahn MMF, Sakai D. A phase 1b study of transforming growth factor-beta receptor I inhibitor galunisertib in combination with sorafenib in Japanese patients with unresectable hepatocellular carcinoma. Invest New Drugs. 2019;37(1):118–26.

    Article  CAS  PubMed  Google Scholar 

  16. Cheng AL, Kang YK, He AR, Lim HY, Ryoo BY, Hung CH, Sheen IS, Izumi N, Austin T, Wang Q, et al. Safety and efficacy of tigatuzumab plus sorafenib as first-line therapy in subjects with advanced hepatocellular carcinoma: A phase 2 randomized study. J Hepatol. 2015;63(4):896–904.

    Article  CAS  PubMed  Google Scholar 

  17. Bitzer M, Horger M, Giannini EG, Ganten TM, Wörns MA, Siveke JT, Dollinger MM, Gerken G, Scheulen ME, Wege H, et al. Resminostat plus sorafenib as second-line therapy of advanced hepatocellular carcinoma - The SHELTER study. J Hepatol. 2016;65(2):280–8.

    Article  CAS  PubMed  Google Scholar 

  18. Li J, Yu T, Sun J, Ma M, Zheng Z, Kang W, Ye X. Comprehensive integration of single-cell RNA and transcriptome RNA sequencing to establish a pyroptosis-related signature for improving prognostic prediction of gastric cancer. Comput Struct Biotechnol J. 2024;23:990–1004.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  19. Kulkarni A, Anderson AG, Merullo DP, Konopka G. Beyond bulk: a review of single cell transcriptomics methodologies and applications. Curr Opin Biotechnol. 2019;58:129–36.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  20. Choi JH, Lee BS, Jang JY, Lee YS, Kim HJ, Roh J, Shin YS, Woo HG, Kim CH. Single-cell transcriptome profiling of the stepwise progression of head and neck cancer. Nat Commun. 2023;14(1):1055.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  21. Wichmann G, Rosolowski M, Krohn K, Kreuz M, Boehm A, Reiche A, Scharrer U, Halama D, Bertolini J, Bauer U, et al. The role of HPV RNA transcription, immune response-related gene expression and disruptive TP53 mutations in diagnostic and prognostic profiling of head and neck cancer. Int J Cancer. 2015;137(12):2846–57.

    Article  CAS  PubMed  Google Scholar 

  22. Lohavanichbutr P, Méndez E, Holsinger FC, Rue TC, Zhang Y, Houck J, Upton MP, Futran N, Schwartz SM, Wang P, et al. A 13-gene signature prognostic of HPV-negative OSCC: discovery and external validation. Clin Cancer Res. 2013;19(5):1197–203.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  23. Gribov A, Sill M, Lück S, Rücker F, Döhner K, Bullinger L, Benner A, Unwin A. SEURAT: visual analytics for the integrated analysis of microarray data. BMC Med Genomics. 2010;3:21.

    Article  PubMed  PubMed Central  Google Scholar 

  24. Aran D, Looney AP, Liu L, Wu E, Fong V, Hsu A, Chak S, Naikawadi RP, Wolters PJ, Abate AR, et al. Reference-based analysis of lung single-cell sequencing reveals a transitional profibrotic macrophage. Nat Immunol. 2019;20(2):163–72.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  25. Cohen M, Giladi A, Gorki AD, Solodkin DG, Zada M, Hladik A, Miklosi A, Salame TM, Halpern KB, David E, et al. Lung Single-Cell Signaling Interaction Map Reveals Basophil Role in Macrophage Imprinting. Cell. 2018;175(4):1031-1044.e1018.

    Article  CAS  PubMed  Google Scholar 

  26. Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, Paulovich A, Pomeroy SL, Golub TR, Lander ES, et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U S A. 2005;102(43):15545–50.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  27. Hänzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics. 2013;14:7.

    Article  PubMed  PubMed Central  Google Scholar 

  28. Bindea G, Mlecnik B, Tosolini M, Kirilovsky A, Waldner M, Obenauf AC, Angell H, Fredriksen T, Lafontaine L, Berger A, et al. Spatiotemporal dynamics of intratumoral immune cells reveal the immune landscape in human cancer. Immunity. 2013;39(4):782–95.

    Article  CAS  PubMed  Google Scholar 

  29. Ward-Kavanagh LK, Lin WW, Šedý JR, Ware CF. The TNF Receptor Superfamily in Co-stimulating and Co-inhibitory Responses. Immunity. 2016;44(5):1005–19.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  30. Janakiram M, Chinai JM, Zhao A, Sparano JA, Zang X. HHLA2 and TMIGD2: new immunotherapeutic targets of the B7 and CD28 families. Oncoimmunology. 2015;4(8):e1026534.

    Article  PubMed  PubMed Central  Google Scholar 

  31. Chrétien S, Zerdes I, Bergh J, Matikas A, Foukakis T. Beyond PD-1/PD-L1 Inhibition: What the Future Holds for Breast Cancer Immunotherapy. Cancers (Basel). 2019;11(5):628.

    Article  PubMed  Google Scholar 

  32. Karasaki T, Nagayama K, Kuwano H, Nitadori JI, Sato M, Anraku M, Hosoi A, Matsushita H, Morishita Y, Kashiwabara K, et al. An Immunogram for the Cancer-Immunity Cycle: Towards Personalized Immunotherapy of Lung Cancer. J Thorac Oncol. 2017;12(5):791–803.

    Article  PubMed  Google Scholar 

  33. Hoshida Y, Brunet JP, Tamayo P, Golub TR, Mesirov JP. Subclass mapping: identifying common subtypes in independent disease data sets. PLoS ONE. 2007;2(11):e1195.

    Article  PubMed  PubMed Central  Google Scholar 

  34. Yang W, Soares J, Greninger P, Edelman EJ, Lightfoot H, Forbes S, Bindal N, Beare D, Smith JA, Thompson IR, et al. Genomics of Drug Sensitivity in Cancer (GDSC): a resource for therapeutic biomarker discovery in cancer cells. Nucleic Acids Res. 2013;41(Database issue):D955-961.

    CAS  PubMed  Google Scholar 

  35. Geeleher P, Cox N, Huang RS. pRRophetic: an R package for prediction of clinical chemotherapeutic response from tumor gene expression levels. PLoS ONE. 2014;9(9):e107468.

    Article  PubMed  PubMed Central  Google Scholar 

  36. Seashore-Ludlow B, Rees MG, Cheah JH, Cokol M, Price EV, Coletti ME, Jones V, Bodycombe NE, Soule CK, Gould J, et al. Harnessing Connectivity in a Large-Scale Small-Molecule Sensitivity Dataset. Cancer Discov. 2015;5(11):1210–23.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  37. Finotello F, Mayer C, Plattner C, Laschober G, Rieder D, Hackl H, Krogsdam A, Loncova Z, Posch W, Wilflingseder D, et al. Molecular and pharmacological modulators of the tumor immune contexture revealed by deconvolution of RNA-seq data. Genome Med. 2019;11(1):34.

    Article  PubMed  PubMed Central  Google Scholar 

  38. Becht E, Giraldo NA, Lacroix L, Buttard B, Elarouci N, Petitprez F, Selves J, Laurent-Puig P, Sautès-Fridman C, Fridman WH, et al. Estimating the population abundance of tissue-infiltrating immune and stromal cell populations using gene expression. Genome Biol. 2016;17(1):218.

    Article  PubMed  PubMed Central  Google Scholar 

  39. Elmusrati A, Wang J, Wang CY. Tumor microenvironment and immune evasion in head and neck squamous cell carcinoma. Int J Oral Sci. 2021;13(1):24.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  40. Zheng Z, Xu Y, Shi Y, Shao C. Neutrophils in the tumor microenvironment and their functional modulation by mesenchymal stromal cells. Cell Immunol. 2022;379:104576.

    Article  CAS  PubMed  Google Scholar 

  41. Hedrick CC, Malanchi I. Neutrophils in cancer: heterogeneous and multifaceted. Nat Rev Immunol. 2022;22(3):173–87.

    Article  CAS  PubMed  Google Scholar 

  42. Willson JA, Arienti S, Sadiku P, Reyes L, Coelho P, Morrison T, Rinaldi G, Dockrell DH, Whyte MKB, Walmsley SR. Neutrophil HIF-1α stabilization is augmented by mitochondrial ROS produced via the glycerol 3-phosphate shuttle. Blood. 2022;139(2):281–6.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  43. Korbecki J, Simińska D, Gąssowska-Dobrowolska M, Listos J, Gutowska I, Chlubek D, Baranowska-Bosiacka I. Chronic and Cycling Hypoxia: Drivers of Cancer Chronic Inflammation through HIF-1 and NF-κB Activation: A Review of the Molecular Mechanisms. Int J Mol Sci. 2021;22(19):10701.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  44. van Vlerken-Ysla L, Tyurina YY, Kagan VE, Gabrilovich DI. Functional states of myeloid cells in cancer. Cancer Cell. 2023;41(3):490–504.

    Article  PubMed  PubMed Central  Google Scholar 

  45. Giese MA, Hind LE, Huttenlocher A. Neutrophil plasticity in the tumor microenvironment. Blood. 2019;133(20):2159–67.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  46. Xiong S, Dong L, Cheng L. Neutrophils in cancer carcinogenesis and metastasis. J Hematol Oncol. 2021;14(1):173.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  47. Huang X, Nepovimova E, Adam V, Sivak L, Heger Z, Valko M, Wu Q, Kuca K. Neutrophils in Cancer immunotherapy: friends or foes? Mol Cancer. 2024;23(1):107.

    Article  PubMed  PubMed Central  Google Scholar 

  48. Barry ST, Gabrilovich DI, Sansom OJ, Campbell AD, Morton JP. Therapeutic targeting of tumour myeloid cells. Nat Rev Cancer. 2023;23(4):216–37.

    Article  CAS  PubMed  Google Scholar 

  49. Hirschhorn D, Budhu S, Kraehenbuehl L, Gigoux M, Schröder D, Chow A, Ricca JM, Gasmi B, De Henau O, Mangarin LMB, et al. T cell immunotherapies engage neutrophils to eliminate tumor antigen escape variants. Cell. 2023;186(7):1432-1447.e1417.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  50. Xue R, Zhang Q, Cao Q, Kong R, Xiang X, Liu H, Feng M, Wang F, Cheng J, Li Z, et al. Liver tumour immune microenvironment subtypes and neutrophil heterogeneity. Nature. 2022;612(7938):141–7.

    Article  CAS  PubMed  Google Scholar 

  51. Mahat U, Garg B, Yang CY, Mehta H, Hanna R, Rogers HJ, Flagg A, Ivanov AI, Corey SJ. Lymphocyte cytosolic protein 1 (L-plastin) I232F mutation impairs granulocytic proliferation and causes neutropenia. Blood Adv. 2022;6(8):2581–94.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  52. Szczerba BM, Castro-Giner F, Vetter M, Krol I, Gkountela S, Landin J, Scheidmann MC, Donato C, Scherrer R, Singer J, et al. Neutrophils escort circulating tumour cells to enable cell cycle progression. Nature. 2019;566(7745):553–7.

    Article  CAS  PubMed  Google Scholar 

  53. Poli V, Zanoni I. Neutrophil intrinsic and extrinsic regulation of NETosis in health and disease. Trends Microbiol. 2023;31(3):280–93.

    Article  CAS  PubMed  Google Scholar 

  54. Ng MSF, Kwok I, Tan L, Shi C, Cerezo-Wallis D, Tan Y, Leong K, Calvo GF, Yang K, Zhang Y, et al. Deterministic reprogramming of neutrophils within tumors. Science. 2024;383(6679):eadf6493.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  55. Li S, Cong X, Gao H, Lan X, Li Z, Wang W, Song S, Wang Y, Li C, Zhang H, et al. Tumor-associated neutrophils induce EMT by IL-17a to promote migration and invasion in gastric cancer cells. J Exp Clin Cancer Res. 2019;38(1):6.

    Article  PubMed  PubMed Central  Google Scholar 

  56. Wang Y, Li X, Zhang T, Li F, Shen Y, He Y, You Q, Zhang Y, Zhai J, Yao X, et al. Neutrophils promote tumor invasion via FAM3C-mediated epithelial-to-mesenchymal transition in gastric cancer. Int J Biol Sci. 2023;19(5):1352–68.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  57. Ou B, Liu Y, Gao Z, Xu J, Yan Y, Li Y, Zhang J. Senescent neutrophils-derived exosomal piRNA-17560 promotes chemoresistance and EMT of breast cancer via FTO-mediated m6A demethylation. Cell Death Dis. 2022;13(10):905.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  58. Ruffin AT, Li H, Vujanovic L, Zandberg DP, Ferris RL, Bruno TC. Improving head and neck cancer therapies by immunomodulation of the tumour microenvironment. Nat Rev Cancer. 2023;23(3):173–88.

    Article  CAS  PubMed  Google Scholar 

  59. Bhatia A, Burtness B. Treating Head and Neck Cancer in the Age of Immunotherapy: A 2023 Update. Drugs. 2023;83(3):217–48.

    Article  CAS  PubMed  Google Scholar 

  60. Guo N, Ni K, Luo T, Lan G, Arina A, Xu Z, Mao J, Weichselbaum RR, Spiotto M, Lin W. Reprogramming of Neutrophils as Non-canonical Antigen Presenting Cells by Radiotherapy-Radiodynamic Therapy to Facilitate Immune-Mediated Tumor Regression. ACS Nano. 2021;15(11):17515–27.

    Article  CAS  PubMed  Google Scholar 

  61. Chen Y, Li ZY, Zhou GQ, Sun Y. An Immune-Related Gene Prognostic Index for Head and Neck Squamous Cell Carcinoma. Clin Cancer Res. 2021;27(1):330–41.

    Article  PubMed  Google Scholar 

  62. Liu ZL, Meng XY, Bao RJ, Shen MY, Sun JJ, Chen WD, Liu F, He Y. Single cell deciphering of progression trajectories of the tumor ecosystem in head and neck cancer. Nat Commun. 2024;15(1):2595.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  63. Zhang S, Sun L, Zuo J, Feng D. Tumor associated neutrophils governs tumor progression through an IL-10/STAT3/PD-L1 feedback signaling loop in lung cancer. Transl Oncol. 2024;40:101866.

    Article  CAS  PubMed  Google Scholar 

  64. Zhang Z, Liu F, Li Z, Wang D, Li R, Sun C. Jak3 is involved in CCR7-dependent migration and invasion in metastatic squamous cell carcinoma of the head and neck. Oncol Lett. 2017;13(5):3191–7.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  65. Wang Z, Kirkwood KL, Wang Y, Du W, Lin S, Zhou W, Yan C, Gao J, Li Z, Sun C, et al. Analysis of the effect of CCR7 on the microenvironment of mouse oral squamous cell carcinoma by single-cell RNA sequencing technology. J Exp Clin Cancer Res. 2024;43(1):94.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  66. Zhao ZJ, Liu FY, Li P, Ding X, Zong ZH, Sun CF. CCL19-induced chemokine receptor 7 activates the phosphoinositide-3 kinase-mediated invasive pathway through Cdc42 in metastatic squamous cell carcinoma of the head and neck. Oncol Rep. 2011;25(3):729–37.

    CAS  PubMed  Google Scholar 

  67. Ragland SA, Criss AK. From bacterial killing to immune modulation: Recent insights into the functions of lysozyme. PLoS Pathog. 2017;13(9):e1006512.

    Article  PubMed  PubMed Central  Google Scholar 

  68. Taube JM, Young GD, McMiller TL, Chen S, Salas JT, Pritchard TS, Xu H, Meeker AK, Fan J, Cheadle C, et al. Differential Expression of Immune-Regulatory Genes Associated with PD-L1 Display in Melanoma: Implications for PD-1 Pathway Blockade. Clin Cancer Res. 2015;21(17):3969–76.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  69. Zou D, Wang Y, Wang M, Zhao B, Hu F, Li Y, Zhang B. Bioinformatics analysis reveals the competing endogenous RNA (ceRNA) coexpression network in the tumor microenvironment and prognostic biomarkers in soft tissue sarcomas. Bioengineered. 2021;12(1):496–506.

    Article  CAS  PubMed  Google Scholar 

  70. Zhang P, Zhao Y, Xia X, Mei S, Huang Y, Zhu Y, Yu S, Chen X. Expression of OLR1 gene on tumor-associated macrophages of head and neck squamous cell carcinoma, and its correlation with clinical outcome. Oncoimmunology. 2023;12(1):2203073.

    Article  PubMed  PubMed Central  Google Scholar 

  71. Yang G, Xiong G, Feng M, Zhao F, Qiu J, Liu Y, Cao Z, Wang H, Yang J, You L, et al. OLR1 Promotes Pancreatic Cancer Metastasis via Increased c-Myc Expression and Transcription of HMGA2. Mol Cancer Res. 2020;18(5):685–97.

    Article  CAS  PubMed  Google Scholar 

  72. Sweetwyne MT, Murphy-Ullrich JE. Thrombospondin1 in tissue repair and fibrosis: TGF-β-dependent and independent mechanisms. Matrix Biol. 2012;31(3):178–86.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  73. Omatsu M, Nakanishi Y, Iwane K, Aoyama N, Duran A, Muta Y, Martinez-Ordoñez A, Han Q, Agatsuma N, Mizukoshi K, et al. THBS1-producing tumor-infiltrating monocyte-like cells contribute to immunosuppression and metastasis in colorectal cancer. Nat Commun. 2023;14(1):5534.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  74. You Y, Tian Z, Du Z, Wu K, Xu G, Dai M, Wang Y, Xiao M. M1-like tumor-associated macrophages cascade a mesenchymal/stem-like phenotype of oral squamous cell carcinoma via the IL6/Stat3/THBS1 feedback loop. J Exp Clin Cancer Res. 2022;41(1):10.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  75. Chun-On P, Hinchie AM, Beale HC, Gil Silva AA, Rush E, Sander C, Connelly CJ, Seynnaeve BKN, Kirkwood JM, Vaske OM, et al. TPP1 promoter mutations cooperate with TERT promoter mutations to lengthen telomeres in melanoma. Science. 2022;378(6620):664–8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  76. Wang Q, Liu Y, Li Z, Tang Y, Long W, Xin H, Huang X, Zhou S, Wang L, Liang B, et al. Establishment of a novel lysosomal signature for the diagnosis of gastric cancer with in-vitro and in-situ validation. Front Immunol. 2023;14:1182277.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  77. Wu L, Liu Y, Deng W, Wu T, Bu L, Chen L. OLR1 Is a Pan-Cancer Prognostic and Immunotherapeutic Predictor Associated with EMT and Cuproptosis in HNSCC. Int J Mol Sci. 2023;24(16):12904.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  78. Wang B, Zhao H, Zhao L, Zhang Y, Wan Q, Shen Y, Bu X, Wan M, Shen C. Up-regulation of OLR1 expression by TBC1D3 through activation of TNFα/NF-κB pathway promotes the migration of human breast cancer cells. Cancer Lett. 2017;408:60–70.

    Article  PubMed  Google Scholar 

  79. Liang M, Zhang P, Fu J. Up-regulation of LOX-1 expression by TNF-alpha promotes trans-endothelial migration of MDA-MB-231 breast cancer cells. Cancer Lett. 2007;258(1):31–7.

    Article  CAS  PubMed  Google Scholar 

  80. Li C, Zhang J, Wu H, Li L, Yang C, Song S, Peng P, Shao M, Zhang M, Zhao J, et al. Lectin-like oxidized low-density lipoprotein receptor-1 facilitates metastasis of gastric cancer through driving epithelial-mesenchymal transition and PI3K/Akt/GSK3β activation. Sci Rep. 2017;7:45275.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  81. Luo H, Ikenaga N, Nakata K, Higashijima N, Zhong P, Kubo A, Wu C, Tsutsumi C, Shimada Y, Hayashi M, et al. Tumor-associated neutrophils upregulate Nectin2 expression, creating the immunosuppressive microenvironment in pancreatic ductal adenocarcinoma. J Exp Clin Cancer Res. 2024;43(1):258.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

Download references

Funding

This study received no funding.

Author information

Authors and Affiliations

Authors

Contributions

HYC and ZKL provided direction and guidance throughout the preparation of this manuscript. HYC and ZKL wrote and edited the manuscript. HYC and ZKL reviewed and made significant revisions to the manuscript. HYC, ZKL, YCY, JLK, SBG, XCS, LKHF, and YKL collected and prepared the related papers. All authors read and approved the final manuscript.

Corresponding authors

Correspondence to Daoke Yang or Yingjuan Zheng.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Cui, H., Li, Z., Liu, Y. et al. Integrated analysis of single-cell RNA-seq and bulk RNA-seq unravels the molecular feature of tumor-associated neutrophils of head and neck squamous cell carcinoma. BMC Cancer 25, 821 (2025). https://doiorg.publicaciones.saludcastillayleon.es/10.1186/s12885-025-14179-9

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doiorg.publicaciones.saludcastillayleon.es/10.1186/s12885-025-14179-9

Keywords