Altered energy metabolism and metabolic gene expression associated with increased metastatic capacity identified in MDA-MB-231 ceII Iine variants

Yan Tu, Cameron N. Johnstone, James G. Ryall, Guillermo H. López-Campos, Christine R.Keenan,0,, Alastair G. Stewart

1Department of Pharmacology and Therapeutics, the University of Melbourne, Parkville, Victoria 3010, Australia.

2Cancer Research Division, Peter MacCallum Cancer Centre, East Melbourne, Victoria 3000, Australia.

3Olivia Newton-John Cancer Research Institute, Heidelberg, Victoria 3084, Australia.

4School of Cancer Medicine, La Trobe University, Bundoora, Victoria 3086, Australia.

5Department of Biochemistry & Molecular Biology, Monash University, Clayton, Victoria 3800, Australia.

6Department of Pathology, the University of Melbourne, Parkville, Victoria 3010, Australia.

7Department of Physiology, the University of Melbourne, Parkville, Victoria 3010, Australia.

8Health and Biomedical Informatics Centre, Melbourne Medical School, the University of Melbourne, Victoria 3010, Australia.

9Current address: The Wellcome-Wolfson Institute for Experimental Medicine, Queen's University, Belfast BT9 7BL, UK.

10Current address: The Walter and Eliza Hall Institute of Medical Research, Parkville, Victoria 3010, Australia.

11Current address: Department of Medical Biology, The University of Melbourne, Parkville, Victoria 3010, Australia.

Abstract Aim: Despite current advances in therapies and the gradual decline in breast cancer-related mortality, metastasis remains a major therapeutic challenge for treatment. Energy reprogramming is now recognized to be an important part of tumorigenic processes, but its relevance in metastatic dissemination has yet to be elucidated.Methods: Using the MDA-MB-231HM.LNm5 cell line, a novel, highly metastatic variant line derived from TN human breast adenocarcinoma MDA-MB-231 line, alteration in growth and energy metabolisms associated with enhanced metastatic potential were described. Glycolysis and oxidative phosphorylation (OXPHOS) was characterized using the seahorse XF analyzer. Whole transcriptome sequencing (RNA-seq) and quantitative real-time PCR was used to ascertain expression differences in metabolic genes.Results: We observed reduced proliferation, and an elevation of both glycolytic and OXPHOS metabolism in the highly metastatic daughter line. The elevated metabolic rate is only partially reflected by transcript levels of relevant metabolic regulators. Heightened mitochondrial respiration is potentially underpinned by increased expression mitochondrial electron transport chain components. However, increased glycolysis was not underpinned by up-regulation of metabolic genes encoding enzymes participating in glycolysis.Conclusion: Our results indicate breast tumour cells with elevated metastatic propensity are more metabolic active. We also identified differentially expressed metabolic genes, such as IDH2, that may play a part in the metastatic process beyond energy reprogramming.

Keywords: Breast cancer, energy reprogramming, cancer metabolism, metastasis, RNA-seq

INTRODUCTION

The majority of breast cancer-related deaths are not caused by the primary tumor itself, but are due to the results of metastasis to vital organs[1]. Although only a small percentage of patients are initially diagnosed with late stage or metastatic breast cancer, the 5-year survival for these patients is 25% compared with 99%for patients diagnosed with localized disease[2]. In addition, current prognostic markers are unable to accurately predict the risk of metastasis development and approximately 30% of patientsfirst diagnosed with earlier-stage breast cancer will eventually develop recurrent metastatic disease[3]. Therefore, despite current advances in therapies and the gradual decline in breast cancer-related mortality[4], the diagnosis and management of metastatic disease remains a major therapeutic challenge for breast cancer treatment.

The dysregulation of cellular energetics is now regarded as one of the hallmarks of cancer[5]. The metabolic phenomenon describing increased glycolytic capacity in cancer cells, known as “the Warburg effect”[6],stimulated decades of research directed towards the characterization of the reprogramming of energy metabolism during cellular transformation and its role in tumor development. The Warburg effect emerged as just one component of global changes in energy metabolism occurring in both cancer cells and the tumor microenviroment[7,8]. Additionally, an increasing number of studies suggest that metabolic reprogramming plays an important role not only in the process of malignant transformation, but also in the growth and survival of tumor cells within a hostile environment, such as the often limited nutrient and oxygen supply in solid tumours[9-11]. However, despite the significant number of studies that investigated the metabolic programming of primary cancer cells, less is known about metabolic alterations in the context of metastatic disease, especially in breast cancer.

Comparison of breast cancer cell lines panel reveals that cell lines with molecular subtypes associated with more aggressive disease progression exhibit an overall increase in energy metabolic processes, including glycolysis and oxidative phosphorylation (OXPHOS)[12-15]. Studies using metastases derived from the same primary tumour reported more puzzling metabolic changes. In a xenograft model using circulating tumor cells from a breast cancer patient, a proteomic comparison between parental cells and cells that metastasized to the brain demonstrated up-regulation in enzymes involved in both glycolysis and mitochondrial respiration pathways[16]. Moreover, compared to primary tumour cells, circulating tumour cells derived from 4T1 mouse mammary tumors exhibited elevated expression of mitochondrial respiration pathway genes, but not glycolytic genes, while lung metastasis from the same primary tumour revealed modest metabolic change[17].Consistent with thisfinding, increased OXPHOS, were observed with increased metastatic potential in several metastatic cell line variants derived from the same primary breast cancer[13,18]. Thesefindings provide evidence that energy reprogramming may be an important feature of the complex process of breast cancer metastasis, but also raise the question of whether the metabolic profiles of metastatic cells vary depending on the stage of metastasis and site of distant metastasis.

To gain a better understanding of the metabolic changes underlying the process of breast cancer metastasis,we characterized a highly metastatic variant line of the commonly used triple-negative human breast ad-enocarcinoma cell line MDA-MB-231[19], and compared cellular and metabolic alterations. The MDA-MB-231HM.LNm5 cell line is a highly angiogenic and metastatic variant of the MDA-MB-231 cell line derived by in vivo passaging whereby spontaneous secondary lesions were isolated and expanded ex vivo[20-23]. We recently demonstrated that the metastatic ability of MDA-MB-231HM.LNm5 line is highly elevated compared to the parental MDA-MB-231 cells. In a metastasis model involving surgical resection of the primary tumour, NSG immune-deficient mice bearing the HM.LNm5 line exhibited primary tumour recurrence, as well as significant lung, liver, spleen, lymph and spine metastasis. By comparison, no metastatic lesions were detected in secondary organs of MDA-MB-231-innoculated mice[23].

In this study, the metabolic profiles of MDA-MB-231HM.LNm5 were compared to the parental MDAMB-231 cell line using the Extracellular Flux (XF) Analyzer thus enabling simultaneous measurement of the two major cellular energy-producing pathways, glycolysis and OXPHOS. We then used whole transcriptome sequencing (RNA-seq) and quantitative real-time PCR (RT-qPCR) to ascertain expression differences in metabolic genes that were associated with enhanced breast cancer metastatic ability.

METHODS

Cell culture

The MDA-MB-231 human breast adenocarcinoma cell line was purchased from ATCC[19]. MDA-MB-231HM cells[20,21]were kindly provided by ZM Shao and ZL Ou (Breast Cancer Institute, Fudan University, Shanghai,China). MDA-MB-231HM.LNm5 cells were derived as described below. All lines were maintained in phenol red - containing RPMI1640 (Invitrogen) supplemented with 10% (v/v) heat-inactivated fetal calf serum (FCS),2 mmol/L L-glutamine, 1% (v/v) non-essential amino acids, 5% (v/v) sodium pyruvate, 100 U/mL penicillin, 15 mmol/L HEPES buffer, and 0.2% (v/v) sodium bicarbonate (Sigma-Aldrich). Cells were maintained at 37 °C in 5% CO2and passaged every 4-5 days.

Generation of reporter gene tagged MDA-MB-231 variants

To facilitate optical imaging of tumors in vivo, both parental MDA-MB-231 and MDA-MB-231HM cells were transduced with retrovirus encoding tdTomatofluorescent protein, selected with Blasticidin S and bulk sorted for tdTomato expression byflow cytometry (FACSAria, Beckton Dickinson), as previously described[24].Both populations were also transduced with retrovirus encoding Firefly luciferase and selected using puromycin[24]. Parental MDA-MB-231 cells were additionally transduced with retrovirus encoding enhanced GFP (encoded by the pFBneoGFP plasmid, a kind gift from Hiroshi Nakagawa, University of Pennsylvania),selected using G418, and bulk sorted for GFP expression usingflow cytometry (FACSAria, Beckton Dickinson). MDA-MB-231HM.LNm5 cells were isolated from a spontaneous axillary lymph node metastasis that developed from a reporter gene tagged MDA-MB-231HM primary inguinal mammary tumour in a BALB/c-SCID mouse.

Cellular proliferation

Cells seeded at 105cells/cm2into 24-well plastic plates were allowed to adhere overnight and were then rendered quiescent by incubation in serum-free medium containing 0.25% (v/v) bovine serum albumin for 24 h before re-exposing to FCS (5%, 10%) for 48 h. Viable cells were identified by the trypan blue (0.06% w/v)exclusion[25]and enumerated (blinded) manually using a haemocytometer chamber.

Cellular proliferation using Resazurin (Alamar Blue)

Cell proliferation was also assessed using the Resazurin dye method by measuring reduction of the redox dye resazurin to resorufin[26]. Cells were seeded and treated as described in previous section, and were then incubated with Resazurin reagent containing 1.5% (w/v) Resazurin, 0.25 (w/v) methylene blue, 2.9% (w/v)potassium hexacyanoferrate (III) and 4.22% (w/v) potassium hexacyanoferrate (II) trihydrate for 2 h at 37 °C in 5% CO2. Resorufin formation was measuredfluorometrically (excitation 570 nm; emission 620 nm) using a FlexStationII (Molecular Devices, Sunnyvale, CA, USA). Results are expressed in relativefluorescent units.

Extracellular flux assay

The extracellular acidification rate (ECAR)and oxygen consumption rate (OCR) were measured in realtime using the XF24 extracellular flux bioanalyzer (Seahorse Bioscience, Agilent Technologies Australia,Mulgrave, Vic, Australia), as described previously[27]. Briefly, MDA-MB-231 or MDA-MB-231HM.LNm5 cells were seeded at a concentration of 100,000 cells/well in RPMI medium the day before the assay. One hour before the start of the metabolic profiling assay the medium was changed to XF Base medium (Seahorse Bioscience) supplemented with sodium pyruvate (1 mmol/L), D-glucose (25 mmol/L) and adjusted to pH7.4.To determine glycolytic parameters, ECAR was measured at baseline and after injection of oligomycin(5 μmol/L) [Supplementary Figure 1A]. To determine respiration parameters, OCR was measured at baseline and after injection of Oligomycin (5 μmol/L), [carbonyl cyanide-4 (trifluoromethoxy) phenylhydrazone(FCCP), 1 μmol/L] and a combination of antimycin A and rotenone (2.5 μmol/L each). Parameters of mitochondrial respiration were measured according to the XF cell Mito Stress test user manual [Supplementary Figure 2A].

RT-qPCR

RNA samples were isolated from 3 or more independent experiments. Total RNA was isolated using TRIzol®reagent (Life technologies). RNA (100 ng) was reverse-transcribed using the High Capacity cDNA Reverse Transcription Kit (Invitrogen, Mulgrave, Vic, Australia). Reactions of 5 µL total volume were performed using a Mastercycler® Pro (Eppendorf, Hamburg, Germany). cDNA (1 ng) was used for real-time PCR using iTaq™ universal SYBR® Green Supermix (Bio-Rad, Gladesville, NSW, Australia) and an ABI Prism 7900HT sequence detection system (Applied Biosystems), as described previously[22]. Gene expression was normalized to 18S ribosomal RNA using the 2-ΔCtmethod[28]. Specificity of the primer sets was confirmed by dissociation curve analysis. Primer sequences are listed in Supplementary Table 1.

RNA-seq library preparation and sequencing

RNA was pooled from 9 individual cultures of independent culture passage of MDA-MB-231HM.LNm5 cell line and the parental MDA-MB-231 cell line. Total RNA was isolated using TRIzol® (as above). RNA-seq libraries were constructed from 500 ng total RNA using NEBNext Ultra RNA library prep kit for Illumina(#E7530) with NEBNext Poly(A)mRNA Magnetic Isolation Module (#E7490). Prior to library preparation,RNA was confirmed to be of high quality (RNA integrity number ＞ 8) by Agilent Bioanalyzer 2100 analysis.Paired end 2 × 50 bp rapid sequencing was performed on an Illumina HiSeq 2500 (Melbourne Translational Genomic Platform, University of Melbourne). Raw data wasfiltered by removing reads with adaptor sequences, reads where the percentage of unknown bases is greater than 10%, and reads considered to be of low quality (where bases with quality ≤ 5 constitute greater than 50% of base reads) to obtain “clean reads”.All subsequent analyses are based on “clean reads” only.

RNA-Seq data analysis

FASTQfiles werefirst analysed using FASTQC software (http://www.bioinformatics.babraham.ac.uk/projects/download.html) before proceeding with an integrated sequence trimming and alignment step against the UCSC hg19 human reference genome downloaded from Illumina’s iGenomes (https://support.illumina.com/sequencing/sequencing_software/igenome.html) using Rsubread package (v 1.20)[29]. Reads that were aligned to annotated coding regions of the genome were counted using the “featureCounts” feature from Rsubread[30]. These counts were subsequently normalized using the trimmed mean of M-value method[31]and transformed into counts per million (CPM) to describe gene expression level. As a single replicate per condition was used, we assigned a biological coefficient of variation of 0.3 to proceed with the pairwise comparison analyses for the detection of differentially expressed genes using EdgeR software[32].

Gene Ontology analysis and gene list extraction

A list of all 1,158 nuclear-encoded mitochondrial genes was obtained from the MitoCarta2.0 human inventory[33]. Genes associated with the processes of glycolysis (canonical glycolysis GO:0061621; glycolytic process GO:0006096; positive regulation of glycolytic process GO:0045821; negative regulation of glycolytic process GO:0045820; regulation of glycolytic process GO:0006110) and tricarboxylic acid (TCA) cycle (GO:0006099)were extracted from the gene ontology consortium website using the AmiGO Gene Ontology browser (http://amigo.geneontology.org/amigo)[34]. Genes contained within the mitochondrial respiratory chain complexes[The Hugo Gene Nomenclature Committee (HGNC) family ID: 639; Complex I GO:006120; complex II GO006121; complex III GO:006122; complex IV GO:006123; complex V GO:006124] as well as mitochondrial respiratory chain complex assembly factors (HGNC family ID:645) were extracted from HGNC data base(http://www.genenames.org) under the “gene family” browser[35]. Genes with CPM ＜ 1 were excluded from analysis.

Statistical analyses

All statistical analyses were conducted using Prism v6.0 software (Graph Pad, San Diego, CA, USA). Results are expressed at mean ± SEM from n independent experiments (performed on separate days on cells from a different passage) and analysed as grouped data. Cell proliferation data are expressed as a percentage of unstimulated control cell number. Two-way ANOVA with repeated measures with Bonferroni’s post hoc test was performed to ascertain statistical significance. For XF analyser profiles and qRT-PCR analysis, signi ficance was determined by paired two-tailed Student’s t-test. P ＜ 0.05 was considered a significantly difference for all analyses.

RESULTS

The MDA-MB-231HM.LNm5 cell line shows slower serum-induced proliferation in vitro than the parental MDA-MB-231 cell line

MDA-MB-231HM.LNm5 cells exhibit lower migratory and invasive capabilities compared to the parental MDA-MB-231 cells, despite enhanced metastatic potential[22]. Here we see decreased MDA-MB-231HM.LNm5 cell proliferation compared to the parental cells, as measured by either cell enumeration or resazurin dye reduction. FCS (5% or 10%) increased the number of both parental MDA-MB-231 and MDA-MB-231HM.LNm5 cells [Figure 1A]. However, the number of cells resulting from 48 h of proliferation was significantly reduced in MDA-MB-231HM.LNm5 cells compared to parental cells [Figure 1A]. The commonly used Resazurin “proliferation” assay demonstrated a similar percentage increase in MDA-MB-231 cell number at 5%and 10% FCS [Figure 1B], whereas the FCS response of MDA-MB-231HM.LNm5 cell was barely detectable.The differences in the outcomes of experiments using the two different methodologies not only illustrate the limitation of Resazurin use for assessment of cell proliferation, but also emphasize that the indirect measurement of cell number using metabolically converted substrates without independent verification is prone to generate incorrect conclusions.

The metastatic line MDA-MB-231HM.LNm5 is more metabolically active than the parental MDAMB-231 cells line

Although the conversion of resazurin to resorufin is widely used as a “proliferation” assay, an estimate of mitochondrial metabolic activity could be extracted by calculating resorufin production per cell, as resazurin undergoes enzymatic reduction in the mitochondria to generate thefluorescent resorufin product[36]. Unstimulated MDA-MB-231HM.LNm5 cells have significantly higher basal mitochondrial activity compared to the parental cells, a difference that was not observed in the presence of serum [Figure 1C]. However, as both cytosolic and microsomal enzymes have the ability to reduce resazurin[37], we sought a more precise method for quantification of the metabolic changes potentially associated with enhanced metastatic phenotype.

The XF bioanalyser facilitates real time measurement of the two major energy-producing pathways in the cell, namely glycolysis and OXPHOS. The ECAR, is a measure of glycolysis, and is determined by the net production and extrusion of protons into the culture medium as a result of the conversion of glucose to pyruvate and subsequently to lactate plus H+. Simultaneously, OXPHOS is measured by calculating the OCR.

Figure 1. Proliferation of the highly metastatic MDA-MB-231HM.LNm5 (MDA-231HM.LNm5) and parental MDA-MB-231_ATCC (MDA-231) cell lines. Cell growth in the presence and absence of fetal calf serum (FCS) (5% or 10%) after 24 h serum starvation was measured by enumeration of viable cells (tryplan blue exclusion) (A) or the resazurin fluorometric method (B). Mitochondrial metabolic activity (C)was also determined [total relative fluorescent units (RFU) divided by the total number of cells]. Data is presented as mean ± SEM, n = 4.Two-way ANOVA with Bonferroni’s post hoc test was used to test for statistical significance ★P ＜ 0.05; ★★P ＜ 0.01; ★★★P ＜ 0.001 vs. MDA-231 control; ^P ＜ 0.05; ^^^P ＜ 0.001 vs. MDA-231HM.LNm5 control. BSA: bovine serum albumin; NS: not significant

Measurement of OCR and ECAR baseline conditions in the absence of glutamine and lipids showed a near doubling of OXPHOS (OCR) and an approximately 25% increase in glycolysis (ECAR) in the MDA-MB-231HM.LNm5 cells compared to the parental MDA-MB-231 cells [Figure 2A]. Both cell lines were challenged to their maximum glycolytic and respiratory capacity by treatment with oligomycin and FCCP, respectively. Oligomycin inhibits ATP production by inhibiting the mitochondrial ATP synthase (complex 5).This subsequently triggers any cellular energy production that was occurring by respiration to shift to glycolysis, thus revealing the maximum glycolytic rate [Supplementary Figure 1A]. FCCP, on the other hand,is an uncoupling agent that disrupts the mitochondrial membrane potential and stimulates the respiratory chain to operate at maximum capacity [Supplementary Figure 2A]. Compared to the parental MDA-MB-231 cell line, MDA-MB-231HM.LNm5 cells showed higher maximum glycolytic and marginally higher respiratory capacity [Figure 2B].

To ensure that the metabolic alteration observed was independent of the exogenousfluorescent proteins and luciferase in the cells the assay was repeated in reporter gene-free MDA-MB-231 and MDA-MB-231HM.LNm5 lines. Basal OCR and ECAR were increased in the MDA-MB-231HM.LNm5 lines compared to the parent, albeit the differences were less striking. Maximum OCR was also higher in the MDA-MB-231HM.LNm5 lines, while maximum ECAR remained similar [Supplementary Figure 3].

Figure 2. The metabolic phenotype of the highly metastatic MDA-MB-231HM.LNm5 (MDA-231HM.LNm5) and parental MDA-MB-231(MDA-231) cell lines. Baseline oxygen consumption rate (OCR) and extracellular acidification rate (ECAR) calculated in the presence of glucose were combined to generate the basal phenogram (A), while maximum ECAR (Oligomycin treated) and maximum OCR[carbonyl cyanide-4 (trifluoromethoxy) phenylhydrazone (FCCP) treated] were combined to generate the maximum phenogram (B).The glycolytic profile shows ECAR readings in the presence of glucose (11 mmol/L) from which the basal glycolytic rate was calculated(C) [Supplementary Figure 1]. Following injection of oligomycin (5 μmol/L), both maximum glycolytic rate (G) and glycolytic reserve(H) could be determined. The respiration profiles of both cell lines (D) show OCR readings in the presence of glucose and following subsequent addition of oligomycin (5 μmol/L), FCCP (1 μmol/L) and finally antimycin plus rotenone (2.5 μmol/L of each). This procedure allows the quantification of basal mitochondrial respiration (E) and ATP production (F), as well as maximal mitochondrial respiration,spare respiration capacity, proton leak and non-mitochondrial respiration [Supplementary Figure 2]. Data are presented as mean ± SEM,n = 5-7. The student’s t-test was used to test for statistical significance. NS: not significant; ★P ＜ 0.05, ★★P ＜ 0.01

MDA-MB-231HM.LNm5 cells display enhanced glycolytic reserve, mitochondrial respiration and ATP synthesis

The increase in ECAR in the presence of oligomycin not only demonstrates the maximum glycolytic rate,but also shows the glycolytic reserve [Supplementary Figure 1A]. MDA-MB-231HM.LNm5 cells showed a larger increase in ECAR compared to the MDA-MB-231 parental line following oligomycin treatment[Figure 2D], revealing higher maximum glycolytic rates and reserves [Figure 2G and H]. Similarly, the difference between maximum OCR and basal OCR allows calculation of the spare respiratory capacity, which did not differ between the two cell lines [Supplementary Figure 2B].

Both mitochondrial and non-mitochondrial respiration contributed to the basal and maximum OCR. The combination of rotenone, a complex I inhibitor, and antimycin A, a complex Ill inhibitor, shut down mitochondrial respiration completely, leaving respiration driven by processes outside the mitochondria only.MDA-MB-231HM.LNm5 cells showed significantly higher mitochondrial-dependent basal respiration[Figure 2E] and similar mitochondrial and non-mitachondrial -dependent maximum respiration rates compared to the parental cells [Supplementary Figure 2C and D].

The two processes that control basal mitochondrial respiration, ATP production and proton leak, can be probed with the blockade of ATP synthase using oligomycin. Measuring the reduction in OCR upon addition of oligomycin revealed significantly higher mitochondrial ATP synthesis in MDA-MB-231HM.LNm5 cells compared to parental MDA-MB-231 cells [Figure 2F], but unchanged proton leak-driven respiration[Supplementary Figure 2E].

Gene expression analysis of energy metabolism pathways

In order to associate the observed metabolic changes with specific genetic or epigenetic alterations, wefirst selected several genes encoding enzymes that participate in glycolysis and the TCA cycle that were documented to contribute to altered metabolism in cancer cells[38]. RT-qPCR analysis of glucose transporter type 1 [solute carrier family 2 member 1 (SLC2A1)], hexokinase 2, fructose-2,6-biphosphatase 3, muscle pyruvate kinase 2, pyruvate dehydrogenase kinase 1, cytosolic isocitrate dehydrogenase-1 (IDH1), succinate dehydrogenase complex subunits C and D and fumarate hydratase mRNA showed similar expression levels between the MDA-MB-231HM.LNm5 cells and parental cells [Figure 3].

To produce an unbiased analysis, the whole transcriptome of each cell line was then deep-sequenced using RNAseq, and the expression of genes involved in key pathways of energy metabolism was compared, including those influencing glycolysis and mitochondrial respiration. Gene expression level was expressed as CPM and the expression level of gene sets was compared by calculating the ratio between two cell lines using the MDA-MB-231 parental line as the denominator. Transcript per million was also compared and yields similar ratio (data not shown). The comparison of all mitochondrial genes showed a symmetrical distribution of expression around a log-fold change of 0, indicating no predominant direction of effect, although some genes were dysregulated between the two cell lines [Figure 4A]. Genes encoding enzymes directly involved in glycolysis were expressed at lower levels in MDA-MB-231HM.LNm5 cells compared to the parental cells [Figure 4B, Supplementary Table 2]. In particular, hexokinase domain containing 1 (HKDC1), encoding the hexokinase isoform HKDC1 which catalyzes the rate-limiting and obligatoryfirst step of glucose metabolism[39],was significantly down-regulated (Table 1 log2FC = -6.64). However, the majority of genes involved in regulating glycolytic processes showed unaltered expression between the two cell lines. The most differentially expressed genes were those that were down-regulated in metastatic cells [Figure 4B, Table 1, Supplementary Table 3], including MLX interacting protein-like (MLXIPL, log2FC = -6.73), encoding a leucine zipper transcription factor of the Myc/Max/Mad superfamily, and FBP1 (log2FC = -5.36), encoding the gluconeogenesis regulatory enzyme fructose-1,6-bisphosphatase-1. Reduced expression of these genes in MDA-MB-231HM.LNm5 was confirmed by RT-qPCR [Figure 5].Expression of TCA cycle genes was similar between two cell lines with the exception of IDH2 (mitochondrial isocitrate dehydrogenase), which was expressed at one- fifth the levels of the parental MDA-MB-231 cells (Figure 4C, log2FC = -2.39, Supplementary Table 4). This down-regulation was confirmed by RT-qPCR [Figure 5].

Figure 3. qRT-PCR analysis of metabolic regulatory genes, glycolysis and tricarboxylic acid (TCA) cycle genes in MDA-MB-231HM.LNm5 (MDA-231HM.LNm5) and parental MDA-MB-231 (MDA-231) cell lines. Expression of hypoxia inducible factor 1 alpha subunit(HIF1a), solute carrier family 2 member 1 (SLC2A1)/GLUT1, hexokinase 2 (HK2), fructose-2,6-biphosphatase 3 (PFKFB3), pyruvate kinase, muscle (PKM2), pyruvate dehydrogenase kinase 1 (PDK1), cytosolic isocitrate dehydrogenase-1 (IDH1), succinate dehydrogenase complex subunits C and D (SDHC, SDHD), and fumarate hydratase (FH) were measured by SYBR-green qRT-PCR relative to 18S rRNA levels. Results are presented as mean ± SEM, n = 7-8. The student’s t-test was used to test for statistical significance. Not significant (NS),compared to MDA-MB-231 cells

The electron transport chain (ETC) in mitochondria is a key site for oxidative phosphorylation and is the major energy source used to produce ATP. The aforementioned XF mitochondrial stress test quantitatively probes this process. Expression levels of allfive complexes were higher in metastatic daughter line compared to parental line, while genes belonging to complex II and III showed the greatest up-regulation. The expression of ubiquinol-cytochrome C reductase complex III chaperone (BCS1L), encoding a ubiquinol-cytochrome C reductase complex III chaperone, was the most strikingly elevated of all the ETC genes (Figure 4D, log2FC= 1.71, Supplementary Table 5).

DISCUSSION

Figure 4. RNA-seq analysis of metabolic gene expression alteration between the MDA-MB-231HM.LNm5 and parental MDA-MB-231 cell lines. Expression level of all mitochondrial genes (MitoCarta 2.0) were compared (A), as well as genes involved in key processes of energy metabolism, including glycolysis (B) (glycolytic process: GO: 0061621 & 0006096; positive regulator of glycolytic process: GO:0045821; negative regulator of glycolytic process: GO: 0045820; regulation of glycolytic process: GO: 0006110), tricarboxylic acid (TCA)cycle (C) (GO: 0006099), and the electron transport chain (D) (mitochondrial respiratory chain complexes: HGNC family ID: 639 &mitochondrial respiratory chain complex assembly factors HGNC family ID: 645). The log2FC (y-axis) is derived from counts per million(CPM) values for MDA-MB-231HM.LNm5 divided by CPM values for MDA-MB-231, where a positive FC value represents up-regulation in the MDA-MB-231HM.LNm5 cells and a negative value represents down-regulation. Genes with a CPM value of ＜ 1 across both samples were not included. FC: fold change; ENO3: enolase 3; HKDC1: hexokinase domain containing 1; MLXIPL: MLX interacting proteinlike; FBP1: fructose-1,6-bisphosphatase-1; PRKAG2: protein kinase, AMP-activated, gamma 2 non-catalytic subunit; IDH2: isocitrate dehydrogenase-2; BCS1L: ubiquinol-cytochrome C reductase complex III chaperone

MDA-MB-231 human breast cancer cells, originally derived from the pleural effusion of a patient with metastatic dissemination[19], exhibit a gene expression signature predicting poor-prognosis[40]. Although this line has been used to study breast cancer metastasis, and despite deriving from tumors with metastatic capability in the original patient, the MDA-MB-231 cell line often displays poor spontaneous metastatic ability when used in immuno-compromised mice, including BALB/c nude and NOD.SCID strains[41]. The MDA-MB-231HM.LNm5 cell line, on the other hand, provides a much better model for the study of breast cancer metastasis in vivo. Compared to the parental MDA-MB-231 cells, this lymph node metastasis-derived sub-line not only shows aggressive spontaneous metastasis, but also mimics the organ tropism of metastatic human breast cancer, with spontaneous metastasis to lung, liver, spleen and sentinel lymph node[22,23].

Figure 5. qRT-PCR analysis of candidate differentially expressed metabolic genes between MDA-MB-231HM.LNm5 (MDA-231HM.LNm5) and parental MDA-MB-231 (MDA-231) cell lines. Results are presented as mean ± SEM, n = 9. Expression in MDA-231 was set to 1.The student’s t-test was used to test for statistical significance. NS: not significant; ★P ＜ 0.05 compared to MDA-MB-231 parental cells

Real time bioenergetics assessment revealed an elevated glycolytic rate and oxidative phosphorylation in MDA-MB-231HM.LNm5 cells compared to the parental line, suggesting that the more metastatic line offers greater energy plasticity. This increased metabolic capacity reflects a composite of both energy demands for energy production used in macromolecule biosynthesis and metabolism and could be a result of an increased energy requirement accompanying the acquisition of metastatic potential. Interestingly, we showed that this enhanced metastatic ability was associated with reduced in vitro migratory and proliferative phenotype[22,23].

Enhanced proliferative rate has long been considered as a hallmark of tumor cells, which is the basis for conventional chemotherapy[5]. Early molecular profiling studies of human breast tumors revealed that increases in proliferative gene signatures (for example genes directly associated with cell cycle progression) were associated with worse clinical outcome[42,43]. However, evidence also shows migratory, and thus invasive phenotype and proliferative phenotype are not expressed simultaneously in breast cancer. Indeed, breast cancer subpopulations with elevated metastatic activity are not more proliferative than their parental population[44].Recentfinding revealed MDA-MB-468 cells with reduced E-cadherin (inducing EMT) were more migratory,invasive and less proliferative[45]. Others showed positive correlation between bone marrow metastasis and the levels of circulating but non-proliferating breast cancer cells[46]. Furthermore, the correlation between breast cancer cell lines extracted from tumours of various disease stages and their growth rate indicate that proliferation decreases with disease progression[47]. These observations, together with our own, support the phenomenon known as the “migration/proliferation dichotomy”[48]or a “go or grow” mechanism[49], where cell motion and proliferation appear to be mutually exclusive phenotypes.

The inverse relationship observed between cell proliferation and metastatic ability may be explained by the cancer stem cells theory, where quiescent/slowly dividing cells exhibit increased tumorigenic potential[50-52].In addition to slow growth rate, these quiescent stem cells are also relatively resistant to current chemotherapy and radiotherapy treatments[53], show increased metastatic ability through the epithelial-to-mesenchymal transition[54]and potentially explain the inter-tumoural heterogeneity and therapeutic failure seen in metastatic breast cancer[55].

Table 1. Differentially expressed metabolic genes in MDA-MB-231HM.LNm5 (MDA-231HM.LNm5) and parental MDA-MB-231(MDA-231) cell lines, as detected by RNA-seq

Speculation can be made on other biological capabilities requiring higher cellular energy that contribute to increased metastatic potential, including the ability to resist cell death (especially in the circulation), induce angiogenesis, and evade immune destruction[5]. Emerging evidence suggests that some key cellular energetics regulators and processes can also be linked to the induction of angiogenesis[56], the triggering of cancer cell death[57], and shaping the immune micro-environment in the tumor stroma[58]. However, the nature of the relationship between these biological processes and cancer metabolism phenotype has been largely unexplored and warrants further study.

Our results show that the increased glycolysis in the MDA-MB-231HM.LNm5 cells was not underpinned by up-regulation of metabolic genes encoding enzymes participating in glycolysis. On the contrary, glycolytic genes were expressed at a comparatively lower level in the metastatic daughter line. Interestingly, reductions in HKDC1 and MLXIPL expression have been reported to be associated with reduced glucose uptake[59,60],although we did not observe any change in expression levels of any of the major glucose transporters such as GLUT1 (SLC2A1). Protein post-translational modification (PTM) is a key mechanism of regulation in signal transduction pathways. Studies have shown that up-regulated glycolysis can be influenced through diverse PTMs including phosphorylation, acetylation, glycosylation and oxidation of glycolytic enzymes as well as other signaling mediators (reviewed[61]). It is not unlikely that the observed elevation of glycolytic activity in the metastatic cells was the result of PTMs and gene expression were lowered as compensating mechanism.Further studies would need to be carried out to investigate whether proteomic changes are correlated with transcriptomic observations.

The XF mitochondrial stress test revealed that the elevated oxidative phosphorylation observed in the metastatic cells is independent of leaky mitochondria and is mainly explained by the enhanced production of ATP. The result further suggests a higher energy demand in the metastatic MDA-MB-231HM.LNm5 line compared to the parental line. Additionally, we found increased expression of allfive complexes of the mitochondrial electron transport chain, which are the mediators of oxidative phosphorylation. Although this elevation was modest in magnitude, it may be sufficient to shift the entire metabolic profile of the cells.

In addition to the XF analyzer, metabolic status could also be measured by a variety of assays such as direct measurements of various metabolic enzymes, substrates, or ATP as surrogates of total energy metabolism.Although these metabolic assays each have their limitations and are mostly single-point measurements, it would have added valuable verification of our XF observation.

IDH2 expression was significantly reduced in MDA-MB-231HM.LNm5 cells while IDH1 levels remain unchanged. Interest in this family of enzymes in relation to cancer biology arose from reports of recurring mutations in IDH1 and IDH2 genes in several cancers including colorectal cancer and gliomas[62]. The functionality of these mutants and their impact on cancer progression has been the focus of many studies. Currently,inhibitors of mutant IDH1 and IDH2 are in Phase I/II clinical trials for both solid and myeloid tumors. In breast cancer, IDH gene mutations are detected at a frequency of less than 5%[63]. Compared to the substantial focus on mutant forms of IDH, little is known about the role of wild-type IDH1, and even less of wildtype IDH2, in cancer progression and metastasis. Hepatocellular carcinoma patients with reduced levels of IDH2 in tumors were at increased risk of metastatic progression and showed worse prognosis[64]. Similarly,in osteosarcoma, IDH2 levels were inversely correlated with pathological grade and metastasis[65]. The suggestion from these correlative observations, that wild-type IDH2 suppresses metastatic processes, is further supported by our data. In addition, ourfindings suggest that the mechanism by which IDH2 may inhibit metastasis is independent of cellular energy pathways.

Our transcriptomicfindings warrant further studies that directly investigate the role of the abovementioned DEGs in metastatic behaviors of breast cancer cells. Knockdown and/or ectopic overexpression of genes of interest found in our study, such as BCS1L or IDH2, in the metastatic MDA-MB-231HN.LNm5 and/or nonmetastatic MDA-MB-231 cells may reveal the relationship between these genes and metastatic phenotypes including metabolic reprogramming. Moreover, related animal experiments involving the manipulation of the expression of these genes of interest would further characterize their contribution in breast cancer growth and progression.

We acknowledge the limitation of having carried out the metabolic and transcriptomic studies in cultured cells. The clinical relevance of human cell line models has been questioned. Indeed, there is not always a linear correlation between in vitro proliferation or motility and spontaneous metastatic capacity in vivo, as other cellular phenotypes, influencing intravasation, extravasation and survival in the circulation (among others) also play a role. However, to determine precise ECAR and OCR measurement in vivo would be technically challenging. Future studies involving metabolic and transcriptomics analysis of tumour cells isolated in situ are required.

In conclusion, until recently, metabolic reprogramming in the context of metastatic dissemination has been largely unexplored in breast cancer. In the present study, a model of spontaneous metastatic breast cancer was used to identify metabolic alterations involved in breast cancer progression. The highly metastatic MDAMB-231HM.LNm5 line displayed higher glycolytic activity and elevated oxidative phosphorylation compared to the parental MDA-MB-231 line, despite reduced proliferative ability. We also showed that this enhanced metabolic rate is only partially reflected by transcript levels of relevant metabolic regulators. Consideration of protein translation, and post-translational modifications, may provide further insight into the molecular alterations underlying the elevated glycolysis and oxidative phosphorylation in cells with higher metastatic capacity. Characterization of the metabolic changes correlated to enhance metastatic potential would deepen knowledge of metastatic mechanisms, and could facilitate the development of new strategies for therapeutic interventions and clinical management of patients with metastatic breast cancer.

DECLARATIONS

Authors’ contributions

Derived the MDA-MB-231HM.LNm5 line, generated the reporter gene tagged MDA-MB-231 and MDA-MB-231HMLNm5 lines: Johnstone CN

Conducted the Seahorse XF assays and contributed to the analysis and interpretation of the data: Ryall JG

Conducted the RNA-seq library preparation: Keenan CR

Analysed the FASTQfiles from RNA sequencing: López-Campos GH

Conducted the majority of the experiments in the study, performed data analysis and interpretation, and drafted the article: Tu Y

Contribution to conception and design of the study: Tu Y, Johnstone CN, Stewart AG

Major contributor to the conception and design of described work, contributed to writing and editing of the manuscript, and will be the guarantor for this article: Stewart AG

Availability of data and materials

The datasets generated during and/or analysed during the current study are available from the corresponding author on reasonable request.

Financial support and sponsorship

This work was supported by the National Health and Medical Research Council of Australia (1059655,1023185).

Conflicts of interest

All authors declared that there are no conflicts of interest.

Ethical approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Copyright

© The Author(s) 2018.