Chemotherapy induces dynamic immune responses in breast cancers that impact treatment outcome

Abstract

To elucidate the effects of neoadjuvant chemotherapy (NAC), we conduct whole transcriptome profiling coupled with histopathology analyses of a longitudinal breast cancer cohort of 146 patients including 110 pairs of serial tumor biopsies collected before treatment, after the first cycle of treatment and at the time of surgery. Here, we show that cytotoxic chemotherapies induce dynamic changes in the tumor immune microenvironment that vary by subtype and pathologic response. Just one cycle of treatment induces an immune stimulatory microenvironment harboring more tumor infiltrating lymphocytes (TILs) and up-regulation of inflammatory signatures predictive of response to anti-PD1 therapies while residual tumors are immune suppressed at end-of-treatment compared to the baseline. Increases in TILs and CD8+ T cell proportions in response to NAC are independently associated with pathologic complete response. Further, on-treatment immune response is more predictive of treatment outcome than immune features in paired baseline samples although these are strongly correlated.

Fig. 1. Study design and consort diagram.

a Study design diagram. AC Adriamycin (Doxorubicin) + Cyclophosphamide, T Taxol (Docetaxel), H Herceptin (Trastuzumab), NAC neoadjuvant chemotherapy, pCR pathologic complete response. (b) Consort diagram for patient enrollment, sample collection and data generation. Tissues were not collected for nine patients who withdrew consent. *55 patients were excluded from NGS data generation due to tissue collection failure (33), follow-up failure (7), withdrawal of consent (6), non-breast cancer pathology (4), distant progression and exclusion from surgery (3), NAC termination (1) and expiration (1). Two non-pCR patients had paired T1, T3 samples. WTS whole transcriptome sequencing.

Fig. 2. Differential expression patterns across treatment times.

Expression patterns of representative pathways (a) and genes (b) mapped to the three DE gene clusters. GSVA scores in a were centered and scaled to z-scores for each pathway geneset. Statistical significance of DE was determined using linear mixed-effects regression analysis (LMER) and shown as sign*−log10(p-value) based on the sign of the t-statistics. Source data are provided as a Source Data file.

Fig. 3. Comparing differential expression patterns across subtypes.

a Aggregate expression patterns of three DE gene clusters in the overall cohort and subtypes. Gray lines represent individual gene expressions and red lines represent the summarized scores (GSVA) for each cluster. The ER+/HER2+ subtype was excluded due to a lack of T3 samples. “ER+ combined” includes ER+ and ER+/HER2+ subtypes. “HER2+ combined” includes ER+/HER2+ and HER2+ subtypes. Error bars represent the standard deviation of the scaled log2TPM value at each time point. Sample sizes used in the analysis are the following. All: n = 112 (T1), n = 88 (T2), n = 27 (T3); ER+: n = 26 (T1), n = 18 (T2), n = 7 (T3); ER+ Combined: n = 49 (T1), n = 36 (T2), n = 8 (T3); HER2+: n = 19 (T1), n = 17 (T2), n = 5 (T3); HER2+ Combined: n = 42 (T1), n = 35 (T2), n = 6 (T3); TN: n = 44 (T1), n = 35 (T2), n = 14 (T3). Source data are provided as a Source Data file. b Aggregate expression patterns of pathways mapped to the three DE clusters in the overall cohort and subtypes. Gray lines represent individual gene signatures (GSVA scores) and red lines represent the averages for each sample group. Error bars represent the standard deviation of the scaled GSVA score at each time point. Sample sizes used to derive statistics are the same as in Fig. 2a. Source data are provided as a Source Data file.

Fig. 4. Tumor infiltrating lymphocyte density increased during treatment then decreased below baseline at surgery time.

TIL density distributions over three time points were compared between tumors from pCR and RD patients in the overall cohort (a), TN (b), and non-TN subtypes (c). d Distribution of stromal TIL scores over time in TN tumors stratified into pCR and RD groups. Asterisks indicate statistical significance based on linear mixed effects regression (LMER) adjusting for tumor purity and subtype as covariates. *0.01 < p < 0.05; **0.001 < p < 0.01; ***p < 0.001. See Supplementary Table 2 for exact p-values. For all box-and-whisker plots, the box is bounded by the first and third quartile with a horizontal line at the median and whiskers extend to the maximum and minimum value. Source data are provided as a Source Data file. e, f Multiplex immunofluorescence (IF), H&E and PD-L1 IHC images showing the same regions of tumors taken at T1 and T2 from a TNBC patient—BR294 (e) and a HER2+ patient—BR308 (f). Markers include CD45RO, CD3, CD4, CD8, PD-L1, and Pan-CK. Chromogenic IHC and multiplex IF assays were only performed once on tumor biopsy samples following assay optimization. Scale bars in the lower-left corner of the micrographs show 200 μm.

Fig. 5. Neoadjuvant chemotherapy induced dynamic changes in tumor associated immune states.

a Integrated clustering of immune signatures and most abundant immune cell fractions. b Sankey plot of immune state changes from T1 to T2 and from T2 to T3. The text denotes the number of samples in each immune state and the percentage of samples that switch from one immune state to another. Source data are provided as a Source Data file.

Fig. 6. Comparing immune states at on-treatment vs. baseline.

a Integrated map of immune signatures and immune cell fractions for paired T2 and T1 samples. Immune state classifications were based on T2 samples. b, c Scatterplots showing CYT scores in paired T1 vs. T2 samples in the overall cohort (b) and different subtypes (c). Source data are provided as a Source Data file. d, e Scatterplots showing TIL densities in paired T1 vs. T2 samples in the overall cohort (d) and different subtypes (e). Statistical significance of correlation was determined using the Spearman method. Shaded error band around the regression line indicates the 95% confidence interval. Source data are provided as a Source Data file.

Fig. 7. Impact of immune response on treatment outcome.

a Forest plot showing the associations of different immune cell fractions vs. NAC response combined T1 and T2 samples. The x-axis shows the log odds ratio of % cell fractions in pCR vs. RD. Asterisks indicate statistical significance based on multiple regression adjusting for subtype as a covariate: *p = 0.038, **p = 0.008, ***p = 0.0008. b Independent associations between clinical factors (columns) and molecular features (rows) associated with pCR. Value in each cell represents –log10(p-value) of the association between each feature and clinical factor after adjusting for the confounding effect of other factors. Features were ranked in descending order by % variable usage as determined by bootstrapping analysis, the percentage of runs in which the elastic net model selected the variable to predict pCR status. Source data are provided as a Source Data file. c Comparison of variable importance estimated for T1 and T2 immune features based on bootstrapping analyses. Statistical significance was determined using two-sided Wilcoxon rank sum test. The box is bounded by the first and third quartile with a horizontal line at the median and whiskers extend to the maximum and minimum value. Source data are provided as a Source Data file. d Distributions of intraepithelial TIL levels between T1 and T2 samples and in pCR vs. RD cohorts. Asterisks indicate statistical significance calculated by using one-sided Chi-squared test: ***p = 4.06E−5. e H&E images of tumor biopsy sections taken from the same TNBC patient (OB-16-0006) showing an increased abundance of intraepithelial TILs between T1 baseline and T2 on-treatment. Sub-panel ii shows the tumor region annotated in subpanel i at a higher zoom with the intraepithelial lymphocytes highlighted in green. H&E staining was only performed once on tumor biopsy samples.