The research strategy is presented in Figure 1 . RNA sequencing data were downloaded from the National Central of Biology Information Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/), including GSE43837, GSE14017, and GSE14018 (11). GSE43837 contains RNA sequence for 19 BCBMs and 19 BCs, GSE14017 contains RNA sequence for 15 BCBMs and 14 extracranial metastases (BCEMs), and GSE14018 contains RNA sequence for 7 BCBMs and 29 BCEMs. Microarray annotation information was used to match probes with corresponding genes. The median expression value was calculated out for the gene matched with more than one probe. We first performed the immune and metabolic analysis on BCBMs and BCs of GSE43837 and then performed a similar analysis on BCBMs and BCEMs of GSE14017 and GSE14018, respectively.

We utilized the Estimation of Stromal and Immune cells in Malignant Tumor tissues using Expression data (ESTIMATE) and Microenvironment Cell Populations-Counter (MCP-counter) R package to characterize immune infiltration in samples. ESTIMATE can infer the proportion of immune cells and stromal cells in tumor samples using gene expression (12). However, ESTIMATE cannot identify the distinct immune cell populations in heterogeneous tissues. In contrast, MCP-counter can quantify the absolute abundance of eight immune cells in heterogeneous tissues using transcriptome data (13).

WGCNA R package were used to construct a weight co-expression network associated with immune infiltration (14). First, based on the Pearson’s correlation value between paired genes, the expression levels of individual transcripts were converted into a similarity matrix. Next, we picked a proper soft threshold power that can increase strong correlations and decrease weak correlations between genes. The adjacency matrix was then converted into a topological overlap matrix when the soft threshold power β = 6. Then, the gene set was divided into several modules with similar expression patterns. Module–trait associations referred to the correlation between the module eigengene and the immune infiltration.

Differentially expressed genes (DEGs) in different groups were identified using edgeR package (15). Specifically, edgeR adjusts gene expression according to different sequencing depths as represented by varying libraries. The Log2 fold-change (Log2FC) is an estimate of the log2 ratio of expression in a cluster to other clusters. A value of 1.0 indicates twofold greater expression in the cluster of interest. The exact test that adapted for the negative binomially distributed counts was chosen to judge the significance for DEGs. Adjusted p-values or false discovery rate (FDR) was determined by the default Benjamini–Hochberg (BH) correction in edgeR. For selecting the top features in a dataset, FDR < 0.05 and fold change (FC) > 1.5 were set as the cutoff criteria.

The Database for Annotation, Visualization and Integrated Discovery (DAVID) v6.8 (https://david.ncifcrf.gov/summary.jsp) database integrates biological data and functional annotation tools to provide systematic and comprehensive biological function annotations for large-scale gene or protein lists. It was used to identify enriched Gene Ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways of DEGs (16). p-values are determined by the Fisher’s exact test in DAVID. Adjusted p-values were determined by BH correction in DAVID. For selecting significant pathways, p-value < 0.05 and FDR < 0.25 were set as the cutoff criteria.

We crossed the genes co-expressed with immune infiltration determined by the WGCNA package with DEGs to obtain gene signatures related to immune infiltration in the cluster of interest. The gene expression values for those signatures were then averaged to form the Immune metagene (IM-metagene). Specific genes were indicated in Supplementary Table S4 .

GSEA determines whether an a priori defined set of genes has statistically significant difference in expression under two different biological conditions (17). GSEA software 3.0 downloaded from the Broad Institute was used for enrichment analysis for our datasets. The gene set of “c2.cp.kegg.v7.1.symbols.gmt”, which summarizes and represents specific, well-defined KEGG metabolic pathways, was downloaded from the Molecular Signatures Database (http://software.broadinstitude.org/gsea/msigdb/index.jsp). The normalized enrichment score (NES) represented the degree of enriched KEGG pathways in cluster of interest. p-values corresponding to each NES were determined by the Fisher’s exact test (1,000 permutations) in GSEA. Adjusted p-values were determined by BH correction in GSEA. For selecting significant pathways, FDR < 0.25 was set as the cutoff criteria.

We crossed core genes in OXPHOS enrichment in interested cluster determined by GSEA with DEGs to obtain signatures related to OXPHOS enrichment in the cluster of interest. The gene expression values for those signatures were then averaged to form the OXPHOS metagene (OP-metagene). Specific genes were indicated in Supplementary Table S7 .

Kaplan–Meier (KM) plotter [Breast cancer] is an online survival analysis tool that can assess the prognostic function of 22,277 genes in breast cancer patients using microarray data (http://kmplot.com/analysis/index.php?p=background) (18). All KM plots were displayed using the “auto select best cutoff” parameter. Relapse-free survival (RFS), overall survival (OS), and distant metastasis-free survival (DMFS) were selected as the endpoints. Hazard ratio (HR) was considered significant when log rank p-value < 0.05. The corresponding 95% confidence intervals (95% CI) were also displayed on all KM plots.

MDA-MB-231 (human breast cancer cell line) cells were purchased from Procell Life Science & Technology Co. Ltd. MDA-MB-231 cells were cultured in Dulbecco’s modified Eagle’s medium (DMEM; BasalMedia, cat. no. L110KJ) supplemented with 10% fetal bovine serum (FBS; gibco, cat. no. A3160801) and 1% penicillin–streptomycin (BasalMedia, cat. no. S110JV) in a 95% humidified incubator containing 5% CO2 at 37°C.

Cell proliferation assay. MDA-MB-231 cells (4 × 10³) were seeded on 96-well plates. After the cells adhered to the wall, the cells were treated with 1.0 µM oligomycin [Oligo(1.0)] and incubated in a 5% CO₂ incubator at 37°C; 10 µl Cell Counting Kit-8 (CCK8; APEXBIO, cat. no. K1018) solution was then added into each well at 0 h, 12 h, 24 h, and 48 h, respectively, and cultured for 2 h. Next, the 96-well plates were put on the enzyme-linked immunoassay instrument and shaken for 2 s. The absorbance was measured at 460 nm. The growth rate was calculated as follows: Growth rate of Control = ABS(OD value of Control − mean(OD value of Control group))/mean(OD value of Control group), Growth rate of Oligo(1.0) = ABS(OD value of Oligo(1.0) − mean(OD value of Control group))/mean(OD value of Control group).

Cell apoptosis assay. Annexin V and propidium iodide (PI) (BD pharmingen, cat. no. 556547) were used to stain the cells cultured in medium. FSC-H and SSC-H of flow cytometry were used to detect single cells. The percentage of annexin V⁻/PI⁻ cells was used to represent the cell viability.

Scratch assay. MDA-MB-231 cells (1 × 10⁶) were seeded on six-well plates. When the cell confluence reached 95%, the fused cells were scratched along the pore diameter with a sterile 200-µl pipette tip and then washed five times with PBS to remove floating cells and debris. The medium in each well was replaced with serum-free medium containing 1.0 µM oligomycin. The wound healing was observed at 0 and 48 h, and photos were taken under a microscope.

We used edgeR package to identify DEGs between BCBMs and BCs of GSE43837. A total of 539 DEGs were identified, of which 394 protein-coding genes were upregulated and 145 protein-coding genes were downregulated in BCBMs compared with BCs, respectively (FDR < 0.05, FC > 1.5; Supplementary Table S1 ).

Then, we performed immune analysis on BCBMs and BCs of GSE43837. We utilized the ESTIMATE and MCP-counter R packages to characterize differences in immune cell infiltration between BCBMs and BCs (GSE43837). ESTIMATE is a tool used to infer tumor purity and immune infiltration from gene-expression data that were originally validated in 11 cancer types (12). However, ESTIMATE can only assess the overall immune status of the tumor. On the contrary, MCP-counter can calculate the specific infiltration of T cells, CD8 T cells, cytotoxic lymphocytes, B lineage, NK cells, monocytic lineage, myeloid dendritic cells, and neutrophils in tumors based on gene expression (13). Together, ESTIMATE assessed that the immune score of BCBMs was lower than that of BC, although there was no statistical significance (p = 0.1542; Figure 2A ); MCP-counter estimated that the infiltration of eight immune cells in BCBMs was also lower than that of BCs; in particular, the infiltration of B lineage (p < 0.05; Figure 2B ) and myeloid dendritic cells (p < 0.05; Figure 2B ) in BCBMs was significantly lower than that of BCs. As the immune infiltration is lower in BCBMs compared with BCs, and the expression of PDL1 and PTEN has been confirmed to be related to tumor immune infiltration in previous studies (19, 20), we also compared the expression of PDL1 and PTEN between BCBMs and BCs. The RNA expression of PDL1 was not different in BCBMs compared with BCs (p = 0.1328; Figure 2C ). The RNA expression of PTEN in BCBMs was lower than that of BCs at the limit of significance (p = 0.0571; Figure 2D ).

We used Weighted Correlation Network Analysis (WGCNA) R package to search for genes related to the immune infiltration of GSE43837. WGCNA R package is an effective tool that can be used to mine hub modules with similar expression patterns related to clinical traits (14). To build a scale-free network, we picked β = 6 (scale-free R ² = 0.86) as the soft-thresholding power ( Figure 3A ). Then, those genes were classified into 16 modules ( Figure 3B ). As previous immune infiltration analysis identified B lineage and myeloid dendritic cell infiltration significantly decreased in BCBMs compared with BCs and the turquoise module had the highest correlation with B lineage (r = 0.71, p = 7e−07; Figure 3C ) and myeloid dendritic cells (r = 0.49, p = 0.002; Figure 3C ), the turquoise module was identified as a hub module significantly related to the immune infiltration of GSE43837 samples. To obtain the core immune signatures associated with BCBMs, we crossed the genes in turquoise module with DEGs. In total, we obtained 30 immune signatures (KRTAP4-9, BNC2, GUCA2B, BMP15, MDGA2, OTOP2, OSBP2, ZNF768, NUDT18, ABRA, KRT37, RHOC, COL8A1, GJA8, WFDC10B, GOLIM4, ASCC2, KITLG, ACOT4, BARX1, KCNC3, C6orf163, ACHE, HSD17B4, BATF3, CD1B, ZNRF4, C1orf158, OR2H2, and VCX2; Figure 3D and Supplementary Table S4 ), all of which were downregulated in BCBMs compared with BCs ( Figure 4F ). The gene expression values for all those signatures were then averaged to form the Immune metagene (IM-metagene).

Survival analysis for IM-metagene was performed in the KM plotter [Breast cancer] database. This was done to determine whether the expression of IM-metagene is related to the biological malignant behavior of breast cancer and whether IM-metagene can be used as a prognostic indicator for patients with breast cancer. KM plotter [Breast cancer] showed a significant decrease of RFS (HR = 0.7, log rank p = 5.3e−06), OS (HR = 0.66, log rank p = 0.011), and DMFS (HR =0.66, log rank p = 0.012) with lower expression of IM-metagene in patients with breast cancer ( Figures 3E–G ). According to molecular classification, breast cancer is divided into three subtypes: luminal epithelial type (luminal type), HER2 overexpression (HER2+) type, and basal-like type. Basal-like type molecules are expressed as ER(−)/PR(−)/HER2(−), which is equivalent to triple-negative breast cancer. Different breast cancer subtypes could vary for the prognosis and adjuvant treatments. Further exploring the relationship between IM-metagene and the prognosis of patients with breast cancer subtypes, we did not identify significant correlation between the expression of IM-metagene and the prognosis of patients with different breast cancer subtypes.

To explore the biological and metabolic features of BCBMs, we used the DAVID tool to analyze the enrichment of GO and KEGG pathways of DEGs. DAVID is an online tool to analyze the biological function and the enrichment of KEGG pathways using gene lists. However, biological regulation is a progressive relationship; small changes in upstream genes may lead to obvious changes in downstream genes. If you use a set threshold to screen DEGs and then perform function/pathway enrichment analysis (GO/KEGG) directly, some gene information will be lost, which may result in missing significant biological and metabolic pathways. Therefore, we performed GSEA in the GSE43837 dataset. GSEA does not require a fixed threshold to filter genes. It is a method based on all-gene expression analysis and avoids the shortcomings of traditional enrichment analysis methods. Because there were not many DEGs in BCBMs compared with BCs, if the threshold was set to FDR < 0.05, a lot of GO terms will be missed. Therefore, we set the screening conditions as p-value < 0.05, FDR < 0.25. The top GO terms for BCBMs included protein folding (CCT3, LRPAP1, TRAP1, LMAN2L, NFYC, TBCC, DNAJB2, GNAO1, MLEC, ERP27, CCT7, CRYAB, PPIA, PFDN5, SIL1, and AARS; Figure 4A and Supplementary Table S2 ) and negative regulation of canonical Wnt signaling pathway (PSMB6, PSMB4, PSMB2, FRZB, HDAC1, DDIT3, PSMD2, UBC, PSMB1, KREMEN2, SOX9, and PFDN5; Figure 4A and Supplementary Table S2 ). Interestingly, GO analysis also showed that overexpressed genes in BCBMs compared with BCs were mainly enriched in cellular components related to oxidative phosphorylation (OXPHOS), such as mitochondria, mitochondrial matrix, mitochondrial inner membrane, proton transport ATP synthase complex, and catalytic core F (1) (p < 0.05, FDR < 0.25; Figure 4A and Supplementary Table S2 ); the top GO terms for BCs include cell–cell signaling (CCR1, CXCL10, GJB2, CXCL9, FGFBP1, CCL8, SH2D1A, IHH, and BARX1) and collagen catabolic process (MMP12, MMP11, COL3A1, MMP13, and MMP1) and other extracellular pathways (p < 0.05, FDR < 0.25; Figure 4A and Supplementary Table S2 ). KEGG pathway analysis using DAVID database showed that upregulated genes were enriched in Alzheimer’s disease and downregulated genes were enriched in cytokine–cytokine receptor interaction (p < 0.05, FDR < 0.25; Figure 4B and Supplementary 2) in BCBMs compared with BCs. GSEA detected the significant enrichment of Parkinson’s disease and OXPHOS (FDR < 0.25; Figure 4C and Supplementary Table S5 ) in BCBMs compared with BCs using the c2.cp.kegg.v7.1.symbols.gmt gene sets. Considering that the enrichment of Parkinson’s disease may be due to the contamination of the surrounding brain tissue, and the GO analysis identified many cell components related to OXPHOS enriched in BCBMs compared with BCs, our next step was mainly focused on OXPHOS ( Figure 4D ).

To obtain the core signatures related to OXPHOS enrichment of BCBMs, we crossed DEGs with 49 core genes in OXPHOS enrichment in BCBMs compared with BCs determined by GSEA ( Supplementary Table S6 ), and obtained six signatures (COX6B1, UQCRFS1, COX4I1, NDUFV1, ATP6V0A1, and NDUFA9; Figure 4E and Supplementary Table S7 ), respectively. All the six OXPHOS signatures were upregulated in BCBMs compared with BCs ( Figure 4F ). The gene expression values for all those signatures were then averaged to form the OXPHOS metagene (OP-metagene).

Next, we performed a series of survival analyses in patients with breast cancer using microarray data in the KM Plotter [Breast cancer] database. KM plotter [Breast cancer] is an analytical database that can be used to determine whether gene expression is statistically related to the prognosis of breast cancer patients. This was done to determine whether OP-metagene was higher in more biologically aggressive tumors and whether it has value as predictive biomarkers for disease progression in patients. Remarkably, the analysis identified significant decrease of RFS (HR = 1.83, log rank p < 1e−16), OS (HR = 1.67, log rank p = 3.1e−06), and DMFS (HR = 1.33, log rank p = 0.0041) with higher expression of OP-metagene in breast cancers ( Figure 4H ). Then, we continued to explore whether OP-metagene is significantly associated with the prognosis of different breast cancer subtypes. The results showed that the higher the expression of OP-metagene, the shorter the RFS of breast cancer patients with luminal A (HR = 1.72, log rank p = 2.7e−10), luminal B (HR = 1.97, log rank p = 2.7e−12), basal-like (HR=1.52, log rank p = 0.0023), and HER2+ (HR = 1.45, log rank p = 0.06, in the edge of significance) ( Figure 4I ) subtypes, suggesting that OP-metagene may be a better biomarker for predicting disease progression in patients with breast cancer than IM-metagene. We also tried to explore the cause for the enrichment of OXPHOS in BCBMs compared with BCs. As a previous study has reported that PGC1A mediates mitochondrial biosynthesis and OXPHOS in cancer cells to promote metastasis (21), we compared the RNA expression of PGC1A in BCBMs and BCs, but there was no difference between the two clusters (p = 0.5204; Figure 4G ).

Exploratory immune and metabolic analysis was performed on BCBMs and extracranial metastases of GSE14017 and GSE14018, respectively, including lung metastases, bone metastases, and liver metastases. ESTIMATE-identified immune scores decreased significantly in BCBMs compared with BCEMs of GSE14017 (p = 4.064e−05; Figure 5A ) and GSE14018 (p = 0.0202; Figure 5C ) using RNA sequence. Specifically, MCP-counter-identified T cells (p < 0.05), Cytotoxic lymphocytes (p < 0.05), Monocytic lineage (p < 0.01), and Neutrophils (p < 0.05) infiltration decreased significantly in BCBMs compared with BCEMs of GSE14017 ( Figure 5B ), and there were more significant decreases in T cells (p < 0.01), Cytotoxic lymphocytes (p < 0.001), Myeloid dendritic cells (p < 0.01), and Neutrophils (p < 0.001) infiltration in BCBMs compared with BCEMs of GSE14018 ( Figure 5D ). However, MCP-counter failed to detect the infiltration of CD8 T cells of GSE14018. Exploratory PTEN RNA expression comparison was also performed on a small cohort RNA-seq of BCBMs versus BCEMs (GSE14017 and GSE14018). PTEN RNA expression was significantly decreased in BCBMs compared with BCEMs of GSE14017 (p = 0.0292; Figure 5E , left) and GSE14018 (p = 0.0014; Figure 5E , right). We also performed GSEA in BCBMs versus BCEMs. No upregulated pathway was discovered in BCBMs compared with BCEMs. Some immune-related pathways were discovered downregulated in BCBMs compared with BCEMs in GSE14017 and GSE14018 (FDR < 0.25, Figure 5F and Supplementary Table S8 ), such as antigen processing and presentation, leishmania infection, jak stat signaling pathway, and nod-like receptor signaling pathway.

Next, we investigated whether increased OXPHOS utilization is functionally important for metastasis or only represents a response to the interplay between metastatic tumor cells and the brain microenvironment. We used oligomycin, an inhibitor of mitochondrial F(1)F(o)ATPase, to inhibit OXPHOS in MDA-MB-231 cells. Schematic diagram of experimental plan for determining the effect of oligomycin treatment on MDA-MB-231 cells was presented in Figure 6A . We first explored the effect of oligomycin on the proliferation and viability of MDA-MB-231cells using cell proliferation curve and flow apoptosis assays. The CCK8 proliferation curve showed that there was no significant difference in proliferation (p > 0.05; Figure 6B ) and growth rate (p > 0.05; Figure 6C ) of MDA-MB-231 cells in the control group and Oligo(1.0) group within 48 h. Flow cytometry analysis of cells stained with Annexin V and PI showed that compared with the control group, Oligo(1.0) did not reduce cell viability or increase cell apoptosis after 48 h (p > 0.05; Figure 6D ). These results were consistent with previous studies that cancer cells can switch between glycolysis and OXPHOS to adapt to the environment (22). Then, we explored the effect of oligomycin on the metastatic potential of MDA-MB-231 cells using scratch, migration, and invasion assays. The scratch assay showed that the healing rate of scratch wounds significantly decreased in the Oligo(1.0) group compared with the control group (p < 0.01; Figure 6E ). Migration and invasion assays showed that the migrated (p < 0.001; Figure 6F ) and invaded cells (p < 0.05; Figure 6F ) significantly decreased in the Oligo(1.0) group compared with the control group. Together, these assays confirmed that MDA-MB-231 cells can switch between glycolysis and OXPHOS to adapt to the environment and had stronger migration and invasion potential in the case of OXPHOS metabolism.