LDRT combined with anti-PD-1 promotes tumor regression in head and neck squamous cell carcinoma through tumor metaprogram evolution and immune cell reprogramming

TL;DR

Low-dose radiotherapy combined with anti-PD-1 therapy improves tumor regression in head and neck cancer by altering tumor and immune cell behavior.

Contribution

This study reveals how LDRT enhances αPD-1 efficacy through tumor and immune cell reprogramming in HNSCC.

Findings

LDRT and αPD-1 together reduce tumor growth and reshape the tumor microenvironment.

Combination therapy increases CD8⁺ T-cell clonal expansion and reduces CD4⁺ Treg-Ctla4 frequency.

An Isg15high macrophage state is linked to CD4⁺ Treg-Ctla4 and is reduced by LDRT-containing regimens.

Abstract

Anti–PD-1 (αPD-1) therapy shows limited efficacy in head and neck squamous cell carcinoma (HNSCC), partly due to an immunosuppressive tumor microenvironment (TME) and insufficient infiltration of effector T cells. Low-dose radiotherapy (LDRT) has been suggested to modify the TME, potentially enhancing the effects of αPD-1.

We profiled murine HNSCC treated with LDRT (1 Gy) and/or αPD-1 using single-cell RNA sequencing (scRNA-seq) and single-cell TCR sequencing (scTCR-seq), complemented by orthogonal functional assays and in situ validation (spatial transcriptomics and multiplex immunofluorescence).

LDRT and αPD-1 synergistically suppressed tumor growth and remodeled the TME. Single-cell analyses showed malignant programs shifting from proliferative states toward stress- and interferon (IFN)-responsive states. Integrated scTCR-seq analysis revealed enhanced clonal expansion of tumor-reactive CD8⁺ T-cell subsets and reduced the frequency of CD4⁺ Treg-Ctla4 under combination therapy. TCR-based metrics (including STARTRAC) indicated higher migration potential and limited intratumoral clonal overlap for CD4⁺ Treg-Ctla4. An Isg15high macrophage state (Macro-Isg15) was positively associated with CD4⁺ Treg-Ctla4 frequency and was more prominent under αPD-1 monotherapy. Ligand–receptor inference highlighted a Macro-Isg15–associated CXCL14–CXCR4 chemokine axis enriched under αPD-1 monotherapy and attenuated in LDRT-containing regimens, suggesting reduced chemokine-driven support for CD4⁺ Treg-Ctla4 accumulation/recruitment.

These findings support an integrated framework in which LDRT augments αPD-1 efficacy by reprogramming malignant states and reshaping immune cell dynamics, including attenuation of immunosuppressive features within the TME. This preclinical evidence provides a rationale for translational evaluation of LDRT-based combination strategies in HNSCC.

The online version contains supplementary material available at 10.1186/s13046-026-03660-3.

Figures

Single-cell sequencing reveals that LDRT in combination with αPD-1 reprograms the tumor microenvironment in HNSCC. a Experimental design and multi-modal profiling workflow in the SCC7 syngeneic model. Tumors received Control (IgG), LDRT (1 Gy × 3 fractions on days 1, 3, and 5), αPD-1 (10 mg/kg on days 1, 3, and 5), or the combination. Tumors were collected at indicated time points for IHC, NanoString, and flow cytometry, and at endpoint (day 21) for 10x Genomics scRNA-seq and scTCR-seq. Created with BioRender. b Representative photographs of excised tumors at day 21 from SCC7 (left) and MOC1 (

LDRT-driven tumor-cell state transition and immunogenic modulation underlying synergy with αPD-1. a UMAP plot showing the distribution of different tumor cell clusters (C1-C9). b GeneNMF-based metaprogram analysis identified four tumor cell metaprograms. c Expression levels of each metaprogram across different tumor cell clusters. d UMAP density plots of tumor cell clusters in different treatment groups, showing cell density distribution. e UMAP of tumor-cell states defined by metaprogram-associated features, annotated as SIRTC, ITC, and PTC. f Monocle3 trajectory analysis revealing differenti

c). These cytokine changes are consistent with an overall shift toward a more immune-active microenvironment. Flow cytometry further showed a marked increase in the frequency of intratumoral CD8⁺ T cells (Fig. [3](#Fig3)e-f, Supplementary Fig. 3b), whereas the frequency of CD4⁺ T cells remained relatively unchanged (Supplementary Fig. S3d). Notably, the proportions of CD3⁺CD8⁺ IFN-γ^+^ cytotoxic T cells (CTLs) and CD3⁺CD8⁺ TNF-α^+^ CTLs were also significantly increased in SCC7 tumors treated with LDRT + αPD-1 (Fig. [3](#Fig3)g). Together, these findings suggest that LDRT + αPD-1 increases the

Analysis of clonal expansion and state transitions in CD8 + and CD4 + T cells. a UMAP plot showing clonal expansion of T cell clusters, with different colors representing different clone sizes. b Distribution of clone sizes in different T cell clusters, showing the frequency of clonal expansion. c Heatmap of clonal overlap between different T cell clusters, revealing shared clonality among these populations. d STARTRAC analysis calculating expansion (expa), migration (migr), and transition (tran) indices for different T cell clusters to quantify their dynamic behaviors in the TME. e Scatter pl

TME cell subset composition dynamics and monocyte-macrophage state polarization under LDRT combined with αPD-1 treatment. a UMAP plots illustrating the distribution of tumor cells, T cells, NK cells, fibroblasts, endothelial cells, and myeloid cells. b CIBERSORTx was used to predict relative fractions of all cell clusters in TCGA-HNSC, followed by Pearson correlation analysis (*n* = 500). c Heatmap displaying proportional changes in monocyte-macrophage clusters based on the ratio of observed cell counts to expected cell counts (Ro/e), calculated using a chi-square test, across treatment groups

Spatial transcriptomics analysis of human HNSCC validates the spatial proximity of PTC-like tumor cells, Macro-Isg15, and Treg cells. a–c Representative spatial maps showing cell2location-inferred spot-level abundances of PTC-like tumor cells (C6), Macro-Isg15 macrophages, and CD4_Treg_Ctla4 cells in a human HNSCC spatial transcriptomics section. d–f Quantification of spatial proximity using SpottedPy based on spot-level cell2location abundance matrices. Getis–Ord Gi* (local spatial autocorrelation; neighborhood defined by a KNN/neighbor-number parameter) was applied to classify spots as hotsp

Regulation of ligand-receptor interactions among different cell types in the TME by LDRT combined with αPD-1 treatment. a Heatmaps of NicheNet analysis showing regulatory patterns between tumor cells and Macro-Isg15 among different treatment groups. b Heatmaps showing changes in ligand and receptor expression across different cell types, based on z-score normalization. c Representative GO biological process (BP) enrichment analysis illustrating the functional enrichment of predicted target genes expressed in Macro-Isg15 and tumor cells. d Heatmaps of NicheNet analysis showing regulatory patter

Full text

Background

Head and neck squamous cell carcinoma (HNSCC) is a prevalent malignancy, ranking sixth in global cancer incidence [1, 2]. The advent of immunotherapy, especially immune checkpoint inhibitors like anti-PD-1 (αPD-1) therapy, has marked a significant shift in the treatment paradigm for HNSCC. The KEYNOTE-048 clinical trial established the critical importance of anti-PD-1 therapy as a first-line treatment for recurrent and metastatic HNSCC [3]. At the same time, the CheckMate-141 clinical trial also established the efficacy of anti-PD-1 therapy in recurrent metastatic HNSCC [4]. However, the efficacy of these therapies is not universal, it is only effective in approximately 15% of patients with HNSCC [5], underscoring a sritical need to overcome therapeutic resistance.

The tumor microenvironment (TME) is a critical factor contributing to the resistance of HNSCC to immunotherapy. For instance, tumor-associated macrophages (TAMs) and regulatory T cells (Tregs) suppress the cytotoxic activity of T lymphocytes and promote tumor growth and metastasis by secreting immunosuppressive cytokines and modulating local inflammation. Additionally, tumor cells within the TME secrete immunosuppressive mediators and express immune checkpoint ligands, leading to the downregulation of human leukocyte antigen (HLA) expression [6, 7]. These findings underscore the importance of addressing both the immunosuppressive effects of cells within the TME and the genetic characteristics of tumor cells when designing immunotherapeutic strategies for HNSCC patients.

To overcome this resistance, combinatorial therapeutic strategies are being explored. Low-dose radiation therapy (LDRT), typically defined in immuno-oncology as non-ablative low per-fraction doses (e.g., 0.5–2 Gy per fraction), delivered as one or multiple fractions depending on the regimen [8, 9], has emerged as a promising approach to modulate the TME and enhance immune responses. LDRT has been shown to promote the reprogramming of immune cells within the TME, suggesting a synergistic potential when combined with immunotherapies. For instance, preclinical studies have indicated that LDRT (0.5–2 Gy) can reprogram immune-desert tumors by inducing M1 macrophage polarization and recruiting effector T cells [8], overcoming the inhibitory stroma [10], and reversing resistance to immunotherapy (Herrera et al., 2022). Therefore, LDRT effectively enhances the antitumor effects of immunotherapy [8–12].

The potential of LDRT to synergize with immune checkpoint inhibitors is an area of growing interest. In fact, our recent phase II clinical trial (NeoRTPC02) demonstrated that neoadjuvant LDRT combined with chemoimmunotherapy leads to a remarkably high pathological complete response rate (60.9%) in patients with resectable HNSCC [13]. However, the biological mechanisms underlying this powerful synergy remain incompletely understood, motivating us to investigate the mechanistic basis in preclinical models. To position our preclinical regimen relative to neoadjuvant clinical schedules, we used a low per-fraction dose of 1 Gy delivered in three fractions (days 1, 3, and 5) in the SCC7 model, consistent with the 1 Gy per-fraction dosing used in NeoRTPC02. We utilized single-cell RNA sequencing (scRNA-seq) and single-cell T cell receptor sequencing (scTCR-seq) of SCC7 mouse models and in vivo experiments to reveal that the combination of LDRT and αPD-1 treatment may lift the limitations on the efficacy of immunotherapy by mediating dynamic changes in the metaprogram of HNSCC tumor cells, promoting the differentiation of tumor-reactive CD8^+^ T cells, and inhibiting signal crosstalk between Macro-Isg15 and CD4⁺ Treg-Ctla4. This study provides a new perspective on the synergistic inhibition of HNSCC tumor progression by combining LDRT with αPD-1.

Methods

Mouse and cell lines

Female C3H mice, aged 6 weeks, were purchased from SPF (Beijing) Biotechnology Co., Ltd. All animal procedures were carried out in compliance with the guidelines established by the Institutional Animal Care and Use Committee at Southern Medical University. The mice were housed under specific pathogen-free (SPF) conditions with a 12-hour light/dark cycle, maintained at 18–22 °C, and with 50–60% humidity.

The SCC7 cell line, obtained from the National Collection of Authenticated Cell Cultures, was cultured in RPMI-1640 medium (Gibco) supplemented with 10% fetal bovine serum (FBS) (Gibco), and incubated at 37 °C in a 5% CO2 atmosphere.

Animal model

After one week of adaptive rearing, 10^6 SCC7 cells were subcutaneously injected into the dorsal region of C3H mice to establish the SCC7 tumor model. When the average tumor volume reached 100 mm³, the mice were randomly assigned to one of four treatment groups: IgG (control), low-dose radiation therapy (1 Gy locally delivered to the tumor, LDRT), anti-mouse PD-1 (BioLegend, 124314, 10 mg/kg), and a combination of anti-mouse PD-1 with LDRT, administered every 2 days for a total of 3 cycles. Tumor growth was monitored every three days using calipers, and tumor volume was calculated using the formula: Volume = (width)² × length/2. Mice were euthanized when the tumor volume reached 2000 mm³.

LDRT

Radiation therapy is directly applied to the tumor. Before the radiation therapy, each mouse was anesthetized and shielded with a lead chamber, exposing only the tumor that was to receive the irradiation. In brief, the mice were anesthetized with a dose of 10 mg/kg of pentobarbital prior to receiving ionizing radiation. Subsequently, local low-dose radiotherapy (LDRT) was delivered using an X-ray irradiator (Rad Source, Rs2000 Pro 225). Tumor-bearing mice were positioned to ensure localized irradiation, and the beam was collimated to a 4 × 4 cm field to cover the tumor area while minimizing exposure to surrounding tissues. The prescribed dose was 1 Gy per fraction (dose rate, 0.5 Gy/min; tube current, 4.3 mA) for a total of 3 fractions.

Flow cytometry analysis

Single-cell suspensions for flow cytometry analysis were prepared as previously described. Tumor samples were dissociated in PBS containing 1 mg/mL Collagen IV (Corning, 354233) and 1 mg/mL DNase I (ThermoFisher, EN0525) at 37 °C for 1 h. After digestion, the tissues were passed through a 70 μm mesh filter and washed with PBS to generate uniform single-cell suspensions. Dead cells were labeled using the LIVE/DEAD Fixable Aqua Dead Cell Stain Kit (ThermoFisher, L34957), following the manufacturer’s protocol. Fc receptors were blocked by incubating the cells with monoclonal antibodies for 15 min at 4 °C. For extracellular staining, the cells were incubated with fluorochrome-conjugated antibodies in 100 µL staining buffer (BioLegend, 420201) at 4 °C for 30 min. Following staining, the cells were washed twice and then fixed with Permeabilization Wash Buffer (BioLegend, 421002) before intracellular staining. For intracellular cytokine detection, cells were restimulated with 2 µL of PMA (BioLegend, 423303) in vitro. Flow cytometry sorting was conducted on a cytometer, and data analysis was performed using FlowJo software.

Antibodies listed for cell-marker staining were purchased from BioLegend: anti-mouse CD45 Alexa Fluor (147715), anti-mouse CD3 FITC (100203), anti-mouse CD8 PE (100707), anti-mouse Granzyme B PE/Cyanine7 (396409), anti-mouse CD11b FITC (101205), anti-mouse F4/80 APC (1231115), anti-mouse CD206 PE (141705), anti-mouse CD86 BV605 (105037), anti-mouse IFN-γ BV421 (505829), anti-mouse anti-TNF-α APC (506307).

Enzyme-linked immunosorbent assay (ELISA)

To measure cytokine levels, supernatants from Control, αPD-1, LDRT, and LDRT + αPD-1 treated samples were analyzed using ELISA kits. These kits specifically quantified mouse IFN-γ (Proteintech, KE10094), IL-12 (Proteintech, KE10068), and TNF-alpha (Proteintech, KE10002), following the manufacturer’s instructions. In brief, 100 µL of each sample was added to antibody-coated wells and incubated for 1 h at 37 °C. After the incubation, wells were washed three times with 350 µL of 1X wash buffer. TMB Development Solution (100 µL) was then added to each well, and the plates were incubated in the dark at 37 °C for 10 min. To halt the reaction, 100 µL of stop solution was added, and the optical density (OD) was subsequently measured at 450 nm with a microplate reader.

Tissue dissociation and preparation of single-cell suspensions

Mouse subcutaneous tumor tissue samples were collected and processed into single-cell suspensions. The tissues were placed into a sterile RNase-free culture dish with calcium- and magnesium-free PBS (1×) on ice and minced into approximately 0.5 mm² pieces. After washing with PBS to remove blood and adipose tissue, enzymatic dissociation was performed using a solution containing 0.35% collagenase IV, 2 mg/mL papain, and 120 U/mL DNase I. The mixture was incubated in a 37 °C water bath with shaking at 100 rpm for 20 min. To stop digestion, PBS containing 10% fetal bovine serum (FBS) was added, followed by gentle pipetting. The suspension was then passed through a 70 μm cell strainer, centrifuged at 300 g for 5 min at 4 °C, and resuspended in PBS with 0.04% BSA. Red blood cells were lysed using red blood cell lysis buffer (MACS 130-094−183, 1×), and dead cells were removed using the Miltenyi Dead Cell Removal Kit (MACS 130-090−101). Cell viability was assessed by trypan blue exclusion, and cells with viability greater than 85% were counted and adjusted to a concentration of 700–1200 cells/µL.

Single-cell library construction and sequencing

The Chromium Single-Cell 3’ kit (V3) from 10x Genomics was used to process single-cell suspensions, which were then loaded onto the Chromium Controller to capture approximately 5000 cells. Following the manufacturer’s instructions, cDNA synthesis and library construction were carried out. Sequencing of the libraries was performed on an Illumina NovaSeq 6000 platform using paired-end reads (150 bp read length) at a depth of at least 20,000 reads per cell.

V(D)J library construction and sequencing

Single-cell suspensions for TCR sequencing were processed with the Chromium Single-Cell 5’ V(D)J kit (10x Genomics). Gel bead-in-emulsion (GEM) technology was used to partition cells and barcode RNA. After cell lysis and barcoding within GEMs, reverse transcription was performed to generate barcoded cDNA. Full-length TCR transcripts were enriched by PCR amplification, followed by enzymatic fragmentation and size selection for V(D)J library construction. Libraries were sequenced using paired-end reads on an Illumina NovaSeq 6000 platform to capture both TCR α and β chains at sufficient depth, ensuring comprehensive coverage of V(D)J regions.

Bulk RNA sequencing (bulk RNA-seq)

Transcriptomic expression data, along with patient clinical and survival information for head and neck cancer (HNSC), were retrieved from the UCSC Xena database. The retrieved mRNA expression data were provided as raw count matrices.

Initial processing and quality control of single-cell RNA-seq datasets

CellRanger from 10× Genomics was used to process the raw scRNA-seq data. This pipeline performed the following steps: alignment, quantification, basic filtering, and quality control. The result was an initial gene expression matrix based on the mouse reference genome GRCh39. For downstream analysis and quality control, the R package Seurat was applied [14]. High-quality cells were selected by applying the following criteria: (i) 500–6000 genes expressed per cell, (ii) mitochondrial gene UMI counts representing less than 20% of the total UMI counts, and (iii) Genes detected in a minimum of three cells.

Downstream analysis of single-cell RNA-seq datasets

Raw counts were normalized using the NormalizeData function in Seurat. The FindVariableFeatures function was used to identify the top 2000 highly variable genes across cells. Dimensionality reduction was carried out using PCA after scaling the data with the ScaleData function. Harmony was used to correct batch effects across samples, ensuring comparability. Clustering was performed at a resolution of 0.6, and the results were visualized with UMAP embedding. Finally, based on specific marker genes, cells were classified into six major types.

CopyKat analysis of aneuploidy status in tumor cells

Using CopyKat (version 1.0.6), aneuploidy in single cells was inferred, which facilitated the identification of malignant cells through their copy number variation (CNV) profiles [15]. Raw gene expression matrices obtained from single-cell RNA-seq were used as input for CopyKat, which applies a Bayesian framework to predict chromosomal gains and losses. To differentiate between tumor and non-tumor cells, CNV scores were computed for each cell, and cells with high aneuploidy levels were classified as malignant. Default parameters were used for CopyKat.

Trajectory analysis

To identify cells with the highest differentiation potential, Cytotrace was applied to the single-cell RNA-seq dataset [16]. Cytotrace is a computational tool that predicts cellular differentiation states from gene expression profiles. The algorithm was used to score each cell based on its differentiation potential, allowing us to identify clusters with the highest potential for further trajectory analysis.

To explore differentiation trajectories within malignant subtypes, Monocle3 was applied to tumor cells [17]. The Seurat-processed tumor-cell count matrix and metadata were converted into a cell_data_set (CDS) object with gene annotation. Cells were processed using the standard Monocle3 workflow (preprocess_cds, reduce_dimension, learn_graph, and order_cells). The cluster with the highest differentiation potential identified by CytoTRACE2 was specified as the root state for pseudotime ordering. Trajectories were visualized using plot_cells and colored by malignant subtype, state, or pseudotime.

In conjunction with trajectory analysis, RNA velocity was calculated to understand the dynamic processes underlying cell differentiation [18]. RNA velocity values for each gene in tumor cells were computed using the velocyto.R package. RNA velocity vectors were subsequently projected into the UMAP embedding to visualize potential future states of cells and further clarify cellular transitions and differentiation pathways.

TCR analysis

For TCR analysis, the Cell Ranger V(D)J pipeline was employed to assemble and annotate paired TCR sequences. Clonotypes were defined based on shared CDR3 nucleotide sequences, and clonotype diversity was visualized through rarefaction analysis and clonotype frequency plots.

ScRepertoire analysis

scRepertoire was used to analyze T cell receptor (TCR) diversity and clonotype dynamics from scTCR-seq data [19]. The package was utilized to characterize TCR clonotypes, assess clonotype expansion, and investigate the distribution of TCR clones across different cell clusters. The analysis was performed to gain insights into the adaptive immune response in the TME, with particular emphasis on clonotype sharing and expansion between tumor and immune compartments.

STARTRAC analysis

For clonotype trajectory and interaction analysis, STARTRAC was applied to evaluate T cell expansion, migration [20], and transition among different cell states. STARTRAC metrics, including expansion (STARTRAC-expansion), migration (STARTRAC-migration), and transition (STARTRAC-transition), were computed to assess the dynamics of T cell populations.

CIBERSORTx Deconvolution analysis

CIBERSORTx was applied to analyze the relative fractions of TME cells in bulk RNA-seq samples using digital deconvolution [21]. The algorithm predicts the proportions of immune and stromal cell types in the TME. TCGA-HNSC data were used as input, and the reference matrix was constructed from data obtained from our own scRNA-seq HNSCC cohort. The algorithm was conducted to provide an estimation of immune cell fractions and their potential associations with clinical features. CIBERSORTx’s default parameters were used for the deconvolution, with 100 permutations to ensure the robustness of the predictions.

UCell macrophage function scoring

UCell was used to quantify macrophage functional scores derived from gene expression profiles [22]. For this analysis, gene sets representing macrophage activation and specific functional markers were selected to assess their activity in the tumor microenvironment.

Spatial transcriptomics dataset and processing

To investigate the spatial distribution of key cell subtypes, we obtained a public spatial transcriptomics (Visium) dataset of 10 HNSCC samples from the Zenodo data repository (doi: 10.5281/zenodo.14284038). ST spots were filtered based on standard quality control metrics (e.g., nUMI > 500, nGene > 300). Gene expression data for each sample was then normalized using SCTransform in Seurat (v4.3), and all 10 samples were integrated using Harmony (v0.1.1) to correct for batch effects.

Spatial deconvolution using cell2location

To estimate cell-type abundances at spot resolution, we applied cell2location using the annotated scRNA-seq dataset generated in this study as the reference atlas. The model was used to infer spot-level abundance patterns for key cell states, including PTC-like tumor cells (C6), Macro-Isg15, and CD4_Treg_Ctla4, which were then used for downstream spatial analyses.

Spatial Co-localization and proximity analysis

To quantify spatial proximity between inferred cell-type patterns, we used SpottedPy on the spot-level cell2location abundance maps. We first identified spatially enriched (“hotspot”) and sparse (“coldspot”) areas for each cell type/state using the Getis–Ord Gi* statistic as implemented in SpottedPy. We then measured proximity by computing the nearest-distance distributions between hotspots of two cell states and summarizing them by the within-sample median distance. Finally, paired Wilcoxon signed-rank tests (paired by sample) were used to compare hotspot- versus coldspot-based distances to the target hotspots, assessing whether hotspots showed significantly closer spatial proximity.

NicheNet analysis

NicheNet was applied to identify key ligands that may regulate the expression of target gene sets by leveraging prior knowledge of ligand-receptor interactions, gene regulatory networks and signaling pathways. Target gene sets were defined as upregulated genes within each treatment group (Control, LDRT, αPD-1, LDRT + αPD-1) based on significance thresholds of LogFC ≥ 0.5 and p-value ≤ 0.05 [23]. Ligands included in the analysis were those expressed in at least 5% of cells in one or more cell types. Subsequently, ligand-receptor activity analysis was conducted among Macro-Isg15, CD4⁺ Treg-Ctla4, and tumor cells to investigate intercellular communication dynamics across the treatment groups.

SCENIC analysis

SCENIC was utilized to infer gene regulatory networks (GRNs) from single-cell data, allowing the prediction of transcription factors and their regulatory elements [24]. The transcription factor annotation utilized the datasets “mm10_10kbp_up_10kbp_down_full_tx_v10_clust.genes_vs_motifs.rankings.feather” and “motifs-v10nr_clust-nr.mgi-m0.001-o0.0.tbl.” GRNBoost2 was employed to identify co-expression modules and potential transcription factors along with their target genes. RcisTarget analyzed these modules to detect enriched transcription factor binding sites, enabling the construction of transcription factor-target (TF-target) networks based on binding site enrichment. Each transcription factor and its directly associated target genes were defined as a regulon. AUCell scores were calculated to measure regulon activity at the single-cell level. Ultimately, SCENIC identified variations in transcription factor and regulon activities across different cell subclusters, providing insights into regulatory dynamics within the single-cell data.

Tissue immunofluorescence

Paraffin-embedded sections were processed for immunofluorescence as described for immunohistochemistry. Sections were incubated with primary antibodies against FOXP3 (1:200, 22228-1-AP, Proteintech), CTLA4 (1:100, YT5933, Immunoway), CD68 (1:200, 25747-1-AP, Proteintech), and ISG15 (1:500, 15981-1-AP, Proteintech). After incubation with fluorophore-conjugated secondary antibodies (G1255, Servicebio) for 1 h in the dark, sections were washed three times with PBST, counterstained with DAPI (C1002, Beyotime), mounted with an anti-fade reagent, and imaged on a Zeiss microscope. Images were analyzed using TissueGnostics.

Statistical analysis

Statistical analysis was performed using GraphPad Prism 9.5. For comparisons involving more than two groups, one-way or two-way analysis of variance (ANOVA) was applied, followed by Tukey’s multiple comparisons test as indicated in the figure legends. For two-group comparisons, unpaired two-tailed Student’s t-test was used for parametric analyses, whereas the Wilcoxon rank-sum test (Mann–Whitney U test) was used for non-parametric analyses as specified in the figure legends. Paired comparisons were performed using the two-sided Wilcoxon signed-rank test where applicable. Significance for Ro/e-based enrichment/depletion was assessed by the chi-square test, and correlation analyses were performed using Pearson correlation. Statistical significance was defined as p < 0.05 (*p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001).

Results

LDRT + αPD-1 therapy remodels the tumor microenvironment and synergistically promotes tumor regression in HNSCC

To systematically evaluate the antitumor potential of low-dose radiotherapy (LDRT) combined with αPD-1 therapy, we employed the regimen in two syngeneic HNSCC mouse models (SCC7 and MOC1), assigning animals to four groups: IgG (control), LDRT, αPD-1, or the combination. The overall study design, dosing/irradiation schedule, and the downstream single-cell profiling workflow are summarized in Fig. 1a. Briefly, LDRT was delivered at 1 Gy per fraction on days 1, 3, and 5 (three fractions total), and αPD-1 was administered at 10 mg/kg on days 1, 3, and 5 (three doses total). Longitudinal monitoring showed that, in both SCC7 and MOC1 models, LDRT + αPD-1 consistently delayed tumor growth more effectively than either monotherapy and reduced endpoint tumor burden; representative tumors and quantitative measurements were concordant (Fig. 1b–d). Throughout treatment, body weights remained stable in both models, supporting favorable tolerability (Supplementary Fig. S1a–b). Together, these cross-validating results from two independent HNSCC models demonstrate robust, reproducible, and well-tolerated synergistic antitumor activity of LDRT combined with αPD-1 blockade.

Fig. 1. Single-cell sequencing reveals that LDRT in combination with αPD-1 reprograms the tumor microenvironment in HNSCC. a Experimental design and multi-modal profiling workflow in the SCC7 syngeneic model. Tumors received Control (IgG), LDRT (1 Gy × 3 fractions on days 1, 3, and 5), αPD-1 (10 mg/kg on days 1, 3, and 5), or the combination. Tumors were collected at indicated time points for IHC, NanoString, and flow cytometry, and at endpoint (day 21) for 10x Genomics scRNA-seq and scTCR-seq. Created with BioRender. b Representative photographs of excised tumors at day 21 from SCC7 (left) and MOC1 (right) models across the four treatment groups (n = 5 mice/group). c Tumor growth curves for SCC7 (left) and MOC1 (right) under the four treatment conditions (n = 5 mice/group; mean ± SEM). Statistical comparisons over time were performed using two-way ANOVA with Tukey’s multiple comparisons test (significance shown in the plots). d Endpoint tumor burden quantification corresponding to the models in (c) (n = 5 mice/group; mean ± SEM). For each tumor model, statistical significance was assessed by one-way ANOVA followed by Tukey’s multiple comparisons test (significance is indicated in the plot). e UMAP visualization of major cell types identified from the integrated single-cell dataset. f Cell-type composition across treatment groups. The stacked bar plot shows the fraction of each major cell type per group. The heatmap displays Ro/e values (observed-to-expected ratio) indicating relative enrichment/depletion of each cell type across treatment groups, with significance assessed by a chi-square test. g Density distributions of major cell types across treatment conditions

To delineate the cellular basis of this synergy, we performed scRNA-seq and scTCR-seq on SCC7 tumors from the four treatment groups (Control, LDRT, αPD-1, and LDRT + αPD-1; n = 12 tumors total). After quality control, 81,430 cells were retained for downstream analyses. Using Harmony for batch correction and clustering, we identified six major cell types based on canonical marker genes: epithelial cells, myeloid cells, fibroblasts, endothelial cells, T lymphocytes and NK cells (Fig. 1e; Supplementary Fig. S1c, Supplementary Data 1). We next compared treatment-associated shifts in cellular composition: αPD-1 monotherapy exhibited a relatively higher epithelial fraction accompanied by reduced myeloid and T/NK fractions, whereas the LDRT + αPD-1 group showed a decreased epithelial fraction with concomitant increases in myeloid and T/NK fractions (Fig. 1f). Consistently, group-wise density distributions further supported an immune-enriched remodeling of the tumor microenvironment under the combination regimen (Fig. 1g). Collectively, these changes indicate that adding LDRT to αPD-1 reshapes the tumor microenvironment toward an immune-enriched composition, in line with the enhanced tumor growth control observed in vivo.

LDRT + αPD-1 mediates dynamic evolution of tumor cell metaprograms in HNSCC

To further explore tumor cell heterogeneity under LDRT combined with αPD-1 treatment, we first subset epithelial cells and applied CopyKAT to infer large-scale CNV profiles from scRNA-seq, thereby classifying malignant (aneuploid) tumor cells and non-malignant (diploid) cells. Based on this CNV-based classification, 29,044 aneuploid malignant tumor cells were retained for downstream analyses (Supplementary Fig. 2a-b). Dimensionality reduction and clustering analyses identified nine clusters (C1-C9) (Fig. 2a). Given the pronounced heterogeneity of tumor cells, we applied non-negative matrix factorization (NMF) to extract interpretable tumor metaprograms that capture recurrent transcriptional programs across cells. After filtering out low-quality related metaprograms, we identified four tumor cell metaprograms: Oxphos or stress (MP1), Cell cycle (MP2), Interferon (MP3), and RNA Transcription Regulation (MP4) (Fig. 2b).

Fig. 2LDRT-driven tumor-cell state transition and immunogenic modulation underlying synergy with αPD-1. a UMAP plot showing the distribution of different tumor cell clusters (C1-C9). b GeneNMF-based metaprogram analysis identified four tumor cell metaprograms. c Expression levels of each metaprogram across different tumor cell clusters. d UMAP density plots of tumor cell clusters in different treatment groups, showing cell density distribution. e UMAP of tumor-cell states defined by metaprogram-associated features, annotated as SIRTC, ITC, and PTC. f Monocle3 trajectory analysis revealing differentiation directions of tumor cells. g Differential expression between C9 (SIRTC) and C6 (PTC-enriched state) and corresponding GO/KEGG enrichment analyses. h Kaplan–Meier survival curves illustrating the association between the signature scores of C6 and C9 and overall survival in TCGA-HNSC patients (n = 500). Group differences were assessed by the log-rank test. i Flow cytometry histograms showing surface PD-L1 and MHC class I molecules (H-2Kb and H-2Db) in MOC1 cells with or without LDRT. j Quantification of mean fluorescence intensity (MFI) fold changes for PD-L1, H-2Kb, and H-2Db in MOC1 cells. Statistical significance was determined using two-way ANOVA followed by Sidak’s multiple comparisons test (n = 6 independent samples per group; data are presented as mean ± SEM; ***p < 0.001). k CCK8 assays showing cell growth curves of SCC7 and MOC1 cells with or without LDRT. Statistical significance was determined using two-way ANOVA followed by Sidak’s multiple comparisons test (n = 6 independent samples per group; data are presented as mean ± SEM; ***p < 0.001). l Colony formation assays showing clonogenic survival of SCC7 and MOC1 cells with or without LDRT. Statistical significance was determined using two-way ANOVA followed by Sidak’s multiple comparisons test (n = 6 independent samples per group; data are presented as mean ± SEM; p < 0.001). m Quantification of γH2AX foci per nucleus at indicated time points after irradiation in SCC7 and MOC1 cells. n Representative immunofluorescence images of γH2AX foci (with DAPI) at indicated time points (0, 0.5, 2 and 4 h) after irradiation in SCC7 and MOC1 cells. Scale bars = 50 μm. Data are presented as mean ± SEM from at least 3 independent biological replicates. For each replicate, 20 random fields were captured and analyzed for γH2AX foci quantification. Statistical significance was determined using two-way ANOVA followed by Sidak’s multiple comparisons test ( p < 0.001 relative to the 0 h control group within each cell line)

Subsequently, we calculated scores for each cluster using UCell based on the top genes from each metaprogram. Results showed that C9 had high scores for both Oxphos or stress (MP1) and Interferon (MP3), and was highly enriched in the LDRT and LDRT + αPD-1 groups, suggesting that C9 represents a stress- and interferon-responsive tumor state enriched under LDRT-containing conditions. Therefore, we named this cluster Stress-Immune Responsive Tumor Cells (SIRTC). Interestingly, we found that clusters C6 and C7 had high Cell cycle (MP2) scores and were highly enriched in the αPD-1 group, suggesting that they may exhibit reduced sensitivity to αPD-1 treatment, maintaining high proliferative potential even under αPD-1 treatment (Fig. 2c-d, Supplementary Fig. 2c). We named these clusters Proliferative Tumor Cells (PTC). The remaining cells were named Intermediate Tumor Cells (ITC) (Fig. 2e, Supplementary Fig. 2c, Supplementary Data 1). Accordingly, SIRTC was defined as cluster C9, PTC as clusters C6–C7, and ITC as clusters C1–C5 and C8 for downstream analyses.

We also used Monocle3 and RNA velocity to analyze the dynamic evolution of tumor cells under LDRT + αPD-1 treatment. First, using CytoTRACE2, we evaluated the differentiation potential of tumor cells. Results showed that PTC had the highest differentiation potential, while SIRTC had the lowest, suggesting that PTC may represent a relatively less differentiated, progenitor-like transcriptional state (Supplementary Fig. 2d-e). Therefore, we used cells with the highest differentiation potential as the starting point and found that tumor cells differentiated in two directions. Branch 1 extended from C7 to C9, indicating differentiation towards a state of high oxidative stress and interferon-related gene expression. Branch 2 extended from C7 toward C6, consistent with a trajectory toward a more proliferative, cell-cycle–dominant PTC-like state (Fig. 2f, Supplementary Fig. 2f).

To define molecular features accompanying this transition, we focused on the two trajectory endpoints and compared C9 (SIRTC) with C6 (PTC-enriched proliferative state). Differential expression highlighted interferon-stimulated and immune-associated genes in C9, including Gbp2, Gbp2b, Iigp1, and Zbp1, together with upregulation of Cd274 (PD-L1) (Fig. 2g, Supplementary Data 2). GO and KEGG enrichment analyses indicated that the C9/SIRTC state was enriched for antigen processing and presentation pathways together with p53 signaling, autophagy, and DNA damage–induced apoptotic programs, collectively reflecting a stress-responsive cellular state, whereas the C6/PTC state was dominated by nuclear division, cell-cycle progression, and DNA replication programs (Fig. 2g). At the hallmark level, PTC was enriched for proliferative programs (E2F targets, G2M checkpoint, MYC targets, mitotic spindle, and DNA repair), while SIRTC was enriched for immune stress programs (TGF-β signaling, IL2–STAT5 signaling, inflammatory response, IFN-α/IFN-γ responses, and TNFα signaling via NFκB), consistent with their proliferative versus stress/IFN-responsive transcriptional states (Supplementary Fig. 2 g, Supplementary Data 4). In TCGA-HNSC, CIBERSORTx-based deconvolution suggested that higher C6-like signature scores were significantly associated with worse survival, whereas higher C9-like scores showed a trend toward better survival but did not reach statistical significance (Fig. 2h).

We performed SCENIC transcription factor analysis to understand the transcriptional regulatory networks of these clusters. Results showed that in the PTC tumor cell subgroup, transcription factors E2F1, MYBL2, ZFP64, and RAD21 exhibited significant activity, suggesting that these factors may promote the proliferative and immune-evading characteristics of this subgroup. E2F1, a key regulator of the cell cycle, had its stability enhanced by m5C modification, which may drive tumor cell proliferation [25]. MYBL2 exhibited high activity in the tumor microenvironment, promoting tumor cell proliferation while enhancing immune evasion through the induction of immunosuppressive macrophage infiltration, thereby increasing resistance to anti-PD-1 therapy [26]. Additionally, ZFP64 drove immunosuppressive remodeling of the tumor microenvironment via the PKCα/ZFP64/CSF1 axis, further exacerbating resistance to immunotherapy [27]. RAD21, through interaction with the YAP/TEAD4 complex, inhibited interferon signaling, promoting immune evasion and possibly enhancing resistance to anti-PD-1 therapy [28]. In contrast, the SIRTC tumor cell subgroup showed significant activity of transcription factors IRF9 and GABPA, suggesting that LDRT may induce SIRTC cells into a stress and immune-responsive state. IRF9, as a key regulator of the interferon signaling pathway, was upregulated and activated interferon-responsive genes, which may limit the proliferative potential of SIRTC [29]. Additionally, GABPA regulated the proliferative properties of SIRTC by activating TGFBR2 expression, and this pathway activation may further inhibit SIRTC proliferation [30]. These results indicate that, compared with the PTC subgroup, SIRTC exhibits a specific stress response under LDRT treatment, which may suppress its proliferative activity (Supplementary Fig. 2 h).

Finally, we performed in vitro assays to evaluate the direct impact of LDRT on tumor-cell immunogenic features and radiation responses. Flow cytometry in MOC1 cells showed that LDRT increased surface PD-L1 and MHC class I molecules (H-2Kb and H-2Db) (Fig. 2i–j). Consistently, LDRT suppressed cell growth in both SCC7 and MOC1 cells in CCK8 assays and reduced clonogenic survival (Fig. 2k–l). Moreover, LDRT induced a robust DNA damage response, as evidenced by increased γH2AX foci and representative immunofluorescence images across time points (Fig. 2m–n). Collectively, these data indicate that LDRT (alone or in combination with αPD-1) is associated with a shift of tumor cells from a proliferative, cell-cycle–dominated PTC state toward a stress/interferon-responsive SIRTC state accompanied by enhanced antigen-presentation features and increased DNA-damage responses.

LDRT + αPD-1 increases tumor-reactive CD8⁺ T-cell representation and reduces CD4⁺ Treg-Ctla4 frequency in HNSCC

To comprehensively explore the role of T cells during LDRT + αPD-1 treatment, we isolated lymphocytes (T & NK cells) and performed re-clustering, resulting in 11 distinct transcriptional states of T and NK cell subpopulations: naive T cells (Tn-Tcf7), CD8⁺ effector memory T cells (CD8⁺ Tem-Ccl4), CD8⁺ tissue-resident memory T cells (CD8⁺ Trm-Klrd1), CD8⁺ effector T cells (CD8⁺ Teff-Gzmk), CD8⁺ exhausted T cells (CD8⁺ Tex-Havcr2), CD8⁺ proliferative T cells (CD8⁺ Mki67), CD4⁺ follicular helper T cells (CD4⁺ Tfh-Cxcl13), CD4⁺ regulatory T cells (CD4⁺ Treg-Ctla4), CD4⁺ proliferative T cells (CD4⁺ Mki67), natural killer cells (NK), and proliferative natural killer cells (NK-Mki67) (Fig. 3a-b, Supplementary Data 1).

Fig. 3LDRT combined with αPD-1 modulates functional states of tumor-infiltrating T cells. a UMAP analysis showing 11 reclustered T and NK cell clusters. b Heatmap showing the expression (z-score) of signature genes in different T and NK cell clusters. c Heatmap showing the distribution of T and NK cell clusters based on the ratio of observed cell counts to expected cell counts (Ro/e), calculated using a chi-square test. The top bar plot illustrates the cell composition. d Cytokine secretion levels (e.g., IFN-γ, TNF-α, and IL-12) across treatment groups, measured by ELISA. Data are presented as mean ± SEM. Statistical significance determined by one-way ANOVA with Tukey’s multiple comparisons test. * p < 0.05, ** p < 0.01, **p < 0.001, **** p < 0.0001. n = 3 biological replicates per group. e Representative flow cytometry plots of CD8 (CD8-PE, y-axis) versus CD4 (CD4-APC, x-axis) among gated [live singlets → CD45⁺CD3⁺] cells across treatment groups, with quantification of [CD8⁺ and/or CD4⁺] frequencies shown on the right (n = 4). f Quantitative analysis of CD45 + CD3+ CD8 + T cells proportions in SCC7 tumors across treatment groups (n = 5 per group). Data are presented as mean ± SEM. Statistical significance determined by one-way ANOVA with Tukey’s multiple comparisons test. * p < 0.05, ** p < 0.01, *** p < 0.001. g Quantitative analysis of CD3 + CD8+ IFN-γ CTLs and CD3 + CD8+ TNF-α CTLs proportions in SCC7 tumors across treatment groups (n = 5 per group). Data are presented as mean ± SEM. Statistical significance determined by one-way ANOVA with Tukey’s multiple comparisons test. * p < 0.05, ** p < 0.01, *** p < 0.001. h Violin plots showing the expression levels of Ctla4 in CD4 + cells across different treatment groups. Pairwise group comparisons were performed using the Wilcoxon rank-sum test (Mann–Whitney U). ( p < 0.05, ** p < 0.01, *** p < 0.001, **** p < 0.0001; ns, not significant)

We then compared the relative frequencies of these subpopulations across treatment groups. The frequencies of CD8⁺ Tem-Ccl4, CD8⁺ Trm-Klrd1, and CD8⁺ Teff-Gzmk—subsets implicated in antitumor immunity—were reduced in the αPD-1 group but rebounded in the LDRT + αPD-1 group. In contrast, the frequency of CD4⁺ Treg-Ctla4, a subset associated with immunosuppression, was highest in the αPD-1 group and decreased in the LDRT + αPD-1 group (Fig. 3c, Supplementary Fig. 3a).

Compared with LDRT or αPD-1 alone, ELISA results showed that the levels of IFN-γ, TNF-α, and IL-12 were elevated in the LDRT + αPD-1 group (Fig. 3d). To validate these findings in an independent tumor model, we conducted parallel analyses in MOC1 tumors. Consistent with SCC7 results, combination therapy significantly elevated IFN-γ, TNF-α, and IL-12 levels in the tumor microenvironment as measured by ELISA (Supplementary Fig. 3c). These cytokine changes are consistent with an overall shift toward a more immune-active microenvironment. Flow cytometry further showed a marked increase in the frequency of intratumoral CD8⁺ T cells (Fig. 3e-f, Supplementary Fig. 3b), whereas the frequency of CD4⁺ T cells remained relatively unchanged (Supplementary Fig. S3d). Notably, the proportions of CD3⁺CD8⁺ IFN-γ^+^ cytotoxic T cells (CTLs) and CD3⁺CD8⁺ TNF-α^+^ CTLs were also significantly increased in SCC7 tumors treated with LDRT + αPD-1 (Fig. 3g). Together, these findings suggest that LDRT + αPD-1 increases the abundance and effector activation of intratumoral CD8⁺ T cells.

αPD-1 resistance may be related to the presence of other immune checkpoints in T cells, such as CTLA-4, TIGIT, and LAG-3 [31]. By comparing immune checkpoint expression across different treatment groups, we found that Ctla4 expression in CD4⁺ T cells was upregulated in the αPD-1 group compared with the control group, which may be associated with the reduced efficacy of αPD-1 treatment. Interestingly, only Ctla4 expression in the LDRT and LDRT + αPD-1 groups was significantly lower than in the αPD-1 group (Fig. 3h, Supplementary Fig. 3e-f), suggesting that LDRT may inhibit Ctla4 expression in CD4⁺ T cells, thereby enhancing the antitumor effect of αPD-1.

In summary, these results indicate that adding LDRT to αPD-1 reshapes the intratumoral immune landscape by increasing tumor-reactive CD8⁺ T-cell abundance and effector cytokine production while reducing the frequency of CD4⁺ Treg-Ctla4, consistent with the superior tumor growth control observed under combination therapy.

LDRT + αPD-1 reshapes T cell immune landscape as revealed by scTCR-seq

To trace the dynamic evolution of T-cell clonoyptes, we utilized scRepertoire to analyze the clonotypes of CD8⁺ and CD4⁺ T cells during LDRT + αPD-1 treatment. We observed distinct distribution of T cell clonotypes on the UMAP, with an increased clone size in CD8⁺ T cell subpopulations. Notably, the total clone size of large clones was highest in CD8⁺ Tem-Ccl4, CD8⁺ Trm-Klrd1, and CD8⁺ Teff-Gzmk, which were highly enriched in the LDRT + αPD-1 group. This further confirmed that these three CD8⁺ T cell subpopulations underwent clonal expansion under tumor antigen stimulation, playing a key role in antitumor immunity, and are critical effector cells in the LDRT + αPD-1 treatment. In contrast, the clone size of CD4⁺ T cells decreased (Fig. 4a-b, Supplementary Data 5).

Fig. 4. Analysis of clonal expansion and state transitions in CD8 + and CD4 + T cells. a UMAP plot showing clonal expansion of T cell clusters, with different colors representing different clone sizes. b Distribution of clone sizes in different T cell clusters, showing the frequency of clonal expansion. c Heatmap of clonal overlap between different T cell clusters, revealing shared clonality among these populations. d STARTRAC analysis calculating expansion (expa), migration (migr), and transition (tran) indices for different T cell clusters to quantify their dynamic behaviors in the TME. e Scatter plot of proliferation frequency versus expansion index, showing the relationship between proliferative activity and clonal expansion in different T cell clusters. f Trajectory map derived from RNA velocity analysis, showing the relationships between T cell clusters during state transitions and highlighting the developmental trajectories among different clusters

To further investigate the lineage relationships among T cell phenotypes, we quantified the fraction of clonotypes shared between T cell clusters. High clonotype overlaps were found among CD8⁺ T cell subpopulations, as well as some clonotype overlaps with the Tn-Tcf7 subpopulation, suggesting that some CD8⁺ T cell subpopulations may have differentiated from Tn-Tcf7. Among CD4⁺ T cell subpopulations, CD4⁺ Tfh-Cxcl13 showed clonotype overlaps with CD8⁺ T cell subpopulations and Tn-Tcf7, indicating potential cooperative roles in the immune response, consistent with reports that low-dose radiotherapy can promote TLS formation, in which CXCL13⁺ Tfh-like programs contribute to organized antitumor immunity [32]. Notably, CD4⁺ Treg-Ctla4 and CD4⁺ Mki67 exhibited the lowest clonotype overlaps with Tn-Tcf7, suggesting that CD4⁺ Treg-Ctla4 and CD4⁺ Mki67 may not have differentiated from Tn-Tcf7 in the tumor microenvironment but rather migrated from the periphery into the tumor microenvironment (Fig. 4c).

To further elucidate T cell clonal expansion potential, tissue migration, and state changes, we performed STARTRAC analysis, integrating scRNA-seq and scTCR-seq data. The results showed that CD8⁺ Tem-Ccl4, CD8⁺ Trm-Klrd1, CD8⁺ Teff-Gzmk, and CD8⁺ Mki67 exhibited clonal expansion and sustained proliferation, with a high degree of clonotype sharing and low migration levels, further supporting their differentiation from Tn-Tcf7 in response to tumor antigens, demonstrating high activity and antigen reactivity. RNA velocity results also confirmed this. In contrast, CD8⁺ Tex-Havcr2 showed high clonal sharing but low migration, clonal expansion, and proliferation, suggesting it may be in a state of functional exhaustion with weak reactivity to tumor antigens. As previously mentioned, CD4⁺ Treg-Ctla4 and CD4⁺ Mki67 exhibited low clonal expansion, low proliferation, and low clonal sharing, but high migration levels, further confirming that they migrated from the periphery rather than differentiating from Tn-Tcf7 in the tumor microenvironment (Fig. 4d-f).

CD4⁺ Treg-Ctla4 frequency correlates with Macro-Isg15

We performed re-clustering and annotation of the remaining myeloid cells and fibroblasts. The myeloid cells were annotated into monocytes, eight macrophage subpopulations (Macro-Spp1, Macro-Mki67, Macro-Lgmn, Macro-Isg15, Macro-Il1b, Macro-Fn1, Macro-Fcer1g, Macro-Arg1), and three dendritic cell (DC) subpopulations (cDC1-Clec9a, cDC2-Cd1a, cDC3-Lamp3) (Fig. 5a, Supplementary Data 1). Fibroblasts were annotated into three clusters (CAF-Mki67, CAF-Col5a5, CAF-Cd74). Endothelial cells, due to their small numbers, were not further clustered.

Fig. 5TME cell subset composition dynamics and monocyte-macrophage state polarization under LDRT combined with αPD-1 treatment. a UMAP plots illustrating the distribution of tumor cells, T cells, NK cells, fibroblasts, endothelial cells, and myeloid cells. b CIBERSORTx was used to predict relative fractions of all cell clusters in TCGA-HNSC, followed by Pearson correlation analysis (n = 500). c Heatmap displaying proportional changes in monocyte-macrophage clusters based on the ratio of observed cell counts to expected cell counts (Ro/e), calculated using a chi-square test, across treatment groups. The top bar plot illustrates the relative enrichment of these clusters. d UCell scores showing M2 markers, antigen presentation capacity, pro-angiogenic ability, and T cell attraction potential of different monocyte-macrophage clusters. e Representative flow cytometry contour plots showing CD206 expression (x-axis, CD206-PE; y-axis, SSC-A) in F4/80⁺CD11b⁺ myeloid cells from SCC7 tumors across the indicated treatment groups (Control, LDRT, αPD-1, and LDRT + αPD-1). Percentages denote the frequency of CD206⁺ cells within the gated F4/80⁺CD11b⁺ population. The bar plot (right) summarizes the proportion of F4/80⁺CD11b⁺CD206⁺ (M2-like) cells (%), with each dot representing one biological replicate (mean ± [SEM/SD]; n = [per group]). Statistical significance was determined by one-way ANOVA with Tukey’s multiple comparisons test. * p < 0.05, ** p < 0.01, *** p < 0.001. f Representative flow cytometry contour plots showing CD86 expression (x-axis, CD86–Brilliant Violet 605; y-axis, SSC-A) in F4/80⁺CD11b⁺ myeloid cells from SCC7 tumors across the indicated treatment groups. Percentages denote the frequency of CD86⁺ cells within the gated F4/80⁺CD11b⁺ population. The bar plot (right) summarizes the proportion of F4/80⁺CD11b⁺CD86⁺ (M1-like) cells (%), with each dot representing one biological replicate (mean ± SEM; n = 4). Statistical significance was determined by one-way ANOVA with Tukey’s multiple comparisons test * p < 0.05, ** p < 0.01, *** p < 0.001. g Trajectory map of monocyte-macrophage clusters based on RNA velocity analysis, revealing their developmental trajectories and potential functions in the TME, showing differentiation fate relationships

To investigate cellular contexts associated with CD4⁺ Treg-Ctla4 abundance—and to identify myeloid programs that co-vary with this Treg signal in larger patient cohorts—we used our scRNA-seq data to derive reference signatures and performed CIBERSORTx deconvolution in TCGA-HNSC, followed by Pearson correlation analysis. The results revealed that the estimated fraction of CD4⁺ Treg-Ctla4 showed the strongest correlation with the estimated fraction of the myeloid subpopulations Macro-Isg15 and Macro-Fcer1g (Fig. 5b, Supplementary Fig. 4a-b). Furthermore, the expression of ISG15 and FCER1G is positively correlated with the estimated fraction of CD4⁺ Treg-Ctla4 (Supplementary Fig. 4c-d), suggesting a potential link between the migration of CD4⁺ Treg-Ctla4 and these macrophage subpopulations.

Monocyte-macrophages play a critical role as mediators in the tumor TME, particularly in modulating T cell immune responses against tumors. To gain a deeper understanding of the dynamic changes in macrophage subpopulations during LDRT + αPD-1 treatment, we used the Ucell R package to characterize the functional features of different macrophage subpopulations through macrophage-related pathways and compared their relative frequencies across different treatment groups. The results showed that the Macro-Lgmn subpopulation had the highest M2 macrophage score, suggesting that it represents M2 macrophages with an immunosuppressive phenotype, consistent with previous studies [33, 34]. Intergroup comparisons revealed that Macro-Lgmn frequency was higher in the LDRT group than in the control group, whereas it significantly decreased in both the αPD-1 and LDRT + αPD-1 groups. In contrast, the Macro-Il1b subpopulation exhibited the lowest M2 macrophage score and the highest antigen presentation score, suggesting that it may represent M1 macrophages with antitumor immune functions. Compared to the control group, Macro-Il1b frequency did not increase in the LDRT group but significantly increased in both the αPD-1 and LDRT + αPD-1 groups (Fig. 5c-d, Supplementary Data 6).

Flow cytometry analysis further validated these findings. Compared to other groups, the LDRT + αPD-1 group exhibited the highest M1 macrophage frequency and the lowest M2 macrophage frequency (Fig. 5e-f). To assess whether these immunomodulatory effects on macrophage polarization were consistent across tumor models, we extended the analysis to MOC1 tumors. Similar to the results observed in SCC7 tumors, the proportion of F4/80⁺CD11b⁺CD86⁺ M1-like macrophages was significantly increased, while F4/80⁺CD11b⁺CD206⁺ M2-like macrophages were markedly decreased in the LDRT + αPD-1 group compared with the monotherapy and control groups (Supplementary Fig. 4e-g). These findings indicate that the combination therapy modulates macrophage phenotypes by decreasing M2 macrophage frequency and increasing M1 macrophages, thereby enhancing the antitumor immune response.

Notably, we observed that the Macro-Isg15 displayed chemokine/immune-interaction–associated transcriptional features (as reflected by the T cell recruitment–related score). Compared to the control group, its relative frequency was higher in the αPD-1 group but decreased in the LDRT and LDRT + αPD-1 groups, further supporting an association in which Macro-Isg15 enrichment co-occurs with higher CD4⁺ Treg-Ctla4 abundance under αPD-1 monotherapy, while LDRT-containing regimens are associated with attenuation of this macrophage–Treg coupling. Meanwhile, the Macro-Fcer1g subpopulation had a low T cell recruitment score. Intergroup analysis revealed that while Macro-Fcer1g frequency was higher in the αPD-1 group compared to the control group, no significant difference was observed between the αPD-1 and LDRT + αPD-1 groups. This suggests that LDRT treatment may not influence its frequency (Fig. 5c-d, Supplementary Fig. 4 h). RNA velocity analysis also indicated that Macro-Isg15 and Macro-Fcer1g originated from monocytes but followed two distinct differentiation directions, representing their cellular heterogeneity (Fig. 5g).

In summary, our data suggest that LDRT + αPD-1 not only modulates macrophage polarization but is also associated with attenuation of an ISG15-high macrophage–Treg axis observed under αPD-1 monotherapy, consistent with reduced immunosuppressive features in the TME.

Spatial transcriptomics and paired mIF support a Macro-Isg15–Treg niche in PTC-like regions

Our CIBERSORTx-based deconvolution analysis in TCGA-HNSC revealed a positive correlation between the estimated fraction of CD4_Treg_Ctla4 cells and the Macro-Isg15 macrophage subset (Fig. 5b). To provide in situ spatial support for this association, we analyzed a public human HNSCC spatial transcriptomics cohort (n = 10). Using our scRNA-seq data as reference, we performed cell2location-based deconvolution to infer per-spot abundances of PTC-like tumor cells (C6), Macro-Isg15, and CD4_Treg_Ctla4. Representative spatial abundance maps showed that C6-enriched regions tend to co-occur with Macro-Isg15- and CD4_Treg_Ctla4-enriched regions within the same tissue context (Fig. 6a–c).

Fig. 6. Spatial transcriptomics analysis of human HNSCC validates the spatial proximity of PTC-like tumor cells, Macro-Isg15, and Treg cells. a–c Representative spatial maps showing cell2location-inferred spot-level abundances of PTC-like tumor cells (C6), Macro-Isg15 macrophages, and CD4_Treg_Ctla4 cells in a human HNSCC spatial transcriptomics section. d–f Quantification of spatial proximity using SpottedPy based on spot-level cell2location abundance matrices. Getis–Ord Gi* (local spatial autocorrelation; neighborhood defined by a KNN/neighbor-number parameter) was applied to classify spots as hotspots, coldspots, or not significant for each cell state. For each sample, nearest-distance distributions to target hotspots were computed and summarized as within-sample medians, followed by paired comparisons across samples (n = 10). Violin plots show the distribution of within-sample median distances, with paired points indicating per-sample medians. d Median nearest distance from C6 hotspot versus C6 coldspot spots to the nearest Macro-Isg15 hotspot spots. * p < 0.05, ** p < 0.01, *** p < 0.001. e Median nearest distance from C6 hotspot versus C6 coldspot spots to the nearest CD4_Treg_Ctla4 hotspot spots. * p < 0.05, ** p < 0.01, *** p < 0.001. f Median nearest distance from Macro-Isg15 hotspot versus Macro-Isg15 coldspot spots to the nearest CD4_Treg_Ctla4 hotspot spots. * p < 0.05, ** p < 0.01, *** p < 0.001. g Representative multiplex immunofluorescence (mIF) images from an independent paired HNSCC cohort (pre-treatment biopsy vs. post-treatment specimen; n = 5). DAPI (blue) was used as nuclear counterstain; CD68 (magenta), ISG15 (yellow), FOXP3 (green), and CTLA4 (red) are shown as indicated. Scale bars, 50 μm (main images) and 10 μm (insets). h–i Paired quantification of (h) FOXP3/CTLA4 double-positive Treg cells and (i) CD68/ISG15 double-positive macrophages per unit area in matched pre- versus post-treatment tissues (paired dots represent matched samples from the same patient; n = 5). Two-sided Wilcoxon signed-rank test; (j) Representative multi-channel mIF images illustrating in situ proximity between CD68/ISG15 double-positive macrophages and FOXP3/CTLA4 double-positive Treg cells in paired specimens

We then quantified this relationship by comparing within-sample nearest-distance distributions between hotspot and coldspot regions. Across samples, C6 hotspot spots were significantly closer to Macro-Isg15 hotspots than C6 coldspot spots (Fig. 6d, p = 0.0117), and likewise were closer to CD4_Treg_Ctla4 hotspots (Fig. 6e, p = 0.0020). Notably, Macro-Isg15 hotspots themselves were also significantly closer to CD4_Treg_Ctla4 hotspots than Macro-Isg15 coldspots (Fig. 6f, p = 0.0488). Together, these spatial proximity patterns provide in situ evidence that PTC-like tumor regions are spatially coupled with Macro-Isg15 enrichment and adjacent Treg hotspots, supporting a localized immunosuppressive niche in human HNSCC.

To independently evaluate these observations in human tissues, we analyzed an independent paired clinical cohort of five HNSCC patients with pre-treatment biopsies and post-treatment specimens obtained under a neoadjuvant low-dose radiotherapy–based combination regimen (two 21-day cycles; IMRT 1 Gy per fraction on days 1, 2, 8, and 15 of each cycle, 8 fractions/8 Gy total; combined with tislelizumab, albumin-bound paclitaxel, and cisplatin). In multiplex immunofluorescence, Treg cells were defined by FOXP3 and CTLA4 co-staining, and ISG15 macrophages were defined by CD68 and ISG15 co-staining. Representative images and paired quantification indicated that both Treg infiltration density (FOXP3/CTLA4 double-positive cells) and ISG15 macrophage infiltration density (CD68/ISG15 double-positive cells) were reduced after treatment (Fig. 6g–i, Supplementary Data 7).

Building on the density-based readouts readouts, we further examined the spatial relationship between CD68/ISG15 double-positive macrophages and FOXP3/CTLA4 double-positive Treg cells in situ. Representative multi-channel images suggested a weakened spatial association, and nearest-neighbor distance–based quantification showed an increased separation between CD68/ISG15⁺ macrophages and FOXP3/CTLA4⁺ Treg cells in post-treatment samples (Fig. 6j).

LDRT + αPD-1 remodels the communication network between CD4⁺ Treg-Ctla4, Macro-Isg15, and malignant cells

Given that Macro-Isg15 enrichment co-occurs with higher CD4⁺ Treg-Ctla4 abundance in our data, and that tumor cells can shape macrophage states through cytokine/chemokine signaling [35–37], we used NicheNet to infer intercellular signaling links among Macro-Isg15, Macro-Fcer1g, CD4⁺ Treg-Ctla4, and tumor cells. First, we compared the communication differences between Macro-Isg15 and Macro-Fcer1g across the four treatment groups. The results showed that, in the αPD-1 group, the ligand Mif was upregulated in tumor cells, and NicheNet prioritized Mif–Cxcr4 and Mif–Cd74 as putative axes associated with Macro-Isg15 enrichment/activation in the αPD-1 group, whereas these inferred interactions were weaker in LDRT-containing groups (Fig. 7a, Supplementary Fig. 5a). The heatmap also indicated that the expressions of Mif, Cxcr4, and Cd74 were cell-specific (Fig. 7b). Furthermore, we found that these interactions regulated downstream target genes, affecting the functions of Macro-Isg15, and GO enrichment analysis revealed that upregulated target genes were associated with monocyte proliferation (Fig. 7c).

Fig. 7. Regulation of ligand-receptor interactions among different cell types in the TME by LDRT combined with αPD-1 treatment. a Heatmaps of NicheNet analysis showing regulatory patterns between tumor cells and Macro-Isg15 among different treatment groups. b Heatmaps showing changes in ligand and receptor expression across different cell types, based on z-score normalization. c Representative GO biological process (BP) enrichment analysis illustrating the functional enrichment of predicted target genes expressed in Macro-Isg15 and tumor cells. d Heatmaps of NicheNet analysis showing regulatory patterns between Macro-Isg15 and CD4⁺ Treg-Ctla4 among different treatment groups. e Heatmaps of ligand and receptor expression based on z-score normalization, showing changes across different cell types in each treatment group. f Representative GO biological process (BP) enrichment analysis illustrating the functional enrichment of predicted target genes expressed in Macro-Isg15 and CD4⁺ Treg-Ctla4

Similarly, using Macro-Fcer1g as a control group, we performed Nichenet analysis between Macro-Isg15 and CD4⁺ Treg-Ctla4. The results showed that, compared to other groups, Cxcl14 in Macro-Isg15 was upregulated in the αPD-1 group, and NicheNet highlighted Cxcl14–Cxcr4 as a putative chemotactic axis associated with a higher CD4⁺ Treg-Ctla4 signal. In parallel, Cd86 was upregulated in both Macro-Isg15 and Macro-Fcer1g in the αPD-1 group, and NicheNet al.so prioritized the Cd86–Ctla4 interaction, consistent with established co-inhibitory signaling between myeloid cells and T cells [38–40]. Notably, both inferred axes were attenuated in the LDRT + αPD-1 group (Fig. 7d, Supplementary Fig. 5b). The heatmap results further demonstrated that Cd86 was highly expressed in myeloid cells, whereas Ctla4 was predominantly expressed in lymphoid cells (T & NK), providing additional evidence of the regulatory influence of myeloid cells on lymphoid cells (Fig. 7e). Enrichment analysis of downstream target genes showed an upregulation in leukocyte chemotaxis-related pathways (Fig. 7f), further supporting the possibility that Macro-Isg15 recruits CD4⁺ Treg-Ctla4.

In conclusion, integrating ligand–receptor inference with multi-layer evidence suggests that, under αPD-1 pressure, tumor-derived MIF–CXCR4/CD74 signaling aligns with enrichment of ISG15-high macrophage states and is accompanied by reinforced CXCL14–CXCR4 and CD86–CTLA4 interactions, consistent with a Treg-permissive immunosuppressive niche. Notably, LDRT-containing regimens are associated with attenuation of this axis, in agreement with the observed cell-state shifts and the in situ spatial/protein readouts. Collectively, these results narrow the mechanistic explanation to a focused set of axes and provide a clear direction for how LDRT reconditions the HNSCC microenvironment to improve PD-1 blockade efficacy.

Discussion

Anti-PD-1 therapy in HNSCC has limited efficacy, with only around 15% of patients responding, primarily because of the immunosuppressive environment within the TME [5]. The TME is enriched with Tregs and TAMs, which reduce the effectiveness of anti-PD-1 therapy through various mechanisms [6, 7]. This study investigates the potential of LDRT to modulate the TME and enhance immune responses, particularly in combination with αPD-1 therapy. Using scRNA-seq and scTCR-seq, we identified the mechanisms by which LDRT enhances immunotherapy efficacy by modulating tumor cell states and reprogramming immune cells. Our findings provide new insights and strategies to enhance the therapeutic response of HNSCC patients to αPD-1 therapy.

Our study demonstrates that LDRT combined with αPD-1 alters tumor cell states primarily by regulating their evolution during treatment. Using scRNA-seq and NMF, we identified distinct tumor cell populations and observed that LDRT + αPD-1 induces a stress and immune-responsive state in tumor cells. We identified a distinct subset of SIRTC that was markedly overrepresented in the LDRT and combination therapy groups, characterized by high expression of interferon-related genes. This likely results from LDRT-induced stress, driving these cells into a more immune-active state. In contrast, proliferative tumor cells (PTC) were relatively enriched under αPD-1 monotherapy and displayed elevated cell-cycle activity, a pattern that may reflect tumor adaptation under immune pressure. Together, these observations support the view that LDRT can shift malignant programs away from a proliferative-dominant state. To provide orthogonal support beyond transcriptomic state inference, our in vitro assays indicate that LDRT directly increases tumor-cell immunogenic features and activates canonical radiation responses. In MOC1 cells, LDRT upregulated surface PD-L1 together with MHC class I molecules (H-2Kb/H-2Db), consistent with enhanced antigen presentation potential and a tumor cell state that may be more “visible” to cytotoxic T cells when PD-1 blockade is present. In parallel, LDRT suppressed tumor-cell growth and clonogenic survival and induced a robust DNA damage response (γH2AX). Collectively, these functional readouts align with the in vivo observation that LDRT-containing regimens enrich stress/IFN-responsive tumor programs and disfavor cell-cycle–dominant proliferative states.

Cell trajectory analysis revealed different evolutionary paths for SIRTC and PTC, with SIRTC transitioning towards an immune-responsive state, while PTC maintained high proliferative potential. Consistent with this observation, PTC was enriched in genes related to the cell cycle (Hallmarks G2 Checkpoint, E2F-Target, Mitotic spindle) and DNA repair, which may indicate higher intrinsic radiosensitivity [41]. Transcription factor analysis further supported these findings, showing significant activation of transcription factors such as IRF9 and GABPA in SIRTC, indicating that these cells may limit their proliferative capacity and enhance immune responses by activating interferon and stress-related pathways [29, 30]. Conversely, transcription factors such as E2F1 and MYBL2 were active in the PTC population, likely promoting their proliferative characteristics and contributing to treatment resistance [25–30].

In addition, our study shows that LDRT markedly enhances the effectiveness of anti-PD-1 therapy in HNSCC by promoting the clonal expansion of antitumor CD8⁺ T cell subsets, while being accompanied by a reduction in CD4⁺ Treg-Ctla4 abundance, consistent with decreased recruitment/migration inferred from TCR-based analyses. In particular, we noted a marked increase in various CD8⁺ T cell subsets, including effector memory T cells (CD8⁺ Tem-Ccl4), tissue-resident memory T cells (CD8⁺ Trm-Klrd1), and effector T cells (CD8⁺ Teff-Gzmk) in the combination group. These cell types play crucial roles in antitumor immunity [42, 43].

Importantly, these preclinical observations are clinically consistent with our Phase II NeoRTPC02 trial in HNSCC [13], where an LDRT-containing regimen (1 Gy per fraction) combined with αPD-1 therapy was associated with increased intratumoral CD8⁺ T-cell infiltration (IHC, p = 0.0028) and enrichment of IL7R⁺ effector-memory–like CD8⁺ subsets (scRNA-seq). While the clinical data primarily reflect immune infiltration and the murine data highlight clonal expansion dynamics, the convergent immune remodeling patterns across species support the translational relevance of LDRT-induced enhancement of CD8⁺ antitumor responses.

In contrast, immunosuppressive Treg (CD4⁺ Treg-Ctla4) showed a significant decrease in frequencyand reduced CTLA-4 expression in the LDRT + αPD-1 group. This finding is important. As confirmed by our scTCR-seq and STARTRAC analyses, CD4⁺ Treg-Ctla4 show limited clonal overlap and high migration, suggesting they mainly originate from peripheral migration. LDRT may reduce the recruitment of these immunosuppressive cells (potentially via the Macro-Isg15 axis discussed later), allowing effector T cells to fully exert their functions within the tumor.

Interestingly, our study shows that LDRT + αPD-1 not only alters T cell landscape but also significantly impacts specific subsets of myeloid cells [9–12], particularly the interaction between the Macro-Isg15 subset and CD4⁺ Treg-Ctla4. We re-clustered and annotated myeloid cells and fibroblasts, identifying multiple subsets. Among them, the Macro-Isg15 subset exhibited transcriptional features consistent with chemokine-associated immune-interacting programs. Our pearson correlation analysis identified a significant positive correlation between CD4⁺ Treg-Ctla4 abundance and Macro-Isg15 abundance, supporting a close association between Macro-Isg15 states and CD4⁺ Treg-Ctla4 accumulation, consistent with a microenvironment permissive for Treg enrichment. Moreover, recent studies have highlighted the role of IFN-regulated tumor-associated macrophages (IFN-TAMs), including those expressing ISG15, which resemble immunosuppressive macrophages that suppress immune responses through mechanisms such as regulatory T cell (Treg) recruitment [44]. This further supports the notion that Macro-Isg15 states are associated with an immunosuppressive microenvironment, consistent with co-occurrence of Treg-Ctla4 features. Crucially, our human in situ analyses provide complementary spatial support for the macrophage–Treg association implied by the single-cell data. In a public HNSCC spatial transcriptomics cohort (n = 10), PTC-like tumor regions (C6) showed significant spatial proximity to Macro-Isg15-enriched hotspots, and Macro-Isg15 hotspots were also proximal to CD4_Treg_Ctla4 hotspots, supporting a localized niche in which proliferative tumor programs co-occur with ISG15-high macrophage states and Treg features.

Extending these observations to protein-level human tissue readouts, multiplex immunofluorescence in an independent paired neoadjuvant cohort (pre-treatment biopsy vs. post-treatment specimen, n = 5) treated with an LDRT-containing combination regimen showed reduced densities of FOXP3/CTLA4 double-positive Treg cells and CD68/ISG15 double-positive macrophages after treatment. Moreover, nearest-neighbor distance–based quantification suggested attenuated in situ proximity between ISG15 macrophages and Treg cells post-treatment. Together, the spatial transcriptomics and paired mIF data strengthen the clinical relevance of an ISG15-macrophage–Treg niche and support the notion that LDRT-containing regimens can disrupt this immunosuppressive microenvironmental coupling. In contrast, the Macro-Il1b subset exhibited higher antigen-presenting capabilities and significantly increased frequency in the LDRT + αPD-1 group, suggesting that combination therapy may also enhance antitumor immunity by promoting M1 macrophage frequency [44]. Together, these findings indicate that LDRT may reshape macrophage states under immune checkpoint blockade, with potential implications beyond the irradiated lesion. Notably, previous studies have shown that low-dose radiotherapy combined with anti–PD-L1 can induce an abscopal effect; moreover, selective depletion of tumor-associated macrophages (TAMs) in distant, non-irradiated tumors abolishes this abscopal response, suggesting that this systemic therapeutic benefit is largely TAM-dependent and involves TAM remodeling [45].

To further clarify how LDRT modulates the tumor–macrophage–Treg circuit observed in our data, we used Nichenet analysis to explore the signaling interactions among macrophages, CD4⁺ Treg-Ctla4, and tumor cells. The results suggested that, under αPD-1 treatment, tumor-derived MIF signaling may contribute to the enrichment and/or activation of Macro-Isg15 states through the Mif–Cxcr4 and Mif–Cd74 axes. However, these interactions were weakened in the LDRT group, suggesting that LDRT might mitigate this immunosuppressive effect.

Further inference indicated that Macro-Isg15 states under αPD-1 treatment were associated with higher Cxcl14 and Cd86 signaling, which could plausibly support Treg-Ctla4 recruitment/activation through Cxcl14–Cxcr4 and Cd86–Ctla4 interactions; these signals were attenuated in the LDRT + αPD-1 group. Heatmap analysis showed that Cd86 was primarily expressed in myeloid cells, while Ctla4 was mainly expressed in lymphoid cells, further supporting the regulatory role of myeloid cells on lymphoid cells. Enrichment analysis of downstream target genes also supported the idea that Macro-Isg15 states are associated with an immunosuppressive TME, consistent with increased Treg-Ctla4 accumulation.

In summary, our data support a model in which LDRT + αPD-1 is associated with attenuation of an ISG15-high macrophage–Treg circuit, thereby relieving immunosuppressive features within the TME. These findings provide new insights into how LDRT can improve immunotherapy efficacy by modulating interactions between macrophage subsets and Treg cells. Our results suggest that enrichment of Macro-Isg15 states is linked to an immunosuppressive milieu under αPD-1 therapy, in part through co-occurrence with Treg-Ctla4 features, whereas LDRT is associated with attenuation of this macrophage–Treg coupling. Moreover, LDRT’s influence on tumor cell states is crucial for boosting immune responses. Collectively, our multi-layer evidence—from tumor-state remodeling and TCR-informed T-cell dynamics in mice, to orthogonal in vitro functional validation and in situ human spatial and mIF readouts—converges on a coherent framework in which LDRT reconditions the HNSCC microenvironment toward a more therapy-permissive state under PD-1 blockade. This framework highlights two actionable levers: shifting malignant programs away from cell-cycle–dominant proliferative states and weakening an ISG15-high macrophage–Treg suppressive niche that co-occurs with reduced effector function. In this context, the immune remodeling patterns observed in NeoRTPC02 [13]. are not merely consistent but strongly supportive of the translational relevance of this mechanism-guided combination strategy. Moving forward, functional perturbation of the inferred axes and prospective validation in patient cohorts will be essential to lock down causality and define biomarker-driven stratification for rational deployment of LDRT plus anti–PD-1 in HNSCC.

Our study has limitations. First, although our in vitro tumor-cell experiments provide direct evidence for tumor-intrinsic effects of LDRT (including upregulation of PD-L1 and MHC-I [H-2Kb/H-2Db], induction of a γH2AX-marked DNA damage response, and suppression of proliferation and clonogenic survival; Fig. 2i–n), we did not perform compartment-separation validations (e.g., matched-dose irradiation of purified macrophages/stromal cells, or in vivo approaches such as bone-marrow chimeras or compartment-restricted genetic models). Therefore, we cannot yet decisively quantify the dominant compartment in which LDRT acts or the relative contribution of tumor-intrinsic versus immune/non-tumor effects. Second, while the NeoRTPC02 clinical study provides translational support for our findings, the key mechanistic chain in this work still relies primarily on murine models and computational inference, and the clinical regimen differs from the animal setting in dose/fractionation and combination components; thus, extrapolation should be made cautiously. Third, the spatial transcriptomics analysis was based on public human datasets and supports spatial associations of cell-state programs/signatures rather than direct cell–cell contact or causality. Fourth, key inferred axes (e.g., MIF–CXCR4/CD74, CXCL14–CXCR4, and CD86–CTLA4) were derived from ligand–receptor inference and correlative analyses; functional perturbation experiments will be needed to establish causality. Finally, other components of the TME not deeply profiled here may also contribute to combination efficacy and warrant further investigation.

Conclusion

In summary, our study demonstrates that low-dose radiotherapy acts as a potent sensitizer for anti-PD-1 therapy in head and neck squamous cell carcinoma by orchestrating a comprehensive remodeling of both the tumor cell states and the immune microenvironment. Through integrated single-cell transcriptomic and TCR sequencing, we have elucidated that LDRT drives a critical shift in tumor metaprograms, transitioning malignant cells from a highly proliferative state (PTC) to an “immunogenic” stress-responsive state (SIRTC) characterized by enhanced antigen presentation and interferon signaling. We observed that an Isg15^high^ macrophage state was positively correlated with CD4^+^ Treg-Ctla4 frequency, particularly under αPD-1 monotherapy. Ligand–receptor inference further identified a Macro-Isg15-associated CXCL14–CXCR4 axis that was enriched during αPD-1 monotherapy but significantly attenuated in LDRT-containing regimens. These findings suggest that LDRT helps overcome resistance by reducing the chemokine-driven recruitment and accumulation of immunosuppressive Treg cells.

Supplementary Information

Supplementary Material 1.

Supplementary Material 2.

Supplementary Material 3.

Supplementary Material 4.

Supplementary Material 5.

Supplementary Material 6.

Supplementary Material 7.