The study was approved by the Stomatology Hospital of Air Force Medical University, and all patients participated in this study had signed the informed consent. Twenty patients with HNSCC who underwent surgery between January 2023 and January 2024 provided primary tumor tissues. The diagnosis of HNSCC was based on WHO Classification of Head and Neck Tumors (5th edition) and the TNM staging system (8th edition, UICC). The collected tissues were fixed in 10% neutral-buffered formalin and embedded in paraffin for subsequent pathological examination and staining.

Gene expression datasets were obtained from the Gene Expression Omnibus (GEO) repository under accession numbers GSE65858 (bulk RNA-seq), GSE150430 (single-cell RNA-seq), and GSE123813 (single-cell RNA-seq). Fifteen primary tumor samples from GSE150430 were included in this study ( Table 1 ). RNA-seq FPKM expression profiles, overall survival (OS) data, and clinical annotations for HNSCC were retrieved from the National Cancer Institute (NCI)’s Genomic Data Commons (GDC).

Fifteen primary samples from GSE150430 were processed using the R software package and underwent quality control in GEO. Two thousand highly variable genes were identified using the “FindVariableFeatures” algorithm. Principal component analysis (PCA) was conducted on these HVGs, with the top 50 principal components retained for subsequent analyses. TAMs associated marker genes were identified using the “FindAllMarkers” function(P<0.05). WGCNA constructed co-expression modules (minimum size=30 genes) through soft thresholding. Finally, module eigengenes were intersected with TAMs signature genes identified through single-cell analysis to pinpoint M2 TAMs related genes.

The relative abundance of M1 and M2 macrophages was quantified using the XCell algorithm. Samples were stratified into high-risk and low-risk groups based on macrophage infiltration levels, applying a median cut-off value derived from the R package XCell. Kaplan-Meier survival analysis with log-rank testing was performed to evaluate the correlation between macrophage infiltration density and overall survival (OS) in HNSCC patients.

Univariate Cox regression was performed to identify M2 TAMs-related prognostic genes based on the survival curve (P<0.05). The R package was employed to construct a LASSO Cox regression model to identify prognostic factors. A risk score model was developed by weighting key prognostic factors with LASSO regression coefficients to predict survival.

Based on corresponding scores, fifteen samples were classified into high-risk group and low risk group and survival curve were visualized using the Kaplan-Meier method with the log-rank test. The receiver operating characteristic (ROC) curve was adapted to evaluate the predictive performance of the scoring system, and the area under the curve (AUC) was visualized with the R package time ROC. Univariate and multivariate Cox regression analyses were performed to evaluate the independent prognostic value of the risk score.

The half-maximal inhibitory concentration (IC50) values for training set samples were estimated using the Phenotype algorithm implemented in the R package Predict (v1.2.3), with drug sensitivity data sourced from the Genomics of Drug Sensitivity in Cancer (GDSC) database (version 2.0; PMID: 22000000).

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were performed using the R package clusterProfiler (v4.0.1) to functionally annotate the signature genes. According to the enrichment analysis result, differences in immune function between the high-risk group and low-risk group were compared using Gene Set Variation Analysis (GSVA) and Gene Set Enrichment Analysis (GSEA). Immune cell infiltration across subgroups was compared using the Wilcoxon test and the ssGSEA (single-sample gene-set enrichment analysis) algorithm. The relative abundance of 28 immune cell subsets (e.g., activated CD8+ T cells, dendritic cells, and macrophages) within the tumor microenvironment (TME) were quantified.

Tissue sections were fixed with 4% paraformaldehyde in PBS for 20 minutes at room temperature. Next, the membranes were blocked with 1% BSA for 2 hours at room temperature. The membranes were then incubated with primary antibodies against CD163 and iNOS at 4°C overnight. The results were observed using a laser scanning confocal microscope.

Sections and TMAs were stained with or incubated with primary antibodies using the Elivision™ Plus Polymer HRP immunohistochemistry kit (Maxim, Fujian, China). The following antibodies were used: anti-FCGBP (ab121199, Abcam, 1:500), anti-GPR35 (ab150635, Abcam, 1:300), GIMAP7 polyclonal antibody (Proteintech, 1:500), WIPF1 polyclonal antibody (Proteintech, 1:500), RASGEF1B polyclonal antibody (Proteintech, 1:300), p47 phox polyclonal antibody (Proteintech, 1:400), HEXB polyclonal antibody (Proteintech, 1:400), IL-10 monoclonal antibody (Proteintech, 1:500), CLECL1 monoclonal antibody (Proteintech, 1:500),GIMAP5 polyclonal antibody (AtaGenix,1:800),and CLECL1 polyclonal antibody(AtaGenix,1:800). The score of each section was classified into 0–4 by the ImageJ software based on the intensity and the positive rate of stained cells.

All experimental protocols were approved by the Institutional Animal Care and Use Committee (IACUC) of the State Key Laboratory, Air Force Military Medical University. Female BALB/c nude mice (6-week-old, n=20) were subcutaneously inoculated with SCC9 (1–5 × 10⁷cells/mouse) into the left forelimb. One-week post-inoculation, mice were randomized into four groups: the experimental group first (n=5) received lipopolysaccharide (LPS) (20 mg/kg, 200 μL), group second (n=5) was administered recombinant IL-4 (20 mg/kg, 200 μL), while the control group third (blank) (n=5) and forth(n=5) with (SCC+Normal Saline,200 μL). LPS and IL-4 treatments were administered to the mice one week after tumor cell injection, when the subcutaneous tumorigenesis model was successfully established with a tumor volume (V) > 100 mm³ (to avoid the impact of early intervention on tumor formation). The experiment was terminated when the tumor volume (V) < 1500 mm³, and all mice were euthanized on day 27 in accordance with the ethical norms for animal experiments. Following confirmation of tumor formation, injections were administered three times at 48-hour intervals. Tumor tissues were harvested for immunohistochemical (IHC) analysis of CD68 (pan-macrophage marker), CD163 (M2 macrophage marker), and iNOS (M1 macrophage marker) expression. Prognostic significance was evaluated based on tumor weight as a key measure of tumor progression in HNSCC.

Statistical analyses were conducted using R software (version 4.1.2). The Wilcoxon rank-sum test was applied to compare differences between two groups, while the Kruskal-Wallis’s test was used for comparisons involving multiple groups. Survival curves for prognostic analysis were generated using the Kaplan-Meier method, and the log-rank test was used to determine the significance of differences. In graphical representations, significance levels were denoted as follows: ns (not significant, P > 0.05), *(P < 0.05), ** (P < 0.01), and *** (P < 0.001).

TAMs infiltration was more abundant in stage III/IV tumors in stage I/II based on the expression of CD68 in immunofluorescence ( Figures 1A, B ). Survival curve based on TCGA data confirmed that M2 macrophages are strongly associated with poor prognosis in malignancies, while M1-type macrophages associated with better prognosis ( Figures 1C, D ). Mouse subcutaneous tumor model was constructed to reveal the infiltration situation ( Figure 1E ). The weight and volume of SCC9 with IL-4 which means M2 TAMs rich infiltration group were higher than SCC9 with LPS which means M1-type TAMs rich infiltration, control and blank (P<0.05) ( Figures 1F, G ). Immunofluorescence staining also showed that M2 TAMs infiltration was richer based on the expression of CD163 in SCC9 with IL-4 group than in SCC9 with LPS (P<0.05) ( Figures 1H, I ).

WGCNA was employed to screen M2 macrophages-related genes to explore their association with the prognosis of HNSCC. As shown in Supplementary File 1 , the WGCNA results identified M2 macrophage-related genes in HNSCC and revealed 25 optimal modules. Among these 25 optimal modules, the blue module, which exhibited the highest Pearson correlation coefficient, was selected for downstream analysis and contained 778 genes.

From 46,001 single-cell transcriptomes, the top 2,000 highly variable genes (HVGs), including CRNN, CRCT1, and HLA-DRA, were selected for further analysis ( Figure 2A ). Uniform Manifold Approximation and Projection (UMAP) visualization was used to display the top 50 principal components (PCs) and reveal distinct cellular clusters. The Harmony algorithm was applied to correct batch effects ( Figure 2B ). Cell-type-specific marker genes were identified using the Find All Markers function and the top five markers per cluster were visualized ( Figure 2C ). In total, 1,208 TAMs-specific signature genes were identified.

The intersection of the 1208 TAMs-specific signature genes from single-cell database and the 778 M2 macrophage-associated genes from bulk database resulted in 259 candidate M2 TAMs-associated genes ( Figure 3A ). Functional enrichment analysis revealed that key prognostic factors were significantly enriched in immune processes such as T cell proliferation and lymphocyte proliferation ( Figure 3B ). Twenty-nine genes associated with clinical prognosis were identified through univariate cox analysis and were furthered narrowed down to 11 key prognostic genes by LASSO Cox regression analysis ( Figures 3C–E ). The HNSCC samples were classified into two subgroups (cluster1 and cluster2), based on the expression of the 11 key factors ( Figures 3F, G , Supplementary File 2 ). Prognostic survival rates in cluster 2 showed a significantly worse prognosis ( Figures 3H ). The subsequent analysis involved examining the expression levels of these 11 key prognostic genes to establish a prognosis-related signature.

The gene coefficients from the linear combination of 11 key prognostic factors were used to define the prognostic signature for each patient, as shown in Table 2 .

According to the median value, cases were classified into high-risk group and low-risk group. Kaplan-Meier survival analysis and log-rank tests revealed that patients in the high-risk group had a significantly poorer prognosis in GSE65858 (P<0.02, Figure 4A ). Based on multivariate and univariate Cox regression analysis performed with clinical features and prognosis signature, demonstrated a consistent trend in predicting prognosis. These analyses also confirmed that the prognostic signature was an independent prognostic factor (HR=1.65, P-value=0.01, Supplementary Figures S2B ). The data from GSE65858 further supported the prognosis signature to be an independent prognostic factor ( Figures 2C D ). Patients in the high-risk group had significantly worse overall survival rate in TCGA (P<0.001, Figure 4B ). A nomogram was plotted, incorporating clinical factors such as stage, sex, lymph vascular invasion (LVI), and perineural invasion (PNI) to provide a more comprehensive survival prediction ( Figure 4C ). The calibration curve ( Figures 4F–H ) and decision curve analysis ( Figure 4E ) demonstrated its reliability of the model. The concordance index (c-index) analysis demonstrated that the prognostic signature exhibited higher accuracy than other clinicopathological indicators ( Figure 4D ). Next, the distribution of prognosis signature among clinical pathology characteristics was analyzed.

The distribution of the prognostic signature across clinical pathological characteristics, showed that the proportion of patients with terminal cancer was higher in the high-risk group ( Figure 5A ). However, there was no significant difference between the two groups in terms of gender ( Figure 5B ), or lymphovascular invasion (LVI) ( Figure 5C ). The proportion of patients with perineural invasion (PNI) in the high-risk group was significantly higher in those with advanced stage disease ( Figure 5D ).

The risk scores of patients with PNI were elevated in those with advanced disease ( Figure 5E ). However, no significant difference in risk scores was observed between genders ( Figure 5F ). Furthermore, there was no significant difference in risk scores between patients with or without LVI ( Figure 5G ). Risk scores were significantly higher in patients with PNI ( Figure 5H ). This malignant biological behavior, which appeared frequently in the high-risk group, indicated a poor prognosis.

In different clinical pathological characteristic subgroups, the signatures showed that the prognosis in patients with advanced stage ( Figure 5M ), men ( Figure 5N ), without PNI ( Figure 5K ), without LVI ( Figure 5K ), female ( Figure 5L ), and those with PNI ( Figure 5P ) were poorer. Patients in the high-risk group generally had worse prognosis. On the contrary, patients in the early stage ( Figure 5I ), in women ( Figure 5J ), showed no significant difference in prognosis between high- and low-risk groups. ( Figure 5O ). These findings suggest that the prognostic signature can effectively predict the poor prognosis of HNSCC.

GSEA analysis revealed that drug metabolism pathways were significantly activated in the high-risk group ( Figure 6A ). Conversely, immune-related biological processes, including the activation and proliferation of B cells and T cells, were markedly enhanced in the low-risk group ( Figures 6B, C ). The immune cell scores were compared before and after treatment. Th17 cells showed a significantly increase in cell scores following immunotherapy in patients who responded to immunotherapy, which can be interpreted as an increase in cellular activity. ( Figures 6D, E ). Activated B cells, immature B cells, and natural killer T cells were significantly more abundant in the low-risk group, whereas the proportion of CD56+ natural killer cells was lower than that in the high-risk group ( Figure 6G ). Subsequently, comparison of immune checkpoint-related gene expression revealed significantly higher expression of CD276 in the high-risk group ( Figure 6F ). The distribution of biomarkers in the high-risk group, as determined by the signature we constructed, can serve as a better indicator for assessing poor prognosis.

Based on the expression levels of the 11 key prognostic factors and gene coefficients in a linear combination, the prognostic signature of twenty patients was evaluated, as shown in Table 3 .

Tissue sections from these patients with HNSCC were examined immunohistochemically, and their expression scores were analyzed ( Figure 7A ). The obtained scores were incorporated into constructed signatures to generate predictions ( Table 4 ). The gene expression levels were consistent with the key genes predicted by LASSO regression for the prognostic factors ( Figures 7B–D ).

According to the constructed prognostic signature, we calculated the risk score for patients in the training set and stratified them into high-risk and low-risk group based on the median risk score ( Figure 7E ). The proportion of patients under 60 years old was higher in the high-risk group and lower in the low-risk group ( Figure 7F ). Furthermore, the high-risk group also had a significantly higher number of stage III/IV patients compared to stage I/II patients ( Figure 7G ). Gender did not show a significantly difference ( Figure 7H ). Additionally, low-risk group had significantly more patients without invasion compared to those with invasion ( Figure 7I ). These findings suggest that the prognosis signature has the potential to predict the prognosis effectively in patients with HNSCC.