Analysis of Pre- and Posttreatment Tissues from the SWOG S0800 Trial Reveals an Effect of Neoadjuvant Chemotherapy on the Breast Cancer Genome

The S0800 (NCT00856492) clinical trial was a three-arm neoadjuvant (i.e., preoperative) study that randomized patients with stage II–III breast cancer to either (i) weekly nab-paclitaxel and bevacizumab followed by dose-dense doxorubicin/cyclophosphamide (ddAC), (ii) nab-paclitaxel followed by ddAC, or (iii) ddAC followed by nab-paclitaxel. The study included both estrogen receptor (ER)-positive and ER-negative patients. The trial demonstrated that bevacizumab increased pathologic complete response (pCR, defined as complete eradication of all invasive cancer from the breast and lymph nodes) from 21% to 36%; P = 0.019) but chemotherapy sequence in the non-bevacizumab arms did not influence efficacy (1). Pretreatment (i.e., baseline) core needle biopsy and posttreatment surgically resected tissues were prospectively collected for biomarker analysis. We previously reported that high baseline tumor-infiltrating lymphocyte (TIL) count and programmed death ligand-1 (PD-L1) protein expression in stromal cells were associated with higher pCR rates in all treatment arms and that TIL counts, but not PD-L1 expression, decreased significantly after treatment (2). We also examined mRNA expression of 750 immune-related genes corresponding to 14 different immune cell types and a broad range of immune functions in matched pre- and posttreatment samples. At baseline, in addition to higher TIL counts and PD-L1 expression, high expression of chemoattractant cytokines (e.g., CCL21 and CCL19) and cytotoxic T-cell markers were also associated with higher pCR rate, whereas high expression of stromal genes (e.g., VEGFB, TGFB3, PDGFB, FGFR1, and IGFR1), mast cell, and myeloid inflammatory cell metagenes, stem cell–related genes (CD90, WNT11, and CTNNB1) and CX3CR1, and IL11RA were higher in cancers that did not achieve a pCR (3). In posttreatment residual cancer samples, most immune gene expression decreased but IL6, CD36, CXCL2, and CD69 expression increased compared with baseline. The goal of this analysis was to perform whole-exome sequencing (WES) and assess the somatic mutation landscape of the tumors before and after neoadjuvant chemotherapy.

Of the 215 patients accrued to the S0800 trial, 134 patients had pretreatment and 63 patients had posttreatment formalin-fixed, paraffin-embedded (FFPE) tissues with informed consent for research, including 60 patients with paired tissues. Patients who had any viable residual invasive cancer after chemotherapy, regardless of clinical response, were categorized as residual disease (RD). Twenty-nine pretreatment samples (22 RD and 7 pCR) and nine posttreatment samples with greater than 10% tumor cell content were available for WES (Supplementary Fig. S1). Demographic and disease characteristics of the WES population and use of tissues are shown in Table 1. The current biomarker study was conducted in accordance with U.S. Common Rule of human subject research. All patients signed informed consent including permission for biomarker analysis of their tissues. The analysis was conducted with approval by the NCI and the Yale University Human Investigations Committee (i.e., institutional review board).

DNA was isolated from 5–7 μm FFPE tissue sections by AllPrep RNA/DNA FFPE Extraction Kit (Qiagen) and PreCR DNA Repair Kit (New England Biolabs). Genomic DNA (1 μg) was sheared to a mean fragment length of 220 bp using the Covaris E210 instrument, purified by Magnetic AMPure XP Beads (Beckman Coulter) and labeled with 6-base barcode during PCR amplification. Exomes were captured using the IDT xGen Exome Research Panel v1.0. Libraries were sequenced on Illumina HS4000 Illumina instrument using 74-base pair paired-end reads by multiplexing four tumor samples per lane to sequence to a median coverage of 174×, 98% of exonic bases passing 30× coverage. Matching normal tissues were not available for DNA sequencing. Sequence data are deposited under dbGAP accession number phs001883.v1.p1.

Reads were filtered by Illumina CASAVA 1.8.2 software, and aligned to the human reference genome (GRCh37) by Burrows-Wheeler Aligner v0.7.5a and PCR duplicates were removed by MarkDuplicates algorithm. We performed local realignment around putative and known insertion/deletion (INDEL) sites using RealignerTargetCreator (Genome Analysis Toolkit: GATK v3.1.1) and applied base quality recalibration using GATK. We used MuTect v.1.1.4 and Strelka v.1.0.14 to identify somatic single-nucleotide variants (SNV) and INDELs, respectively.

WES data from seven posttreatment biopsies from patients that experienced a pCR, and did not have any cancer cells, were pooled to serve as reference normal cohort for somatic variant calling by MuTect and Strelka as described previously (4, 5). Samples were combined after being down-sampled to 14% of their total reads to retain a similar library size where each of the seven samples were represented equally. To minimize misclassifying germline variants as somatic events, we filtered variants with total coverage <20 and variants that were present in at least five of the breast normal samples in The Cancer Genome Atlas. We also excluded from further analysis, variants that we considered likely to be germline because they were listed in any of the following databases: dbSNP, ESP6500, 1000Genome, or Exac01. In addition, we filtered out a variant found in both the pre- and posttreatment biopsies of a patient if either called variant has a minor allele frequency of >0.40 or total coverage <20, and we recognize that this step likely has also removed genuine somatic mutations. We considered all remaining variants potentially somatic. Recurrent (N ≥ 5 cases) annotated variants in COSMIC v64 and Clinvar (http://www.ncbi.nlm.nih.gov/clinvar/) were white listed. A variant was designated as high functional impact (HFI) if it was either an indel with a predicted deleterious effect (frameshift deletion, frameshift insertion, stop gain, or stop loss) or if more than 3 of 5 functional predictors (SIFT, PolyPhen, LRT, MutationTaster, and PhyloP) predicted a deleterious effect. For the majority of the analysis, mutations were aggregated across the 50 biological hallmark pathways downloaded from MSigDB (6).

Mutational signature weights for the Cosmic signatures were estimated using the deconstructSigs R package (7). Signatures with zero weight were discarded, and Wilcoxon rank-sum was used to test signature weights between groups. Bonferroni correction was used to correct for the number of nonzero signatures tested.

Copy number–neutral regions were identified using SynthEX with the panel of pCR residual disease samples used in SNV mutation calling (8). Briefly, SynthEX identifies a representative panel of normal copy number DNA segments whose read depth ratios across the exome compared with the tumor sample have the least variance. This method assumes that the majority of the exome even in cancer samples is in a copy-neutral state. We use the SynthEX K-nearest neighbor approach with a K of 3. Segments with a designated integer copy number between two and three were considered copy neutral. All mutations found in copy-neutral regions that are shared in both pre- and posttreatment biopsies were collected to analyze variant allele frequency (VAF) changes.

To correct for changes in tumor cellularity due to therapy and sampling variation that would influence VAF changes between pre- and posttreatment samples, we assessed changes in VAF as follows: for each patient, mutations were ordered by the change in VAF between pre- and posttreatment mutations and binned by 10% percentiles. We then examined the overall selective pressure exerted by treatment on a given hallmark pathway across the study cohort by looking at the median percentile bin membership across all HFI mutations that effect gene members of that pathway. P values were obtained through resampling (without replacement) of VAF deltas for each mutation within a patient. This resampling was performed 1,000 times for empirical P value calculations.

RNA was isolated from 5-μm FFPE sections using the Qiagen RNeasy FFPE kit and 100 ng total RNA was hybridized to the NanoString PanCancer IO 360 code to quantify the expression of 750 immune-related and 20 housekeeping genes as described previously (3). The nSolver 2.6 software was used to normalize expression values and immune cell type scores were defined as the average log-transformed expression value of the cell type–specific gene sets (9). A total TIL signature was calculated as the average of all log (2) normalized cell type expression scores (excluding dendritic, regulatory T cells, and mast cells) following the manufacturer's recommendation.

Pretreatment samples (n = 29) were sequenced to a median depth of coverage of 174×, with 98% of exonic bases passing 30× coverage. A median of 340 somatic mutations were detected per sample, and a median of 54 mutations were annotated as HFI in each sample.

Overall, patients that experienced a pCR had a similar number of mutations to those who had RD (median 225 mutations, range 69–436 in pCR vs. median 354, range 210–612 in RD, Wilcoxon rank sum, P = 0.079). No correlations were detectable between total TIL abundance and overall mutation rate for the 23 patients with both data types available (Spearman rho = 0.15; P = 0.49). No individual genes demonstrated significant differences in HFI mutation frequency between the pCR and RD cohorts (Supplementary Fig. S2).

Several known biological pathways were affected by HFI mutations in the majority (>60%) of the cohort, including the p53 pathway, E2F target pathway, mTOR signaling pathway, myogenesis, apical junction, mitotic spindle, and complement pathway (Fig. 1). However, no pathways showed significantly higher mutation frequency in cases with pCR after adjusting for multiple testing.

We also examined mutational signatures present in the pretreatment cohort by deconvoluting the frequency of the 96 different possible trinucleotide substitutions against known signatures of mutation patterns (ref. 10; Fig. 2). Notably, 34% of the samples had detectable presence of cosmic signature 3, a signature closely associated with BRCA-mediated homologous recombination deficiency. Contribution of signature 3 to overall mutation signature spectrum was significantly higher in patients who achieved pCR with chemotherapy (median weight 24%, range 0%–38% in pCR vs. median weight 0%, range 0%–19% in RD, Wilcoxon rank sum, Bonferroni-adjusted P = 0.007).

Because of the diminished tumor cellularity of cancers following neoadjuvant therapy, comparing pre- and posttreatment sequencing results is challenging. Often large-scale differences in tumor cellularity are seen, which is the goal of therapy (11). However, the most interesting pre-and posttreatment genomic comparisons may be those that involve tumors with substantial amounts of residual invasive cancer. These cancers demonstrated treatment resistance and the higher posttreatment tumor cellularity also facilitates comparisons between pre- and posttreatment genomes. In our matched cohort, we could identify only nine of 29 samples with greater than 10% posttreatment tumor cellularity (Supplementary Fig. S3), and these were used to compare the tumor genomic landscape before and after therapy. Posttreatment tumors were sequenced to a median depth coverage of 164×, with 98% of exonic bases passing 30% coverage. A median of 283 somatic mutations were detected per sample (not significantly different from pretreatment mutation load), with a median of 74 HFI mutations. Even at the pathway level, pre- or posttreatment exclusive mutations were sparse, with few recurrent pathway alterations detectable (Fig. 3).

We also examined the change in VAF of HFI mutations between paired baseline and posttreatment samples to assess whether selection pressure is evident for any variant. We rank-ordered VAF changes of all HFI and non-HFI mutations detected in copy-neutral regions in each paired sets of samples. The percentile delta in ranking order between pre- and posttreatment samples reflects the magnitude of VAF change and suggests emergence (or depletion) of a cancer cell clone harboring the variant under treatment pressure. For example, a percentile delta of 90 would indicate a mutation being in the top 10% largest VAF increases from pretreatment VAF to posttreatment VAF suggesting that the mutation mediates treatment resistance because it is being selected for over most other variants observed in that cancer. This rank-based comparison was performed at variant, gene, and pathway levels. There were no recurrent highly selected individual variants in this small cohort. We also failed to find significantly selectively mutated single genes. At pathway level, “E2F targets” pathway mutations (5/9 patients affected, pathway size = 200 genes; Fig. 3) were significantly enriched in posttreatment residual cancer (median percentile delta = 80; permutation P = 0.027; Tables 2 and 3). Mutations in “G₂–M checkpoint” pathway (4/9 patients affected, pathway size = 200 genes; Fig. 3) also showed a significant VAF enrichment in postchemotherapy tissues, (median percentile delta = 80; permutation P = 0.048; Tables 2 and 3). We also observed a statistically significant VAF depletion of HFI mutations in the “myogenesis pathway” in residual cancer suggesting the cell clones harboring these variants were eradicated by chemotherapy (median percentile delta = 15; permutation P = 0.021; Tables 2 and 3).

We examined whole-exome sequences from breast cancer tissues obtained before and after 20 weeks of chemotherapy with or without bevacizumab. We could not identify any recurrent variant or gene-level mutation at baseline that was associated with pathologic response to therapy. This is not entirely surprising considering the small sample size and known sparsity of recurrent mutations, except TP53 and PIK3CA, in breast cancer. Even in substantially larger studies, PIK3CA mutations were not associated with response to preoperative anthracycline- and taxane-based chemotherapies in breast cancer (12). On the other hand, loss-of-function mutations in TP53 have been reported to occur more frequently in patients with pCR, but this association is attributed to the greater frequency of TP53 mutations in triple-negative breast cancers and is not statistically significant after adjustment for ER status (13–15). Other mutations are too infrequent for statistical analysis in currently available small to moderate sized studies. Individually rare mutations can be aggregated at pathway level to reflect genomic disturbance of a biological process. Few studies addressed association between mutations at biological pathway level and response to neoadjuvant chemotherapy in breast cancer. In a previous study in HER2- amplified breast cancers, we found that while no single-gene mutations were predictive of pCR to neoadjuvant chemotherapy and HER2-targeted therapy (lapatinib or trastuzumab), pathway-level mutations were statistically significantly associated with pCR to lapatinib (RhoA pathway) and resistance to trastuzumab (PI3K-related gene network; ref. 4). In this study, we also observed that several biological pathways were affected by HFI mutations (Fig. 1) but no pathways showed significantly higher mutation frequency in cases with pCR after adjusting for multiple testing, which may be due to the limited power of this study caused by the small sample size.

Different mutational processes result in different and unique patterns of single-nucleotide alterations that can be sorted into mutational signatures. We examined whether any of the known mutation signatures are associated with greater chemotherapy sensitivity and found that signature 3, that is associated with BRCA-mediated homologous recombination deficiency, was significantly higher in patients who achieved pCR even after adjusting for multiple comparisons. This observation is consistent with findings of the GeparSepto trial (NCT01583426), that used a similar taxane- and anthracycline-based neoadjuvant chemotherapy regimen and also reported significantly higher pCR rate in tumors with mutation signature 3 (16). Several studies demonstrated that patients with germline BRCA1 or BRCA2 mutations have increased pCR rates after anthracycline-based chemotherapy regimens compared with BRCA wild-type cancers, even among triple-negative breast cancers (17, 18). A BRCA deficiency transcriptional signature is also associated with higher chemotherapy sensitivity, even in the absence of germline BRCA mutations (19). Collectively, these results indicate that disturbances in BRCA-related DNA repair functions in early stage breast cancer confer increased chemotherapy sensitivity.

We also examined the differences in DNA sequence alterations in paired pre- and posttreatment tissues. No single-variant, or gene-level mutations were enriched in posttreatment samples, however, we observed statistically significant enrichment of mutations in the “E2F targets” and “G₂–M checkpoint” pathways in residual cancer samples. Cancers that survived neoadjuvant chemotherapy frequently had alterations in these pathways, but different genes were affected in different patients. The E2F targets include a coexpression network of genes that regulate cell-cycle progression and are targeted by the E2F transcription factors. The G₂–M checkpoint pathway genes also regulate cell proliferation (6). While it has been extensively documented that breast cancers with high proliferation rate also have greater chemotherapy sensitivity compared with less proliferative tumors (20, 21), sustained high proliferation after chemotherapy in residual cancers is consistently associated with worse overall survival (22, 23). Our data suggest that some of the sustained proliferation after chemotherapy may be due to mutations in regulatory genes acquired during treatment or selected for by the treatment. We also detected significant depletion of mutations in the “myogenesis pathway” suggesting the cells harboring these variants were effectively eradicated by therapy.

This study has important limitations. Our sample size is small, and this limits the power of our analyses. The particularly small sample size of paired pre- and posttreatment tissues highlight the challenge of analyzing clinical trial material when highly effective chemotherapies are used that result in major changes in tumor cellularity. Also, because matched normal tissues were not collected during the trial, we had to apply aggressive variant filtering and used the cases with pCR as the “normal cohort” for variant calling. Previous studies have shown that using this cohort-normal approach to variant calling can retain an acceptable level of specificity in mutation calling, but at the cost of sensitivity (5). Despite these limitations, our results are among the first WES efforts to assess genomic changes during neoadjuvant chemotherapy in breast cancer. Overall, the results suggest that genomic disturbances in BRCA-related DNA repair mechanisms, reflected by a dominant mutational signature 3, confer increased chemotherapy sensitivity and cancers that survive neoadjuvant chemotherapy frequently have alterations in cell-cycle–regulating genes.

R.L. Powles and C. Hatzis are employees/paid consultants for Bristol-Myers Squibb. L. Pusztai reports receiving commercial research grants from AstraZeneca, Merck, H3Bio, Clovis, Novartis, and Roche Genenetech, and other commercial research support from AstraZeneca, Seattle Genetics, and Bristol-Myers Squibb. No potential conflicts of interest were disclosed by the other authors.

The content is solely the responsibility of the authors and does not necessarily represent the official views of the NIH.

Conception and design: R.L. Powles, Z. Nahleh, C. Hatzis, L. Pusztai

Development of methodology: R.L. Powles, V.B. Wali, L. Pusztai

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): V.B. Wali, W.E. Barlow, Z. Nahleh, L. Pusztai

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): R.L. Powles, X. Li, W.E. Barlow, C. Hatzis, L. Pusztai

Writing, review, and/or revision of the manuscript: R.L. Powles, V.B. Wali, X. Li, W.E. Barlow, Z. Nahleh, A.M. Thompson, A.K. Godwin, L. Pusztai

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): R.L. Powles, A.M. Thompson, A.K. Godwin, L. Pusztai

Study supervision: Z. Nahleh, A.M. Thompson, L. Pusztai

Research reported in this article was supported by the NCI of the NIH under Award numbers CA180888, CA180819, and CA180826; and in part by Genentech (Roche), Abraxis BioScience (Celgene), Helomics, The HOPE Foundation, a Susan Komen Foundation Leadership Award (to L. Pusztai), and grants from the Breast Cancer Research Foundation (to L. Pusztai and C. Hatzis).

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.