Skip to main content
  • Submit
  • About
    • Editorial Board
    • PNAS Staff
    • FAQ
    • Accessibility Statement
    • Rights and Permissions
    • Site Map
  • Contact
  • Journal Club
  • Subscribe
    • Subscription Rates
    • Subscriptions FAQ
    • Open Access
    • Recommend PNAS to Your Librarian
  • Log in
  • My Cart

Main menu

  • Home
  • Articles
    • Current
    • Special Feature Articles - Most Recent
    • Special Features
    • Colloquia
    • Collected Articles
    • PNAS Classics
    • List of Issues
  • Front Matter
  • News
    • For the Press
    • This Week In PNAS
    • PNAS in the News
  • Podcasts
  • Authors
    • Information for Authors
    • Editorial and Journal Policies
    • Submission Procedures
    • Fees and Licenses
  • Submit
  • About
    • Editorial Board
    • PNAS Staff
    • FAQ
    • Accessibility Statement
    • Rights and Permissions
    • Site Map
  • Contact
  • Journal Club
  • Subscribe
    • Subscription Rates
    • Subscriptions FAQ
    • Open Access
    • Recommend PNAS to Your Librarian

User menu

  • Log in
  • My Cart

Search

  • Advanced search
Home
Home

Advanced Search

  • Home
  • Articles
    • Current
    • Special Feature Articles - Most Recent
    • Special Features
    • Colloquia
    • Collected Articles
    • PNAS Classics
    • List of Issues
  • Front Matter
  • News
    • For the Press
    • This Week In PNAS
    • PNAS in the News
  • Podcasts
  • Authors
    • Information for Authors
    • Editorial and Journal Policies
    • Submission Procedures
    • Fees and Licenses

New Research In

Physical Sciences

Featured Portals

  • Physics
  • Chemistry
  • Sustainability Science

Articles by Topic

  • Applied Mathematics
  • Applied Physical Sciences
  • Astronomy
  • Computer Sciences
  • Earth, Atmospheric, and Planetary Sciences
  • Engineering
  • Environmental Sciences
  • Mathematics
  • Statistics

Social Sciences

Featured Portals

  • Anthropology
  • Sustainability Science

Articles by Topic

  • Economic Sciences
  • Environmental Sciences
  • Political Sciences
  • Psychological and Cognitive Sciences
  • Social Sciences

Biological Sciences

Featured Portals

  • Sustainability Science

Articles by Topic

  • Agricultural Sciences
  • Anthropology
  • Applied Biological Sciences
  • Biochemistry
  • Biophysics and Computational Biology
  • Cell Biology
  • Developmental Biology
  • Ecology
  • Environmental Sciences
  • Evolution
  • Genetics
  • Immunology and Inflammation
  • Medical Sciences
  • Microbiology
  • Neuroscience
  • Pharmacology
  • Physiology
  • Plant Biology
  • Population Biology
  • Psychological and Cognitive Sciences
  • Sustainability Science
  • Systems Biology
Research Article

The harsh microenvironment in early breast cancer selects for a Warburg phenotype

View ORCID ProfileMehdi Damaghi, View ORCID ProfileJeffrey West, Mark Robertson-Tessi, View ORCID ProfileLiping Xu, View ORCID ProfileMeghan C. Ferrall-Fairbanks, Paul A. Stewart, Erez Persi, Brooke L. Fridley, View ORCID ProfilePhilipp M. Altrock, Robert A. Gatenby, Peter A. Sims, View ORCID ProfileAlexander R. A. Anderson, and View ORCID ProfileRobert J. Gillies
PNAS January 19, 2021 118 (3) e2011342118; https://doi.org/10.1073/pnas.2011342118
Mehdi Damaghi
aDepartment of Cancer Physiology, Moffitt Cancer Center and Research Institute, Tampa, FL 33612;
bDepartment of Oncologic Sciences, Morsani College of Medicine, University of South Florida, Tampa, FL 33612;
  • Find this author on Google Scholar
  • Find this author on PubMed
  • Search for this author on this site
  • ORCID record for Mehdi Damaghi
Jeffrey West
cDepartment of Integrated Mathematical Oncology, Moffitt Cancer Center and Research Institute, Tampa, FL 33612;
  • Find this author on Google Scholar
  • Find this author on PubMed
  • Search for this author on this site
  • ORCID record for Jeffrey West
Mark Robertson-Tessi
cDepartment of Integrated Mathematical Oncology, Moffitt Cancer Center and Research Institute, Tampa, FL 33612;
  • Find this author on Google Scholar
  • Find this author on PubMed
  • Search for this author on this site
Liping Xu
aDepartment of Cancer Physiology, Moffitt Cancer Center and Research Institute, Tampa, FL 33612;
  • Find this author on Google Scholar
  • Find this author on PubMed
  • Search for this author on this site
  • ORCID record for Liping Xu
Meghan C. Ferrall-Fairbanks
cDepartment of Integrated Mathematical Oncology, Moffitt Cancer Center and Research Institute, Tampa, FL 33612;
  • Find this author on Google Scholar
  • Find this author on PubMed
  • Search for this author on this site
  • ORCID record for Meghan C. Ferrall-Fairbanks
Paul A. Stewart
dBiostatistics and Bioinformatics Shared Resource, Moffitt Cancer Center and Research Institute, Tampa, FL 33612;
  • Find this author on Google Scholar
  • Find this author on PubMed
  • Search for this author on this site
Erez Persi
eNational Center for Biotechnology Information, National Library of Medicine, NIH, Bethesda, MD 20894;
  • Find this author on Google Scholar
  • Find this author on PubMed
  • Search for this author on this site
Brooke L. Fridley
fDepartment of System Biology, Columbia University Irving Medical Center, New York, NY 10032;
  • Find this author on Google Scholar
  • Find this author on PubMed
  • Search for this author on this site
Philipp M. Altrock
cDepartment of Integrated Mathematical Oncology, Moffitt Cancer Center and Research Institute, Tampa, FL 33612;
  • Find this author on Google Scholar
  • Find this author on PubMed
  • Search for this author on this site
  • ORCID record for Philipp M. Altrock
Robert A. Gatenby
aDepartment of Cancer Physiology, Moffitt Cancer Center and Research Institute, Tampa, FL 33612;
  • Find this author on Google Scholar
  • Find this author on PubMed
  • Search for this author on this site
Peter A. Sims
fDepartment of System Biology, Columbia University Irving Medical Center, New York, NY 10032;
gDepartment of Biochemistry and Molecular Biophysics, Columbia University Irving Medical Center, New York, NY 10032
  • Find this author on Google Scholar
  • Find this author on PubMed
  • Search for this author on this site
Alexander R. A. Anderson
cDepartment of Integrated Mathematical Oncology, Moffitt Cancer Center and Research Institute, Tampa, FL 33612;
  • Find this author on Google Scholar
  • Find this author on PubMed
  • Search for this author on this site
  • ORCID record for Alexander R. A. Anderson
Robert J. Gillies
aDepartment of Cancer Physiology, Moffitt Cancer Center and Research Institute, Tampa, FL 33612;
  • Find this author on Google Scholar
  • Find this author on PubMed
  • Search for this author on this site
  • ORCID record for Robert J. Gillies
  • For correspondence: Robert.Gillies@moffitt.org
  1. Edited by Robert H. Austin, Princeton University, Princeton, NJ, and approved December 9, 2020 (received for review June 17, 2020)

  • Article
  • Figures & SI
  • Info & Metrics
  • PDF
Loading

Significance

Glucose is converted to energy through “fermentation” or “oxidation.” Generally, if oxygen is available, cells will oxidize glucose to CO2 because it is more efficient than fermentation, which produces lactic acid. But Warburg noted that cancers ferment glucose at a “remarkable” rate even if O2 is available! This “Warburg Effect” is still misunderstood because it doesn’t make sense that a cell would ferment glucose when it could get much more energy by oxidizing it. The current paper goes to the heart of this problem by defining the microenvironmental conditions that exist in early cancers that would select for a Warburg Effect. This is important because such cells are much more aggressive and like to lead to cancers that are lethal.

Abstract

The harsh microenvironment of ductal carcinoma in situ (DCIS) exerts strong evolutionary selection pressures on cancer cells. We hypothesize that the poor metabolic conditions near the ductal center foment the emergence of a Warburg Effect (WE) phenotype, wherein cells rapidly ferment glucose to lactic acid, even in normoxia. To test this hypothesis, we subjected low-glycolytic breast cancer cells to different microenvironmental selection pressures using combinations of hypoxia, acidosis, low glucose, and starvation for many months and isolated single clones for metabolic and transcriptomic profiling. The two harshest conditions selected for constitutively expressed WE phenotypes. RNA sequencing analysis of WE clones identified the transcription factor KLF4 as potential inducer of the WE phenotype. In stained DCIS samples, KLF4 expression was enriched in the area with the harshest microenvironmental conditions. We simulated in vivo DCIS phenotypic evolution using a mathematical model calibrated from the in vitro results. The WE phenotype emerged in the poor metabolic conditions near the necrotic core. We propose that harsh microenvironments within DCIS select for a WE phenotype through constitutive transcriptional reprogramming, thus conferring a survival advantage and facilitating further growth and invasion.

  • DCIS
  • Warburg Effect
  • tumor evolution
  • clonal selection
  • adaptation

Ductal carcinomas in situ (DCIS) of the breast are a heterogeneous group of neoplastic lesions confined to the lumens of breast ducts. In early intraductal cancers, hyperplasia forces cells to grow toward the ductal lumens, which moves cells further from their supplying blood vessels that are restricted to the surrounding stroma (Fig. 1A) (1). As a consequence, these cells are significantly nutrient deprived. Hyperplastic tissue in DCIS can be >0.16 mm thick, which is larger than the diffusion distance of oxygen in tissues, rendering the periluminal areas of DCIS hypoxic (2, 3). This lack of oxygen would induce glucose fermentation due to a Pasteur effect, and the resulting production of lactic acid would make the periluminal areas of DCIS profoundly acidic. This has been verified following identification (4) and subsequent validation (5) of membrane-associated Lamp2b as a marker for acid adaptation, which is abundant in the periluminal cells of DCIS. These microenvironmental properties of hypoxia, acidity, and nutrient deprivation exert strong selection pressure on cancer cell survival, and the metabolic adaptations subsequently feed back to the microenvironment, creating a dynamically changing landscape. Over many years in this environment, cells adapt and emerge with flexible, aggressive, and dedifferentiated phenotypes (6).

Fig. 1.
  • Download figure
  • Open in new tab
  • Download powerpoint
Fig. 1.

Early DCIS conditions can select for glycolytic phenotype. (A) Schematic of early and late DCIS progression. (B) Multiplex IHC staining of DCIS patient sample with markers of glycolysis (Glut1 [green]), acid adaptation (LAMP2b [orange]), hypoxia (CA9 [purple]), lactate production (MCT4 [cyan]), vasculature marked (CD138 [red]), and nuclei (DAPI [blue]). (C) Lactate production rate of clones grown out from cells selected under conditions of being selected through multiple rounds of the following conditions: unfed for 1 mo (UF); low glucose (G); low glucose, oxygen, and pH (GOP); low oxygen and pH (OP); and growth in rich media (Control). (D) Seahorse glycolytic rate assay to measure ECAR and OCR following addition of glucose. Basal glycolysis was higher in UF cells (E), but compensatory glycolysis showed no difference between control clones and overall UF clones (F). (G) UF clones have higher WE phenotype (expressed as ECAR/OCR ratio) than control.

The most prominent metabolic hallmark to emerge from DCIS selection is the Warburg Effect (WE) phenotype, which is defined as aerobic glycolysis, where glucose is fermented to lactic acid even in the presence of adequate oxygen, contributing to the acidity of the ductal microenvironment (1). A WE is commonly observed in aggressive cancers (7, 8) and has been exploited clinically with 18fluoro-2-deoxy-d-glucose positron emission tomography scans as a diagnostic marker of tumor stage and is prognostic of cancer outcome (9). Despite its almost ubiquitous expression in cancers, the causes and consequences of a WE remain a mystery. There have been dozens of mechanisms proposed, yet none have been proven. We have previously proposed that these conditions (hypoxia, acidosis, or nutrient deprivation) would select for cells with WE phenotype. In an initial study, cells were selected with periodic hypoxia (16 h 0.2% O2 and 8 h 21% O2 for 50 cycles). Multiple clones were derived from surviving cells, and these were shown to be pan-therapy resistant and had an E-cadherin to N-cadherin switch and a loss of p53, with a moderate increase in aerobic glycolysis that was not sustained (10). In a subsequent study, we adapted cells to growth in acidic conditions, and this selected for a number of important phenotypes, including anchorage-independent growth, yet it did not select for cells with aerobic glycolysis, although cells adapted to acid pH did ferment glucose more rapidly at a low pH compared to nonadapted cells (4, 11, 12). Vogelstein’s group has shown that nutrient deprivation, specifically limiting (0.2 mM) glucose, promoted the outgrowth of pancreatic cancer cells that express mutant k-ras and a WE phenotype in mixed starting cultures, although de novo selection was not shown (13).

Hence, we hypothesize that, if conditions in DCIS select for constitutive aerobic glycolysis, it may involve a complex and dynamic interplay between the multiple factors of hypoxia, acidosis, and nutrient deprivation. To test this, we subjected low-glycolytic breast cancer cells to a series of these combined selection pressures over a period of many months. At endpoints, individual clones were isolated and characterized for their metabolic and transcriptomic profiles. The resulting selected clones were enriched for populations that constitutively expressed an aerobic –glycolytic (WE) phenotype. Transcriptomic analyses identified a number of relevant factors that could account for constitutive glycolysis, including SP1/KLF4 and NFκB. KLF4 expression was validated on selected clones using Western blots and immunocytochemistry (ICC). Tissue microarray (TMA) and whole mount staining of DCIS patients showed increased expression of KLF4 in DCIS samples, when compared to adjacent normal, as well as a relatively elevated expression at the core of each DCIS where the most selective environment exists. NFκB was also validated to be overexpressed at the protein level in vitro, with a significant increase in the transcriptionally active phosphorylated form, which was associated with increased HK2 expression. Knockdown of NFκB-related p65 reversed the WE in highly glycolytic clones.

We further investigated the emergence of WE phenotypes in areas under harsh selective pressures by adapting a previously published mathematical model of tumor metabolism and growth, informed by empirical data (14, 15). The model simulates a tumor growing in a homeostatic tissue, initialized within a ductal structure with diffusion-limited nutrients. Different tumor phenotypes were allowed to evolve due to selection, and multiple simulations showed that the selection of a WE phenotype occurred in the harshest conditions near the periluminal necrotic core. The model was calibrated to the in vitro results presented herein, and simulations under different conditions suggested that different modes of selection can be in action, depending on cellular turnover and the specific microenvironmental conditions. In particular, the harsh conditions had bottleneck-like selection events, whereas the less-harsh conditions tended to show phenotypic drift.

Thus, we conclude that the microenvironmental conditions existing in DCIS are sufficient, with time, to select for cells with a WE phenotype. In this particular case, the switch to a WE phenotype is related to KLF4 as a phenotypic switch and/or NFκB expression as a survival strategy. This study unravels the role of harsh microenvironmental selection pressures in driving activation of pathways, controlled by key transcription factors, that lead to the WE phenotype and subsequent cancer progression.

Results

Harsh Microenvironments, Similar to Early DCIS Conditions, Select for Clones with Higher Aerobic Lactate Production Rate.

In early carcinogenesis, intraductal hyperplasia leads to significant alterations in the physical microenvironment, especially in periluminal cells that are far (>0.16 mm) from their blood supplies, leading to a highly selective microenvironment of hypoxia, acidosis, and severe nutrient deprivation (1, 6, 16). This suggests that the periluminal cells should be oxygen deprived, which is consistent with increased expression of hypoxia-inducible factor client proteins, such as CA9 and Glut1, in periluminal areas of late-stage DCIS (17, 18). As proof of principle, we performed multiplexed immunohistochemistry (mIHC) on our DCIS stage patient whole mount samples for markers of high glycolysis (Glut1), hypoxia-induced acid production (CA9), non-hypoxia-induced acid production (MCT4), and acid resistance (LAMP2b) (Fig. 1B). Our results illustrate that all these conditions exist inside the DCIS ducts individually or in combination. To better understand the impact of these conditions in DCIS breast cancer and their correlation to the WE phenotype, we subjected low-glycolytic Michigan Cancer Foundation-7 (MCF7) cells to a range of selection forces such as acidity (pH 6.7), hypoxia (1% O2), low glucose (0.1 mM), and combinations thereof, reflecting increasing levels of stress (i.e., low glucose (G); low oxygen and pH [OP]; and low glucose, oxygen, and pH [GOP]). Additionally, as an extreme condition, we selected cells by placing them in a flask and not replenishing the media for 4 wk (unfed [UF]) which caused >95% of the cells to die (Fig. 1C). We excluded acidosis and hypoxia alone as selection pressures, because previous results showed that these conditions alone do not strongly select for a WE phenotype (4, 10). Each of these harsh microenvironments resulted in significant cell death, followed by regrowth under rich microenvironmental conditions (neutral pH, 21% O2, and 5.8 mM glucose). This process was repeated multiple (two to six) times with flasks regaining confluence, typically within 4 wk, before reexposure to harsh conditions. After the final outgrowth, we isolated individual clones (>20 per condition) both from controls that were continuously grown in a rich microenvironment and those that were selected to survive in harsh microenvironments (G, OP, GOP, and UF) by seeding individual cells in 96-well plates, which were then regrown under rich microenvironmental conditions. These clones were then expanded in individual T25 flasks, which were then harvested for freezing and for metabolic profiling for rates of lactate production and glucose consumption under normal culture conditions, as the first sign of WE phenotype, using colorimetric kits. Fig. 1C shows the lactate production rates (LPRs, in nanomole per minute per milligram of protein) for individual clones from the four harsh and one rich (control) conditions. These data demonstrate that harsh environmental conditions preferentially select for clones with increased rates of aerobic lactate production; specifically, the UF and low-GOP conditions had the greatest number of high LPR clones (Fig. 1C). To relate our finding of high LPR to the WE, we further measured lactate production and glucose consumption rates at the same time using a multianalyte system (YSI 2900, Yellow Springs, OH) in 96-well plates. For these studies, we used the three UF clones (UF1, UF9, and UF18) with highest lactate production rates and three clones from the rich-microenvironment MCF7 cells with low LPR. Results shown in SI Appendix, Fig. S1 confirmed the higher LPR, observed by colorimetric assays in the harsh- compared to rich- microenvironment clones.

To confirm the WE phenotype of the harsh microenvironmentally selected clones, compared to parental MCF7 cells, we performed the Seahorse glucose stress test (GST) assay to measure both basal and maximal glycolytic capacity of cells as well as their respiratory capacity (see Methods for more details) (Fig. 1D). We found all theUF clones had higher basal glycolysis rate compared to control clones (Fig. 1E), although their compensatory glycolysis was not different in general (Fig. 1F). However, compared to their parental MCF7 clones, all UF clones showed an increase in the ratio of extracellular acidification to oxygen consumption (ECAR/OCR), which is a measure of the WE phenotype (Fig. 1G). These results indicated that the harsh microenvironmental conditions similar to those found in early DCIS select for a WE (aerobic glycolytic) phenotype; more specifically, the combination of low glucose, low oxygen, and low pH or starvation provide the greatest selective pressure for a WE phenotype.

Clonal Evolution under DCIS Microenvironmental Conditions.

These data have demonstrated that harsh microenvironmental conditions selected for cells with increased rates of aerobic glycolysis. Furthermore, and slightly counterintuitively, these selected cells maintained their WE phenotype even after being placed in abundant nutrients and oxygen conditions for multiple generations (i.e., exceeding 20 passages). To investigate the mechanisms leading to this stable (“hardwired”) phenotypic switch, we performed RNA sequencing (RNA-seq) analyses of our harsh and rich microenvironment selected single-cell–derived clones (see Methods). Briefly, the harsh (UF, GOP, OP, and G) and rich clones were plated in 96-well plates and grown to confluence, from which RNA was extracted and sequenced using pooled library amplification for transcriptome expression (PLATE-seq) (19) (see Methods). After filtering, 12,568 genes were used for further analysis. Unsupervised clustering of the RNA-seq data identified five distinct groups that corresponded to each of the microenvironmental conditions (Fig. 2A and SI Appendix, Fig. S2 A and B). Principal component analysis showed generally good segmentation for the different microenvironmental conditions (Fig. 2B and SI Appendix, Fig. S3). It is notable that the UF cluster was readily segmented from the rest of the cells, suggesting that this condition, which more accurately reflects the in vivo situation, adds selection pressures beyond those imposed by the metabolic selections of G, OP, and GOP. Furthermore, there was some overlap between the parental unselected and some of the selected (G and OP) clones, suggesting some clonal heterogeneity in the parental population or original phenotype recovery due to the highly plastic nature of MCF7 cells (15).

Fig. 2.
  • Download figure
  • Open in new tab
  • Download powerpoint
Fig. 2.

RNA sequencing analysis of selected clones reveals the molecular mechanism of switch to WE phenotype. (A) Heatmap showing the top 500 most variable genes, grouped by selection condition. A preliminary analysis of RNA-seq data was performed by linearly regressing gene expression data with lactate production rate and filtered for significantly correlated or anti-correlated genes. Unsupervised clustering (SI Appendix, Fig. S3) of these data showed that the 1,000 most highly correlated and anticorrelated genes clustered within selection condition. (B) Principal component analysis of gene expression data showed separation of the UF and GOP groups compared to control. (C) Phenotype evolution trajectory alignment of single-clone RNA-seq for evolving breast cancer cell populations. Cell fate analysis with Palantir was applied to the single-clone RNA-seq dataset to determine differentiation potential from an initial, unselected, parental lineage to selected, phenotypic terminal states of G, OP, and UF. UMAP projections were used to visualize the high-dimensional dataset and known identity of each clones was colored on the UMAP projections. Unselected clones were indicated in red, UF clones were indicated in purple, G clones were indicated in green, OP clones were indicated in mint, and GOP clones were indicated in blue.

To determine which genes were associated with the WE phenotype, the gene expression data were linearly regressed against LPR using the limma and voom R packages (20, 21). A total of 676 genes had a significant association with LPR (Padj < 0.1), with 388 having a positive association (correlation coefficient > 0) and 288 having a negative association (correlation coefficient < 0) (Fig. 2A and Dataset S1).

To study phenotypic evolution under the different microenvironmental selection pressures on control MCF7 cells, Palantir analysis (SI Appendix, Supplementary Data1) was applied to the single-clone RNA-seq dataset of unselected (parental), UF, G, OP, and GOP clones to detect evolutionary trajectories of these clones and alignment along pseudotime (Fig. 2C and SI Appendix, Figs. S4 and S5). All clones started in the rich microenvironment of the parental phenotype, indicated by orange points in Fig. 2C. Aligning clones along pseudotime revealed three distinct terminal phenotypes: UF, G, and OP. Interestingly, the GOP phenotype lay along the trajectory of the UF terminus. Analysis of the differentiation potential of the clones showed that those aligned to the earliest pseudotime (dark blue in SI Appendix, Fig. S4) also had the highest differentiation potential, indicating that they are most likely to evolve to one of the terminal phenotypes over pseudotime. Likewise, especially for the UF and G phenotypes, pseudotime was near 1 (yellow in SI Appendix, Fig. S4), indicating that these clones had the lowest differentiation potential and that they were at their terminal phenotypic states.

Mathematical Modeling Shows the WE Phenotype Is Rapidly Selected in Harsh Conditions.

To investigate the emergence of the WE phenotype in more detail, we extended a mathematical model of tumor metabolism (14, 15) to simulate the experiments presented herein. The extensions to the model that simulate the in vitro portions of this work are provided in the Methods, Eq. 1, and Tables 1 and 2. First, we applied our previously published model to simulate DCIS development in vivo. These results indicated that the WE phenotype primarily emerged from the most metabolically depleted area of a simulation, namely far from blood vessels and adjacent to the necrotic core. Fig. 3 A–D shows representative examples of these simulations (refer to SI Appendix, Fig. S6 for the schematic of the model and SI Appendix, Fig. S7 for more information about the simulated barcoding data in in Fig. 3B). The WE phenotype is pink, and after 100 time increments (“days”), this phenotype began to emerge in the center of the duct where glucose and oxygen were highly depleted and the pH was acidic. This suggests that the harsh heterogeneous conditions of a tumor growing within a duct (or other similarly poorly vascularized region) would select for WE phenotypes.

View this table:
  • View inline
  • View popup
Table 1.

In vivo parameters

View this table:
  • View inline
  • View popup
Table 2.

In vitro parameters

Fig. 3.
  • Download figure
  • Open in new tab
  • Download powerpoint
Fig. 3.

Computational modeling of the emergence of the WE phenotype. (A–D) A previously published two-dimensional hybrid discrete-continuum homeostatic cancer metabolism model (42, 43) shows the evolution of acid resistance (blue to green) and WE (blue to pink) phenotypes over time. The model simulates growth into ductal structure (A) where increased acidity in the center of the duct promotes acid resistance phenotypes (blue). After depletion of glucose, the WE phenotype emerges in harsh conditions near the center of the duct, on the edge of the necrotic core. (B) A Muller plot shows phenotypic selection and lineage over time. Vertical axis indicates size of clones, colored by its acid-resistance/WE phenotype shown in C. Final distributions of oxygen, glucose, and acid are shown in D. (E–I) An in vitro version of the model simulated for identical conditions as Fig. 2A confirms that WE phenotype (E) emerges in harsh conditions (“unfed”). Furthermore, these cells have enhanced efficiency in producing ATP from nutrients (F). Model simulated barcode proportions are shown for three timepoints: 0 d (G), 60 d (H), and 120 d (I). Barcodes are colored by average final phenotype with dead clones colored in black. Control and glucose-depleted conditions have low turnover, leading to slowed evolution. OP and UF conditions have increased turnover, selecting for WE phenotypes. Parameters are as follows (Eq. 1): S = 0.08, ko = 0.005, kg = 0.3, kg0 = 2.5, and V0 = 0.93.

We then calibrated the model to directly simulate our in vitro experiments and, the results of this simulation are shown in Movie S1. In short, to recapitulate an in vitro environment, blood vessels were removed from the simulation and nutrient concentrations and pH levels were changed globally across the “flask.” Cells were plated at the same seeding density as in the experiments and were allowed to adapt to their particular conditions through phenotypic drift as in the original in vivo model. Fig. 3 E–I shows that the WE phenotype was selected primarily in the “UF” case, when glucose and pH levels dropped significantly after 14 d without media change. The harsh conditions toward the end of the UF period induced rapid turnover that enabled faster phenotypic drift. Furthermore, these cells exhibited increased adenosine triphosphate (ATP) efficiency (Fig. 3 E and F), which is useful for survival in low-glucose conditions (SI Appendix, Fig. S8).

The barcoding plots in Fig. 3 G–I show how phenotypic selection changes through the period of the simulation for each of the five conditions. The colors correspond to the phenotypes of Fig. 3C, with pink cells having a WE phenotype. The UF condition (bottom panel) showed rapid turnover of the population, driven by an early bottleneck, which quickly drove adaptation to the WE phenotype due to the severely depleted glucose and the highly acidic microenvironment. In contrast, the other conditions showed less turnover, even though some had low glucose and/or acidosis. The adaptation was slower and was primarily aimed at mitigating acidosis rather than becoming glycolytic. Notably, UF cells with a WE phenotype emerged from only a few of the initial cells. This is in contrast to our metabolic profiling of the unselected clones, which showed that one clone had a slightly elevated lactate production rate, alluding to the likely preexistence of this phenotype (Fig. 1C).

The Role of Transcription Factors in Selection of WE Phenotype.

To investigate whether specific transcription factors were involved in the transcriptional switch for generating a WE phenotype, we used the list of 388 significantly and positively associated genes for enrichment analysis by oPOSSUM and Enrichr (22, 23). oPOSSUM “single site analysis” with “human” was selected with default parameters. The top oPOSSUM hit was KLF4 (Z-score = 54.053) (Fig. 4A). As a test of the oPOSSUM analysis, we investigated whether cells from UF clones had high nuclear KLF4 expression using ICC. Fig. 4B shows nuclear localization and higher expression of KLF4 in UF clones compared to their parental MCF7 cells. KLF4 plays a role in early development and promotes a stem-like phenotype (24⇓–26). Consistent with this, we also observed increased RNA expression of genes associated with stemness (SI Appendix, Table S1).

Fig. 4.
  • Download figure
  • Open in new tab
  • Download powerpoint
Fig. 4.

Clinical validation of KLF4 expression in breast cancer patients’ samples. (A) Enrichment analysis of 388 genes positively associated with LPR using oPOSSUM revealed KLF4 as the top regulator of LPR genes. (B and C) ICC analysis of UF and parental MCF7 and DCIS validates the higher expression and nuclear localization of KLF4 in UF cells. (D) Harsh condition in DCIS selects for aggressive cells that can invade other organs. One million of MCF7 parental, UF9, and UF18 clones were injected to tail vein of NSG mice and looked for possible metastasis. The figure shows the metastasis in the lung of the UF group. (E) There was only one metastasis in the control group of 12 mice compared to 100% metastasis in both UF cells. (F) TMA analysis of 204 breast biopsies of Moffitt Cancer Center patients for KLF4 expression shows higher expression of this protein in DCIS compared to adjacent normal tissues. The expression of KLF4 stays high with higher grade of breast cancer. (G) Representative images of the TMA analysis done in A. The box is zoomed in the center of the duct to show spatial correlation of KLF4 in DCIS samples. Cells with high KLF4 expression are located to the center of the duct in DCIS samples that proves our hypotheses that harsh condition selects the cells with glycolytic phenotype through a transcriptional factor switch.

We also queried the Enrichr “Genome Browser Position Weight Matrices” pathway, which contains a list of genes and associated binding motifs from transcription factors, and identified an NFκB-associated list as significantly enriched (Padj = 0.04) (Fig. 4C and SI Appendix, Fig. S9). To investigate this, we performed Western blotting on UF clones and parental MCF7. SI Appendix, Fig. S10 shows that the total amount of NFκB was slightly, but significantly, higher in the UF clones relative to control clones (C clones). In contrast, Western blot of phospho-NFκB (p-NFκB) shows that the p-NFκB, which promotes nuclear localization, was significantly elevated in UF clones compared to that present in the C clones (SI Appendix, Fig. S10A). To further investigate this, we performed ICC on UF clones and parental control MCF7 and found higher expression of p-NFκB in nuclei of UF cells (Fig. 4D). To investigate the role of NFκB in promoting glycolysis in these systems, three separate anti-p65 small interfering RNAs (siRNAs) were prepared and tested for efficacy against UF18 cells and were all able to effect knockdown at different doses. Using the p65siRNA-C, we were able to optimally knockdown expression in UF18 and UF9 clones (SI Appendix, Fig. S10B). In UF18 cells, knockdown of p65 significantly reduced aerobic glycolysis (SI Appendix, Fig. S10C) to the levels observed in nonselected cells. In addition to NFκB, other metabolically relevant proteins that were observed to be different between UF and C clones were pAkt and the NFκB client, HK2 (SI Appendix, Fig. S10A), which were also consistent with increased aerobic glycolysis in the UF cells. These results propose the prosurvival activity of NFκB in our UF cells, which has also been shown for acid-exposed cells in sarcoma (27).

To determine which of the cells expressing a WE were more aggressive compared to nonselected parental clones, we injected cells from two UF clones and one unselected clone into the tail vein of NOD scid gamma (NSG) mice and observed that 100% of both UF cell groups formed metastases, while only one mouse (8.3%) showed metastasis with MCF7 parental cells (Fig. 4E).

KLF4 expression in clinical breast tumor samples.

To further investigate KLF4 expression following selection of WE phenotype in DCIS lesions, we then interrogated expression of KLF4 in breast tumor samples from the Moffitt Cancer Center TMA collection. KLF4 was selected for further clinical validation over NFκB based on its nature of switching phenotype in stem cells and early embryonic cells (25) as well as its previously shown role in regulating glycolysis in stem cells and pancreatic cancer (28, 29). Breast TMAs from Moffitt breast cancer patients were used to study the expression level of KLF4 at different stages. TMA4 contains 204 biopsy cores including adjacent normal samples, DCIS, invasive ductal carcinoma (IDC) with no metastasis, IDC with local metastasis, and lymph node metastasis core biopsy samples. Staining of TMA4 with KLF4 antibody showed significantly increased expression of KLF4 in DCIS samples compared to normal breast tissue. The KLF4 expression remained high in IDCs and metastasis samples (Fig. 4 F and G and SI Appendix, Fig. S11). To relate the role of KLF4 to selection of the WE phenotype in DCIS, we performed spatial analysis of KLF4 in our DCIS samples and observed that multiple sites had very high expression of KLF4 at the center of the duct, where access to nutrient resources is severely restricted, and thus exerts a strong selection pressure with regard to increased acidosis and decreased oxygen (Fig. 4G).

Discussion

The WE phenotype is associated with progression and aggressiveness of cancers and is defined by a high glycolytic rate in the absence or presence of oxygen (aerobic glycolysis). Most cancer cells reprogram their metabolism in favor of aerobic glycolysis despite the presence of plentiful oxygen in their microenvironment. This observation was first reported by Otto Warburg and is thus referred to as the “Warburg Effect” (30⇓⇓–33). We and others have observed this high glycolysis rate in tumors using positron emission tomography (34). We also know that most cancer cells in hypoxic environments (Pasteur effect) compensate for the low ATP yield of glycolysis by overexpression of glucose transporters, such as Glut1 (35). The driving theory for why the WE takes place in cancer is that the high rate of glycolysis benefits cancer cells by increasing ATP production. It also provides many intermediates that are used in subsidiary metabolic pathways for de novo synthesis of nucleotides, amino acids, lipids, and NADPH that are required for cancer cell survival and proliferation. However, none of these cellular regulations individually are enough to hardwire the WE phenotype in cells because they can be altered based on microenvironmental conditions. At the center of individual DCIS, the harsh microenvironment consists of low glucose, low oxygen, and high acidity. Therefore, we hypothesized that there are some biological controls or switches at the genome, transcriptome, or epigenome level that initiate and control the WE phenotype.

Previous research has shown evidence of driver mutations, such as p53 and KRAS, up-regulating the WE phenotype in different cancer types; however, none of these mutation-driven WE phenotypes are consistently present in all cancer cells. Furthermore, there are cancer cells with a WE and no known driver gene mutations. This suggests that there may be mutation-independent drivers of this phenomenon (i.e., the microenvironment). To test our hypothesis, we probed the transcriptome of single selected clones under different microenvironmental conditions recapitulating the environments found in DCIS. Using single clones over single cells had the benefit of allowing us to measure the derived diversity and heterogeneity of a single cell’s progeny over time. Surprisingly, we found a highly variable transcriptome among the clones across all of the selection conditions, which may have been lost at the single-cell resolution. Using sophisticated transcriptome analysis of oPOSSUM and Enrichr, we discovered the transcription factor KLF4 controls all of the LPR genes. KLF4 was previously identified as one of the essential factors for induced pluripotent stem cell development (25). KLF4 was previously reported to regulate WE phenotype (28, 29, 36), although none of these studies connected the KLF4 expression or activation to microenvironmental conditions as evolutionary selection pressures. Here, we have shown that KLF4-induced WE is connected to the microenvironment of cancer cells in DCIS lesions. Open questions still remain regarding the heterogeneity of KLF4 expression in selected clones as well as clinical samples. This might imply redundant mechanisms, such as NFκB, that we also uncovered as a mechanism to maintain the WE phenotype or co-opt adjacent cancer cells (37).

The emergence of a glycolytic phenotype following multiple rounds of harsh selection may be a general phenomenon. For example, selection of immortalized breast epithelial cells by starvation does result in emergence of clones with a WE in two different cell lines. Both of these are associated with increased nuclear KLF4. However, it is entirely possible that other mechanisms exist to generate a WE phenotype and that these may be observed in other cellular systems.

While selected clones clearly can survive and thrive under harsh conditions, it may seem paradoxical that the selected clones can outcompete the nonselected under nutrient-rich conditions, where their growth rates are identical. As the selected clones with a WE are more metastatic, we expect that their greater fitness is derived from their ability to migrate and invade and not necessarily out-proliferate the nonselected clones, although all the clones have the same growth rate as parental ones and some have higher. In the context of avascular DCIS, the first step in our models (1) is the loss of a requirement for basement membrane (BM) attachment, which forces these cells via hyperplasia toward the perilumen, where they are subject to harsh environmental selection. Although those near the BM do have access to a richer environment, they are not under selection and therefore do not adopt an aggressive WE phenotype. Notably, as we have shown in prior works, once cells are selected for a WE, they tend to relocate to the invasive edge (7, 38). This has also been observed in our extended mathematical models that indicated that the WE phenotype primarily evolved in the most metabolically depleted area, namely far from blood vessels and adjacent to the necrotic core. SI Appendix, Fig. S12 shows the evolutionary trajectory of cells in the presence and absence of vascular renewal, clearly implicating emergence of an invasive, aggressive phenotype under the harshest microenvironmental conditions. While mutagenicity of the harsh environmental conditions is likely (6), our current model is based on the WE phenotype emerging from preexisting clones with a fast turnover in the harsh microenvironment with some phenotypic drift. Indeed, metabolic profiling of unselected clones shows at least one with a relatively high preexisting glycolytic rate (Fig. 3G). This is nonetheless an unsettled question and will be a major focus of future research in this area.

There is growing evidence that the stromal microenvironment of preinvasive breast cancer actively participates in DCIS–IDC transition (39, 40) and “reactive” stroma can often be observed in DCIS before the emergence of microinvasion, where progression to invasion can be promoted by cancer-associated fibroblasts (41⇓–43). Furthermore, macrophage infiltration in the periductal stroma surrounding DCIS is also correlated with early dissemination and recurrence (44, 45). Relevant to the current work, it is also possible that the metabolic sequelae (acidity, hypoxia) emanating from the DCIS may induce stromal reactivity (13).

Finally, these results are paradoxical to our notion that cells under very low nutrient conditions should reduce their demand and energy expenditures based on the energy availability. Our findings suggest that the WE phenotype may be more efficient than previously assumed since we show that the WE phenotype is a highly regulated and controlled energy consumption source. Our results also illuminate the evolutionary trajectory of the WE phenotype driven by microenvironment selection pressures. We observed that transcription factors can activate the WE phenotype under appropriate environmental conditions that can both select for the WE phenotype and facilitate its hardwired statues. The activation of transcription factors such as KLF4 and NFκB may serve as an adaptive mechanism for cancer cells to switch to fitter phenotypes (WE) that can withstand the harsh environmental selective forces found in early DCIS lesions.

Methods

Cell Culture and In-Vitro Clonal Selections.

MCF7 cells were acquired from American Type Culture Collection (ATCC, Manassas, VA, 2007 to 2010) and were maintained in Roswell Park Memorial Institute 1640 (Life Technologies, catalog number 11875-093) supplemented with 10% fetal bovine serum (HyClone Laboratories). Growth medium was further supplemented with 25 mmol/L−1 each of Pipes and Hepes and the pH adjusted to 7.4 or 6.7. Cells were tested for mycoplasma contamination and authenticated using short tandem repeat DNA typing according to ATCC’s guidelines.

Western Blotting.

Selected and nonselected MCF7 cells were grown with the same number of passages and used for whole-protein extraction. Lysates were collected radioimmunoprecipitation assay buffer containing 1× protease inhibitor mixture (P8340; Sigma-Aldrich). A total of 20 mg of protein per sample were loaded on polyacrylamide–sodium dodecyl sulfate gels, which later were electrophoretically transferred to nitrocellulose membranes. Membranes were incubated with primary antibodies against rabbit monoclonal KLF4 (1:1,000, ab215036 Abcam), NFκB (1:1,000, #8242 Cell), SignalingHK (1:1,000, #2867s Cell Signaling), p-AKT (1:1,000, #4060s Cell Signaling), Tubulin (1:1,000, #3873 Cell Signaling), and GAPDH (1:4,000, anti-rabbit; Santa Cruz Biotechnology).

siRNA Transfection.

Three unique 27-mer RELA human siRNA oligo duplexes (SR304030A, B, and C) were obtained from Origene (SR304030). Universal scrambled negative control siRNA duplex (SR30004, Origene) was used as a nontargeting control for this study. Cells were seeded in a 6-well plate and reached 70 to 80% confluence before transfection. Cells were transfected with the negative control siRNA or p65-targeting siRNA according to standard protocols using lipofectamine RNAiMAX transfection reagent (13778030, Themo Fisher).

Immunofluorescence.

Selected and nonselected MCF7 cells cultured with the same number of passages were rinsed with phosphate buffered saline (PBS), fixed in cold Methanol:Acetone (1:1) for 10 min and further permeabilized by 0.5% triton 100 and then blocked with 5% bovine serum albumin in PBS. Samples were incubated with KLF4 rabbit monoclonal primary antibody (1:100; ab 215036 Abcam) and secondary Alexa-Fluor 488 anti-rabbit (1:1,000) antibody. Coverslips were mounted using ProLong Gold Antifade Reagent (Life Technologies) and images were captured with a Leica TCS SP5 (Leica) confocal microscope.

Glycolytic and Oxygen Consumption Rate Measurements (Seahorse).

Glycolytic rate of MCF7 and selected MCF7 cancer cells was measured using Seahorse XF96 extracellular flux analyzer and a glycolysis rate kit (Seahorse Biosciences). OCR and ECAR of cancer cells were determined by seeding them on XF96 microplates in their growth medium until they reached over 90% confluence. In these studies, seeding started with 20,000 cells (80% of well area). Measurements were determined 24 h later when the cells reached the 90% confluence. One hour before the Seahorse measurements, culture media were removed and cells were washed three times with PBS followed by addition of base medium (nonbuffered Dulbecco's Modified Eagle Medium supplemented with 25 mM glucose) or our nonbuffered only-glucose-containing solution. For glycolytic rate measurements, mitochondria inhibitors including rotenone (1 μM) and antimycin A (1 μM) were injected after basal measurements of ECAR and OCR of the cells under treatment to stop the mitochondrial acidification. 2-deoxy-glucose (100 mM) was added next to bring down glycolysis to basal levels. Finally, data were normalized for total protein content of each well using the Bradford protein assay (Thermofisher). Seahorse measurements were performed with four to six technical replicates, and these experiments were repeated four times.

Solutions for seahorse experiments.

A total of 2 mM Hepes, 2 mM MES, 5.3 mM KCl, 5.6 mM NaPhosphate, 11 mM glucose, 133 mM NaCl, 0.4 mM MgCl2, and 0.42 mM CaCl2 were titrated to the given pH with NaOH. For reduced Cl− experiments, 133 mM NaCl was replaced with 133 NaGluconate, and MgCl2 and CaCl2 were raised to 0.74 and 1.46 mM, respectively, to account for gluconate-divalent binding. The amount of dilute HCl or NaOH added to medium to reduce pH to target level was determined empirically.

Respiratory capacity is a measure of the maximum rate of O2 consumption and mitochondrial electron transport in a cell (46). Glycolytic capacity is the maximum rate of glucose conversion to pyruvate and/or lactate by a cell. Glucose breakdown to two lactates produces two protons, allowing for the capability of indirect measurement of glycolytic rate using the extracellular acidification (46). Compensatory glycolysis is the maximum possible rate of glycolysis in cells following inhibition of oxidative phosphorylation with rotenone/antimycin. The WE phenotype (“Warburgness”) can be expressed as the ratio of glycolysis (ECAR) to oxidative phosphorylation (expressed as the OCR) from the GST.

RNA-Seq.

High-multiplexed library preparation for RNA-seq (PLATE-seq) was performed as described previously (19). Briefly, we captured poly-adenylated mRNA from cell lysates using a 96-well plate with oligo(dT) grafted to the inner walls of each well (Qiagen). Next, we eluted the poly-adenylated mRNA and reverse transcribed using 96 different barcoded oligo(dT) primers (Integrated DNA Technologies). Following exonuclease digestion of excess primers, the barcoded cDNA libraries were pooled for second-strand synthesis and Illumina library construction. We sequenced the resulting pooled and barcoded 3′-end libraries on an Illumina NextSEq. 500.

Metabolic Profiling.

Cells were seeded in a 24-well plate in the growth medium containing 10% FBS under standard culture condition. Once cells reached 90% confluence, the growth media were removed, and cells were washed twice in PBS and incubated in 2% serum and phenol-red free medium for 24 h. Then, the media were collected for lactate production assay. L-(+)-Lactate was measured with a colorimetric kit (BioVision, K627-100) according to the manufacturer’s instructions. Absorbance (optical density 450 nm) for each sample was background corrected with the culture medium (2% FBS) collected from the well without growing cells. Final data of the lactate production rates were normalized to the protein amount per well.

Lactate and glucose concentration measurement was also done by a YSI machine and followed their protocol (YSI 2900 multianalyte system [YSI, Yellow Springs, OH]).

RNA Sequencing Data Analyses.

Bioinformatics processing and statistical analysis of RNA-seq data.

Paired-end PLATE-seq data have the sample-identifying barcode sequences in read 1 and transcript sequences in read 2. We first aligned read 2 to the hg19 human genome with University of California Santa Cruz Genome Browser (https://genome.ucsc.edu/) known genes transcriptome annotation using STAR39. Next, we demultiplexed the aligned reads based on the barcode sequence in read 1. Finally, we quantified the number of uniquely aligned reads associated with each gene in each sample using the featureCounts (47).

We next filtered genes and only kept those with >2 counts per million in at least five samples, resulting in 12,568 genes for further analysis. To account for differences in library size between the samples, trimmed mean of M values (TMM) normalization was applied, followed by data transformation using the mean–variance relationship estimated on the observed log count data as implemented in the R package voom (20, 48). This results in approximately normally distributed count data for each gene, thus allowing for standard normal theory methods to be applied. We determined there that no batch effects were present using principal component analyses. Association of gene expression with continuous measure of LPR was completed with linear regression models using the limma package (21). Genes with false discovery rate (FDR) q-values < 0.10 were considered significant (49, 50). Gene list enrichment was performed using oPOSSUM (http://opossum.cisreg.ca/oPOSSUM3/) and Enrichr (https://maayanlab.cloud/Enrichr/) (22, 23).

Evolutionary Trajectory Analysis.

Alignment of cells along their evolutionary trajectories from the parental, unselected lineage to several selected states was performed using Palantir (51), a recently published trajectory-detection algorithm for single cell RNA-seq data. Here, with single-clone RNA-seq, we had complete RNA-seq expression per clone and did not need to impute any missing data, as is done with single-cell RNA-sequencing datasets. Palantir models cell fate choices as a continuous probabilistic process over pseudotime, estimating the probability of a cell in an intermediate state to reach a terminal state (here: G, OP, and GOP). Palantir calculates differentiation potential of a given cell, leveraging the entropy over branch probabilities. Differentiation potentials near 1 correlate with earlier in the pseudotime lineage, which in this case corresponds to the parental lineage (unselected) and indicates cells with the highest potential of becoming a different phenotype over time. The high-resolution data allows for mapping of gene expression trends and dynamics over this pseudotime, which can be interpreted as how gene expression changes as the populations were driven from the parental lineage (unselected) to alternate terminal trajectories of G and OP phenotypes. Visualization of this dataset was performed using uniform manifold approximation and projection (UMAP) projections (52, 53) of the high-dimensional dataset and further analyses were overlaid onto this representation. Python code implementing Palantir on this single-clone dataset is available in SI Appendix.

Mathematical Modeling.

We used the mathematical model described in ref. 14 and extended in ref. 15 as a starting point. An interaction network and decision tree for the model are shown in SI Appendix, Fig. S6. For the in vivo simulations in Fig. 3 A–D, we set up an initial condition of a duct, as in (15). Vascular was initialized with normal vascular density outside the duct and no vessels within. For the in vitro simulations in Fig. 3 E–I, we altered the model as follows: vasculature was removed, and concentrations of oxygen, glucose, and protons (pH) were considered to be well-mixed and therefore had a global value across the simulation domain. No diffusion was necessary, and metabolic reaction rates for glucose were calculated per cell and then summed across the entire population for each time step. This lowered the concentration of glucose over time; the pH was calculated via the metabolic equations of the modelATP=(P−Vooo+ko)(gg+kG)+(Vooo+ko)(gg+kG0)+PS[1],

where P is the cell’s glycolytic phenotype while g and o are glucose and oxygen concentrations, respectively. Cells in the model were shown to survive in the UF conditions well after glucose was depleted, suggesting that a secondary survival effect was in operation. This could be due to glutamine in the media, autophagic response, etc. To account for this behavior in the model, we added a term to the equation for ATP production under the hypothesis that this behavior emerges in concert with the WE phenotype. The term is a simple linear scaling with the glycolytic phenotype (pG) of a given cell, kS pG, added to the ATP production derived from normal metabolism (SI Appendix, Fig. S9). We fit the parameter kS based on the dynamics of the population and metabolites seen in the experimental system.

Replating was mimicked in the simulation by restoring the nutrient and pH values to their initial conditions every 3 or 14 d, as per the experimental protocol for the five different conditions. Simulations were implemented using the “Hybrid Automata Library” framework (54) and barcoding visualized using the EvoFreq package in R (55). Parameters for in vivo and in vitro models are below in Tables 1 and 2 (14, 15).

Statistical Analysis.

Bioinformatics processing and statistical analysis of RNA-seq data.

Primary analysis and de-multiplexing are performed using Illumina’s CASAVA software, resulting in demultiplexed FASTQ files for subsequent analysis by the mapping software and aligner. These data will then be checked with fastqc program for quality assessment. Then cutadapt will be used to trim off adaptor contaminant sequences and low-quality bases at the ends. Reads pairs with either end too short (<25 bps) will be discarded from further analysis. Fastqc will be used again to examine characteristics of the sequencing libraries after trimming and to verify its efficiency. Next, trimmed and filtered reads will be aligned to the hg19 human transcriptome using spliced alignment transcripts to a reference (STAR) (56), followed by gene abundance estimation completed using RNA-seq by expectation maximization (57), as this approach accounts for reads mapping to multiple locations.

We next filtered genes and only kept those with >2 counts per million in at least five samples, resulting in 12,568 genes for further analysis. To account for differences in library size between the samples, TMM normalization was applied, followed by data transformation using the mean–variance relationship estimated on the observed log count data as implemented in the R package voom (20, 48). This results in approximately normally distributed count data for each gene, thus allowing for standard normal theory methods to be applied. We determined there that no batch effects were present using principal component analyses. Association of gene expression with continuous measure of LPRs was completed with linear regression models using the limma package (21). Genes with FDR q-values < 0.10 were considered differentially expressed (49, 50). Gene list enrichment was performed using oPOSSUM and Enrichr (22, 23).

Data Availability.

All study data are included in the article and supporting information.

Acknowledgments

We gratefully acknowledge funding from both the Physical Sciences Oncology Network at the National Cancer Institute (grant U54CA193489) and the Cancer Systems Biology Consortium (grant U01CA232382) as well as support from the Moffitt Center of Excellence for Evolutionary Therapy. This work was also supported partly by R01 grant R01CA077575. This work has been also supported in part by the Biostatistics and Bioinformatics Shared Resource at the H. Lee Moffitt Cancer Center & Research Institute, a National Cancer Institute-designated Comprehensive Cancer Center (P30-CA076292).

Footnotes

  • ↵1To whom correspondence may be addressed. Email: Robert.Gillies{at}moffitt.org.
  • Author contributions: M.D., R.A.G., and R.J.G. designed research; M.D., J.W., M.R.-T., L.X., P.A.S., and A.R.A.A. performed research; M.D. contributed new reagents/analytic tools; M.D., M.C.F.-F., P.A.S., E.P., B.L.F., P.M.A., and P.A.S. analyzed data; and M.D. and R.J.G. wrote the paper.

  • The authors declare no competing interest.

  • This article is a PNAS Direct Submission.

  • This article contains supporting information online at https://www.pnas.org/lookup/suppl/doi:10.1073/pnas.2011342118/-/DCSupplemental.

  • Copyright © 2021 the Author(s). Published by PNAS.

This open access article is distributed under Creative Commons Attribution-NonCommercial-NoDerivatives License 4.0 (CC BY-NC-ND).

References

  1. ↵
    1. R. A. Gatenby,
    2. R. J. Gillies
    , Why do cancers have high aerobic glycolysis? Nat. Rev. Cancer 4, 891–899 (2004).
    OpenUrlCrossRefPubMed
  2. ↵
    1. J. A. Bertout,
    2. S. A. Patel,
    3. M. C. Simon
    , The impact of O2 availability on human cancer. Nat. Rev. Cancer 8, 967–975 (2008).
    OpenUrlCrossRefPubMed
  3. ↵
    1. R. H. Thomlinson,
    2. L. H. Gray
    , The histological structure of some human lung cancers and the possible implications for radiotherapy. Br. J. Cancer 9, 539–549 (1955).
    OpenUrlCrossRefPubMed
  4. ↵
    1. M. Damaghi et al
    ., Chronic acidosis in the tumour microenvironment selects for overexpression of LAMP2 in the plasma membrane. Nat. Commun. 6, 8752 (2015).
    OpenUrlCrossRefPubMed
  5. ↵
    1. N. Rohani et al
    ., Acidification of tumor at stromal boundaries drives transcriptome alterations associated with aggressive phenotypes. Cancer Res. 79, 1952–1966 (2019).
    OpenUrlAbstract/FREE Full Text
  6. ↵
    1. R. J. Gillies,
    2. D. Verduzco,
    3. R. A. Gatenby
    , Evolutionary dynamics of carcinogenesis and why targeted therapy does not work. Nat. Rev. Cancer 12, 487–493 (2012).
    OpenUrlCrossRefPubMed
  7. ↵
    1. R. J. Gillies,
    2. R. A. Gatenby
    , Adaptive landscapes and emergent phenotypes: Why do cancers have high glycolysis? J. Bioenerg. Biomembr. 39, 251–257 (2007).
    OpenUrlCrossRefPubMed
  8. ↵
    1. R. J. Gillies,
    2. R. A. Gatenby
    , Hypoxia and adaptive landscapes in the evolution of carcinogenesis. Cancer Metastasis Rev. 26, 311–317 (2007).
    OpenUrlCrossRefPubMed
  9. ↵
    1. M. V. Liberti,
    2. J. W. Locasale
    , The Warburg effect: How does it benefit cancer cells? Trends Biochem. Sci. 41, 211–218 (2016).
    OpenUrlCrossRefPubMed
  10. ↵
    1. D. Verduzco et al
    ., Intermittent hypoxia selects for genotypes and phenotypes that increase survival, invasion, and therapy resistance. PLoS One 10, e0120958 (2015).
    OpenUrlCrossRefPubMed
  11. ↵
    1. M. Damaghi,
    2. R. Gillies
    , Phenotypic changes of acid adapted cancer cells push them toward aggressiveness in their evolution in the tumor microenvironment. Cell Cycle 16, 1739–1743(2016).
    OpenUrl
  12. ↵
    1. M. Damaghi,
    2. R. J. Gillies
    , Lysosomal protein relocation as an adaptation mechanism to extracellular acidosis. Cell Cycle 15, 1659–1660 (2016).
    OpenUrlCrossRefPubMed
  13. ↵
    1. J. Yun et al
    ., Glucose deprivation contributes to the development of KRAS pathway mutations in tumor cells. Science 325, 1555–1559 (2009).
    OpenUrlAbstract/FREE Full Text
  14. ↵
    1. M. Robertson-Tessi,
    2. R. J. Gillies,
    3. R. A. Gatenby,
    4. A. R. Anderson
    , Impact of metabolic heterogeneity on tumor growth, invasion, and treatment outcomes. Cancer Res. 75, 1567–1579 (2015).
    OpenUrlAbstract/FREE Full Text
  15. ↵
    1. A. Ibrahim-Hashim et al
    ., Defining cancer subpopulations by adaptive strategies rather than molecular properties provides novel insights into intratumoral evolution. Cancer Res. 77, 2242–2254 (2017).
    OpenUrlAbstract/FREE Full Text
  16. ↵
    1. R. A. Gatenby,
    2. R. J. Gillies
    , A microenvironmental model of carcinogenesis. Nat. Rev. Cancer 8, 56–61 (2008).
    OpenUrlCrossRefPubMed
  17. ↵
    1. C. C. Wykoff et al
    ., Expression of the hypoxia-inducible and tumor-associated carbonic anhydrases in ductal carcinoma in situ of the breast. Am. J. Pathol. 158, 1011–1019 (2001).
    OpenUrlCrossRefPubMed
  18. ↵
    1. C. C. Wykoff et al
    ., Hypoxia-inducible expression of tumor-associated carbonic anhydrases. Cancer Res. 60, 7075–7083 (2000).
    OpenUrlAbstract/FREE Full Text
  19. ↵
    1. E. C. Bush et al
    ., PLATE-seq for genome-wide regulatory network analysis of high-throughput screens. Nat. Commun. 8, 105 (2017).
    OpenUrlCrossRefPubMed
  20. ↵
    1. C. W. Law,
    2. Y. Chen,
    3. W. Shi,
    4. G. K. Smyth
    , voom: Precision weights unlock linear model analysis tools for RNA-seq read counts. Genome Biol. 15, R29 (2014).
    OpenUrlCrossRefPubMed
  21. ↵
    1. M. E. Ritchie et al
    ., Limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 43, e47 (2015).
    OpenUrlCrossRefPubMed
  22. ↵
    1. M. V. Kuleshov et al
    ., Enrichr: A comprehensive gene set enrichment analysis web server 2016 update. Nucleic Acids Res. 44, W90–W97 (2016).
    OpenUrlCrossRefPubMed
  23. ↵
    1. A. T. Kwon,
    2. D. J. Arenillas,
    3. R. Worsley Hunt,
    4. W. W. Wasserman
    , oPOSSUM-3: Advanced analysis of regulatory motif over-representation across genes or ChIP-Seq datasets. G3 (Bethesda) 2, 987–1002 (2012).
    OpenUrlCrossRef
  24. ↵
    1. B. Ye et al
    ., Klf4 glutamylation is required for cell reprogramming and early embryonic development in mice. Nat. Commun. 9, 1261 (2018).
    OpenUrlCrossRef
  25. ↵
    1. D. Wei et al
    ., KLF4 is essential for induction of cellular identity change and acinar-to-ductal reprogramming during early pancreatic carcinogenesis. Cancer Cell 29, 324–338 (2016).
    OpenUrlCrossRefPubMed
  26. ↵
    1. L. Wang,
    2. F. Shen,
    3. J. R. Stroehlein,
    4. D. Wei
    , Context-dependent functions of KLF4 in cancers: Could alternative splicing isoforms be the key? Cancer Lett. 438, 10–16 (2018).
    OpenUrlCrossRef
  27. ↵
    1. S. Avnet et al
    ., Acid microenvironment promotes cell survival of human bone sarcoma through the activation of cIAP proteins and NF-κB pathway. Am. J. Cancer Res. 9, 1127–1144 (2019).
    OpenUrl
  28. ↵
    1. K. Minami et al
    ., MiR-145 negatively regulates Warburg effect by silencing KLF4 and PTBP1 in bladder cancer cells. Oncotarget 8, 33064–33077 (2017).
    OpenUrlCrossRef
  29. ↵
    1. M. Shi et al
    ., A novel KLF4/LDHA signaling pathway regulates aerobic glycolysis in and progression of pancreatic cancer. Clin. Cancer Res. 20, 4370–4380 (2014).
    OpenUrlAbstract/FREE Full Text
  30. ↵
    1. O. Warburg
    , The chemical constitution of respiration ferment. Science 68, 437–443 (1928).
    OpenUrlFREE Full Text
  31. ↵
    1. O. Warburg
    , On respiratory impairment in cancer cells. Science 124, 269–270 (1956).
    OpenUrlCrossRefPubMed
  32. ↵
    1. O. Warburg
    , On the origin of cancer cells. Science 123, 309–314 (1956).
    OpenUrlFREE Full Text
  33. ↵
    1. O. Warburg,
    2. F. Wind,
    3. E. Negelein
    , The metabolism of tumors in the body. J. Gen. Physiol. 8, 519–530 (1927).
    OpenUrlFREE Full Text
  34. ↵
    1. J. B. Bomanji,
    2. D. C. Costa,
    3. P. J. Ell
    , Clinical role of positron emission tomography in oncology. Lancet Oncol. 2, 157–164 (2001).
    OpenUrlCrossRefPubMed
  35. ↵
    1. R. J. Gillies,
    2. I. Robey,
    3. R. A. Gatenby
    , Causes and consequences of increased glucose metabolism of cancers. J. Nucl. Med. 49 (suppl. 2), 24S–42S (2008).
    OpenUrlAbstract/FREE Full Text
  36. ↵
    1. M. Chen,
    2. J. Zhang,
    3. J. L. Manley
    , Turning on a fuel switch of cancer: hnRNP proteins regulate alternative splicing of pyruvate kinase mRNA. Cancer Res. 70, 8977–8980 (2010).
    OpenUrlAbstract/FREE Full Text
  37. ↵
    1. A. R. Maximilian et al
    ., Mix & Match: Phenotypic coexistence as a key facilitator of solid tumour invasion. Bull. Math. Biol. 82, 15 (2019).
    OpenUrl
  38. ↵
    1. M. C. Lloyd et al
    ., Darwinian dynamics of intratumoral heterogeneity: Not solely random mutations but also variable environmental selection forces. Cancer Res. 76, 3136–3144 (2016).
    OpenUrlAbstract/FREE Full Text
  39. ↵
    1. D. C. Sgroi
    , Preinvasive breast cancer. Annu. Rev. Pathol. 5, 193–221 (2010).
    OpenUrlCrossRefPubMed
  40. ↵
    1. M. Hu,
    2. K. Polyak
    , Microenvironmental regulation of cancer development. Curr. Opin. Genet. Dev. 18, 27–34 (2008).
    OpenUrlCrossRefPubMed
  41. ↵
    1. K. O. Osuala et al
    ., Il-6 signaling between ductal carcinoma in situ cells and carcinoma-associated fibroblasts mediates tumor cell growth and migration. BMC Cancer 15, 584 (2015).
    OpenUrlCrossRefPubMed
  42. ↵
    1. C. Strell et al
    ., Impact of epithelial-stromal interactions on peritumoral fibroblasts in ductal carcinoma in situ. J. Natl. Cancer Inst. 111, 983–995 (2019).
    OpenUrl
  43. ↵
    1. M. Hu et al
    ., Regulation of in situ to invasive breast carcinoma transition. Cancer Cell 13, 394–406 (2008).
    OpenUrlCrossRefPubMed
  44. ↵
    1. X. Y. Chen et al
    ., Higher density of stromal M2 macrophages in breast ductal carcinoma in situ predicts recurrence. Virchows Arch. 476, 825–833 (2020).
    OpenUrl
  45. ↵
    1. N. Linde et al
    ., Macrophages orchestrate breast cancer early dissemination and metastasis. Nat. Commun. 9, 21 (2018).
    OpenUrlCrossRef
  46. ↵
    1. S. A. Mookerjee,
    2. D. G. Nicholls,
    3. M. D. Brand
    , Determining maximum glycolytic capacity using extracellular flux measurements. PLoS One 11, e0152016 (2016).
    OpenUrlCrossRefPubMed
  47. ↵
    1. Y. Liao,
    2. G. K. Smyth,
    3. W. Shi
    , featureCounts: An efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics 30, 923–930 (2014).
    OpenUrlCrossRefPubMed
  48. ↵
    1. M. D. Robinson,
    2. A. Oshlack
    , A scaling normalization method for differential expression analysis of RNA-seq data. Genome Biol. 11, R25 (2010).
    OpenUrlCrossRefPubMed
  49. ↵
    1. J. D. Storey,
    2. R. Tibshirani
    , Statistical significance for genomewide studies. Proc. Natl. Acad. Sci. U.S.A. 100, 9440–9445 (2003).
    OpenUrlAbstract/FREE Full Text
  50. ↵
    1. J. D. Storey
    , A direct approach to false discovery rates. J. R. Stat. Soc. B 64, 479–498 (2002).
    OpenUrlCrossRef
  51. ↵
    1. M. Setty et al
    ., Characterization of cell fate probabilities in single-cell data with Palantir. Nat. Biotechnol. 37, 451–460 (2019).
    OpenUrlCrossRefPubMed
  52. ↵
    1. E. Becht et al
    ., Dimensionality reduction for visualizing single-cell data using UMAP. Nat. Biotechnol., doi:10.1038/nbt.4314 (2018).
    OpenUrlCrossRefPubMed
  53. ↵
    1. L. McInnes,
    2. J. Healy,
    3. N. Saul,
    4. L. Großberger
    , UMAP: Uniform Manifold approximation and projection. J. Open Source Softw. 3, (2018).
  54. ↵
    1. R. R. Bravo et al
    ., Hybrid Automata library: A flexible platform for hybrid modeling with real-time visualization. PLoS Comput. Biol. 16, e1007635 (2020).
    OpenUrlCrossRef
  55. ↵
    1. C. D. Gatenbee,
    2. R. O. Schenck,
    3. R. R. Bravo,
    4. A. R. A. Anderson
    , EvoFreq: Visualization of the evolutionary frequencies of sequence and model data. BMC Bioinformatics 20, 710 (2019).
    OpenUrlCrossRef
  56. ↵
    1. A. Dobin et al
    ., STAR: Ultrafast universal RNA-seq aligner. Bioinformatics 29, 15–21 (2013).
    OpenUrlCrossRefPubMed
  57. ↵
    1. B. Li,
    2. C. N. Dewey
    , RSEM: Accurate transcript quantification from RNA-seq data with or without a reference genome. BMC Bioinformatics 12, 323 (2011).
    OpenUrlCrossRefPubMed
PreviousNext
Back to top
Article Alerts
Email Article

Thank you for your interest in spreading the word on PNAS.

NOTE: We only request your email address so that the person you are recommending the page to knows that you wanted them to see it, and that it is not junk mail. We do not capture any email address.

Enter multiple addresses on separate lines or separate them with commas.
The harsh microenvironment in early breast cancer selects for a Warburg phenotype
(Your Name) has sent you a message from PNAS
(Your Name) thought you would like to see the PNAS web site.
CAPTCHA
This question is for testing whether or not you are a human visitor and to prevent automated spam submissions.
Citation Tools
The harsh microenvironment in early breast cancer selects for a Warburg phenotype
Mehdi Damaghi, Jeffrey West, Mark Robertson-Tessi, Liping Xu, Meghan C. Ferrall-Fairbanks, Paul A. Stewart, Erez Persi, Brooke L. Fridley, Philipp M. Altrock, Robert A. Gatenby, Peter A. Sims, Alexander R. A. Anderson, Robert J. Gillies
Proceedings of the National Academy of Sciences Jan 2021, 118 (3) e2011342118; DOI: 10.1073/pnas.2011342118

Citation Manager Formats

  • BibTeX
  • Bookends
  • EasyBib
  • EndNote (tagged)
  • EndNote 8 (xml)
  • Medlars
  • Mendeley
  • Papers
  • RefWorks Tagged
  • Ref Manager
  • RIS
  • Zotero
Request Permissions
Share
The harsh microenvironment in early breast cancer selects for a Warburg phenotype
Mehdi Damaghi, Jeffrey West, Mark Robertson-Tessi, Liping Xu, Meghan C. Ferrall-Fairbanks, Paul A. Stewart, Erez Persi, Brooke L. Fridley, Philipp M. Altrock, Robert A. Gatenby, Peter A. Sims, Alexander R. A. Anderson, Robert J. Gillies
Proceedings of the National Academy of Sciences Jan 2021, 118 (3) e2011342118; DOI: 10.1073/pnas.2011342118
Digg logo Reddit logo Twitter logo Facebook logo Google logo Mendeley logo
  • Tweet Widget
  • Facebook Like
  • Mendeley logo Mendeley
Proceedings of the National Academy of Sciences: 118 (3)
Table of Contents

Submit

Sign up for Article Alerts

Article Classifications

  • Biological Sciences
  • Evolution

Jump to section

  • Article
    • Abstract
    • Results
    • Discussion
    • Methods
    • Data Availability.
    • Acknowledgments
    • Footnotes
    • References
  • Figures & SI
  • Info & Metrics
  • PDF

You May Also be Interested in

Surgeons hands during surgery
Inner Workings: Advances in infectious disease treatment promise to expand the pool of donor organs
Despite myriad challenges, clinicians see room for progress.
Image credit: Shutterstock/David Tadevosian.
Setting sun over a sun-baked dirt landscape
Core Concept: Popular integrated assessment climate policy models have key caveats
Better explicating the strengths and shortcomings of these models will help refine projections and improve transparency in the years ahead.
Image credit: Witsawat.S.
Double helix
Journal Club: Noncoding DNA shown to underlie function, cause limb malformations
Using CRISPR, researchers showed that a region some used to label “junk DNA” has a major role in a rare genetic disorder.
Image credit: Nathan Devery.
Steamboat Geyser eruption.
Eruption of Steamboat Geyser
Mara Reed and Michael Manga explore why Yellowstone's Steamboat Geyser resumed erupting in 2018.
Listen
Past PodcastsSubscribe
Birds nestling on tree branches
Parent–offspring conflict in songbird fledging
Some songbird parents might improve their own fitness by manipulating their offspring into leaving the nest early, at the cost of fledgling survival, a study finds.
Image credit: Gil Eckrich (photographer).

Similar Articles

Site Logo
Powered by HighWire
  • Submit Manuscript
  • Twitter
  • Facebook
  • RSS Feeds
  • Email Alerts

Articles

  • Current Issue
  • Special Feature Articles – Most Recent
  • List of Issues

PNAS Portals

  • Anthropology
  • Chemistry
  • Classics
  • Front Matter
  • Physics
  • Sustainability Science
  • Teaching Resources

Information

  • Authors
  • Editorial Board
  • Reviewers
  • Librarians
  • Press
  • Site Map
  • PNAS Updates

Feedback    Privacy/Legal

Copyright © 2021 National Academy of Sciences. Online ISSN 1091-6490