Stress‐Programmed Immune Niches Fuel TNFR2+ Treg Activation and Drive Neoadjuvant Chemotherapy Resistance in Breast Cancer
Zhibo Shao, Xuliren Wang, Han Zhu, Ling Yao, Qi Zhang, Zihan Zhai, Xinyi Lv, Xinyuan Xia, Yi Zhang, Xinya Lu, Xinyi Quan, Jiawen Yang, Xinwei Li, Jieyi Dong, Zhi‐ming Shao, Ruoxi Wang, Bingqiu Xiu, Jiong Wu, Sheng Chen
TL;DR
Stress-programmed immune cells in breast cancer promote chemotherapy resistance by activating Tregs through TNFα–TNFR2 signaling, and blocking this pathway could improve treatment outcomes.
Contribution
The study identifies a stress-programmed immune module and its role in chemotherapy resistance via TNFα–TNFR2 signaling in breast cancer.
Findings
Stress-programmed immune states are linked to elevated TNFα signaling and chemotherapy resistance in breast cancer.
Blocking the TNFα–TNFR2 axis suppresses tumor growth without toxicity, suggesting a new therapeutic strategy.
Single-cell data reveals coordinated heat shock gene expression across immune cells in nonresponders to chemotherapy.
Abstract
The tumor microenvironment (TME) harbors diverse immune cell states that shape therapeutic outcomes in breast cancer. Here, we identify a conserved stress‐programmed cellular module as a key responder to neoadjuvant therapies in breast cancer, characterized by coordinated heat shock gene expression across multiple immune cells, based on single‐cell transcriptomic data from neoadjuvant chemotherapy‐treated patients. We discover that this multicellular program enhances the effector fate of regulatory T cells (Tregs) via chronic and TME‐wide TNFα signaling, compromising the efficacy of neoadjuvant chemotherapy. TNFα signaling, typically considered an antitumor cytokine, is paradoxically elevated in nonresponders both pre‐ and post‐treatment, with a particularly prominent TNFα–TNFR2 interaction. Blocking this axis, with or without chemotherapy, significantly suppresses tumor growth without…
Genes, proteins, chemicals, diseases, species, mutations and cell lines named across the full text — each resolved to its canonical identifier and authoritative record.
Click any figure to enlarge with its caption.
FIGURE 1
FIGURE 2
FIGURE 3
FIGURE 4
FIGURE 5
FIGURE 6
FIGURE 7- —National Natural Science Foundation of China10.13039/501100001809
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Taxonomy
TopicsCancer Immunotherapy and Biomarkers · Cancer, Stress, Anesthesia, and Immune Response · Immune cells in cancer
Introduction
1
Immune infiltrates within the tumor microenvironment (TME) plays a decisive role in shaping the response to clinical outcomes of cancer treatment [1, 2, 3]. However, immune cell subsets in TME are phenotypically and functionally diverse and their characteristics determine the effectiveness and potential side effects of anticancer therapies [4, 5]. Even within the same immune cell type, distinct cellular states can lead to markedly different functions [6, 7, 8]. The heterogeneity of immune cellular states is tightly linked to treatment response, particularly in the setting of neoadjuvant chemotherapy [9, 10]. However, the complexity of immune cell states within the TME and their influence on treatment response remain to be fully elucidated.
Heat shock proteins (HSPs) are a family of highly conserved molecular chaperones that are markedly upregulated in response to cellular stress [11, 12, 13]. HSPs are highly expressed in cancers and closely associated with tumorigenesis and progression [14, 15, 16]. Recently, among the diverse immune cell states, a subset of immune cells in the TME can acquire stress‐related programs characterized by heat shock gene expression. A pan‐cancer T atlas encompassing over 300 000 cells has delineated a distinct T cell subset marked by elevated expression of heat shock‐related genes, which has been implicated in resistance to neoadjuvant immunotherapy across diverse tumor types [17]. Another large‐scale pan‐cancer single‐cell analysis of NK cells demonstrated that those with elevated heat shock program exhibit diminished cytotoxic potential and increased expression of immune checkpoint molecules, enriched within tumors as tumor‐associated NK cells [18]. In B cells, a subset marked by elevated heat shock program activity has been identified and is associated with unfavorable tumor prognosis [19, 20]. Although several studies have identified immune cell states dominated by heat shock programs and associated with adverse tumor outcomes, their connection and underlying mechanisms remain unclear.
In this study, we comprehensively analyzed single‐cell RNA sequencing data from breast cancer patients undergoing neoadjuvant chemotherapy and uncovered a stress‐programmed immune module defined by heat shock gene expression. This module spans multiple immune cell types, including CD8⁺ T cells, CD4⁺ T cells, B cells, and myeloid cells. Moreover, heat shock–enriched immune subtypes rarely occur in isolation, as different subtypes tend to be strongly correlated. Correlation analysis of immune niche revealed that the stress‐programmed immune niche is positively associated with Treg‐enriched environments and negatively correlated with cytotoxic‐enriched niches. The TNFα–TNFR2^+^Treg axis is a central pathway mediating the immunosuppressive effects of the stress‐programmed immune cellular module.
Result
2
A Single‐Cell Atlas of Breast Cancer TME During Neoadjuvant Chemotherapy
2.1
To investigate the impact of neoadjuvant chemotherapy, heterogeneity of immune cell states and their interactions on the TME of breast cancer, we collected 11 breast cancer samples from 7 patients who received 6 cycles of neoadjuvant chemotherapy based on taxane and platinum for single‐cell sequencing (Figure 1A). Among these patients, four provided single‐cell sequencing data from both pretreatment biopsies and post‐treatment surgical samples. The remaining three patients had only pretreatment biopsy samples. Three were diagnosed with tripe‐negative breast cancer, while four were classified as estrogen receptor‐positive. Based on the pathological response evaluated from surgical pathology slides, we categorized the patients into two groups: responsive (n = 4) and nonresponsive (n = 3) (Figure 1B and Table S1).
The dynamic single‐cell landscape of breast cancer patients during neoadjuvant chemotherapy. (A) Schematic representation of the patient cohort and experimental design. (B) Clinical characteristics and treatment outcomes of enrolled patients. (C) Dotplot illustrating the expression of specific markers across major cell types. (D) Uniform manifold approximation and projection (UMAP) plot showing integrated major celltypes. Major cell types, including Epithelial cells, Smooth muscle cells, Endothelial cells, Fibroblasts, T/NK cells, B cells, Myeloid cells, and Plasma cells, are identified in distinct colors. (E) UAMP plot showing Immune subtype clusters. The major compartment identity of each cell is shown in the inset on the right.
The single‐cell sequencing data were processed through a strict quality control (QC) that included the removal of doublets and cells with high RNA contamination, in addition to standard quality control (Figure S1A,B). Finally, we integrated 68 534 cells that met quality control criteria and categorized them into 9 major cell types (Figure 1D). In brief, epithelial cells, fibroblasts, endothelial cells, smooth muscle cells, and immune cells—specifically mast cells, B cells, plasma cells, myeloid cells, and T/NK cells—were classified in this study based on the expression patterns of well‐established marker genes (Figure 1C). The uniform manifold approximation and projection (UMAP) embeddings showed no apparent segregation by timepoint (Figure S1C), patient identity (Figure S1D), or treatment response (Figure S1E), suggesting successful integration and the absence of overt batch effects in downstream clustering. All sample contained 9 cell types, and the cell count distribution was balanced, indicating a comprehensive and robust sampling of the breast cancer microenvironment during neoadjuvant chemotherapy (Figure S1F).
Immune cell subsets are essential components of the tumor immune microenvironment, and their transcriptional heterogeneity critically shapes their functional states and regulatory roles in tumor progression and therapeutic response [21, 22, 23]. To further examine TME heterogeneity and its impact on breast cancer patients’ outcomes at a finer resolution, we conducted precise clustering of the immune cells. To this end, we constructed a hierarchical framework encompassing all major immune cell types, guided by comprehensive literature curation of transcriptional evidence supporting the definition of 31 distinct immune cell states (Figure 1E and Figures S1G and S2). For example, our annotation revealed four and six states in the CD4^+^ and the CD8^+^T cell compartments, respectively (Figures S1G and S2B). Moreover, a subset of T cells coexpressing high levels of CD4 and CD8—termed double‐positive (DP) T cells—has been previously characterized as a naïve‐like population [24, 25]. Intriguingly, we annotated stress‐responsive subsets across multiple immune lineages, including B cells (B_c3, B_c4), CD4⁺ T cells (CD4T_c2), CD8⁺ T cells (CD8T_c5), and macrophages (Macro_c5_HSPA1B), suggesting a convergent stress‐programmed state emerging in distinct cell populations. Together, this high‐resolution single‐cell dataset provides a valuable resource for understanding how neoadjuvant chemotherapy reshapes immune states in the breast cancer TME.
Stress‐Programmed Immune Cell States are Enriched in Nonresponders and Predict Poor Treatment Outcomes
2.2
Intertumoral TME heterogeneity established during tumorigenesis and dynamic variations during neoadjuvant chemotherapy often portends diverse patient outcomes [26], prompting us to examine the dynamic changes of immune cell subsets during treatment. We first constructed a comprehensive cell‐state hierarchy covered all immune major cell types using R o/e value [6] and compared cell‐state proportions across treatment responses and treatment phases, identifying a wide range of cell‐state change (Figure 2A).
*Stress‐programmed immune subtype states predict poor outcomes to neoadjuvant chemotherapy. (A) Top: Unsupervised hierarchical clustering of defined TME cell states. Bottom: Heatmap showing Ro/e values of immune cell clusters across three comparisons: (1) nonresponders (NR) versus responders (R) at baseline, (2) surgery versus baseline in responders, and (3) surgery versus baseline in nonresponders. Clusters with high HSPA1B expression are highlighted in red. (B) Shared top‐ranking cell populations across the three differential abundance comparisons described above. (C) Boxplots comparing the stress score of CD4+ T cells, CD8+ T cells, B cells, and macrophage clusters between baseline and surgery samples in responders and nonresponders. Statistical comparisons were performed using two‐sided Wilcoxon tests. (D) Schematic of the 4T1 mouse tumor implantation model and treatment regimen. (E) Representative images of tumors from PTX‐resistant and PTX‐effective mice at endpoint. (F) Quantification of residual tumor weight on day 27 in resistant versus effective groups. Data are shown as mean ± SEM; statistical analysis was performed using a two‐tailed unpaired t‐test. (G) Tumor growth curves of individual mice in resistant and effective groups. (H) Difference in the proportion of HSP70high cells among CD45⁺ immune cells between resistant and effective samples according to flow cytometry results. (J) Correlation of the percentage of CD45⁺HSP70⁺ immune cells and residual tumor burden across individual mice. (K) Dataset from Hengstler et al. (n = 91) show molecular subtype distribution and elevated stress‐programmed scores in nonresponders (NR) relative to responders (R) (p < 0.001). (L) Dataset from Korde et al. (21 breast cancer patients receiving neoadjuvant chemotherapy) shows subtype distribution and changes in stress‐programmed scores from baseline to cycle 2 (C2) in nonresponders (NR) and responders (R) (p = 0.074 and 0.69, respectively). ns, not significant; *p < 0.05; **p < 0.01; ***p < 0.001; ***p < 0.0001.
To focus on the most convergent signals, we identified four overlapped cell states across the three comparisons‐namely heat shock‐responsive CD8⁺ T cells (CD8T_c5_HSPA1B), CD4⁺ T cells (CD4T_c2_HSPA1B), and two stress‐associated B cell subsets (B_c3_HSPA1B_1 and B_c4_HSPA1B_2)—all of which were more abundant in treatment‐naive nonresponders and remained elevated at surgery, while showing lower levels in responders at baseline (Figure 2B), reflecting a conserved stress‐response program across immune lineages. This shared stress signature prompted us to validate its link to treatment response [17, 18, 19]. Given that, we first performed AUCell scoring, based on a curated heat stress gene set, to quantify the stress‐related programs across immune cell subsets (Figure S3A). Consistent with our earlier findings, stress scores increased in nonresponders after treatment and declined in responders, across CD4⁺ T cells, CD8⁺ T cells, B cells, and macrophages (Figure 2C), with directionally consistent trends despite limited significance in B cells from responders and CD8⁺ T cells from nonresponders. Although not among the initially identified cell states, stress signatures in macrophages exhibited similar trends, warranting their inclusion in downstream analysis. Second, we observed significantly higher expression of HSPA1A and HSPA1B in CD4⁺ T cells, CD8⁺ T cells, macrophages, and B cells from nonresponders compared to responders at both baseline and surgery, in both our in‐house cohort (Figure S3B) and the dataset from Zhang et al. (Figure S3C). These findings link stress‐related immune states to resistance, warranting further study.
To further test this observation, we conducted an in vivo mouse experiment (n = 12) in which mice received five doses of paclitaxel (PTX) at 200 µg per mouse, administered via intravenous injection every 3 days (Figure 2D). At the experimental endpoint, mice were stratified into resistant and effective groups based on residual tumor weight, to investigate the relationship between post‐treatment tumor burden and heat stress signaling within immune cells (Figure 2E). Residual tumor weight was significantly higher in the resistant group compared to the effective group (Figure 2F), confirming the dichotomy in therapeutic response. Tumor growth was suppressed in the effective group but continued to progress in the resistant group despite PTX treatment (Figure 2G). We next assessed HSP70 expression in tumor‐infiltrating immune cells from individual mouse tumor samples using flow cytometry. HSP70, a major stress‐inducible chaperone encoded by HSPA1A and HSPA1B, reflects activation of the heat stress response at the protein level [27, 28, 29]. Across samples, we observed that mice with larger residual tumor weight exhibited a higher proportion of HSP70^high^ immune cells (Figure 2H–I). In addition, correlation analysis revealed a strong positive association between the proportion of HSP70 ^high^ immune cells and residual tumor weight in the mice cohort (Figure 2J), suggesting a positive association between immune stress signaling and treatment resistance.
To determine the clinical impact of stress‐programmed immune cell states, we analyzed HSPA1B expression‐based scores across immune subsets in the TCGA breast cancer cohort. Kaplan–Meier survival analysis revealed that high HSPA1B scores in CD4⁺ T cells, CD8⁺ T cells, B cells, and macrophages were each significantly associated with worse overall survival (Figure S3D). To validate the predictive effect of the stress‐programmed signatures on neoadjuvant chemotherapy response, we analyzed two phase II clinical cohorts of breast cancer patients: the dataset reported by Jan G. Hengstler et al. [30], which includes 91 patients treated sequentially with docetaxel followed by the FEC regimen (5‐fluorouracil, epirubicin, and cyclophosphamide), and the dataset reported by Larissa A. Korde et al. [31], comprising 21 paired tumor samples collected before and after docetaxel + capecitabine (TX) treatment. The analysis revealed that stress‐programmed scores were significantly enriched in nonresponders but tended to decrease after treatment in responders, although not reaching statistical significance; this pattern was not observed in nonresponders, and both cohorts encompassed patients across all major molecular subtypes of breast cancer (Figure 2K,L). These results suggest that the marked changes in stress programmed signatures within immune subsets may underlie the adverse outcomes of neoadjuvant treatment.
Trajectory Characteristics of Stress‐Programmed B Cells and Macrophages
2.3
Previously reported work has characterized stress‐programmed T cells [17], in which stress‐programmed CD8⁺ T cells occupied a distinct trajectory, separate from both effector and exhausted CD8⁺ T cell populations, suggesting that this represents a unique and previously underexplored immune cell state. Consistently, we also observed stress‐programmed CD4⁺ and CD8⁺ T cells within breast cancer tissue using multiplex immunofluorescence (Figure S4A). However, stress‐programmed states in B cells and macrophages remain insufficiently explored, prompting us to explore their trajectory dynamics.
To characterize the differentiation trajectory of stress‐programmed B cells, we first confirmed their presence in tumor tissues through multiplex immunofluorescence staining (Figure 3A). Then, monocle3 trajectory analysis revealed two primary differentiation trajectories paths of B cells, originating from a shared naïve‐like states (Figure 3B). Within the shared naïve‐like branch point, B_c1_TCL1A and B_c5_CD70 accounted for 34.6% and 29.7% of the population, respectively. Downstream, Path 1 was predominantly composed of the stress‐programmed subset B_c3_HSPA1B_1 (41.4%), while Path 2 was mainly driven by the interferon‐associated subset B_c2_IFN (35.5%) (Figure 3C). To further characterize the differentiation trajectories, we projected the expression of B cell activation, HLA, interferon (IFN), stress related gene signatures onto the monocle3‐defined UMAP embedding, each of which followed its expected dynamic pattern along pseudotime (Figure 3E). The two main paths imply divergent cell fates, which are supported by expression kinetics of related gene signatures along the inferred pseudotime axis (Figure 3D).
Pseudotime trajectory and stress‐programmed differentiation states in B cells and macrophages. (A) Representative multiplex immunofluorescence images of stress‐programmed B cells. (B) Monocle 3 trajectory analysis of B cells cell differentiation revealing two main divergent trajectories. Cells are color coded for their corresponding pseudotime. (C) Distribution of B cell subclusters across Naive, Path 1, and Path 2 states. (D) Dynamics of signature scores for B cell activation, HLA, interferon (IFN), and Stress response along Path 1 and Path 2 pseudotime. (E) UMAP visualization of B cell subclusters generated by Monocle3 and their corresponding signature scores for B cell activation, HLA, IFN, and Stress pathways. (F) Representative multiplex immunofluorescence images of stress‐programmed macrophages. (G) Monocle 3 trajectory analysis of monocyte/macrophage cell differentiation revealing one main trajectories. Cells are color coded for their corresponding pseudotime. (H) Ridge plots displaying the distribution of inferred pseudotime across monocyte/macrophage subclusters. (I) Distribution of monocyte/macrophage subclusters across three differentiation phases and the dynamics of Stress, M1, and M2 signature scores across pseudotime trajectory. (J) UMAP visualization of monocyte/macrophage subclusters generated by Monocle3 and their corresponding signature scores for Stress, M1, and M2 pathways.
As for macrophages, we similarly validated the presence of stress‐programmed subsets within the tumor microenvironment (Figure 3F). Monocle3 trajectory analysis revealed a single predominant differentiation path in macrophages, beginning with Mono_c0_AQP9 and Mono_c1_CXCL3 and progressing toward the Macro_c0_SPP1 state, with the stress‐programmed subset (Macro_c5_HSPA1B) occupying an intermediate position along the trajectory (Figure 3G and H). To elucidate the role of stress signaling in macrophage differentiation, we projected M1 phenotype, M2 phenotype, and stress‐related gene signatures onto the Monocle3‐defined UMAP embedding, which demonstrated that stress programs were initiated early and sustained throughout the differentiation continuum (Figure 3J). We defined three differentiation phases along the macrophage trajectory based on Monocle3 analysis: Phase 1 featured low stress signatures and moderate M1 signatures, Phase 2 exhibited strong M1 polarization along with elevated M2 markers, suggesting a transitional state, and Phase 3 was marked by elevated M2 and diminished M1 expression. These phases corresponded to distinct subsets—Mono_c0_AQP9 (23.8%) and Mono_c1_CXCL3 (29.6%) in Phase 1, Marco_c1_CCL18 (51.8%) in Phase 2, and Marco_c0_SPP1 (79.9%) in Phase 3—indicating a transition from early activation to terminal immunoregulatory states (Figure 3I).
These findings extend the concept of stress‐programmed states beyond T cells and underscore their broader relevance across immune lineages.
Stress‐Programmed Cell States Constitute a Core Coordinated Cellular Programs in the Tumor Immune Ecosystem
2.4
Since stress‐programmed states occur in multiple immune subsets, we conducted cell‐type co‐occurrence analyses across tumor patients to examine the interplay between stress‐programmed states and the tumor immune microenvironment. Macro_c6_TAM and DP_T cells were excluded due to their limited representation and functional redundancy (Figure S2B,S2C). In result, eight covarying multicellular modules (CMs) were identified (Figure 4A). With these CMs, we found that CM1 and CM2 were enriched with immune‐activating populations, including CXCL13 expressing T cells [32] and NK cells [33, 34], while CM5 was dominated by immunosuppressive regulatory CD4^+^T cells, and CM8 showed a marked aggregation of stress‐programmed immune states (Figure 4A). To delve deeper, CM1 and CM2 exhibited elevated expression of cytotoxic genes such as GZMK and GZMB, along with enhanced “Effector Signature” and chemokine signaling. In contrast, CM5 showed high expression of immunosuppressive markers, such as FOXP3 and IL2RA and was enriched in pathways related to immune checkpoints and Treg‐mediated suppression (Figure 4B). CM8 was devoid of both effector and checkpoint signatures but was instead characterized by prominent stress responses and mild TNFα–NF‐κB activity (Figure 4B and Figure S4C), reflecting a stress‐programmed but immune‐inactive cellular modules.
*The diversity of coordinated cellular modules in the TME and the distinct correlations of each module with CM8. (A) Heatmap displaying pairwise Pearson correlations of cell‐type proportions across tumor samples. Cellular modules are delineated by colored boxes at the top. (B) Dot plot showing the expression of specific marker genes and the expression levels of specific signature scores for samples categorized into cellular modules. (C) Cross‐cohort correlation of cellular modules with CM8 from the in‐house cohort (x‐axis) and the TCGA BRCA cohort (y‐axis). Mutually correlated with CM8 are marked with asterisks. (D) Pearson correlation plots depicting the relationships between CM8 niche fractions and those of CM1, CM2, and CM5 in both the in‐house and TCGA BRCA cohorts. (E) Kaplan–Meier survival curves of TCGA BRCA patients stratified by dominant module (CM1, CM2, or CM5). The significance of differences was performed by the log‐rank test, multiple testing correction (Benjamini–Hochberg FDR adjustment) was applied. (F) Left: Spatial transcriptomic analysis showing regional distributions of effector and CM8 signatures in two representative tumor regions (#1 and #2). Right: Boxplots indicate significantly reduced effector signature scores in high‐CM8 areas. Wilcoxon test. (G) Distribution of CM8 and effector signature scores along a scaled spatial distance axis. The x‐axis reflects distance from the tumor core, with lower values indicating proximity and higher values representing more distal regions. ns, not significant; *p < 0.05; **p < 0.01; ***p < 0.001; ***p < 0.0001.
Building on these observations, we further explored the function of CM8 within the immune context of tumors, we conducted covariation analyses with other modules via BayesPrism [35] deconvolution using both in‐house single‐cell data and the TCGA BRCA cohort (Figure S4B). The correlation patterns between CM8 and CM1, CM2, and CM5 were consistently observed in both our in‐house single‐cell dataset and the TCGA BRCA cohort (Figure 4C). Specifically, CM1 and CM2 were negatively correlated with CM8, whereas CM5 showed a positive correlation, as illustrated by the corresponding scatter plots (Figure 4D). Furthermore, survival analysis based on CM classification showed that patients enriched for CM5 had significantly worse overall survival than those with CM1 or CM2 (Figure 4E). These findings collectively position CM8 as a distinct, stress‐programmed cellular module, supporting its potential role in shaping an immune‐inert tumor microenvironment at both single‐cell and bulk transcriptomic levels.
Given the absence of spatial context in single‐cell and bulk analyses, we next leveraged spatial transcriptomics [35, 36] to examine the in situ distribution and niche associations of stress‐programmed cell states. We first analyzed a published spatial transcriptomics dataset [37] and used Cottrazm [38] package to define malignant, boundary, and nonmalignant regions within tumor spatial microenvironment (Figure S4D). Intriguingly, we observed a spatially inverse relationship between stress and effector immune signatures. Specifically, CM8 signature scores peaked in the tumor core and gradually declined toward the periphery, while effector signatures exhibited the opposite trend (Figure 4G). Consistently, analysis of an integrated single‐cell dataset comprising normal breast tissue [39], peripheral blood, and breast tumors [40] revealed a progressive increase in HSPA1A and HSPA1B expression in immune cells along the normal‐to‐tumor axis (Figure S4E,F). To refine the spatial interplay between CM8 and effector programs, we mapped their respective signatures across spatial coordinates and found that regions with high CM8 signature scores consistently exhibited diminished effector activity (Figure 4F), reinforcing their spatially antagonistic relationship. Together, these integrative analyses define CM8 as a spatially confined, immunologically inert cellular module, enriched for stress signaling yet devoid of effector activity, and underscore its potential role in shaping an immune‐excluded tumor niche.
TNFα Signaling Shapes the Intercellular Network of Stress‐Programmed Cellular Modules
2.5
To further elucidate how stress‐programmed cellular modules (CM8) functionally interacts with its surrounding microenvironment, we next examined its cell–cell communication landscape.
First, we employed CellChat [41] to decipher reciprocal ligand–receptor interactions quantitatively. Importantly, we noticed that the outgoing signaling from stress‐programmed cell states was comparable to that of other cells, whereas their incoming signaling was markedly reduced (Figure S5A,B), potentially reflecting an immune‐inert or communication‐impaired state that limits their responsiveness to external cues. Then, with NicheNet [42] analysis, we are able to identify a set of four genes (defined as core ligands) that were predicted as the potential ligands to drive the phenotype of the stress‐programmed cellular modules (Figure 5A). Of note, among the four core ligands, TGFB1 emerged as a canonical immunosuppressive regulator, known to modulate cancer‐associated fibroblast (CAF) activation and epithelial‐to‐mesenchymal transition (EMT) [43, 44, 45, 46]. TNF, a pleiotropic proinflammatory cytokine, can promote immune clearance but may also facilitate tumor progression under chronic inflammatory conditions [46, 47, 48]. ITGB2 encodes a critical integrin subunit essential for immune cell adhesion and migration [49]. These results suggested a role for CM8 in orchestrating leukocyte dynamics. Among them, TNF stood out as the principal mediator orchestrating the immunological role of CM8. The expression of TNF was highest in cell states within CM8 (Figure 5C), whereas the other ligands did not exhibit such specificity (Figure S5C). Moreover, the expression of TNF and HSPA1B displayed notable colocalization on the UMAP embedding (Figure 5D). To examine cell–cell communication involving specific stress‐programmed cell states, we analyzed their interactions with selected immune and stromal populations. We found that TNF signaling emerged as the most consistently engaged pathway, with stromal cells (e.g., epithelial cells and fibroblasts) interacting mainly via the TNF–TNFRSF1A axis, while immune cells preferentially engaged the TNF–TNFRSF1B axis (Figure 5B).
*Stress‐programmed cellular modules converge on a shared TNF signaling axis. (A) Venn diagram illustrating the overlap of the top 30 predicted ligands identified by NicheNet analysis. (B) Bubble heatmap showing the cell–cell communication of selected ligand–receptor pairs among specific cell types. Dot size indicates p value, colored by communication probability. (C) Ranking of immune cell states based on median TNF expression. (D) UMAP plots of TNF and HSPA1B expression across T cells, B cells and Myeloid. (E) Flow cytometry validation of TNFαhigh percentage and mean fluorescence intensity of TNFα in HSP70⁺ versus HSP70− immune cell subsets. Left panels show analysis in tumor samples from breast cancer patients (n = 6 for CD8T and B cells, n = 5 for CD4T); right panels show results from the 4T1 mouse tumor model (n = 9 for CD8T and B cells, n = 8 for CD4T). Statistical comparisons were performed using paired two‐tailed t‐test. ns, not significant; *p < 0.05; **p < 0.01; ***p < 0.001; ***p < 0.0001.
Next, we performed flow cytometry on both human breast cancer tissues and 4T1 xenograft tumors in mice to validate this observation. As anticipated, within B cells, CD4⁺ T cells, and CD8⁺ T cells, the HSP70^high^ subsets exhibited a higher proportion of TNFα^high^ cells and elevated mean fluorescence intensity (MFI) of TNFα (Figure 5E), a trend similarly observed in the 4T1 xenograft model (Figure 5E). In addition, multiplex immunofluorescence analysis revealed that HSPA1B⁺ immune cells exhibited higher TNFα expression intensity (Figure S5E). However, this pattern was absent in macrophages (Figure S5D), potentially due to their uniformly elevated stress signature (Figure 3J,I), which may obscure TNFα‐associated heterogeneity. Collectively, our analyses construct a cell–cell interaction network in stress‐programmed and point to TNFα as the key communication “hub” shaping the stress‐programmed ecosystem.
TNFα–TNFR2+ Treg Axis Underpins Adverse Treatment Response
2.6
To elucidate how TNFα, a pleiotropic cytokine with both pro‐ and antitumor functions [46, 47, 48], shapes the immune landscape of breast cancer patients undergoing neoadjuvant chemotherapy, we systematically investigated TNFα signaling dynamic across eight CMs spanning different treatment phases and clinical response categories (Figure 6A). Beyond, reaffirming CM8 as the principal source of TNFα signaling, we further observed that, after treatment, TNFα signaling strength in CM5 and CM1 was significantly higher in nonresponders than in responders (Figure 6A). Moreover, we found that the TNF–TNFRSF1B axis predominated within the TNFα signaling landscape. Notably, overall TNF signaling activity was stronger in nonresponders and showed a less pronounced decline following treatment compared to responders—a trend similarly observed for TNFRSF1B‐mediated signaling (Figure 6C). Additionally, in a neoadjuvant‐treated breast cancer cohort derived from a clinical trial comparing paclitaxel plus PD‐L1 inhibitor versus paclitaxel monotherapy [40], TNF expression in CD4⁺ T cells, CD8⁺ T cells, and B cells remained markedly elevated in nonresponders compared to responders, regardless of treatment phase (Figure 6B).
*Association of TNFα signaling in the tumor microenvironment with treatment efficacy and functional Treg expansion. (A) TNF signaling pathway networks inferred across co‐occurrence modules (CM1–CM8) at baseline and postsurgery in responders and nonresponders. Edge thickness represents overall communication strength between modules. (B) Violin plots showing TNF expression levels in CD4⁺ T cells, CD8⁺ T cells, and B cells from responders and nonresponders at baseline and surgery, based on the Zhang et al. dataset. Black dots represent the mean TNF expression. (C) Left: Relative contribution of individual ligand–receptor pairs to overall TNF signaling. Right: Bar plots comparing overall TNF signaling and TNFRSF1B‐specific signaling contributions across groups. (D) Boxplots showing TNFRSF1B signaling intensity across specific immune cell clusters in responders and nonresponders before and after treatment. Statistical significance was assessed by Wilcoxon test. (E) Scatter plot showing the alterations in TNFRSF1B signaling intensity for individual immune cell clusters in responders (R, y‐axis) and nonresponders (NR, x‐axis) between surgery and baseline. (F) UMAP plot showing the distribution of TNFRSF1B‐positive cells (pink) across the Treg landscape. (G) Violin plots comparing selected gene marker expression between TNFRSF1B‐positive and ‐negative Treg cells. Black dots represent the mean expression for each gene marker. (H) Violin plots showing TNFRSF1B expression in Tregs from responders and nonresponders in both the in‐house cohort and the Zhang et al. dataset. Statistical significance was assessed by Wilcoxon test. (I) Treg cells from three patients were treated with IgG, TNFα, or TNFα plus anti‐TNFR2 antibody, and CCR8 expression was measured by flow cytometry. Statistical significance was assessed by one‐way ANOVA. ns, not significant; *p < 0.05; **p < 0.01; ***p < 0.001; ***p < 0.0001.
To better pinpoint the cellular targets of TNFα signaling, we selected representative cell subsets from CM1, CM2, and CM5, three modules that showed the strongest correlation with CM8 and the most pronounced changes in TNFα signaling, to assess dynamic shifts in TNFα activity at the cell‐type level. We found that TNFRSF1B signaling intensity significantly declined postneoadjuvant therapy within the selected cellular modules in responders, whereas no substantial change was observed in nonresponders (Figure 6D). Notably, CD4T_c3_Treg exhibited the most pronounced decrease in TNFRSF1B activity among responders (Figure 6E), highlighting this subset as a potential key targets of TNFα in patients undergoing neoadjuvant chemotherapy.
TNFRSF1B is not only a key receptor promoting Treg proliferation and function, but also a crucial mechanism for maintaining immune tolerance in inflammatory and tumor settings [50, 51], making it a promising target for cancer therapy [52]. To investigate how TNFRSF1B influences Treg function in tumors, we defined TNFRSF1B‐positive and TNFRSF1B‐negative Treg populations using our in‐house Treg dataset (Figure 6F). Intriguingly, we found that TNFRSF1B‐positive Tregs exhibited a transcriptional profile more consistent with effector Tregs, characterized by elevated expression of canonical effector Treg markers [53, 54] (TNFRSF18 and TNFRSF4), core Treg identity genes (FOXP3 and IL2RA), and immune checkpoint molecules (TIGIT and CTLA4) (Figure 6G). To validate this finding, Treg cells isolated from three healthy donors were treated with TNFα or TNFα in combination with an anti‐TNFR2 antibody. We found that TNFα stimulation markedly upregulated the effector marker CCR8 [53] on Tregs, whereas blockade of TNFR2 effectively reversed this upregulation, indicating a critical role of the TNFα–TNFR2 axis in promoting Treg activation and function (Figure 6I). Additionally, we observed that TNFRSF1B expression on Tregs was significantly higher in nonresponders compared to responders, a trend consistently observed in both our in‐house dataset and the previously referenced public cohort [40] (Figure 6H).
Combining TNFR2 Blockade With PTX Offers an Effective Therapeutic Strategy
2.7
Our systematic analyses highlighted the abundant enrichment of stress‐programmed cell states and TNFα–TNFR2^+^ Treg axis in nonresponders, which may be a key contributor to clinical benefit from neoadjuvant chemotherapy. Although TNFR2 blockade confers antitumor immunity [55, 56], the pleiotropic and context‐dependent role of TNFα in breast cancer [57] obscures the mechanistic clarity of this pathway under neoadjuvant therapy. Based on this context and our preceding analyses, we were prompted to explore the potential of the TNFα–TNFR2⁺ Treg axis as a sensitization target for neoadjuvant chemotherapy in breast cancer patients.
To this end, we first employed the 4T1 murine breast cancer model by orthotopically implanting tumor cells into the mammary fat pad of BALB/c mice. Using anti‐TNFα (Bio X Cell, BE0058), anti‐TNFR2 (Bio X Cell, BE0247), and paclitaxel (PTX; MCE, B0015), we established seven treatment arms in the 4T1 murine breast cancer model: (1) control (IgG), (2) PTX alone, (3) anti‐TNFα alone, (4) anti‐TNFR2 alone, (5) PTX + anti‐TNFα, (6) PTX + anti‐TNFR2, and (7) PTX + anti‐TNFα + anti‐TNFR2 (Figure 7A). Importantly, while all three agents—PTX, anti‐TNFα, and anti‐TNFR2—exhibited significant antitumor effects when administered individually, anti‐TNFα demonstrated efficacy comparable to PTX, whereas anti‐TNFR2 monotherapy achieved superior outcomes; moreover, the combination of PTX and anti‐TNFR2 produced the most pronounced tumor suppression, even outperforming the triple‐agent regimen (Figure 7B and Figure S6A,B). There were no significant differences in body weight among mice across treatment arms, indicating the tolerability and safety of the therapeutic regimens (Figure S6C). We also observed that the proportion of TNFR2⁺ Tregs tended to increase following PTX monotherapy, albeit without statistical significance (Figure 7C,D). When used individually, anti‐TNFR2 exhibited the strongest capacity to deplete TNFR2⁺ Tregs, surpassing the limited effect of anti‐TNFα, and this depletion was similarly observed when these regimens were combined with PTX (Figure 7C,D). Focusing on the IgG, PTX, anti‐TNFR2, and PTX plus anti‐TNFR2 groups, we observed that PTX monotherapy significantly increased the frequency of PD‐1⁺ Tregs, whereas both anti‐TNFR2 alone and its combination with PTX led to marked reductions in this subset (Figure 7E,F). To further explore the clinical relevance of these findings, we employed a coculture system of patient‐derived organoids (PDOs) and paired PBMCs to evaluate the efficacy of PTX combined with anti‐TNFR2 therapy. We found that the combination of PTX and anti‐TNFR2 markedly enhanced tumor cell killing compared with either treatment alone (Figure 7G,H and Figure S6E). Thus, the depletion of TNFR2⁺ Tregs by combining PTX and anti‐TNFR2 represents a promising intervention strategy for breast cancer patients undergoing neoadjuvant chemotherapy. Furthermore, we observed that while inhibiting the TNFα/TNFR2 pathway effective, it specifically upregulates immunosuppressive signals like TGF‐β and IL‐10, but not PD‐1 (PDCD1) and PDL1 (CD274) (Figure S6D). These observations suggest that further investigations are merited to determine whether the combination therapy can enhance neoadjuvant chemotherapy outcomes and confer sustained survival advantages.
*Blockade of TNFR2 combined with PTX reduces breast cancer burden. (A) Schematic of the 4T1 tumor implantation and treatment regimen in BALB/c mice. Mice were randomized on day 7 and treated with paclitaxel (PTX, 200 µg) intraperitoneally on days 10, 13, and 17. Antibodies targeting TNFR2 and/or TNFα (200 µg each) were administered concurrently with PTX. (B) Left: Tumor growth curves across different treatment groups. Right: Pairwise comparisons of tumor volume on day 20, shown as color‐coded heatmaps of mean value of tumor volumes and corresponding –log10(p) values (t‐test). Specifically, blue denotes a smaller mean in the column group compared to the row group, and red denotes a smaller mean in the row group. (C) Flow cytometry plots showing TNFR2 expression on tumor‐infiltrating FOXP3⁺CD25⁺CD4⁺ cells of 4T1 implanted tumor under each treatment. (D) Quantification of TNFR2⁺ Treg proportions across groups. Left: Boxplot showing the proportions of TNFR2⁺ Tregs under all treatment conditions. Right: Heatmaps of pairwise comparisons of the proportions of TNFR2⁺ Tregs similar to (B). (E) Representative immunofluorescent images showing FOXP3 (green), PD‐1 (red), and nucleus (blue). White arrows mark FOXP3+PD‐1+ cells. Scale bars, 50 mm. (F) Quantification of FOXP3⁺PD‐1⁺ Tregs from immunofluorescence analysis. Statistical comparisons were performed using unpaired t‐test. (G) Representative images of patient‐derived organoids (PDOs) after treatment with IgG, paclitaxel (PTX), anti‐TNFR2 antibody, or the combination of PTX and anti‐TNFR2. Organoids were stained with live (green) and dead (red) cell markers and imaged by fluorescence microscopy. (H) Representative flow cytometry plots and quantification of dead tumor cells in PDOs from three patients under indicated treatments. Statistical comparisons were performed using one‐way ANOVA ns, not significant; *p < 0.05; **p < 0.01; ***p < 0.001; ***p < 0.0001.
Discussion
3
Our study delineates an immunological mechanism contributing to neoadjuvant chemotherapy response heterogeneity in breast cancer, centering on the enrichment of stress‐programmed immune cell states. Through an integrated high‐resolution single‐cell profiling, we uncovered that nonresponders harbor a conserved transcriptional stress program across B cells, macrophages, and T cells. Although stress‐programmed states in B cells, T cells, and NK cells have been identified in independent studies [17, 18, 19], these investigations were largely confined to individual immune cell types and did not systematically explore how such states interact within the broader immune microenvironment. Moreover, prior work primarily only established correlations with clinical outcomes, leaving the underlying mechanisms poorly understood. These limitations motivated us to conduct a more in‐depth examination.
In this study, we demonstrate for the first time that stress‐programmed immune cells are not isolated entities but rather coexist and interact across tumors, collectively forming a coordinated multicellular module. Furthermore, we mechanistically elucidate that this stress‐programmed module sustains chronic TNFα secretion, which in turn signals through TNFR2 on Tregs to maintain their immunosuppressive effector functions [58, 59]—a process that likely contributes to adverse clinical outcomes. Together, our findings extend previous observations by providing a more comprehensive understanding of stress‐programmed immune states and highlight the TNFα–TNFR2⁺ Treg axis as a potential therapeutic target in breast cancer. Although TNFα has traditionally been considered a key mediator of antitumor immunity, recent studies have highlighted its central role in promoting tumor stemness and immune evasion through both tumor‐intrinsic and immunological mechanisms [57, 60, 61, 62]. The net effect of TNFα on the TME may heavily rely on the location, timing, and duration of exposure. Depending on the tissue context, tumor progression status, and neoadjuvant chemotherapy regimen, these parameters may vary substantially across breast cancer patients [63]. In this context, prolonged and chronic TNFα stimulation is particularly relevant in an immunosuppressive tumor microenvironment in patients undergoing neoadjuvant chemotherapy, as such patients often present with more advanced disease stages and extensive tumor burden [64].
In early single‐cell research, stress‐related gene expression was often dismissed as a technical artifact [65], underscoring a long‐overlooked need for their systematic investigation. In our data, stress‐programmed cellular modules represent a distinct immune ecosystem enriched in tumor core regions, driving widespread immunoediting of the TME through sustained TNFα signaling, which is noteworthy. Although our study advances previous work by revealing the coordinated nature of stress‐programmed cell states and their downstream target subsets, the mechanisms underlying their emergence remain poorly understood. We speculate that tumor cells [66, 67], as key regulators of the TME, may orchestrate a broad reprogramming of tumor‐infiltrating lymphocytes through direct or indirect mechanisms, driving immune cells into a stress‐programmed state and impairing their effector functions—reflecting a critical and unmet need for further investigation.
Despite the lack of direct mechanistic evidence, our study hints at a role for stress‐programmed cellular modules in promoting immune evasion in breast cancer by fostering the expansion of TNFR2^+^ Treg, accentuating their potential relevance for clinical translation. Accordingly, blockade of the TNFα–TNFR2^+^ Treg axis effectively suppressed breast tumor growth across multiple treatment arms, with the combination of anti‐TNFR2 and PTX yielding the most pronounced tumor regression. Our strategy centers on the role of stress‐programmed cellular modules within the immune ecosystem and aims to eliminate their downstream immunosuppressive effects, thereby securing a net gain in antitumor immunity. Additionally, stress‐programmed immune cell states have been observed across multiple cancer types [17, 18, 19]. The extent to which the influence of stress‐programmed immune cell states on the tumor microenvironment can be extrapolated to other malignancies with distinct immunological contexts remains a key question for future research.
In summary, through comprehensive analysis of well‐controlled neoadjuvant‐treated patient samples combined with in vivo validation in preclinical models, we (1) revealed that nonresponders are characterized by a high abundance of stress‐programmed immune cell states; (2) uncovered the coordinated multicellular nature of these stress‐responsive programs across the tumor immune microenvironment; and (3) supported TNFα–TNFR2⁺ Treg targeting as a strategy to boost neoadjuvant chemotherapy response. Future investigations should aim to elucidate the underlying mechanisms driving the emergence of stress‐programmed immune cell states in tumors and translate these insights into therapeutic strategies to advance the neoadjuvant treatment landscape for breast cancer.
Limitations of the Study
3.1
Despite providing important insights into stress‐programmed immune niches and the TNFα–TNFR2 axis in breast cancer therapy resistance, our study has several limitations. The relatively small sample size of our neoadjuvant‐treated patient cohort may limit the generalizability of our findings across broader patient populations and breast cancer subtypes. Although our integrated analyses highlight the TNFα–TNFR2 axis as a central communication hub driving immune suppression, the mechanistic underpinnings of how stress‐programmed cell states emerge and propagate within the tumor microenvironment remain incompletely understood. Our conclusions are primarily supported by correlative data from single‐cell transcriptomics and in vivo tumor models; therefore, further mechanistic validation using genetic or pharmacological perturbation in lineage‐specific settings is necessary.
Additionally, while we demonstrated the therapeutic potential of TNFR2 blockade in murine models, preclinical systems cannot fully recapitulate the cellular complexity and immune heterogeneity of human breast tumors. Translation into clinical practice will require cautious assessment of safety, especially considering the pleiotropic roles of TNFα in immune regulation. Future studies incorporating longitudinal clinical samples will be essential to validate and extend our findings toward therapeutic application.
Experimental Section
4
Patient Cohorts and Sample Collection
4.1
We collected core needle biopsy samples prior to treatment and surgical specimens post‐treatment from breast cancer patients who received the taxane and platinum based neoadjuvant chemotherapy regimen at Fudan University Shanghai Cancer Center between 2019 and 2020. The patient received 4–6 cycles (21 days per cycle) of taxane and platinum based neoadjuvant therapy (Paclitaxel + Carboplatin‐based chemotherapy). We collected a total of 11 breast cancer tissue samples from 7 breast cancer patients, including paired core needle biopsy samples and surgical specimens from 4 of these patients. Additional tumor specimens were collected for flow cytometry and mIHC. All tumor samples were collected in strict adherence to clinical practice guidelines. The clinical information of all patients is presented in Table S1.
Single‐Cell RNA‐seq Library Preparation and Sequencing
4.2
Single‐cell RNA‐seq libraries were prepared using the Chromium Single Cell 3’ Reagent Kit v3.1 (10x Genomics) according to the manufacturer's protocol. Briefly, single cells were washed once with Phosphate‐Buffered Saline (PBS) containing 0.04% bovine serum albumin (BSA), and then resuspended in PBS with 0.04% BSA to reach a final concentration of 500–1200 cells/mL, determined by the Rigel S2 cell counter (Countstar). Around 10 000 cells were encapsulated into droplets to form nanoliter‐scale Gel Beads in Emulsion (GEMs). Reverse transcription was performed in a C1000 Touch Thermal Cycler (Bio‐Rad) under the following conditions: 53°C for 45 min, 85°C for 5 min, followed by a hold at 4°C. After reverse transcription and cell barcoding, the emulsions were broken, and cDNA was purified using Cleanup Mix with DynaBeads and SPRIselect reagent (Thermo Fisher Scientific), followed by PCR amplification.
For RNA‐seq library preparation, the amplified cDNA was fragmented, end‐repaired, and underwent double‐sided size selection. PCR amplification was performed using sample‐specific indexing primers, followed by an additional round of size selection. The libraries were then purified and evaluated for quality. Single‐cell RNA libraries were sequenced using an Illumina NovaSeq 6000 platform with 150 bp paired‐end reads.
Multiplex Immunofluorescence
4.3
Multiplex immunofluorescence staining was performed on formalin‐fixed, paraffin‐embedded (FFPE) sections of human breast cancer tissues and 4T1 implanted tumor following the manufacturer's instruction of Absin 7‐color IHC kit (cat abs50037, Absin). Tissue sections with a thickness of 4 µm were obtained from paraffin‐embedded blocks and affixed to glass slides. Deparaffinization was carried out using xylene, followed by rehydration through a descending ethanol gradient (100%, 95%, and 70%). Antigen retrieval was achieved by heating the slides in citric acid buffer (pH 9.0, Panovue) at boiling temperature in an oven for 30 min. After a 30‐min preincubation with blocking buffer at room temperature, the sections were sequentially incubated with primary antibodies. For human breast cancer sections, antihuman CD4 (Abcam, ab288724), antihuman CD8 (Abcam, ab217344), antihuman CD20 (Abcam, ab78237), antihuman CD68 (Abcam, ab303565), antihuman CD45(Abcam, ab10558), antihuman TNFa (Abcam, ab220210), and antihuman HSPA1B (Abclonal, A19317) antibodies were used. For 4T1 implanted tumor sections, antimouse FOXP3 (Abcam, ab215206) and antimouse PD‐1 (Abcam, ab214421) antibodies were applied. Following TBST washes, the sections were incubated with horseradish peroxidase (HRP)‐conjugated secondary antibodies, followed by application of the tyramide signal amplification (TSA) working solution. In addition, between each round of antibody incubation, antigen retrieval was performed by microwave heating in a water bath for 20 min. After labeling all the above antigens, nuclei were stained with DAPI. The stained slides were then scanned using the Mantra System (PerkinElmer). Quantitative analysis was performed using Halo software.
Flow Cytometry
4.4
The fresh human breast cancer tissue and fresh 4T1 implanted tumor tissue was cut into small pieces on ice, digested with digestion buffer (cat abs9482, Absin) on a shaker at 37°C for 2 h, and then filtered through a strainer to prepare a single‐cell suspension as per the manufacturer's instructions. For human breast cancer tissues, We prepared three different flow cytometry panels to detect TNFα signaling in various immune cell populations. After incubation with Fixable Viability Dye (BioLegend) and surface protein antibodies, including B cells panel:CD45‐AF700 (BD Pharmingen), CD3‐FITC (Absin), CD19‐APC (BD Pharmingen); T cells panel: CD45‐AF700 (BD Pharmingen), CD3‐FITC (Absin), CD4‐APC (Absin), CD8‐percp‐cy5.5 (Absin); Myeloid panel: CD45‐AF700 (BD Pharmingen), CD3‐FITC (Absin), CD11B‐percp‐cy5.5 (BD Pharmingen), CD68‐FITC (BD Pharmingen), the cells were fixed and permeabilized with a Fixation/Permeabilization Solution Kit (BD Biosciences). Then these fixed cells were subsequently incubated with HSP70‐PE (Novus) and TNFa‐PE‐CY7 (BD Pharmingen) in 13 perm buffer for 1 h at 4 C. For Treg cells and patient‐derived organoids (PDOs), the following antibodies were used for flow cytometry: CD4–PB450 (BioLegend), CD25–PerCP‐Cy7 (BD Biosciences), FOXP3–PE (eBioscience), CD45–AF700 (BioLegend), and EPCAM–FITC (BioLegend). For 4T1 cc implanted tumor tissue, two flow cytometry panels were used. Panel 1 included CD45‐AF700 (BD Pharmingen), CD3‐RB705 (BD Pharmingen), CD4‐FITC (BD Pharmingen), CD25‐PE‐Cy7 (BD Pharmingen), FOXP3‐AF647 (BD Pharmingen), and TNFR2‐PE (BioLegend). Panel 2 included CD45‐AF700 (BD Pharmingen), CD3‐RB705 (BD Pharmingen), CD4‐FITC (BD Pharmingen)∖CD8‐FITC (BD Pharmingen)∖B220‐FITC (BD Pharmingen), HSP70‐PE (Novus), and TNFα‐APC (BD Pharmingen). FOXP3, TNFR2, HSP70‐PE, and TNF‐α were detected as intracellular markers, using the same staining protocol and sequence as previously described. The fluorescence intensity was determined using a Beckman Coulter Cytoflex and downstream analysis was performed by FlowJo. The gating strategies for different panels are illustrated in (Figures S7 and S8).
Single‐Cell RNA‐seq Data Preprocessing and Quality Control
4.5
FASTQ files were processed using Cellranger (v3.0.0) to obtain 10x count matrix for downstream analysis with the GENCODE v38 human genome as the reference. In this version of Cellranger, only exonic reads would be counted to the final UMI count.
The matrix of read counts per gene per sample was further analyzed by the Seurat (v 4.3.0) package. First, to obtain high‐quality cells, we used DoubletFinder [68] (v2.0.3) to detect and remove droplets containing more than one cell. Second, cells that satisfy any of the following criteria were excluded from downstream analysis: (1) < 500 expressed genes or > 10 000 expressed genes, (2) >100 000 nCount_RNA, (3) > 20% UMIs of mitochondria genes, (4) > 50% UMIs of ribosome genes, and (5) > 5% UMIs of hemoglobin‐associated genes. Finally, we used decontX [69] (v 1.0.0) to remove cells with potential RNA contamination. Cells with a contamination score greater than 20% were defined as contaminated and excluded from downstream analyses. After these quality control steps, we further removed cells expressing multiple cell markers based on existing biological knowledge to obtain the final clean matrix. At last, we preserved a total of 68 534 cells.
Single‐Cell RNA‐seq Data Integration and Annotation
4.6
After the preprocessing and quality control of each sample as described above, we selected the top 2000 highly variable genes for PCA‐based dimension reduction after normalizing the gene expression matrix using the NormalizeData function. Harmony package was used to remove sample‐level batch effects after calculating the principal components and before conducting UAMP‐based dimension reduction. We considered the batch effect correction successful when no cell clusters were dominated by a single sample.
For cell type annotation, we first identified 9 major cell types at a resolution of 0.6 based on the expression distribution of well‐established marker genes. In brief, epithelial cells are marked by EPCAM and KRT18, smooth muscle cells are marked by TAGLN and ACTA2, Endothelial cells are marked by VWF and PLVAP, Fibroblasts are marked by COL1A1 and COL1A2, T/NK cells are marked by CD3D and NKG7, B cells are marked by CD19 and CD79A, Myeloid cells are marked by AIF1 and CD163, Plasma are marked by JCHAIN and MZB1 and Mast cells are marked by KIT. Next, for each major cell type, we subsetted and reprocessed the count matrices using the workflow mentioned above. In this section, clustree was used to identify appropriate resolution for each major cell type. We then identified the subcluster through an extensive search of the literature after examining the expression distribution of gene markers across cell clusters and named them by selection of one specific signature gene, taking into account both the observed expression patterns and established knowledge.
Calculation of Ro/e Value
4.7
The Ro/e value, as defined by Zhang et al. [6]., measures the ratio of the observed to expected cell numbers for a particular combination of cell clusters and treatment stages. To assess the distribution patterns of distinct cell subpopulations before and after treatment, as well as between responders and nonresponders patients, we computed the Ro/e value for each subpopulation across different conditions. A cluster was considered enriched in a specific condition when Ro/e exceeded 1.
Calculation of Signature Score at the Single‐Cell Level
4.8
R package AUCell (v1.22.0) were used to calculate the signature score of a specific gene set. AUCell is a ranking‐based tool for assessing gene set activity at the single‐cell level, independent of gene expression units or normalization methods. We first ranks genes for each cell using the AUCell_buildRanking function, creating a ranked expression matrix that orders genes according to their expression levels within each cell. Then, the AUCell_calcAUC function is applied to compute the Area Under the Curve (AUC) scores for each gene set. Detailed gene sets used in this study were listed in Table S2.
Definition of HSPA1B Score in Immune Cell Subsets and Its Prognostic Significance
4.9
For precise evaluation of HSPA1B score in immune cells while minimizing confounding effects from cellular subtypes, its mRNA expression were normalized across canonical immune cell populations. For example, to define the HSPA1B score in CD8⁺ T cells, HSPA1B expression was normalized by dividing it by the product of CD8A and CD8B expression levels, thereby accounting for the abundance of the CD8⁺ T cell population. Based on this, we defined the CD4T HSPA1B score, CD8T HSPA1B score, B HSPA1B score, and Macro HSPA1B score. This normalization and denoising approach was adapted from Eynav et al. [70] and Gao et al. [53].
To evaluate the prognostic relevance of the HSPA1B score within specific immune cell subsets, TCGA breast cancer samples were classified into high‐ and low‐expression groups for each of the four scores—CD4T HSPA1B score, CD8T HSPA1B score, B cell HSPA1B score, and Macro HSPA1B score—using the median value as the cutoff. This stratification enabled survival analyses to assess associations between HSPA1B expression and clinical outcomes. Kaplan–Meier survival curves were generated using the survival package (v3.5.5) and survminer package (v0.4.9).
Definition of Stress‐Programmed Score
4.10
A similar strategy to that used for defining the HSPA1B score in immune cell subsets was applied to evaluate the overall Stress‐Programmed score. To more precisely assess the global stress response, the stress score was calculated by ssGSEA based on the expression of HSPA1A, HSPA1B, DNAJB1, DNAJB6, and DNAJA1, and normalized across canonical immune cell populations.
Monocle3 Trajectory
4.11
Monocle3 was applied to reconstruct the cellular differentiation trajectories of B cells and myeloid cell subsets. Specifically, cell clusters were first grouped into well‐separated partitions using the cluster_cells function, and a principal graph was then fitted within each partition using the learn_graph function with default parameters. The resulting principal graph, visualized as “skeleton lines” on the UMAP plot, represents inferred differentiation paths. Based on prior biological knowledge, the clusters enriched in B_c0 and B_c5 were designated as the root for the B cell trajectory, while those enriched in Mono_c0 and Mono_c1 were used as the root for the myeloid cell trajectory.
Correlation and Clustering Analysis of Cell Subpopulations
4.12
To calculate the correlation coefficient between cell subpopulations, we first assessed the proportion of each subcluster with the 4 major immune cell types across individual patients. Subsequently, we utilized the Hmisc (v5.1.0) package to calculate the correlations between these cell subpopulations, aiming to distinguish different immune cell niches. Subsequent clustering and visualization were carried out using the ComplexHeatmap (v2.16.0) package.
TCGA RNA‐seq Data Deconvolution
4.13
TCGA RNA‐seq expression matrix and clinical information were downloaded and processed using the TCGAbiolink (v 2.28.4) package. To precisely map annotated cell subpopulations from single‐cell data onto bulk RNA‐seq data, we employed the BayesPrism [35] (v2.2.2) package for the deconvolution of cell subtype and gene expression profiles using count matrix from TCGA BRCA cohorts. In this practice, to remove unnecessary noise, we focused on the deconvolution of immune subpopulations in the bulk data. Therefore, we subsetted immune cell subpopulations from the integrated single‐cell data and then performed deconvolution using the BayesPrism package with default parameters.
Cellular Communication Analysis
4.14
To identify potential ligands that drive the phenotype of selected cell clusters with NicheNet [42] (v2.1.0), we first calculated the differentially expressed genes of each cluster. In this practice, the top 400 DEGs with p_val_adj ≤ 0.05 and absolute value of avg_log2FC > 0.25 were selected for interest. All expressed genes in the integrated single‐cell data were used as background genes. Expressed ligands from sender cells and receptors form receiver cells were identified, and only significantly expressed pairs were included in the analysis. Using the NicheNet R package, we calculated ligand activity scores to evaluate the influence of ligands on the expression of target genes in receiver cells. Then, the top 30 predicted ligands ordered by ligand activities were used to find the mutual ligand in the stress related immune cell clusters. The analysis of specific cellular pathways in subsequent steps was conducted using the CellChat [41] (v1.6.1) R package.
Spatial Transcriptomics Data Analysis
4.15
Spaceranger (v 2.1.1) was used to align the raw data to the GRCh38 reference genome and the downstream analysis was conducted by Seurat package. Initial quality control filtering spots with fewer than 200 detected genes, and genes with fewer than 10 read counts or those expressed in fewer than three spots were removed. Normalization and scaling were performed using Seurat's standard workflow, followed by dimension reduction with PCA and clustering to identify cell populations. Spatial distribution of clusters was visualized using UMAP and spatial features plots.
To investigate the spatial distribution dynamics of immune cell subsets, we utilized UCell package to assign scores to these populations across spatial spots. Notably, we implemented a filtering process based on the classification of broader cell types to ensure accuracy. For instance, spots were excluded if they displayed a stress‐responsive CD4T score greater than 0, while the overall CD4T cell score remained below 0.1. Each spot containing the corresponding cell type was evaluated for its expression using established and literature‐referenced cell markers (Table S2). Each 10x Visium sample was analyzed separately, with each sample undergoing the aforementioned workflow individually.
Identification of Tumor Core Regions and Calculation of Spatial Distance
4.16
Tumor core regions and boundary regions were identified using spatial transcriptomics data with the Cottrazm [38] (v 0.1.1) package. Malignant spots were defined based on the highest CNV scores, indicating regions with substantial copy number variations and the highest proportion of malignant cells. These regions were further validated by mapping the spatial distribution of CNV‐enriched clusters.
Spatial distance calculations were performed using the semla [71] (v 1.1.6) package to determine the proximity of each spatial spot relative to the tumor core, with the tumor edge region serving as the boundary. This was achieved by calculating the Euclidean distance between each spot and the center of the identified malignant core. Distances were measured in relation to the tumor edge, allowing us to assess spatial gradients in gene expression and cell type distribution from the edge toward the tumor core.
RNAseq Analysis
4.17
Raw RNA‐seq reads in FASTQ format were subjected to quality control and adapter trimming using fastp [72] (v0.23.2). Clean reads were aligned to the mouse reference genome (mm10) using STAR [73] (v2.7.10a) with default parameters. Gene‐level read counts were quantified using FeatureCounts [74] (v2.0.3) from the Subread package based on the Gencode M25 annotation. Gene expression levels were subsequently normalized as TPM (Transcripts Per Million) for visualization.
Cell Cultures
4.18
Cell line 4T1 were obtained from American Type Culture Collection (ATCC) (Manassas, VA). 4T1 were cultured in high‐glucose Dulbecco's Modified Eagle's Medium (DMEM) (Cytiva, USA) supplemented with 10% Fetal Bovine Serum (FBS) (FBS, ExCell, China). All cells were incubated in a constant‐temperature incubator at 37°C with 5% CO_2_.
Human Tregs Isolate and Culture
4.19
Peripheral blood mononuclear cells (PBMCs) were isolated from peripheral blood of health donors using Ficoll‐Paque PLUS (GE Healthcare) density gradient centrifugation. CD4⁺CD25⁺ regulatory T cells (Tregs) were purified from PBMCs using the CD4⁺CD25⁺ Regulatory T Cell Isolation Kit, human (Miltenyi Biotec, Cat. no. 130‐091‐301) according to the manufacturer's instructions.
Purified Tregs were cultured in RPMI 1640 medium supplemented with 10% fetal bovine serum (FBS), 1% penicillin–streptomycin, 50 µM 2‐mercaptoethanol, 100 U/mL recombinant human IL‐2 (STEMCELL Technologies, Cat. no. 78067), and 5 ng/mL recombinant human TGF‐β1 (STEMCELL Technologies, Cat. no. 78036) at 37°C with 5% CO_2_. For activation, cells were stimulated with anti‐CD3/CD28‐coated beads (25 µL/mL; STEMCELL Technologies, Cat. no. 10981) according to the manufacturer's instructions. Where indicated, recombinant human TNFα (UA Bioscience, Cat. UA040018) and anti‐TNFR2 antibody (Thermo Fisher, Cat. MA1‐24723) were added at the indicated concentrations prior to functional assays. After treatment, cells were collected for flow cytometry analysis of CCR8 expression.
In Vivo Mouse Models
4.20
All animal procedures were performed in compliance with institutional ethical standards and were approved by the Animal Ethics Committee of Fudan University Shanghai Cancer Center (protocol number: FUSCC‐IACUC‐2024645). Mice were maintained in a specific pathogen‐free (SPF) environment with controlled temperature (18–23°C), humidity (40%–60%), and a 12‐h light/dark cycle. For each study, animals were randomly allocated to different treatment groups; no blinding was applied.
Female BALB/c mice (6–8 weeks old) were used for all in vivo experiments. 4T1 cells were resuspended in D‐PBS (Gibco) and combined with Matrigel Matrix (Corning, 356234) at a 1:1 ratio, yielding a final concentration of 5 × 10^6^ cells/mL. A total of 5 × 10^5^ cells were then subcutaneously injected into the fourth mammary fat pad of each mouse. Two animal cohorts were constructed for experimental purposes. In the first cohort, the mice were intravenously injected with PTX (200 µg per mouse) on day 10, 13, 17, 20, and 24, as illustrated in Figure 2D. In the second cohort, mice were randomly assigned to seven treatment groups: (1) Control (IgG), (2) PTX alone, (3) anti‐TNFα, (4) anti‐TNFR2, (5) PTX + anti‐TNFα, (6) PTX + anti‐TNFR2, and (7) PTX + anti‐TNFα + anti‐TNFR2. PTX (200 µg per mouse) was administered intravenously via the tail vein, while anti‐TNFα, anti‐TNFR2, and IgG (each at 200 µg per mouse) were delivered via intraperitoneal injection. All treatments were performed on days 10, 13, and 17 after tumor cell inoculation, as illustrated in Figure 7A. When the tumor was palpable, its long diameter (L) and short diameter (W) were measured every 3 days with a calliper, the tumor volume was determined using the formula: 1/2 * L * W^2^. The mice were euthanized by CO_2_ inhalation and the tumors were weighed.
Coculture of Patient‐Derived Organoids (PDOs) and PBMCs
4.21
Paired PDOs and PBMCs were established from the same patient to preserve the autologous tumor‐immune context. PDOs were generated from freshly resected tumor tissues as previously described, while PBMCs were isolated from matched peripheral blood samples using Ficoll‐Paque PLUS (GE Healthcare) density gradient centrifugation.
Mature PDOs were dissociated into small clusters using TrypLE Express (Thermo Fisher Scientific) and re‐embedded in Matrigel (Corning, Cat. 356 231) in 48‐well plates. PBMCs were activated for 24 h with anti‐CD3/CD28‐coated beads (25 µL/mL; STEMCELL Technologies, Cat. 10 981) in RPMI 1640 medium supplemented with 10% FBS, 1% penicillin–streptomycin, 100 U/mL recombinant human IL‐2 (STEMCELL Technologies, Cat. 78 067), and 5 ng/mL recombinant human TGF‐β1 (STEMCELL Technologies, Cat. 78 036).
Activated PBMCs were then cocultured with autologous PDOs at an effector‐to‐target (E:T) ratio of 5:1 in a mixed medium composed of 50% RPMI 1640 and 50% organoid basal medium. Where indicated, paclitaxel (PTX), and anti‐TNFR2 antibody (Thermo Fisher, Cat. MA1‐24723) were added. Cocultures were maintained for 1 week at 37°C with 5% CO_2_. After coculture, PDOs were collected for live/dead cell staining using Calcein‐AM (green, live cells) and Ethidium homodimer‐1 (red, dead cells) and dissociated for flow cytometry analysis.
Statistical Analysis
4.22
Detailed data analysis procedures can be found in the figure legends and the Experimental Section. The Wilcoxon test assessed statistical significance, with survival analysis utilizing the Kaplan–Meier method. Correlation analysis employed Pearson correlation, and a p‐value of less than 0.05 was considered statistically significant.
Author Contributions
Z.S. and X.W. conducted the research, analyzed the single‐cell data, and wrote the manuscript. H.Z. and L.Y. aided with animal experiment management and tumor harvesting in mice. B.X. and R. W. provided conceptual guidance and contributed to the overall logic of the manuscript. J.W., S.C., and Z.S. initiated the study, organized, designed, and revised the manuscript. Q.Z., Z.Z., X.L., X.X., Y.Z., X.L., X.Q., J.Y., X.L., and J.D. participated in data collection, sample processing, or manuscript preparation. Z.S. provided technical support. All authors have read and approved the final version of the manuscript.
Funding
This work was supported by grants from the National Natural Science Foundation of China (82373373) and National Natural Science Foundation of China (82203097).
Conflicts of Interest
The authors declare no conflicts of interest.
Supporting information
Supporting File: advs73624‐sup‐0001‐SuppMat.docx.
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1T. F. Gajewski , H. Schreiber , and Y. X. Fu , “Innate and Adaptive Immune Cells in the Tumor Microenvironment,” Nature Immunology 14, no. 10 (2013): 1014–1022, 10.1038/ni.2703.24048123 PMC 4118725 · doi ↗ · pubmed ↗
- 2D. S. Chen and I. Mellman , “Elements of Cancer Immunity and the Cancer–Immune Set Point,” Nature 541, no. 7637 (2017): 321–330, 10.1038/nature 21349.28102259 · doi ↗ · pubmed ↗
- 3M. Binnewies , E. W. Roberts , K. Kersten , et al., “Understanding the Tumor Immune Microenvironment (TIME) for Effective Therapy,” Nature Medicine 24, no. 5 (2018): 541–550, 10.1038/s 41591-018-0014-x.PMC 599882229686425 · doi ↗ · pubmed ↗
- 4Y. Zhang and Z. Zhang , “The History and Advances in Cancer Immunotherapy: Understanding the Characteristics of Tumor‐Infiltrating Immune Cells and Their Therapeutic Implications,” Cellular & Molecular Immunology 17, no. 8 (2020): 807–821, 10.1038/s 41423-020-0488-6.32612154 PMC 7395159 · doi ↗ · pubmed ↗
- 5E. Papalexi and R. Satija , “Single‐Cell RNA Sequencing to Explore Immune Cell Heterogeneity,” Nature Reviews Immunology 18, no. 1 (2018): 35–45, 10.1038/nri.2017.76.28787399 · doi ↗ · pubmed ↗
- 6L. Zhang , X. Yu , L. Zheng , et al., “Lineage Tracking Reveals Dynamic Relationships of T Cells in Colorectal Cancer,” Nature 564, no. 7735 (2018): 268–272, 10.1038/s 41586-018-0694-x.30479382 · doi ↗ · pubmed ↗
- 7E. Azizi , A. J. Carr , G. Plitas , et al., “Single‐Cell Map of Diverse Immune Phenotypes in the Breast Tumor Microenvironment,” Cell 174, no. 5 (2018): 1293–1308.e 1293, 10.1016/j.cell.2018.05.060.29961579 PMC 6348010 · doi ↗ · pubmed ↗
- 8D. Lambrechts , E. Wauters , B. Boeckx , et al., “Phenotype Molding of Stromal Cells in the Lung Tumor Microenvironment,” Nature Medicine 24, no. 8 (2018): 1277–1289, 10.1038/s 41591-018-0096-5.29988129 · doi ↗ · pubmed ↗