Changes in Peripheral and Local Tumor Immunity after Neoadjuvant Chemotherapy Reshape Clinical Outcomes in Patients with Breast Cancer

This article is featured in Highlights of This Issue, p. 5541

Combination of conventional chemotherapy with an immunotherapeutic targeting programmed death-ligand 1 (PD-L1), atezolizumab, was recently approved for the treatment of patients with metastatic triple-negative breast cancer (TNBC) based on results from a phase III clinical trial (1). Furthermore, addition of the anti–programmed death-1 (PD-1) mAb pembrolizumab to neoadjuvant chemotherapy (NAC) can significantly improve TNBC pathologic complete response (pCR) rates (2). Thus, existing clinical data indicate that chemotherapy combinations with immunotherapy demonstrate enhanced efficacy compared with chemotherapy alone. However, these results suggest a growing need to better understand how chemotherapy modulates the tumor immune microenvironment (TIME).

High levels of stromal tumor-infiltrating lymphocytes (sTIL) in the pretreatment biopsy are predictive of pCR in patients with TNBC treated with NAC (3). In NAC-treated patients with TNBC with residual disease (RD) at surgery or in untreated primary TNBC tumors, higher sTILs in the resected tumor also confer improved prognosis (4–7). However, the immunomodulatory effect of chemotherapy on sTILs in patients and how chemotherapy influences the TIME are poorly understood. In a study of patients with non–small cell lung carcinoma, tumors treated with NAC had higher expression of PD-L1 and higher density of CD3⁺ T cells, suggesting NAC may be immunomodulatory in some tumor types (8). In addition, in a small breast cancer study, post-NAC natural killer (NK) cells and IL6 expression were associated with better response. However, this study was limited by inclusion of only a small number of TNBCs (n = 4; ref. 9). Intriguingly, a larger study of patients with breast cancer (n = 60) showed that pre-NAC sTILs and higher pre-NAC expression of cytotoxic T-cell markers and cytokines were associated with higher pCR rate. However, this study did not examine long-term outcomes in patients with RD, for whom prognosis is worse than those with pCR and only included a limited number of TNBC samples (n = 13; ref. 10). Furthermore, whether immunomodulatory effects of NAC are locoregional or systemic (i.e., able to be detected in the peripheral blood) is unknown.

To address this gap in knowledge, we examined expression patterns of immune-related genes before and after NAC in a series of 83 breast tumors, including 44 TNBCs, from patients with RD. As patients with pCR generally experience excellent outcomes, we chose to focus on patients with RD, who may benefit from additional therapies. Changes in gene expression patterns in the TIME were tested for association with recurrence-free (RFS) and overall survival (OS). T-cell receptor sequencing was performed on a subset (n = 15) of tumors. In addition, in 4 patients undergoing NAC, PD-1–high (PD-1^HI; top 20% expressors of CD8 or CD4 cells) and PD-1–negative (PD-1^NEG) CD3⁺CD8⁺ peripheral blood mononuclear cells (PBMC) were profiled using single-cell RNA sequencing (scRNAseq) and multiplexed cytokine secretion assays. Finally, we used an scRNAseq-derived signature of activated cytolytic cells to measure immune activation in the peripheral blood of two cohorts of patients: a Vanderbilt cohort consisting of 34 patients after NAC (collected within 2 weeks prior to surgery) and 24 untreated patients, and a cohort from the Dana-Farber Cancer Institute (DFCI; Boston, MA) consisting of 30 hormone receptor–positive, HER2⁻ patients treated with NAC (including bevacizumab) as part of a clinical trial (11). We then tested the association of this signature with surgical outcome (pCR or RD burden) and postsurgical cancer recurrence.

Intriguingly, higher expression of the gene signature at the time of surgery was associated with higher disease burden [i.e., in those with RD who experienced disease recurrence within 3 years following surgery or those with the highest residual cancer burden (RCB)]. Thus, peripheral cytotoxic activity, guarded by immune checkpoints, may reflect ongoing micrometastatic and primary disease burden, and could be a useful biomarker for disease recurrence and possibly immune checkpoint inhibitor benefit.

Three cohorts of patients were combined for the tumor-profiling study. All included patients received neoadjuvant therapy and had RD and matched pretreatment tissue was required for inclusion. All but 2 estrogen receptor–positive (ER⁺) patients received cytotoxic chemotherapy as part of their regimen. Four patients (ER⁺) received courses of hormone therapy as part of their neoadjuvant regimen, two of which were in conjunction with cytotoxic chemotherapy. Five patients (HER2⁺) received HER2-directed therapy as part of their neoadjuvant regimen.

For the “Peru” cohort, clinical characteristics and molecular analysis of the patients (n = 48 with matched pretreatment tissue) were previously described at the Instituto Nacional de Enfermedades Neoplásicas (12). Clinical and pathologic data were retrieved from medical records under an institutionally approved protocol (INEN IRB 10–018). For the “VICC” cohort, which included PBMC and whole blood analyses, clinical and pathologic data were retrieved from medical records under an institutionally approved protocol (VICC IRB 030747). For the DARTMOUTH cohort, patient samples were collected under a protocol approved by the Dartmouth College Institutional Review Board and the waiver of the subject consent process was IRB approved (IRB 28888). Metadata for the primary cohort of patients is presented in Supplementary Table S1. For the peripheral blood study, two cohorts of patients were used (summarized in Table 2). For the Vanderbilt cohort, all blood was collected within 14 days preceding definitive surgery. Metadata for the Vanderbilt cohort is presented in Supplementary Table S7. For the DFCI cohort (11), all blood was collected in the interlude between completion of NAC and definitive surgery.

sTILs were analyzed using full face hematoxylin and eosin sections from pre-NAC diagnostic biopsies or post-NAC RD surgical specimens. Samples were scored according to the International TILs Working Group Guidelines (13–15). The predefined cut point of 30% (4) was used for all survival analyses.

Gene expression and gene set analysis on pre- and post-NAC formalin-fixed tissues were performed using the nanoString Pan-Cancer Immunology panel (770 genes) according to the manufacturer's standard protocol. Data were normalized according to positive and negative spike-in controls, then endogenous housekeeper controls, and transcript counts were log transformed for downstream analyses. Normalized linear data are presented in Supplementary Table S2. Gene sets were calculated by summing the log₂-transformed normalized NanoString counts for all genes contained in a given gene set (Supplementary Table S3). Samples were simultaneously assayed for PAM50 molecular subtyping using a custom-designed 60-gene Elements panel (Supplementary Table S4). Briefly, 10 μm sections of diagnostic biopsies or residual tumors were used for RNA preparation [Promega Maxwell 16 RNA formalin-fixed, paraffin-embedded (FFPE)] and 50 ng of total RNA >300 nt (assayed on an Agilent Tapestation 2200 Bioanalyzer) was used for input into nCounter hybridizations for Pan-Cancer Immunology panels or 500–1,000 ng RNA for PAM50 analysis. Data were normalized according to positive and negative spike-in controls, then endogenous housekeeper controls, and transcript counts were log transformed for downstream analyses. Subtype prediction was performed in R using the genefu package.

For the 8-gene signature analysis in whole blood, a custom NanoString Elements was constructed to measure the gene expression levels of PDCD1, NKG7, LAG3, GZMH, GZMB, GNLY, FGFBP2, HLA-DRB5, and HLA-G, as well as 3 normalization control genes (PTPRC, RPL13a, TBP). Probe design is shown in Supplementary Table S9. RNA was isolated from whole blood (Promega Maxwell 16 Simply RNA Blood) and 100–200 ng was used for input into the nCounter analysis. Data were normalized as above. Linear normalized data are presented in Supplementary Table S8 for the Vanderbilt Cohort.

On day 1, cryopreserved PBMCs were thawed and resuspended in complete RPMI media with IL2 (10 ng/mL) at a density of 1–5 × 10⁶ cells/mL. Cells were recovered at 37°C, 5% CO₂, overnight. Plates were prepared by coating with anti-human CD3 (10 μg/mL in PBS, 200–300 μL/well) in a 96-well flat-bottom plate at 4°C, O/N. On day 2, nonadherent cells for each sample were collected and viability was confirmed, with dead cell depletion by Ficoll. For each sample, where sufficient, volume was split in half for each of the following negative isolations: with one half of cells from each sample, CD4 T cells were isolated with CD4⁺ negative isolation kit following Miltenyi protocol (130-096-533); with the other half of cells from each sample, CD8 T cells were isolated with CD8⁺ negative isolation kit following Miltenyi protocol (130-096-495). The PD-1⁺ and PD-1⁻ subsets were from isolated CD4 or CD8 T cells by staining with PE-conjugated anti–PD-1 antibody using the manufacturer's protocol (Miltenyi, 130-096-164) as follows: (i) stain each subset with 10 μL stain: 100 μL Robosep buffer for every 1 ×10⁷ total cells; (ii) incubate at 4°C for 10 minutes; (iii) rinse cells by adding 1–2 mL of Robosep and C/F at 300 × g for 10 minutes; (iv) aspirate supernatant and resuspend cell pellets in 80 μL buffer per 1 × 10⁷ total cells. PD-1⁺ cells were then isolated with anti-PE microbeads following the manufacturer's protocol (Miltenyi, 130-097-054). Cells were resuspended in complete RPMI media at a density of 1 × 10⁶/mL and seeded into wells of the CD3-coated 96-well flat-bottom plate with soluble anti-human CD28 (5 μg/mL). Plates were incubated at 37°C, 5% CO₂ for 24 hours. On day 3, supernatants (100 μL/well) were collected from all wells and stored at −80°C for population assays. T cells were collected and stained with Brilliant Violet cell membrane stain and AlexaFluor-647–conjugated anti-CD8 at room temperature for 20 minutes, rinsed with PBS and resuspended in complete RPMI media at a density of 1 × 10⁶/mL.

Approximately 30 μL of cell suspension was loaded into the IsoCode Chip and incubated at 37°C, 5% CO₂ for an additional 16 hours. Protein secretions from approximately 1,000 single cells were captured by the 32-plex antibody barcoded chip and analyzed by fluorescence ELISA-based assay (16–21). Polyfunctional T cells that cosecreted 2+ cytokines per cell were evaluated by the IsoSpeak software across the five functional groups: Effector: Granzyme B, TNFα, IFNγ, MIP1α, Perforin, TNFβ; Stimulatory: GM-CSF, IL2, IL5, IL7, IL8, IL9, IL12, IL15, IL21; Chemoattractive: CCL11, IP-10, MIP-1β, RANTES; Regulatory: IL4, IL10, IL13, IL22, sCD137, sCD40L, TGFβ1; and Inflammatory: IL6, IL17A, IL17F, MCP-1, MCP-4, IL1β.

TP53 gene sequencing was performed using either the Foundation Medicine assay as reported previously (12) or using the SureMASTR TP53 sequencing assay (Agilent). For the latter, purified DNA from FFPE breast tumor sections were amplified and sequenced according to the manufacturer's standard protocol. Samples were sequenced to a depth of approximately 10,000 and mutations were called using the SureCall software (Agilent). Mutation allele frequency was set at 5% and only likely functional mutations (early stops, frameshift deletions, and known recurrent hotspot single-nucleotide variation mutations) were selected for sample annotation.

TCR sequencing and clonality quantification was assessed in FFPE samples of breast cancer specimens or PBMCs. For FFPE tissue, DNA or RNA was extracted from 10 μm sections using the Promega Maxwell 16 FFPE DNA or FFPE RNA kits and the manufacturer's protocol. For PBMCs, PD-1^HI and PD-1^NEG CD8⁺ T cells were sorted by FACS from samples isolated from EDTA collection tubes and processed using a Ficoll gradient. At least 100 K cells were collected, centrifuged, and utilized for RNA purification. TCRs were sequenced using survey level immunoSEQ (DNA; Adaptive Biotechnologies) and the Immunoverse assay (RNA; ArcherDX), as described previously (22, 23). Sequencing results were evaluated using the immunoSEQ analyzer version 3.0 or Archer Immunoverse analyzer. CDR3 sequences and frequency tables were extracted from the manufacturers' analysis platforms and imported into R for analysis using the Immunarch package (https://immunarch.com; ref. 24) in R. Shannon entropy, a measure of sample diversity, was calculated on the clonal abundance of all productive TCR sequences in the data set. Shannon entropy was normalized by dividing Shannon entropy by the logarithm of the number of unique productive TCR sequences. This normalized entropy value was then inverted (1 − normalized entropy) to produce the “clonality” metric. TCRβ clonotypes and metadata based on the primary cohort (Supplementary Fig. S5; Adaptive) are included as Supplementary Dataset 1. TCRβ clonotypes and metadata based on prospectively collected peripheral blood and tumor from Fig. 4 (Archer) are included as Supplementary Dataset 2.

PBMCs were isolated from EDTA collection tubes, processed using a Ficoll gradient, and cryopreserved in 10% DMSO 90% FBS. Upon thaw, whole live PBMCs were prepared by depleting dead cells using a dead cell removal Kit (Miltenyi, catalog no.: 130–090–101). PD-1^HI and PD-1^NEG CD8⁺ T cells were sorted by FACS from PBMCs. Each sample (targeting 5,000–10,000 cells/sample) was processed for single cell 5′ RNAseq utilizing the 10x Chromium system. Libraries were prepared using P/N 1000006, 1000080, and 1000020 following the manufacturer's protocol. The libraries were sequenced using the NovaSeq 6000 with 150 bp paired end reads. RTA (version 2.4.11; Illumina) was used for base calling and analysis was completed using 10× Genomics Cell Ranger software v2.1.1. Data were analyzed in R using the filtered h5 gene matrices in the Seurat (25, 26) package. Briefly, samples were merged, and all cells were scored for mitochondrial gene expression (a marker of dying cells) and cell-cycle genes to determine phase. Data were transformed using SCTransform, regressing against mitochondrial gene expression and cell-cycle phase. Dimensional reduction was performed using Harmony (22). Missing values were imputed using the RunALRA function in the SeuratWrappers package (23). For the whole PBMC single-cell data, cell types were assigned to individual cells using SingleR (27). BlueprintEncodeData was used as a reference (28, 29).

All statistical analyses were performed in R or Graphpad. For survival curves (RFS and OS), the log-rank statistic was reported (or trend-test for >2 groupings). For gene-level analysis, nominal P values were calculated using a Cox proportional hazards model, and then an adjusted P value (q value) was calculated on the basis of the FDR (30). All single-cell statistical analyses were calculated in R using the Seurat package (25). Shared Nearest Neighbors were calculated using the Harmony reduction, and clusters were identified at a resolution of 0.3, defining 5 total clusters. UMAP was performed for visualization, and missing values were imputed using ALRA (23). Visualization and graph generation was performed in R. Some heatmaps were made using the complex heatmap package in R (31). P-value cut-offs displayed on plots correspond to “ns” equals P > 0.05, * equals 0.01 < P ≤ 0.05, ** equals 0.001 < P ≤ 0.01, *** equals 0.0001 < P ≤ 0.001, **** equals P ≤ 0.0001. Code used to generate figures can be accessed at https://github.com/MLAxelrod/Immunologic_changes_with_chemotherapy_inTNBC.

We procured matched archived pretreatment (diagnostic biopsy) and posttreatment (RD surgical specimen) tumor specimens from a series of 83 patients, including 44 patients with TNBC. Importantly, we refined our study to include only patients who had RD at surgery for analysis, thereby excluding patients who achieved pCR. This was a purposeful selection strategy, as patients with pCR usually experience good outcomes, and we instead chose to focus on patients with RD for whom additional risk stratification could identify those who are most likely to benefit from additional therapies. Metadata for the patients, including treating institution, molecular subtype (PAM50), RFS, OS, TP53 mutation status, and other molecular and clinical data are summarized in Table 1. Individual patient-level data are available in Supplementary Table S1.

As sTILs have been described and rigorously validated in breast cancer as both a prognostic factor (in surgical specimens for postsurgical outcomes, particularly in TNBC and HER2⁺ cancers), and a predictive factor (in diagnostic biopsies for benefit from NAC), we first asked whether these findings were consistent with our study cohort. Using the published cutoff (30%) of sTILs (4), we found that higher abundance of sTILs in the post-NAC RD in patients with TNBC (n = 44) was significantly prognostic for both RFS (log-rank P = 0.019) and OS (P = 0.05; Fig. 1A). Representative histology of high (>30%) and low (<30%) sTILs are shown in Supplementary Fig. S1. Interestingly, pre-NAC sTILs in the diagnostic biopsy were not prognostic for outcomes in patients with TNBC (Supplementary Fig. S2A), presumably due to the selection strategy of including only patients who lacked pCR. Consistent with prior literature that the prognostic and predictive effect of sTILs differs by breast cancer subtype (32), neither pre-NAC nor post-NAC sTILs were prognostic for OS when considering our entire cohort (n = 83). However, post-NAC sTILs were prognostic for RFS (P = 0.031) in the whole cohort (Supplementary Figs. S2B and S3A). This effect seems primarily driven by TNBC tumors as post-NAC sTILs are not prognostic of either RFS or OS in non-TNBCs (Supplementary Fig. S3B). Neither pre-NAC nor post-NAC sTILs were prognostic in either ER⁺ or HER2⁺ patients only (Supplementary Fig. S4). This may be limited by our sample size as sTILs have been previously shown to be associated with longer RFS in HER2⁺ patients with cancer and shorter OS in luminal/HER2⁻ patients with cancer (32). Stratifying patients with TNBC by whether sTILs were qualitatively increased or decreased/equivocal in the surgical resection compared with the diagnostic biopsy did not provide any prognostic capability in this cohort (Supplementary Fig. S5A). Interestingly, most patients, regardless of clinical subtype, had a decrease in sTILs over the course of NAC (Supplementary Fig. S5B). Importantly, our cohort includes only patients with RD, and there is currently not a validated method for quantifying sTILs in pCR where, by definition, a tumor is no longer present. Thus, in our cohort, abundance of sTILs has the strongest prognostic effect for the post-NAC surgical resection specimen in TNBC tumors with an incomplete response to NAC. These findings, consistent with both retrospective studies and analyses from randomized controlled trials, prompted us to perform more detailed molecular studies aimed at understanding how NAC influences the TIME.

To measure transcriptional changes occurring in the TIME induced by NAC, we performed gene expression profiling for a series of 770 immune-related genes using nanoString (Pan-Cancer Immunology panel), before and after NAC in the entire cohort (n = 83). Transcriptional patterns and hierarchical clustering for all data primarily segregated tumors based on receptor status (ER/PR/HER2) and/or molecular subtype, with most luminal/hormone receptor–positive tumors appearing in the first cluster, most HER2⁺ tumors in the second cluster, and most basal-like/TNBC tumors in the third cluster (Fig. 1B). Examining the data as the change in gene expression for each gene after NAC in a patient-matched fashion (Δ expression; post-NAC minus pre-NAC) yielded similar patterns, with a trend of most patients with TNBC having generalized decreased immune gene expression patterns after NAC (Fig. 1C). To test for effects driven by differences in sampling (i.e., pre-NAC samples are biopsies whereas post-NAC samples are surgical resection samples) or treating institution, we performed principal component analyses. We did not detect significant clustering by time of sampling or cohort (Supplementary Fig. S6).

While the TIME change in sTIL abundance did not prognosticate outcome in patients with TNBC, we asked whether changes in individual immune-related genes are associated with outcome. We performed iterative Cox proportional hazards models, using the delta (Δ) of each gene (post-NAC minus pre-NAC) in an independent univariate analysis, for both RFS and OS. All analyses are reported using a nominal P value as well as an FDR (Benjamini–Hochberg method) q value for association with RFS or OS. After correction for FDR (q < 0.10), upregulation of 11 genes were associated with improved RFS, while upregulation of only one gene was significantly associated with worse RFS (CDH1, which encodes e-cadherin) in our TNBC cohort. Interestingly, e-cadherin is known to interact with killer cell lectin-like receptor G1 (KLRG1), an inhibitory receptor expressed by memory T cells and NK cells (33). In contrast, upregulation of a larger number of genes was associated with improved OS (n = 189) or reduced OS (n = 15) at FDR q < 0.10 (Fig. 2A). These genes are listed in Supplementary Table S5. Kaplan–Meier visualization examples of strongly prognostic genes (negative prognostic: CDH1; positive prognostic: CD70) reinforced the prominent association of TNBC disease outcomes with changes in immune gene expression during NAC (Fig. 2B). Conversely, no changes in immune-related gene expression were significantly associated with RFS or OS in patients with non-TNBC at q < 0.10 (Supplementary Fig. S7A).

Dimensional reduction through collapsing individual genes into pathways or defined functions can improve interpretation of high-dimensional data. Thus, we collapsed the gene expression data into bioinformatically categorized immune signatures (sum scores, defined as the summation of the log₂ expression values for all genes in a category). Organization of the data in this manner and testing the signatures (n = 70) for association with RFS and OS yielded a surprising finding: all significant (q < 0.10) gene sets (n = 40 for RFS and n = 60 for OS; listed in Supplementary Table S6) identified in this analysis were associated with good outcome (Fig. 3A). Many of the top-scoring gene sets were associated with T cells, including “T-cell polarization,” “T-cell immunity,” and “T-cell activation,” among others. Although manual inspection revealed some overlap in these gene sets, they were largely composed of signature-exclusive genes (Supplementary Table S3). Kaplan–Meier visualization examples of strongly prognostic gene sets (“NK-cell functions” and “T-cell activation”) reinforced the considerable association of changes in immune gene sets during NAC with outcomes (Fig. 3B). Interestingly, no gene or gene set was significantly associated with RFS or OS in patients with non-TNBC at q < 0.10 (Supplementary Fig. S7A and S7B). When the non-TNBC group was separated into ER⁺ and HER2⁺ groups, no gene or gene set was significantly associated with outcome (Supplementary Fig. S8). Thus, these data suggest that NAC, exclusively in TNBC, could promote immunologic activity leading to improved outcomes in a subset of patients. However, these effects may be related to factors beyond TNBC biology, as hormone receptor–positive and HER2⁺ patients receive additional endocrine or HER2-directed therapy in the adjuvant setting, complicating associations with RFS and possibly OS. Nonetheless, immune-related signatures, particularly those derived from T cells, appeared to be strongly associated with improved outcomes in TNBC.

A robust T-cell response is characterized by oligoclonal expansion of antigen-specific T cells. Therefore, we next asked whether clonality of T cells in the TIME was altered during NAC. In a subset of samples (n = 15; 8 TNBC, 7 non-TNBC), we performed TCRβ chain sequencing using the ImmunoSeq assay to estimate the number of unique T-cell clones (diversity), and the presence of expanded T-cell clones in the TIME before and after NAC. Given the breadth of sTIL fractions observed among breast tumors as well as caveats associated with comparison of samples derived from diagnostic core needle biopsies versus surgical resections, we first verified that the number of productive T cells was associated with estimation of sTILs determined on adjacent sections. We detected a strong association between these parameters (R² = 0.6; P < 0.0001; Supplementary Fig. S9A), raising confidence in the assay results. In this sample set, NAC did not universally alter productive clonality (Supplementary Fig. S9B), a measurement of the number of times the same (productive) TCRβ sequence is represented in the sample, which is a descriptor of T-cell clonal expansion. When stratified by breast cancer subtype, there was no significant change in productive clonality with NAC (one-sample t test). However, TNBC tumors demonstrated a qualitative trend toward decreased clonality after NAC, while non-TNBC tumors trended toward increased clonality after NAC. The difference between these two subgroups approached significance (P = 0.054; two-sample t -test; Supplementary Fig. S9C). There was no association of change in clonality with change in sTILs, suggesting that changes in sTIL abundance after NAC are not necessarily due to expansion of existing clones (Supplementary Fig. S9D).

To further explore changes in T-cell clonality and function in response to chemotherapy, we prospectively collected PBMCs from 4 patients with breast cancer (including 2 TNBCs) before and after NAC (Fig. 4A). In addition, the post-NAC RD (or in one case, pCR residual scar) was analyzed in tandem. On the basis of previous findings demonstrating that tumor-reactive T cells are enriched in the CD8⁺ PD-1^HI population of peripheral T cells (34), we purified CD4⁺ and CD8⁺ cells from each sample by FACS, further stratifying by PD-1^NEG and PD-1^HI status (gating scheme shown in Supplementary Fig. S10). Using a functional fluorescence ELISA-based assay of cytokine (32-plex antibody barcoded chip; Supplementary Fig. S11A) secretion following CD3/CD28 stimulation, we determined that PD-1^HI peripheral T cells had functional capacity, secreting multiple cytokines following activation, and these effects were particularly pronounced in CD8⁺ T cells (Supplementary Fig. S11B). In 2/2 patients with TNBC, the percentage of “polyfunctional” PD-1^HI CD8⁺ T cells—those capable of expressing multiple cytokines after TCR stimulation—were increased following NAC (Fig. 4B). In contrast, 2/2 ER⁺ patients with breast cancer experienced a drop or stasis in the functionality of the PD-1^HI CD8⁺ population of cells following NAC (Fig. 4B). Of note, the patient with ER⁺HER2⁺ disease has a near complete loss of T-cell functionality after NAC. Cytokines produced by individual PD-1^HI CD8⁺ cells in patients with TNBC were primarily effector (e.g., Granzyme B, IFNγ, MIP-1α, TNFα, and TNFβ) and chemoattractive (MIP-1β) cytokines (Fig. 4C). PD-1^HI CD4⁺ T cells also produced primarily effector cytokines including IFNγ and TNFα (Supplementary Fig. S11C).

PD-1^HI CD8⁺ and PD-1^NEG CD8⁺ T cells from pre-NAC and post-NAC blood (except patient 4, for whom a sufficient pre-NAC sample was not available) were also analyzed by TCR sequencing. While the number of detected T cells was consistent among all samples, the clonotypes detected (unique TCRs) were considerably lower in PD-1^HI CD8⁺ T cells (Fig. 4D). This suggests that there are more repetitive sequences detected in the PD-1^HI population, indicating clonal expansion. Consistent with this observation, the proportion of the overall TCR repertoire occupied by expanded clonotypes (large or hyperexpanded clonotypes consisting of greater than 0.1% or 1% of the total repertoire, respectively) was substantially higher in the PD-1^HI than in PD-1^NEG CD8⁺ T-cell fractions (Fig. 4E). We additionally sequenced the TCR repertoire in the post-NAC RD, although the number of T cells sequenced in these samples were limited because of fixation of tissue and small T-cell abundance as a function of total RNA in the bulk samples, and thus should be interpreted with caution. Nonetheless, we found that the similarity (Jaccard index, normalized to size of repertoire detected) of tumor-infiltrating TCRs in the post-NAC sample was universally more similar to the PD-1^HI CD8⁺ peripheral TCR repertoires, compared with the PD-1^NEG CD8⁺ repertories (Supplementary Fig. S12). This suggests that the PD-1^HI peripheral compartment is enriched for similarity to TILs relative to the PD-1^NEG peripheral compartment.

Next, we utilized scRNAseq to describe the post-NAC peripheral PD-1^HI CD8⁺ T-cell populations at the time of surgery in the blood of 2 patients with TNBC: one with RD (Pt. 1) and one with matrix-producing metaplastic TNBC who experienced pCR (Pt. 4). Uniform Manifold Approximation and Projection (UMAP) analysis was performed on Harmony-normalized samples to adjust for intersample technical variation, and we stratified cells based on 5 clusters identified through the Louvain algorithm (Fig. 5A and B). A heatmap of the top 10 most differentially expressed genes by cluster is presented in Fig. 5C. Although the composition of the cells was largely similar, we identified one cluster (cluster 0) which was enriched in Pt. 4. Examination of genes differentially expressed in this cluster of cells suggested a cellular identity concordant with that of highly cytotoxic memory (TBX21-expressing) T cells, which had an abundance of MHC-I (HLA-A/B/C) and MHC-II (e.g., HLA-DRA, HLA-DRB5) family member expression as well as expression of cytolytic and immune checkpoint genes [e.g., LAG3, FCRL6 (35), and higher transcriptional expression of PDCD1; Fig. 5C and D]. Verification of the pattern of expression of key cytolytic and killer-identity genes [GNLY (granulysin), GZMB (granzyme B), and FGFBP2 (killer-secreted protein 37)] showed that these genes were almost exclusively expressed in cluster 0 (Supplementary Fig. S13A). This analysis also demonstrated purity-of-sort in that all clusters expressed CD8A and PDCD1, but not CD4 (Supplementary Fig. S13B).

These data led us to propose two competing hypotheses: (i) Cluster 0 genes, reflective of cytolytic CD8⁺ T cells, are a positive prognostic factor reflective of robust antitumor immunity as evidenced by their enrichment in the patient with metaplastic TNBC with pCR, or (ii) Cluster 0 genes are reflective of ongoing disease including the micrometastatic component that cannot be sampled from the primary tumor. The second hypothesis is supported by the observations that pCR is less prognostic of RFS and OS in metaplastic disease (36, 37) and that rates of recurrence following chemotherapy are higher for metaplastic disease than nonmetaplastic TNBC (36, 38, 39). Interestingly, matrix-producing metaplastic breast cancer (Pt. 4) has been shown to be associated with pCR to NAC, but often can still recur despite pCR (36, 40). Follow-up for this individual patient was immature at the time of reporting, and thus recurrence, and therefore presence of micrometastatic disease at the time of sampling, cannot be ruled out.

To gain a deeper understanding of cluster 0 genes, we performed additional scRNAseq on post-NAC whole PBMCs from 2 patients with TNBC with pCR (Pt. 4, described above, and Pt. 5, not used in any prior analysis). UMAP was used for dimensionality reduction on Harmony-normalized samples. SingleR was used to computationally assign cell-type annotations (Fig. 5E). Analysis of the top 5 differentially expressed transcripts by cluster (Supplementary Fig. S14A) and expression of key cell-type identity markers (Supplementary Fig. S14B) validated annotations by SingleR. Low-confidence cell-type annotations were collapsed into the category “other” and make up a minor fraction of cells (Supplementary Fig. S15A). Seven genes enriched in cluster 0 (PDCD1, NKG7, LAG3, GZMH, GZMB, GNLY, FGFBP2) and one gene strongly deenriched in Pt. 4 (HLA-G) were chosen as an 8-gene signature for downstream applications, including validation in a larger cohort of patients (Supplementary Fig. S15B and S15C). Expression of this 8-gene score (PDCD1 + NKG7 + LAG3 + GZMH + GZMB + GNLY + FGFBP2 – HLA-G) was the highest in CD8⁺ T cells (as expected, given the derivation from PD-1^HI CD8⁺ T cells) and NK cells (Fig. 5F). This finding was similar in both patients (Supplementary Fig. S16A and S16B). Strong expression of this signature in post-NAC whole PBMCs led us to test whether this signature was predictive of response in archived whole blood samples from patients with breast cancer.

To determine whether cytolytic signatures representative of cluster 0 genes were associated with disease outcome, we evaluated archived whole blood from a series of 58 patients with breast cancer. All samples were collected within 14 days preceding surgical resection for primary breast cancer, with 34 samples having received NAC, in addition to 24 samples from untreated patients. As described above, a series of eight genes enriched in cluster 0 (PDCD1, NKG7, LAG3, GZMH, GZMB, GNLY, FGFBP2, HLA-DRB5), one gene enriched in Pt. 1 (RCB-II) over Pt. 4 (pCR; HLA-G; Supplementary Fig. S15B and S15C), and three normalization control genes (PTPRC, RPL13a, TBP; ref. 41) were selected for a 12-gene custom NanoString gene expression analysis. HLA-G has been described as an immune checkpoint which can dampen anti-tumor immune responses (42–44), and thus we expected HLA-G expression to be inversely correlated with the other selected genes, as is the case in our scRNAseq dataset. One of these genes performed poorly (HLA-DRB5), likely due to frequent polymorphisms in the gene leading to highly variable probe binding and was therefore omitted from further analysis. Information on the presence of pCR/RD at surgery, ER/PR/HER2 status, and clinical follow-up (recurrence at 1,000 days after surgery for patients with RD) was collected (Table 2; Supplementary Table S7).

Nearly, all tested genes demonstrated a pattern supporting the hypothesis that gene expression in whole blood is associated with ongoing disease, being highest (or lowest in the case of HLA-G) in untreated patients (who have ongoing tumor burden by virtue of not having received therapy prior to surgery) and those with RD compared with those with pCR. However, it is important to note that there is a high degree of heterogeneity in the untreated patients, possibly reflecting heterogeneity in tumor size, strength of antitumor immune response at baseline, or antitumor response following systemic therapy. Furthermore, among patients with RD, higher expression (or lower in the case of HLA-G) tended to be observed in patients who had early recurrences in the first 3 years following surgery (and thus may have had micrometastatic disease at the time of surgery). Several of these genes (FGFBP2, GNLY, PDCD1, LAG3, and NKG7) were also significantly differentially expressed or approached statistical significance across the outcome groups (Kruskal–Wallis test). Comparisons were particularly striking between the group of patients with RD who experienced early disease recurrence and the group with pCR following NAC (post hoc Dunn test; Fig. 6A). A composite score of PDCD1 + NKG7 + LAG3 + GZMH + GZMB + GNLY + FGFBP2 – HLA-G also demonstrated statistically significant associations with presence of ongoing disease (Fig. 6B). Interestingly, expression levels of these genes did not always correlate with one another, indicating heterogeneity in their expression patterns and some degree of independence (Fig. 6C). Trends in gene expression were similar for patients with TNBC and non-TNBC, but in-depth subgroup analyses were limited by sample size. Thus, peripheral antitumor immunity in blood may be a useful measure of persistent residual primary or micrometastatic disease and could identify patients likely to benefit from additional therapy.

To extend these findings, we evaluated gene expression in a second cohort of archival blood from patients with breast cancer treated at DFCI (Boston, MA; n = 30), with blood collected in the interlude between completing NAC and definitive surgery. Notably, this cohort had differing baseline characteristics (summarized in Table 2). All of the patients were hormone receptor positive (a subtype known to have a lower pCR rate than patients with TNBC), HER2⁻, and all were treated with bevacizumab. Intriguingly, expression of the 8-gene composite score was also correlated with enhanced disease burden in this cohort, being highest in patients with RCB III (regardless of recurrence) and lowest in patients with RCB 0/I/II who did not have a breast cancer recurrence within 3 years (Supplementary Fig. S17). These findings suggest that peripheral blood gene expression may be predictive of response in diverse groups of patients with breast cancer.

With the recent approval of immunotherapy in combination with chemotherapy in mTNBC, and promising early results in the neoadjuvant setting, an improved understanding of how chemotherapy reshapes antitumor immunity, both in the tumor and in the peripheral compartment, is needed. Perhaps, the most widely studied marker to approximate antitumor immunity is sTILs, which are predictive of improved NAC response when measured in the primary untreated tumor and are prognostic of good outcomes in the RD of patients lacking pCR. Moreover, sTILs have primarily been a useful biomarker in TNBC as opposed to hormone receptor–positive cancers, but this inference is complicated by the routine use of endocrine-targeted agents in the adjuvant setting, that can impact postsurgical outcomes. However, this likely confounder leaves space for a biological contribution, as numerous differences between molecular and clinical subtypes exist. Even when considering only TNBC, quantification of sTILs is an imperfect biomarker, which does not precisely inform on the immunobiology of the tumor.

In this study, we present a molecular analysis of the TIME in response to NAC in 83 patients with breast cancer, specifically focusing on patients lacking a pCR, as these patients have worse outcomes. We found that changes in tumor immunity seem to be most prevalent in TNBC, often resulting in decreases in expression of immune-related genes. However, an upregulation of immune-related gene expression in tumors following NAC was associated with a strikingly improved outcome after surgery, specifically in TNBC. Of these genes, those involved in cytotoxic effector cells were among the most robustly associated with outcome. Furthermore, we found that cytokines expressed by PD-1^HI CD8⁺ T cells in the peripheral blood were increased dramatically in patients with TNBC following NAC.

Analysis of TCR clonotypes infiltrating into tumors suggested that chemotherapy may preferentially increase the recruitment of new T-cell clones into the tumor, rather than expanding the T cells already present. This effect was consistent with that observed in the peripheral blood, where PD-1^HI CD8⁺ T cells, while highly clonal compared with PD-1^NEG cells, did not substantially change in clonality during NAC; these observations reflect a lack of clonal expansion in response to NAC, as we found in the TIME.

Assessment of peripheral blood represents a unique opportunity to monitor antitumor immunity through minimally invasive means. Using scRNAseq, we identified a population of cytolytic effector T cells in blood that expressed elevated levels of exhaustion/checkpoint genes. A gene expression signature derived from this population was used to test the hypothesis that these highly cytolytic, but potentially exhausted cells may be reflective of an ongoing disease process, and therefore a peripheral approximation of disease burden. This hypothesis was confirmed in two validation sets totaling 88 patients with breast cancer and serves as a proof-of-principle for the use of this signature as a possible biomarker of outcome.

There are several limitations of our study. We intentionally chose to focus on patients lacking a pCR for the tumor study, as these patients have the worst outcomes after surgery, and because of the caveats and difficulty in defining or assessing antitumor immunity after chemotherapy in surgical specimens that lack tumor tissue. Future studies will need to be done to determine whether alterations in immune-related genes are seen over the course of NAC in patients experiencing pCR. Furthermore, our sample size was limited in our exploration of peripheral blood PD-1^Hi CD8⁺ T cells by scRNAseq to 2 patients. Using this method, we identified a unique cluster of cytotoxic cells and derived a gene signature. We used further scRNAseq to identify that this gene signature is expressed predominately by circulating CD8⁺ T cells and NK cells. Therefore, expression of this gene signature in whole blood is not limited to PD-1^Hi CD8⁺ T cells. We tested expression of this gene signature in whole blood of two different cohorts of patients with breast cancer and found that expression was highest in those with RD who recurred or had the largest residual cancer burden. This provides proof-of-concept that peripheral blood gene signatures may predict response, but optimization of signature genes is likely needed to improve prediction.

Importantly, there has been a paucity of studies looking at the effect of chemotherapy on peripheral blood (45), with little data on disease outcomes. This study provides a novel assessment and framework for an improved understanding of how chemotherapy alters antitumor immunity both in the TIME and the peripheral compartment. These data represent a unique opportunity to better understand patient populations most likely to benefit from the addition of immunotherapy to chemotherapy, particularly in the neoadjuvant setting. Furthermore, our findings demonstrating the association of expression of key cytolytic and immune-activation genes in the peripheral blood with presence of RD and recurrence represent a possible biomarker platform. Peripheral gene expression signatures may identify high-risk populations with potentially exhausted T cells and either primary or micrometastatic disease who may benefit from additional immunotherapeutic strategies.

M.L. Axelrod is listed as a coinventor on a provisional patent application on methods to predict therapeutic outcome using blood-based gene expression patterns, that is owned by Vanderbilt University Medical Center, and is currently unlicensed. M.J. Nixon reports grants from NIH/NCI and DoD during the conduct of the study. W.J. McDonnell reports personal fees and other from 10× Genomics (currently a shareholder and employee of 10× Genomics, which occurred after publication) outside the submitted work. S. Loi reports nonfinancial support from Seattle Genetics, Pfizer, Novartis, BMS, Merck, AstraZeneca, and Roche-Genentech [consultant (nonremunerated)]; other from Aduro Biotech, Novartis, GlaxoSmithKline, Roche-Genentech, and G1 Therapeutics [consultant (remunerated, paid to institution)]; and other from Novartis, Bristol Meyers Squibb, Merck, Roche-Genentech, Puma Biotechnology, Pfizer, Eli Lilly, and Seattle Genetics [research funding (to institution)] outside the submitted work; and is supported by the National Breast Cancer Foundation of Australia Endowed Chair and the Breast Cancer Research Foundation, New York. V.M. Jansen reports personal fees from Eli Lilly and Company (past employee) and Mersana Therapeutics (current employee) outside the submitted work. S.M. Tolaney reports grants and personal fees from Genentech/Roche, Merck, Eli Lilly, Novartis, AstraZeneca, Nektar, Nanostring, Pfizer, Immunomedics, Bristol-Myers Squibb, Eisai, Sanofi, Odonate, and Seattle Genetics (research funds to institution; honorarium for consulting/advisory board participation); grants from Cyclacel and Exelixis (research funds to institution); personal fees from Puma, Celldex, Daiichi-Sankyo, Silverback Therapeutics, G1 Therapeutics, Abbvie, Athenex, OncoPep, Kyowa Kirin Pharmaceuticals, Samsung Bioepsis Inc., and CytomX (honorarium for consulting/advisory board participation) outside the submitted work. I. Krop reports grants from Genentech (to institution) during the conduct of the study; personal fees from Genentech/Roche; grants from Pfizer (to institution); personal fees from Bristol Meyers Squibb, Daiichi-Sankyo, Macrogenics, Context Therapeutics, Taiho Oncology, Seattle Genetics, Novartis (for DSMB participation), Merck (for DSMB participation), Celltrion, and AstraZeneca outside the submitted work. A.C. Garrido-Castro reports other from Genentech (institutional funding from Genentech to Dana-Farber Cancer Institute for the clinical trial under which samples included in the submitted work were obtained) during the conduct of the study. T.W. Miller reports grants from Takeda Pharmaceuticals outside the submitted work. I.A. Mayer reports grants and personal fees from Pfizer and Genentech (research grant and advisory board participation) and personal fees from Novartis, Lilly, AstraZeneca, GSK, Abbvie, Seattle Genetics, Immunomedics, Macrogenics, Eisai, and Puma (advisory board participation) outside the submitted work. J.M. Balko reports grants from Genentech/Roche and Incyte Pharma, personal fees from Novartis Inc (consulting), and grants from BMS outside the submitted work, and is listed as a coinventor on a provisional patent application on methods to predict therapeutic outcome using blood-based gene expression patterns, that is owned by Vanderbilt University Medical Center and is currently unlicensed. No potential conflicts of interest were disclosed by the other authors.

M.L. Axelrod: Conceptualization, data curation, formal analysis, investigation, visualization, writing-original draft. M.J. Nixon: Conceptualization, data curation, formal analysis, investigation. P.I. Gonzalez-Ericsson: Formal analysis, writing-review and editing. R.E. Bergman: Data curation. M.A. Pilkinton: Investigation, writing-review and editing. W.J. McDonnell: Formal analysis, writing-review and editing. V. Sanchez: Investigation. S.R. Opalenik: Investigation, project administration. S. Loi: Resources, writing-review and editing. J. Zhou: Formal analysis, methodology. S. Mackay: Formal analysis, methodology. B.N. Rexer: Resources, writing-review and editing. V.G. Abramson: Resources, writing-review and editing. V.M. Jansen: Resources, writing-review and editing. S. Mallal: Resources, writing-review and editing. J. Donaldson: Resources, data curation, writing-review and editing. S.M. Tolaney: Resources, data curation, writing-review and editing. I.E. Krop: Resources, writing-review and editing. A.C. Garrido-Castro: Resources, writing-review and editing. J.D. Marotti: Resources, writing-review and editing. K. Shee: Resources, writing-review and editing. T.W. Miller: Resources, writing-review and editing. M.E. Sanders: Formal analysis, writing-review and editing. I.A. Mayer: Resources, data curation, writing-review and editing. R. Salgado: Resources, writing-review and editing. J.M. Balko: Conceptualization, formal analysis, supervision, funding acquisition, visualization, writing-original draft.

Funding for this work was provided by Susan G. Komen Career Catalyst Grant CCR14299052 (to J.M. Balko), NIH/NCI R00CA181491 (to J.M. Balko), NIH/NCI SPORE 2P50CA098131–17 (to J.M. Balko), Department of Defense Era of Hope Award BC170037 (to J.M. Balko), and the Vanderbilt-Ingram Cancer Center Support Grant P30 CA68485. Additional funding was provided by NIH T32GM007347 (to M.L. Axelrod) and F30CA236157 (to M.L. Axelrod); F30CA216966 (to K. Shee); R01CA200994 (to T.W. Miller) and R01CA211869 (to T.W. Miller); and Dartmouth College Norris Cotton Cancer Center Support Grant P30CA023108 (to T.W. Miller).

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.