Spatial and genomic profiling of residual breast cancer after neoadjuvant chemotherapy unveil divergent fates for each breast cancer subtype

Abstract

Residual cancer burden (RCB) is a strong prognostic marker after neoadjuvant chemotherapy (NAC) in breast cancer (BC), yet some BCs defy their predicted outcomes. Using single-cell spatial transcriptomics and genomic profiling, we investigate mechanisms underlying divergent fates of BCs with high RCB across subtypes. In triple-negative BC (TNBC), CXCL9+ macrophage-CD8⁺ T cell interactions via chemokines and interferon-gamma signaling promote favorable outcomes, while SPP1+ macrophage-cancer cell interactions driven by hypoxia signaling correlate with poor prognosis. In non-TNBC, the extent of basal-like cancer cells and their proximity to scarce immune cells are linked to prognosis. Additionally, tumor-intrinsic features-such as homologous recombination deficiency in hormone receptor (HR)-positive cancers and structural variations, including extrachromosomal ERBB2 DNA in human epidermal growth factor receptor 2 (HER2)-positive cancers-predict worse outcomes. This study highlights distinct genomic and microenvironmental strategies governing BC subtype-specific fates after NAC.

Keywords: breast cancer; extrachromosomal DNA; intrinsic subtype; neoadjuvant treatment; residual cancer burden; spatial transcriptomics; whole-genome sequencing.

Graphical abstract

Figure 1

Spatial profiling identifies distinct myeloid cell subtypes and spatial features associated with dismal prognosis in post-NAC residual TNBC (A) Flightpath plot showing InSituType cell identification with posterior probability. (B) Dotplot of signature marker expression by cell type. (C) Uniform manifold approximation and projection (UMAP) of reclustered myeloid cells. (D) Dotplot of signature markers for myeloid cells. (E) Milo cellular neighborhood differential abundance plots (false discovery rate [FDR] < 0.1) showing myeloid composition changes (left); beeswarm plot of log-fold changes in neighborhoods by subcluster (right). (F) Boxplots of CXCL9+/SPP1+ Mac ratios. (G) FOV example showing spatial cell classification relative to margin. (H) Comparison of spatial distributions of myeloid cells. (I) Heatmap of myeloid-other cell type co-localization within 50 μm. (J and K) Boxplots showing the co-localization of macrophages with CD8⁺ T cells (J) and cancer cells (K). (L and M) Pearson correlation analyses between the proportions of CXCL9+ Macs (L) or SPP1+ Macs (M) and the counts of other cell types in FOV. (N) Correlation between cancer cell hypoxia signatures and minimal distance to specified cell types. (O) Visualization of the cancer hypoxia signature and the spatial distribution of CD8⁺ T cells and SPP1+ Macs. (P) Correlation between the proportions of CXCL9+ and SPP1+ in each FOV. Boxplots show medians and interquartile ranges; scatterplots show Pearson’s R and regression lines with 95% confidence intervals.

Figure 2

Enriched interaction between hypoxia-driven SPP1+ Macs and tumor niches associated with dismal prognosis in residual TNBC post-NAC (A) CellCharter cluster stability across varying cluster numbers (left) and a representative spatial distribution of niches derived from clusters (right). (B) Heatmap of cell type enrichment across identified niches. (C and D) Boxplots comparing niche proportions (C) and cancer cell hypoxia signatures across tumor-enriched niches (D) between ND and D groups. (E) Differential neighborhood enrichment analysis between tumor-enriched (source) and TME-enriched (target) niches in ND vs. D. Asterisks: p < 0.05; two-sided tests. (F) Bar plot showing information flow differences in inferred networks between ND and D, with significant signaling pathways ranked by variation. (G) Chord diagram of SPP1 signaling interactions between cell types. (H) Bubble plot of significant interactions from SPP1+ Macs to other cells. For all the boxplots, data are presented as median with interquartile range.

Figure 3

Spatial and molecular interactions between CXCL9+ Macs and CD8⁺ T cells (A and B) Scatterplots of correlations between average IFNG (A) or IFNA1 (B) expression levels and CXCL9+ Mac abundance per FOV. (C) Dotplot showing IFNG, IFNA1, and receptor expression across cell types. (D) Correlation analysis of CD8⁺ T cell abundance with average CCL/CXCL chemokine gene expression per FOV. Spearman’s rank correlation estimates (circles) are shown with 95% confidence intervals. (E) Boxplot of spatial autocorrelation using bivariate Moran’s I between CD8A and chemokine genes. Data are presented as median with interquartile range. (F) Dotplot showing CCL/CXCL chemokine gene expression across cell types. (G) Chord plot of ligand-receptor pairs in CCL pathways mediating interactions received by CD8⁺ T cells. (H) CellChat communication probability analysis of upregulated ligand-receptor interactions in CCL pathways received by CD8⁺ T cells in ND vs. D FOVs (p < 0.01).

Figure 4

Validation of the prognostic significance of CXCL9+ Macs and CD8⁺ T cells using bulk RNA-seq (A) Violin plots comparing CXCL9+ Macs, SPP1+ Macs, and CD8⁺ T cell proportions by prognosis (Wilcoxon rank-sum test). Data are presented as median and interquartile range. (B–D) Scatterplots correlating SPP1+ and CXCL9+ Macs with CD8⁺ T cell proportions (B), CXCL9+ Mac proportion with survival duration (C), and SPP1+ Mac proportion with survival duration (D). (E) Kaplan-Meier survival curves grouped by SPP1+ Mac and CXCL9+ Mac proportions (groups determined by median values). (F) Forest plots of Cox multivariable regression for overall survival (left) and disease-free survival (right). Hazard ratios shown as rectangles with 95% confidence intervals. ∗p < 0.05; ∗∗p < 0.01; ∗∗∗p < 0.001. For scatterplots, R represents Pearson’s correlation coefficient; solid lines show fitted linear regression values with 95% confidence intervals.

Figure 5

Spatial profiling reveals distinct cell type distributions and immune-tumor interactions in residual non-TNBC post-NAC based on patient outcomes (A) Flightpath plot of InSituType cell identification showing posterior probability. (B) Dotplot of signature marker expression by cell type. (C) UMAP of reclustered myeloid cells. (D) Dotplot of signature markers for myeloid cells. (E) Radar plot of normalized M1/M2 polarization, angiogenesis, and phagocytosis scores across macrophage subtypes. (F and G) Bar plots of log₂ fold changes in cell abundance between ND and D groups: all cell types (F) and TME (G), adjusted for intrinsic subtype. ∗p < 0.05. (H) Heatmap of average normalized co-localization scores between cancer subtypes (rows) and TME cells (columns) using 50 μm radius. (I) Heatmap of Pearson correlations between cancer gene signatures (rows) and TME cells (columns) based on minimum distance to cancer cells. Only adjusted p < 0.05 shown; |R| > 0.1 annotated. (J) Dotplots of cell type enrichment by niche. (K) Differential neighborhood enrichment analysis between tumor-enriched (source) and TME-enriched (target) clusters in ND vs. D non-TNBC. (L and M) Representative FOVs showing a C0-dominant region in D non-TNBC (L) and a region with minimal C0 in ND cases (M). (N) Heatmaps comparing co-localization scores between cancer subtypes and TME cells in ND vs. D groups (t tests). ∗p < 0.05; ∗∗p < 0.01.

Figure 6

Basal-like subtype and homologous recombination deficiency as determinants of dismal prognosis in residual HR+HER2− BC post-NAC (A) Intrinsic subtype prediction in HR+HER2− bulk RNA-seq dataset. (B) Somatic mutation profiles of HR+HER2− BCs by prognosis (Fisher’s test). COSMIC Cancer Gene Census v.90 genes are highlighted in red. (C and D) Violin plots of SBS1 and SBS3 proportions (C) and HRD scores (D) by prognosis. (E) HRD score and SBS3 proportion comparison between basal and non-basal intrinsic subtypes. (F) Correlation between SBS3 proportion and HRD score in basal and non-basal subtypes. R: Pearson’s correlation coefficient; solid line: fitted linear regression with 95% confidence intervals. (G and H) Kaplan-Meier survival curves stratified by SBS3 proportion (G) and HRD score (H) using median cutoff. (I) Kaplan-Meier survival curves for basal vs. non-basal intrinsic subtypes. (J) Forest plot of Cox multivariable regression for disease-free survival. Hazard ratios are shown as rectangles with 95% confidence intervals. ∗p < 0.05; ∗∗p < 0.01; ∗∗∗p < 0.001. For all the violin plots, central lines represent the median; boxes indicate the interquartile range.

Figure 7

Comparative analysis of genomic and amplicon characteristics in HER2+ BC based on patient outcomes (A) Somatic mutation profiles of HER2+ BCs by prognosis. COSMIC Cancer Gene Census v.90 genes are highlighted in red. (B) Lollipop plot showing ERBB2 somatic mutation loci in D HER2+ group. (C) Violin plots of HRD score and mutational signature proportion differences between D and ND HER2+ groups. (D) Distribution of oncogene-carrying focal amplifications with characteristics (interchromosomal, intrachromosomal, and copy counts) in HER2+ patients (N = 18) grouped by prognosis and ecDNA presence. (E) Distribution of ecDNA-bearing patient classification in HER2+ subtype between D and ND conditions. (F–J) Comparison of amplicon features between D and ND conditions: (F) chromosome counts, (G) amplicon size, (H) oncogene counts, (I) amplicon copy number, and (J) interval counts. ecDNA + ChrAmp, all amplicons; ecDNA Only, ecDNA amplicons; ChrAmp Only, ChrAmp amplicons. p value is from two-sided Wilcoxon rank-sum test. For all the violin plots, central lines represent the median; boxes indicate the interquartile range.