The Role of Proliferation in Determining Response to Neoadjuvant Chemotherapy in Breast Cancer: A Gene Expression–Based Meta-Analysis

In spite of remarkable progress in targeted therapeutics in breast cancer, including hormonal and HER2-directed therapies, chemotherapy continues to play a critical role in the management of over one third of all women with breast cancer (1–3). Neoadjuvant chemotherapy (NAC) is associated with improved rates of conversion from mastectomy to breast conservation and also offers prognostic information: in a large meta-analysis, patients who achieved pathologic complete response (pCR) following NAC had improved disease-free and overall survival relative to those with residual disease (4). In addition, NAC provides a window through which to evaluate the intrinsic chemosensitivity of breast cancer, as patients receive no treatments prior to NAC and there is a discrete endpoint—pCR or residual disease (RD).

Gene expression signatures are sets of genes whose expression patterns associate with a specific biological process or phenotype. Over 15 years ago, Perou and colleagues used gene expression profiling to describe molecular subtypes within breast cancer that initially correlated with overall survival and ultimately were shown to reflect sensitivity to NAC (5–8). Studies of gene expression signatures by Desmedt and colleagues over the past decade revealed the importance of proliferation and other biologic processes, including immune activation, in breast cancer prognosis and response to NAC (9–11). Recent prospective data from the TAILORx and MINDACT studies demonstrate the importance of signature-based genomic predictors (OncotypeDX and MammaPrint, respectively) in clinical decision-making for breast cancer (12, 13). While these studies and many others demonstrate the utility of gene expression signatures, most evaluate only one or a small number of signatures, making cross-study comparison of signatures difficult due to sample and methodological differences. Systematically interrogating multiple individual gene expression signatures could provide insight into the complex network of pathways and processes involved in chemosensitivity and resistance (14–16).

In this patient-level meta-analysis, we assembled publicly available mRNA microarray gene expression data from over 1,400 breast cancer biopsies obtained prior to NAC, the largest dataset of its type to date. We collected study-reported response to NAC defined as pCR versus RD as well as patient clinicopathologic characteristics. With this large clinically annotated dataset, we interrogated pathways involved in response to chemotherapy by evaluating and comparing 125 published breast cancer–related gene expression signatures implicated in chemosensitivity and resistance. We demonstrate that proliferation accounts for the vast majority of the predictive power of gene signatures within each subset of breast cancer and identify proliferation-independent signatures predictive of response to NAC, including immune signatures in ER⁺/HER2⁻ breast cancers and signaling pathway signatures in triple-negative breast cancers (TNBC).

Using Preferred Reporting Items for Systematic Reviews and Meta-Analyses (PRISMA) guidelines, two gene expression repositories—the National Center for Biotechnology Information Gene Expression Omnibus (NCBI GEO; ref 17) and the European Molecular Biology Laboratory/European Bioinformatics Institute (ArrayExpress; ref. 18)—were queried for available gene expression datasets; date of last query was May 1, 2015. Using the search terms “breast cancer" and (neoadjuvant OR chemotherapy OR residual), we identified 22 studies with 2,407 patient records (Fig. 1). Only studies with microarray gene expression data from Affymetrix (U133 and U133Plus2.0 arrays) and Agilent platforms were included to limit cross-platform variability. Individual samples were excluded if biopsy was obtained after NAC, if the patient received HER2-directed therapy (e.g., trastuzumab), or if pathologic response data were not available. Duplicate records were removed by careful review of GEO annotation and any records with Pearson correlation coefficient of one. In total, 1,419 records from 17 studies were included for analysis (Table 1; Supplementary Table S1; NCBI GEO accession numbers GSE8465, GSE16446, GSE18728, GSE19697, GSE20194, GSE20271, GSE21974, GSE21997, GSE22093, GSE22226, GSE22358, GSE22513, GSE23988, GSE25066, GSE28796, and GSE32646). Raw Affymetrix microarray data were processed with RMA normalization and Agilent microarray data was subjected to Lowess-normalization and log₂ transformation. Samples were batch-median centered by study to allow cross-study comparison (19, 20). To determine intrinsic breast cancer subtype, the “Bioclassifier” package from Parker and colleagues (21) was used with no additional calibration parameters given that each dataset was relatively balanced with respect to ER status by pathology: 49% and 58% ER-positive for the Affymetrix and Agilent datasets, respectively (Table 2; ref. 7). Receptor status by gene expression was determined using the genefu package in R (22). Triple-negative breast cancer (TNBC) subtype was determined using the TNBCtype tool (23) after re-normalizing data within TNBCs per recommended procedures (24, 25). The ESTIMATE package was used to calculate estimated tumor purity using R version 2.15.3 (26).

Patient-specific data were obtained from NCBI GEO and ArrayExpress, including response to NAC, age, grade, pathologic receptor status (ER, PR, and HER2), and chemotherapy regimen (Table 2). Response was categorized as pCR or RD based on study-specific definitions of endpoint. Most studies defined response by breast and axilla, though four studies reported breast-only response (Table 1 and Supplementary Table S1). Pathologic receptor status was “positive” or “negative” based on study-specific definitions. For ER or PR immunohistochemistry (IHC), 0 was defined as “negative” and 2–3 defined as “positive,” while IHC of 1 was considered indeterminate. A FISH HER2/CEP17 ratio of greater than 2.0 was defined as positive. Chemotherapy was patient-specified in 15 studies and based on the associated publication in two studies.

Gene expression signatures were compiled from published studies describing a variety of processes implicated in chemosensitivity or resistance (11, 15) and from curated databases, including MSigDB (27) and GeneSigDB (28). We included only signatures that were defined from mammary or breast cancer cell lines, breast cancer samples, or had been previously validated in breast cancer expression datasets. Signatures were categorized based on the proposed phenotype characterized by that signature (Fig. 2A). To establish a single value for each signature, we calculated the mean of the gene expression values identified as “up” in each signature based on the original publication. As a sensitivity analysis, we evaluated those signatures that had both “up” and “down” genes (32/125 signatures; 25.6%) and calculated the mean of the “down” genes subtracted from the mean of the “up” genes. Bidirectional signature scores were highly correlated with scores calculated using only “up” genes (median Pearson r of 0.7412). Based on these data, we pursued subsequent analyses using only “up” genes. We also included “random” signatures of varying sizes (5, 10, 50, 100, 500 genes) composed of genes selected randomly from the transcriptome, termed “RANDOM.GENOME” or genes randomly selected from the list of genes present in any signature, termed “RANDOM.SAMPLE.” Signature nomenclature, categorization, reference/PMID, and full gene lists are provided (Supplementary Table S2) as well as code for analyses in R (Supplementary Data).

MCF-10A mammary epithelial cells were cultured in MatriGel as described previously, and total mRNA was previously isolated from triplicate samples and hybridized onto Affymetrix gene chips U133A and U133B (29). Raw gene expression data were RMA processed and then batch median centered. The proportion of proliferating cells was previously determined by detecting total DNA content via FACS at each time point (30).

To determine proliferation-associated versus proliferation-independent signatures, we performed a proliferation adjustment using the method proposed by Venet and colleagues (31). Briefly, a linear model was constructed using the proliferation signature value for each sample fitted to the expression of each gene using the “lm” function in R (version 3.1.3). Each expression measurement was then substituted by the sum of its residual and mean expression across the dataset. This approach was performed using two independent proliferation signatures, the 11-gene PAM50 proliferation index (32) and the 131-gene PCNA signature (31). Venn diagrams were created using BioVenn (33).

To determine if we could improve predictive performance by combining signatures, we developed a proliferation-immune ‘meta-signature’ consisting of 18 genes—11 from the PAM50 proliferation signature (BIRC5, CCNB1, CDC20, NUF2, CEP55, KNTC2, MKI67, PTTG1, RRM2, TYMS, and UBEC2) and 7 from the GeparSixto immune activation signature (CXCL9, CCL5, CD8A, CD80, CXCL13, IGKC, and CD21; refs. 32, 34). To evaluate the performance of this meta-signature, we combined ER⁺/HER2⁻ cases from both Affymetrix and Agilent datasets (total n = 642) and we randomly partitioned the combined ER⁺/HER2⁻ dataset into equal training/validation subsets (n = 321) and then calculated the association of the “Proliferation-Immune” meta-signature and the other 125 signatures via t test. We repeated this approach for 1,000 iterations and then calculated an average FDR P value for each signature as a measure of consistent performance.

All microarray data processing and statistical analyses were performed in R version 3.1.3. Contrasts in patient and tumor characteristics were evaluated using Pearson χ² tests. Correlation among gene expression signatures was calculated using Pearson's coefficient, and hierarchical clustering was performed using average linkage. The association of signatures to continuous and categorical factors was evaluated using Student t test and analysis of variance, respectively. The association of individual genes significantly associated with NAC response was identified using limma package and by t test (35). All calculations of association with response were multiple-testing corrected using the Benjamini–Hochberg procedure for false discovery rate.

This meta-analysis includes microarray expression data from breast cancer biopsies obtained prior to NAC from 1,419 breast cancer samples in 17 studies (Fig. 1, Table 1, and Supplementary Table S1). The data were divided into two independent datasets based on microarray platform for validation purposes: an Affymetrix Dataset (12 studies with 1,033 total samples) and an Agilent Dataset (5 studies with 386 total samples). The two datasets were relatively balanced with respect to age, grade, pathologic ER and HER2 receptor status, and clinical subtype (Table 2). HER2⁺ breast cancers were relatively underrepresented relative to general incidence as we excluded those patients who received HER2-directed therapy, while TNBCs were relatively overrepresented given that patients with TNBC are more likely to receive NAC. The chemotherapeutic agents that these patients received are similar to those agents that are currently administered in most centers. Clinicopathologic characteristics known to be associated with response to NAC—including tumor grade, receptor status, receptor-based subtype, and PAM50 status—were all significantly associated with response to NAC in both the Affymetrix and Agilent datasets (all P < 0.05; Supplementary Table S3).

We evaluated 125 gene expression signatures representing the expression patterns of a diverse set of pathways and cellular processes implicated in chemosensitivity and resistance in the 1,033 breast cancers samples in the Affymetrix dataset (Fig. 2A). Pairwise correlation of the calculated values for each published signature in every sample and subsequent unsupervised clustering demonstrates that independently derived signatures reveal marked redundancy (Fig. 2B). The largest cluster is composed of 46 signatures (46/125 signatures; 36.8%) with five well-established proliferation signatures at the center of the cluster (10, 15, 24, 31, 32). The 46 signatures in this group each demonstrate a strong correlation to the 11-gene PAM50 proliferation signature (ref. 32; Pearson r > 0.40 for all associations; P < 1 × 10⁻⁴⁵) and will be considered “proliferation-associated” by correlation. Although only 11 signatures were designed to represent proliferation or cell cycle (Fig. 2A), the pairwise correlation data suggest that a large proportion of published signatures, irrespective of described signature target nonetheless correlate highly with proliferation. Several other distinct clusters were evident, centering around ER and HER2-related signatures, mesenchymal-related signatures, and immune-related signatures (Fig. 2B). This supports prior evidence of significant redundancy among published gene expression signatures (11, 36).

To explore gene expression–based proliferation signatures in greater detail, we evaluated proliferation gene expression signatures in a 3D culture model of immortalized normal mammary epithelial cells that undergo a programmed proliferation arrest, transitioning from a high proliferation state (30%–35% in G₂–S phase on days 2–5) to a low proliferation state (approximately 5% in G₂–S phase days 12–15; Fig. 3A). We assessed the proportion of proliferating cells daily over a 15-day time course and compared this with five well-established proliferation signatures, calculated from RNA collected in triplicate each day (Fig. 3A; 10, 15, 24, 31, 32). The five proliferation signatures show a strong intersignature correlation (Pearson r > 0.97 for all associations) and each has a strong correlation with proportion of proliferating cells (Pearson r > 0.88 for all associations). These in vitro results indicate that the five proliferation signatures, derived from patient tumor expression data, reproducibly measure cellular proliferation of mammary epithelial cells.

Tumor biopsies often capture nontumor cells that could affect the proliferation score, for example low proliferative fibroblasts or high proliferative immune cells. We calculated tumor purity for each sample using a method that infers the fraction of stromal and immune cells in tumor samples, shown to correlate with tumor purity based on DNA copy number (26). Estimated tumor purity demonstrated no relationship with the 11-gene PAM50 proliferation score (32) across all samples (R² = 0.013; Fig. 3B).

We evaluated the five proliferation signatures described above in 1,419 breast tumors. More highly proliferative tumors were associated with a greater proportion of pCR to NAC in both the Affymetrix and Agilent datasets (Fig. 3C and D; ANOVA P = 1.6 × 10⁻¹³ and P = 3.2 × 10⁻¹¹, respectively). Within individual breast cancer subgroups, the highest proliferative quartile of tumors via the 11-gene PAM50 proliferation index were associated with a greater proportion of pCR in all receptor-based (ER⁺/HER2⁻, HER2⁺, TNBC; Fig. 3E) and PAM50 (Luminal A, Luminal B, Her2-like, Basal-like, Normal-like; Fig. 3F) subgroups. This association was highly statistically significant by stratified χ²: P = 3.0 × 10⁻⁸ and P = 0.003, respectively (Supplementary Table S4A). As a sensitivity analysis, we also compared proliferation above/below median; higher proliferation remained significantly associated with a greater proportion of pCR for both receptor-based subgroups (P = 0.0004) and PAM50 subgroups (P = 0.047; Supplementary Table S4B). Additionally, we evaluated proliferation and pCR rate among the TNBC subtypes from Lehmann and colleagues using proliferation above/below median (24). Five of the TNBC subtypes demonstrated higher pCR rates for highly proliferative tumors (basal-like 2, immunomodulator, luminal AR, mesenchymal, unselected) while two did not (basal-like 1, mesenchymal stem-like). Given the limited number of tumors within each TNBC subtype, it is not surprising that the trend for proliferation correlating with pCR within the TNBC subtypes did not reach statistical significance (stratified χ² P = 0.22; Supplementary Table S4B).

To evaluate the performance of proliferation signatures versus individual genes, we evaluated three individual genes commonly associated with proliferation: MKI67, PCNA, and AURKA. Comparing each individual proliferation gene with each of the five proliferation signatures demonstrates robust correlation across all comparisons (median Pearson r = 0.79; range, 0.59–0.88; Supplementary Fig. S1A). We stratified the expression of each gene into quartiles calculated across subtypes or within subtypes as an attempt to approximate a tissue-based biomarker (e.g., 0, 1+, 2+, 3+). Each individual marker of proliferation was significantly associated with pCR; however, the 11-gene PAM50 signature outperformed (lower stratified χ² P value) each individual gene in comparisons where quartiles are defined within subtypes. In analysis across subtypes, the 11-gene PAM50 signature outperformed MKI67 and AURKA and was comparable to PCNA (P = 4.3 × 10⁻⁵ vs. 1.3 × 10⁻⁵, respectively), suggesting that expression-based proliferation signatures have the potential to be more robust biomarkers than individual genes (Supplementary Fig. S1B and S1C).

We evaluated the significance of association of each of the 125 gene expression signatures with response to NAC (pCR vs. RD; Supplementary Table S5A–S5H). The signatures associated with response were highly overlapping between Affymetrix and validation Agilent datasets: among all breast cancers, 88.9% (48/54) of statistically significant (FDR P < 0.05) signatures in the smaller Agilent dataset were also significant in the Affymetrix dataset. To investigate the potential impact of median center normalization by study, multivariate analysis of signatures with study as a covariate demonstrated that very few signatures that shifted from significant to not significant (or vice-versa): 1.6% in ER⁺/HER2⁻, 5.6% in TNBC, and 4.8% in basal-like (Supplementary Fig. S2A–S2C).

Among receptor-based ER⁺/HER2⁻ breast cancers (Fig. 4A), four of the five most highly associated signatures were proliferation signatures (14, 24, 31, 32) and the fifth a STAT1 signature (37); all five signatures were positively correlated to pCR. For TNBC (Fig. 4B), the five signatures most highly associated with response included a signature of E2F1 activity (14), an expression-based predictor of chemosensitivity (38), a PI3K signature (14), a PTEN deletion signature (39), and the “mature luminal” mammary lineage signature (40); the “mature luminal” signature was positively correlated with RD while the other four signatures with pCR. Thirty signatures were predictive of response in both ER⁺/HER2⁻ and TNBCs, while 16 signatures were predictive in both triple-negative and basal-like breast cancers (Supplementary Fig. S3A). Within each of the PAM50 subgroups, multiple signatures had nominal P values less than 0.05 for association with response to NAC (Supplementary Table S5D–S5H). However, after correcting for multiple testing at a false-discovery threshold of 0.05, signatures remain statistically significant only in basal-like (17 signatures) and normal-like (2 signatures) subgroups. The ranked statistical significance, direction of association, and grouping for all signatures is provided in Supplementary Table S5A–S5H.

To investigate the contribution of proliferation relative to other processes, we compared pCR rates of the five “pure” proliferation signatures with the most highly significant five signatures in the proliferation cluster (“proliferation-associated”), five signatures that were significantly associated with response but anticorrelated to proliferation, and five immune signatures that were “proliferation-independent.” When all breast cancers are stratified into quartiles by signature score, “pure” proliferation and “proliferation-associated” signatures revealed a very similar pCR increase from approximately 10% pCR at the lowest signature level to approximately 35% pCR at the highest signature level with a parallel rise in average proliferation score (Supplementary Fig. S3C and S3D). The five proliferation anticorrelated signatures demonstrate a decline in pCR rate with a parallel—but less dramatic—decline in proliferation score (Supplementary Fig. S3E). Five immune signatures reflected a pCR rate increased with higher immune signature with minimal change in average proliferation score (Supplementary Fig. S3F).

We evaluated proliferation-associated versus proliferation-independent signatures by two distinct approaches: (i) correlation to the 11-gene PAM50 signature and (ii) linear modeling and residualization to adjust total gene expression data for proliferation (31). For the correlation approach, we considered signatures in the proliferation cluster as “proliferation-associated” (Fig. 2B; Pearson r > 0.40 to the 11-gene PAM50 proliferation signature, P < 1 × 10⁻⁴⁵), with the remaining signatures “proliferation-independent.” In the linear modeling approach, those signatures that were significantly associated with response after proliferation normalization were considered “proliferation-independent.” As a sensitivity analysis, we evaluated linear-modeling–based normalization using two separate proliferation signatures, the 11-gene PAM50 proliferation index (32) and the 131-gene PCNA proliferation signature (31). Nonproliferation signatures were largely unaffected by normalization via either signature (Supplementary Fig. S2D).

Using the correlation approach, “proliferation-associated” signatures accounted for the majority of the signatures that achieved statistical significance: 26/42 (61.9%) in ER⁺/HER2⁻ breast cancers, and 38/60 (63.3%) in TNBCs (Fig. 4A and B, red bars). Among signatures that were not “proliferation-associated” by correlation, in ER⁺/HER2⁻ breast cancers, most were related to immune activation (10/16; 62.5%; refs. 10, 14, 20, 24, 34, 37, 41), while in TNBC the signature phenotypes were diverse, including signaling pathways, immune, DNA damage, and mesenchymal phenotype.

Using the linear modeling approach, in ER⁺/HER2⁻ breast cancer nearly all signatures significantly associated with response were “proliferation-associated” (38/42; 90.5%). The four “proliferation-independent” signatures were related to immune activation: a signature of immune genes predictive in the GeparSixto trial (34), an immune signature from TNBC subtypes (24), and two STAT1 signatures (Fig. 4C; refs. 10, 37). As a sensitivity analysis, we also evaluated ER⁺ tumors defined by gene expression (ESR1-positive/ERBB2-negative) and found only three “proliferation-independent” signatures, which were likewise immune-related and overlapping with those in ER⁺/HER2⁻ (10, 24). In luminal A or B breast cancers, 8 of the 21 signatures with nominal P values less than 0.05 for association with response to NAC were likewise immune signatures (Supplementary Table S5D and S5E). We repeated our approach using the 131-gene PCNA signature–normalized data and the results demonstrated an overlapping list of proliferation-independent immune activation signatures (Supplementary Fig. S2E). At the individual gene level, we found that >95% of individual genes significantly associated with response for both ER⁺/HER2⁻ and ESR1⁺/ERBB2⁻ were proliferation-associated using either proliferation normalized dataset (Supplementary Fig. S2F and S2G).

Among TNBCs, 63.3% (38/60) of signatures significantly associated with pCR were “proliferation-associated” using the linear modeling approach. The 22 (36.7%) “proliferation-independent” signatures comprised a variety of processes, including DNA damage/chromosomal instability/mutation, signaling pathways (including PI3K, AKT, Ras, and EGFR), “stemness” signatures, prognostic/predictive signatures, and two signatures of randomly selected genes (Supplementary Table S5C). Signatures made up of random genes from across the genome have been shown to be prognostic in breast cancer previously (31). At the individual gene level, the proportion of significant genes that were “proliferation-independent” in TNBCs was 22.5% to 32.4%, similar to the proportion of signatures (Supplementary Fig. S2F and S2G). Additional sensitivity analyses in PAM50-defined basal-like breast cancer reveal a comparable proportion and overlapping “proliferation-independent” signatures and genes associated with response (Fig. 4D; Supplementary Fig. S2F and S2G).

Based on evidence that immune signatures are a proliferation-independent predictor of pCR in ER⁺/HER2⁻ disease, we developed a proliferation-immune “meta-signature” as a proof-of-concept application of these data. This proliferation-immune “meta-signature” consists of 18 genes—11 that comprise the PAM50 proliferation signature and 7 that comprise the GeparSixto immune activation signature (10, 32). We combined ER⁺/HER2⁻ cases from both Affymetrix and Agilent datasets (total n = 642) and iteratively evaluated the performance of the proliferation-immune “meta-signature” in randomly assigned, equally sized (n = 321) training/validation subsets. Over 1,000 iterations, the proliferation-immune “meta-signature” signature demonstrated a stronger association with pCR than any other individual signature in ER⁺/HER2⁻ breast cancers (Supplementary Fig. S3G).

We compiled gene expression data from 1,419 breast cancer biopsies obtained prior to NAC excluding patients who received HER2-directed therapy, the largest data set of its kind, and calculated 125 gene expression signatures for each sample. We demonstrate that proliferation signatures are highly associated with pCR in all subgroups of breast cancer. These proliferation signatures reflect a broad range of proliferation in an in vitro model of mammary epithelial cell proliferation. We further demonstrate that proliferation differences account for the majority of variation in chemosensitivity in ER⁺/HER2⁻ breast cancer and remains important—but less so—in TNBCs. Importantly, we provide an approach to distinguish proliferation-related from proliferation-independent gene expression signatures and identify immune-related signatures as an important proliferation-independent predictor of sensitivity to NAC in ER⁺/HER2⁻ breast cancer. An 18-gene combined proliferation-immune “meta-signature” performed better than all other signatures in predicting response to NAC in ER⁺/HER2⁻ breast cancers.

Proliferation has long been implicated as a critical mediator of sensitivity to chemotherapy (42, 43) as well as a predictor of long-term outcome (22, 31, 36), but the extent to which proliferation accounts for chemosensitivity relative to other processes is not well understood. Our findings reinforce the overwhelming importance of proliferation in response to chemotherapy in tumors: even within each PAM50 subgroup, which are defined in part by 11 proliferation-related genes, higher proliferative tumors are more likely to achieve pCR. The association of proliferation with chemosensitivity varies based on the subgroup, more prominent in higher proliferative subgroups (e.g., basal-like more than luminal A). These data support multivariate analyses in a prior gene expression study of an overlapping dataset (11) and work from Prat and colleagues on the PAM50 subtypes and Prosigna assay in NAC (44, 45). A recent study retrospectively evaluated proliferation and PAM50 subtype in CALGB 9741, a clinical trial evaluating adjuvant dosing strategies of doxorubicin, cyclophosphamide, and paclitaxel in node-positive breast cancer patients (46). A trend toward better prognosis for dose-dense chemotherapy with higher proliferation score suggests that the proliferation–chemosensitivity association may also be affected by dosing. Most prognostic multigene expression assays used clinically are already heavily weighted toward proliferation-related genes, including OncotypeDX, Mammaprint, and Prosigna (7, 47, 48). As novel signatures are defined, reported, or brought into clinical practice, demonstrating predictive or prognostic significance above and beyond existing proliferation signatures could serve a benchmark for identifying novel, clinically relevant signatures.

In this study, using a linear modeling approach to identify proliferation-associated signatures and genes we find that the majority of the predictive capacity for response to neoadjuvant chemosensitivity within ER⁺ breast cancer is associated with proliferation: 90% of signatures and >95% of individual genes significantly associated with response to NAC can be attributed to proliferation. The striking degree to which proliferation accounts for predictive capacity of gene expression is supported by stratifying into the ER⁺ PAM50 subtypes, luminal A and B, which are distinguished primarily based on proliferation. No signatures remained statistically significant for association with response to NAC in either luminal A or B after multiple-testing correction. Within TNBCs, proliferation remains important but accounts for a lesser extent of variation in response: approximately 65% of signatures and individual genes associated with response to NAC are proliferation associated by our linear modeling approach. The smaller proportion of “proliferation-associated” signatures in TNBCs may reflect less variation in proliferation among TNBCs relative to other breast cancer subtypes (24). The differences between ER⁺/HER2⁻ and TNBC, as well as among PAM50 subtypes, suggest that subgroup-specific genomic predictors are more likely to provide optimal prediction of chemosensitivity and resistance.

A crucial next step is to understand what affects chemosensitivity outside of proliferation. Our data suggest that immune activity plays a crucial role in response to NAC as multiple distinct immune signatures are associated with response to NAC, particularly in ER⁺/HER2⁻ breast cancers, an observation that may have been masked in historical analyses by effects of proliferation (16). These signatures appear to primarily reflect immune activation, as they are comprised of activating cytokines, chemokines, and immune response genes (Supplementary Table S1). In transcriptional analyses of tumor biopsies, it is difficult to discern whether these signatures reflect tumor cell–intrinsic expression differences or, alternately, transcriptional programs of infiltrating lymphocytes. Tumor-infiltrating lymphocytes (TIL), evidence of immune activation, have been shown to be a predictive marker of response to NAC primarily in HER2⁺ and TNBC breast cancer, but to date the association has been less robust in ER⁺/HER2⁻ breast cancers (49, 50). We are currently evaluating the underlying cause or driver of these immune signatures.

From a clinical perspective, predictive markers of pCR are needed as we increasingly consider patients for NAC. This study supports expanded use of expression-based proliferation signatures, given our data regarding suggesting enhanced predictive capacity of proliferation signatures relative to individual proliferation gene expression, two prior studies demonstrating greater prognostic value of the 11-gene proliferation score than tissue-based Ki-67 staining (32, 51), and ongoing debate regarding the reproducibility of immunohistochemical proliferation markers in clinical care (52–55). The combination of a proliferation signature and an immune activation signature into a single “meta-signature” demonstrated a higher predictive capacity than any other individual signature in ER⁺/HER2⁻ breast cancers. More robust predictors for ER⁺/HER2⁻ breast cancers would be particularly useful to patients with “intermediate” scores in proliferation-based assays like OncotypeDX or MammaPrint.

Our findings implicate many processes other than proliferation in neoadjuvant chemosensitivity in TNBC and basal-like breast cancer. Sequencing of paired TNBC samples before and after NAC revealed few recurrent chemoresistance-related mutations, suggesting that mutations alone do not explain chemoresistance, potentially implicating transcriptional differences (56). Proliferation-independent signatures predictive of response in our analysis include DNA damage/chromosomal instability, signaling pathway activation (PI3K, AKT, Ras, and EGFR), and “stemness” phenotypes. Signaling pathway activation signatures may appear paradoxical, as these are known to induce proliferation, but this analysis suggests that proliferation-independent components of these pathways also drive an association with response to chemotherapy in TNBCs. The phosphoinositide 3-kinase (PI3K) pathway is particularly of interest in TNBCs given evidence of activating PIK3CA mutations in the luminal AR subtype as well as pathway activation frequently by PTEN loss in nonluminal TNBCs (24, 57, 58). Multiple immune signatures were associated with pCR in TNBCs; however, none remained statistically significant after proliferation adjustment. Despite the absence of a robust association between proliferation and immune signatures in TNBCs (Pearson r² < 0.08 for all associations, data not shown), we are investigating potential links between DNA damage, proliferation, and immune signatures as a possible explanation. Collectively, these data add to our growing understanding of the particularly complex nature of chemosensitivity in TNBC and suggest that stratifying TNBCs based on gene expression signatures may improve prediction of response to NAC in addition to existing expression-based tools/subgroups that rely heavily on proliferation.

Although pooled analyses of clinical trials have not validated the proportion of patients who achieve a pCR as a surrogate endpoint for event-free or overall survival, NAC remains valuable in clinical practice and as a research platform. Predictive markers of pCR are needed as we increasingly consider patients for NAC (4). These data support consideration of the expanded use of expression-based proliferation signatures, possibly in conjunction with immunohistochemical proliferation markers (e.g., mitotic index, Ki-67, or PCNA; refs. 53–55). In ER⁺/HER2⁻ disease, incorporating markers of immune activation may provide more robust gene expression–based predictors of response and help discern likelihood of response for those with “intermediate” proliferation. In TNBC, this study provides candidate cellular processes that could be integrated into a unified tool or stratification approach to improve chemotherapy response prediction in this challenging subgroup.

Our study should be interpreted in the context of several limitations. The inclusion criteria of the 17 trials varied and definition of pCR was primarily, but not exclusively, both breast and axillae. In addition, we defined response as pCR versus RD but analyzing by Miller-Payne (59) or RCB (60) may add information, particularly as there is evidence that RCB-1 may have similar long-term outcomes as RCB-0 (60). Although the chemotherapy agents used in the studies included are similar to those we use for standard of care today, there was variability in regimen and dosing. We focused our analysis exclusively on response to NAC and cannot extrapolate to the adjuvant or metastatic settings, although those analyses are ongoing. There are also general limitations to gene expression data, including failure to capture protein level effects (e.g., apoptotic/BH3 proteins).

In this gene expression-based meta-analysis, proliferation differences account for the majority of predictive capacity of gene expression signatures in neoadjuvant chemosensitivity for ER⁺/HER2⁻ breast cancers and, to a lesser extent TNBCs. In ER⁺/HER2⁻ breast cancer, proliferation-independent signatures predictive of response were all related to immune activation. In TNBC, sensitivity and resistance to NAC is complex and associated with a diverse set of proliferation-associated and proliferation-independent signatures.

No potential conflicts of interest were disclosed.

Conception and design: D.G. Stover, W.T. Barry, L.M. Selfors

Development of methodology: D.G. Stover, W.T. Barry, L.M. Selfors

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): D.G. Stover

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): D.G. Stover, J.L. Coloff, W.T. Barry, J.S. Brugge, L.M. Selfors

Writing, review, and/or revision of the manuscript: D.G. Stover, J.L. Coloff, W.T. Barry, J.S. Brugge, E.P. Winer, L.M. Selfors

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): L.M. Selfors

Study supervision: J.S. Brugge, L.M. Selfors

Financial support was provided by Susan G. Komen for the Cure PDF14299961 (D.G. Stover), Breast Cancer Alliance FS45312 (J.S. Brugge), and Breast Cancer Research Foundation (J.S. Brugge).

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.