EGR1 regulates cellular metabolism and survival in endocrine resistant breast cancer

About 70% of all breast cancers are estrogen receptor alpha positive (ER+; ESR1). Many are treated with antiestrogens. Unfortunately, de novo and acquired resistance to antiestrogens is common but the underlying mechanisms remain unclear. Since growth of cancer cells is dependent on adequate energy and metabolites, the metabolomic profile of endocrine resistant breast cancers likely contains features that are deterministic of cell fate. Thus, we integrated data from metabolomic and transcriptomic analyses of ER+ MCF7-derived breast cancer cells that are antiestrogen sensitive (LCC1) or resistant (LCC9) that resulted in a gene-metabolite network associated with EGR1 (early growth response 1). In human ER+ breast tumors treated with endocrine therapy, higher EGR1 expression was associated with a more favorable prognosis. Mechanistic studies showed that knockdown of EGR1 inhibited cell growth in both cells and EGR1 overexpression did not affect antiestrogen sensitivity. Comparing metabolite profiles in LCC9 cells following perturbation of EGR1 showed interruption of lipid metabolism. Tolfenamic acid, an anti-inflammatory drug, decreased EGR1 protein levels and synergized with antiestrogens in inhibiting cell proliferation in LCC9 cells. Collectively, these findings indicate that EGR1 is an important regulator of breast cancer cell metabolism and is a promising target to prevent or reverse endocrine resistance.


INTRODUCTION
Resistance to endocrine therapy is a major clinical problem for the management of estrogen receptor positive (ER+) breast cancers. ER+ tumors comprise 70% of all breast cancer cases. Antiestrogens such as tamoxifen, ICI 182,780/faslodex/fulvestrant (ICI) or aromatase inhibitors (AI) are widely used endocrine therapies but little is known about the complex cellular pathways that contribute to resistance [1,2]. Deregulation of metabolic pathways regulated by oncogenes such as MYC has been implicated in endocrine resistant breast cancer [3][4][5].
However, to understand the systems-level changes in endocrine resistance, biologically relevant interactions between genes and metabolites need to be identified and validated. Using paired cell lines that are either sensitive or resistant to antiestrogens we generated and integrated data from transcriptomics (microarray analysis) and metabolomics (GC/MS and UPLC/MS). Within our genemetabolite integrated model, we selected to further study the role of EGR1 (early growth response 1), a gene that is known to be deregulated in some cancers [6].
EGR1 is an immediate-early gene induced by estrogen, growth factors, or stress signals, and can exhibit both tumor Research Paper www.impactjournals.com/oncotarget suppressor and promoter activities. A nuclear phospho-protein and transcription factor that can promote cell proliferation and cell death [7,8], EGR1 can be induced by external stimuli, with its induction being either transient or sustained. Such a diverse array of functions is achieved through differential regulation of EGR1 expression and its selection of target genes [7]. A highly conserved DNA-binding domain on EGR1 targets the GC-rich consensus sequence GCG (G/T) GGGCG. Transcriptional activity of EGR1 is further regulated by NAB-1 and NAB-2 (NGF-I A-binding proteins) [9,10]. EGR1 can promote growth of some hormone regulated cancers including prostate cancer [11]. EGR1 mediated signaling is important for the normal development of female reproductive organs [12] but its precise role in breast cancer remains unclear. The EGR1 gene is deleted in some ER-negative breast tumors [13]. In ER+ breast cancer cells, EGR1 is induced by estrogen treatment following raf-1 kinase activation [14] and is inhibited with acquired resistance to ICI [15]. EGR1 has been associated with sensing cellular glucose levels [16,17], in fatty acid metabolism and inflammation [18] in various cells. In this study, we investigated the role of EGR1 in endocrine resistance to validate an integrated model consisting of differentially expressed genes and metabolites in endocrine resistant breast cancer cells. Decreased levels of EGR1 in ER+ breast cancer cells and human tumors correlated with decreased sensitivity to antiestrogens. However, sustained inhibition of EGR1 with siRNA or tolfenamic acid (TOLE) suppressed the growth of endocrine resistant breast cancer cells and interacted synergistically with both 4-hydroxytamoxifen (hereafter referred to as TAM), major active metabolite of tamoxifen, and ICI in inhibiting cell proliferation. EGR1 inhibition also disrupted fatty acid metabolism in endocrine resistant cells. Overall, our gene-metabolite integrated model suggests a novel role for EGR1 in regulating cellular metabolism in endocrine resistant breast cancer.

Analysis and integration of transcriptomics and metabolomics data from endocrine sensitive and resistant breast cancer cell lines
Comparing LCC1 and LCC9 cell lines yielded 4,010 unique differentially expressed genes (DEGs), 46 unique m/z values for metabolites from the Metabolomics Shared Resources Core (MSRC) analysis conducted at Georgetown University Medical Center and 12 identified metabolites from the Metabolon analysis at an FDR threshold of 0.01 for the first two analyses, respectively; 0.05 for the third analysis. (Figure 1A and 1B). A heatmap was generated using all differentially expressed genes. Data visualization by principal component analysis (PCA) analysis was performed for both the MSRC and Metabolon metabolome analyses. While there was inherent noise, LCC1 and LCC9 cells were clearly separated in both plots. Figure 1C shows the integrated network of differentially expressed genes and putative metabolites in LCC9 cells (endocrine resistant) using 300 DEGs, 46 unique m/z values for metabolites from the MSRC analysis and 11 identified metabolites from the Metabolon analysis. The genes and metabolites present in the integrated network are also presented in Table 1. Based on our metabolomics analysis, glutamate and prostaglandin levels were significantly higher in LCC9 compared with LCC1 cells (log2 FC=1.522, p-value = 0.00011, q-value = 0.00574; respectively log2 FC=1.047, p-value = 0.00019, q-value = 0.00759). Estradiol and endothelin along with three genes (PTGS1, PTGS2 and GRM7) were added to the network based on predictions from the STITCH database.
Results from the pathway analyses are presented in Table 2. Table 2A represents results using the significant metabolites, Table 2B using the top 300 genes, and Table 2C using both metabolites and genes. In particular, we noted the D-glutamine and D-glutamate metabolism pathway (Table 2A, p-value = 0.001, q-value = 0.052), several key signaling pathways (Table 2B, and prostaglandin synthesis and regulation (Table 2C, p-value < 0.001, q-value = 0.002 when combining genes and metabolites).

Decreased EGR1 expression correlates with decreased responsiveness to antiestrogens in human breast tumors
To determine whether EGR1 expression was associated with disease free survival, we used publicly available gene expression datasets (see Methods; Figure 2; Table 3 [20] (Figure 2A and 2B). Furthermore, in Miller et al., GSE20181 [21], pre-treatment versus 90 days post-treatment comparisons for treatment with the aromatase inhibitor Letrozole showed significantly increased levels of EGR1 expression (p < 0.0001) only in the responder group ( Figure 2C).

EGR1 regulates cell proliferation in both endocrine sensitive and resistant breast cancer cell lines
To elucidate the role of EGR1 in endocrine responsiveness, we inhibited (RNAi) or overexpressed the EGR1 cDNA in LCC1 and LCC9 cells, followed by www.impactjournals.com/oncotarget . Figure 3A and 3B shows Western blot data that confirm the EGR1 protein levels with knockdown (EGR1-siRNA) or overexpression (EGR1-cDNA). Figure 3C and 3D show graphs of EGR1 protein levels normalized to actin protein levels from three independent experiments where LCC1 and LCC9 cells were transfected with either EGR1 siRNA, cDNA, or their respective controls. Under control conditions (control siRNA or empty vector/EV), the EGR1 protein levels were 1.25-fold higher in LCC1 cells compared with that in LCC9 cells. Thus, we confirmed the model prediction in Figure 1C that endogenous EGR1 expression is higher in LCC1 cells and lower in LCC9 cells. In EGR1 siRNA transfected cell, EGR1 protein levels decreased 2.5-and 3.8-fold in LCC1 and LCC9 cells, respectively, compared with control siRNA transfected cells. In EGR1 cDNA transfected cell, EGR1 protein levels increased 1.4-fold and 2.5-fold in LCC1 and LCC9 cells, respectively, compared to that in empty vector (EV) transfected cells. Independent of antiestrogen treatment, transfection with EGR1-siRNA significantly reduced cell proliferation in both cell lines within 48 h compared with control-siRNA ( Figure 3D). To determine whether EGR1 siRNA changed cell viability, we studied changes in apoptosis and necrosis in LCC1 and LCC9 cells transfected with either control or EGR1 siRNA for 48 h. Figure 3E shows significant decrease in cell viability in both LCC1 and LCC9 cells transfected with control siRNA compared with EGR1 siRNA. Transfection with EGR1-cDNA did not initially affect LCC1 or LCC9 proliferation ( Figure 3F). However, 5-days post-transfection, EGR1-cDNA transfected LCC1 and LCC9 cells each exhibited a significant decrease in proliferation compared with their respective EV-transfected controls ( Figure 3G). LCC1 cells, at 5-days post-transfection with EGR1-cDNA and treated with ICI, showed a modestly additive growth inhibition relative to EV control cells. At 5-days post-transfection with EGR1-cDNA and E2 treatment, LCC1 cells showed a significant decrease in E2 response compared to their EV controls. Thus, some basal level of EGR1 protein expression may be essential for the survival of both endocrine sensitive and resistant cells, whereas changes beyond this base level in sensitive cells determines their responsiveness to E2.

EGR1 regulates fatty acid metabolism in endocrine resistant breast cancer cells
To determine the role of EGR1 in affecting cell metabolism in LCC9 cells, we transfected cells with either the EGR1-siRNA or control siRNA for 48 h, or with the EGR1 cDNA (EGR1-cDNA) and the respective EV controls. Five biological replicates were used for each group. Metabolomics analysis was performed by Metabolon. Data analysis followed the same steps as for the LCC1/LCC9 comparison. Negative correlation between the fold-changes from the two experiments indicated good global agreement between the knockdown and overexpression approaches ( Figure 4). Applying an FDR cut off = 0.1 yielded 18 metabolites for the siRNA experiment; 15 of these metabolites had HMDB IDs (Table 4). Pathway analysis of these 15 metabolites was done using MetaboAnalyst [22] and IMPaLA: Integrated Molecular Pathway Level Analysis [23] (Table 5). No metabolites reached statistical significance for the EGR1 overexpression analysis. Several metabolites that were significantly regulated following EGR1 knockdown including 7-hydroxycholesterol and acetyl CoA, had fold-changes in the opposite direction when EGR1 was overexpressed.
Results from the siRNA and cDNA experiments suggest that disrupted EGR1 expression may have a subtle impact on the metabolic profile of LCC9 breast cancer  Figure 1C. The metabolite names are putative metabolites, so they could be one of the many annotations obtained. The metabolites marked with * were from Metabolon and experimentally validated.
cells. However, since cell proliferation was significantly reduced in with EGR1-siRNA (at 48 h; Figure 3B) and with EGR1-cDNA (at 5-days; Figure 3D), these seemingly subtle metabolic changes at 48 h post-transfection likely underestimate their ability to affect cell phenotype.
Pathway analysis using the significant metabolites from the siRNA experiment is presented in Table 5. Two of the pathways implicate fatty acid metabolism (biosynthesis of unsaturated fatty acids, transport of fatty acids). Fatty acids are a critical source of energy for mitochondrial oxidation and cellular ATP generation. Silencing of EGR1 was accompanied by high levels of glycerol and multiple monoacylglycerols such as 1-myristoylglycerol that may reflect an increase in complex lipid hydrolysis. Consequently, long chain fatty acids such as palmitoleate and medium chain fatty acids including caprylate and heptanoate were also elevated compared to control-siRNA cells. Fatty acid availability may ultimately alter mitochondrial β-oxidation.

TOLE down-regulates EGR1 and sensitizes resistant cells to antiestrogens
Since EGR1 is an essential regulator of cell survival, we tested the effect of TOLE, a nonsteroidal anti-inflammatory drug (NSAID) that induced cell     death in an EGR1-dependent manner in colorectal cancer cells [24]. Western blot analysis of whole cell lysates from LCC1 cells ( Figure 5A, Figure 5C). Together, these data indicate that antiestrogen treatment can differentially affect EGR1 levels in endocrine sensitive versus resistant cells. Furthermore, TOLE can downregulate EGR1 levels and re-sensitize endocrine resistant cells to antiestrogens.

DISCUSSION
Molecular adaptations that lead to drug resistance in cancer cells are largely dependent on cellular context and the nature of the stress signal. To determine the pathways that promote endocrine resistance, we integrated gene expression data with metabolite concentrations studying only those signals that were significantly changed in resistant breast cancer cells (LCC9) compared with sensitive cells (LCC1). The resulting model implicated EGR1 as a gene that is downregulated in endocrine resistant cells and proposed an increased activation of the glutamine and arachidonic pathways ( Figure 1C; Tables 2 and 3). Interestingly, pathway analysis of top DEGs showed that a number of pathways associated with immune response were significantly altered in endocrine resistant LCC9 cells. It is still unclear how immune response genes are regulated in ER+ breast cancer cells and tumors. In breast cancer patients, depression in cellular immunity was associated with resistance to endocrine therapy [25]. In vivo models of mammary tumors suggest a role of immune-associated genes in antiestrogen resistance [26,27]. The signaling interactions between cancer cells and the tumor microenvironment remain to be elucidated. In this study, to validate our gene-metabolite integration model, we tested whether alteration of EGR1, which is down-regulated in LCC9 cells [15], changed endocrine responsiveness or survival in resistant cells. We also asked whether gene expression data from human tumors treated with endocrine therapy showed a correlation between higher EGR1 levels and a more favorable prognosis (Figure 2A). EGR1 levels are high in some prostate cancers [28,29], Wilm's tumors [30], and melanoma cells bearing oncogenic B-RAF mutation [31] compared to normal tissue. An array of stress stimuli including radiation, chemotherapy, or hypoxia can alter EGR1 levels and the nature of a response is determined by whether EGR1 transcriptionally up-or down-regulates specific target genes [6,7,32]. While the precise role of EGR1 in cell survival remains unclear, disruption of endogenous levels of EGR1 can either inhibit growth or promote tumor   [33,34]. In prostate cancer cell lines, EGR1 inhibition decreased cell growth, induced apoptosis [35], and decreased expression of the pro-inflammatory chemokine interleukin 8 (IL8) in a NF-κB-dependent pathway [36]. EGR1 has been implicated in the acquisition of resistance to hormone therapy particularly through its role in the androgen receptor (AR) pathway [37,38]. In breast cancer, the role of EGR1 remains ambiguous.
In an ER+, carcinogen-induced (7,12-dimethylbenz(a) anthracene; DMBA), rat mammary tumor model, EGR1 levels in tumors were reduced relative to normal mammary tissues but then increased with TAM treatment [39]. In ER-negative breast cancer, the EGR1 gene is frequently deleted [13]. Disruption of the ER signaling pathway can affect EGR1 levels in ER+ breast cancer cells and tumors [14,39]. Interestingly, overexpression of EGR1 dramatically reduced E2-mediated proliferation in these cells. In a model to identify topological and temporal effects of E2 regulatory networks in MCF7 cells, EGR1 was identified as a mediator of some of the late responses to E2 [40]. Thus, EGR1 levels in breast cancer cells may be closely regulated by a functional ER pathway.
A role for EGR1 in regulating cellular metabolic pathways has been reported in several disease models including cancer. Several cellular metabolic pathways were altered with perturbation of EGR1 levels in LCC9 cells (Table 5), particularly molecules associated with fatty acid metabolism. EGR1 and the lipogenic enzyme fatty acid synthase (FASN) are elevated in tissues adjacent to prostate cancer; this relationship is used as a predictive marker of recurrence [41]. FASN activation is required for estrogen-mediated signaling in ER+ breast cancer cells [42] . EGR1 is an immediate-early prostaglandin E2 (PGE2) target gene that can mediate eicosanoid regulation of genes involved in the immune and inflammatory responses [43]. Together, these findings highlight the critical role of EGR1 in fatty acid metabolism.
In LCC9 cells, combining TOLE and antiestrogens synergistically inhibited cell growth ( Figure 5B). TOLE can inhibit synthesis of prostaglandins and has been used for treating migraines [44]. In a triple-negative breast cancer xenograft model, treatment with TOLE resulted in a significant reduction in tumor volume over 5 weeks compared to treatment with vehicle alone [45] The pathway analysis was performed using tool http://impala.molgen.mpg.de/ . The pathways enriched at q-value < 0.1 are shown.
inhibit cell growth in cancer cells through cyclooxygenaseindependent pathways including inhibition of ErbB2 expression [46], activation or ATF3 [47], or induction of NSAID-activated gene-1 (NAG-1) [48] and EGR1 [24]. Contrary to the latter study, in our breast cancer cells, treatment with TOLE inhibited EGR1 protein levels in LCC9 cells and synergized with ICI and TAM ( Figure 5B). Thus, cellular context is likely a key determinant of the outcomes of EGR1 action. Overall, we present a comprehensive model incorporating the differential expression/concentration of genes and metabolites that may interact to determine breast cancer cell fate in the response to select endocrine therapies. We tested our model by further elucidating the role of EGR1 in endocrine sensitive and resistant breast cancer cells. Our data suggest that although EGR1 levels are significantly reduced in endocrine resistant cells, it is an essential driver of cell survival and metabolic pathways such as fatty acid metabolism. Targeting EGR1 with TOLE may be an effective therapeutic strategy in some endocrine resistant breast cancers. Furthermore, EGR1 levels in human breast tumors may be useful as a favorable prognosis marker in ER+ breast cancer.

Cell proliferation and viability
For determination of cell density, cells were plated in 96-well plates at 5 × 10 3 cells/well. At 24 h, cells were treated with specified drugs for 48 h (or otherwise indicated). After treatment, media were removed, and plates were stained with a solution containing 0.5% crystal violet and 25% methanol, rinsed, dried overnight, and re-suspended in citrate buffer (0.1 M sodium citrate in 50% ethanol). Intensity of staining, assessed at 570 nm and quantified using a V Max kinetic microplate reader (Molecular Devices Corp., Menlo Park, CA), is directly proportional to cell number [50,51]. For assessing cell viability and cell death (apoptosis and necrosis), cells were treated for 48 h, and stained with an Annexin V-fluorescein isothiocyanate and propidium iodide, respectively (Trevigen, Gaithersburg, MD).

Transfections with EGR1 siRNA or cDNA
Cells were plated at 60-80% confluence. EGR1 (10 nM of 3 unique 27mer siRNA duplexes; Origene, Rockville, MD) or their respective control siRNA, were transfected using the RNAiMAX (Invitrogen) transfection reagent. For EGR1 overexpression, EGR1 cDNA (catalog #SC128132) was purchased from Origene and transfected with TransIT-2020 (Mirus). Cells were lysed at 48 h posttransfection and subjected to Western blot analysis or cell number assay as described above. Antiestrogens, 100 nM ICI or TAM, or vehicle (0.02% ethanol) was added to the transfected cells at 24 h and treatment was allowed for the time-points indicated. For 5-day treatments, cells were retreated with the indicated drugs in fresh cell culture media at day 3.

Generation and integration of transcriptomic and metabolomic data from LCC1 and LCC9 cells Transcriptome data
We obtained and analyzed gene expression and untargeted metabolomics data from antiestrogen sensitive (LCC1) or antiestrogen resistant (LCC9). Microarray analysis was performed using three biological replicates from LCC1 and three biological replicates from LCC9 using Affymetrix HG U133 Plus 2.0 microarray at our Genomics and Epigenomics Shared Resources. Briefly, total RNA was extracted using the RNeasy kit (Qiagen, Valencia, CA, USA). RNA labeling and hybridization were performed according to the Affymetrix protocol for one-cycle target labeling. For each experiment, fragmented cRNA was hybridized in triplicates to Affymetrix GeneChip HG-U95 arrays (Affymetrix, Santa Clara, CA). Affymetrix data analysis included pre-processing of the probe-level Affymetrix data (CEL files).

Metabolomics data
Metabolomics analysis was two-part: we sent cell samples to both our in-house MSRC and to Metabolon Inc. MSRC samples were five biological replicates from each of the two groups, each with two technical replicates. Metabolon samples were six biological replicates from each of the two groups. LC-MS was used to analyze the MSRC samples; both LC-MS and GC-MS were used by Metabolon.
For the MSRC protocol, metabolite extraction was performed as described by Sheikh et al. [52]. Briefly, the residual pellet was resuspended in 200 μL of solvent A (98% water, 2% ACN and 0.1% formic acid) for Ultra-performance liquid chromatography-electro-spray ionization quadrupole-time-of-flight mass spectrometry Data were acquired in centroid mode from 50 to 850 m/z in MS scanning. Centroided and integrated mass spectrometry data from the UPLC-TOFMS was processed to generate a multivariate data matrix using MarkerLynx (Waters).
For Metabolon, samples were prepared using the automated MicroLab STAR ® system (Hamilton Company, Reno, NV). A recovery standard was added prior to the first step in the extraction process for quality control (QC) purposes. Proteins were precipitated with methanol to remove protein, dissociate small molecules bound to protein or trapped in the precipitated protein matrix, and to recover chemically diverse metabolites. The resulting extract was divided into five fractions: one for analysis by UPLC-MS/ MS with positive ion mode electrospray ionization, one for analysis by UPLC-MS/MS with negative ion mode electrospray ionization, one for LC polar platform, one for analysis by GC-MS, and one sample was reserved for backup. Samples were placed briefly on a TurboVap ® (Zymark) to remove the organic solvent. For LC, the samples were stored overnight under nitrogen before preparation for analysis. For GC, each sample was dried under vacuum overnight before preparation for analysis. The LC/MS portion of the platform was based on a Waters ACQUITY ultra-performance liquid chromatography (UPLC) and a Thermo Scientific Q-Exactive high resolution/accurate mass spectrometer interfaced with a heated electrospray ionization (HESI-II) source and Orbitrap mass analyzer operated at 35,000 mass resolution. The MS analysis alternated between MS and datadependent MS2 scans using dynamic exclusion, and the scan range was from 80-1000 m/z. The samples destined for analysis by GC-MS were were analyzed on a Thermo-Finnigan Trace DSQ fast-scanning single-quadrupole mass spectrometer using electron impact ionization (EI) and operated at unit mass resolving power. The scan range was from 50-750 m/z. Raw data files are archived and extracted as described below. The scope of the Metabolon LIMS system encompasses sample accessioning, sample preparation and instrumental analysis and reporting and advanced data analysis. It has been modified to leverage and interface with the MSRC information extraction and data visualization systems, as well as third party instrumentation and data analysis software.

Metabolite quantification and data normalization
For the MSRC metabolomics data, peaks were detected and quantified using the estimated area-under-thecurve of the LC/MS signals via the CentWave algorithm [53]. For the Metabolon data, an in-house software was used for detection and integration of peaks [54]. For studies spanning multiple days, a data normalization step was performed to correct variation resulting from instrument inter-day tuning differences. Essentially, each compound was corrected in run-day blocks by registering the medians to equal one (1.00) and normalizing each data point proportionately (termed the 'block correction'). For studies that did not require more than one day of analysis, no normalization was necessary, other than for data visualization. In certain instances, biochemical data were normalized to an additional factor (such as cell counts, total protein as determined by Bradford assay, osmolality) to account for differences in metabolite levels due to differences in the amount of material present in each sample.

Data analysis and integration of transcriptomics and metabolomics
Gene expression raw data were processed using the RMA algorithm as implemented in the R affy package [55], followed by a log2-transformation and analysis using a moderated t-test via the limma approach implemented in the limma package in Bioconductor [56]. MSRC metabolomics data were processed using the XCMS approach [53, 57,58]. Internal controls were used to standardize raw values. Intensity values were then log2transformed and quantile normalized [59], to avoid infinite value produced by 0 expression level during the log2transform phase, 1e-06 was introduced to replace all 0 in the expression profile and the two technical replicates were averaged. Putative isotopes were identified via the CAMERA package [60] and higher-weight isotopes removed. Finally, results were analyzed using limma along with surrogate variable analysis [61,62] to remove sources of variability unrelated to the comparison of interest. Metabolon data were processed by Metabolon as described above and log2-transformed, followed by the use of limma. Log2 transformation was performed for both data types to meet the normality assumptions for t-tests. All statistical tests were adjusted for multiple testing using the Benjamini-Hochberg approach to control the false discovery rate (FDR). Results are presented as both the p-values and the q-values (transformed p-values used for FDR control).
Metabolites significant at an FDR of < 0.01 and < 0.05 for the MSRC set and the Metabolon set, respectively, were annotated using the HMDB and Metlin databases [63,64]. Pathway analyses for the genes significant at FDR of < 0.01 were performed using Pathway Studio (http://www.elsevier.com/online-tools/pathwaystudio) and Enrichr [65] that allowed us to perform pathway analysis using the KEGG [66] and Reactome [67] databases. Pathway analyses for the top metabolites significant at an FDR of < 0.01 for our MSRC analysis and < 0.05 for the Metabolon analysis were performed using MetaboAnalyst [22] and IMPaLA: Integrated Molecular Pathway Level Analysis [23], respectively.
The top 300 DEGs according to their q-value were chosen for integration, along with the metabolites obtained after the filtering described above. Knowledge-based integration of these gene expression and metabolomics results was performed using STITCH [68]. This created a network of genes and metabolites using text mining and database curation sources. Based on the information obtained from these sources, each connection was given a confidence value. For our network, connections with atleast medium confidence were selected [68,69]. STITCH lists up to 50 articles that support each connection in the network. For each article, the PubMed ID and abstract was provided and the information therein was compared for consistency/relevance to the model solution. STITCH added nodes (genes or metabolites) to the predicted network to indicate an indirect connection between two input nodes. A network was created using the nearest 5 and nearest 10 interactions of the input genes and metabolites. The final network obtained in STITCH was downloaded, and formatted using Cytoscape [70] to overlay fold change values from our experiments and to visualize the network.
For EGR1 inhibition and overexpression experiments, five biological replicates were considered and the analysis was performed by Metabolon. The data analysis followed the same steps as that described above for the LCC1/LCC9 comparison. A relaxed FDR cut off of 0.1 was applied for each experiment. Since transfection with EGR1 siRNA in LCC9 cells resulted in significant reduction in cell proliferation at 48 h, the small metabolic changes at this time-point are likely sufficient to affect cell phenotype. The metabolites showing significant differences between the control group and the inhibition or overexpression group were considered in pathway analyses.

Estimates of relapse-free survival and EGR1 gene expression levels from public gene expression datasets
Publicly available datasets for gene expression from human ER+ breast cancer tumors that were treated with endocrine therapy were obtained: GSE17705 (Symmans et al., treated with Tamoxifen) [19], GSE6532, ER+ samples on GPL96 platform (Loi et al.,

Statistical analysis for cell proliferation experiments
Statistical analyses were performed using the Sigmastat software package (Jandel Scientific, SPSS, Chicago, IL). Where appropriate, cell growth under different conditions were compared using ANOVA with a post hoc t-test for multiple comparisons. Differences were considered significant at p ≤ 0.05. Nature of interaction between TOLE, TAM and ICI was defined by measuring the R-index (RI). The RI values were obtained by calculating the expected cell survival (S exp ; the product of survival obtained with drug A alone and the survival obtained with drug B alone) and dividing S exp by the observed cell survival in the presence of both drugs (S obs ). S exp /S obs > 1.0 indicates a synergistic interaction [71].