Impact of Tumor Microenvironment and Epithelial Phenotypes on Metabolism in Breast Cancer Free

Highly proliferative tumor cells undergo fundamental changes in metabolism and nutrient usage to survive and progress (1), and metabolic transformation appears to be necessary for sustained proliferation (2). Much of the metabolic transformation is glucose-dependent, with invasive cancers exhibiting increased aerobic glycolysis (3) via the “Warburg Effect” (4). The Warburg effect occurs when there is a metabolic shift toward glycolysis, with increased cellular production of biomass, especially amino acids and nucleic acids. While the Warburg phenomenon has been investigated for more than 85 years, the metabolic interactions between stroma and epithelium are not well characterized, despite the important role of stroma in breast cancer biology (5–12).

The gap in our understanding of stromal–epithelial interactions was recently illuminated when it was shown that epithelial tumor cells induce oxidative stress in the normal stroma (9), in turn, leading to activation of NFκB and hypoxia-inducible-1α (HIF-1α) pathways in cancer-associated fibroblasts (CAF; ref. 11). Concomitant changes in inflammation, autophagy, mitophagy, and aerobic glycolysis are induced in the stroma, which then produces energy-rich metabolites (such as lactate and pyruvate) that are secreted and used by epithelial cells to generate ATP (11, 13). This bidirectional communication, with aerobic glycolysis in stroma fueling cancer growth is referred to as “Reverse Warburg Effect” (14–16), and its discovery established an important role for stroma in altering metabolism (9–11). However, differences in metabolomics by stromal–epithelial interaction and according to breast cancer subtype have not been widely investigated.

Breast tumors show substantial heterogeneity (17) with at least 5 distinct subtypes: luminal A, luminal B, HER2-positive, basal-like, and claudin-low breast cancers (18–20). In addition to their unique cell-autonomous characteristics [including unique mutation patterns, chemotherapy responses, and cellular phenotypes (21–23)], these subtypes have distinct interactions with the stroma. The interaction of basal-like breast cancer cells with fibroblasts leads to an increase in inflammatory cytokines and migratory behavior, whereas luminal cells in culture with fibroblasts show distinct gene expression and altered proliferation (12, 24). This evidence of subtype-specific interactions with surrounding stroma underscores the importance of studying epithelial and stromal factors simultaneously when evaluating tumor metabolism.

In the present study, we apply metabolomics and radiotracer metabolic studies to simultaneously investigate the role of tumor gene expression and stromal–epithelial interactions in tumor metabolism. Tissue samples were used to conduct gene expression analyses and samples were classified using multiple published signatures: the PAM50 signature, claudin-low signature, and in vivo wound response signature (17, 20, 25). These gene expression phenotypes were then compared with metabolomics classes identified by analysis of 379 metabolites. Complementary coculture experiments were carried out to evaluate subtype-specific metabolic interactions of stromal and epithelial cells experimentally. Our results illustrate that compared with luminal cancer cells or luminal-derived CAFs, basal-like cells and basal-like CAFs play a strong role in determining substrate metabolism in culture. Strong associations between tissue metabolic phenotype and tissue gene expression also suggest the importance of stroma in metabolic cancer microenvironments.

Thirty-one fresh-frozen breast tumor samples, 6 breast tissue samples from reduction mammoplasty, and 5 metastatic samples were obtained under Institutional Review Board (IRB)-approved protocols at the University of North Carolina at Chapel Hill (UNC-CH; Chapel Hill, NC). Patients were treated according to the standard of care dictated by disease stage, estrogen response (ER), and HER2 status. Table 1 describes patient demographics and tumor characteristics. All tissues were handled by snap freezing immediately after surgery, and RNA was isolated using a protocol as described by Hu and colleagues (26). RNA integrity was determined using the RNA 6000 Nano LabChip Kit and Agilent 2100 Bioanalyzer.

Microarrays experiments were carried out as described (26). Briefly, labeled cRNA was generated using Agilent's Low RNA Input Linear Amplification Kit. Cy5-labeled experiment samples were combined with Cy3-labeled reference (Stratagene Universal Human Reference spiked with 1:1,000 with MCF-7 RNA and 1:1,000 with ME16C RNA to increase expression of breast cancer genes) and hybridized to 4 × 44 K Agilent whole genome arrays or 244 K Agilent custom arrays. For 2 cases with duplicate arrays, the intraclass correlation coefficients (ICC) were 0.597 and 0.886, suggesting moderate-to-strong agreement between replicates. All tumor gene expression data are publicly available through the Gene Expression Omnibus (GSE6128 and GSE6130).

For each tumor, a single sample was analyzed for metabolites. Tissue (100 mg) was homogenized in a volume of water at 4 μL per mg of sample in a 2-mL cryovial with two 3-mm zirconium oxide beads, one 3:8 steel ballcone, and one 1:8 steel ballcone on the GenoGrinder at 1,000 strokes/min for 5 minutes. A 100 μL aliquot of the homogenate was placed in a 96-well deepwell plate along with two 3-mm zirconium oxide grinding beads per well. The plate was placed in the chiller on the robot, and a 400 μL volume of EtOAc/EtOH (1:1) with 2.5 μg/mL 2-fluorophenylglycine and 25 μg/mL d2-maleic acid and tridecanoic acid was added to each well. The plate was capped with a plate mat and samples were shaken on GenoGrinder at 675 strokes/min for 2 minutes. The plate was centrifuged on a Beckman GS-6R centrifuge at 3,200 rpm for 5 minutes at 4°C and then placed back on the Hamilton LabStar robot and the centrifugates were transferred to a second plate. A 200-μL volume of MeOH was placed in each well of the first plate, and the plate was shaken and centrifuged as before. The centrifugates were transferred as before and the same process was repeated using 200 μL volume of MeOH/H2O (3:1) and then repeated using 200 μL volume of DCM/MeOH (1:1). The samples were mixed by pipetting up and down several times and a 225 μL aliquot of each was transferred to each of two 250-μL autosampler vial inserts. All extracts were placed in the Zymark TurboVap 96 for concentration under nitrogen streams for 25 minutes. A second 225 μL aliquot of each sample was transferred to the appropriate insert, and the samples were placed in the Labconco Centrivap Concentrator for drying for 6.5 hours for drying. Samples were transferred to the Labconco Freezone 6 lyophilizer for further overnight drying. For liquid chromatography/mass spectrometry (LC/MS), a 40-μL volume of MeOH with 20 μg/mL D10 benzophenone was added to each sample. Samples were shaken with a Lab-Line Titer plate shaker for 5 minutes on setting 7 and then a 60 μL volume of 0.1% formic acid was added containing 5 μg/mL d3 leucine, 5 μg/mL dl-4-chlorophenylalanine, 5 μg/mL 4-bromo-dl-phenylalanine, and 0.5 μg/mL amitryptyline.

For gas chromatography/mass spectrometry (GC/MS), samples were derivatized to a final volume of 50 μL for GC/MS analysis using equal parts bistrimethyl-silyl-trifluoroacetamide and solvent mixture acetonitrile:dichloromethane:cyclohexane (5:4:1) with 5% triethylamine at 60°C for 1 hour. Three types of controls were analyzed in concert with the experimental samples: samples generated from pooled experimental samples served as technical replicates throughout the data set, extracted water samples served as process blanks, and a cocktail of standards spiked into every analyzed sample allowed instrument performance monitoring. Experimental samples and controls were randomized across the platform run. Raw data are included in Supplementary Table S1.

Metabolite levels with signal intensity greater than 10 dpi in both channels and at least 80% present data were selected and the red/green Lowess-normalized ratios for each gene were log₂-transformed. Missing data were imputed using k-nearest neighbors' imputation (with k = 10). This resulted in a complete dataset consisting of 379 metabolites, and each row/metabolite was median-centered. Data were analyzed by unsupervised significance analysis of microarray (SAM; ref. 27), and significant metabolites were clustered across all samples using average-linkage hierarchical cluster analysis, results were visualized using Java Treeview (28). SAM uses the variance structure of the data to compute expected distributions. Small fold changes were statistically significant by SAM analysis due to low variation between samples within a class. To test whether there was a significant trend in metabolite level for classes of metabolites [amino acids, carbohydrates/sugars, nucleic acids, and the tricarboxylic acid (TCA) cycle metabolites], we estimated the β-value and P value corresponding to linear regression of metabolite level on an ordinal variable equal to 1 for normal tissue, 2 for cluster 1 samples, 3 for cluster 2 samples, and 4 for metastases. Statistical tests of the trend were 2-tailed and SE was calculated (SAS version 9.2). Because these clusters emerged from SAM analyses where multiple comparisons were adjusted, no adjustment for multiple comparisons was made in the trend analyses. The objective of these analyses is to show whether there is evidence of a monotonic increase or decrease in metabolite level.

Using the probe set common to both array platforms (25), probes with signal intensity greater than 10 dpi in both channels and at least 80% present data were selected and the red/green Lowess-normalized ratios for each gene were log₂-transformed. Duplicate microarrays corresponding to the same patient sample were combined by averaging. Missing data were imputed using k-nearest neighbors' imputation (with k = 10). Gene annotations from each dataset were mapped to UniGene cluster IDs (UCIDs, Build 161) using the SOURCE database, and multiple occurrences of a UCID were collapsed by taking the median value for that ID within each experiment and platform. This resulted in a complete dataset consisting of about 2,800 genes, and each experiment was standardized to N(0, 1) and each row/gene was median-centered.

To identify genes that significantly changed by metabolic cluster (cluster 1 or 2), data were analyzed by 2-class SAM (27). Significant genes were clustered across all samples using average-linkage hierarchical cluster analysis, and results were visualized using Java Treeview (28). Significant genes were evaluated for ontologic enrichment using Ingenuity Pathway Analysis (IPA), with Benjamini–Hochberg (B–H) multiple testing correction. Significant functions and pathways were defined as those with B–H P < 0.05.

To characterize the biologic phenotypes of the metabolic clusters, gene expression in each sample was compared with an existing in vivo breast-derived wound response signature (25). The wound response signature was generated by comparing cancer-adjacent stroma-rich breast tissue to breast tissue from reduction mammoplasty, which identified an active wound signature in the cancer-adjacent tissue. The median-centered gene expression profile of each patient was evaluated for correlation with this signature by calculating Pearson correlation coefficients, using the method of Creighton and colleagues (29). Briefly, vectors corresponding to the genes in the wound response signature were constructed, with 1 assigned to upregulated genes and −1 assigned to downregulated genes. Sample arrays were filtered to retain only genes with interquartile range (IQR) of at least 0.8. A Pearson correlation coefficient was calculated comparing this standard vector to the vector of median-centered gene expression for each patient. Patients were classified as positive if the Pearson correlation coefficient was greater than zero and negative if the coefficient was less than zero. The association between metabolic cluster and wound response score was evaluated using a 2-tailed Fisher exact test (conducted in SAS version 9.2).

Primary CAFs were isolated from breast tissue of patients undergoing breast surgery for primary invasive breast carcinoma at UNC Hospital. Tissue specimens were procured under an IRB-approved protocol (LCC 0913) by the Lineberger Cancer Center Tissue Procurement Facility and stored in 10 mL of Dulbecco's Modified Eagle's Media (DMEM)/F12 with 10% FBS on ice until processing (up to 6 hours). Tissue was minced and transferred to a 15-mL conical tube with 9 mL of suspension media: DMEM/F12 (GIBCO) supplemented with 10% FBS (GIBCO), 1% Pen/Strep (Invitrogen), 2.5 μg/mL Amphotericin B solution (Invitrogen), 300 U/mL collagenase (Sigma), and 100 U/mL hyaluronidas (Sigma). Collagenase digestion occurred overnight at 37°C, and the sample was then centrifuged for 30 seconds at 100 × g. A top layer of hydrolyzed fat was discarded and supernatant was centrifuged at 200 × g for 3 minutes. Supernatant was transferred to a new tube and centrifuged at 400 × g for 5 minutes. The supernatant was discarded, and the remaining fibroblast pellet was resuspended in 7 mL of suspension media, as described earlier, and transferred to a T25 flask. Cells were grown at 37°C in a 5% CO₂, with media changed every 3 to 5 days and split into a T75 at 80% confluency.

SUM149 (basal-like) and MCF7 (luminal) breast cancer cell lines were obtained from American Type Culture Collection (30). Cell lines were maintained 37°C and 5% CO₂ in DMEM/F12 supplemented with 10% FBS and 50 units/mL penicillin/streptomycin. Cell lines were tested for mycoplasma by the Tissue Culture Facility at Lineberger Cancer Center.

Two types of cocultures were conducted to model the tumor metabolic microenvironment. First, we conducted direct cocultures, defined as a coculture where the 2 cell types are grown in physical contact, in the same well. Direct cocultures were used for all of the glucose assays plated at a physiologically relevant 1:2 ratio of cancer cells to fibroblasts in 24-well plates and were maintained for 96 hours, changing the media at 48 hours, before assays were conducted. Second, we conducted indirect cocultures or transwell cultures, where fibroblasts and cancer cells are grown separated by a membrane but in contact via soluble factors. Indirect cocultures were solely used to calculate individual cell growth rates for interacting stromal and epithelial cells. These studies were conducted by seeding 1 of the 2 cell types on the insert layer of Corning Transwell plates with 0.4-mm pore polycarbonate membranes, whereas the other cell type was grown in the bottom of the well. At 96 hours, 48 hours after media change, cells were harvested from top and bottom wells and counted using a Coulter Counter (Beckman Coulter, Inc.) to determine the ratio of cells. Identical starting stromal–epithelial ratios and analysis time points were used for both direct and indirect cultures. Two luminal and 1 basal-like primary CAF cell lines were used for all in vitro experiments, and all experiments were repeated in triplicate. For hepatocyte growth factor (HGF) inhibition, cells were treated with 0.5 μg/mL of anti-HGF antibody (Abcam) at time of plating, and media were not changed during the course of the experiment (48 hours).

Cells were plated for direct coculture at 1.0 × 10⁵ epithelial cells per well and 2.0 × 10⁵ fibroblasts/well into 24-well plates and allowed to grow for 96 hours, with media changed at 48 hours. Cells that were treated with anti-HGF antibody or sham were assayed at 48 hours. The cells were washed twice with 0.5 mL 37°C 1× PBS and then incubated for 10 minutes at 37°C in 0.5 mL/well KRH buffer (136 mmol/L NaCl, 4.7 mmol/L KCl, 1.25 mmol/L CaCl₂, 1.25 mmol/L MgSO₄, 10 mmol/L HEPES, pH 7.4) with 0.5 μCi/mL [³H]2-deoxy-d-glucose (PerkinElmer). Cells were washed twice with 0.5 mL iced 1× PBS with 20 mmol/L d-glucose and then lysed in 0.5 mL lysis buffer (0.025% SDS, 1% Triton X-100 in 1× PBS). Lysates were centrifuged at 20,000 × g for 5 minutes at 4°C, and the protein concentration was determined with a Bicinchoninic Acid (BCA) assay (Thermo Fisher Scientific). The rest of the lysate was then added to a scintillation vial filled with 4.5 mL EcoScint H (National Diagnostics) for scintillation counting. Because cell composition changes during the time of coculture (due to different growth rates for cancer cells and fibroblasts), expected metabolic responses for cocultures were computed using cell numbers at time of measurement and using measured activity levels. For each of n cocultures, the following formula was used to calculate expected metabolic activity level:

where m = metabolic activity and a = fraction of epithelial cells as estimated by cell counts for both cell types at a given time point after the start of coculture. Expected activity computed by this formula was compared with observed activity, and fold change was computed as a function of expected over observed metabolic activity. Statistical analyses of these data were conducted with 2-tailed t tests (e.g., all MCF7 cocultures vs. all SUM149 cocultures, or all LCAF cocultures vs. BCAF cocultures, etc.).

Cells were plated for direct coculture at 1.0 × 10⁵ epithelial cells per well and 2.0 × 10⁵ fibroblasts per well into 24-well plates and allowed to grow for 96 hours, with media changed at 48 hours. The cells were washed twice with 0.5 mL 37°C 1× PBS and then incubated for 120 minutes at 37°C in 0.5 mL/well KRH buffer (136 mmol/L NaCl, 4.7 mmol/L KCl, 1.25 mmol/L CaCl₂, 1.25 mmol/L MgSO₄, 10 mmol/L HEPES, pH 7.4) with 2 μCi/mL [¹⁴C] uniformly labeled d-glucose (PerkinElmer). Following incubation, 400 μL of KRH media from each well was placed into the upper well of an oxidation plate with 200 μL 1 mol/L NaOH in the lower well. The oxidation plate was sealed with a gasket and lid and the media well acidified by injecting 100 μL 70% perchloric acid. The oxidation plate was placed on an orbital shaker for 1 hour and then 150 μL of the NaOH from each well was placed into a scintillation tube with 4.5 mL EcoScint H for scintillation counting. The data are expressed as DPM/mg of protein/minute. Cells from this assay were washed twice with 0.5 mL iced 1× PBS with 20 mmol/L d-glucose and then lysed in 0.5 mL lysis buffer (30% KOH saturated with Na₂SO₄). Lysates were heated to 95°C for 10 minutes and centrifuged at 20,000 × g for 5 minutes at 4°C. The protein concentration was determined with a BCA assay (Thermo Fisher Scientific). To the remaining lysate, an equal volume of 100% ethanol was added to precipitate glycogen and samples stored at −20°C overnight. Samples were centrifuged at 20,000 × g for 5 minutes at 4°C and the supernatant removed. The pellets were washed in 1 mL of 70% ethanol and centrifuged at 20,000 × g for 5 minutes at 4°C. Supernatants were removed and pellets air-dried overnight. Pellets were then resuspended in 0.5 mL water and transferred to a scintillation vial with 4.5 mL EcoScint H for scintillation counting. The data are expressed as DPM/mg of protein/h fold change of observed/expected. For each of these cocultures, the expected activities of glucose oxidation and glycogen synthesis pathways were computed as described earlier and used to compute fold change relative to expected. Statistical analyses of these data were conducted as described earlier for glucose uptake.

Cells were plated at 1.0 × 10⁵ epithelial cells per well and 2.0 × 10⁵ fibroblasts per well into 24-well plates and allowed to grow for 48 hours, media were collected and stored in aliquots at −80°C and assayed for lactate content. Briefly, fresh NAD+ was added to assay buffer (175 mmol/L hydrazine sulfate, 68 mmol/L glycine, 2.9 mmol/L EDTA, 11.3 mmol/L NAD+, pH 9.5) just before assay. In a 96-well plate, 200 μL of assay buffer and 40 μL of lactate standard or media sample were combined in each well. The plate was then preread at 340 nm followed by the addition of 1U of lactate dehydrogenase in a 10 μL volume of water. The plate was mixed and read at 340 nm. The lactate concentration in each sample was determined from the standard curve. Computed expected levels of lactate were conducted as described for other metabolic assays earlier. Statistical analyses of these data were conducted as described earlier for glucose uptake.

Cells were harvested from culture, and protein was isolated and quantitated. Lysates were denatured by boiling with β-mercaptoethanol, and 30 μg of protein was electrophoresed on a 4% to 20% Tris-HCl Criterion precast gel (Bio-Rad) and transferred to a Hybond-P membrane (Amersham Biosciences) by electroblotting. The blots were probed with antibodies against the GLUT1 (Abcam) and β-actin (Cell signaling). Blots were washed 3 times with Tris-buffered saline supplemented with 0.1% Tween and then were probed with ECL anti-mouse IgG horseradish peroxidase (HRP)-linked whole antibody from rabbit (Amersham-GE Healthscience). Blots were rewashed, and detection was by enhanced chemiluminescence Western blotting detection system (Amersham-GE Healthcare). Relative GLUT1 protein concentration was quantified using ImageJ software, pixel intensity was used to measure the protein band of GLUT1 and divided by the intensity of the β-actin band.

Conditioned media from direct 1:1 cocultures (48 hours) was analyzed on a RayBio Human Cytokine Antibody Array 5 (80; Raybiotech) designed to detect 80 cytokines and chemokines. These glass arrays were used according to manufacturer protocol to measure cytokine and chemokine expression in the conditioned media from direct cocultures. Briefly, slides were blocked by incubation with blocking buffer at room temperature for 30 minutes and incubated with 100 μL of the sample at room temperature for 90 minutes. Membranes were washed and incubated with biotin-conjugated antibodies overnight at 4°C. Finally, the membranes were washed and incubated with fluorescent dye–conjugated streptavidin at room temperature for 2 hours. After final washing, slides were dried by centrifugation at 1,000 rpm for 3 minutes. Fluorescent signal was detected on a laser scanner using a cy3 (green) channel (excitation frequency 532 nm). Data for each cytokine were normalized to positive controls on the same slide to estimate relative protein expression. Each monoculture or direct coculture was analyzed in duplicate.

RNA was isolated from cultured cells using RNeasy Mini kit (Qiagen) following manufacturers' instructions. One microgram of total RNA was reverse transcribed into cDNA using iScript cDNA synthesis kit (Bio-Rad). Quantitative PCR (qPCR) was carried out with a dilution of cDNA equivalent to 100 ng RNA in 18 μL of master mix [10 μL SsoFast 2X Probes Supermix (Bio-Rad), 0.5 μL 18S-VIC, and 0.5 μL gene-specific Assay-On-Demand-FAM (Applied Biosystems), 7 μL water] was used in each well of the qPCR 96-well plate. The following primer/probe sets (Applied Biosystems) were used: GLUT1, 18S. Amplification conditions were as follows: 1 cycle of 95°C for 1 minute; 40 cycles of 95°C for 5 seconds, 60°C for 20 seconds. The best linear fit equation generated by the amplicon standard curve was used to determine the number of copies of GLUT1 in the cDNA generated. Relative transporter expression was normalized to the expression of 18s. Statistical analyses of these data were conducted with a 2-tailed t test (e.g., cocultures vs. monocultures, or cocultures with and without HGF antibody).

To identify metabolic subgroups of tumors, unsupervised hierarchical clustering was conducted on 379 metabolites across 31 breast tumor and 6 normal breast tissue samples. Patient characteristics for these 37 patients are described in Table 1. Two main clusters resulted (dendrogram shown in Fig. 1A and heatmap shown in Fig. 1C), one of which (cluster 1) was primarily normal tissue and less aggressive tumors (90% of luminal A tumors; 100% of normal breast tissues; Fig. 1B, colored boxes). The other cluster (cluster 2) was enriched for more aggressive tumor subtypes and metastatic tumors (60% of basal-like tumors; 66.7% of luminal B tumors; 100% of metastases). Cluster 2 tumors had significantly larger tumor size (P = 0.047) than those in cluster 1. While associations were not statistically significant, cluster 2 tumors were also substantially more likely to be ER-negative (Fig. 1B, gray scale boxes) and have poorly differentiated tumor grade.

Given the importance of stromal cells in altering glucose metabolism of cancer cells in previous literature (10, 31, 32), we hypothesized that metabolite subtypes would correlate with stroma-driven gene expression in the human tissue specimens. Using a published wound response signature identified in cancer-adjacent stroma-rich breast tissue when compared with breast tissue from reduction mammoplasty (25), we observed an association between metabolic cluster and wound response, with a higher percentage of tumors in cluster 2 (the more aggressive metabolic group) expressing activated wound response signature (Fig. 1D, P = 0.054). Of the tumors in cluster 1, 39% were positively correlated with the wound response signature, whereas 79% of tumors in cluster 2 were positively correlated with the signature. Interestingly, the association between metabolic subtype and wound response signature was stronger than the association with either ER status or breast tumor subtype (Table 1).

Four well-defined groups of metabolites defined the differences between metabolic clusters: amino acids, sugars, nucleic acids, and metabolites involved in the tricarboxylic citric acid (TCA) cycle. Principal component analysis (PCA) using metabolites in each of these 4 metabolite classes separated the cluster 1 and 2 tumors and delineated normal from metastatic samples (Fig. 2). The association between metabolite levels and metabolite cluster was strongest for amino acids (Fig. 2A; PC1, 85.7%; PC2 4.3%), followed by nucleic acids (Fig. 2B; PC1, 51.6%; PC2, 24.3%), TCA cycle (Fig. 2C; PC1, 54.2%; PC2, 17.2%), and carbohydrates/sugars (Fig. 2D; PC1, 32.8%; PC2, 13.6%), suggesting the relative importance of these types of metabolites in driving clustering and aggressive tumor subtype. Considering individual metabolites in each of these 4 categories, there were 19 amino acids that increased between normal samples and cluster 1 and between cluster 1 and non-metastatic tumors in cluster 2 (Table 2). Many of these metabolites increased further in metastatic tumors from different sites (skin, lung, liver, or brain). Several individual metabolites showed a trend from normal to cluster 1 tumors to cluster 2 non-metastatic tumors, including glucose (decrease), glucose-6-phosphate (increase), lactate (increase), and ribulose-5-phosphate (increase). Malic acid and fumaric acid, both metabolites in the TCA cycle, increased from normal to cluster 1 tumors to cluster 2 non-metastatic tumors, inversely proportional to citric acid levels. In summary, many metabolites that are crucial for glucose metabolism and the citric acid cycle were present at lower or higher levels in the more aggressive tumor cluster (cluster 2) than in cluster 1. A diagram illustrating the pathways affected by this metabolic shift is presented in Fig. 3.

Previous experimental studies focused on specific metabolic hallmarks of Warburg phenomenon recently identified a reverse Warburg effect, where stromal cells act as important drivers of the metabolic phenotype of cancers (10, 11, 14, 15, 33). Given this previous literature and the correlation between metabolic phenotype and wound response genomic signature, we conducted cocultures that modeled basal-like and luminal metabolic microenvironments to assess whether fibroblast characteristics, cancer cell characteristics, or interactions between the 2 were important in tumor metabolism. Basal-like microenvironments were modeled with SUM149 basal-like breast cancer cells and with basal-like CAFs (BCAFs). Luminal microenvironments were modeled with MCF7 luminal breast cancer cells and luminal CAFs (LCAFs). Results showed a complex interplay of stromal and epithelial cells in determining metabolic phenotype.

Consistent with the observation that more aggressive subtypes have more extreme metabolic phenotypes, SUM149 (basal-like) cells displayed higher glucose uptake than MCF7 (luminal) cells in monoculture (Supplementary Fig. S1A). However, SUM149 cells had even greater glucose uptake when grown with BCAFs (Fig. 4A), showing a 2.9-fold increase in observed relative to expected (P < 0.001). BCAFs also increased MCF7 glucose uptake (1.9-fold higher than expected, P = 0.001). In contrast, LCAFs suppressed glucose uptake in MCF7 cocultures but had no effect in coculture with SUM149 cells. Overall, BCAF cocultures had a 2.8-fold increase in glucose uptake when compared to LCAF cocultures (P = 0.04).

Once glucose is taken up, utilization includes glucose oxidation, glycogen synthesis, and lactate production, which are endpoints of oxidative metabolism, storage, and glycolysis, respectively. Basal-like cells had higher glucose oxidation in monoculture (Supplementary Fig. S1). However, glucose oxidation was suppressed in all luminal coculture conditions (MCF7 cells or LCAFs present, Fig. 4B). Increased glucose oxidation among SUM149 cells occurred only when these cells were cocultured with BCAFs (2.1-fold higher glucose oxidation than SUM149s with LCAFs, P < 0.001). This is interesting, suggesting that both stromal and epithelial factors contribute to glucose oxidation.

Glucose oxidation in the cell is balanced by non-oxidative glucose utilization, including glycogen synthesis. Stored cellular glycogen can promote cell survival in conditions of hypoxia (33). Glycogen synthesis was increased in all cocultures relative to monocultures (Fig. 4C); however, the increase in glycogen synthesis was most noticeable among SUM149 cocultures.

Finally, we assessed lactate as a measure of anaerobic glycolysis. Lactate production was higher in basal-like cancer cells and CAFs than in luminal cells in monoculture (Supplementary Fig. S1A); however, all coculture models showed about a 1.5-fold higher-than-expected lactate production that was not subtype-driven (Supplementary Fig. S1B). Considering all of these metabolites, coculture of basal-like epithelial cells with basal-like fibroblasts revealed the strongest phenotype of elevated glucose metabolism including glycolysis, oxidation, and storage.

Given the importance of stromal cells in altering glucose metabolism of cancer cells in previous literature (10, 32, 33) and our coculture data, we hypothesized that the expression of glucose transporters may be regulated by tumor–stromal interactions. GLUT1–4 are class I integral membrane proteins responsible for the transport of glucose. We focused on GLUT1 because it is responsible for basal levels of glucose uptake in all cells and glucose levels inversely regulate its expression (34). GLUT1 has also been shown to be the predominant glucose transporter in human breast carcinomas and mouse mammary carcinomas (35, 36). GLUT1 protein levels were higher in the more aggressive SUM149 monocultures than in MCF7 monocultures (Fig. 5A). Coculture conditions further increased GLUT1 protein levels in both cell types, showing a role for stroma in the regulation of glucose metabolism (Fig. 5A). To identify soluble factors and specific pathways regulating this change in GLUT1, cytokine protein arrays were conducted. These results indicated significantly higher levels of HGF in coculture conditions (Fig. 5B). Consistent with previous data in muscle (37), we hypothesized that HGF was playing a regulatory role in GLUT1 expression. We were able to completely block GLUT1 RNA expression by inhibiting HGF levels in the media, confirming our hypothesis (Fig. 5C). To further confirm this mechanism and its downstream effects on metabolism, we compared levels of glucose uptake in stromal–epithelial cocultures with and without HGF inhibition. Our results show 65.6% decrease in glucose uptake (P = 0.055) with HGF inhibition in SUM149:RMF cocultures (Fig. 5D).

Gillies and Gatenby (38) have argued that adaptations to resource scarcity are fundamental in the evolution of carcinogenesis and have shown that a variety of pathways dysregulate aerobic glycolysis in tumors (39). Consistent with metabolic adaptation during carcinogenesis, metabolite profiles can distinguish cancer from non-cancer (40). However, while it has been increasingly recognized that cancers evolve into distinct breast cancer subtypes, few studies have evaluated evolution of metabolic differences between breast cancer subtypes (39, 41–44). Our study shows that tumors differ from normal samples in their metabolomic profiles but extended this to describe metabolomic heterogeneity within tumors. High-resolution quantitative profiling of metabolites from normal and tumor tissue identified 2 metabolic subgroups that were associated with aggressive tumor characteristics. Interestingly, although these metabolic subtypes did not strongly correlate with established PAM50 or claudin-low tumor gene expression subtypes. While larger studies are warranted to allow better powered analysis of how intrinsic subtype affects metabolism, an important next step was to understand whether other genomic signals correlate with metabolic subgroups.

We hypothesized that interactions between tumor epithelium and stroma can be detected in genomic data and that these interactions contribute to the evolution of distinct metabolic microenvironments. Previous genomic studies have emphasized the importance of host–tumor interactions in progression (45), suggesting possible pathways to test in association with metabolomic subtype. From microarray studies, a wound response is tumor promoting in the presence of initiated cells (46, 47) and wound response gene signatures are prevalent in tumor-adjacent, stroma-enriched tissue (25). Thus, we tested an in vivo wound response signature and found its expression to be strongly associated with metabolic subtype. These results show that integrated analyses of metabolite and gene expression data can identify phenotypically distinct groups of breast cancers.

Other recent articles have conducted integrated analyses of metabolomics and gene expression datasets, emphasizing identification of genes associated with individual metabolite levels or identification of metabolic phenotypes associated with specific tumor characteristics (48–51). A large study of metabolite profiles across 289 tumor samples identified groups of tumors corresponding to hormone receptor status and grade, however, that study did not seek to identify associated genomic changes beyond these tumor characteristics (49). Borgan and colleagues analyzed metabolite profiles of 46 mostly luminal A breast cancers (48) and linked metabolic heterogeneity within the luminal A breast cancers to gene expression differences, including differences in extracellular matrix (ECM) genes. Their observation of the importance of ECM genes in metabolism is consistent with our finding that metabolism was strongly correlated with stromal signature and stromal–epithelial interactions. Also consistent with our findings, Borgan and colleagues noted heterogeneity within the luminal A subtype, showing that intrinsic subtype alone does not determine metabolic phenotype (48). In our study, we aimed to include a more diverse sample set with 6 different subtypes of breast cancer and an approximately even distribution of ER⁺ and ER⁻ tumors. Our power to conduct comprehensive assessments of subtype-specific changes was limited due to sample size, but our results do suggest heterogeneity of response within classes defined by tumor subtype.

Few integrated, observational studies of gene expression and metabolomics data, including our own dataset, have been of sufficient size to evaluate both tumor and stromal characteristics. To address this limitation, we combined our tissue-based observations with well-established cell culture–based models that mimic the tumor interactions with microenvironment (5) to confirm the pathway changes seen in our gene expression and metabolomics data. While we were unable to measure all metabolites, we selected the metabolism of glucose for more focused investigation because it is central to many of the Warburg-like changes observed in the full metabolic profile. Our results show that metabolic phenotype is a complex interplay between tumor characteristics and the surrounding stromal biology. Basal-like breast cancer epithelial cells more readily take up glucose from surrounding tissue than luminal breast cancer cells. Furthermore, basal-like CAFs also stimulate marked increases in glucose uptake, even in luminal breast cancer cells. Conversely, luminal CAFs produce no change or can slightly decrease glucose uptake and oxidation in basal-like breast cancer cells. Glucose oxidation depended on both stromal and epithelial characteristics, whereas glycogen synthesis appears to be less sensitive to the differences between luminal CAFs and basal-like CAFs, correlating most strongly with epithelial cell content. This latter pattern of metabolite production in coculture suggests an acquired capability—increased glycogen storage in response to microenvironmental signals—and raises the question as to whether increased glycogen synthesis is a hallmark of basal-like breast cancers.

While primary CAFs used in this study show that there is interindividual variation in the effect of CAFs, to establish central trends for basal-like versus luminal stroma will require greater numbers of cell lines. However, these results underscore the importance of understanding variation in tumor stroma. Characteristics of fibroblasts, such as aging or senescence phenotypes, can also alter metabolism as shown in recent studies (52). Our results confirm that understanding which metabolites are sensitive to stromal factors and which are dominated by epithelial characteristics is necessary if metabolic processes are to be targeted for cancer treatment or prevention.

Future metabolic-targeted treatment demands understanding of not just the metabolomics phenotypes but the genomic signals/pathways that drive these phenotypes. We conducted cytokine arrays to identify factors induced in coculture that could explain the shift to greater glucose uptake. Having observed large fold changes in HGF secretion, we conducted a literature search that suggested GLUT1 expression is downstream of HGF-cMET signaling in liver and muscle (53, 54). Therefore, we hypothesized that the same pathway may be active in breast cancer. Previous literature shows a role for HGF/cMET in breast cancer aggressiveness. MET receptor protein tyrosine kinase regulates cell motility and invasion (55, 56) and is stimulated by HGF (57, 58). Together, the HGF-MET pathway has been shown to regulate stromal–epithelial interactions in multiple cancers (59–63). In breast cancer, HGF/c-MET signaling promotes tumorigenesis (64), increases metastasis (65, 66), and mediates drug resistance in most aggressive breast cancers (67–69). Our results show that HGF regulates GLUT1 expression, which in turn regulates glucose uptake. Inhibition of exogenous HGF completely blocks GLUT1 expression and decreases glucose uptake. Given the success of small-molecule inhibitors of c-MET in the clinic, the HGF/c-MET pathway may also be a target for the regulation of tumor cell metabolism.

In summary, reciprocal interactions between cancer epithelial cells and the surrounding microenvironment have an established impact on tumor growth (5, 6, 46) and a broad range of other metabolic and signaling processes (26, 43, 70). Previous studies by Castello-Cros et al. have shown extensive evidence in vitro for the stromal role in Warburg metabolism through matrix remodeling (2), stromal autophagy (11), and stromal–epithelial lactate exchange (71). Our results further show that both fibroblast and epithelial characteristics modulate specific metabolic phenotypes in vitro and pairs one of these changes, glucose uptake, with specific targetable gene expression changes (HGF/GLUT1). This study shows that integration of multiple data types in human tumor studies, together with in vitro experimental studies that dissect heterotypic interactions, can yield important advances in understanding the complex metabolic and genomic interactions during tumor evolution.

No potential conflicts of interest were disclosed.

Conception and design: H.A. Brauer, C.M. Perou, M.A. Troester

Development of methodology: H.A. Brauer

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): L. Makowski, P. Casbas-Hernandez, L.J. Lang, A.J. Freemerman, C.M. Perou

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): H.A. Brauer, L. Makowski, K.A. Hoadley, L.J. Lang, E. Romàn-Pèrez, M. D'Arcy, C.M. Perou, M.A. Troester

Writing, review, and/or revision of the manuscript: H.A. Brauer, L. Makowski, K.A. Hoadley, L.J. Lang, C.M. Perou, M.A. Troester

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): L.J. Lang, M. D'Arcy

Study supervision: M.A. Troester

The project was funded by the National Institutes of Environmental Health Sciences and NCI via Breast Cancer and the Environment Research Program (BCERP) U01-ES019472 and R01-CA138255 and via a Breast SPORE program (P50-CA58223). L. Makowski was supported by UNC University Cancer Research Fund, National Institute of Diabetes, Digestive and Kidney Disease (P30DK056350), National Institute of Alcohol Abuse and Alcoholism (AA017376), and National Institutes of Environmental Health Sciences and NCI (ES019472).

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.