Cellular reprogramming during anti-PD-1 and chemotherapy treatment in early-stage primary hormone receptor-positive breast cancer

Abstract

The efficacy of immune checkpoint inhibitors combined with chemotherapy varies among breast cancer subtypes and is particularly less effective in hormone receptor-positive (HR + ) breast cancers. Here, we analyze pre-, on-, and post-treatment biopsies from 20 female patients with stage II-III HR+ breast cancer who participated in a clinical trial of neoadjuvant chemo-immunotherapy with nab-paclitaxel and pembrolizumab. Through single-nucleus RNA and ATAC sequencing of these tumor biopsies, we identified gene expression metaprograms (MPs) associated with differential therapy responses. Here we show that favorable responders exhibit increased activity in pathways related to tumor state transition, T cell effector functions, and pro-inflammatory macrophage states. Unfavorable responders demonstrate increased tumor estrogen signaling and immunosuppressive tumor-immune interactions. In this work, we highlight the interplay between tumor and microenvironmental cells in treatment naïve and exposed HR+ breast cancers and reveal that pivotal shifts in tumor cell, macrophage, and T cell states may mediate response to chemo-immunotherapy.

Fig. 1. The transcriptomic and epigenetic landscape of tumor microenvironment in primary hormone receptor-positive breast cancer at single cell resolution.

A The study design (Created in BioRender. Miler-jones, L. (2025) https://BioRender.com/fv15ktr). B Genomic and clinical overview of the hormone receptor-positive (HR+) breast cancer sample. Each column represents a tumor sample. Tumors are ordered by residual cancer burden (RCB) response (RCB 0–I or RCB II-III) and within each subgroup ordered by treatment arm. “Stage” denotes the breast cancer stage. “BluePrint” refers to the molecular subtype assessed by MammaPrint. The presence of bulk RNA-seq and bulk WES data is indicated by a black box for each tumor. The type of 10x Genomics assay used for each tumor is denoted by different colors. An additional four cycles of neoadjuvant Adriamycin/cyclophosphamide (AC) are marked by a black box for each patient. Somatic mutations in genes frequently mutated in breast cancer and common copy number alterations are displayed for each tumor. Sample indicators connect samples from the same patients. C UMAP representation of transcriptional (left) and epigenetic landscape (right) of sequenced cells. Top: UMAPs of broad cell-type annotation. Bottom: Transcriptionally based UMAPs displaying canonical marker gene expression. D Proportion of cellular compartments in each sample. Samples are colored by the RCB response. Source data are provided as a Source Data file.

Fig. 2. Tumor intrinsic gene signatures associated with combination therapy tesponse.

A Heatmap displaying pairwise Jaccard similarity indices among robust non-negative matrix factorization (rNMF) programs based on their top 50 genes. Programs are clustered into nine metaprograms (MPs); MPs, timepoint and patient information are labeled at the top. The central scatterplot illustrates Pearson correlation of rNMFs with cellular complexity. B Heatmap depicting gene membership within MPs, with rows as top representative genes, columns as MPs, and functional annotations on the right. C Distribution of ER-I-related MP7 and EMT-III-related MP11 signature level in tumor cells from pretreatment biopsies. P values from two-sided Wald tests on linear mixed model coefficients (patient as random effect; Methods), without multiple-comparison adjustment. Data presented as median with nested quantile ranges and boxes narrowing toward distribution tails to show extreme values. D Baseline comparison of the relative abundance of ER-I-related MP7 and EMT-III-related MP11 states in the tumor population between favorable responders (R; N(pt) = 3) and unfavorable responders (NR; N(pt) = 9), with p values from two-sided Mann–Whitney–Wilcoxon test. Data presented as median with interquartile range (first and third quartiles). E Dotplot of scaled transcription factor expression (color) and target region enrichment scores (dot size). Each row shows a gene regulatory network (GRN), named after the transcription factor and its target region; columns represent four tumor cell states. F PROGENy-inferred activity in ER-I (MP7) and EMT-III (MP11) tumor cell states. Barplots (left) show pathway activity scores relative to other tumor states. X-axis (right) shows gene weights in the p53 (top) and TGFβ (bottom) pathways; Y-axis shows t-values from differential expression comparing MP7 or MP11 to other states. G ER-I (MP7) and EMT-III (MP11) abundance in tumors with or without TP53 mutations from pretreatment biopsies. P values from two-sided Mann–Whitney–Wilcoxon test. Data presented as median with interquartile range (first and third quartiles). H ER-I-related MP7 and EMT-III-related MP11 gene signatures relative to the overall tumor signature in bulk RNA-Seq data from TCGA primary HR+ breast cancer. Patients grouped by TP53 mutation status, with significance tested by two-sided Mann–Whitney–Wilcoxon test. Data presented as median with interquartile range (first and third quartiles). Source data are provided as a Source Data file.

Fig. 3. Reprogramming of CD8+ T cells and their molecular signatures during monotherapy and combination therapy.

A Heatmap illustrating average gene expression from three metaprograms (MPs) among three distinct CD8 + T cell states. Rows represent gene names. Top: cell state, patient, and timepoint information. B Differential cytokine activity between cytotoxic and exhausted CD8 + T cell states, analyzed using the two-sided Mann–Whitney–Wilcoxon test. The Y-axis displays negative log10 p values; the X-axis shows t-values. C Boxplot depicting distribution of predictive IL-15 cytokine activity across three CD8 + T cell states. Cytotoxic, native, and exhausted CD8 + T cell sates were detected in 38, 38, and 35 samples, respectively. Each point represents a sample, with p values calculated using the two-sided Mann–Whitney–Wilcoxon test. Data presented as median with interquartile range (first and third quartiles). D Violin plots illustrating shift in distribution of three MP gene signatures from baseline to monotherapy in favorable and unfavorable responders with biopsies at both timepoints. Top: CD8 + T cells during chemotherapy vs. pretreatment CD8 + T cells. Bottom: CD8 + T cells during Pembrolizumab treatment vs. pretreatment CD8 + T cells. P values were calculated via two-sided Wald tests on linear mixed model coefficients (Methods: Differential gene signature analysis), with patient ID modeled as a random effect. No adjustment made for multiple comparison. E Violin plots displaying shift in distribution of three MP gene signatures from monotherapy to combination therapy in favorable and unfavorable responders with biopsies at both timepoints. P values calculated via two-sided Wald tests on linear mixed model coefficients (Methods: Differential gene signature analysis), with patient ID modeled as a random effect. No adjustment made for multiple comparison. F Violin plots showing shift in distribution of three MP gene signatures from combination therapy to post-combination therapy (with AC treatment) in favorable and unfavorable responders with biopsies at both timepoints. P values calculated via two-sided Wald tests on linear mixed model coefficients (Methods: Differential gene signature analysis), with patient ID modeled as a random effect. No adjustment made for multiple comparison. RCB residual cancer burden. Source data are provided as a Source Data file.

Fig. 4. Reprogramming of macrophages and their molecular signatures during monotherapy and combination therapy.

A Heatmap displaying average gene expression from six metaprograms (MPs) among five macrophage cell states. Rows represent gene names. Top: cell state, patient, and timepoint information. B Pairwise Pearson correlation among the six MP gene signatures. P values calculated with two-sided Student’s t-distribution. C Boxplot illustrating distribution of six MP gene signatures in macrophages during monotherapy, comparing 2802 macrophage cells from two favorable responders and 1990 macrophage cells from seven unfavorable responders. P values were calculated via Wald tests on linear mixed model coefficients (Methods: Differential gene signature analysis), with patient ID modeled as a random effect. No adjustment made for multiple comparison. Data presented as median with multiple nested quantile ranges (50–75%, 75–87.5%, 87.5–93.75%, 93.75–96.875%, and 96.875–100%), with boxes narrowing toward distribution tails to show extreme values. D Violin plots showing shift in distribution of six MP gene signatures from baseline to monotherapy in favorable and unfavorable responders with biopsies at both timepoints. Top: Macrophages during chemotherapy vs. pretreatment macrophages. Bottom: Macrophages during Pembrolizumab treatment vs. pretreatment macrophages. P values were calculated via two-sided Wald tests on linear mixed model coefficients (Methods: Differential gene signature analysis), with patient ID modeled as a random effect. No adjustment made for multiple comparison. E Violin plots depicting shift in distribution of six MP gene signatures from monotherapy to combination therapy in favorable and unfavorable responders with biopsies at both timepoints. P values calculated via two-sided Wald tests on linear mixed model coefficients (Methods: Differential gene signature analysis), with patient ID modeled as a random effect. No adjustment made for multiple comparison. F Violin plots illustrating shift in distribution of six MP gene signatures from combination therapy to post-combination therapy (with AC treatment) in favorable and unfavorable responders with biopsies at both timepoints. P values calculated via two-sided Wald tests on linear mixed model coefficients (Methods: Differential gene signature analysis), with patient ID modeled as a random effect. No adjustment made for multiple comparison. RCB residual cancer burden. Source data are provided as a Source Data file.

Fig. 5. Tumor cellular interactions with CD8+ T cells and macrophages.

A Rank plots displaying baseline ligand-receptor (L–R) interactions between tumors and CD8 + T cells, sorted by their MultiNicheNetR prioritization score. Top: Predominant L–R interactions in favorable responders. Bottom: Predominant L–R interactions in unfavorable responders. Annotations highlight L–R interactions involving immune checkpoint receptors on CD8 + T cells that rank within the top 50 for each group. (Top illustration was created in BioRender. Miler-jones, L. (2025) https://BioRender.com/e6aheln). B Boxplots illustrating the distribution of three L–R interaction metrics for tumors expressing ligands that interact with the CTLA-4 receptor on CD8 + T cells, comparing favorable and unfavorable responders. Left: Probability of producing both ligand and receptor; Middle: Probability of producing the ligand; Right: Probability of producing the receptor. Data were presented as median with interquartile range (first and third quartiles). C Rank plots displaying baseline ligand-receptor (L–R) interactions between tumors and macrophages, sorted by their MultiNicheNetR prioritization score. Top: Predominant L–R interactions in favorable responders. Bottom: Predominant L–R interactions in unfavorable responders. Annotations highlight the top five L–R interactions for each group. (Top illustration was created in BioRender. Miler-jones, L. (2025) https://BioRender.com/e6aheln). D Boxplots illustrating the distribution of three L–R interaction metrics for tumors expressing ligands that interact with the receptors on macrophage, comparing favorable and unfavorable responders. Left: Probability of producing both ligand and receptor; Middle: Probability of producing the ligand; Right: Probability of producing the receptor. Data are presented as median with interquartile range (first and third quartiles). E Heatmap illustrating the changes in tumor-CD8 + T cell L–R interactions from pretreatment (BL) to on-combination therapy (W7D1). The top five altered L–R interactions for each group are displayed. F Heatmap depicting the changes in tumor-macrophage L–R interactions from pretreatment (BL) to on-combination therapy (W7D1). The top five altered L–R interactions for each group are shown. RCB residual cancer burden. Source data are provided as a Source Data file.