Mitochondria-Related Core Genes in Breast Cancer: Potential Diagnostic, Prognostic Markers and Therapeutic Targets

Background: Mitochondria-nuclear cross talk and mitochondrial retrograde regulation are involved in the genesis and development of breast cancer (BC). Therefore, mitochondria can be regarded as a promising target for BC therapeutic strategies. In the present study, we aimed to seek the potential biomarkers of BC diagnosis and prognosis and also the molecular therapeutic targets from the perspective of mitochondrial dysfunction. Methods: The microarray data of mitochondria-related encoding genes of BC were downloaded from GEO including GSE128610 and GSE72319. GSE128610 was treated as test set and validation sets consisted of GSE72319 and TCGA, which were used for identifying mitochondria-related differential expressed genes (mrDEGs). GO analysis and KEGG mapping of validated mrDEGs was performed by STRING and KEGG mapper, respectively. PPI network was constructed by STRING. MCODE was applied to modeling prediction of PPI and hub mrDEGs were screened by cytohubba. The online tool of Kaplan Meier plotter was adopted to analyze the correlation between hub mrDEGs and the overall survival of BC patients. Results: A total of 23 up-regulated and 71 down-regulated mrDEGs were identied and validated. Enrichment analyses indicated that mrDEGs were associated with several biological processes, including cell migration, cell surface receptor signaling pathway, cell differentiation and cell communication. Meanwhile, these genes were mapped in PI3K-ALT pathway, TGF-beta pathway, evading apoptosis pathway and resistance to chemotherapy pathway. Moreover, 9 hub mrDEGs were identied. Among them, up-regulated hub mrDEGs contained FN1, BGN, EFNA3, COL5A2 and SEMA3F, and down-regulated hub mrDEGs comprised RHOQ, SEMA3A, NRP1 and DDR2. Overall survival analysis suggested that the up-regulated FN1 and SEMA3F and down-regulated DDR2 conferred to poor BC prognosis. Conclusion: The 9 hub mrDEGs, especially FN1, SEMA3F and DDR2, were likely to regulate mitochondrial function and might be novel biomarkers of BC diagnosis and prognosis as well as the therapeutic targets.


Background
Mitochondria, the only extranuclear organelle carried with genetic material, plays an important role in carcinogenesis through its communication and retrograde regulation of nucleus (1). The Reactive Oxygen Species in mitochondria were suggested to promote proliferation, migration and apoptosis of tumor cells (2). The mitochondria in BC cells could exert retrograde regulation of nucleus by transmitting signal to them, facilitating the bidirectional communication between each other (3), and making mitochondria an anticancer drug target for tumor. In addition, mitochondria from non-cancer cell lines has been shown to suppress multiple carcinogenic pathways and reverse the carcinogenic properties of tumor cells under the same nuclear background, including cell proliferation, viability in hypoxia, anti-apoptosis property, resistance to anticancer drugs, invasion, colony formation, and enhancing the response of tumor cells to treatment (4). These ndings emphasize that mitochondria has critical regulatory roles for cancer cell property. The correction of mitochondrial function is a promising target of anticancer therapy (4).
Breast cancer (BC) is the most common cancer and the leading cause of cancer-related death in female (5,6). The exploration of potential biomarkers for early diagnosis and therapeutic targets of BC has important scienti c signi cance and application value. In recent years, it remains rare about the differential analysis of mitochondria-related encoding genes in BC. Hence, the model construction for predicting early BC diagnosis and prognosis via bioinformatics would greatly bene t the identi cation of potential mitochondrial diagnostic biomarkers and therapeutic targets for BC. In the present study, two microarray datasets of mitochondria-related genes in BC were collected from Gene Expression Omnibus (GEO), of which one served as test set and the other served as validation set. Then the mitochondriarelated differential expressed genes(mrDEGs)were screened out and validated by TCGA database. Our study aims to focus on mrDEGs and seek potential diagnostic and prognostic biomarkers as well as the molecular therapeutic targets of BC from the perspective of mitochondrial dysfunction.

Data collection
To identify mrDEGs involved in BC genesis and development, two datasets (GSE128610 and GSE72319) were collected from Gene Expression Omnibus (GEO) (https://www.ncbi.nlm.nih.gov/gds/). GSE128610 contained three BC cell lines (MDA-MB-468) and three BC-free epithelial cell lines (MCF10A). GSE72319 was composed of three triple-negative BC cell lines (SUM159) and three benign BC cell lines (A1N4). Both of them adopted transmitochondrial cybrid system (Cybrid), which is well acknowledged in mitochondrial function research currently. For the Cybrid model, the nucleus in experimental and control groups were both replaced by other cells' to eliminate the interference of nuclear encoding genes in mitochondrial function research (7,8). BC data was downloaded from The Cancer Genome Atlas (TCGA, https://www.cancer.gov/about-nci/organization/ccg/research/structural-genomics/tcga), including 1112 BC tissues and 113 normal breast tissues.

Data processing
In this research, GSE128610 was treated as test set, while GSE72319 and TCGA data respectively served as the rst and second validation set. To identify mrDEGs, the original microarray datasets of GSE128610 and GSE72319 were analyzed using GEO2R, and TCGA data was processed by edgeR and SangerBox.
The screening criteria were set at |logFC|>=1, P < 0.05 and P. adjust < 0.05. GSE72319 and TCGA were successively used to verify mrDEGs in GSE128610 through the "MATCH function".

GO enrichment analysis and KEGG mapping
Gene ontology (GO) analysis was performed for validated mrDEGs by Search Tool for the Retrieval of Interacting Genes (STRING), including cellular component, molecular function and biological process.. P < 0.05 was considered as statistical signi cance. KEGG mapper (https://www.genome.jp/kegg/mapper.html) was applied to Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway map of mrDEGs.
2.4 Protein-protein interaction (PPI) network construction and modeling analysis PPI network of validated mrDEGs was constructed by STRING. The cut-off value of Interaction score was set at 0.4, and PPI network was visualized. Subsequently, the classical models were screened out by Molecular Complex Detection (MCODE) plug-in of Cytoscape_3.7.2 based on the criteria of score > = 3 and nodes > = 3. The function enrichment analysis of single model was performed by STRING, and P < 0.05 was considered as cut-off criteria. In addition, hub mrDEGs were selected by cytohubba plug-in of Cytoscape_3.7.2 according to MCC > = 6.

Survival analysis based on hub mrDEGs
The correlation between hub mrDEGs and the overall survival of 1402 BC cases was analyzed applying Kaplan Meier plotter (http://kmplot.com/analysis/). The hazard ratios with 95% con dence intervals and P values of log rank test were calculated and displayed in the gure. P < 0.05 was regarded to be statistically signi cant.

The identi cation and validation of mrDEGs in BC
The ow diagram was shown in Fig. 1. GSE128610 and GSE72319 datasets were analyzed online by GEO2R. The mrDEGs in BC were identi ed based on the cut-off criteria of |logFC|>=1, P < 0.05 and P. adjust < 0.05. We found out 1756 up-regulated and 3225 down-regulated mrDEGs in GSE128610.
Subsequently, 251 up-regulated and 1162 down-regulated mrDEGs in GSE128610 were validated by GSE72319. After initial validation, they were further veri ed by TCGA data, and then 23 up-regulated and 71 down-regulated mrDEGs were nally identi ed (Table 1).
mrDEGs: mitochondria-related differential expressed genes 3.2 GO function enrichment and KEGG pathway map for mrDEGs GO enrichment analysis of 94 validated mrDEGs was performed by STRING. The results were presented in Table 2 and supplementary Table 1. The mrDEGs in BC were shown to play roles in cancer-related biological processes, such as neural crest cell migration involved in autonomic nervous system development, regulation of cell migration, cell surface receptor signaling pathway and cell differentiation.
However, no signi cant molecular function or cellular component was observed.

mrDEGs-related PPI network construction and Modeling analysis
PPI network of 94 mrDEGs in BC was constructed by STRING database, with a total of 94 nodes and 60 edges (Fig. 3a). Three models were found to meet the criteria of score > = 3 and nodes > = 3 by MOCODE plug-in of Cytoscape software (Fig. 3b). GO and KEGG enrichment analyses were carried out for the three models respectively (Table 3 and Supplementary Table 2). The results showed that model 1 was mainly associated with the structure and function of nerve cells, model 2 was involved in extracellular matrix, and model 3 took parts in tissue development and gene expression.

Discussion
In this study, the microarray data of mitochondria-related genes in BC collected from GEO were used to identify the mrDEGs further validated in TCGA. Moreover, GO enrichment analysis and KEGG pathway mapping for validated mrDEGs were performed to explore the potential function of mrDEGs in breast carcinogenesis. Based on that, we constructed the PPI network, discovered the hub mrDEGs, and analyzed the correlation between hub mrDEGs and the overall survival of BC patients to investigate the in uence of hub mrDEGs on BC prognosis.
Mitochondria play a critical role in multiple cell processes, and mitochondrial dysfunction may affect the occurrence and development of BC. The initiation and metastasis of BC could be altered by regulating the genetic background of mitochondria, making mrDEGs potential therapeutic targets (9). Here, we utilized multiple bioinformatics tools to analyze the microarray data of mitochondria-related genes, and found out 23 up-regulated and 71 down-regulated mrDEGs in BC. They were closely associated with mitochondrial dysfunction in breast carcinogenesis. GO enrichment analysis demonstrated that 94 mrDEGs were enriched in cancer-related biological processes, such as neural crest cell migration involved in autonomic nervous system development, regulation of cell migration, cell surface receptor signaling pathway, cell differentiation and regulation of cell communication. These biological processes conformed to tumor cell properties, including unlimited cell proliferation, cell invasion and migration, and reduced intercellular adhesion (10), suggesting that mrDEGs were tightly linked to breast carcinogenesis. KEGG pathway mapping showed that mrDEGs might participate in cancer-related regulation pathways, including PI3K-ALT pathway, TGF-beta pathway, evading apoptosis and resistance to chemotherapy. The relationship of them with BC could be listed as follows: (1) Inhibiting PI3K-ALT pathway may induce mitochondria-mediated cell apoptosis of BC (11); (2) Ligand-dependent or cell-autonomous activation of the TGF-β pathway in stromal cells could induce metabolic reprogramming, enhance oxidative stress, mitochondrial autophagy and aerobic glycolysis, and decrease Cav-1, which can spread to adjacent broblasts and maintain BC cell growth (12); (3). Regarding the well-known property of unlimited proliferation in BC cells, the GSTs gene mapped in evading apoptosis pathway could regulate cell apoptosis by its interaction with various protein partners (13); (4). MITF, a differential gene identi ed in our research, is able to enhance mitochondrial oxidative phosphorylation (14). It has been reported that enhanced mitochondrial oxidative phosphorylation may induce the resistance to chemotherapy of BC cells (15), thus these genes could be related to drug resistance of BC cells. Overall, the validated mrDEGs mentioned above might be enriched in the pathways for BC progression through regulating mitochondrial function.
PPI network analysis indicated that three interaction networks could be classical models to predict BC occurrence. Model 1 consisted of SEMA3F, EFNA3, SEMA3A and NRP1, which were mainly associated with the structure and function of nerve cells. EFNA3 was induced by HIF under anoxic conditions, and then Ephrin-A3 protein encoded by EFNA3 was aberrantly accumulated to promote the metastasis of BC cells (16). Model 2 was composed of BGN, DDR2 and COL5A2, which were mainly involved in extracellular matrix of cells. COL5A2 related to extracellular matrix remodeling was up-regulated during ductal carcinoma in situ developed to invasive ductal carcinoma, leading to BC progression (17). Model 3 included ZEB1, VIM and FN1, which participated in tissue development and gene expression. ZEB1 increased the expression of vascular endothelial growth factor (VEGF) via paracrine to stimulate angiogenesis in BC (18). ZEB1 also promoted epithelial mesenchymal transformation (EMT), proliferation and migration of BC (19). All the three models with different function took their parts in BC progression.
In our study, 9 hub mrDEGs were screened out based on MCC method, including up-regulated FN1, BGN, EFNA3, COL5A2 and SEMA3F as well as down-regulated RHOQ, SEMA3A, NRP1 and DDR2. Next, Kaplan Meier plotter was used to analyze the association between hub mrDEGs and the overall survival of BC patients. The results showed that only up-regulated FNA and SEMA3F and down-regulated DDR2 suggested poor BC prognosis (P < 0.05) with the potential to be a signi cant biomarker. FN1 has been demonstrated to be up-regulated in BC epithelial cells without mitochondria DNA (20). FN1 was also a core gene of mrDEGs network and its encoded bronection distributed in BC cell matrix affecting tumor progression (20). Meanwhile, FN1 could regulate EMT of BC cells (21) and might be one of the key genes in BC invasion and migration (22). Therefore, FN1 has the potential to be a diagnostic biomarker and a molecular therapeutic target of BC. During the initial of breast neoplasms, SEMA3F could prevent tumor cells from spreading and attaching to stromal cells for further implantation (23). Additionally, SEMA3F may interact with NRP1 or NRP2 receptor to suppress the metastasis and invasion of BC cells (24). Hence, it could be speculated that up-regulated SEMA3F is likely to be a capable biomarker for BC detection. DDR2 was activated by brillar collagen to regulate the synthesis of extracellular matrix and wound healing (25), exerting important roles in microenvironment. DDR2 was involved in hypoxia-induced cancer metastasis by accelerating migration, invasion and EMT of BC cells (26), which might be a potential molecular therapeutic target. Above all, mitochondria-related hub genes may function in different stages of BC. Further investigations on these genes would help to elucidate BC etiology from the perspective of mitochondrial dysfunction, and thus to identify diagnostic and prognostic biomarkers and also molecular targets for BC targeted therapy.

Conclusions
In summary, we employed bioinformatics analyses to discover mrDEGs in BC. Then enrichment analyses for these genes were carried out and three interaction networks were constructed to serve as classical models for predicting breast carcinogenesis. We also selected 9 hub mrDEGs including FN1, BGN, EFNA3, COL5A2, SEMA3F, RHOQ, SEMA3A, NRP1 and DDR2. Among them, FN1, SEMA3F and DDR2 were closely associated with BC prognosis with the potential to be prognostic biomarkers and therapeutic targets.

Declarations
Ethics approval and consent to participate Not applicable Consent for publication Not applicable Availability of data and materials The data that support the results of this manuscript are available from the corresponding author upon reasonable request.

Competing interests
All authors disclose no con icts of interest that might bias their work.

Funding
This work is supported by the National Key R&D Program of China (Grant #2018YFC1311600) and partly by grants from the Natural Science Foundation of Liaoning Province in China (201705040987).

Authors' contributions
Yuan Yuan and Qian Xu conceived and designed this study. Li-rong Yan and Ang Wang were responsible for the data analysis and performed data interpretation. Li-rong Yan wrote the paper. Qian Xu, Zhi Lv and Yuan Yuan revised the manuscript. Flow diagram of bioinformatics analysis. mrDEGs: mitochondria-related differential expressed genes.