A pan-cancer longitudinal atlas uncovers a Macro_BIRC5-Treg axis driving immune checkpoint blockade resistance in head and neck squamous cell carcinoma

Introduction

Immune checkpoint blockade (ICB) has fundamentally transformed the treatment landscape of multiple malignant solid tumors over the past decade. Despite these advances, durable clinical responses are achieved in only a subset of patients, highlighting the persistent challenge of primary and acquired resistance to immunotherapy (1). The magnitude of this challenge varies substantially across cancer types: while melanoma patients achieve objective response rates of 33–45% with anti- programmed cell death protein 1 (PD-1) therapy, head and neck squamous cell carcinoma (HNSCC) demonstrates markedly lower response rates of only 13–16% in platinum-refractory disease and approximately 23% in routine clinical practice (2-5). Increasing evidence indicates that such heterogeneous therapeutic outcomes are not solely determined by tumor-intrinsic factors but are critically shaped by the tumor microenvironment (TME) (6).

While cytotoxic T lymphocytes have long been considered the principal mediators of ICB efficacy, recent studies have underscored the pivotal role of tumor-infiltrating myeloid cells (TIMs)—one of the most abundant immune populations in the TME—in orchestrating anti-tumor immunity and regulating treatment response (7,8). TIMs encompass diverse cell types, including macrophages, monocytes, and dendritic cells (DCs), which exhibit remarkable plasticity and functional diversity (9), and actively shape the balance between immune activation and suppression within the TME.

In immunotherapy-resistant cancers, suppressive myeloid populations can impair anti-tumor immunity through several interconnected mechanisms. First, tumor-associated macrophages (TAMs) and myeloid-derived suppressive cells can attenuate T-cell activation through inhibitory ligand-receptor pathways, including programmed death-ligand 1 (PD-L1)/PD-1 and major histocompatibility complex (MHC)-I-LILRB signaling. In particular, engagement of MHC-I by LILRB1 has been shown to suppress macrophage-mediated phagocytosis, indicating that innate immune checkpoint pathways may also contribute to tumor immune escape (10). Second, suppressive myeloid cells can remodel the cytokine and chemokine milieu of the TME. For example, CCL22 produced by tumor cells and TAMs promotes the recruitment of regulatory T cells (Tregs), thereby fostering a locally immune-privileged microenvironment (11). Third, myeloid cells can impose metabolic constraints on effector T cells through arginine depletion, nitric oxide production, reactive oxygen species generation, and nutrient competition. Consistently, arginase I produced by mature myeloid cells in the TME has been shown to inhibit T-cell receptor expression and antigen-specific T-cell responses (12). These results indicate that myeloid-mediated immune suppression is not driven by a single pathway, but by coordinated inhibitory, chemokine-dependent, and metabolic programs.

Myeloid-mediated immunosuppression is particularly prominent in HNSCC, which is characterized by extensive infiltration of immunosuppressive myeloid populations, including TAMs, myeloid-derived suppressor cells (MDSCs), and tumor-associated neutrophils (TANs) (13-16). This profoundly immunosuppressive TME not only limits baseline T cell function but also creates a barrier to effective ICB therapy. Understanding how these myeloid populations are dynamically remodeled during ICB treatment in HNSCC, and how this compares to cancers with higher response rates such as melanoma, is critical for developing effective combination strategies.

Accordingly, the classical M1/M2 polarization model, which categorizes macrophages into pro-inflammatory or immunosuppressive states, is increasingly recognized as insufficient to capture the complex and context-dependent phenotypes observed in vivo (17,18). Single-cell RNA sequencing (scRNA-seq) technologies have substantially expanded our understanding of myeloid heterogeneity, revealing a continuum of transcriptional and functional states shaped by tumor types, tissue context, and therapeutic exposure (19-21). Beyond their intrinsic functions, TIMs serve as key regulators of adaptive immunity through dynamic interactions with T cells. Specific myeloid subsets can suppress T cell activity via inhibitory ligands, metabolic competition, and the recruitment or activation of Treg, whereas other subsets—particularly specialized DCs—are essential for antigen presentation and the initiation of effective anti-tumor T cell responses during ICB (22-24). These results highlight the importance of defining functionally distinct TIM subsets and understanding how their states and interactions evolve during therapy.

Despite these advances, critical gaps remain. Most existing studies rely on cross-sectional analyses of single tumor types, limiting the ability to distinguish shared, therapy-associated myeloid programs from cancer-specific features (20). More importantly, our understanding of how TIMs are dynamically remodeled during ICB treatment—and how those changes relate to clinical response—remains incomplete. Longitudinal single-cell analyses tracking the same patients before and after treatment are still scarce, particularly across multiple cancer types (25). As a result, it remains unclear which myeloid states truly emerge in response to therapeutic pressure and which intercellular communication programs drive resistance to ICB.

To address these limitations, we assembled a longitudinal pan-cancer single-cell atlas comprising 82 patients across basal cell carcinoma (BCC) (26), clear cell renal cell carcinoma (ccRCC) (27,28), HNSCC (29,30), and skin cutaneous melanoma (SKCM) (31), with pre- and post-treatment tumor samples annotated with precise clinical response outcomes. To capture therapy-induced myeloid remodeling, we applied scCURE, a computational framework that distinguishes myeloid cells undergoing significant transcriptional shifts during treatment from stable background populations (32).

Using this approach, we systematically delineated divergent myeloid remodeling programs related to ICB response and resistance. We found that effective therapy is characterized by coordinated metabolic and functional reprogramming of myeloid cells, whereas non-responders (NRs) exhibit persistent immunosuppressive and metabolically constrained states. Notably, we identified distinct resistance-associated macrophage populations, including Macro_C1QC and a hyper-oxidative Macro_BIRC5 subsets, which promote immunosuppressive networks through extensive crosstalk with Treg. These results have particular relevance for HNSCC, where the Macro_BIRC5-Treg axis may represent a key therapeutic vulnerability. Together, our study provides a longitudinal and pan-cancer view of therapy-driven myeloid remodeling and identifies specific myeloid-Treg interaction axes that may be exploited to overcome resistance to ICB.

Methods

Patient cohorts and data collection

This study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments. To construct a longitudinal pan-cancer single-cell atlas suitable for analyzing therapy-associated myeloid remodeling, we systematically screened publicly available scRNA-seq datasets of ICB-treated tumors, including Yost et al. (26) (GEO: GSE123814), Bi et al. (27) (Single Cell Portal: SCP1288), Au et al. (28) (EGA: EGAD00001008166), van der Leun et al. (29) (GEO: GSE232240), Luoma et al. (30) (GEO: GSE200996), and Sade-Feldman et al. (31) (GEO: GSE120575). Eligible cohorts were required to meet the following criteria: (I) availability of longitudinal tumor biopsies collected before treatment (Pre) and during or after ICB treatment (Post); (II) availability of single-cell transcriptomic data with sufficient sequencing depth and cell quality for downstream myeloid subset analysis; and (III) clear patient-level clinical response annotations, such as complete response/partial response (CR/PR) versus stable disease/progressive disease (SD/PD) according to Response Evaluation Criteria in Solid Tumors (RECIST) or study-specific response criteria, enabling the classification of patients into responders (REs) and NR. Cohorts with insufficient myeloid cell recovery, inadequate sequencing depth, or incomplete response annotation were excluded.

Based on these criteria, we retrospectively collected scRNA-seq data covering four cancer types: HNSCC, ccRCC, SKCM, and BCC. The final aggregated cohort encompassed 129 samples from 82 patients, enabling multi-dimensional stratification according to treatment timepoint and clinical response.

The four included cancer types are all clinically relevant to ICB treatment but represent distinct immunological and tissue contexts. HNSCC is directly relevant to oral and maxillofacial oncology and shows heterogeneous responses to PD-1/PD-L1-based immunotherapy (33). ccRCC is also an important ICB-treated malignancy and is characterized by prominent immune infiltration, vascular remodeling, and a myeloid-rich TME (34). SKCM is one of the most established solid tumor models for checkpoint blockade, particularly in the context of PD-1- and CTLA-4-targeted immunotherapy (35). BCC provides an additional epithelial skin tumor context in which PD-1 blockade, such as cemiplimab, has demonstrated clinical activity in advanced or difficult-to-treat disease after hedgehog pathway inhibitor therapy (36). Therefore, these cohorts were not intended to exhaustively represent all malignancies, but rather to provide complementary ICB-treated tumor contexts with distinct tissue origins, response patterns, and myeloid landscapes.

In this study, “Pre” samples were defined as tumor biopsies collected before ICB treatment, whereas “Post” samples referred to tumor biopsies collected during or after ICB treatment (at least two doses) according to the original study annotations. Because the included cohorts were derived from different clinical studies, the exact post-treatment sampling time varied across datasets. Detailed sample-level metadata was recorded in Table S1.

External bulk transcriptomic datasets were used for survival analysis and clinical correlation validation. The Cancer Genome Atlas (TCGA) cohorts, including TCGA-HNSCC, TCGA-ccRCC, TCGA-SKCM, TCGA-lung adenocarcinoma (LUAD), and TCGA-lung squamous cell carcinoma (LUSC), were used to assess tumor-associated expression patterns, clinical stage association, and TME scores. In addition, previously published ICB-treated bulk cohorts were used for immunotherapy-related survival validation, including the FoyJP2022 HNSCC cohort (37), Braun2020 ccRCC cohort (38-41), OAK2017 and POPLAR2017 non-small cell lung cancer (NSCLC) cohorts (42), IMvigor210 urothelial carcinoma cohort (43), and Liu 2019 and Riaz 2017 SKCM cohorts (44,45). Detailed information on cohort name, cancer type, treatment regimen, sample size, data source is provided in Table S2.

scRNA-seq data preprocessing and integration

Publicly available scRNA-seq expression matrices and matched metadata were obtained from the original studies. When available, processed gene-cell count matrices generated by the original authors were used directly for downstream analysis. Raw read alignment and initial gene-cell matrix generation were therefore not repeated in this study and were performed as described in the corresponding original publications.

To ensure comparability across cohorts, standardized downstream preprocessing was performed using the Seurat R package (version 4.4.0). Cells were filtered according to the following quality control criteria: detected gene count <200 or >6,000, and mitochondrial gene expression proportion >15%. After quality control, data were normalized using the LogNormalize method, and highly variable genes (HVGs) were identified for downstream dimensionality reduction.

Principal component analysis (PCA) was then performed based on the selected HVGs. To minimize technical variation across patients, samples, and independent datasets, batch correction was performed using Harmony (version 0.1.1), with cohort- and sample-level identifiers used to account for inter-sample and inter-dataset variation. The Harmony-corrected embeddings were subsequently used for Uniform Manifold Approximation and Projection (UMAP) visualization, nearest-neighbor graph construction, and graph-based clustering. Cell clusters were manually annotated as major immune and stromal cell lineages according to canonical marker gene expression.

scCURE analysis and myeloid subset identification

To precisely delineate myeloid cell subsets undergoing significant transcriptomic alterations during ICB therapy, we applied the scCURE algorithm framework. The core principles and procedures were as follows:

Data preparation and dimensionality reduction

Myeloid lineage cells were extracted from all samples, and their preprocessed expression matrices were merged. PCA was performed on the merged dataset, and the top five principal components were selected as input for subsequent modeling.

Gaussian mixture model (GMM) construction

We fitted GMMs separately for pre-treatment and post-treatment cell subsets using the expectation-maximization algorithm. The model posits that each cell subset consists of K Gaussian components, each representing a latent cellular substate or functional state. The optimal K value was determined via leave-one-out cross-validation on training data, guided by predictive performance.

Discrimination of changed versus steady cells

The Kullback-Leibler (KL) divergence was systematically calculated between all pairs of pre- and post-treatment Gaussian components. Cells belonging to component pairs that were mutual nearest neighbors based on KL divergence were classified as “Steady cells”, while the remaining cells were categorized as “Changed cells”.

Myeloid subset re-clustering

The identified “Changed” myeloid cells from all samples were integrated and subjected to unsupervised clustering using standard Seurat workflows (resolution =0.8).

Cluster quality assessment

We utilized the Ratio of Global Unshifted Entropy (ROGUE) R package, an entropy-based universal metric, to evaluate the cluster purity of each defined cell subset. This entropy-based metric quantifies cluster purity, with values approaching 1 indicating minimal intra-cluster heterogeneity. To ensure that the defined cellular states represented distinct and stable phenotypes rather than technical artifacts or transient mixtures, we assessed ROGUE scores across multiple stratifications, including individual patients, treatment timepoints, and clinical response groups.

Quantification of cellular distribution preference [ratio of observed to expected (Ro/e) analysis]

To rigorously quantify the distributional bias and preferential enrichment of myeloid subsets across different treatment timepoints, we calculated the Ro/e cell numbers. This analysis was performed separately for RE and NR cohorts to deconvolute response-specific dynamic shifts.

We utilized the calTissueDist function from the sscVis R package with the calculation method based on the chi-square test. The expected cell numbers for each subset were estimated based on the null hypothesis of a random distribution across timepoints. The Ro/e index was then defined as the ratio of the observed cell count to this expected value. An Ro/e ratio greater than 1 signifies a preferential enrichment of a specific cell type at a given timepoint, whereas a ratio less than 1 indicates depletion.

Differential expression and functional enrichment analysis

Differentially expressed genes (DEGs) between cell subsets or treatment conditions were identified using the FindAllMarkers functions in Seurat. The significance threshold was set at an adjusted P value <0.05 and |log₂ (fold change)| >0.5. Gene Ontology (GO) biological process and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were performed on identified DEGs or specific gene modules using the clusterProfiler R package.

Trajectory analysis

Single-cell trajectory analysis was conducted separately on RE and NR myeloid cells identified by scCURE using the Monocle2 R package (version 2.3.1). Cellular developmental trajectories were constructed based on HVGs, with the state primarily composed of pre-treatment cells defined as the pseudotime root. Monocle differential expression testing function was employed to identify genes changing significantly along the pseudotime trajectory, which were then organized into distinct dynamic expression modules via hierarchical clustering.

Metabolic pathway scoring

Metabolic pathway scoring was performed using the AddModuleScore function in the Seurat R package. Curated metabolic gene sets were obtained from the MSigDB and KEGG pathway databases, including pathways related to oxidative phosphorylation (OXPHOS), glycolysis, fatty acid metabolism, cholesterol metabolism, and amino acid metabolism. The normalized expression matrix was used as input, and a module score was calculated for each pathway in each single cell. This score represents the relative activity of a given metabolic pathway based on the average expression of its member genes after subtracting the expression of control gene sets with comparable expression levels.

The metabolic scores were calculated for each myeloid cell and compared across myeloid subsets, cancer types, treatment timepoints, and clinical response groups.

Cell-cell communication analysis

Cell-cell communication analysis was performed using the LIANA R package and the CellChat R package. The normalized expression matrices and cell type annotations were used as input. The analysis was based on known ligand-receptor interactions between sender and receiver cell populations.

LIANA was first used to identify ligand-receptor interactions between myeloid subsets and between myeloid cells and T cells. Filtering was performed based on aggregate rank and cell-specific interaction frequency. Aggregate rank reflects the overall confidence of each ligand-receptor interaction, with lower values indicating higher-ranked interactions. Cell-specific interaction frequency reflects the frequency of ligand-receptor interactions between specific sender and receiver cell populations. Ligand-receptor interactions with an aggregate rank <0.01 were retained, and interactions with cell-specific interaction frequencies below the 75th percentile of the overall frequencies within each sample condition were discarded.

Separate communication networks were constructed for post-treatment REs and NRs. Interaction frequencies, interaction strengths, and source-target specificities were compared between the two groups to identify resistance-associated communication programs. CellChat was further used to calculate communication probabilities and visualize global signaling networks. Network-level metrics, including interaction number, interaction strength, and information flow, were used to identify dominant communication hubs and signaling pathways associated with ICB resistance.

Transcription factor (TF) activity analysis

TF activity analysis was performed using the DoRothEA regulon resource and the VIPER R package. The goal of this analysis was to identify TFs associated with myeloid cell states across different cancer types, treatment timepoints, and clinical response groups. DoRothEA provides curated TF-target regulons with assigned confidence levels, and VIPER estimates TF activity scores based on the coordinated expression patterns of TF target genes.

The normalized expression matrix was used as input. Regulons with confidence levels A, B, and C were retained for downstream analysis. VIPER was then applied to calculate TF activity scores for each cell or each myeloid subset. To identify TFs associated with therapy-related myeloid remodeling, TF activity scores were compared across myeloid subsets, Pre/Post treatment groups, and RE/NR groups. TFs showing subset-specific or response-associated activity patterns were considered candidate regulators of myeloid functional states.

Gene signature construction and survival analysis

For key myeloid subsets, specific gene signatures were derived from the top five most discriminant marker genes. In external bulk RNA-seq clinical cohorts, the abundance of subsets was estimated by calculating the average expression of these signature genes. Patients were stratified into high-expression and low-expression groups based on the optimal cutoff of the signature score. Survival curves were generated using the Kaplan-Meier method, and differences in overall survival (OS) between groups were assessed using the log-rank test. Immune and stromal scores for tumor samples were calculated based on bulk transcriptomic data using the ESTIMATE algorithm.

Spatial transcriptomics (ST) analysis

Publicly available ST datasets for HNSCC and SKCM were utilized for validation. Following standard preprocessing of ST data, single-cell defined myeloid subsets were mapped to spatial locations using deconvolution methods (SPOTlight). Spatial co-localization was assessed by calculating the correlation of different cell type scores across spatial spots.

Statistical analysis

All statistical analyses were performed using R (version 4.3.2). Comparisons between two groups were conducted using the Wilcoxon rank-sum test or t-test, while comparisons among multiple groups used the Kruskal-Wallis test or one-way analysis of variance (ANOVA). Correlation analyses were performed using Pearson or Spearman correlation coefficients. A P value <0.05 was considered statistically significant.

All data analyzed in this study were obtained from publicly available databases [Gene Expression Omnibus (GEO), TCGA, and Database of Genotypes and Phenotypes (dbGaP)]. No new human or animal experiments were conducted, and therefore institutional ethical approval and informed consent were not required.

Results

Construction of a pan-cancer myeloid cell atlas and dissection of functional heterogeneity

To systematically evaluate the heterogeneity of TIMs across cancers, we established a unified single-cell analytical framework, as illustrated in Figure 1A. This framework was applied to single-cell transcriptomic datasets from four cancer types—BCC, ccRCC, HNSCC, and SKCM—encompassing a total of 82 patients. All samples possessed complete information regarding treatment response status (RE/NR) and treatment timepoints (Pre-treatment, Pre/Post-treatment, Post). After strict quality control and filtration to remove low-quality cells, we applied batch correction to integrate the datasets. We manually annotated eight common major cell subsets based on canonical markers, primarily identifying CD4⁺ T cells, CD8⁺ T cells, DCs, mast cells, monocyte/macrophage cells, natural killer (NK) cells, plasma cells, and Treg (Figure S1A-S1E) (46).

Figure 1 Construction of a pan-cancer myeloid cell atlas and dissection of functional heterogeneity. (A) Schematic workflow of the study design. Overview of the analytical framework employed to delineate the dynamic landscape of TIMs under ICB. (B) UMAP embeddings of the integrated TIMs, colored by clinical response, treatment timepoint, cancer type, and annotated cell subtype. (C) Stacked bar plots displaying the relative frequencies of myeloid subsets across four cancer types, stratified by therapeutic response (NR vs. RE, top) and treatment timing (Pre vs. Post, bottom). (D) Heatmap showing the Ro/e ratios of each myeloid subcluster in pre- and post-treatment samples from NRs. (E) Dot plot displaying the enriched biological pathways for each myeloid cluster. BCC, basal cell carcinoma; ccRCC, clear cell renal cell carcinoma; cDC, conventional dendritic cell; COVID-19, coronavirus disease 2019; HNSCC, head and neck squamous cell carcinoma; ICB, immune checkpoint blockade; NR, non-responder; pDC, plasmacytoid dendritic cell; PPAR, peroxisome proliferator-activated receptor; RE, responder; Ro/e, ratio of observed to expected; SKCM, skin cutaneous melanoma; TIM, tumor-infiltrating myeloid cell; UMAP, Uniform Manifold Approximation and Projection.

Because the intrinsic heterogeneity of the TME often obscures treatment-induced transcriptional programs, we applied scCURE to distinguish myeloid cells undergoing therapy-related remodeling from relatively stable background populations. Briefly, scCURE uses statistical modeling to compare pre- and post-treatment cellular states and classifies cells as “Changed” or “Steady” according to the degree of transcriptional shift during treatment. After integrating and re-clustering these “Changed” myeloid cells, we identified a total of 13 distinct myeloid subsets, comprising seven macrophage subsets (e.g., Macro_SPP1, Macro_C1QC, Macro_FOLR2), three monocyte subsets (Mono_DUSP6, Mono_FOLR3, Mono_TREM1), and three DC subsets (cDC_CD1C, cDC_CLEC9A, pDC) (Figure S1E). UMAP visualization displayed the distributional patterns of these subsets across different cancer types, treatment stages, and clinical response groups (Figure 1B). The establishment of this atlas laid a solid foundation for subsequent in-depth analyses of the functional characteristics, developmental trajectories, intercellular communication, and clinical relevance of these “Changed” myeloid subsets influenced by ICB treatment.

To assess the biological robustness of the subsets, we used ROGUE, an entropy-based metric where a score near one denotes high transcriptional homogeneity (Figure S1F). Notably, the results showed that most subsets maintained high scores across different treatment timepoints and response states, indicating that despite the intrinsic plasticity, the myeloid subsets we defined represent relatively stable and distinct functional phenotypes.

After establishing the pan-cancer myeloid cell atlas, we further explored the distributional characteristics of these subsets within different clinical contexts. We compared the distribution of myeloid cell subsets between each pair of stratified groups to assess dynamic variations by cell proportion analysis (Figure 1C). In the Pre group, Macro_SPP1 was the predominant cell type present in Melanoma. In the Post group, Mono_TREM1 emerged as the most abundant population in melanoma, whereas cDC_CD1C maintained its dominance in BCC. Notably, Mono_FOLR3 constituted a major proportion in HNSCC within the Pre group but showed a marked reduction in the Post group. Macro_VSIG4 consistently constituted the major cell type in ccRCC across both Pre and Post groups, suggesting a stable lineage dependency in this cancer type. Comparing the response groups, Macro_C1QC was highly enriched in the NR group of ccRCC, whereas its proportion was substantially lower in the RE group. Conversely, Macro_HAVCR2 was specifically enriched in the RE group, while significantly reduced in the NR group. Mono_TREM1 and Macro_VSIG4 displayed increased proportions in the RE group across multiple cancer types, suggesting a potential association with favorable therapeutic outcomes. The relative abundance of these subsets exhibited dynamic changes with treatment timepoints and clinical response groups, suggesting that specific myeloid populations may be closely related to therapeutic intervention and clinical outcomes.

We next quantified the preferential distribution of myeloid subsets across different clinical response groups and treatment timepoints based on the Ro/e analysis—a metric utilized to evaluate the specific enrichment of cell subsets within distinct groups (Figure 1D, Figure S1G). The results showed that, in the NR group, Macro_C1QC and Macro_BIRC5 exhibited elevated Ro/e ratios in post-treatment samples (Ro/e =1.7), suggesting a preferential enrichment of these subsets following ineffective therapy. Conversely, Macro_HAVCR2 and Mono_DUSP6 displayed a marked trend of depletion. Strikingly, a divergent pattern was observed in the RE group, where Macro_VSIG4 was specifically accumulated, whereas Macro_C1QC and Macro_MERTK showed universally decreased enrichment scores. Collectively, these results highlight that ICB treatment profoundly remodels the myeloid cell compositional structure within the TME, particularly promoting the accumulation of specific subsets.

Finally, we characterized the functionality of these 13 identified cell types (Figure 1E). For instance, Macro_SPP1 was primarily engaged in “Viral protein interaction with cytokine and cytokine receptor” and “Mineral absorption” pathways, indicating a cytokine-active, metabolically specialized macrophage phenotype related to tissue remodeling and immune microenvironment modulation, whereas Macro_C1QC and Macro_FOLR2 were predominantly enriched in “Staphylococcus aureus infection” and “Leishmaniasis” pathways, suggesting a role in pathogen response and immune defense. Additionally, Macro_VSIG4 was uniquely linked to “Phagosome”, indicating a robust capacity for debris clearance and immune modulation. In contrast, conventional DC (cDC) subsets (cDC_CD1C, cDC_CLEC9A) showed strong enrichment for “Antigen processing and presentation”, consistent with their canonical role in adaptive immunity priming.

This analysis collectively highlights notable variations in the composition of TIMs among different cancer types, emphasizing the heterogeneity of these cells across various cancers, treatment statuses, and response groups.

ICB therapy drives divergent differentiation trajectories and functional remodeling of myeloid cells in REs versus NRs

To elucidate the dynamic remodeling of TIMs and uncover the mechanisms underlying therapeutic response, we employed Monocle to construct pseudotime trajectories, analyzing the evolutionary patterns between REs and NRs. Trajectory analysis revealed that in both R and NR patients, myeloid cells evolved from a pre-treatment state towards a post-treatment state (Figure 2A,2B). To elucidate the molecular basis driving these distinct trajectories, we tracked the dynamic expression of key functional gene sets. In REs, we observed a coordinated metabolic and immune reprogramming. As pseudotime progressed from pre- to post-treatment, key glycolytic enzymes (PKM, PGK1, ENO1) were sharply upregulated, indicating a switch in cellular metabolism towards aerobic glycolysis. The upregulation of mitochondrial electron transport chain genes, including MT-ND1, MT-ND2, MT-ND4, MT-CO1, and MT-CO2, indicates enhanced OXPHOS and mitochondrial respiratory activity, reflecting a high-energy, functionally active state of myeloid cells that supports antigen processing, inflammatory responses, and effector functions. Concurrently, the expression of key immunoinhibitors, HAVCR2 (TIM-3), and the metabolic inhibitory enzyme, IDO1, significantly decreased, marking the alleviation of immunosuppressive brakes (Figure 2C). In contrast, in the NR group, the induction of other key glycolytic genes (such as HK2, ENO1, SLC2A1) was noticeably weaker, indicating incomplete metabolic reprogramming. More critically, HAVCR2 and IDO1 sustained high expression levels throughout the entire trajectory in NRs, even rising in the later stages of treatment, suggesting failure to overcome the immunosuppressive environment. Together with this persistent immunosuppression, ICB treatment was related to a coordinated downregulation of mitochondrial electron transport chain genes, including MT-ND1, MT-ND2, MT-ND4, MT-CO1, and MT-CO2, indicating impaired OXPHOS. This metabolic collapse suggests that myeloid cells in NRs fail to sustain the energetic demands required for effective antigen presentation, inflammatory signaling, and immune effector functions under ICB pressure (Figure 2D).

Figure 2 ICB therapy drives divergent differentiation trajectories and functional remodeling of myeloid cells in REs versus NRs. (A,B) Monocle trajectories of myeloid cells from the RE (A) and NR (B) groups. Cells are colored by pseudotime (top left), treatment timepoint (top right; Pre vs. Post), cell subtype (bottom left), and trajectory state (bottom right). (C,D) Kinetic plots illustrating the expression changes of selected metabolic and immune-related genes along the pseudotime trajectory in REs (C) and NRs (D). (E) Heatmaps displaying the expression dynamics of trajectory-dependent genes along the pseudotime axis for the RE group. Genes are clustered into four distinct functional modules. The color gradient represents scaled gene expression values (Z-scores). (F,G) Dot plots illustrating the top enriched GO biological processes (F) and KEGG pathways (G) for genes within Cluster 1 of the RE group trajectory. Dot size corresponds to the gene count, and the color gradient indicates statistical significance. (H) Heatmaps displaying the expression dynamics of trajectory-dependent genes along the pseudotime axis for the NR group. Genes are clustered into four distinct functional modules. The color gradient represents scaled gene expression values (Z-scores). (I,J) Dot plots illustrating the top enriched GO biological processes (I) and KEGG pathways (J) for genes within Cluster 1 of the NR group trajectory. Dot size corresponds to the gene count, and the color gradient indicates statistical significance. Ag, antigen; cDC, conventional dendritic cell; ECM, extracellular matrix; GO, Gene Ontology; ICB, immune checkpoint blockade; IgA, immunoglobulin A; KEGG, Kyoto Encyclopedia of Genes and Genomes; MHC, major histocompatibility complex; NR, non-responder; pDC, plasmacytoid dendritic cell; RE, responder; TNF, tumor necrosis factor.

To dissect the transcriptional programs driving therapeutic outcomes, we performed unsupervised clustering analysis on trajectory-dependent genes independently for the RE and NR groups, identifying four distinct functional modules (Figure 2E). To explore the biological functions related to this dynamic expression pattern, we conducted GO and KEGG pathway functional enrichment analyses. Specifically, the characteristics of genes contained in Cluster 1 of the RE group showed high expression at the beginning of the pseudotime trajectory, followed by a gradual and significant downregulation towards the trajectory terminus. GO biological process analysis showed a striking enrichment of functional terms critical for initiating adaptive immune responses. The most significantly enriched biological processes were mainly centered on antigen presentation-related functions, including “antigen processing and presentation of peptide or polysaccharide antigen via MHC-II”, “MHC-II protein complex assembly”, and “positive regulation of leukocyte cell-cell adhesion” (Figure 2F). Corroborating these results, KEGG pathway analysis indicated that Cluster 1 genes were significantly enriched in “Antigen processing and presentation” and “Phagosome” pathways (Figure 2G). These data indicate that the gene signature defined by Cluster 1 represents a highly functional state of myeloid cells possessing efficient antigen uptake and processing capabilities, as well as the ability to perform MHC-II-restricted antigen presentation to CD4⁺ T cells. During ICB therapy, this is a prerequisite for initiating effective anti-tumor immune responses.

We further examined the remaining gene modules (Clusters 2, 3, and 4) that typically showed upregulation along the pseudotime trajectory in the RE group. These modules represent functional states where myeloid cells were successfully induced and enhanced in the context of an effective response to ICB therapy. In GO biological process analysis of Cluster 2, the significant enrichment of “response to molecule of bacterial origin”, “response to lipopolysaccharide”, as well as “positive regulation of cytokine production” and “chemotaxis” (Figure S2A) indicates a transition of myeloid cells from a quiescent state to a highly alert state, acquiring the ability to recruit other immune cells. KEGG pathway analysis further confirmed this, with core inflammatory signal hubs such as the “NF-κB signaling pathway”, “TNF signaling pathway”, and “Toll-like receptor signaling pathway” all showing high enrichment (Figure S2B). These results strongly indicate that successful ICB therapy breaks immune tolerance in the TME in REs, reactivating the potent intrinsic proinflammatory mechanisms of myeloid cells. Cluster 3 mainly reflects the execution and regulation functions of innate immunity. GO and KEGG analyses emphasized the key roles of myeloid cells in clearance and defense, including enrichment of pathways such as “phagocytosis”, “lysosome”, “complement and coagulation cascades”, and “cell killing” (Figure S2C,2D). In addition, the appearance of multiple viral infection-related pathways suggests a defense state similar to a strong interferon response. The functional enrichment results of the Cluster 4 module were almost entirely dominated by terms related to protein synthesis and energy metabolism. GO analysis showed extreme activity in “ribosome biogenesis”, “cytoplasmic translation”, and “rRNA processing” (Figure S2E), while KEGG analysis highlighted “OXPHOS” and “proteasome” pathways (Figure S2F). This indicates that myeloid cells responding to therapy are in a high-synthesis, high-energy-consumption state, providing necessary materials and an energy basis for generating large amounts of cytokines, chemokines, and performing phagocytic functions by significantly enhancing ribosome function and mitochondria-dependent OXPHOS.

In the NR group, trajectory analysis delineated a distinct transcriptional program indicative of therapeutic failure. We observed a progressive and coordinated downregulation of genes indispensable for acute inflammation and adaptive immunity priming (Clusters 1 & 2) (Figure 2H). Specifically, Cluster 1—enriched in tumor necrosis factor (TNF) and NF-κB signaling as well as leukocyte chemotaxis—showed high initial expression that was rapidly dampened, suggesting an inability to sustain the recruitment of effector cells (Figure 2I,2J). Concurrently, Cluster 2 exhibited a marked decline in antigen processing and presentation and adaptive immune response pathways. The loss of these critical functions, coupled with the downregulation of response to virus signatures, implies that myeloid cells in NRs suffer from a compromised capacity to sense tumor antigens and prime cytotoxic T cells under ICB pressure (Figure S2G,S2H).

These results indicate that successful ICB therapy drives profound and directed differentiation and functional remodeling of myeloid cells in REs, characterized by coordinated metabolic reprogramming, alleviation of immune suppression checkpoints, and comprehensive enhancement of proinflammatory, effector, and antigen presentation functions. In contrast, myeloid cells in NRs failed to complete this remodeling process, remaining locked in a metabolically restricted and immunosuppressive state.

Macro_C1QC reveals prognostic relevance of myeloid remodeling in specific cancer types

The intricate dynamics of cell-cell interaction within the TME involve not only TIMs but also other immune cell populations, particularly T cells. To further elucidate the association between intercellular crosstalk and therapeutic response, we performed a systematic ligand-receptor interaction analysis on post-treatment samples. Specifically, we extended the scCURE workflow to CD8⁺ and CD4⁺ T cell lineages (Figure S3A,S3B). Within each lineage, we identified all “Changed cells”, followed by cross-sample integration and re-clustering. In total, we analyzed T cells derived from four different cancer types and annotated nine distinct subsets of CD4⁺ T cells and 11 subsets of CD8⁺ T cells based on their distinctive molecular signatures. Subsequently, we utilized LIANA to investigate the intercommunication among TIMs, as well as the reciprocal interactions between TIMs and CD4⁺/CD8⁺ T cells in the Post group. Strikingly, the Macro_C1QC subset exhibited the highest frequency of intercellular interactions within the NR group (Figure 3A). This prominence in interaction strength remained evident when specifically analyzing crosstalk between myeloid cells and either CD4⁺ or CD8⁺ T cells (Figure 3B). Notably, we identified a repertoire of inhibitory ligand-receptor pairs that were specifically enhanced or exclusively present in the NR group, including the MHC-I-LILR axes (HLA-A/B/C–LILRA1/3, B2M–LILRB1/2/3), complement-LILR interactions (C1QA/B–LILRB3), and APOE signaling pathways (APOE–LRP1/LDLR/SCARB1) (Figure 3C). These results suggest that within the TME of NR patients, a cellular communication network centered on Macro_C1QC and mediated by inhibitory axes, such as MHC-I-LILR, is significantly reinforced. This amplified crosstalk likely impedes the effective activation of anti-cancer T cells by concurrently disrupting antigen presentation and delivering inhibitory signals.

Figure 3 Specific enrichment and prognostic value of Macro_C1QC in ccRCC. (A) Scatter plot depicting the FC in signaling strength for interactions occurring between myeloid populations, comparing post-treatment NRs versus REs. The dot size corresponds to the absolute log10FC, and the color indicates the specific myeloid cell type. (B) Scatter plots illustrating the FC in interaction strength derived from myeloid subsets (sources/senders) to CD8⁺ T cells (targets/receivers, left) and CD4⁺ T cells (targets/receivers, right) between the post-NR and post-RE groups. (C) Bubble plot characterizing specific ligand-receptor signaling axes originating from Macro_C1QC and targeting various cell types within the post-treatment NR cohort. The dot size reflects statistical significance (interaction specificity), and the color gradient represents the magnitude of expression for the ligand-receptor pair. (D) Kaplan-Meier curves for OS in an independent cohort of ICB-treated ccRCC patients. Patients were stratified into High (red) and Low (blue) expression groups based on the Macro_C1QC gene signature. (E) Kaplan-Meier survival analysis assessing the association between Macro_C1QC signature expression and overall survival in the TCGA ccRCC cohort. ccRCC, clear cell renal cell carcinoma; cDC, conventional dendritic cell; FC, fold change; ICB, immune checkpoint blockade; NR, non-responder; OS, overall survival; pDC, plasmacytoid dendritic cell; RE, responder; TCGA, The Cancer Genome Atlas.

Building upon our analysis of preferential distribution across clinical strata, we observed that Macro_C1QC was markedly enriched in ccRCC patients. Notably, within the NR group, this subset significantly increased after immunotherapy. To validate the clinical consequences of Macro_C1QC-mediated immune suppression, we evaluated the prognostic value of its gene signature in ccRCC patients. In an independent cohort of ccRCC patients treated with ICB, high expression of the Macro_C1QC signature was significantly related to shortened OS (Figure 3D). This result was further corroborated in the larger TCGA ccRCC cohort, where a high Macro_C1QC score similarly predicted significantly poorer clinical outcomes (Figure 3E).

These findings indicate that some myeloid subsets related to treatment resistance may show prognostic relevance in selected cancer types. In particular, Macro_C1QC appears to be preferentially enriched in ccRCC and is related to inhibitory intercellular communication and unfavorable clinical outcomes. However, this pattern does not indicate that Macro_C1QC serves as a major resistance mechanism in HNSCC. Therefore, we next focused on Macro_BIRC5, another macrophage subset enriched in NRs, which showed broader features related to treatment resistance and greater relevance to the HNSCC-focused narrative of this study.

The Macro_BIRC5 subset promotes ICB resistance via hyper-oxidative metabolism and establishment of a synergistic immunosuppressive network

After identifying Macro_C1QC as a key immunosuppressive myeloid subset enriched in NRs, we further focused on another myeloid population that is specifically expanded in the NR group, Macro_BIRC5 (Figure 1D). This subset exhibits pronounced therapy-associated dynamic changes across multiple cancer types, yet displays functional characteristics that are clearly distinct from those of Macro_C1QC. To avoid cell-number-driven bias, we quantified the proportion of Macro_BIRC5 cells at the patient level in post-treatment samples (Figure 4A). We observed that the proportion of Macro_BIRC5 was significantly higher in post-treatment NRs (post-NR) compared with REs (post-RE) (P=0.006), whereas this population was nearly absent in post-R patients. KEGG pathway enrichment analysis revealed that genes highly expressed in Macro_BIRC5 were concentrated in pathways related to cell junctions (“Tight junction”, “Gap junction”) and cell survival and stress (“Apoptosis”, “p53 signaling pathway”) (Figure 4B), suggesting unique anti-apoptotic capabilities and tissue remodeling potential for this subset.

Figure 4 The Macro_BIRC5 subset promotes ICB resistance via hyper-oxidative metabolism and establishment of a synergistic immunosuppressive network. (A) Boxplot showing the proportion of Macro_BIRC5 cells among myeloid populations in post-treatment NRs (post_NR) and REs (post_R). Each dot represents an individual patient. (B) Dot plot illustrating the top enriched KEGG pathways for genes upregulated in the Macro_BIRC5 subset. Dot size represents the gene count, and the color gradient indicates statistical significance. (C) Boxplots displaying the expression levels of BIRC5 in tumor versus matched normal tissues across TCGA HNSCC, ccRCC, SKCM, LUAD and LUSC cohorts. (D) Kaplan-Meier survival analysis for OS in three independent immunotherapy cohorts (SKCM, ccRCC, and NSCLC). Patients were stratified based on the Macro_BIRC5 gene signature expression. (E) Boxplots comparing Immune Scores and Stromal Scores between Macro_BIRC5-high and -low expression groups in the SKCM cohort. (F) Violin plots illustrating the correlation between Macro_BIRC5 signature scores and clinical tumor stages in TCGA LUAD and ccRCC cohorts. (G) Heatmap displaying the relative activity scores of major metabolic pathways across distinct myeloid subsets. The color gradient represents scaled metabolic pathway activity scores (Z-scores). (H) Metabolic disparity in therapeutic response. Violin plots comparing OXPHOS activity scores in Macro_BIRC5 cells between post-treatment NRs and REs. (I) Scatter plot showing the correlation between OXPHOS activity score and immunosuppressive score in Macro_BIRC5 cells. (J) Heatmap displaying the inferred regulatory activity of key TFs across myeloid subsets. The color gradient represents scaled TF activity scores (Z-scores). (K) Bubble plot characterizing specific ligand-receptor interactions originating from Macro_BIRC5 and targeting other immune cell types. (L,M) Network plot (L) and interaction strength heatmap (M) visualizing the intercellular crosstalk among key resistance-associated myeloid subsets. *, P<0.05; **, P<0.01. BIRC5, baculoviral IAP repeat containing 5; ccRCC, clear cell renal cell carcinoma; HNSCC, head and neck squamous cell carcinoma; ICB, immune checkpoint blockade; KEGG, Kyoto Encyclopedia of Genes and Genomes; LUAD, lung adenocarcinoma; LUSC, lung squamous cell carcinoma; NR, non-responder; NSCLC, non-small cell lung cancer; OS, overall survival; OXPHOS, oxidative phosphorylation; RE, responder; SKCM, skin cutaneous melanoma; TCGA, The Cancer Genome Atlas; TF, transcription factor; TPM, transcripts per million; Treg, regulatory T cell.

To evaluate the clinical relevance of Macro_BIRC5, we examined the relationship between the expression levels of its signature gene set, TME characteristics, and prognosis. Firstly, Macro_BIRC5 gene signatures expression was significantly higher in tumor tissues compared to normal tissues across multiple TCGA cancer cohorts (SKCM, ccRCC, LUSC) (Figure 4C). More importantly, across independent clinical cohorts receiving ICB therapy, high Macro_BIRC5 gene signature scores were related to shorter OS in SKCM, ccRCC, and NSCLC (Figure 4D; P=0.006, P<0.001, and 0.009, respectively). In SKCM cohorts, the Macro_BIRC5-high expression group exhibited lower Immune Scores, indicating a close association with immunosuppressive microenvironment characterized by reduced immune cell infiltration (Figure 4E). Furthermore, in ccRCC and LUAD datasets, the Macro_BIRC5 score increased with advancing tumor stage, further suggesting its potential role in tumor progression (Figure 4F).

Next, we analyzed the metabolic profile of Macro_BIRC5. A heatmap based on metabolic pathway gene set scores showed that Macro_BIRC5 exhibited high activity in pathways related to both OXPHOS and glucose metabolism, a metabolic state significantly distinct from other myeloid subsets (Figure 4G), indicating a distinct hyper-metabolic state. Notably, in post-ICB treatment samples, Macro_BIRC5 cells from NRs displayed higher OXPHOS activity than those from REs (Figure 4H), suggesting that metabolic reprogramming may be part of its resistance mechanism. To further investigate the functional relevance of the metabolic state of Macro_BIRC5, we assessed the relationship between OXPHOS activity and immunosuppressive gene expression at the single-cell level. Pearson correlation analysis revealed a statistically significant but modest positive correlation between OXPHOS activity and the immunosuppressive score in Macro_BIRC5 cells (R=0.27, P=0.005), suggesting that enhanced oxidative metabolism is related to increased immunosuppressive potential (Figure 4I).

To explore how this metabolic state may be linked to transcriptional regulation, we inferred TF activities using the DoRothEA regulon resource coupled with the VIPER algorithm. Macro_BIRC5 exhibited a distinct transcriptional regulatory profile characterized by the specific hyperactivation of E2F3, TFDP1, SREBF2, NR2F6, and KLF4. Among these, E2F3 and TFDP1 are well-known regulators of the cell cycle, aligning with the proliferative characteristics of this subset, while SREBF2 suggests a reprogramming of cholesterol metabolism (Figure 4J). These results suggest that the hyper-oxidative state of Macro_BIRC5 is accompanied by transcriptional programs related to proliferation and metabolic adaptation. Crucially, our analysis highlighted a robust immunosuppressive program in Macro_BIRC5, orchestrated by NR2F6 and KLF4. NR2F6 acts as an intracellular immune checkpoint that suppresses pro-inflammatory cytokine transcription, while KLF4 is a well-established inducer that reinforces the immunosuppressive phenotype. Together with the positive association between OXPHOS and immunosuppressive scores, these TF activity patterns support a model in which Macro_BIRC5 couples oxidative metabolic remodeling with immunosuppressive transcriptional regulation.

Next, we utilized LIANA to explore the specific signals emanating from Macro_BIRC5 and their potential immunoregulatory functions (Figure 4K). Given that Macro_BIRC5 preferentially accumulates in NRs post-ICB treatment and is closely related to poor prognosis, we considered it crucial to understand its functional interactions within this specific resistant microenvironment. Therefore, to precisely capture active signaling pathways driving or maintaining immune escape under therapeutic pressure, we specifically selected cells from the “Post-NR” group for communication network analysis. This focused analytical strategy allowed us to isolate key communication events mediated by Macro_BIRC5 that are most active in the context of therapeutic failure. Notably, the analysis revealed strong directional communication links between Macro_BIRC5 and Treg. Macro_BIRC5 interacts with receptors on the membrane of Treg (such as TSHR) by expressing various ligands (such as FN1, POMC, ADM), and also interacts through the PSEN1-NOTCH2 axis. Given that Treg are a key immunosuppressive population in the TME, the active signals sent to them by Macro_BIRC5 suggest a role in enhancing Treg function, thereby consolidating the immune-tolerant microenvironment.

These results further connect the metabolic and transcriptional features of Macro_BIRC5 with its Treg-directed communication pattern. The elevated OXPHOS activity may support the high-energy demand of this proliferative and stress-adaptive macrophage state, while activated TFs such as E2F3 and TFDP1 are consistent with cell-cycle activity, and SREBF2 suggests accompanying lipid/cholesterol metabolic remodeling. In parallel, the activation of NR2F6 and KLF4 indicates that this metabolically active state is related to an immunosuppressive transcriptional program. This program may provide a permissive context for the expression of Treg-interacting ligands, including FN1, POMC, ADM, and PSEN1, which engage Treg-associated receptors such as TSHR and NOTCH2. Therefore, the Macro_BIRC5-Treg communication pattern is consistent with a coordinated resistance-associated phenotype in which oxidative metabolism, immunosuppressive transcriptional regulation, and Treg-oriented signaling are linked.

Finally, to explore the broader immune landscape related to Macro_BIRC5, we extended the scCURE framework to other immune lineages, including CD8⁺ T cells, CD4⁺ T cells, NK cells and B/plasma cells (Figure S3). Within each lineage, we identified therapy-driven “Changed cells”. Subsequently, we integrated and re-clustered these cells to resolve functionally refined subsets. We systematically identified cell subsets preferentially enriched in post-treatment NRs based on the Ro/e ratio; their gene signatures were significantly related to poor clinical prognosis in at least three cancer types. Ultimately, we determined that CD8T_RRM2, Treg_1, Plasma_Conventional, Macro_BIRC5, and Macro_SPP1 constitute key resistance-associated cell subsets (Figure S4). We constructed a critical communication network among these five subsets. To decipher the cooperative interactions among these subsets, we reconstructed their communication network using CellChat. Analysis of both network topology (Figure 4L) and interaction frequency (Figure 4M) consistently identified Macro_BIRC5 as the central communication hub. This underscores its role as a dominant mediator of intercellular crosstalk, driving the coordination of this immunosuppressive microenvironment.

The Macro_BIRC5-Treg axis is validated as resistance mechanism in HNSCC

Given the pan-cancer identification of Macro_BIRC5 as a resistance-associated subset with extensive Treg interactions, we conducted focused analyses to validate its clinical relevance specifically in HNSCC.

We first examined the expression pattern of Macro_BIRC5 signature genes in HNSCC clinical samples. Analysis of the TCGA-HNSCC cohort revealed that Macro_BIRC5 gene signatures were significantly elevated in tumor tissues compared to normal tissues (Figure 5A; P=0.03). Furthermore, we observed a progressive increase in Macro_BIRC5 signature scores with advancing tumor stage (Figure 5B; P=0.04), suggesting that the accumulation of this macrophage subset may contribute to HNSCC disease progression. To assess the relationship between Macro_BIRC5 and the immune landscape of HNSCC, we evaluated its association with TME characteristics. Consistent with our pan-cancer observations, HNSCCs with high Macro_BIRC5 signature expression exhibited lower Immune Scores and higher Stromal Scores in an HNSCC cohort treated with ICB therapy (Figure 5C), suggesting that this signature may be linked to an immunosuppressive, stromal-enriched microenvironment with reduced ICB efficacy. Critically, in this independent cohort, high Macro_BIRC5 signature expression was significantly related to shortened OS (Figure 5D; P=0.02), establishing its prognostic value in ICB-treated HNSCC patients.

Figure 5 The Macro_BIRC5-Treg axis is validated as resistance mechanism in HNSCC. (A) Boxplots displaying the expression levels of Macro_BIRC5 in tumor versus normal tissues in TCGA HNSCC cohort. (B) Violin plots illustrating the correlation between Macro_BIRC5 signature scores and clinical tumor stages in TCGA HNSCC cohort. (C) Boxplots comparing Immune Scores and Stromal Scores between Macro_BIRC5-high and -low expression groups in an independent immunotherapy cohort of HNSCC. (D) Kaplan-Meier survival analysis for OS in the HNSCC cohort. Patients were stratified based on the Macro_BIRC5 gene signature expression. (E) Spatial feature plots visualizing the deconvolution-inferred abundance of the Macro_BIRC5 subset across four representative HNSCC tissue sections. The color gradient indicates the spatial score intensity. The spatial resolution of the 10× Genomics Visium platform is approximately 55 μm per capture spot. (F,G) Scatter plots with linear regression lines depicting the correlation between BIRC5 expression and microenvironment scores (Stromal/Immune) across spatial spots in (F) HNSCC and (G) SKCM datasets. Pearson correlation coefficients and P values are indicated. (H) Overlay of Macro_BIRC5 signature scores and FOXP3 expression on a single HNSCC tissue section. The arrows indicate representative regions showing spatial co-localization of Macro_BIRC5 and FOXP3 signals. The spatial resolution of the 10× Genomics Visium platform is approximately 55 μm per capture spot. (I) Smoothed expression profiles illustrating the spatial congruency of scaled Macro_BIRC5 and FOXP3 values along the defined tissue coordinates. *, P<0.05. BIRC5, baculoviral IAP repeat containing 5; FOXP3, forkhead box P3; HNSCC, head and neck squamous cell carcinoma; N, normal tissue; OS, overall survival; SKCM, skin cutaneous melanoma; T, tumor; TCGA, The Cancer Genome Atlas.

The spatial architecture of the TME fundamentally shapes therapeutic response, where the physical juxtaposition of myeloid cells with specific lymphocyte subsets often dictates the efficacy of anti-tumor immunity. Having inferred Macro_BIRC5 as a central communication hub from scRNA-seq, we leveraged ST data across multiple cancer types to rigorously validate its spatial distribution and physical interacting partners in situ.

We then mapped the spatial distribution of Macro_BIRC5 using deconvolution-based scoring. As shown in Figure 5E, Macro_BIRC5 signatures were significantly enriched within the annotated tumor regions of HNSCC samples, confirming its tumor-resident nature. Next, to assess the broader microenvironmental context of this population, we correlated Macro_BIRC5 signature expression with tumor purity indices in HNSCC and SKCM cohorts. Notably, BIRC5 expression was positively correlated with Stromal Scores in both HNSCC and SKCM cohorts (HNSCC: R=0.13, P=0.04; SKCM: R=0.32, P<0.001), but negatively correlated with Immune Scores (HNSCC: R=−0.15, P=0.02; SKCM: R=−0.16, P<0.001; Figure 5F,5G). These results suggest that higher BIRC5 expression may reflect a stromal-enriched and relatively immune-low TME.

To further spatially resolve the intercellular crosstalk identified by scRNA-seq, we mapped the distribution of Macro_BIRC5 relative to Treg. We observed a striking physical juxtaposition of these two populations across multiple tissue sections in HNSCC. Quantitative trend analysis further revealed a highly synchronized spatial distribution pattern, where the scaled abundance of Macro_BIRC5 exhibited strong spatial congruency with FOXP3 expression (Figure 5H,5I). This spatial co-localization provides robust physical evidence for the formation of a localized Macro_BIRC5-Treg immunosuppressive axis, fostering their local cooperative interactions.

These spatially resolved results are highly congruent with our scRNA-seq results regarding cellular interactions, metabolic features, and clinical prognosis. Collectively, they indicate that in HNSCC, Macro_BIRC5 is not only enriched in tumor domains with low immune infiltration but also actively reinforces the immune-tolerant microenvironment through precise spatial co-localization with Treg, thereby promoting resistance to ICB therapy.

Discussion

TIMs are increasingly recognized as pivotal regulators of the TME and key determinants of immune checkpoint blockade efficacy. However, resolving their fine-grained heterogeneity and therapy-induced plasticity has remained challenging, largely due to the confounding effects of steady-state signals in conventional single-cell analyses. In this study, we addressed this limitation by applying scCURE, a novel computational framework that distinguishes therapy-driven “Changed” myeloid cells from homeostatic “Steady” populations. Applied to a longitudinal pan-cancer atlas of 82 patients with pre- and post-ICB samples across four malignancies, scCURE enabled precise interrogation of bona fide transcriptional remodeling under therapeutic pressure. Our analysis uncovered three major advances: (I) markedly divergent evolutionary trajectories of TIMs in REs versus NRs; (II) a previously unrecognized Macro_BIRC5-Treg axis that orchestrates resistance across cancer types; and (III) specific validation of this axis as resistance mechanism and therapeutic predictive marker in HNSCC.

A central discovery of this work is the divergent remodeling trajectories that distinguish clinical outcomes. In REs, TIMs executed a coordinated metabolic switch toward enhanced aerobic glycolysis and OXPHOS, accompanied by downregulation of immunosuppressive checkpoints (e.g., HAVCR2, IDO1) and upregulation of proinflammatory and antigen-presentation programs. This reprogramming mirrors the energetic and functional demands of activated effector T cells, suggesting that successful ICB requires myeloid cells to transition into a supportive, high-energy state. In contrast, NRs exhibited arrested differentiation, persistent immunosuppressive markers, and metabolic collapse—indicative of a “locked” state resistant to therapeutic reprogramming. These results have implications for HNSCC, where baseline myeloid dysfunction is already pronounced. In HNSCC tumors, myeloid cells are heavily infiltrated and exhibit constitutive immunosuppressive features including high expression of arginase-1, IDO1, and PD-L1, which collectively suppress T cell receptor expression and antigen-specific T cell responses (13,14). The failure of myeloid metabolic reprogramming observed in our NR cohort suggests that in HNSCC, therapeutic strategies must not only block inhibitory checkpoints but also actively restore myeloid metabolic fitness to overcome the deeply entrenched immunosuppressive state. This may explain why single-agent ICB achieves limited efficacy in HNSCC and supports the rationale for metabolic combination approaches. These results implicate metabolic interventions (e.g., promoting glycolysis or inhibiting lipid/OXPHOS pathways) as promising strategies to unlock myeloid plasticity and overcome resistance.

We further identified two resistance-associated macrophage subsets with distinct mechanisms. Macro_C1QC, preferentially enriched in non-responding ccRCC patients, functioned as a hub for inhibitory crosstalk, notably via MHC-I-LILR axes (10), and its signature robustly predicted poor survival in independent cohorts. Most notably, we characterized Macro_BIRC5 as a proliferative, hyper-oxidative subset that drives pan-cancer resistance. As a central hub in the resistance network, it sustained extensive crosstalk with Treg through multiple ligand-receptor pairs (e.g., FN1-TSHR, PSEN1-NOTCH2) (47,48). In HNSCC specifically, we observed that Macro_BIRC5 gene signature expression was significantly elevated in tumor tissues compared to normal tissues, increased progressively with tumor stage, and strongly predicted poor OS in ICB-treated patients. ST analysis of HNSCC specimens validated the tumor-specific enrichment of Macro_BIRC5 and its physical co-localization with Treg within the TME, providing direct spatial evidence for this resistance-associated cellular interaction in HNSCC. The consistent association of Macro_BIRC5 signatures with adverse prognosis, advanced stage, and reduced immune infiltration across multiple cancers positions this subset as both a predictive biomarker and a high-priority therapeutic target (49).

Importantly, these features of Macro_BIRC5 may be viewed as an integrated resistance-associated program rather than as separate observations. Elevated OXPHOS may support the energetic demands of a proliferative and stress-adaptive macrophage state, whereas increased activity of E2F3, TFDP1, and SREBF2 is consistent with cell-cycle progression and lipid/cholesterol metabolic adaptation. In parallel, activation of NR2F6 and KLF4 may link this metabolically active state to immunosuppressive macrophage polarization. Within this transcriptional context, Macro_BIRC5 expressed multiple ligands predicted to communicate with Treg, including FN1, POMC, ADM, and PSEN1, which were linked to Treg-associated receptors such as TSHR and NOTCH2. Thus, the Macro_BIRC5-Treg axis may reflect a coordinated metabolic-transcriptional-communication program in which oxidative metabolic remodeling supports an immunosuppressive macrophage state that reinforces Treg-centered immune tolerance. However, the causal hierarchy among metabolism, TF activity, and ligand-receptor communication remains to be experimentally determined.

From a therapeutic perspective, targeting the metabolic reprogramming of myeloid-Treg interactions may represent a feasible but challenging strategy to overcome ICB resistance. The feasibility of this approach is supported by increasing evidence that myeloid cells are therapeutically targetable in solid tumors and that metabolic programs can shape the suppressive functions of both macrophages and tumor-infiltrating Treg. In our study, Macro_BIRC5 exhibited elevated OXPHOS activity, immunosuppressive TF programs, and Treg-directed ligand-receptor interactions, suggesting several potential intervention points. These may include attenuating oxidative metabolic programs in suppressive macrophages, disrupting Macro_BIRC5-Treg communication axes such as FN1-TSHR or PSEN1-NOTCH2, or combining myeloid metabolic modulation with ICB to restore anti-tumor immunity.

However, several challenges should be considered. First, metabolic pathways such as OXPHOS, glycolysis, and lipid metabolism are broadly used by tumor cells and immune cells, making cell-type-specific targeting difficult and raising the possibility of systemic toxicity or unintended impairment of effector immune cells. Second, because myeloid cells and Treg-associated communication pathways are both highly heterogeneous, plastic, and involved in immune homeostasis, non-selective metabolic inhibition or broad blockade of these interactions may impair beneficial antigen-presenting/phagocytic myeloid subsets and increase immune-related adverse effects. Finally, the current study is based mainly on computational inference and spatial association, and functional experiments are required to determine whether disrupting Macro_BIRC5 metabolism or its communication with Treg can enhance ICB efficacy. Thus, future therapeutic development should prioritize cell-state-specific markers, selective delivery systems, and rational combination strategies that selectively attenuate immunosuppressive myeloid-Treg communication while preserving productive anti-tumor immunity.

Although our results were validated across independent clinical cohorts and spatially corroborated, several limitations remain. First, this study was based on retrospective integration of publicly available scRNA-seq datasets; therefore, a prospective sample size calculation was not applicable. Although the final cohort of 82 patients and 129 tumor samples is comparable to recent pan-cancer single-cell studies of TIMs during immunotherapy, the sample size remains limited for definitive cancer type-specific clinical stratification (50). Future studies using larger, prospectively collected longitudinal cohorts are required to further validate the robustness and clinical generalizability of the identified myeloid remodeling programs.

Conclusions

In conclusion, this longitudinal pan-cancer single-cell study reveals distinct therapy-associated remodeling patterns of TIMs during ICB therapy, with particular relevance to HNSCC. We identified divergent myeloid trajectories between REs and NRs, characterized by coordinated metabolic and functional reprogramming in REs and persistent immunosuppressive remodeling in NRs. We further identified Macro_BIRC5 as a resistance-associated macrophage state marked by hyper-oxidative metabolism, immunosuppressive transcriptional programs, and Treg-directed communication, together suggesting that the Macro_BIRC5-Treg axis may contribute to the formation of an immune-suppressive TME. These results provide a myeloid-centered framework for understanding ICB resistance and suggest that the Macro_BIRC5-related program may have potential value for predicting unfavorable immunotherapy response, particularly in HNSCC. Selective modulation of suppressive myeloid metabolic programs or myeloid-Treg communication may represent a potential direction for future combination strategies.

Acknowledgments

None.

Footnote

Peer Review File: Available at https://fomm.amegroups.com/article/view/10.21037/fomm-2026-1-0005/prf

Funding: This work was supported by the National Natural Science Foundation of China (grant Nos. 82230090 and 32570782 to X.Z.), the Shanghai Committee of Science and Technology (grant No. 21DZ2292000), the Shanghai Jiao Tong University Trans-med Awards Research (grant No. 20250203), and the Explorers Program of Shanghai (grant No. 25TS1403700).

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://fomm.amegroups.com/article/view/10.21037/fomm-2026-1-0005/coif). The authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. This study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments. All data analyzed in this study were obtained from publicly available databases (GEO, TCGA, and dbGaP). No new human or animal experiments were conducted, and therefore institutional ethical approval and informed consent were not required.

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Kraehenbuehl L, Weng CH, Eghbali S, et al. Enhancing immunotherapy in cancer by targeting emerging immunomodulatory pathways. Nat Rev Clin Oncol 2022;19:37-50. [Crossref] [PubMed]
2. Ferris RL, Blumenschein G Jr, Fayette J, et al. Nivolumab vs investigator's choice in recurrent or metastatic squamous cell carcinoma of the head and neck: 2-year long-term survival update of CheckMate 141 with analyses by tumor PD-L1 expression. Oral Oncol 2018;81:45-51. [Crossref] [PubMed]
3. Mehra R, Seiwert TY, Gupta S, et al. Efficacy and safety of pembrolizumab in recurrent/metastatic head and neck squamous cell carcinoma: pooled analyses after long-term follow-up in KEYNOTE-012. Br J Cancer 2018;119:153-9. [Crossref] [PubMed]
4. Kim H, Kwon M, Kim B, et al. Clinical outcomes of immune checkpoint inhibitors for patients with recurrent or metastatic head and neck cancer: real-world data in Korea. BMC Cancer 2020;20:727. [Crossref] [PubMed]
5. Cohen EEW, Soulières D, Le Tourneau C, et al. Pembrolizumab versus methotrexate, docetaxel, or cetuximab for recurrent or metastatic head-and-neck squamous cell carcinoma (KEYNOTE-040): a randomised, open-label, phase 3 study. Lancet 2019;393:156-67. [Crossref] [PubMed]
6. Aliazis K, Christofides A, Shah R, et al. The tumor microenvironment's role in the response to immune checkpoint blockade. Nat Cancer 2025;6:924-37. [Crossref] [PubMed]
7. Goswami S, Anandhan S, Raychaudhuri D, et al. Myeloid cell-targeted therapies for solid tumours. Nat Rev Immunol 2023;23:106-20. [Crossref] [PubMed]
8. Engblom C, Pfirschke C, Pittet MJ. The role of myeloid cells in cancer therapies. Nat Rev Cancer 2016;16:447-62. [Crossref] [PubMed]
9. Pan Y, Yu Y, Wang X, et al. Tumor-Associated Macrophages in Tumor Immunity. Front Immunol 2020;11:583084. [Crossref] [PubMed]
10. Barkal AA, Weiskopf K, Kao KS, et al. Engagement of MHC class I by the inhibitory receptor LILRB1 suppresses macrophages and is a target of cancer immunotherapy. Nat Immunol 2018;19:76-84. [Crossref] [PubMed]
11. Curiel TJ, Coukos G, Zou L, et al. Specific recruitment of regulatory T cells in ovarian carcinoma fosters immune privilege and predicts reduced survival. Nat Med 2004;10:942-9. [Crossref] [PubMed]
12. Rodriguez PC, Quiceno DG, Zabaleta J, et al. Arginase I production in the tumor microenvironment by mature myeloid cells inhibits T-cell receptor expression and antigen-specific T-cell responses. Cancer Res 2004;64:5839-49. [Crossref] [PubMed]
13. Mito I, Takahashi H, Kawabata-Iwakawa R, et al. Comprehensive analysis of immune cell enrichment in the tumor microenvironment of head and neck squamous cell carcinoma. Sci Rep 2021;11:16134. [Crossref] [PubMed]
14. Elmusrati A, Wang J, Wang CY. Tumor microenvironment and immune evasion in head and neck squamous cell carcinoma. Int J Oral Sci 2021;13:24. [Crossref] [PubMed]
15. Bhat AA, Yousuf P, Wani NA, et al. Tumor microenvironment: an evil nexus promoting aggressive head and neck squamous cell carcinoma and avenue for targeted therapy. Signal Transduct Target Ther 2021;6:12. [Crossref] [PubMed]
16. Li K, Shi H, Zhang B, et al. Myeloid-derived suppressor cells as immunosuppressive regulators and therapeutic targets in cancer. Signal Transduct Target Ther 2021;6:362. [Crossref] [PubMed]
17. Kumar Jha P, Aikawa M, Aikawa E. Macrophage Heterogeneity and Efferocytosis: Beyond the M1/M2 Dichotomy. Circ Res 2024;134:186-8. [Crossref] [PubMed]
18. Yunna C, Mengru H, Lei W, et al. Macrophage M1/M2 polarization. Eur J Pharmacol 2020;877:173090. [Crossref] [PubMed]
19. Jing ZQ, Luo ZQ, Chen SR, et al. Heterogeneity of myeloid cells in common cancers: Single cell insights and targeting strategies. Int Immunopharmacol 2024;134:112253. [Crossref] [PubMed]
20. He D, Wang D, Lu P, et al. Single-cell RNA sequencing reveals heterogeneous tumor and immune cell populations in early-stage lung adenocarcinomas harboring EGFR mutations. Oncogene 2021;40:355-68. [Crossref] [PubMed]
21. Mu Q, Zhang H, Shi Y, et al. Single-cell transcriptomic profiling reveals the heterogeneity of epithelial cells in lung adenocarcinoma lymph node metastasis and develops a prognostic signature. Front Immunol 2025;16:1637625. [Crossref] [PubMed]
22. He ZN, Zhang CY, Zhao YW, et al. Regulation of T cells by myeloid-derived suppressor cells: emerging immunosuppressor in lung cancer. Discov Oncol 2023;14:185. [Crossref] [PubMed]
23. Wang Y, Huang T, Gu J, et al. Targeting the metabolism of tumor-infiltrating regulatory T cells. Trends Immunol 2023;44:598-612. [Crossref] [PubMed]
24. MacNabb BW, Tumuluru S, Chen X, et al. Dendritic cells can prime anti-tumor CD8+ T cell responses through major histocompatibility complex cross-dressing. Immunity 2022;55:982-997.e8. [Crossref] [PubMed]
25. Kirschenbaum D, Xie K, Ingelfinger F, et al. Time-resolved single-cell transcriptomics defines immune trajectories in glioblastoma. Cell 2024;187:149-165.e23. [Crossref] [PubMed]
26. Yost KE, Satpathy AT, Wells DK, et al. Clonal replacement of tumor-specific T cells following PD-1 blockade. Nat Med 2019;25:1251-9. [Crossref] [PubMed]
27. Bi K, He MX, Bakouny Z, et al. Tumor and immune reprogramming during immunotherapy in advanced renal cell carcinoma. Cancer Cell 2021;39:649-661.e5. [Crossref] [PubMed]
28. Au L, Hatipoglu E, Robert de Massy M, et al. Determinants of anti-PD-1 response and resistance in clear cell renal cell carcinoma. Cancer Cell 2021;39:1497-1518.e11. [Crossref] [PubMed]
29. van der Leun AM, Traets JJH, Vos JL, et al. Dual Immune Checkpoint Blockade Induces Analogous Alterations in the Dysfunctional CD8+ T-cell and Activated Treg Compartment. Cancer Discov 2023;13:2212-27. [Crossref] [PubMed]
30. Luoma AM, Suo S, Wang Y, et al. Tissue-resident memory and circulating T cells are early responders to pre-surgical cancer immunotherapy. Cell 2022;185:2918-2935.e29. [Crossref] [PubMed]
31. Sade-Feldman M, Yizhak K, Bjorgaard SL, et al. Defining T Cell States Associated with Response to Checkpoint Immunotherapy in Melanoma. Cell 2018;175:998-1013.e20. [Crossref] [PubMed]
32. Zou X, Liu Y, Wang M, et al. scCURE identifies cell types responding to immunotherapy and enables outcome prediction. Cell Rep Methods 2023;3:100643. [Crossref] [PubMed]
33. Masarwy R, Kampel L, Horowitz G, et al. Neoadjuvant PD-1/PD-L1 Inhibitors for Resectable Head and Neck Cancer: A Systematic Review and Meta-analysis. JAMA Otolaryngol Head Neck Surg 2021;147:871-8. [Crossref] [PubMed]
34. Rysz J, Ławiński J, Franczyk B, et al. Immune Checkpoint Inhibitors in Clear Cell Renal Cell Carcinoma (ccRCC). Int J Mol Sci 2025;26:5577. [Crossref] [PubMed]
35. Karlsson AK, Saleh SN. Checkpoint inhibitors for malignant melanoma: a systematic review and meta-analysis. Clin Cosmet Investig Dermatol 2017;10:325-39. [Crossref] [PubMed]
36. Stratigos AJ, Sekulic A, Peris K, et al. Cemiplimab in locally advanced basal cell carcinoma after hedgehog inhibitor therapy: an open-label, multi-centre, single-arm, phase 2 trial. Lancet Oncol 2021;22:848-57. [Crossref] [PubMed]
37. Foy JP, Karabajakian A, Ortiz-Cuaran S, et al. Datasets for gene expression profiles of head and neck squamous cell carcinoma and lung cancer treated or not by PD1/PD-L1 inhibitors. Data Brief 2022;44:108556. [Crossref] [PubMed]
38. Choueiri TK, Fishman MN, Escudier B, et al. Immunomodulatory Activity of Nivolumab in Metastatic Renal Cell Carcinoma. Clin Cancer Res 2016;22:5461-71. [Crossref] [PubMed]
39. Motzer RJ, Rini BI, McDermott DF, et al. Nivolumab for Metastatic Renal Cell Carcinoma: Results of a Randomized Phase II Trial. J Clin Oncol 2015;33:1430-7. [Crossref] [PubMed]
40. Braun DA, Hou Y, Bakouny Z, et al. Interplay of somatic alterations and immune infiltration modulates response to PD-1 blockade in advanced clear cell renal cell carcinoma. Nat Med 2020;26:909-18. [Crossref] [PubMed]
41. Motzer RJ, Rini BI, McDermott DF, et al. Nivolumab plus ipilimumab versus sunitinib in first-line treatment for advanced renal cell carcinoma: extended follow-up of efficacy and safety results from a randomised, controlled, phase 3 trial. Lancet Oncol 2019;20:1370-85. [Crossref] [PubMed]
42. Patil NS, Nabet BY, Müller S, et al. Intratumoral plasma cells predict outcomes to PD-L1 blockade in non-small cell lung cancer. Cancer Cell 2022;40:289-300.e4. [Crossref] [PubMed]
43. Mariathasan S, Turley SJ, Nickles D, et al. TGFβ attenuates tumour response to PD-L1 blockade by contributing to exclusion of T cells. Nature 2018;554:544-8. [Crossref] [PubMed]
44. Liu D, Schilling B, Liu D, et al. Integrative molecular and clinical modeling of clinical outcomes to PD1 blockade in patients with metastatic melanoma. Nat Med 2019;25:1916-27. [Crossref] [PubMed]
45. Riaz N, Havel JJ, Makarov V, et al. Tumor and Microenvironment Evolution during Immunotherapy with Nivolumab. Cell 2017;171:934-949.e16. [Crossref] [PubMed]
46. Mulder K, Patel AA, Kong WT, et al. Cross-tissue single-cell landscape of human monocytes and macrophages in health and disease. Immunity 2021;54:1883-1900.e5. [Crossref] [PubMed]
47. Kuang DM, Peng C, Zhao Q, et al. Activated monocytes in peritumoral stroma of hepatocellular carcinoma promote expansion of memory T helper 17 cells. Hepatology 2010;51:154-64. [Crossref] [PubMed]
48. Benaiges E, Ceperuelo-Mallafré V, Madeira A, et al. Survivin drives tumor-associated macrophage reprogramming: a novel mechanism with potential impact for obesity. Cell Oncol (Dordr) 2021;44:777-92. [Crossref] [PubMed]
49. Yang S, Liu X, Mao S, et al. BIRC5 expression correlates with immunosuppressive phenotype and predicts inferior response to immunotherapy in lung adenocarcinoma. J Thorac Dis 2025;17:5921-35. [Crossref] [PubMed]
50. Li W, Pan L, Hong W, et al. A single-cell pan-cancer analysis to show the variability of tumor-infiltrating myeloid cells in immune checkpoint blockade. Nat Commun 2024;15:6142. [Crossref] [PubMed]