Molecular stratification of early breast cancer identifies drug targets to drive stratified medicine

Many women with hormone receptor-positive early breast cancer can be managed effectively with endocrine therapies alone. However, additional systemic chemotherapy treatment is necessary for others. The clinical challenges in managing high-risk women are to identify existing and novel druggable targets, and to identify those who would benefit from these therapies. Therefore, we performed mRNA abundance analysis using the Tamoxifen and Exemestane Adjuvant Multinational (TEAM) trial pathology cohort to identify a signature of residual risk following endocrine therapy and pathways that are potentially druggable. A panel of genes compiled from academic and commercial multiparametric tests as well as genes of importance to breast cancer pathogenesis was used to profile 3825 patients. A signature of 95 genes, including nodal status, was validated to stratify endocrine-treated patients into high-risk and low-risk groups based on distant relapse-free survival (DRFS; Hazard Ratio = 5.05, 95% CI 3.53–7.22, p = 7.51 × 10−19). This risk signature was also found to perform better than current multiparametric tests. When the 95-gene prognostic signature was applied to all patients in the validation cohort, including patients who received adjuvant chemotherapy, the signature remained prognostic (HR = 4.76, 95% CI 3.61-6.28, p = 2.53× 10−28). Functional gene interaction analyses identified six significant modules representing pathways involved in cell cycle control, mitosis and receptor tyrosine signaling; containing a number of genes with existing targeted therapies for use in breast or other malignancies. Thus the identification of high-risk patients using this prognostic signature has the potential to also prioritize patients for treatment with these targeted therapies.


INTRODUCTION
Despite significant improvements in the treatment of early estrogen receptor positive (ER+) breast cancer, there are ongoing clinical challenges. Targeted anti-endocrine therapies have reduced mortality over the last 30-40 years, 1, 2 but ER+ disease, which comprises 80% of breast cancers, still leads to the majority of deaths from early breast cancer. 3 Multiparametric gene assays are used increasingly to guide clinical treatment decisions. 4 Most prognostic tests provide an estimate of relapse risk following the treatment for ER+ breast cancer, but still lack predictive value for novel targeted treatment options. 2,4 These multiparametric tests, which include OncotypeDx® (Genomic Health Inc.), 5,6 Prosigna ™ (NanoString Technologies, Inc.), [7][8][9] MammaPrint® (Agendia Inc.), 10,11 Breast Cancer Index (BioTheranostics Inc.), 12,13 and EndoPredict (Sividon Diagnostics GmbH), 14 all provide broadly similar clinical utility. 15,16 Although each is derived from RNA abundance studies, there are surprisingly few overlapping genes between different RNA signatures. 17 Prat et al., demonstrated in silico that combined signatures may more accurately predict outcome; leading to greater clinical significance. 18 Nonetheless, despite a decade of development of multiple residual risk signatures, progress towards stratified or targeted medicine has not been markedly accelerated by these tests. None of the existing tests have identified actionable targets which might form the basis for the next generation of stratified medicine approaches.
Using the Tamoxifen and Exemestane Adjuvant Multinational trial (TEAM) pathology cohort, 19 comprised of 3825 hormonereceptor positive (ER+ and/or PgR+) cases and including 477 (13%) HER2-positive cases, we have discovered and validated a novel 95-gene signature of residual risk, which has the potential to markedly improved risk stratification in the context of endocrinetreated patients. Moreover, this gene signature has also revealed potentially druggable targets, thus moving closer to the vision of stratification to targeted therapies for such high-risk patients.

RESULTS
The RNA abundance profiles of all genes were generated for 3825 patients. Of patients who had complete therapy information, 2549 were treated with endocrine therapies alone, while 1275 also received adjuvant chemotherapy. The endocrine-treated only patients were divided into a 576-patient training cohort (n = 67 events), and a 1973-patient validation cohort (n = 253 events), which was used for signature discovery and validation, respectively. To test the prognostic ability of the signature, which was trained and validated in the endocrine-treated patients, to patients who were treated with adjuvant chemotherapy, the signature was then modeled against all patients in the validation cohort and adjusted for adjuvant chemotherapy (n = 3035). The median follow-up in each cohort was 7.51 and 6.21 years, respectively. The clinical characteristics of the endocrine-treated training and validation cohorts are described in Table 1. The clinical characteristics of the entire cohort of 3825 patients are summarized in Supplementary Table 1. High tumor grade, nodal status, pathological size, and HER2 IHC status were univariately prognostic in both training and validation cohorts (Table 1 and  Supplementary Table 1).
Identification and validation of a residual risk signature following endocrine treatment Univariate assessment of the original gene list of 165 genes identified 95 genes which were prognostically significant in the endocrine-treated only patients (Supplementary Table 2, Supplementary Fig. 1). The 95 genes were aggregated into functional modules and used to train a residual risk model. Modeling of the module dysregulation (MDS) scores (as described in the Supplementary Methods) generated from these 95 genes, with and without clinical covariates, resulted in a final refined signature that included nodal status as the only clinical covariate (Fig. 1). This risk model was found to be comparable in the training cohort when the 95-gene signature was used, without clinical covariates (HR high = 4.05, 95% CI 2.25-7.3, p = 3.28 × 10 −6 ; 10-fold cross validation ) and when clinical co-variates such as age, tumor grade, pathological tumor size, and nodal status were included (HR high = 2.74, 95% CI 1.61-4.65, p = 2.06 × 10 −4 ; 10-fold cross validation) ( Supplementary Fig. 2). When dichotomized around the median and applied to the validation set, the resulting 95gene signature was a robust predictor of DFRS following endocrine treatment (HR high = 5.05, 95% CI 3.53-7.22, p = 7.51 × 10 −19 , Fig. 1a). As with the training set, similar results were obtained when all clinical covariates were included in the model of the validation cohort (HR high = 5.56, 95% CI 3.85-8.03, p = 5.75 × 10 −20 , Supplementary Fig. 3). When samples were split into quartiles (Fig. 1b), the signature identified patients were at a very low risk (<5% DRFS at 10 years). The continuous risk scores from this signature were directly correlated with the likelihood of recurrence at 5 (Fig. 1c) and 10-years (Fig. 1d), with a higher risk score associated with a markedly higher likelihood of a metastatic event.
Performance of the 95-gene signature of residual risk in the presence of adjuvant chemotherapy To determine whether the 95-gene residual signature continued to be prognostic amongst patients who also received adjuvant chemotherapy, the model was applied to all patients in the validation cohort (with and without chemotherapy), but stratified to chemotherapy (Supplementary Fig. 4A and B). The results showed that the 95-gene signature was still prognostic in this subset of patients (HR high = 4.7, 95% CI 3.61-6.28, p = 2.53× 10 −28 ). Stratifying according to adjuvant chemotherapy showed no difference in the DRFS between patients defined as low or high risk by the signature (Supplementary Fig. 4C).
Performance of the 95-gene signature of residual risk when adjusted for HER2 status To determine whether the 95-gene residual risk signature remained prognostic in both HER2-positive and HER2-negative patients, the model was applied to patients in the validation cohort who did not receive any additional adjuvant chemotherapy and results stratified by HER2-status ( Supplementary Fig. 5). When the model was applied to all patients and stratified by HER2-status ( Supplementary Fig. 5A), patients identified as low-risk by the 95gene signature, showed no significant difference in DRFS between HER2-positive or HER-2 negative patients (p = 0.78). Similarly, for patients identified as high-risk, we observed no statistically significant difference in DRFS between HER2-positive or HER2negative patients (p = 0.09), although we did observe that HER2positive patients showed a trend for worse outcome. Overall, the signature can differentiate high-risk from low-risk individuals within either HER2-positive (HR = 5.17; 95% CI: 1.25-21.38; p = 0.023) or HER2-negative (HR = 4.75; 95% CI: 3.23-6.97; p = 2.01 × 10 −15 ) patient subsets.
Performance of the 95-gene signature to multiparametric tests Using the NanoString RNA abundance data, risk scores from current multiparametric test were generated and are summarized in Fig. 2a and Supplementary Table 3, along with known prognostic clinical factors. Molecular intrinsic subtyping results are also shown (Fig. 2a). While there exists a common group of  Table 4). When compared to the risk scores generated based on the commercial tests (Fig. 2b, Table 2), our 95-gene signature performed better than these multiparametric tests, with an Area Under the Curve (AUC) of 0.76. The differences in AUC between the commercial tests and the 95-gene risk score were found to be statistically significant (  Supplementary Fig. 6), between patients at low or high risk for recurrence.
Identification of druggable targets in the 95-Gene signature and implications for stratified precision medicine Six significant network modules were identified using the Reactome Functional Interaction (FI) tool, comprising 52 of 95 genes in the signature (Fig. 3a, Table 4). Modules 1, 3, and 4  Table 5) within the modules showed that all genes within Modules 1, 3, and 4 were more highly expressed among patients classed as high-risk (Supplementary Table 5; Wilcoxon rank-sum test). These differences were found to be statistically significant (Supplementary Table 5). As individual modules, they were statistically significant predictors of outcome (Fig. 3b) Using the Integrity Compound Search (Thomson Reuters) for the genes within these modules, a number of targeted compounds were identified as being currently used in the clinic for treatment of breast cancer or other neoplasms; or in phase II and/or phase III development (https://clinicaltrials.gov/) (Fig. 3, Supplementary Table 6). Among these compounds, a number have potential for stratified use in the early luminal breast cancer setting (Table 4) for those deemed high-risk by our classifier. Therefore, these compounds hold potential for repurposing targeted therapies to early luminal breast cancers ( Supplementary  Fig. 7).

DISCUSSION
Relapse following endocrine treatment remains a significant clinical challenge, as more women die following treatment for ER+ disease than for any other breast cancer subtype. 3 Therefore, there is an ongoing need to identify women who are at risk for relapse following anti-endocrine therapy. More importantly, simultaneously identifying targets for future therapeutic intervention and the means to effectively stratify women to such targeted therapies will improve the clinical management of these patients, and potentially reducing their overtreatment, or conversely identifying patients who may be currently undertreated. Using 3825 patients from the TEAM pathology cohort, we are able to derive a signature that both significantly improves risk stratification and identifies genes for which there are drugs currently in use, or under evaluation (https://clinicaltrials.gov/) in other malignancies. These patients could potentially be matched to the specific functional modules within this 95 gene signature (Table 4, Supplementary Table 6). As alluded to by the prognostic capacity of the individual modules (Fig. 3), this approach has the potential to better stratify patients to existing targeted therapies based on the molecular drivers of their cancer, and/or to novel/ putative targets for in vitro validation studies ( Supplementary  Fig. 7). Despite the fact that our commercial risk score and subtype classification was derived based on NanoString RNA abundance  profiling, our work confirmed the recent findings of the UK-OPTIMA prelim trial 17 that most current breast cancer multiparametric risk tests provide broadly equivalent risk information for a population of women with ER+ breast cancers (Supplementary Fig. 6), but can exhibit discordance between tests at the individual patient level (Fig. 2, Supplementary Tables 3 and 4). However, while efforts are being made to identify those patients who may be sensitive to current standard cytotoxic chemotherapies, we must also look to the promise of targeted therapies against the molecular drivers of these high-risk patients as revealed by the genes in the signature. While current multiparametric tests can identify those who may benefit from current adjuvant chemotherapy regimens, none of these tests predicts response to a drug-specific chemotherapy. This challenge is hampered by the identification of driver pathways in addition to the complexities of both global and individual chemotherapeutic response. Using the information generated by this data we envision a model for future prospective clinical trial design through the examination and validation of the drugs targeting the gene modules comprising the 95-gene signature ( Supplementary Fig. 7, Supplementary Table 6). In this way, genes associated with the G2/M checkpoint, as identified in Module 1, such as BIRC5 (Survivin), could be targeted. Indeed YM155, a Survivin suppressor, was recently evaluated in the metastatic breast cancer setting in combination with docetaxel in a phase II, multicenter, open-label, 2-arm study. 20 However, the lack of up-front patient stratification for YM155 benefit likely contributed to the finding of no significant benefit in its addition to docetaxel, thus obscuring the potential benefit of targeting this pathway. While known to be overexpressed in breast cancers, the relatively higher expression of BIRC5 observed among our highrisk patients (q = 1.37× 10 −178 ) (Supplementary Table 5) suggests there is a tipping point of mRNA abundance leading to increased risk. We observed that all genes within Modules 1, 3, and 4 showed a higher expression among patients at higher risk for relapse which were statistically significant (Supplementary Table 5), reflecting the prominent role of cell cycle and proliferation in breast cancer pathogenesis. We see that Module 3 is characterized by pathways involving late mitotic events. The overexpression of CDK1 offers a theranostic target, with the use of Dinacilib or similar molecules, currently under evaluation in phase III trials (Supplementary Table 6). Regaining cell cycle and mitotic checkpoint control is another attractive mechanism for directed therapies, with theranostic targets such as PLK1 (Fig. 3), being treated with inhibitors in the preclinical and clinical setting. 21,22 The regulation of S-phase and DNA replication pathways of Module 4, including CCND1, supports the potential stratification of patients to Palbociclib or other CDK inhibitors (Supplementary Table 6). Findings for the PALOMA-1 trial 23 resulted in approval for Palbociclib (CDK4/6 inhibitor) in combination with Letrozole in the metatstatic breast cancer setting; paving the way for the randomization of high-risk patients with ER+/HER2-cancer and Fig. 3 Signaling modules within the 95-gene residual risk signature. a Summary of REACTOME interactions amongst the genes of the 95-gene residual risk signature. Six major interaction modules comprising 52 genes were identified from the 95-gene residual risk signature. Relationships between genes, between and within modules, are shown by connecting lines. Solid lines with arrows indicate known and direct positive relationships. Solid lines ending in a perpendicular line indicate a known negative regulatory relationship. Dotted lines indicate relationships linked by other genes. Genes with red circles indicate gene targets for which there are known targeted therapies or at phase II/III development based on the Integrity compound search tool (Thompson Reuters) and ClinicalTrials.gov (https://clinicaltrials.gov/). b Kaplan-Meier survival curves (left) for each module are shown, and representing the validation cohort. To the right of each Kaplan-Meier curve are risk score estimates grouped as quartiles with each group compared against Q1. Hazard ratios were estimated using Cox proportional hazards model and significance of survival difference was estimated using the log-rank test. All patients represented are those who only received endocrine treatment Molecular stratification of early breast cancer J Bayani et al residual disease, in the PENELOPE-B trial. While promising in the late and metastatic setting, CDK inhibitor use in the early breast cancer setting has not yet been adequately assessed, nor is there a validated method to stratify patients who would most benefit from this treatment. Interestingly, recent in vitro evidence of synergy between palbociclib with tamoxifen showed resensitization to tamoxifen in ER-resistant cell lines, 24 suggesting that the identification of those who may be ER-resistant, could experience greater benefit with the use combined use of endocrine therapy and a CDK inhibitor. Genes of Module 2 are characterized by receptor tyrosine kinase signaling, apoptosis and control of the cell cycle have druggable targets among the members of ERBBfamily of genes. Anti-HER therapies are effective in ERBB2/HER2positive patients, but crosstalk between other members of the Epidermal Growth Factor Receptor (EGFR)/ERBB family suggest the aberrant expression of members aside from ERBB2/HER2 could justify their use in the absence of HER2 amplification. EGFR inhibitors such as gefitinib and laptinib have shown efficacy in other malignancies, but only moderate success in breast cancer, suggesting that an improved method of patient selection is required to identify those who would benefit the most. Interestingly, 296/342 (86.5%) HER2-enriched-like patients were identified as high-risk by our classifier. However, with 33.8% of HER2enriched-like patients possessing confirmed HER2 gene amplification or protein over-expression, these results suggest some patients may benefit from therapy targeting the ERBB-family and associated pathways. In fact, the 95-gene signature was still prognostic irrespective of ERBB2/HER2-status, in this population of patients that pre-dates the use of anti-ERBB2/HER2 therapies ( Supplementary Fig. 5). Moreover, while ERBB3 and ERBB4 were found to be univariately prognostic and part of the final signature, ERBB2/HER2 expression was not (Supplementary Table 2). This data would suggest that in current clinical practice, a number of our ERBB2/HER2-positive, low-risk patients would have received anti-ERBB2/HER2 therapies, resulting with an outcome potentially no better than ERBB2/HER2-negative patients. With respect to high-risk patients who were also ERBB2/HER2-positive, anti-ERBB2/ HER2 treatment would have some benefit to a subset in this group, but it is clear that there are other molecular drivers of recurrence in this high-risk population. Downstream pathways of ERBB, like PIK3/AKT/mTOR, which was identified as a significant pathway by our analyses, supports the potential use of everolimus in patients identified as high-risk. 25,26 Interestingly, 278/352 (78.9%) of patients identified as Basal-like were classified as highrisk by our gene signature despite being clinically classified as ER+; highlighting the need to recognize the importance of molecular heterogeneity among the hormone receptor-positive cancers, and the implications for novel treatment.
We have demonstrated that a novel 95-gene signature of residual risk, which integrates nodal status, has significantly better clinical utility for early recurrence than the currently available multiparametric tests. The signature appears to remain prognostic for later recurrence, though confirmation of its prognostic ability should be evaluated upon longer follow-up. Unlike these tests, modular analysis of the genes in the signature, have identified several genes and pathways suitable for therapeutic intervention among the high-risk patients. While the clinical potential of novel therapeutic targets identified in this study are being investigated, 27 both additional independent clinical validation of the signature and preclinical validation of the effects of targeting these pathways are required prior to implementation of this approach in a clinical trial. It is clear, however, that there is a need for significant improvement in the targeted selection of patients suitable for new therapies, rather than the randomization of allcomers in future clinical trial design. As we have demonstrated, hormone-receptor positive cancers are molecularly heterogeneous, thus requiring novel treatment strategies (Fig. 2a and  Supplementary Table 3, Supplementary Table 6). Clearly a multiparametric gene signature, as we have shown here is one means of selection, but improved stratification must also include the integration of gene mutational and copy-number status. Therefore we propose that in order to improve the clinical management of women with early hormone-receptor positive breast cancer, future clinical trial design requires a multiparametric test that not only improves identification of high-risk patients, but also improves the selection of patients to existing therapeutics targeting key genes/pathways that underlies the signature.

MATERIALS AND METHODS
The TEAM trial was a multinational, open-label, phase III trial in which postmenopausal women with hormone receptor-positive 19 early breast cancer were randomly assigned to receive exemestane (25 mg) once daily, or tamoxifen (20 mg) once daily for the first 2.5-3 years; followed by exemestane (25 mg) (totaling 5 years of treatment) (Supplementary Fig. 8). Hormone-receptor (ER and PgR) and HER2 status by immunohistochemistry were locally assessed for entry into the trial and then centrally confirmed, 28 and HER2 status was confirmed by immunohistochemistry and fluorescence in situ hybridization. 29 All assessment was performed according to American Society of Clinical Oncology (ASCO)/College of American Physicians (CAP)/ guidelines. [30][31][32] None of the patients received anti-HER2 therapy. This study complied with the Declaration of Helsinki, Institutional Ethics Committee Guidelines, and the International Conference on Harmonisation and Good Clinical Practice guidelines. All patients provided informed consent. DRFS was defined as time from randomization to distant relapse or death from breast cancer. 19 The TEAM trial included a pathology research study comprised of 4736 patients from five countries with an average clinical follow-up of 6.86 years. Power analysis was performed to confirm the study size had 88.6 and 100% power to detect a HR of at least 3.0 in the training and validation cohorts respectively (Supplementary Figure 8B). RNA was available and successfully assayed from 3825 samples. Patients from the UK cohort were assigned as the training cohort (n = 790); while the remaining patients from Germany, Belgium, Netherlands, and Greece comprised a fullyindependent validation cohort (n = 3035). All patients were assayed for mRNA abundance (Supplementary Fig. 8C). To identify a signature of residual risk following endocrine treatment only, the main analyses excluded those patients who received neo-adjuvant and adjuvant chemotherapy. However, analyses of patients who also received adjuvant chemotherapy using the signature trained in the absence of chemotherapy are discussed and included in the Supplementary Results.

RNA extraction and expression profiling
Five 4 µm formalin-fixed paraffin-embedded (FFPE) sections per case were deparaffinised, tumor areas were macro-dissected and RNA extracted using the Ambion® Recoverall ™ Total Nucleic Acid Isolation Kit-RNA extraction protocol (Life Technologies TM , Ontario, Canada). RNA aliquots were quantified using a Nanodrop-8000 spectrophometer (Delaware, USA). All 3825 RNAs extracted from the TEAM pathology cohort were successfully assayed. Probes for each gene were designed and synthesized at NanoString® Technologies (Seattle, Washington, USA); and 250 ng of RNA for each sample were hybridized, processed and analyzed using the NanoString® nCounter® Analysis System, according to NanoString® Technologies protocols. mRNA abundance analysis and survival modeling Raw mRNA abundance count data were pre-processed using the NanoStringNorm R package 33 (v1.1.19) using normalization factors derived from the geometric mean of the top expressing 75 genes. Samples with RNA content |z-score| >6 were flagged and removed as outliers. To assess the performance of the chosen normalization method in this cohort, a combination of 252 preprocessing methods were evaluated. Firstly, each preprocessing method was ranked based on their ability to maximize Euclidean distance of ERBB2 mRNA abundance between HER2-positive and HER2-negative samples. The process was repeated for one million random subsets of HER2-positive and HER2-negative samples for each of the preprocessing schemes. Fifteen replicates of an RNA pool extracted from selected anonymized FFPE breast tumor samples were profiled across multiple batches; and preprocessing methods were ranked based on the inter-replicate variation. A mixed effects linear model was fit and residual Molecular stratification of early breast cancer J Bayani et al estimates were used as an estimate of inter-batch variation (nlme v3.1-117). Cumulative ranks based on these two criteria were calculated using RankProduct and the method chosen was amongst the top 10 preprocessing methods in rank product (Supplementary Data, Supplementary Fig. 9). Univariate survival analysis of preprocessed mRNA abundance data was performed by median-dichotomizing patients into high-expressionand low-expression groups. Clinical variable age was modeled as binary (dichotomized around age 55), while grade and nodal status were modeled as ordinal variables, and pathological size was modeled as a continuous variable.
Network-based signature derivation Feature-selection of genes was first performed based on univariate Cox proportional hazards modeling in the endocrine-treated only training cohort; those with p < 0.25 were retained. These retained genes were used to calculate a "module-dysregulation score" (MDS; Supplementary Data). A multivariate Cox proportional hazards model was then fit on MDSs, along with clinical covariates (age, grade, pathological size and nodal status); a stepwise backward selection approach using Akaike Information Criterion was performed to refine the multivariate model. The final selected model was trained in the training cohort and validated in the fully independent validation cohort (Table 1). DRFS truncated to 10 years was used as an endpoint. Recurrence probabilities were estimated as described in Supplementary Data. All survival modeling was performed on DRFS, in the R statistical environment with the survival package (v2.37-4). Model performances were evaluated through area under the receiver operating characteristic (ROC) curve (AUC, Supplementary Data).

Pathway analyses using reactome
The final gene list was loaded into the Cytoscape Reactome FI plugin in Cytoscape (v3.0.2). Symbols were loaded as a gene set with the 2013 version of the FI network. A FI network was constructed with FI annotations and no linker genes. Spectral clustering and pathway enrichment were computed for each module using the Reactome FI plugin functions.