Molecular signature of response to preoperative radiotherapy in locally advanced breast cancer

Background Radiation therapy is an indispensable part of various treatment modalities for breast cancer. Specifically, for non-inflammatory locally advanced breast cancer (LABC) patients, preoperative radiotherapy (pRT) is currently indicated as a second line therapy in the event of lack of response to neoadjuvant chemotherapy. Still approximately one third of patients fails to respond favourably to pRT. The aim of this study was to explore molecular mechanisms underlying differential response to radiotherapy (RT) to identify predictive biomarkers and potential targets for increasing radiosensitivity. Methods The study was based on a cohort of 134 LABC patients, treated at the Institute of Oncology and Radiology of Serbia (IORS) with pRT, without previous or concomitant systemic therapy. Baseline transcriptional profiles were established using Agilent 60 K microarray platform in a subset of 23 formalin-fixed paraffin-embedded (FFPE) LABC tumour samples of which 11 radiotherapy naïve and 3 post-radiotherapy samples passed quality control and were used for downstream analysis. Biological networks and signalling pathways underlying differential response to RT were identified using Ingenuity Pathways Analysis software. Predictive value of candidate genes in the preoperative setting was further validated by qRT-PCR in an independent subset of 60 LABC samples of which 42 had sufficient quality for data analysis, and in postoperative setting using microarray data from 344 node-negative breast cancer patients (Erasmus cohort, GSE2034 and GSE5327) treated either with surgery only (20%) or surgery with RT (80%). Results We identified 192 significantly differentially expressed genes (FDR < 0.10) between pRT-responsive and non-responsive tumours, related to regulation of cellular development, growth and proliferation, cell cycle control of chromosomal replication, glucose metabolism and NAD biosynthesis II route. APOA1, MAP3K4, and MMP14 genes were differentially expressed (FDR < 0.20) between pRT responders and non-responders in preoperative setting, while MAP3K4 was further validated as RT-specific predictive biomarker of distant metastasis free survival (HR = 2.54, [95%CI:1.42–4.55], p = 0.002) in the postoperative setting. Conclusions This study pinpoints MAP3K4 as a putative biomarker of response to RT in both preoperative and postoperative settings and a potential target for radiosensitising combination therapy, warranting further pre-clinical studies and prospective clinical validation. Electronic supplementary material The online version of this article (10.1186/s13014-018-1129-4) contains supplementary material, which is available to authorized users.


Background
Non-inflammatory locally advanced breast carcinoma (LABC) is a late stage breast cancer presented as a bulky primary chest wall tumour and/or extensive adenopathy including patients with large (> 5 cm), usually inoperable tumours and node positive disease [1]. It is a common presentation worldwide but is of special concern in developing countries with limited breast cancer awareness and efficient population screening programs. For instance, in Serbia over 4500 new breast cancer cases are diagnosed each year, with as many as 30% presenting as late stage initially inoperable LABC.
A multimodal approach including systemic therapy, radiotherapy and surgery is usually applied in the treatment of LABC [2]. Currently, neoadjuvant systemic therapy (CHT) is usually administered to downstage the tumour for breast-conserving surgery, while preoperative radiotherapy (pRT) is often indicated if there is no objective reduction of tumour volume after the administration of neoadjuvant (CHT). Radiation therapy (RT) is frequently used in various modalities for treatment of breast cancer of different stages including LABC, and large meta-analyses of multiple randomised trials demonstrated clear long-term benefit both in terms of locoregional control and reduced mortality in breast cancer patients treated with radiotherapy after breast conserving surgery and after radical mastectomy [3][4][5]. Given the rarity of this treatment modality only a few studies looked into the effects of preoperative radiotherapy (pRT) in combination with breast conserving surgery on locoregional recurrence and overall survival reporting similar results compared to protocols involving neoadjuvant chemotherapy without irradiation [6][7][8][9]. However, the molecular basis of tumour sensitivity to radiotherapy is complex, and at present, there are no conclusive biomarkers in clinical use to predict if a patient will, in fact, benefit from radiotherapy.
The concept of personalised medicine has been successfully implemented in medical oncology for over a decade with several biomarkers approved for clinical use. The same principle could be applied to radiation oncology to achieve better clinical responses to radiotherapy, lower radio-toxicity and avoid overtreatment [10,11]. To achieve this goal there is a clear need to develop biomarkes specific for breast radiotherapy. Hovewer, most of the studies in radiation oncology have been limited to the study of biomarkers not necessarily chosen based on their specificity to radiotherapy. Molecular subtypes have shown limited predictive estimation of RT efficacy and have been potentially confounded by adjuvant systemic therapy [11][12][13]. Several studies aimed to identify a radiosensitivity molecular signature in breast cancer by studying changes in gene and protein expression in response to radiation in cellular and animal model systems. These included determination of cellular radiosensitivity defined by survival fraction at 2 Gy [14], clonogenic doubling time, hypoxic fraction, or clonogenic number [15], with some of these multigene signatures having been validated in retrospective studies in solid tumours including breast cancer [16][17][18].
Although these studies provided valuable insights into cellular radiation response, the breast tissue has a complex microenvironment composed of several interacting cell types and extracellular molecules that may affect tumour response. The ideal model system for researching the breast tumour response to radiation therapy and evaluating the predictive value of markers is the preoperative setting [19,20]. Detection and characterization of biomarkers in the patient's tumour biopsies with known clinical response, before and after radiation therapy could select the group of patients with worse response to radiotherapy, to facilitate the choice of more efficient treatment, avoid overtreatment and consequentially reduce associated healthcare costs.
Leveraging a cohort of LABC patients treated with pRT without neoadjuvant or concomitant CHT (IORS LABC cohort) we determined baseline molecular differences between radio-resistant and radio-responsive breast tumours and identified putative predictive biomarkers of response to pRT. In our previous study on this cohort [21] we have shown that the extent of the clinical response to pRT in LABC is predictive of overall survival. Here, we analysed global gene expression patterns and biological pathways associated with differential response to pRT in a discovery subset of the IORS LABC cohort. Selected biomarkers were validated by an orthogonal assay (qPCR) in an independent validation set of pRT IORS LABC samples, and in the postoperative setting using an external microarray dataset (Erasmus cohort) consisting of breast cancer patients treated only with radiotherapy following surgery.

Patient cohorts IORS LABC cohort
This retrospective cohort included 134 patients with locally invasive non-inflammatory breast cancer (LABC) (93 patients stage III-a, and 41 patients III-b) treated with pRT between 1997 and 2000, delivering 45Gy in 15 fractions every second day alternately to the breast and regional lymph nodes, followed by radical mastectomy and adjuvant chemo and/or hormonal therapy (Fig. 1a). The median follow-up was 74 months (4-216) counted from the breast cancer diagnosis until last check-up or death from any cause. Five-year overall survival was 56% and 5-year disease free survival was 39.2%. Tumour biopsy was taken both prior to radiation treatment to obtain a radiotherapy-naïve sample (series A), and after RT and radical mastectomy before adjuvant treatment (series B) stored as formalin-fixed paraffin-embedded (FFPE) tissue samples. Clinical response to pRT in the breast was defined per RECIST criteria [22]. All patients gave their informed written consent for the use of residual tissue for research. The study was approved by the Institute of Oncology and Radiology Ethical Review Board for human studies. An overview of the patients' clinicopathological characteristics is summarised in Additional file 1: Table S1.

Erasmus breast cancer cohort
For biomarker validation, we have used a previously published clinical data set from the Erasmus cohort (GSE2034 and GSE5327). This cohort includes 344 lymph-node negative breast cancer patients treated at the Erasmus Medical Center (Rotterdam, Netherlands) from 1980 to 1995, who had not received neoadjuvant nor adjuvant systemic treatment [23]. Primary treatment was breast-conserving surgery or modified radical mastectomy, and 87% of the patients received postoperative RT. Early metastasis was defined as distant recurrence in the first 5 years following completion of primary treatment.

RNA extraction from FFPE tissue
Tumour tissue sections were stained with H&E and examined by a pathologist. Tumour tissue was macrodisected prior to RNA extraction to ensure there was > 70% tumour content. Total RNA was extracted from 3 to 5 10 μm thick FFPE tissue sections using RNeasy A B Fig. 1 a Infographic summarizing treatment protocol and sample collection of locally advanced breast cancer patients treated at the Institute of Oncology and Radiology of Serbia between 1997 and 2000 (IORS-LABC cohort). The IORS-LABC cohort included 134 patients who had initial biopsy taken before any treatment, followed by radiotherapy and radical mastectomy. Exceeding tumour material was formalin-fixed and paraffin embedded (FFPE) and stored at room temperature. b Flowchart representing the study outline for sample processing, quality control, data analysis and biomarker validation. FFPE samples were review by a pathologist to select those with > 70% of tumour material, retaining 118 pre-RT biopsy tumour samples (a) and 21 post-RT tumour samples (b). These samples were split into discovery (N A = 22 pre-RT and N B = 21 post-RT) and validation (N A = 96) subsets. After quality control of extracted RNA only 23 samples from the discovery subset were selected for microarray hybridization (N A = 18 and N B = 5), of which only 14 (N A = 11 and N B = 3) passed data quality control. Out of 96 pre-RT samples designated for validation, only 60 had passed RNA quality control and were used for qRT-PCR. Of those, 42 samples passed data quality control and were retained for the downstream analysis FFPE Kit (Qiagen) with an 18-h Proteinase K tissue digestion step. RNA quantity and purity were assessed by BioSpec-nano (Shimadzu Scientific Instruments).

Microarray hybridization and data analysis
Agilent SurePrint G3-Hmn-GE-v.2-8x60K Microarray platform was used for gene expression profiling following Agilent Gene Expression FFPE Workflow. Raw data pre-processing and quality control was performed using R version 3.0.1 and R/Bioconductor packages 'limma' , 'ffpe' and ' ArrayQualityMetrics'. Data were deposited in the GEO database under the GSE101920 accession number. Hierarchical average linkage clustering was performed using Cluster 3.0 and visualised using JavaTreeView [24]. Differential expression analysis was performed using the POMELO II tool applying moderated t-test [25]. The estimated significance levels were corrected for multiple hypotheses testing using Benjamini and Hochberg False Discovery Rate (FDR) adjustment [26]. The ranked target list of the differentially expressed genes was subjected to pathway enrichment analyses using Ingenuity Pathway Analysis software (Ingenuity Systems). Significantly enriched gene networks and canonical pathways based on the curated IPA database KEGG, Biocarta, and Reactome, were identified as previously described [15]. Methods are described in detail in Additional file 2: Supplementary Methods.

Quantitative RT-PCR and statistical analysis
Applied Biosystems High Capacity cDNA Reverse Transcription Kit was used for preparing cDNA from 200 ng RNA. Quantitative RT-PCR was done on an ABI Prism 7300 (Applied Biosystems) using TaqMan® Gene Expression Assays and TaqMan® PreAmp Master Mix Kit (Life Technologies). All qPCR amplicons were designed to be less than 100 bp long and all assays were done in triplicate. Assay ID numbers are shown in Supplementary  Table 4. Each plate included a HeLa cell line as inter-plate calibrator (IPC) and a non-template control (NTC). Average Ct values for each gene were standardised to IPC, dCt values were calculated relative to ACTB as a reference gene. Genes and samples with over 70% of missing values were excluded from further analysis, retaining 8 genes and 42 samples for further analysis. Undetermined values were set to the number of cycles performed (Ct = 45). Data was log2 transformed and differences in gene expression levels between groups were tested using Student t-test and corrected for multiple testing with FDR set to 20%.

Erasmus dataset processing and statistical analysis
Erasmus dataset was downloaded from GEO database and processed using affy R package. Distant metastasis free survival (DMFS) was defined as any distant recurrence within 5 years after the end of treatment. Survival curves were plotted with the Kaplan-Meier method and log-rank test was used to evaluate differences between groups defined by candidate gene expression status. Cutoff Finder was used to determine the optimal cutpoint for gene expression dichotomization based on the log-rank test minimum P-value approach [27]. Hazard ratios were estimated using the Cox proportional hazards model, stratified by RT status in a univariable and multivariable analysis. Pearson's χ 2 test was used to check for unbalanced distribution of clinico-pathologcal varables (ER, PR, T-stage and age categories) in subgroup analysis (Additional file 3: Table  S2). Measures of biological interaction were determined both on additive scale and multiplicative scale [28,29]. Stata command icp was used for calculating 3 different measures of interaction contrast on an additive scale: relative excess risk due to interaction [RERI], attributable proportion [AP] and synergy index [S] as described in [30]. Multiplicative interaction was assessed by including an interaction term with main effects in the Cox proportional hazard model. Statistical calculations were performed using STATA version 11.2 (StataCorp). All reported p-values were two-sided with a 0.05 significance level.

Molecular signature of radiosensitivity in LABC tumour samples
To gain a better understanding of the molecular response to radiotherapy independently of systemic treatment and to identify a baseline transcriptional signature of radiosensitivity between radio-responsive and non-responsive tumours, we have analysed radiotherapy-naïve LABC tumour biopsies and post radiotherapy mastectomy samples by gene expression profiling. A total of 43 FFPE tumour samples (N A = 22, N B = 21) from the patient cohort was randomly selected maintaining balanced group representation of clinical response (CR, PR, SD) and other clinicopathological characteristics (Additional file 4: Table S1). Of those, due to the low concentration or purity of the RNA extracted from FFPE only 23 samples were selected for microarray hybridization. Following stringent microarray quality control 14 samples (N A = 11, N B = 3) were included in subsequent data analysis (Additional file 4: Figure S1), comprising of 8 patients with stable disease designated as non-responders (NR) and 6 patients which experienced either complete (2 pts) or partial clinical response (4 pts) to pRT, classified as responders (R) (Fig. 1b).
Although unsupervised clustering analysis over top 20% most variable transcripts revealed separation between non-responders and responders (Additional file 5: Figure S2), not enough post-RT samples were left after QC to draw relevant conclusions regarding changes in pRT response between pre-RT and post-RT biopsies. Therefore, to study intrinsic differences underlying differential response to radiotherapy in the neoadjuvant setting, we have analysed only pre-RT tumour samples (N A = 11) to look for differences in baseline transcriptional profiles between responders (R) and non-responders-(NR). We identified 192 significantly differentially expressed mRNA transcripts (> 2-fold change and FDR < 0.1), including 89 annotated protein coding genes (PCG) and 78 long non-coding RNAs (lncRNAs) (Additional file 6: Table S3). Of those, only 7 genes were found to be upregulated while the rest of the genes (185) were downregulated in radio-responsive tumours (R) (Fig. 2a). The top 20 differentially expressed PCG are listed in Table 1.
To elucidate which biological processes and signalling pathways are associated with differential response to pRT, we have applied a gene set enrichment analysis using Ingenuity Pathways Analysis software. The bulk of differentially expressed genes were organised in two top scoring gene networks: lipid metabolism, molecular transport and small molecule biochemistry (Network 1) and cell cycle, DNA replication, recombination, and repair (Network 2) (Additional file 7: Figure S3). Specifically, differentially activated canonical pathways (Fisher's test p < 0.05) between pRT responders and nonresponders, included cell cycle control of chromosomal replication, and pathways related to glucose metabolism and de novo NAD biosynthesis (Fig. 2b). Similarly, molecular functions significantly enriched within the list of differentially expressed genes included cellular development, cell growth and proliferation, cell cycle and functions related to cell morphology, movement, assembly and organisation in addition to metabolic processes (Fig. 2c).

Candidate gene validation by qRT PCR
Even after controlling for false discovery rate, due to the small size of the discovery cohort and technical challenges related to FFPE-based microarray hybridization, several associations with the outcome of interest may have occurred due to chance alone. Therefore, we were selected for validation based on the significance threshold, high expression in tumour tissue and plausible biological function. Samples in which over 70% of assays have failed were discarded, retaining 42/60 samples for statistical analysis (Fig. 1b). Considering the cost of missing a potentially interesting gene (false negative) versus the low cost of further external validation to discard any false positive calls, we decided to use a lax False Discovery Rate of 20%. Differential gene expression between responders (N R = 30) and non-responders (N NR = 12) was observed for APOA1, MAP3K4 and MMP14 genes with over 2-fold downregulation in pRT-responsive tumours ( Table 2).

External validation of candidate gene predictive value for radiotherapy response
To evaluate whether expression of APOA1, MAP3K4 and MMP14 genes has an impact on patient survival in the postoperative setting, we analysed distant metastasis-free survival (DMFS) using a well-characterized cohort of 344 lymph node-negative breast patients undergoing surgical treatment, with or without postoperative radiotherapy (Erasmus cohort) (Fig. 1b). The Erasmus cohort was chosen to eliminate potential confounding effects of systemic therapy, and to asses weather these genes are radiation-specific (predictive) or not (simply prognostic) in a stratified analysis. Of the three genes tested in a subgroup analysis, only low MAP3K4 expression (< 7.94) was significantly associated with better DMFS (HR = 2.41 [95%CI:1.37-4.24], p = 0.002) in 282 patients treated with both surgery and RT, but not for those 62 patients  (Fig. 3). After controlling for age, steroid receptor status, T-stage and menopause status as potential confounders in a multivariable analysis MAP3K4 remained an independent predictor of DMFS (HR = 2.54, [95%CI:1.42-4.55], p = 0.002) in RT-treated patients. (Table 3). To check for the presence of biological interaction between MAP3K4 levels and radiotherapy, we calculated hazard ratios for each category combination and summary measures of effect modification on both multiplicative and additive scale. There was no evidence of statistical interaction on a multiplicative scale, while the combined effect of MAP3K4 levels and radiotherapy on additive scale exceeded that of each exposure alone with a 0.91 relative excess risk due to interaction [RERI] (Table 4).

Discussion
Treatment of LABC continues to be challenging with patients being at increased risks of locoregional recurrence, distant metastasis and reduced quality of life. Breast radiotherapy was shown to be effective in the locoregional control and provided benefit in distant metastasis-free survival and for downstaging the tumour in the preoperative treatment of LABC. However, not all patients achieve a satisfying response to radiotherapy. Clinically, a tumour is considered radioresistant when irradiation is unable to reduce its volume or when a recurrence occurs after a possible regression. Thus, it would be beneficial to identify biomarkers predictive of initial response to pRT that could be useful to predict clinical outcome in RT treated patients.
Here we explored gene expression profiles in pre-and post-RT tumour biopsies of LABC samples with  different clinical response to pRT. However, after microarray QC there were too few pre-and post-RT matched samples to be able to draw any statistically significant conclusions regarding changes induced by pRT. Therefore, we focused our further analysis on pre-RT biopsies only, to establish baseline differences in transcriptional profiles between patients achieving either complete or partial clinical response and those who did not respond to pRT. Among 192 significantly differentially expressed transcripts, only 89 corresponded to known protein coding genes and to 78 lncRNA. Although, lncRNAs were shown to have important roles in a broad range of biological processes such as at the level of post-transcriptional processing and transcriptional gene silencing [6,7] their biological function in breast cancer, and especially in relation to radiation therapy remains unknown.
Using IPA analysis, we explored which gene networks and signalling pathways are conferring radiation resistance/sensitivity in LABC tumour samples. Looking at protein interactions, most of the genes were organised in two major networks, one constituted of genes involved in lipid metabolism, molecular transport and small molecule biochemistry and the second related to cell cycle, DNA replication, recombination, and repair. These results further emphasise the importance of the cell cycle and proliferation state on cellular radiosensitivity, in line with the current body of knowledge [31,32]. Interestingly, when focusing on specific canonical pathways in addition to cell cycle control of chromosomal replication, pathways related to glucose metabolism and de novo NAD biosynthesis were significantly overrepresented. This profile potentially reflects a comparatively activated metabolic pathway for de novo synthesis of NAD+ in radio-resistant tumours. In addition to its numerous functions in redox reactions, NAD+ is a sole substrate for the activity of PARP enzyme in the repair of DNA single-strand brakes [33,34]. Several inhibitors of IDO1 gene, a rate-limiting enzyme for the NAD+ de novo synthesis from tryptophan found to be downregulated in radiosensitive LABC tumours, have shown a radiosensitising effect in in pre-clinical studies, making it an attractive target for drug-radiation combination therapy [35][36][37].
The second aim of the study was to identify potential biomarkers of response to radiotherapy. To this end, we selected 10 genes to be validated by an orthogonal assay in an independent subset of LABC tumours. Only 3 genes (APOA1, MAP3K4 and MMP14) were confirmed to be downregulated in radiosensitive tumours in the preoperative setting. These markers were further validated in an independent dataset of breast tumours treated only with mastectomy with or without postoperative radiotherapy. Only MAP3K4 gene was found to independently and specifically predict DMFS in radiotherapy-treated patients. Patients with high levels of MAP3K4 treated with RT had shorter DMFS than those with low MAP3K4 levels, due to an adverse interaction of high MAP3K4 expression with RT. MAP3K4 is a member of first layer of kinases of the MAPK signalling pathway that is activated by a variety of stimuli, including ionizing radiation, to mediate activation of transcription factors controlling differentiation, proliferation, cell growth and survival [38]. Increased expression of MAP3K4 in breast tumours may confer radioresistance through augmented signalling for RT-induced DNA damage repair through G2 arrest thus aiding survival of irradiated cancer cells [39]. Therefore, detecting the levels of MAP3K4 expression in breast tumours, may be useful for the prediction of response to RT for tumour downstaging for breast conserving surgery in LABC. Furthermore, inhibition of Ras-Raf-MEK-ERK cascade was shown to increase radiosensitivity both in in vitro as well as in vivo studies rendering MAPK signalling as an attractive radiosensitising target [40][41][42].