Multistate Gene Cluster Switches Determine the Adaptive Mitochondrial and Metabolic Landscape of Breast Cancer

Abstract

Adaptive metabolic switches are proposed to underlie conversions between cellular states during normal development as well as in cancer evolution. Metabolic adaptations represent important therapeutic targets in tumors, highlighting the need to characterize the full spectrum, characteristics, and regulation of the metabolic switches. To investigate the hypothesis that metabolic switches associated with specific metabolic states can be recognized by locating large alternating gene expression patterns, we developed a method to identify interspersed gene sets by massive correlated biclustering and to predict their metabolic wiring. Testing the method on breast cancer transcriptome datasets revealed a series of gene sets with switch-like behavior that could be used to predict mitochondrial content, metabolic activity, and central carbon flux in tumors. The predictions were experimentally validated by bioenergetic profiling and metabolic flux analysis of 13C-labeled substrates. The metabolic switch positions also distinguished between cellular states, correlating with tumor pathology, prognosis, and chemosensitivity. The method is applicable to any large and heterogeneous transcriptome dataset to discover metabolic and associated pathophysiological states. Significance: A method for identifying the transcriptomic signatures of metabolic switches underlying divergent routes of cellular transformation stratifies breast cancer into metabolic subtypes, predicting their biology, architecture, and clinical outcome.

Figure 1.

Large-scale transcriptome switches identify breast cancer subtypes. A, Gene set switch identification in bicluster samples. MCbiclust starts with a random gene set in a random sample set and performs a greedy search aiming for the sample set, across which, the gene set has the highest correlation. As the absolute value of correlation is maximized during the search, the resulting top correlating samples display two gene sets with sharp anticorrelation (switch) between them. The gene sets and sample sets can be extended to the whole transcriptome and sample set by quantifying correlations (see main text), which provides the basis for gene set enrichment analysis. B, Visualization of individual switches (MB1, MB2, MB3 biclusters) in all METABRIC samples. Two-dimensional distribution by PC1 and ranking index are shown, overlaid with PAM50 classifiers. Dashed lines represent the thresholds to include samples in the forks, calculated as 0.04 × forkscale (see Materials and Methods). C, UpSet matrix and plots of sample intersections between top and bottom positions of each bicluster switch. Unique and overlapping (intersect) samples are shown separately. D, Scheme of switches. MB1 and MB2 share their LF, thus representing a multistate switch, while MB3 is an independent bistate switch, thus the MB1_2 and MB_3 biclusters represent different switch mechanisms. E and F, Interswitch relationships are shown in two-dimensional distribution plots of samples by two biclusters (switches). The axes represent the scale between the UF and LF of each bicluster. The distribution of samples according to their MB1 vs. MB2 forkscale value (see Materials and Methods) is shown in E, while the combinations of MB3 with MB1 or MB2 are shown in F. G, Heatmap of the correlation values (CVs, Pearson r) across the whole transcriptome in the MB1 bicluster; a small set of midrange values were hidden (see gaps in heatmap). Annotations show the locations of mitochondrial and metabolic genes (top), and bicluster genes (bottom). For comparison between the switches, bottom annotation plots show the CV values of each sample in the three biclusters (switches). H, Fold enrichment of mitochondrial and metabolic (with intersections) genes in the anticorrelated gene sets. Heatmap shows enrichment values of each gene set in each position of all switches. Gene set intersections are on the Venn diagram below the heatmap.

Figure 2.

Metabolic transcriptome switches override PAM50 classification. A–D, Visualization of bicluster switches based on CV values of individual genes in the METABRIC breast cancer dataset analysis in the indicated metabolic pathways. Genes are ranked by the MB1 CV values and directly compare the MB1 and MB2 biclusters on the left side of each panel. The ranked MB3 values are shown on the right side of each panel. E–H, Sub-pathway gene clusters group according to the bicluster switches across PAM50 categories. Bicluster forks samples were grouped according to their intrinsic subtype and the average transcript levels of each gene in the indicated pathways were calculated. Hierarchical clustering of genes and sample groups was performed (Pearson distance and ward2 clustering). Group classifications according to the biclusters and intrinsic breast cancer subtypes are shown on the right.

Figure 3.

Oncogenic signaling and cell of origin associates with the biclusters. A, Heatmap of distribution of the top 250 most variable PARADIGM concepts across the bicluster switches in the TCGA-BRCA dataset. The term values were obtained from Berger and colleagues (30) pan-cancer analysis and hierarchically clustered. Annotations show proliferation, tumor purity, switch identity (UF and LF of MB_1, MB_2, or MB_3), Pam50, and histology classifications. Transcriptomic signatures based on PARADIGM concepts are shown on the right. B, Association of mammary gland cell states (37) and transcriptomic switches in the METABRIC dataset. Two-dimensional distribution plots of samples by two biclusters are shown overlaid with cell-of-origin marker GSVA values. C, Scheme of signaling of the multistate switch in the MB1 and MB3 biclusters.

Figure 4.

Mitochondrial function: MB1_UF cells show active functional biogenesis in contrast to the MB1_LF phenotype. A, Two-dimensional distribution (PC1 vs. ranking index, shown for the MB1 bicluster switch) of biopsies obtained for generating PDX samples from the study by Bruna and colleagues (38), overlaid in left panel with Pam50 categories and in the right panel with PDX status indicating whether the sample has grown from the biopsy to PDX status. B, Cell line scoring to sort cell lines to the MB1 switches. Volcano plot of the scores and −log₁₀ of the P values are shown. The cell lines chosen for experimental work are highlighted in red. C, OCR in cell lines grouped by their switch state (MB1_UF/MB1_LF). Reserve, ATP synthesis coupled, and uncoupled respirations were quantified using FCCP and oligomycin treatments. Two-way ANOVA with Sidak multiple comparisons test. D, Average cellular TMRM intensity grouped by cell line switch state (MB1_UF/MB1_LF). Data were obtained from high content microscopy from >1K cells per group. Unpaired Student t test. E, Relative expression of respiratory chain subunits quantified from Western blot analysis. Results are grouped by switch state (MB1_UF/MB1_LF). Two-way ANOVA with Sidak multiple comparisons test. F, High content imaging analysis of mitochondrial biogenesis, structure, and function in MB1_UF and MB1_LF cells. Scale bars, 10 μm. Quantifications of the picoGreen integrated intensity (readout for mtDNA content), mitochondrial number, total volume, and individual mitochondrial volume are shown. Unpaired Student t tests. ns, nonsignificant, P > 0.05; *, P ≤ 0.05; **, P ≤ 0.01; ***, P ≤ 0.001; ****, P ≤ 0.0001.

Figure 5.

Metabolic wiring of the switch positions probed by ¹³C-labeled substrates. A, Nutrient supply and ¹³C-labeling conditions. High and low glutamine (H/L = 1/0.1 mmol/L) levels were combined with the presence and absence of pyruvate (+/−, 1/0 mmol/L), and glucose was kept at 10 mmol/L. B, Glucose uptake, lactate secretion, and fractional enrichment of m + 3 labeled intracellular lactate. Cells were grouped by switch state. Two-way ANOVA with Sidak multiple comparisons test. C, Fractional enrichment of the m+3 pyruvate isotopomer, following ^U13C-glucose labeling for 24 hours, in cells grouped by switch state. Two-way ANOVA with Sidak multiple comparison tests. D, Overall model of the MB1_UF/MB1_LF metabolic switch in central carbon metabolism. E, Relative carbon contribution to citrate following 24 hours labeling with ^U13C-glucose and ^U13C-glutamine, in cells grouped by switch state. Two-way ANOVA with Sidak multiple comparison tests. F, Fractional enrichment of the m+2 (PDH route; left) and m+5 (PDH + PC route; right) isotopologues in citrate following 24 hours labeling with ^U13C-glucose, in cells grouped by switch state. Two-way ANOVA with uncorrected Fisher LSD multiple comparison tests. G, Glutamine uptake, glutamate secretion, and fractional enrichment of malate m+4 (oxidative route) and citrate m+5 (reductive carboxylation) isotopologues following 24 hours labeling with ^U13C- glutamine (right panels). Two-way ANOVA with uncorrected Fisher LSD multiple comparison tests. H, Substrate dependence of basal oxygen consumption rate and total cellular ATP content in cells grouped by switch state. Two-way ANOVA with uncorrected Fisher LSD multiple comparison tests. ns, nonsignificant, P > 0.05; *, P ≤ 0.05; **, P ≤ 0.01; ***, P ≤ 0.00; ****, P ≤ 0.0001.

Figure 6.

Protein expression profiles match transcriptomic switches. A, Western blot analysis of 29 metabolic enzymes across five cell lines representing MB1_UF and MB1_LF switch positions and cell states. Representative blots from n ≥ 3 independent experiments, used for quantification of signals, normalizing density values to either grp75 or β-actin. B, D, and E, Clustering and quantitative analysis of enzyme expression levels obtained from A. Values were separately normalized to either mitochondrial (left) or cytosolic (right) markers, grp75 or β-actin. The comparison of mitochondrial abundance based on the grp75/β-actin ratio is shown at the bottom in B. Heatmaps are based on a pairwise, robust to outliers distance measure of the mean protein expression values. The k-means clustering. D and E, Two-way ANOVA with uncorrected Fisher LSD multiple comparison tests. C, Heatmap and k-means clustering of transcript CV values of the 29 enzymes from the METABRIC dataset (see Figs. 1 and 2, positive values predicting higher expression of genes in MB1_UF and negative values corresponding to higher expression in MB1_LF) and scaled protein expression levels (wCV) averaged across all cell lines, normalized to grp75 (wCV_mito) or β-actin (wCV_cyto). ns, nonsignificant, P > 0.05; *, P ≤ 0.05; **, P ≤ 0.01; ***, P ≤ 0.001; ****, P ≤ 0.0001).

Figure 7.

Pathological and clinical predictions from metabolic gene switches. A, The χ² distribution of discrete histological features across the upper, lower, and non-fork samples of the three biclusters (MB1–MB3). Heatmap shows k-means clustering of χ² test −log₁₀P values, dividing the histological parameters into three groups: high, low, and no association with the bicluster. B, Visualization of histological feature distribution on the MB1 and MB2 switches. Classifiers of histological tumor type (top) and degree of nuclear pleomorphism are overlaid onto two-dimensional distribution plots of TCGA samples. The axes represent the range between the UF and LF of each bicluster. C, Clustering of histological and genetic features of samples belonging to the MB1 and MB2 switches. Fractions of samples for each feature were calculated and clustered with a custom robust to outliers distance measure function and k-means clustering. Genetic features representative of each switch position are highlighted (bold), Pam50 categories (red), and histology types and features (blue and orange, respectively). D, Key genetic and transcriptional architecture of the MB1 and MB2 switches. E–L, Survival analysis of samples in the MB1 switch. Kaplan–Meier survival plots and Cox proportional hazard analysis results from ref. (E–G) and ref. (H–L) are shown, generated by the kmplot.com tool (43). MB1_UF and MB1_LF bicluster (E–J) and OXPHOS (K and L) patterns (Supplementary Table S1) were used to stratify subsets of ER-positive, Her2-negative samples as indicated on the panels.