Establishment and characterization of the treatment cohort
Between 2018 and 2019, we collected primary tumor tissue samples from 83 patients with LARC prior to treatment and performed bulk RNA-seq. All patients received radiotherapy-based regimens, and their tumor regression grades (TRGs) were documented. The clinical characteristics of these patients are summarized in Supplementary Table S1. Among the 49 patients with known microsatellite status, only one had dMMR/MSI-H, while the remainder had pMMR/MSS. Compared with patients with an intermediate treatment response (TRG 2), those with the poorest response (TRG 3) and those with the best response (TRG 0/1) displayed similarly elevated baseline T cell exhaustion scores and MHC class II antigen presentation scores (Fig. 1a), although deconvolution analysis did not indicate exclusive T cell enrichment in either the TRG 0/1 or TRG 3 samples (Supplementary Fig. S1d). Given that the tumor cells in the TRG 0/1 group were largely eradicated by radiotherapy, we focused on understanding why patients in the TRG 3 group had the poorest therapeutic response, which was worse than that in the TRG 2 group wherein minimal to no immune inflammation was detected, despite exhibiting localized immune inflammation prior to treatment. Similarly, bulk RNA-seq was performed on tumor tissue samples from an independent cohort of 114 treatment-naïve patients with LARC32, and among the three pMMR/MSS subtypes classified on the basis of immune inflammation, the IG2 subtype — with intermediate immune inflammation — exhibited a worse trend in prognosis than the IG1 subtype, although the degree of immune inflammation was lowest in patients with the IG1 subtype (Supplementary Fig. S1a, b). At baseline, increased T cell exhaustion and MHC class II antigen presentation are associated with improved responses to ICI therapies33,34,35,36. This raises the question of whether patients in the TRG 3 group may benefit from immunotherapy. However, since ICI-based immunotherapy is not currently a standard first-line treatment for patients with MSS CRC, the aim of our study was to explore the potential benefits of combining radiotherapy with immunotherapy to improve treatment efficacy and to identify potential predictive biomarkers for this therapeutic strategy.
Fig. 1: Single-cell atlas of tumor and peripheral cellular dynamics during the treatment of patients with LARC.
a Boxplots displaying GSVA scores for the T-cell exhaustion gene set (comprising PDCD1, HAVCR2, CTLA4, TIGIT, and CXCL13) and the MHC-II antigen presentation gene set (GO:0019886) across different TRG groups in the CRC bulk RNA-seq cohort (n = 82 after quality control) at baseline. P-values were calculated using the Wilcoxon rank-sum test. b Schematic diagram showing the study design, treatment timeline, sample collection points, and data types analyzed. c UMAP plot depicting the distribution of major cell types across 238,680 cells from all LARC tumor samples. d Heatmap showing the distributions of major cell types between the CR and PR groups at baseline (Pre), post-radiotherapy (In), and post-immunotherapy (Post). The Ro/e index (see Methods for calculation) is represented using the symbols +++ (> 2), ++ (> 1.5), + (> 1), and +/- ( > 0.5). e UMAP plot showing the distribution of major cell types across 243,609 PBMCs from all peripheral blood samples. f Heatmap showing the enrichment of each co-occurrence cluster in the pre-treatment samples.
Herein, we monitored the complete treatment course of 22 patients diagnosed with LARC who received care at Peking University People’s Hospital between January 2023 and August 2024. Among these patients, 20 underwent long-term preoperative radiotherapy combined with anti-PD-1 immunotherapy, while the remaining two received chemotherapy alone without immunotherapy. All tumors were confirmed to be mismatch repair-proficient and microsatellite-stable (pMMR/MSS) on the basis of immunohistochemical (IHC) staining for mismatch repair (MMR) proteins. All patients exhibited tumor regression following treatment. Patients with a tumor regression grade (TRG) of 0 were classified into the complete response (CR) group (n = 7), while the remaining patients were assigned to the partial response (PR) group (n = 15). No treatment-related or immune-related adverse events or perioperative or 30-day postoperative surgical complications were observed. Sphincter-preserving surgery was achieved in 18/20 patients (90%). Detailed clinical information for all patients is summarized in Supplementary Table S2.
Of the 20 patients who received the combined radiotherapy and anti-PD-1 immunotherapy regimen, we obtained consent from 18 patients and collected a range of samples. These included tumor tissue and peripheral blood samples collected before treatment (LARC_Pre and PBMC_Pre, respectively), tumor tissue samples after radiotherapy (LARC_In), and tumor tissue and peripheral blood samples following immunotherapy (LARC_Post and PBMC_Post, respectively). In total, we collected 85 paired fresh samples, which were subsequently subjected to scRNA-seq and scTCR-seq. Our primary analysis focuses on these samples to address two key questions: (1) What pre-treatment factors contributed to the differing clinical responses between patients who achieved a CR and those who achieved a PR, and did these differences diminish or disappear as treatment progressed? (2) Why did patients in the PR group retain residual tumor cells after radiotherapy and immunotherapy, and what are the mechanisms underlying their resistance?
To assess whether our findings in patients with pMMR/MSS CRC extended to immune-responsive tumors, we analyzed two independent cohorts of patients with dMMR/MSI-H CRC who received anti-PD-1 monotherapy. The first cohort consisted of 22 patients for whom processed scRNA-seq data were available in GSE23658131. The second cohort consisted of 19 patients for whom processed scRNA-seq data were available in GSE20550637. Moreover, we performed whole-exome sequencing (WES) on tumor tissue samples from three patients (P8, P9, and P10) at three time points: Before treatment (Pre), after radiotherapy (In), and after immunotherapy (Post). For each of these samples, the WES results from the corresponding PBMC_Post samples were used as controls. We also analyzed a published WES dataset comprising 56 paired LARC tumor tissue samples collected before and after chemoradiotherapy for validation (Fig. 1b).
As mentioned earlier, for the two patients in this study who received chemotherapy alone, we collected tumor and peripheral blood samples pre- and post-chemotherapy for scRNA-seq. Through further integration with scRNA-seq data from 29 rectal cancer patients who received neoadjuvant chemotherapy38, we explored the following question: (3) Were there differences in tumor microenvironment remodeling after radiotherapy and chemotherapy administration, and which approach was more suitable for combination with subsequent immunotherapy?
Single-cell atlas of tumor and peripheral cellular dynamics during treatment of patients with LARC
We first constructed a single-cell atlas of tumor and peripheral cellular dynamics during treatment of patients with LARC. After obtaining the scRNA-seq data for the 85 samples mentioned above, we performed rigorous quality control (see Methods for details) and excluded 5 sequencing samples with low quality. After removing low-quality cells and doublets, we retained a total of 482,289 cells, comprising 238,680 from LARC tissue and 243,609 PBMCs. The cells from LARC tissue were categorized into four major groups: lymphocytes (B cells, plasma cells, CD4+ T cells, CD8+ T cells, γδ T cells and NK cells), myeloid cells (mast cells, plasmacytoid dendritic cells (pDC), conventional dendritic cells (cDC), monocytes, macrophages, and neutrophils), stromal cells (fibroblasts, smooth muscle cells, and endothelial cells), and epithelial cells (normal epithelial and tumor cells), along with a group of proliferating cells (cycling cells) (Fig. 1c; Supplementary Fig. S1c). The effects of different treatment regimens on TME remodeling in LARC tissue were sufficiently pronounced to overshadow the differences in the TME between patients who achieved a CR and patients who achieved a PR at various stages of treatment. Specifically, patients in the CR and PR groups exhibited predominant enrichment of stromal and tumor cells in the TME prior to treatment. Radiotherapy reduced the abundances of stromal and tumor cells and increased the proportion of innate immune cells, particularly myeloid cells, within the TME. This finding is consistent with previous studies showing that radiotherapy can activate the innate immune system by inducing immunogenic cell death (ICD) and releasing tumor-associated antigens39,40. Following immunotherapy, the TME in the CR and PR groups shifted toward a state with greater immune infiltration. There were notable increases in the abundances of T lymphocytes and plasma cells, indicating a transition from a “cold” tumor to a “hot” tumor, a hallmark of immune activation and increased adaptive immune responses (Fig. 1d).
Given that changes in TME cell composition after radiotherapy and immunotherapy were similar between the CR and PR groups, we sought to determine whether pre-existing differences could explain the variation in treatment response. Our analysis revealed that, prior to treatment, there was no significant difference in the proportion of immune cells in LARC tissue between the two groups, except for a slightly greater proportion of neutrophils in the CR group (P = 0.092) (Supplementary Figs. S1e, S2a). This finding suggests that differences in treatment response may have arisen from differences in molecular expression that affect the function of similar cell populations in the two groups. Furthermore, we categorized all the PBMCs into 14 major cell types and 59 subtypes (Fig. 1e; Supplementary Fig. S2b). On the basis of the correlations among cell proportions before treatment, we grouped these 59 subtypes into 10 co-occurrence clusters (Supplementary Fig. S2c). Cluster 7 and Cluster 8 were enriched in the PR group before treatment and were composed primarily of activated lymphocytes, including CD4T_MHCII and tTreg_MHCII cell subtypes expressing high levels of MHC-II molecules41, the CD8T_HAVCR2 cell subtype expressing the exhaustion marker HAVCR2, and plasma cells (Fig. 1f). These findings suggest that pre-existing differences existed between the PR and CR groups prior to treatment, with patients in the PR group potentially exhibiting a state of sustained immune activation, which may be associated with systemic chronic inflammation.
Lower HLA-DQA2 expression in dendritic cells at baseline is associated with better treatment responses in patients with LARC
Since myeloid cells, such as macrophages, monocytes, and DCs, are among the first innate immune populations to accumulate in the TME after radiotherapy, we extracted these cells from patients at all LARC stages and subdivided them into 16 subtypes to study their functions (Fig. 2a; Supplementary Fig. S3a). Similarly, monocytes and DCs from PBMCs were subdivided into 11 subtypes for detailed analysis (Supplementary Fig. S3c). However, the proportions of these subtypes were largely consistent between the CR and PR groups before treatment (P > 0.05), suggesting that no single subtype predominantly contributed to the effects of combined radiotherapy and immunotherapy. We therefore compared the differences in gene expression between these subtypes in the CR and PR groups at each treatment stage. In addition to identifying differentially expressed genes (DEGs) with Padj < 0.05, we incorporated diff_pct (the difference in the proportion of cells expressing a gene between the CR and PR groups) to better evaluate DEGs and mitigate the influence of individual patient variation on avg_log2FC. Compared with those in the CR group, the avg_log2FC and diff_pct of HLA-DQA2 expression were greater before treatment and after radiotherapy in the PR group (Fig. 2b, c).
Fig. 2: Lower baseline HLA-DQA2 expression in dendritic cells is associated with better treatment responses in patients with LARC.
a UMAP projection of macrophages, monocytes, and cDCs from all LARC samples. Each dot represents a single cell and is colored according to its cell type. b, c Volcano plots showing significantly upregulated or downregulated genes (P < 0.05) in the CR group compared with the PR group before treatment (b) and after radiotherapy (c), along with corresponding avg_log2FC and diff_pct values (diff_pct: the difference in the proportion of cells expressing each gene between the two groups). d Diffusion map illustrating the developmental trajectory of DCs, with the colors of each DC subset consistent with those in a. e Ridge plot showing HLA-DQA2 expression in cDCs across post-radiotherapy LARC samples, with sample groupings and total cDC counts indicated on the right. f Violin plot showing the expression of MHC-II molecules, including HLA-DQA2, in cDCs from LARC tissue (blue) and peripheral blood (red) at all treatment stages for each patient in this dataset. CR patients are shown on the left, and PR patients are shown on the right. g Boxplot showing HLA-DQA2 expression in DCs from individual patients in the CR and PR groups; P values were calculated using the Wilcoxon rank-sum test. h Diffusion map illustrating the CytoTRACE scores of individual DCs, with higher scores indicating later positions along the developmental trajectory and arrows denoting the developmental direction. i Density plot showing the distribution of DC subtypes in the low- and high-HLA-DQA2 groups across treatment stages, with coordinates matched to those in d, h. j Boxplots comparing the proportions of the cDC_LAMP3 subtype between the low- and high-HLA-DQA2 groups (left) and between the CR and PR groups (right) at baseline, with each dot representing one patient; P values were calculated using the Wilcoxon rank-sum test. k Boxplot showing changes in cDC_LAMP3 subtype proportions before and after radiotherapy in low-HLA-DQA2 patients, with patient colors matched to those in j; P values were calculated using the Wilcoxon signed-rank test. l Enrichment analysis showing the top five upregulated pathways in the cDC_LAMP3 (red) and cDC_CD1C (blue) subtypes.
Unlike HLA-DQA1, which is polymorphic, HLA-DQA2 is considered a non-polymorphic gene and is expressed primarily in antigen-presenting cells such as DCs42. Therefore, we extracted all DCs (Fig. 2d) from LARC tissues of all patients following radiotherapy and measured HLA-DQA2 expression to exclude the possibility that differences in HLA-DQA2 expression between the CR and PR groups were driven by a single patient. The results revealed that patients P4, P5, P10, and P14 (with a TRG of 0) and patient P9 (with a TRG of 1) exhibited low HLA-DQA2 expression. This pattern was not attributable to a small number of captured DCs, as these patients contributed a greater number of DCs (Fig. 2e). We also observed that patients P1 and P15, both of which had high HLA-DQA2 expression, achieved CR. These findings are consistent with our previous observations from the bulk RNA-seq cohort, where a group of patients with higher pre-treatment immune activity achieved a TRG of 0 (Fig. 1a; Supplementary Fig. S1b). We were then concerned with why a greater number of patients with high HLA-DQA2 expression, compared with those with low expression, do not derive maximal benefit from combined radiotherapy and immunotherapy. Considering that treatment-induced inflammation can sometimes increase the expression of MHC-II molecules, we next assessed the expression of HLA-DQA2 and other MHC-II molecules in DCs from LARC tissues and PBMCs at all treatment stages for each patient. LARC tissues from patients with low HLA-DQA2 expression consistently exhibited this pattern before treatment, after radiotherapy, and after immunotherapy, as did their PBMCs before and after treatment (Fig. 2f). This suggested that low HLA-DQA2 expression was unlikely the result of sequencing dropout in a single sample, as the likelihood of dropout co-occurring across all five samples from the same patient is extremely low. Instead, this likely reflected an intrinsic patient characteristic that may have contributed to a greater likelihood of achieving a clinical CR.
Given the limited sample size of our cohort, we validated our findings using data from two additional Chinese cohorts of dMMR/MSI-H CRC patients who received anti-PD-1 monotherapy. scRNA-seq data from CRC tissue and peripheral blood samples collected before, during, and after treatment were reported for the cohort in the study by Zhang et al.31. In contrast, Deng et al. provided scRNA-seq data from CRC tissue samples collected before and after treatment37. Patients with low HLA-DQA2 expression were more likely to experience favorable outcomes (P = 0.0015) (Fig. 2g; Supplementary Fig. S3d, e). Moreover, data from the Zhang et al. dataset confirmed that low HLA-DQA2 expression in DCs from these patients was not attributable to low cell counts and did not vary by tissue type (CRC vs. peripheral blood) or sampling time point (Supplementary Fig. S3f). This consistency across cohorts and sampling conditions reinforces the reliability of our findings and suggests that low HLA-DQA2 expression may serve as a robust biomarker for predicting favorable responses to immunotherapy.
Although the function of HLA-DQA2 in DCs is not fully understood, recent studies have highlighted its potential role in predicting the outcomes of individuals after SARS-CoV-2 infection. Specifically, increased HLA-DQA2 expression is associated with milder symptoms or even asymptomatic cases after SARS-CoV-2 infection43,44. We performed pseudotemporal trajectory analysis on all the DCs (Fig. 2d) and found that cDC_LAMP3, a mature DC subtype, was located at the end of the developmental trajectory (Fig. 2h). This subset was also enriched in pre-treatment samples from patients with high HLA-DQA2 expression (P = 0.066) (Fig. 2i, j), including those in the CR group (P1 and P15). This explains why there was no significant difference in the proportion of the cDC_LAMP3 subset between the CR and PR groups at baseline. In contrast, a greater proportion of patients with low HLA-DQA2 expression had a higher proportion of the immature DC subtype cDC_CD1C before treatment (P = 0.05) (Supplementary Fig. S3b). However, these cells rapidly matured and became activated after radiotherapy and then enriched after ICI therapy, with functional characteristics such as increased positive regulation of T cell activation, positive regulation of leukocyte cell–cell adhesion, and differentiation. In contrast, the DCs from patients with high HLA-DQA2 expression formed a TME dominated by the cDC_CD1C subtype after combined radiotherapy and immunotherapy and played a greater role in neutrophil activation (Fig. 2i, k, l). This dynamic change in DC behavior may be because although higher HLA-DQA2 expression increases immune responses to viral infections, it may contribute to chronic or dysregulated immune activation in the context of cancer, making it difficult to benefit from immunotherapy and leading to distinct outcomes in these two types of patients.
These findings may also help explain our earlier observation that patients in the PR group paradoxically exhibited stronger baseline immune inflammation in their peripheral blood, which could be driven by an increased antigen presentation capacity associated with high HLA-DQA2 expression (Fig. 1f), and this activated MHC-II antigen presentation pathway may reflect a strong positive correlation with Treg activity among all CD4+ T cells (Supplementary Fig. S3i). Integration of multiple CRC scRNA-seq cohorts revealed that HLA-DQA2 is more highly expressed in DCs within the TME of patients with metastatic CRC, while it is expressed at lower levels in the peripheral blood of healthy individuals (Supplementary Fig. S3g). In contrast, HLA-DQA1 and HLA-DRB1 are highly expressed in the DCs of normal adjacent cancer tissues and are consistently highly expressed in nearly all patients (Supplementary Fig. S3h). These findings address the first question: Low baseline HLA-DQA2 expression is associated with better treatment response, and this hallmark does not disappear as treatment progresses.
Persistent expansion of pre-existing tumor-infiltrating CD8+ T cell clones leads to a suboptimal treatment response
We next systematically delineated the dynamic clonal evolution of CD8+ T cells during sequential radiotherapy and immunotherapy in patients with low HLA-DQA2 expression (CR group) versus those with high expression (PR group). CD8+ T cells from patients with all LARC stages were classified into nine subtypes. Among them, the CD8T_PDCD1 subtype highly expressed exhaustion markers such as HAVCR2, PDCD1, LAG3, TIGIT, CTLA4, and CXCL13 and was therefore defined as a terminally exhausted CD8+ T cell (Tex) subtype (Fig. 3a, d). In contrast, the CD8T_MHCII subtype exhibited high expression of activation markers, including MHC-II molecules and various granzymes, and was thus defined as an activated CD8+ T cell subtype. We reconstructed TCR sequences from all the samples using TCR sequencing data and found that the clonal diversity and richness of the naïve subtypes, CD8T_ANXA1 and CD8T_IL7R, were greatest (Fig. 3b). In contrast, activated CD8+ T cells (CD8T_MHCII) displayed the lowest diversity, which was even lower than that of terminally exhausted CD8+ T cells (CD8T_PDCD1). This may be attributed to CD8T_MHCII being the most significantly expanded subtype in both groups following immunotherapy, representing a highly activated population associated with the immune response in this treatment context. Furthermore, we observed that the Tex subtype was enriched primarily in baseline samples from the high-HLA-DQA2 expression group, whereas the activated CD8T_MHCII subtype following immunotherapy was enriched predominantly in the group with low HLA-DQA2 expression (Fig. 3c). This pattern was also detected at the individual patient level, where HLA-DQA2 expression in DCs inversely correlated with the proportion of CD8T_MHCII cells (R = -0.55; P = 0.044) and positively correlated with that of CD8T_TNFSF9 cells (R = 0.69; P = 0.007) (Supplementary Fig. S4a, b).
Fig. 3: Persistent expansion of pre-existing tumor-infiltrating CD8+ T cell clones drives a suboptimal treatment response.
a UMAP projection of CD8+ T cells from all LARC samples. Each dot represents a single cell and is colored according to its cell type. b Box plots showing TCR clonotype diversity among CD8+ T cell subtypes. The colors correspond to the subsets shown in a. c Heatmap showing the distributions of CD8+ T cells in the low- and high-HLA-DQA2 groups at baseline (Pre), post-radiotherapy (In), and post-immunotherapy (Post), with higher Ro/e index values (see Methods) indicating stronger enrichment tendencies. d Bubble plot showing the expression of marker genes in CD8+ T cell subtypes. e UMAP projection of terminally exhausted CD8+ T (Tex) cells, Tex-shared CD8+ T cells (Tex-S cells; sharing TCR clonotypes with Tex cells), and bystander CD8+ T cells. f Pie charts showing the distribution of CD8+ T cell subtypes among Tex-S cells in the low-HLA-DQA2 (top) and high-HLA-DQA2 (bottom) groups. g Heatmap showing the distributions of Tex, Tex-S, and bystander CD8+ T cells in the low- and high-HLA-DQA2 groups at baseline (Pre), post-radiotherapy (In), and post-immunotherapy (Post). h Heatmap showing the distributions of TCR clonotypes of different sizes across the six groups, with “>50” indicating clonotypes shared by at least 50 CD8+ T cells. i Lollipop plot showing k-mer enrichment in tumor-reactive CD8+ T cells (protective effect) and bystander CD8+ T cells (risk effect). j Heatmap showing the Bray–Curtis (BC) similarity index (upper left) and Morisita index (lower right) of tumor-reactive CD8+ T cells in LARC tissues across the Pre, In, and Post groups with low and high HLA-DQA2 expression. k Boxplot showing paired differences in the Bray–Curtis similarity index of tumor-reactive CD8+ T cells before and after treatment in patients with low and high HLA-DQA2 expression; P values were calculated using the Wilcoxon rank-sum test.
Considering that bystander CD8+ T cells can also undergo expansion following immunotherapy, we focused on tumor-reactive CD8+ T cells, which are thought to share TCR clonotypes with the Tex subset25,45. We defined these non-terminally exhausted tumor-reactive CD8+ T cells as Tex-Shared (Tex-S) cells, which included precursor exhausted CD8+ T (Tpex) cells as well as some early-stage CD8+ T cells (Fig. 3e). Together, Tex and Tex-S cells constitute the tumor-reactive CD8+ T cell population. Although the composition of Tex-S differed between the two groups—primarily because of the enrichment of CD8T_MHCII in the group with low HLA-DQA2 expression and CD8T_TNFSF9 in the group with high HLA-DQA2 expression — CD8T_MHCII remained the main contributor to Tex-S in both groups and likely represents Tpex (Fig. 3f). As expected, the proportion of Tex-S cells increased significantly following immunotherapy in the low-HLA-DQA2 expression group (P < 0.0001), whereas in the high-expression group, the proportion of Tex-S cells remained relatively stable across the baseline, post-radiotherapy, and post-immunotherapy stages, with the abundance of Tex cells decreasing after treatment (Fig. 3g; Supplementary Fig. S4c). Similarly, hyperexpanded clonotypes were also enriched in the low-expression group after therapy (Fig. 3h). Collectively, these results indicate that clonal expansion is more pronounced in patients with low HLA-DQA2 expression.
We next performed an in-depth analysis of paired scTCR-seq data. TCR k-mer analysis was conducted to examine the CDR3 sequence features across different CD8+ T cell populations. Tumor-reactive CD8+ T cells were enriched in short k-mers, such as KL, RA, and EK, which contain hydrophobic and charged residues, reflecting flexible motifs likely involved in tumor antigen recognition. In contrast, bystander CD8+ T cells were enriched in longer k-mers, including NTEAFF, NTEA, and NTEAF, which often contain polar or aromatic residues, suggesting that these motifs are structural or public and lack direct tumor specificity (Fig. 3i) (see Methods).
Additionally, tumor-reactive CD8+ T cells in the group with low HLA-DQA2 expression exhibited lower TCR sequence similarity between pre-treatment and post-immunotherapy time points, as measured by both the Bray‒Curtis similarity index (0.02) and the Morisita index (0.01) (Fig. 3j) (see Methods for the calculations). In contrast, in the group with high HLA-DQA2 expression, tumor-reactive CD8+ T cells after immunotherapy were more likely to be derived from tumor-reactive cells already present in the pre-treatment TME, likely reflecting their prior local expansion and proportional dominance before therapy. This pattern was validated at the individual patient level (P = 0.067) (Fig. 3k). Notably, P1, a CR patient with high HLA-DQA2 expression, displayed a relatively similar TCR composition before and after treatment, whereas P9, a PR patient with low HLA-DQA2 expression, showed markedly different TCR repertoires before and after therapy. These observations suggest that stratification on the basis of HLA-DQA2 expression may better capture the underlying biology, given that numerous factors influencing whether a patient achieves CR or PR cannot be fully accounted for. Interestingly, we observed a similar phenomenon in an independent cohort of triple-negative breast cancer (TNBC) patients treated with combined radiotherapy and immunotherapy46 (Supplementary Fig. S4d).
Tumor-reactive CD8+ T cells derived from PBMCs exhibit an enhanced effector function
As mentioned previously, in the group with high HLA-DQA2 expression, tumor-reactive CD8+ T cells that expanded following immunotherapy originated primarily from the pre-treatment TME, whereas in the low-expression group, post-treatment tumor-reactive CD8+ T cells showed low similarity to pre-treatment TCR clonotypes. We therefore hypothesized that these new TCR clonotypes in the low-expression group were derived primarily from PBMCs. CD8+ T cells from PBMCs were classified into nine subtypes, and paired scTCR-seq analysis allowed us to identify potential tumor-reactive Tex-S CD8+ T cells in the peripheral blood whose composition was similar between the two groups (Fig. 4a, b; Supplementary Fig. S4e). As expected, in the low-HLA-DQA2 expression group, tumor-reactive CD8+ T cells in the post-treatment TME displayed high similarity to Tex-S cells in PBMCs (Bray‒Curtis similarity index = 0.51), even though the baseline similarity was essentially negligible (Bray‒Curtis similarity index = 0.00) (Fig. 4c). This likely reflects the scarcity of exhausted CD8+ T cells and the lack of pre-existing tumor-reactive CD8+ T cells in the TME of the low-expression group before treatment. In contrast, Tex-S cells among PBMCs from the high-expression group remained relatively stable before and after therapy (Fig. 4d).
Fig. 4: Tumor-reactive CD8+ T cells derived from PBMCs exhibit enhanced effector function.
a UMAP projection of CD8+ T cells from all the PBMC samples. b UMAP projection of Tex-shared CD8+ T cells (Tex-S cells; sharing the same TCR clonotypes as Tex cells), and bystander CD8+ T cells from all the PBMC samples. c Heatmap showing the Bray–Curtis similarity indices between tumor-reactive CD8+ T cells from LARC tissues and PBMCs in the low- and high-HLA-DQA2 groups before and after treatment. d Bar plot showing changes in the proportions of tumor-reactive CD8+ T cells in peripheral blood across the Pre and Post groups with low and high HLA-DQA2 expression. e Boxplots showing the expansion (expa), Gini, migration (migr), and transition (tran) indices of Tex, Tex-S, and bystander CD8+ T cells from LARC tissues and PBMCs, with significance assessed for PBMC-derived Tex-S cells and bystander CD8+ T cells. **P < 0.01, ***P < 0.001, ****P < 0.0001, as determined as determined by the Wilcoxon rank-sum test. f Boxplots comparing the migration (migr) and transition (tran) indices of tumor-reactive and bystander CD8+ T cells from LARC tissues and PBMCs between the low- and high-HLA-DQA2 groups. *P < 0.05, as determined as determined by the Wilcoxon rank-sum test. g Pie charts showing the origins of newly emerged tumor-reactive TCR clonotypes in the low-HLA-DQA2 (top) and high-HLA-DQA2 (bottom) groups, with PBMC-derived clonotypes shown in red and Unknown (other tissue sources or not captured by sequencing) clonotypes shown in blue. h Boxplots comparing Le and Gascuel substitution distances (see Methods) between Unknown TCR clones and PBMC clones in the low- and high-HLA-DQA2 groups. ***P < 0.001, as determined as determined by the Wilcoxon rank-sum test. Boxplots and cumulative density plots comparing exhaustion scores (i) and effector scores (j) among Preexisting, Unknown, and PBMC-derived tumor-reactive CD8+ T cells. ****P < 0.0001, as determined as determined by the Wilcoxon rank-sum test and FDR correction.
Next, we applied STARTRAC47 to monitor the cross-tissue migration capacity of tumor-reactive CD8+ T cells before and after therapy in the two groups. Compared with bystander CD8+ T cells, these tumor-reactive Tex-S cells exhibited greater expansion, migration, and transition indices (P < 0.01) in PBMCs (Fig. 4e). Notably, at the individual patient level, Tex-S cells in the low-HLA-DQA2 expression group displayed significantly higher migration and transition indices than those in the high-expression group (P < 0.05) (Fig. 4f), even though their expansion indices were comparable (Supplementary Fig. S4f). These findings suggest that although peripheral CD8+ T cells in the low-expression group highly expressed various naïve markers prior to treatment (Supplementary Fig. S4g), they acquired a stronger ability to migrate into the tissue microenvironment and differentiate into tumor-reactive CD8+ T cells during therapy. In contrast, CD8+ T cells from patients in the high-expression group, despite exhibiting greater baseline peripheral immune inflammation, demonstrated a weaker response and mobilization capacity upon treatment.
On the basis of the origin of tumor-reactive CD8+ T cells, we classified all post-immunotherapy tumor-reactive CD8+ T cells in the TME into three groups: Pre-existing, PBMC-derived, and Unknown. Pre-existing cells are CD8+ T cells carrying TCR clonotypes already present in the pre-treatment TME, whereas PBMC-derived cells are CD8+ T cells carrying TCR clonotypes detected in PBMCs either before or after treatment (Fig. 4g). The Unknown category includes potential contributions from other tissues and TCR clonotypes that were not captured in the pre-treatment TME or PBMCs. To characterize these Unknown clones, we employed the Le and Gascuel substitution model to measure their evolutionary distance from known PBMC clones. Unknown clones in the high-HLA-DQA2 expression group exhibited significantly greater evolutionary distances from PBMCs (P < 0.001) (Fig. 4h), suggesting that they are more likely derived from other tissues or TCR clonotypes not captured in the pre-treatment TME. In contrast, Unknown clones in the low-expression group were more likely to represent TCR clonotypes not detected in PBMCs. We then compared the functional profiles of the three groups of TCR origin and found that, although Pre-existing and PBMC-derived tumor-reactive CD8+ T cells displayed comparable exhaustion scores (Fig. 4i), PBMC-derived tumor-reactive CD8+ T cells exhibited significantly higher effector scores, indicating stronger antitumor activity (P < 0.0001) (Fig. 4j).
Overall, these observations suggest that when the TCR clonotypes that expand significantly after immunotherapy are derived primarily from pre-treatment intratumoral CD8+ T cell pools, particularly those already partially expanded under the chronic state of inflammation, as observed in patients with high HLA-DQA2 expression, they may contribute to an immune-tolerant state. This may be because the tumor antigens recognized by these pre-existing TCR clonotypes have already triggered immune evasion or are suboptimal for eliciting a robust cytotoxic response, thereby limiting their effectiveness in tumor clearance. This mechanism may also explain the poor response of pMMR/MSS CRC patients to immunotherapy alone, in which expansion is limited mainly to ineffective pre-existing CD8+ T cell clones, resulting in an insufficient generation of novel, tumor-reactive populations. We define this phenomenon as clonal entrapment, a state in which, during sequential radiotherapy and immunotherapy, pre-existing local inflammation drives the post-treatment expansion of TCR clonotypes that originate predominantly from the pre-treatment TME, whereas chronic inflammation in PBMCs restricts the migration and transition of potentially tumor-reactive peripheral CD8+ T cells, ultimately resulting in suboptimal treatment responses.
Residual and treatment-resistant tumor cells exhibit increased GDF15 expression
Next, we sought to investigate the second question: whether resistant tumor cells use specific mechanisms to suppress the infiltration of peripheral CD8+ T lymphocytes, thereby exacerbating clonal entrapment. First, we calculated a Tumor_score and a Normal_score for each epithelial cell and subtracted the Normal_score from the Tumor_score (Methods). Moreover, we used the R package inferCNV to infer copy number variations (CNVs) in each epithelial cell. Tumor cells were defined as those exhibiting high tumor scores in conjunction with apparent CNVs (Supplementary Fig. S5a, b). Despite removing batch effects prior to epithelial cell integration, these tumor cells still exhibited substantial inter-patient heterogeneity. To assess the overall chromosomal instability (CIN) status of each sample quantitatively, we integrated a published WES dataset of LARC patients48 and calculated the fraction of the genome altered (FGA) both before and after treatment. Interestingly, tumors that exhibited decreased CIN levels after treatment generally had higher baseline CIN levels prior to therapy (P < 0.001), suggesting a potential adaptive mechanism in which tumor cells with a higher CIN were selectively eliminated, while those with a lower CIN, which are potentially resistant to therapy, persisted (Supplementary Fig. S5c–e). Although higher levels of tumor aneuploidy have been associated with worse prognosis and immune evasion in several cancers, such as melanoma, following immunotherapy49,50, findings suggest that in the context of combined radiotherapy and immunotherapy, patients with lower immunogenicity and greater aneuploidy may actually have better outcomes51.
Gene co-expression modules in tumor cells were then analyzed using Hotspot52. Ultimately, we identified 17 gene expression modules and selected genes from each module with a Z score greater than 150, which reflect the autocorrelation strength between the gene and the module, for enrichment analysis. We displayed the top three significantly enriched pathways to define the function of each module. Five modules (Module_4, Module_8, Module_9, Module_11, and Module_14) were excluded from further discussion because the enriched pathways were not statistically significant (P.adjust > 0.05) or because no pathways were enriched (Fig. 5a; Supplementary Fig. S6). Among the remaining 12 modules, nine were detected in pre-treatment tumor cells only. These included Module_1 (MKI67 and PTTG1) and Module_2 (MCM7 and STMN1), which are associated with cell cycle progression; Module_3 (H2AC11 and H3C10), which is associated with chromatin organization; and the well-characterized Metal module (Module_10), Interferon module (Module_12), and pEMT module (Module_15)53, all of which were no longer detectable after radiotherapy (Fig. 5b).
Fig. 5: Residual and treatment-resistant tumor cells exhibit upregulated GDF15 expression.
a Heatmap showing the 17 coexpression gene modules identified in tumor cells. Each module and its representative genes are labeled. b Bubble plot showing module expression in tumor cells before treatment (Pre), after radiotherapy (In), and after immunotherapy (Post). c GDF15 staining of tumor tissues acquired after immunotherapy from representative patients who achieved a CR (P1 and P5) and patients who achieved a PR (P7 and P20). Scale bars: 700 μm (left), 200 μm (right). d Chord diagram (left) and bubble plot (right) showing cell–cell interactions between tumor cells, DC subtypes, and CD8+ T cell subtypes. e The GDF15-positive area (%) determined by IHC staining was quantified using ImageJ. The data are presented as the mean ± SD. **P < 0.01, two-tailed unpaired Student’s t-test. f Bubble plot showing the expression of Module_13 and its representative genes, GDF15 and BMP2, in tumor cells from public datasets. I–IV represent sampling time points from pre-treatment (I) to post-immunotherapy (IV).
Following the completion of radiotherapy, the expression of two modules, antigen processing and presentation via MHC-II (Module_5) and p53-mediated apoptosis (Module_16), was upregulated in tumor cells. These findings suggest that radiotherapy activated the p53 pathway by inducing DNA damage, triggering programmed cell death in tumor cells and inducing an immune-inflamed response that increased MHC-II expression in tumor cells. This increase in immunogenicity may synergize with subsequent immunotherapy to promote anti-tumor effects. However, following immunotherapy, the expression of these two modules decreased in the remaining tumor cells, whereas Module_13 emerged as the dominant module. Since the expression of Module_13 increased after radiotherapy, we defined it as a drug resistance gene module expressed by tumor cells. Enrichment analysis of genes within this module revealed that tumor cells may adapt to treatment pressure by promoting angiogenesis (placenta blood vessel development), maintaining cell adhesion (hemidesmosome assembly), activating pro-survival signaling pathways (response to peptide hormone), acquiring stem cell-like properties (endoderm formation), and regulating their proliferation (regulation of epithelial cell proliferation) (Fig. 5b; Supplementary Fig. S6).
Here, we found that growth differentiation factor 15 (GDF15) and CEACAM1 were the two proteins with the highest Z-scores among those involved in intercellular communication within Module_13 (Fig. 5a). Notably, cell–cell communication analysis revealed that only GDF15 interacted with immature cDC_CD1C cells as well as multiple CD8+ T cell subtypes (Fig. 5c, d). Recently, GDF15 was identified as a key mediator of immunotherapy failure in patients with non-small cell lung cancer and urothelial cancer54, and it may similarly contribute to treatment resistance in patients with LARC. Tumor cell-derived GDF15 inhibits DC maturation, suppresses T lymphocyte activation and recruitment, and contributes to poor prognosis55. We collected tissue samples from the LARC cohort (n = 20) in our study and from an additional cohort at Peking University Cancer Hospital (n = 10) following radiotherapy combined with immunotherapy. GDF15 was highly expressed in residual tumor cells from the PR group only (P < 0.01) (Fig. 5e). Moreover, we observed a similar pattern in the dataset from Zhang et al.: as treatment progressed, the expression of GDF15 and BMP2 (which is also a member of the transforming growth factor-beta (TGF-β) superfamily) and that of Module_13 steadily increased in the remaining tumor cells, supporting their potential role in the development of resistance to immunotherapy (Fig. 5f).
GDF15 inhibits CD8+ T cell infiltration and reduces the efficacy of combined radiotherapy and immunotherapy
To investigate the role of GDF15 in mediating LARC tumor resistance, we engineered the murine CRC cell line CT26 to overexpress or knockdown GDF15 (CT26-Gdf15-OE and CT26-shGdf15, respectively) (Supplementary Fig. S7a–c). Using these cell lines, we established a subcutaneous xenograft tumor model in mice and administered a long-term radiotherapy regimen combined with immunotherapy (Fig. 6a). Sequential radiotherapy and immunotherapy were effective across all three groups of mice (Fig. 6b, c). However, tumors in the GDF15 overexpression group exhibited the most significant resistance to treatment (P < 0.001), whereas tumors with GDF15 knockdown demonstrated the most pronounced regression after combination therapy (Fig. 6d). Furthermore, immunohistochemistry revealed that, compared with that in the other groups, the GDF15 expression in residual tumor cells in the GDF15-overexpressing group remained high (Supplementary Fig. S7d, e).
Fig. 6: GDF15 inhibits CD8+ T cell infiltration and reduces the efficacy of combined radiotherapy and immunotherapy.
a Schematic of the subcutaneous xenograft tumor model and treatment schedule combining long-term radiotherapy with immunotherapy. b Macroscopic view of subcutaneous xenograft tumors after combined long-term radiotherapy and immunotherapy. c Growth curves of subcutaneous xenografts in mice during treatment. The data are shown as the mean ± SEM. ****P < 0.0001, as determined by the Wilcoxon rank-sum test. d Tumor weights of the subcutaneous xenografts after treatment. The data are shown as the mean ± SEM. **P < 0.01, ***P < 0.001, ****P < 0.0001, two-tailed unpaired Student’s t-test. e CD8+ T-cell migration (%) in Transwell assays using conditioned media from CT26-Gdf15-OE cells, CT26-shGdf15 cells, and their respective controls. f Representative flow cytometry density plots of lower-chamber events (x-axis: SSC-A; y-axis: FITC-A) gated on CFSE+ CD8+ T cells; counting beads were used for normalization. g The CD8A-positive area (%) determined by IHC staining was quantified using ImageJ. The data are shown as the mean ± SD. ns, P > 0.1, *P < 0.1, **P < 0.01, two-tailed unpaired Student’s t-test. h CD8A staining of mouse tumor tissues from the CT26-Gdf15-OE, CT26-shGdf15, or their respective control groups after treatment. Scale bars: 100 μm (left), 25 μm (right). i Heatmap showing correlations between GDF15 expression and CD8+ T cell marker genes, along with related GSVA scores, in the CRC bulk RNA-seq public dataset. ***P < 0.001, Spearman correlation analysis. j Bubble plot showing the expression of GDF15 and Module_13 in tumor cells from each patient in the PR group after treatment. k The mean fluorescence intensity of three cell lines in nine different regions (100× magnification) was quantified using ImageJ. The data are shown as the mean ± SD. *P < 0.1, ****P < 0.0001, two-tailed unpaired Student’s t-test.
Previous studies and our analysis of public CRC datasets revealed that high expression of GDF15 is significantly negatively correlated with CD8+ T cell infiltration in the TME and with the expression of various effector and exhaustion markers in CD8+ T cells (Fig. 6i). Consistent with these findings, our CFSE-based Transwell migration assay demonstrated that tumor-derived GDF15 suppressed CD8+ T-cell chemotaxis (Fig. 6e, f; Supplementary Fig. S7f). To further investigate this, we examined CD8+ T cell infiltration in tumor tissues from different groups of mice following treatment. Compared with the control group, the GDF15 overexpression group exhibited significantly lower CD8⁺ T cell infiltration (P < 0.01) (Fig. 6g, h). Although some mice with GDF15-knockdown tumors exhibited increased CD8+ T cell infiltration, overall, there was no significant difference between the GDF15 knockdown and control groups. One possible explanation is that the GDF15 knockdown group had the lowest tumor burden after treatment, resulting in a reduced inflammatory response in the TME. Alternatively, these tumor cells may retain additional resistance mechanisms. Notably, in the PR group, the TRG 3 patient (P11) presented the highest expression of GDF15 and Module_13, whereas the TRG 1 patient (P09) presented lower expression levels of both (Fig. 6j). These findings suggest that in patients such as P09, other tumor resistance mechanisms or interacting pathways may be involved, which warrants further investigation.
In vitro, compared with overexpressing cells, cells with knockdown of GDF15 presented a significant increase in the expression of the DNA damage marker γ-H2AX following exposure to 4 Gy of radiation (P < 0.0001) (Fig. 6k; Supplementary Fig. S7g), suggesting impaired DNA repair mechanisms. This impaired DNA repair capacity may lead to genomic instability, which is often associated with the accumulation of CNVs. This result is particularly relevant in light of our earlier findings, in which resistant tumor cells exhibited decreased CNV levels and increased GDF15 expression. The reduction in CNV levels in resistant cells may reflect a more stable genome, potentially maintained by GDF15, as tumor cells themselves may also express the GDF15 receptor TGFBR2, similar to CD8+ T cells (Fig. 5d).
To support the role of GDF15 in mediating resistance to radiotherapy, we analyzed two independent public LARC datasets. The GSE209746 dataset includes bulk RNA-seq data from the tumor tissues of 85 MSS LARC patients prior to radiotherapy32. We found that the patients with the poorest treatment response (TRG 3) exhibited the highest baseline expression of GDF15 (Supplementary Fig. S7h). Additionally, recently published single-cell sequencing data on tumor tissues from MSS LARC patients before and after radiotherapy56 revealed that residual tumor cells post-radiotherapy had significantly increased expression of GDF15 (P < 0.0001) (Supplementary Fig. S7i).
In summary, we established high-resolution, longitudinal single-cell transcriptomic atlases of tumor tissues and peripheral blood from LARC patients undergoing sequential radiotherapy and immunotherapy. Our profiling led to the identification that clonal entrapment is a key feature associated with limited therapeutic efficacy in patients with pMMR/MSS tumors — a state characterized by high baseline HLA-DQA2 expression in DCs and GDF15 expression in resistant tumor cells. Specifically, patients with elevated HLA-DQA2 expression exhibit a pre-existing chronic inflammatory microenvironment characterized by the accumulation of mature DCs and terminally exhausted CD8+ T cells. This hyperactivated but functionally compromised state is correlated with the limited presence of the naïve DCs required to present tumor neoantigens. Furthermore, GDF15 expression in resistant tumor cells is associated with reduced infiltration of potential peripheral tumor-reactive CD8+ T cells into the TME. Consequently, the immune response following ICI therapy is observed to be “entrapped”; that is, expansion is confined predominantly to pre-existing intratumoral TCR clonotypes already partially expanded under chronic inflammation, rather than novel, high-effector peripheral clones. This lack of recruitment and expansion of novel peripheral clones indicates restricted TCR repertoire diversification, characterizing the clinical resistance observed in these patients (Fig. 7a, b).
Fig. 7: Single-cell atlas and mechanism of clonal entrapment.
a Schematic diagram of the in-house data used in this study (left) and the public datasets used for validation (right). b Schematic diagram of the specific mechanism of clonal entrapment.
Radiotherapy effectively reduces the abundance of THBS2+ CAFs, enabling immune infiltration during subsequent immunotherapy
Finally, we focused on the third question: The remodeling effects of radiotherapy on fibroblast subtypes and how these differ from those of chemotherapy-induced remodeling, given that chemoimmunotherapy remains one of the most commonly used neoadjuvant treatment strategies57,58. Radiotherapy effectively reduces the proportion of fibroblasts in the TME, which is important because cancer-associated fibroblasts (CAFs) promote the formation of “cold tumors” and contribute to tumor progression59. On the basis of marker gene expression, we classified all the fibroblasts in our dataset into the following 11 subtypes: Two normal fibroblast precursor subtypes (pFB_COL15A1 and pFB_PI16), one myofibroblast subtype (myCAF_WNT2), three inflammatory fibroblast subtypes (iCAF_AREG, iCAF_CXCL5, and iCAF_CXCL10), and five other subtypes (including CAF_SOX6 and CAF_ADAMDEC1, which are considered to be specific to the gastrointestinal microenvironment60) (Supplementary Fig. S8a, b).
Radiotherapy effectively reduced the proportion of or even eliminated the myCAF subtype (Cluster 9, myCAF_WNT2) in the TME, which has been associated with poor prognosis in CRC61,62. We further analyzed the expression of the myCAF marker gene CTHRC1 and the pro-tumorigenic factors WNT2 and THBS2, which are secreted by myCAFs, across the treatment timeline. The expression of these genes significantly decreased in the CR and PR groups after radiotherapy (P < 0.0001); however, baseline expression levels of these factors were higher in the PR group than in the CR group before treatment (Supplementary Fig. S8c). Following immunotherapy, CAFs were largely remodeled toward a normal fibroblast state, characterized by high expression of ADAM28 and ADAMDEC1. These findings were corroborated by data from two public CRC scRNA-seq datasets (SMC and KUL), which included adjacent normal samples and tumor-normal junction samples63 (Supplementary Fig. S9b–d).
Given recent reports that THBS2+ CAFs contribute to chemotherapy resistance in the context of CRC64, we analyzed a scRNA-seq dataset from 29 patients with pMMR/MSS RC treated exclusively with chemotherapy, among which paired tumor samples were collected from 27 patients before and after treatment38. Moreover, we obtained pre- and post-chemotherapy LARC samples from two other patients (P26 and P28). We mapped the CAF subtypes from our dataset to these data and discovered that the previously identified myCAF_WNT2 subtype could be further divided into two distinct populations: myCAF_THBS2 (CTHRC1, THBS2, and INHBA) and myCAF_WNT4 (WNT4 and WNT5A), both of which express myCAF marker genes, including ACTA2, TAGLN, and WNT2 (Supplementary Fig. S8d, e).
Comparative analysis revealed that chemotherapy significantly reduced the proportion of the myCAF_WNT4 subtype (P < 0.0001) but had a minimal effect on the myCAF_THBS2 subtype (P = 0.13). In contrast, radiotherapy significantly reduced the proportion of both subtypes (myCAF_WNT4: P < 0.0001; myCAF_THBS2: P < 0.0001) (Supplementary Fig. S8f, g), suggesting that radiotherapy is more effective than chemotherapy at targeting CAF populations associated with tumor progression and drug resistance, particularly the myCAF_THBS2 subtype, and that this effect was consistently observed in both patients who achieved a CR and those who achieved a PR (Supplementary Fig. S8h, i; Supplementary Fig. S9a). We also compared the effects of radiotherapy and chemotherapy on endothelial cell dynamics. Radiotherapy more effectively induced the expression of CD74 and various MHC-II molecules in endothelial cells, facilitating their conversion into non-professional antigen-presenting cells. The upregulation of adhesion molecules such as MADCAM1 and ICAM1 in endothelial cells after radiotherapy also effectively promoted lymphocyte infiltration into tumor tissue, providing a favorable environment for subsequent immunotherapy (Supplementary Fig. S9e, f).
