The identification of biomarkers predictive of neoadjuvant chemotherapy response in breast cancer patients would be an important advancement in personalized cancer therapy. In this study, we hypothesized that due to similarities between radiation- and chemotherapy-induced cellular response mechanisms, radiation-responsive genes may be useful in predicting response to neoadjuvant chemotherapy. Murine p53 null breast cancer cell lines representative of the luminal, basal-like and claudin-low human breast cancer subtypes were irradiated to identify radiation-responsive genes across subtypes. These murine tumor radiation-induced genes were then converted to their human orthologs, and subsequently tested as a predictor of pathologic complete response (pCR), which was validated on two independent published neoadjuvant chemotherapy datasets of genomic data with chemotherapy response. A radiation-induced gene signature consisting of 30 genes was identified on a training set of 337 human primary breast cancer tumor samples that was prognostic for survival. Mean expression of this signature was calculated for individual samples on two independent published datasets and was found to be significantly predictive of pCR. Multivariate logistic regression analysis in both independent datasets showed that this 30 gene signature added significant predictive information independent of that provided by standard clinical predictors and other gene expression-based predictors of pCR. This study provides new information for radiation-induced biology, as well as information regarding response to neoadjuvant chemotherapy and a possible means of improving the prediction of pCR.
INTRODUCTION
Neoadjuvant chemotherapy has been widely used in recent years as part of the standard of care for locally advanced breast cancer patients. Although neoadjuvant and adjuvant chemotherapy have similar efficacy in terms of disease-free and overall survival rates (1), neoadjuvant chemotherapy has been shown to improve breast-conserving operability in locally advanced breast cancers (2, 3). Another advantage of neoadjuvant chemotherapy is that it allows for the direct and timely observation of tumor treatment response. Twenty to thirty percent of breast cancer patients who receive neoadjuvant chemotherapy achieve a pathologic complete response (pCR) (4), which has been correlated with improved long-term, disease-free and overall survival (1–3, 5–9) and thus is a valuable surrogate end point for survival.
The prediction of neoadjuvant chemotherapy response (i.e., pCR) is an active area of research. Since patients who achieve pCR after neoadjuvant chemotherapy are more likely to experience excellent cancer-free, long-term survival (1–3, 5–9), the accurate prediction of pCR would be of significant value. It would serve to identify those patients who could benefit most from neoadjuvant chemotherapy and identify those unlikely to benefit from therapy and thus be spared treatment-associated toxicities and be selected to receive alternative therapeutic approaches. Rouzier et al. (10) published in 2005, a nomogram consisting of clinical variables including tumor size, histologic grade and estrogen receptor (ER) status to predict pCR to neoadjuvant FAC (fluorouracil, doxorubicin and cyclophosphamide) or T/FAC (paclitaxel, FAC) chemotherapy. This same group reported in another study (11) that intrinsic breast cancer subtype as determined by gene expression profiling of tumors prior to treatment, was associated with pCR. In this study, basal-like and HER2-enriched (HER2-E) intrinsic subtypes were associated with the highest rates of pCR (45%) while luminal A and luminal B tumors showed the lowest rate of pCR (6%). Building on these studies, we sought to determine whether we could identify a novel biomarker that could add predictive information independent of that provided by known clinical variables and other gene expression-based predictors of pCR.
Cellular responses to radiation and chemotherapy involve multiple overlapping pathways. Both therapies involve cytotoxic agents that can ultimately induce apoptosis in target cells (12, 13). Radiation therapy and many chemotherapeutics result in DNA damage, activating multiple overlapping proteins/pathways, such as the ATM and ATR protein kinases, TP53 and TP53-independent pathways including CHK2-mediated signaling (14). Other proteins/pathways induced by both radiation and chemotherapy includes the extracellular signal-regulated kinase (ERK) pathway, Jun-N-terminal kinase (JNK) pathway, p38 MAPK family, NK-kB, AKT, mTOR and the checkpoint proteins (14–18). We hypothesized that due to the similarities in the cellular response mechanisms to radiation and chemotherapy, radiation-responsive genes may be useful in predicting response to neoadjuvant chemotherapy for breast cancer.
MATERIALS AND METHODS
Cell Culture
To identify radiation-responsive genes, we used cell lines derived from the T11, 2225L and 2250L murine mammary tumors described by Herschkowitz et al. (19). Briefly, these tumors were produced by removal and transplantation of 6-week-old BALB/c p53–/– mammary tissues into 3-week-old wild-type BALB/c recipients. This was necessary to circumvent the appearance of other tumor types that occurred with short latency in mice homozygous for p53 loss. As described by Herschkowitz et al. (19), based on murine-intrinsic gene list expression analysis, the T11 tumor had characteristics of the claudin-low human breast carcinoma subtype; the 2225L tumor had characteristics of the basal-like human breast carcinoma subtype; the 2250L tumor had characteristics of the luminal human breast carcinoma subtype. Respective tumor samples stored in liquid nitrogen were thawed and placed into appropriate media to allow cell line formation as follows: for the T11 cell line, the media used was RPMI 1640 containing 10% FBS, 100 U/ml penicillin and 100 μg/ml streptomycin; for the 2225L and 2250L cell lines, media used were HMEC media containing 5% FBS, 100 U/ml penicillin and 100 μg/ml streptomycin and HMEC supplement (Gibco® Life Technologies, Grand Island, NY). Cells were incubated in a humidified incubator at 37°C in 95% air/5% CO2. Genomic DNA was harvested from cell lines using QIAGEN DNeasy Kit (QIAGEN, Valencia, CA), and presence of the p53 null transgene was verified using polymerase chain reaction.
Cell Irradiation and Collection of RNA
Cells were plated in 150 mm dishes and grown until 50% confluence. Cells were then irradiated to a dose of 8 Gy using an RS 2000 irradiator (Rad Source Technologies Inc., Suwanee, GA) operating at a dose rate of 100 cGy/min. Cells were then immediately returned to the incubator and harvested at 4, 8, 12, 24 and 48 h after irradiation, at which point RNA was isolated using the QIAGEN RNeasy Mini Kit. A control nonirradiated RNA sample for each cell line was collected from cells harvested immediately after mock irradiation.
Microarray Experiments
Mouse whole-genome 4×180,000 features microarrays (Agilent Technologies Inc., Palo Alto, CA) were hybridized according to manufacturer's protocol with Cy3-CTP-labeled cRNA from mock-irradiated cells (2 ug/sample) and Cy5-CTP-labeled cRNA from irradiated cells (2 ug/sample), with replicates for a 24 h time point for each cell line. Microarrays were scanned and image files analyzed as described previously (19). All primary microarray data are available from the University of North Carolina (Chapel Hill, NC) Microarray Database ( https://genome.unc.edu) and the Gene Expression Omnibus [GEO; National Center for Biotechnology Information (NCBI)] ( http://www.ncbi.nlm.nih.gov/geo/) with series number GSE48073.
Analysis of Microarray Data to Identify Radiation-Responsive Genes
Data from microarray experiments were calculated as described (19), where we used the Lowess normalized log2R/G ratio. To identify radiation-responsive genes, we used a one-class Significance Analysis of Microarrays (SAM) to identify genes that changed in all time points for all three cell lines (as a single class) relative to the mock-irradiated cells (20). Using a false discovery rate (FDR) of 0%, SAM identified 4,278 radiation-induced genes and 2,689 radiation-repressed genes. Hierarchical cluster analysis was conducted using Cluster 3.0 (21) and results were visualized in Java TreeView (22).
To identify genes that were differentially expressed after irradiation across the cell lines, we used a multi-class SAM analysis to identify genes that changed in all time points after irradiation, but that were different between the cell lines (each cell line as its own class) relative to mock-irradiated cells. Using an FDR of 0%, SAM identified 1,289 genes differentially expressed after irradiation across the cell lines. Hierarchical cluster analysis was conducted using Cluster 3.0 (21) and results were visualized in Java TreeView (22).
Analysis of Human Primary Breast Tumor Microarray Data Using the Radiation-Induced Gene Set
The human primary breast tumor samples used as the training dataset are described in Prat et al. (23); this training dataset comprised a total of 337 tumor samples represented by microarray experiments from breast cancer patients (consisting of 320 breast tumor samples and 17 normal breast samples) heterogeneously treated in accordance with standard of care. To analyze this human training dataset with the radiation-induced murine gene set identified by the above described SAM, the list of 4,278 radiation-induced murine genes identified from the cell lines were converted to their human gene orthologs using the UCSC genome website ( http://genome.ucsc.edu/). A total of 1,964 unique human genes present on our human microarrays were identified. These genes were then used to hierarchically cluster the Prat et al. 337 human breast sample dataset using Cluster 3.0 to identify gene sets/clusters composed of genes that are highly coordinately expressed across these patient samples. We hypothesized that the mean expression of one or more of these gene sets/clusters may be useful in predicting an individual tumor's response to neoadjuvant chemotherapy. To avoid spurious results that may occur with small gene clusters or poorly correlated gene clusters, we limited the analysis to the most correlated gene clusters (node correlation >=0.45) and that were greater than 20 genes in size. This identified 23 non-overlapping gene clusters, where we next calculated a mean gene expression value of each gene cluster, for each of the 337 samples from the Prat et al. dataset (23).
By matching Entrez gene identifiers ( www.ncbi.nlm.nih.gov/gene), microarray data for as many genes as possible for each of the 23 gene clusters identified above was obtained on two independent test microarray datasets (24, 25). Briefly, the Hatzis et al. (24) dataset consisted of microarray data from 473 HER2-negative breast cancer tumors taken prior to sequential taxane and anthracycline-based neoadjuvant chemotherapy. This dataset had associated pCR and distant relapse free survival (DRFS) data. The Gluck et al. (25) dataset consisted of microarray data from 95 HER2-negative breast cancer tumors taken prior to neoadjuvant capecitabine and docetaxel chemotherapy; this dataset had associated pCR data and TP53 mutation status data using the AmpliChip TP53 assay. As done for the Prat et al. dataset, a mean gene expression value of each gene cluster was calculated for each patient in the Hatzis et al. dataset, and we validated the most significant gene cluster on the Gluck et al. dataset.
Statistical and Survival Analysis
Kaplan-Meier survival plots were compared using the Cox-Mantel log-rank test in WinSTAT® (R. Fitch Software, Staufen, Germany) for Excel (Microsoft Corp., Redmond, WA). Univariable and multivariable logistic regression analysis were used to assess the significance of gene cluster mean values to pCR. All statistical tests were two-tailed and P < 0.05 were declared significant. Analyses were performed using JMP Pro 9.0 (SAS Institute, Cary, NC) and R version 2.13.1 (R Development Core Team, Vienna, Austria).
RESULTS
Identification of Radiation-Responsive Genes
To identify radiation-responsive genes, we used cell lines derived from three different mouse p53-deficient tumors, including the T11 (claudin-low), 2225L (basal-like) and 2250L (luminal) tumors described by Herschkowitz et al. (19). A one-class SAM with a false discovery rate (FDR) of 0% identified 4,278 radiation-induced and 2,689 radiation-repressed genes in microarray experiments on the combined T11, 2225L and 2250L cell lines at 4, 8, 12, 24 and 48 h after a single radiation dose of 8 Gy. Hierarchical clustering of the radiation-induced genes is shown in Fig. 1 and that of the radiation-repressed genes is shown in Fig. 2.
Many genes identified as radiation induced in our analysis were previously known to be radiation responsive in humans including CREM, BNIP3, FAS, TNFRSF11B, IFITM1, LGALS3PB, COX7B and SESN1 (26–29), indicating the conservation of the radiation responsiveness of many genes in this p53-deficient mouse model. Using the program DAVID [Database for Annotation, Visualization and Integrated Discovery (30, 31); http://david.abcc.ncifcrf.gov], the gene ontology categories that were significantly over-represented relative to chance in the set of all radiation-induced genes included: “protein catabolic process” (P = 3.40 × 10−9); “apoptosis” (P = 3.45 × 10−5); “cell death” (P = 1.17 × 10−4) and “macromolecule catabolic process” (P = 8.24 × 10−8). Kyoto Encyclopedia of Genes and Genomes [KEGG (32); http://www.genome.jp/kegg/] pathways (which are graphical diagrams representing knowledge on molecular interaction and reaction networks for a wide variety of cellular processes) that were significantly over-represented relative to chance included: “pyruvate metabolism” (P = 3.65 × 10−5); “ubiquitin mediated proteolysis” (P = 0.003); “VEGF signaling pathway” (P = 0.003), “mTOR signaling pathway” (P = 0.004); and interestingly, “p53 signaling pathway” (P = 0.02). Since the cell lines used in our experiments were p53 null, the enrichment for members of the KEGG pathway “p53 signaling pathway” in the radiation-induced gene set suggested that these members known to be activated by p53 [such as Scotin (33) and the sestrins SESN1 and SESN2 (34, 35)] can also be activated independent of p53; Scotin has been confirmed independently by Terrinoni et al. to be able to be activated independent of p53 (36).
Subclusters of radiation-induced genes, as opposed to the complete list, showed enrichment of many distinct ontologies in the various gene clusters. For example, cluster A in Fig. 1, which consisted of genes whose expression was induced at a relatively constant level with respect to time across all cell lines, was significantly enriched for members of the gene ontology categories “mitochondrion” (P = 6.24 × 10−9) and “cellular protein catabolic process” (P = 0.02). In addition, cluster C in Fig. 1, which consisted of genes whose expression increased with time across all cell lines, was enriched for members of gene ontology categories “positive regulation of programmed cell death” (P = 0.002), “antigen processing and presentation” (P = 6.40 × 10−7), “lysosome” (P = 3.71 × 10−5) and “immune response” (P = 2.44 × 10−4). Other clusters of radiation induced genes and the respective gene ontologies enriched in those clusters are shown in Fig. 1.
Gene ontology categories significantly over-represented relative to chance in the set of all radiation-repressed genes included “spliceosome” (P = 3.19 × 10−14), “cell cycle” (P = 2.03 × 10−11), “cell division” (P = 3.33 × 10−7) and “mitosis” (P = 3.98 × 10−5). KEGG pathways significantly over-represented relative to chance in the set of radiation-repressed genes included “ribosome” (P = 5.45 × 10−22), “cell cycle” (P = 5.24 × 10−4) and “spliceosome” (P = 3.78 × 10−19). The enrichment of genes involved in apoptosis in the radiation-induced gene set, and enrichment of genes involved in mitosis in the radiation-repressed gene set, is consistent with prior studies examining radiation-responsive genes in humans (26–29) indicating conservation of the radiation responsiveness of these gene ontology categories/pathways in our mouse model. As for the radiation-induced genes, clusters of radiation-repressed genes showing highly correlated expression were observed (Fig. 2) and gene ontology analysis of these gene clusters showed enrichment. For example, cluster A in Fig. 2, which consisted of genes whose expression was more repressed in the 2250L cell line than in the other cell lines, was enriched for the gene ontology category “Zinc ion binding” (P = 0.01), while cluster G in Fig. 2 that consisted of genes whose repression increased with time across all cell lines, was enriched for the gene ontology categories “cell cycle” (P = 2.93 × 10−5), “RNA processing” (P = 1.20 × 10−21), “cell division” (P = 0.003), “mitosis” (P = 0.02) and “nuclear division” (P = 0.02); the repression of these gene clusters likely reflect reduced cell proliferation, which is a known cellular response to radiation. Other clusters of radiation-repressed genes and the respective gene ontologies enriched in those clusters are shown in Fig. 2.
Differential Response to Radiation among the Cell Lines
To identify genes that were differentially expressed after irradiation between the T11, 2225L and 2250L cell lines, we used a multi-class SAM analysis with each cell line as its own class. Using an FDR of 0%, SAM identified 1,289 genes differentially expressed after irradiation between the cell lines (hierarchical clustering of these genes is shown in Fig. 3), indicating that while the cell lines do share common genes that respond similarly to radiation (as discussed above and shown in Figs. 1 and 2), there are also clear differences in the cell lines' radiation response. These data suggest that the different intrinsic subtypes of breast cancer may have differences in their response to radiation.
Figure 3 shows that the genes found by SAM to be differentially expressed after irradiation between the cell lines fall into distinct gene clusters that differ in terms of how they respond across the three cell lines. For example, cluster A in Fig. 3 (Supplementary Fig. S1; http://dx.doi.org/10.1667/RR13485.1.S1 (10.1667 RR13485.1.S1.pdf)) consists of genes that are repressed in the 2225L basal-like cell line, but that are induced in the other 2 cell lines in response to radiation. Gene ontology categories significantly enriched in this cluster according to DAVID include “response to wounding” (P = 0.005) and “inflammatory response” (P = 0.02). This gene cluster was also significantly enriched for members of the KEGG pathway “MAPK signaling” (P = 0.01). Cluster B in Fig. 3 (Supplementary Fig. 2; http://dx.doi.org/10.1667/RR13485.1.S1 (10.1667 RR13485.1.S1.pdf)) consists of genes observed to be repressed in the T11 claudin-low and 2225L basal-like cell lines, while less repressed in the 2250L luminal cell line in response to radiation. Gene ontology categories significantly enriched in this cluster include “M phase” (P = 0.006), “mitotic cell cycle” (P = 0.01) and “ribosome” (P = 1.25 × 10−4). Cluster C in Fig. 3 (Supplementary Fig. 3; http://dx.doi.org/10.1667/RR13485.1.S1 (10.1667 RR13485.1.S1.pdf)) consists of genes repressed in the 2250L and 2225L cell lines while less repressed or slightly induced in the T11 cell line in response to radiation; gene ontology categories significantly enriched in this cluster include “myofibril assembly” (P = 0.002), “muscle contraction” (P = 0.02) and “muscle cell development” (P = 0.02). Cluster E in Fig. 3 (Supplementary Fig. 5; http://dx.doi.org/10.1667/RR13485.1.S1 (10.1667 RR13485.1.S1.pdf)) consists of genes induced in the 2250L and 2225L cell lines while repressed or less induced in the T11 cell line in response to radiation; gene ontology categories significantly enriched in this cluster include “zinc ion binding” (P = 5.76 × 10−4) and “DNA binding” (P = 0.006).
Cluster D in Fig. 3 (Supplementary Fig. 4; http://dx.doi.org/10.1667/RR13485.1.S1 (10.1667 RR13485.1.S1.pdf)) consists of genes induced in the 2225L basal-like cell line while repressed or less induced in the other 2 cell lines in response to radiation. Gene ontology categories significantly enriched in this cluster include “blood vessel development” (P = 0.003), “positive regulation of mesenchymal cell proliferation” (P = 0.008) and “cell motion” (P = 0.02). Interestingly this gene cluster contained the genes SESN3 and TCF4 found to be regulated by the transcription factor YBX1 by Evdokimova et al. (37). In the experiments done by Evdokimova et al. the enforced expression of YBX1 in noninvasive breast epithelial cells was found to directly activate translation of mRNAs encoding Snail1 and other transcription factors such as SESN3 and TCF4 that are implicated in activation of mesenchymal genes, resulting in induction of an epithelial-mesenchymal transition (EMT) accompanied by enhanced metastatic potential. Other genes in cluster D in Fig. 3 such as DDR1, FLT1, NOTCH3 and PDZRN3 have been identified as candidate YBX1-regulated genes through chromatin immunoprecipitation (ChiP)-on-chip analysis performed by Finkbeiner et al. (38). Also included in cluster D is the gene STAT3, the prosurvival pathway of which has been shown to be engaged by YBX1 to protect cells from apoptosis (39). Another study has suggested that STAT3 may indirectly activate transcription of YBX1 by activating transcription of TWIST, which in turn activates transcription of YBX1 (40). In summary, our multiclass SAM analysis has identified genes differentially expressed after irradiation between the cell lines, suggesting that differences in radiation response exist within the various breast cancer intrinsic subtypes and further, suggests possible biological bases for these differences in radiation response.
Analysis of Human Breast Tumors Using the Radiation-Induced Gene Set
We hypothesized that due to similarities between the cellular response mechanisms to radiation and chemotherapy, expression differences of radiation-induced genes may be useful in predicting response to neoadjuvant chemotherapy. To test this hypothesis, first we hierarchically clustered a training set of 337 human primary breast samples (consisting of 320 breast tumor samples and 17 normal breast samples) published by Prat et al. (23), using the murine radiation-induced gene set of 4,278 genes (shown in Fig. 1) converted into their human orthologs (complete cluster diagram shown in Fig. 4A); this was done to identify gene sets/clusters composed of genes that are highly concordantly expressed within human breast samples.
We hypothesized that the mean expression of one or more of these gene sets/clusters may be useful in predicting an individual tumor's response to neoadjuvant chemotherapy. To avoid spurious results that may occur with small gene clusters, or poorly correlated gene clusters, we limited analysis to the most correlated gene clusters (node correlation >=0.45) and those nodes that were greater than 20 genes in size. This identified 23 non-overlapping gene clusters. We then tested the ability of the mean expression of each of these 23 gene sets to predict response to neoadjuvant chemotherapy by using a dataset published by Hatzis et al. (24), which consisted of microarray data from 473 HER2-negative breast cancer tumors taken prior to sequential taxane and anthracycline-based neoadjuvant chemotherapy; this dataset had associated pCR and distant relapse-free survival (DRFS) data. Using univariate logistic regression modeling, only 4/23 of the tested gene clusters had a significant Bonferroni-corrected P value in predicting pCR. To determine which of these 4 gene clusters could be useful in predicting pCR, we then performed multivariate logistic regression modeling with each of these 4 gene sets in a model that included the standard clinical variables and other gene expression based predictors including those used by Hatzis et al. (24, 41). Only 1 gene cluster (consisting of 30 genes, Fig. 4B and Table 1) was significant in multivariate analysis (Table 3). As seen in Tables 2 and 3, a higher mean expression of this 30 gene set was associated with a higher likelihood of pCR. As shown in Table 2B, this 30 gene set was not only able to significantly predict pCR in univariate logistic regression modeling on the entire Hatzis et al. patient dataset, but was also able to significantly predict pCR within individual clinical subsets of patients, including the clinically relevant triple-negative subset and the more biologically relevant basal-like patient subset.
We next tested the ability of the 30 gene cluster to predict pCR on a second and completely independent test dataset published by Gluck et al. (25), which consisted of microarray data from 95 HER2-negative breast cancer tumors taken prior to neoadjuvant capecitabine and docetaxel chemotherapy. This dataset had associated pCR data and TP53 mutation status data using the AmpliChip TP53 assay. As shown in Table 4, the 30 gene set significantly predicted pCR in univariate logistic regression modeling on the entire dataset as well as on the biologically relevant basal-like patient subset. In multivariable logistic regression modeling (Table 5) that included the standard clinical parameters, intrinsic subtype and TP53 mutation status, the 30 gene set again remained a statistically significant predictor of pCR.
As both the Prat et al. and Hatzis et al. datasets had associated survival data, we also tested the prognostic value of the 30 gene set to predict survival. To achieve this goal, we grouped patients from each dataset into halves or tertiles based on rank order mean expression value of the 30 gene set. Kaplan-Meier analysis showed significant differences among patients when divided into halves or tertiles for distant relapse-free survival (DRFS) in the Hatzis et al. dataset (P = 2.48 × 10−7 and 0.0002, respectively), and for overall survival (OS; P = 1.45 × 10−6 and 1.18 × 10−7, respectively) and relapse-free survival (RFS; P = 2.08 × 10−5 and 7.95 × 10−5, respectively) in the Prat et al. dataset. In the Hatzis et al. dataset, patients with tumors of the basal-like subtype, when grouped into halves or tertiles based on rank order mean expression value of the 30 gene set, showed significant differences in DRFS (Fig. 5), with patients having higher expression of the 30 gene set having better DRFS.
Gene ontology analysis of the 30 gene set using DAVID revealed significant enrichment for the following categories: “pyruvate metabolic process” (genes LDHB, PPARGC1A, PDHA1; P = 0.002); “response to extracellular stimulus” (genes ADORA2B, PPARGC1A, PLA2G4A, VLDLR; P = 0.007); and “generation of precursor metabolites and energy” (genes LDHB, PPARGC1A, PDHA1, SLC25A27; P = 0.017). Interestingly, our 30 gene set was also significantly enriched for members of the KEGG pathway “vascular smooth muscle contraction” (genes PLA2G4A, PRKX, ADORA2B; P = 0.027). When it was compared with other gene expression-based predictors, only one common gene (NFIB) was found between our 30 gene predictor and the predictors published by Hatzis et al. (24). There was no overlap with either the 11-gene proliferation signature developed by Parker et al. (42), the PAM50 gene signature used by Gluck et al. or Prat et al. (23, 25, 43) or the 21 gene recurrence score assay developed by Paik et al. (44), which is known as Oncotype DX.
DISCUSSION
The search for biomarkers that predict pCR to neoadjuvant chemotherapy for breast cancer has been an area of intense research in recent years. pCR is an important clinical end point, as patients who achieve pCR to neoadjuvant chemotherapy are more likely to have improved overall and disease-free survival (3). Gene expression-based analyses have contributed to this research, with the development of several predictive gene sets and assays based on selection of genes that directly correlate with patient/tumor outcomes (24, 45, 46). We took a different approach and selected genes using no knowledge of survival outcomes or response to therapy; these were instead selected on the basis of radiation induction and their natural patterns of expression in primary breast tumors. We objectively identified gene sets/clusters composed of radiation-induced genes that were highly concordantly expressed in human breast tumor samples. We then systematically tested the ability of each of these gene sets to predict pCR on an independent test dataset (24). Only one gene cluster consisting of 30 genes was able to significantly predict pCR in both univariate analysis and multivariate analysis that included standard clinical variables and breast intrinsic subtype. In addition, this 30 gene signature was able to predict pCR in both univariate and multivariate analysis on a second independent test dataset (25), and thus it has shown utility on two different validation data sets.
As determined by gene ontology analysis, our radiation-induced 30 gene signature was significantly enriched for genes involved in pyruvate metabolism, generation of precursor metabolites and energy, suggesting that tumors most sensitive to chemotherapy are those with highest metabolic activity, consistent with the known mechanisms of action of several chemotherapeutics used in the datasets we analyzed, including doxorubicin and capecitabine (47–53). Genes composing our 30 gene signature include lactate dehydrogenase B (LDHB), which is a subunit of the lactate dehydrogenase enzyme, the key glycolytic enzyme catalyzing formation of lactic acid from pyruvate. LDHB has recently been shown to be a downstream target of mTOR critical for oncogenic mTOR-mediated tumorigenesis (54). It was also recently shown to be an essential gene in basal-like/triple-negative breast cancer metabolism (55). Proteomic analysis done by Cortesi et al. (56) showed that high expression of the LDHB protein in tumor interstitial fluid was associated with response to chemotherapy in breast cancer patients, consistent with our results. Interestingly, the downregulation of DAB2 interacting protein (DAB2IP), another member of our 30 gene set and a novel member of the Ras GTPase-activating protein family, was recently reported by Kong et al. (57) to result in resistance to ionizing radiation in prostate cancer cells. Our results suggest that the low expression of DAB2IP is also associated with resistance to chemotherapy in breast cancer.
Not only has our analysis of radiation-responsive genes yielded a useful tool for predicting pCR, it has also yielded information indicating that the different intrinsic subtypes of breast cancer may have distinct biological differences in their response to radiation. As shown in Fig. 3, the genes found to be differentially expressed after irradiation between the murine cell lines representative of the different intrinsic subtypes (claudin-low, basal-like and luminal) fall into distinct clusters. These clusters differ not only in how they respond in the cell lines to radiation but also their biological function, possibly reflecting differential activation or repression of various pathways in the different subtypes in response to radiation. To support this hypothesis, genome wide sequencing of the 2225L basal-like tumor has shown a possible YBX1 activating mutation in this tumor, but not in the T11 or 2250L tumors (unpublished results). Consistent with the possibility of an YBX1 activating mutation in the 2225L tumor, we observed induction of multiple YBX1 known and putative targets in the 2225L basal-like cell line (Fig. 3, gene cluster D; Supplementary Fig. S4: http://dx.doi.org/10.1667/RR13485.1.S1 (10.1667 RR13485.1.S1.pdf)) in response to radiation, which was not seen in the other cell lines. Active/expressed YBX1 appears to be an important feature of the basal-like subtype (58–61). With our observations, this suggests an important role of this gene not only in defining the basal-like subtype but also in the radiation response of this subtype. Experiments are ongoing in our lab to test this hypothesis.
In conclusion, our radiation-induced 30 gene signature shows significant promise as a clinically relevant predictor of pCR in breast cancer patients. We have shown its efficacy in predicting pCR in two independent test datasets. Furthermore we have shown by multivariate analysis that it adds significant information in predicting pCR beyond what is provided by the standard clinical parameters, breast tumor intrinsic subtypes, p53 mutation status and other gene expression-based predictors including those published by Hatzis et al. (24). Our results are hypothesis generating, but based on the analysis presented here, we feel that our predictor warrants further investigation, including validation in prospectively designed clinical trials, to confirm the clinical validity and utility of our 30 gene signature.
SUPPLEMENTARY INFORMATION
Supplementary Fig. S1. http://dx.doi.org/10.1667/RR13485.1.S1 (10.1667 RR13485.1.S1.pdf); Gene cluster A shown Fig. 3 (upper left) significantly enriched for genes involved in wound and inflammatory response as determined by DAVID.
Supplementary Fig. S2. http://dx.doi.org/10.1667/RR13485.1.S1 (10.1667 RR13485.1.S1.pdf); Gene cluster B from Fig. 3 (upper left) significantly enriched for genes involved in M phase, mitotic cell cycle and ribosome as determined by DAVID.
Supplementary Fig. S3. http://dx.doi.org/10.1667/RR13485.1.S1 (10.1667 RR13485.1.S1.pdf); Gene cluster C from Fig. 3 (upper left) significantly enriched for genes involved in muscle contraction, muscle cell development, and myofibril assembly as determined by DAVID.
Supplementary Fig. S4. http://dx.doi.org/10.1667/RR13485.1.S1 (10.1667 RR13485.1.S1.pdf); Gene cluster D from Fig. 3 (upper left) significantly enriched for genes involved in positive regulation of mesenchymal cell proliferation, blood vessel development and cell motion as determined by DAVID.
Supplementary Fig. S5. http://dx.doi.org/10.1667/RR13485.1.S1 (10.1667 RR13485.1.S1.pdf); Gene cluster E from Fig. 3 (upper left) significantly enriched for genes involved in zinc ion and DNA binding as determined by DAVID.
ACKNOWLEDGMENTS
This work was supported by funds from NCI Breast SPORE (P50-CA58223-09A1), R01-CA138255, R01-CA148761 and the Breast Cancer Research Foundation.
REFERENCES
FIG. 1.
Hierarchical clustering analysis of genes determined by one-class significance analysis of microarrays (SAM) to be radiation-induced (4,278 genes) when testing three different mouse mammary tumor cell lines. Colored bars and/or letters indicate various gene clusters as discussed in this article. To the left of the diagram are the gene ontology categories significantly enriched in the corresponding gene clusters indicated by the colored bars and/or letters.
FIG. 2.
Hierarchical clustering analysis of genes determined by one-class significance analysis of microarrays (SAM) to be radiation-repressed (2,689 genes) when testing three different mouse mammary tumor cell lines. Colored bars and/or letters indicate various gene clusters as discussed in the text. To the left of the diagram are the gene ontology categories significantly enriched in the corresponding gene clusters indicated by the colored bars and/or letters.
FIG. 3.
Hierarchical clustering analysis of genes determined by multi-class significance analysis of microarrays (SAM) to be differentially expressed across the three murine-derived mammary tumor cell lines T11, 2225L and 2250L in response to radiation (1,289 genes). Colored bars and/or letters indicate various gene clusters as discussed in this article. To the left of the diagram are the gene ontology categories significantly enriched in the corresponding gene clusters indicated by the colored bars and/or letters.
FIG. 4.
Hierarchical cluster analysis of the Prat et al. (49) 337 primary human breast sample dataset (consisting of 320 breast tumor samples and 17 normal breast samples) using the radiation-induced murine gene list converted into human orthologous genes. Panel A: Scaled down representation of the complete cluster diagram. Panel B: 30 gene set found to be significant in predicting pCR and survival.
FIG. 5.
Kaplan-Meier survival curves of patients with basal-like tumors from the Hatzis et al. dataset (24) divided into (panel A) halves or (panel B) tertiles according to rank order mean expression value of the radiation-induced 30 gene set. P values were calculated using the log-rank test.
TABLE 1
Radiation-Induced 30 Gene Signature List, Which Was Predictive of Pathologic Complete Response (pCR) to Neoadjuvant Chemotherapy in Breast Cancer Patients
TABLE 2
Univariate Logistic Regression Analysis to Predict Pathologic Complete Response (pCR) Using the Radiation-Induced 30 Gene Signature for (Part A) All Patients and (Part B) Select Patient Subsets of the Hatzis et al. (24) Dataset
TABLE 3
Multivariate Logistic Regression Analysis to Predict Pathologic Complete Response for the Hatzis et al. (24) Dataset Using Multiple Possible Predictive Factors Including the Radiation-Induced 30 Gene Signature, the Standard Clinical Parameters and Additional Genomic Biomarkers
TABLE 4
Univariate Logistic Regression Analysis to Predict Pathologic Complete Response (pCR) Using the Radiation-Induced 30 Gene Signature for (Part A) All Patients, and (Part B) Basal-like Patient Subset of the Gluck et al. (22) Dataset
TABLE 5
Multivariate Logistic Regression Analysis to Predict Pathologic Complete Response for the Gluck et al. (22) Dataset Using Multiple Possible Predictive Factors Including the Radiation-Induced 30 Gene Signature, the Standard Clinical Parameters and Other Possible Biomarkers
