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-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 gene-metabolite 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 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.

RESULTS

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.

Figure 1

Analysis and integration of gene and metabolites in LCC1 (sensitive) and LCC9 (resistant) ER+ breast cancer cells

(A) Heatmap: R package Limma was used for microarray analysis of LCC1 versus LCC9 data; significantly different genes were selected (q-value < 0.1, fold change, FC > 2) to plot the heatmap; there were 3-biological replicates. (B) Principal component analysis: PCA analysis performed for the transcriptomics and metabolomics datasets - MSRC and Metabolon. (C) Integration of differentially expressed genes and putative metabolites comparing LCC9 (resistant) and LCC1 (sensitive) cells. Metabolites are shown as rectangular nodes, and genes as ellipses. Orange nodes are over-expressed; blue nodes are under-expressed; darker color represents a higher fold change (FC). Lowest log2 FC = –5.66, highest log2 FC = 5.89. Grey nodes are those added into the network based on prediction by STITCH. Edge thickness increases with the confidence of the connection as predicted by STITCH. Gene-metabolite connections are shown in grey lines, gene-gene connections are shows as purple lines, and metabolite-metabolite connections are shown as gold lines. EGR1 is significantly decreased in LCC9 cells (log2 FC= –2.33).

Table 1

Genes and metabolites from integrated network (Figure 1C)

Name	Expanded Name	Type	Probe/Metabolite ID	log2 fold-change	p-value	q-value
ESR1	ESR1	gene	205225_at	–4.036	2.36E-12	6.91E-09
LGALS3	LGALS3	gene	208949_s_at	–3.232	4.4E-11	3.55E-08
KITLG	KITLG	gene	211124_s_at	–2.666	5.47E-08	4.43E-06
SOCS3	SOCS3	gene	227697_at	–2.575	5.68E-10	1.96E-07
EGR1	EGR1	gene	201694_s_at	–2.333	6.25E-08	4.88E-06
IL1R1	IL1R1	gene	202948_at	–2.331	4.10E-08	3.71E-06
ABAT	ABAT	gene	209459_s_at	–2.219	2.44E-08	2.53E-06
STC2	STC2	gene	203438_at	–2.176	3.24E-10	1.29E-07
ABCG2	ABCG2	gene	209735_at	–2.133	7.61E-10	2.34E-07
LFNG	LFNG	gene	228762_at	–2.121	1.67E-09	4.12E-07
ABCC3	ABCC3	gene	208161_s_at	–2.061	4.81E-09	8.16E-07
ERBB4	ERBB4	gene	214053_at	–1.884	6.01E-09	9.71E-07
EDN1	EDN1	gene	218995_s_at	–1.823	2.36E-09	5.31E-07
SLC12A2	SLC12A2	gene	204404_at	–1.779	3.16E-09	6.38E-07
GHR	GHR	gene	205498_at	–1.734	2.62E-08	2.68E-06
DLK1	DLK1	gene	209560_s_at	1.829	3.2E-09	6.41E-07
CD36	CD36	gene	206488_s_at	1.906	4.86E-08	4.15E-06
LGALS3BP	LGALS3BP	gene	200923_at	2.065	3.38E-09	6.60E-07
AOX1	AOX1	gene	205083_at	2.14	3.54E-10	1.38E-07
SLC7A11	SLC7A11	gene	209921_at	2.231	1.77E-09	4.23E-07
GNAI1	GNAI1	gene	227692_at	2.866	2.59E-10	1.19E-07
EGR3	EGR3	gene	206115_at	2.975	1.21E-08	1.55E-06
RUNX2	RUNX2	gene	232231_at	3.02	4.77E-09	8.15E-07
CYP2B6	CYP2B6	gene	206754_s_at	3.161	1.2E-10	6.69E-08
DUSP4	DUSP4	gene	204014_at	4.272	4.56E-11	3.56E-08
HPGD	HPGD	gene	203913_s_at	5.252	2.16E-10	1.08E-07
GRM7	GRM7	gene	Added to the network
PTGS1	PTGS1	gene
PTGS2	PTGS2	gene
docosapentaeno.	Docosapentaenoic acid (22N-6)	metabolite	HMDB01976	–1.911	0.00001	0.00141
DGLA	8;11;14-Eicosatrienoic acid	metabolite	HMDB02925	–0.938	0.00024	0.00859
Mead_acid	5;8;11-Eicosatrienoic acid	metabolite	HMDB10378	–0.938	0.00024	0.00859
lysine	Lysine	metabolite	HMDB00182	0.872	0.00010	0.00521
pyroglutamate	Pyroglutamic acid	metabolite	HMDB00267	1	0.00012	0.00586
PGE1	Prostaglandin E1	metabolite	HMDB01442	1.047	0.00019	0.00759
PGFM	3,14-dihydro-15-keto PGF2a	metabolite	HMDB04685	1.047	0.00019	0.00759
prostaglandin	Prostaglandin D1	metabolite	HMDB05102	1.047	0.00019	0.00759
NMDA	N-Methyl-D-Aspartic acid	metabolite	HMDB02393	1.094	0.00006	0.00421
dihydrodipicol.	L-2;3-Dihydrodipicolinate	metabolite	HMDB12247	1.094	0.00006	0.00421
LysoPE_(16:0/0.	LysoPE(16:0/0:0)	metabolite	HMDB11503	1.294	0.00017	0.00712
lysoPC	LysoPC (17:0/0.	metabolite	HMDB12108	1.365	0.00002	0.00248
L-valine	L-valine	metabolite	HMDB00883	1.438	0.00013	0.00612
alpha-(methyla.	alpha-(methylamino)isobutyric acid	metabolite	HMDB02141	1.438	0.00005	0.00399
betaine	betaine	metabolite	HMDB00043	1.438	0.00005	0.00399
glutamate	glutamate	metabolite	HMDB03339	1.522	0.00011	0.00574
LysoPC_(17:0/0.	LysoPC(17:0)	metabolite	HMDB12108	1.795	0.00002	0.00248
hydroxybutyrate*	hydroxybutyrate	metabolite	HMDB00710	1.297	0.00079	0.03740
hypotaurine*	Hypotaurine	metabolite	HMDB00965	2.017	0.00001	0.00132
endothelin	Endothelin	metabolite	Added to the network
estradiol	Estradiol	metabolite

The table shows the genes and metabolites from 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.

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).

Table 2A

Pathway analysis of significant metabolites (from MSRC and Metabolon) performed using MetaboAnalyst (http://www.metaboanalyst.ca/), showing pathways with p-value < 0.05

Name of pathway	Number of significant metabolites in pathway / Number of metabolites in pathway	p-value	q-value
Arachidonic acid metabolism	7/62	< 0.001	0.020
D-Glutamine and D-glutamate metabolism	3/11	0.001	0.052
Lysine degradation	5/47	0.003	0.073
Arginine and proline metabolism	6/77	0.005	0.101
Sphingolipid metabolism	3/25	0.015	0.237
Glycerophospholipid metabolism	3/39	0.048	0.640

Table 2B

Pathway analysis of top 300 genes (according to q-value) using Reactome (www.reactome.org), showing pathways with p-value < 0.05

Name of pathway	Number of top genes in pathway/ Number of genes in pathway	p-value	q-value
Translocation of ZAP-70 to Immunological synapse	16/39	< 1.00E-10	< 1.00E-10
Phosphorylation of CD3 and TCR zeta chains	16/44	< 1.00E-10	< 1.00E-10
PD-1 signaling	16/45	< 1.00E-10	< 1.00E-10
Generation of second messenger molecules	16/57	< 1.00E-10	< 1.00E-10
Co-stimulation by the CD28 family	16/96	< 0.001	< 0.001
MHC class II antigen presentation	18/141	< 0.001	< 0.001
Downstream TCR signaling	16/123	< 0.001	< 0.001
Cytokine Signaling in Immune system	41/747	< 0.001	< 0.001
Interferon Signaling	24/291	< 0.001	< 0.001
Interferon gamma signaling	18/176	< 0.001	< 0.001
TCR signaling	16/145	< 0.001	< 0.001
ERBB2 Activates PTK6 Signaling	4/18	< 0.001	0.026
ERBB2 Regulates Cell Motility	4/19	< 0.001	0.029
Interleukin-19, 20, 22, 24	3/9	< 0.001	0.035
Downregulation of ERBB4 signaling	3/10	0.001	0.045
SHC1 events in ERBB2 signaling	4/25	0.002	0.064
Nuclear signaling by ERBB4	5/44	0.002	0.077
Signaling by PTK6	6/80	0.007	0.205
GRB2 events in ERBB2 signaling	3/20	0.008	0.245
Interferon alpha/beta signaling	8/140	0.009	0.252
PTK6 Activates STAT3	2/7	0.009	0.253
PI3K events in ERBB2 signaling	3/22	0.011	0.273
Termination of O-glycan biosynthesis	3/28	0.021	0.450
Growth hormone receptor signaling	3/29	0.022	0.450
Signaling by ERBB2	4/54	0.026	0.450
Activation of anterior HOX genes in hindbrain development during early embryogenesis	6/113	0.030	0.450
Activation of HOX genes during differentiation	6/113	0.030	0.450
Signaling by Interleukins	15/425	0.031	0.450
Constitutive Signaling by Aberrant PI3K in Cancer	5/85	0.032	0.450
Adaptive Immune System	31/1075	0.035	0.450
Immune System	52/1984	0.035	0.450
RA biosynthesis pathway	3/39	0.047	0.450
NCAM signaling for neurite out-growth	11/300	0.047	0.450
ABC-family proteins mediated transport	4/66	0.048	0.450

Table 2C

Pathway analysis of top 300 genes (according to q-value) and significant metabolites using http://impala.molgen.mpg.de/

Name of pathway	Source of pathway	Pathway analysis for genes			Pathway analysis for metabolites			Pathway analysis for genes and metabolites
		Number of top genes in pathway/ Number of genes in pathway	p-value	q-value	Number of significant metabolites in pathway/ Number of metabolites in pathway	p-value	q-value	p-value	q-value
Prostaglandin Synthesis and Regulation	Wikipathways	4/28	< 0.001	0.457	2/9	0.016	0.687	< 0.001	0.002
ABC-family proteins mediated transport	Reactome	4/36	< 0.001	0.457	2/10	0.019	0.746	< 0.001	0.006
Synthesis of Prostaglandins (PG) and Thromboxanes (TX)	Reactome	3/15	< 0.001	0.457	3/34	0.038	1	< 0.001	0.009
Arachidonic acid metabolism	Reactome	3/54	0.029	1	7/78	0.001	0.096	< 0.001	0.011
Transmembrane transport of small molecules	Reactome	15/594	0.006	0.63	9/178	0.014	0.618	< 0.001	0.024
Transport of inorganic cations/anions and amino acids/oligopeptides	Reactome	5/94	0.006	0.63	4/46	0.017	0.73	0.001	0.027
GABA synthesis, release, reuptake and degradation	Reactome	2/20	0.025	1	2/15	0.042	1	0.008	0.180

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) for ER+ human breast tumors treated with endocrine therapy (adjuvant tamoxifen or AI as the only systemic therapy). Kaplan-Meier estimates of relapse-free survival over time (rfs_t) showed that high EGR1 gene expression levels were significantly correlated with favorable prognosis in at least two different datasets where breast tumors were treated with tamoxifen: in Symmans et al., GSE17705 [HR=0.38 (0.21–0.69); p = 0.00083] [19] and Loi et al., GSE6532 (ER+ samples on GPL96 platform) [HR=0.62(0.4–0.95); p = 0.028] [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).

Figure 2

Lower EGR1 levels correlate with lower survival in ER+ breast cancer patients treated with endocrine therapy

(A) and (B) Kaplan-Meier plots were generated using the Symmans et al. and Loi et al. datasets to estimate the number of patients living over time post endocrine treatment (Tamoxifen) with indicated levels of EGR1 expression in their breast tumors; rfs_t (recurrence free survival time) (C) Pre-treatment vs. 90 days post-treatment (Letrozole) comparisons show significantly increased levels of EGR1 expression (p < 0.0001) only in the responder group.

Table 3

Gene expression public dataset for ER+ breast cancer used for correlating EGR1 expression and endocrine response

Dataset	Treatment	Duration	Sample_Size
Symmans et al.	Tamoxifen	5 years	298
Loi et al.	Tamoxifen	N/A	181
Miller et al.	Letrozole	0,10-14,90 day time-point	36 in each time-point

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 treatment with vehicle (ethanol), TAM, ICI or 17β-estradiol (E2). 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.

Figure 3

EGR1 expression regulate cell proliferation and viability in both endocrine sensitive and resistant ER+ breast cancer cells

(A) Western blot of LCC1 and LCC9 cells showing the effect of EGR1 knockdown (EGR1-siRNA) and its respective control (EGR1-control-siRNA) or EGR1 overexpression (EGR1) or its respective control, empty vector (EV). Cells were transfected with siRNA or cDNA plasmid for 72 h. EGR1 protein appeared as a doublet, perhaps due to phosphorylation. Actin was used as a loading control. (B–C) Quantification of EGR1 protein (normalized to actin) following transfection with EGR1-siRNA compared with control siRNA in LCC1 and LCC9 cells show 2.5- and 3.8-fold reduction, respectively, (B) EGR1 protein in LCC1 and LCC9 cells show 1.4- and 2-fold increase, respectively, with EGR1-cDNA compared with EV, (C, D) EGR1 knockdown in both LCC1 and LCC9 cells significantly decreased cell proliferation at 48 h regardless of TAM or ICI treatment (ANOVA, p < 0.001). (E) EGR1 knockdown significantly decreased cell viability in both LCC1 and LCC9 cells (ANOVA, p < 0.01; *p < 0.01 for cell death in EGR1-siRNA versus control-siRNA for respective cells lines) at 48 h. (F) and (G) EGR1 overexpression for 48 h followed by treatment with TAM or ICI for 3-days or 5-days, respectively. While EGR1 overexpression did not change cell proliferation of either LCC1 or LCC9 cell under control or treatment conditions at 3-days, at 5-days, EGR1 transfected LCC1 and LCC9 cells showed significant decrease in cell proliferation compared with respective cells transfected with EV. At 5-day transfection with EGR1 combined with E2 treatment showed a significant decrease in E2 response compared to EV control (ANOVA, p < 0.05).

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.

Figure 4

EGR1 knockdown in endocrine resistant cells disrupt fatty acid metabolism pathway

Correlation between estimated log2 fold changes for the EGR1 knockdown experiment (siEGR1 vs. siCtrlEGR1) and the estimated log2 fold changes for the EGR1 siRNA experiment (EGR1 cDNA vs. EV EGR1). The negative correlation indicates agreement a global agreement between the two approaches, as the direction of change is expected to be different when comparing the knockdown to the overexpression experiments.

Table 4

Metabolites that were significantly altered with EGR1-siRNA knockdown in LCC9 cells, with q-value < 0.1 with EGR1-siRNA versus EGR1-control siRNA in LCC9 cells

Name	HDMB ID	KEGG ID	Results from siRNA experiment			Results from cDNA experiment
			log2 fold-change	p-value	q-value	log2 fold-change	p-value	q-value
stearoyl-arachidonoyl-glycerophosphoinositol (1)			1.818	< 0.001	0.097	–0.006	0.988	0.999
1-arachidonoylglycerophosphoinositol	HMDB61690		2.829	< 0.001	0.097	–0.060	0.865	0.999
7-hydroxycholesterol (alpha or beta)			–3.116	< 0.001	0.097	0.361	0.392	0.999
desmosterol	HMDB02719	C01802	1.037	< 0.001	0.097	–0.399	0.121	0.999
N-acetylglucosamine	HMDB00215	C00140	1.108	0.002	0.097	–0.418	0.216	0.999
5-dodecenoate (12:1n7)	HMDB00529		1.329	0.002	0.097	–0.220	0.522	0.999
N-palmitoyl-sphingosine	HMDB04949		1.415	0.002	0.097	0.141	0.648	0.999
acetyl CoA	HMDB01206	C00024	–3.331	0.002	0.097	0.358	0.750	0.999
dihomo-linolenate (20:3n3 or n6)	HMDB02925	C03242	1.080	0.002	0.097	0.086	0.730	0.999
linoleate (18:2n6)	HMDB00673	C01595	1.354	0.003	0.097	0.016	0.959	0.999
erucate (22:1n9)	HMDB02068	C08316	1.172	0.003	0.097	0.234	0.454	0.999
1-oleoylglycerophosphoserine			1.107	0.003	0.097	–0.062	0.893	0.999
1-myristoylglycerol (1-monomyristin)	HMDB11561	C01885	1.937	0.003	0.100	–0.105	0.775	0.999
nicotinamide ribonucleotide (NMN)	HMDB00229	C00455	1.076	0.004	0.100	0.170	0.642	0.999
phosphopantetheine	HMDB01416	C01134	2.875	0.004	0.100	0.070	0.952	0.999
caprate (10:0)	HMDB00511	C01571	1.152	0.004	0.100	–0.122	0.693	0.999
arachidonate (20:4n6)	HMDB01043	C00219	1.019	0.004	0.100	0.176	0.493	0.999
uridine 5’-diphosphate (UDP)	HMDB00295	C00015	–1.032	0.004	0.100	–0.040	0.945	0.999

Table 5

Pathway analysis on significant metabolites with EGR1-siRNA versus EGR1-control siRNA in LCC9 cells

Name of pathway	Source of pathway	Number of significant metabolites in pathway/ Number of metabolites in pathway	p-value	q-value
Biosynthesis of unsaturated fatty acids - Homo sapiens (human)	KEGG	4/32	< 0.001	0.042
Linoleic acid metabolism - Homo sapiens (human)	KEGG	3/19	< 0.001	0.042
Signal Transduction	Reactome	6/169	< 0.001	0.042
Heparan sulfate/heparin (HS-GAG) metabolism	Reactome	3/21	< 0.001	0.042
Regulation of lipid metabolism by Peroxisome proliferator-activated receptor alpha (PPARalpha)	Reactome	2/5	< 0.001	0.042
Activation of Gene Expression by SREBP (SREBF)	Wikipathways	2/5	< 0.001	0.042
YAP1- and WWTR1 (TAZ)-stimulated gene expression	Wikipathways	2/5	< 0.001	0.042
Glycosaminoglycan metabolism	Wikipathways	3/27	< 0.001	0.042
Leishmaniasis - Homo sapiens (human)	KEGG	2/6	< 0.001	0.042
Circadian Clock	Wikipathways	2/6	< 0.001	0.042
triacylglycerol degradation	HumanCyc	3/29	< 0.001	0.042
Defective SLC26A2 causes chondrodysplasias	Reactome	3/29	< 0.001	0.042
Defective PAPSS2 causes SEMD-PA	Reactome	3/29	< 0.001	0.042
Defective B4GALT7 causes EDS_ progeroid type	Reactome	3/29	< 0.001	0.042
Defective B3GAT3 causes JDSSDHD	Reactome	3/29	< 0.001	0.042
Defective CHSY1 causes TPBS	Reactome	3/29	< 0.001	0.042
Defective CHST3 causes SEDCJD	Reactome	3/29	< 0.001	0.042
Defective CHST14 causes EDS_ musculocontractural type	Reactome	3/29	< 0.001	0.042
Defective B4GALT1 causes B4GALT1-CDG (CDG-2d)	Reactome	3/29	< 0.001	0.042
Defective CHST6 causes MCDC1	Reactome	3/29	< 0.001	0.042
Diseases associated with glycosaminoglycan metabolism	Reactome	3/29	< 0.001	0.042
Glycosaminoglycan metabolism	Reactome	3/29	< 0.001	0.042
Defective EXT2 causes exostoses 2	Reactome	3/29	< 0.001	0.042
Defective EXT1 causes exostoses 1_ TRPS2 and CHDS	Reactome	3/29	< 0.001	0.042
MPS IX - Natowicz syndrome	Reactome	3/29	< 0.001	0.042
MPS I - Hurler syndrome	Reactome	3/29	< 0.001	0.042
MPS II - Hunter syndrome	Reactome	3/29	< 0.001	0.042
MPS IIIA - Sanfilippo syndrome A	Reactome	3/29	< 0.001	0.042
MPS IIIB - Sanfilippo syndrome B	Reactome	3/29	< 0.001	0.042
MPS IIIC - Sanfilippo syndrome C	Reactome	3/29	< 0.001	0.042
MPS IIID - Sanfilippo syndrome D	Reactome	3/29	< 0.001	0.042
MPS IV - Morquio syndrome A	Reactome	3/29	< 0.001	0.042
MPS IV - Morquio syndrome B	Reactome	3/29	< 0.001	0.042
MPS VI - Maroteaux-Lamy syndrome	Reactome	3/29	< 0.001	0.042
MPS VII - Sly syndrome	Reactome	3/29	< 0.001	0.042
Mucopolysaccharidoses	Reactome	3/29	< 0.001	0.042
phospholipases	HumanCyc	3/30	< 0.001	0.044
sphingomyelin metabolism/ceramide salvage	HumanCyc	3/30	< 0.001	0.044
sphingosine and sphingosine-1-phosphate metabolism	HumanCyc	3/36	< 0.001	0.063
the visual cycle I (vertebrates)	HumanCyc	3/36	< 0.001	0.063
Transport of fatty acids	Reactome	2/9	< 0.001	0.070
Transcriptional Regulation of White Adipocyte Differentiation	Wikipathways	2/9	< 0.001	0.070
Regulation of Lipid Metabolism by Peroxisome proliferator-activated receptor alpha (PPARalpha)	Wikipathways	2/9	< 0.001	0.070

The pathway analysis was performed using tool http://impala.molgen.mpg.de/ . The pathways enriched at q-value < 0.1 are shown.

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 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, left panel) show decreased levels of EGR1 at 48 h following treatment with 100 nM TAM or ICI, or 50 µM TOLE. In LCC9 cells (Figure 5A, right panel), treatment with TAM and ICI increased EGR1 protein levels at 72 h compared with control cells. TOLE treatment decreased EGR1 levels in LCC9 cells. In both LCC1 and LCC9 cells, co-treatment with TOLE and either TAM or ICI decreased EGR1 levels. In LCC1 cells, treatment with 100 nM TAM or ICI alone significantly (p < 0.001) inhibited cell number; a combination of 50 µM TOLE and TAM or ICI significantly (p < 0.001) decreased cell number compared with the individual treatments but the interaction was not synergistic (Figure 5B). In LCC9 cells, treatment with TAM or ICI had no effect on cell number compared with vehicle, but TOLE treatment alone (p < 0.001) or in combination with either TAM (RI = 1.31; p < 0.001) or ICI (RI = 1.20; p < 0.001) synergistically reduced the number of cells within 72 h, implying a partial restoration of ICI sensitivity. To confirm that inhibition of cell number with TOLE was mediated through EGR1 downregulation, we treated LCC9 cells that were either transfected with control siRNA or EGR1 siRNA followed by treatment with increasing doses of TOLE (0–50 µM). EGR1 knockdown resulted in a significant decrease in TOLE-mediated inhibition in cell number at 25 and 50 µM (p < 0.05; 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.

Figure 5

TOLE decreased EGR1 protein in both sensitive and resistant cells and re-sensitize resistant cells to antiestrogens

(A) Western blot analysis of LCC1 and LCC9 cells, treated with vehicle, TOLE (50 μM), TAM (100 nM) or ICI (100 nM) or the combination for 72 h. In LCC1 cells, TOLE, TAM or ICI treatment decreased EGR1 protein levels. However, in LCC9 cells, antiestrogens increased but TOLE deceased EGR1 protein levels. Actin was used as a loading control. (B) Cell proliferation was significantly decreased in both LCC1 and LCC9 with treatment with TOLE at 72 h. Combination of TAM or ICI with TOLE did not show a significant interaction in LCC1 cells. However, cell proliferation was synergistically decreased in LCC9 cells treated with TOLE +TAM (RI = 1.31) or ICI+TOLE (RI = 1.20) within 72 h. ANOVA, p < 0.001; *p < 0.001 for specified treatment and cell line compared to vehicle. Dashed line denotes decrease in relative cell proliferation in each cell line with TOLE alone. (C) In LCC9 cells, knockdown of EGR1 with siRNA showed significant decrease in cell proliferation with 25 or 50 µM TOLE in LCC9 cells suggesting that TOLE-mediated EGR1 downregulation contributes to TOLE-induced decrease in cell proliferation in LCC9 cells. ANOVA, p < 0.05; *p < 0.05 for indicated concentration of TOLE in control-siRNA versus EGR1-siRNA.

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 progression [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]. TOLE can inhibit cell growth in cancer cells through cyclooxygenase-independent 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.

MATERIALS AND METHODS

Cell culture and reagents

LCC1 (sensitive) and LCC9 (ICI resistant and TAM cross-resistant) cells were established as previously described [49-51]. Cells were grown in phenol red–free IMEM (Life Technologies, Grand Island, NY; A10488-01) with 5% charcoal-stripped calf serum (CCS); this media contains 2 mM L-glutamine and ∼12 mM glucose. ICI182,780 (ICI) and 4-hydroxytamoxifen (TAM) were obtained from Tocris Bioscience (Ellisville, MO). Tolfenamic acid (TOLE) was purchased from Selleck (Houston, TX). All cells were authenticated by DNA fingerprinting and tested regularly for Mycoplasma infection. All other chemicals were purchased from Sigma-Aldrich.

Cell proliferation and viability

For determination of cell density, cells were plated in 96-well plates at 5 × 103 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 VMax 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 post-transfection 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 re-treated with the indicated drugs in fresh cell culture media at day 3.

Western analyses

Total protein (∼20 μg) was isolated from cells following 72 h treatment or vehicle control (0.02% DMSO or ethanol) for protein analysis as previously described [50, 51]. The following antibodies were used: EGR1 (Cell Signaling, Danvers, MA), and β-actin (Santa Cruz Biotechnology, Santa Cruz, CA).

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 (UPLC-ESI-Q-TOFMS) analysis. Mass spectrometry was performed on a Q-TOF Premier (Waters) operating in either negative-ion (ESI-) or positive-ion (ESI+) electro-spray ionization mode with a capillary voltage of 3200 V and a sampling cone voltage of 20 V in negative mode and 35 V in positive mode. The cone gas flow was 25 L/h, and the source temperature was 120°C. Accurate mass was maintained by introduction of LockSpray interface of sulfadimethoxine (311.0814 [M+H]+ or 309.0658 [M-H]−). 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 data-dependent 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-the-curve 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 log2-transformed and quantile normalized [59], to avoid infinite value produced by 0 expression level during the log2-transform 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/pathway-studio) 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 at-least 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., treated with Tamoxifen) [20] and GSE20181 (Miller et al., treated with Letrozole) [21]. Kaplan-Meier plots were generated using the Symmans et al. and Loi et al. datasets to estimate the number of patients living over time post endocrine treatment with indicated levels of EGR1 expression in their breast tumors. The dataset in Miller et al. was used to generate box plots to show difference in EGR1 expression in responders. Graphs were generated using tools in the R statistical programming language.

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 (Sexp; the product of survival obtained with drug A alone and the survival obtained with drug B alone) and dividing Sexp by the observed cell survival in the presence of both drugs (Sobs). Sexp/Sobs > 1.0 indicates a synergistic interaction [71].