Mitochondria-related core genes and TF-miRNA-hub mrDEGs network in breast cancer

Abstract 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. The present study aimed to construct regulatory network and seek the potential biomarkers of BC diagnosis and prognosis as well as the molecular therapeutic targets from the perspective of mitochondrial dysfunction. Methods: The microarray data of mitochondria-related encoding genes in BC cell lines were downloaded from GEO including GSE128610 and GSE72319. GSE128610 was treated as test set and validation sets consisted of GSE72319 and TCGA tissue samples, intending to identify mitochondria-related differentially expressed genes (mrDEGs). We performed enrichment analysis, PPI network, hub mrDEGs and overall survival analysis and constructed transcription factor (TF)-miRNA-hub mrDEGs network. Results: A total of 23 up-regulated and 71 down-regulated mrDEGs were identified and validated in BC cell lines and tissues. Enrichment analyses indicated that mrDEGs were associated with several cancer-related biological processes. Moreover, 9 hub mrDEGs were identified and validated in BC cell lines and tissues. Finally, 5 hub coregulated mrDEGs, 21 miRNAs and 117 TFs were used to construct TF-miRNA-hub mrDEGs network. MYC associated zinc finger protein (MAZ), heparin binding growth factor (HDGF) and Sp2 transcription factor (SP2) regulated 3 hub mrDEGs. Hsa-mir-21-5p, hsa-mir-1-3p, hsa-mir-218-5p, hsa-mir-26a-5p and hsa-mir-335-5p regulated 2 hub mrDEGs. Overall survival analysis suggested that the up-regulation of fibronectin 1 (FN1), as well as the down-regulation of discoidin domain receptor tyrosine kinase 2 (DDR2) correlated with unfavorable prognosis in BC. Conclusion: TF-miRNA-hub mrDEGs had instruction significance for the exploration of BC etiology. The hub mrDEGs such as FN1 and DDR2 were likely to regulate mitochondrial function and be novel biomarkers for BC diagnosis and prognosis as well as the therapeutic targets.


Introduction
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 breast cancer (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 noncancer 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 therapy [4]. These findings emphasized that mitochondria had critical regulatory roles in cancer cell property. Some reports showed that mitochondrial-related gene expression might lead to human pathogenesis and the correction with mitochondrial function was a promising target for anticancer therapy [4][5][6].
MiRNAs could participate in the whole signal pathways of tumor genesis and progression, including the regulation of mitochondrial function [7,8]. Moreover, miRNAs are also involved in the regulation of Otto Warburg effect, thus affecting tumor progression [8]. Transcription factors (TFs) are regulatory factors at transcriptional level for the progression of breast cancer [9,10], and modulate mitochondria biogenesis and mitochondria-to-nuclear communication [11,12]. The transcription of mRNAs and miRNAs are regulated by transcription factors (TFs) and the expression of TFs is modulated by miRNAs [13], both of which are closely related to mitochondria function. Therefore, it is of great importance to construct regulatory network for 'TF-miRNA-hub mrDEGs' to explore the mitochondria dysfunction of BC.
Breast cancer (BC) is the most common malignancy and the leading cause of cancer-related death in female [14,15]. The exploration of potential biomarkers and regulatory mechanisms for early diagnosis and therapeutic targets of BC has important scientific significance and application values. In recent years, study remains rare about the differential analysis and network regulation mechanism of mitochondria-related encoding genes in BC. Hence, the model and network construction for predicting early BC diagnosis and prognosis by means of bioinformatics would greatly benefit the identification of potential mitochondrial diagnostic biomarkers, therapeutic targets and pathogenic mechanism for BC. In the present study, two microarray datasets of mitochondria-related genes in BC cell lines were collected from Gene Expression Omnibus (GEO), of which one served as test set and the other served as validation set. Then, the mitochondria-related differentially expressed genes (mrDEGs) were screened out and validated with tissue samples in The Cancer Genome Atlas (TCGA) database. Our study aimed to focus on mrDEGs, construct potential TF-miRNA-mrDEGs network and seek potential diagnostic and prognostic biomarkers as well as the molecular therapeutic targets for 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 samples of BC cell lines (MDA-MB-468) and three BC-free samples of BC-free epithelial cell lines (MCF10A). GSE72319 is composed of three BC samples of triple-negative BC cell lines (SUM159) and three BC-free samples of benign BC cell lines (A1N4). Both of them adopted transmitochondrial cybrid system (Cybrid), which was well acknowledged in current mitochondrial function research. For 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 [16,17]. BC data were downloaded from The Cancer Genome Atlas (TCGA, https://www.cancer.gov/tcga), including 1112 BC tissues and 113 normal breast tissues [18].

Data processing
In this research, GSE128610 was treated as test set, while GSE72319 and TCGA data respectively served as the first and second validation set. To identify mrDEGs, original microarray datasets of GSE128610 and GSE72319 were analyzed with GEO2R, and TCGA data were processed with edgeR and SangerBox. GSE72319 and TCGA were successively employed to verify mrDEGs in GSE128610 through 'MATCH function' . The screening and validation criteria were set as |logFC|>= 1, P<0.05 and P. adjust<0.05. FunRich 3.1.3 was used to plot the heatmap to visualize the validated mrDEGs in GSE128610.

The identification and validation of mrDEGs in BC
The flow diagram was shown in Figure 1. GSE128610 and GSE72319 datasets of BC cell lines were analyzed online with GEO2R. MrDEGs in BC were identified 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 with GSE72319. After initial validation, they were further verified with tissue samples in TCGA database, and then 23 up-regulated and 71 down-regulated mrDEGs were finally identified (Table 1 and Figure 2). The criteria for both validation in GSE72319 and TCGA were the same with GSE128610 (|logFC| ≥ 1, P<0.05 and P. adjust<0.05).

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 ( Figure  4A). Three models were found to meet the criteria of score ≥ 3 and nodes ≥ 3 by MOCODE plug-in of Cytoscape software ( Figure 4B). GO and KEGG enrichment analyses were carried out for the three models, respectively (Table 3  and Supplementary Table S2). 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.

Selection and tissue validation of hub mrDEGs
We screened out 9 hub mrDEGs according to MCC ≥ 6 by cytoHubba plug-in of Cytoscape. Among them, up-regulated hub mrDEGs contained FN1, BGN, ephrin A3 (EFNA3), collagen type V alpha 2 chain (COL5A2) and SEMA3F, and down-regulated hub mrDEGs comprised ras homolog family member Q (RHOQ), SEMA3A, NRP1 and DDR2 (Table 4). Consistent results were obtained after validating and visualizing the expression of 9 hub mrDEGs in Ualcan database with tissue samples from TCGA ( Figure 5).

Analysis of TF-hub mrDEGs network
In order to explore the potential regulatory relationship of hub mrDEGs, we predicted the TFs targeting hub mrDEGs with ENCODE database. The result demonstrated that 8 hub mrDEGs were matched except for COL5A2. The Cytoscape software was applied to visualize the TF-hub mrDEGs network, and 167 links among 8 hub mrDEGs and 121 TFs were predicted ( Figure 6). MAZ could regulate 5 hub mrDEGs (e.g. EFNA3, SEMA3A and SEMA3F). SP2, MAX

Discussion
In the present study, the microarray data of mitochondria-related genes in BC cell lines collected from GEO were utilized to identify mrDEGs and further validated in TCGA tissue samples. 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 PPI network, discovered and validated the hub mrDEGs. Furthermore, TF-miRNA-hub mrDEGs network was constructed and correlation of coregulated hub mrDEGs with the overall survival of BC patients was analyzed to investigate the influence of coregulated hub mrDEGs on BC prognosis. Mitochondria function plays 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 [28]. Here, multiple bioinformatics tools were adopted to analyze the microarray data of mitochondria-related genes, and found out 23 up-regulated and 71 down-regulated mrDEGs in BC lines and tissue samples. 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 [29], 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. Their relationship with BC could be listed as follows: (1) Inhibiting PI3K-AKT pathway may induce mitochondria-mediated cell apoptosis of BC [30]; (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, Table 5 Hub mrDEGs of coregulated by miRNAs and TFs   TF  Hub mrDEGs  Gene counts  miRNA  Hub mrDEGs  Gene counts   MAZ  EFNA3, NRP1, RHOQ  3  hsa-mir-218-5p  NRP1, FN1  2   HDGF  DDR2, EFNA3, RHOQ  3  hsa-mir-26a-5p  NRP1, FN1  2   SP2  EFNA3, NRP1, RHOQ  3  hsa-mir-335-5p  NRP1,  Abbreviations: mrDEGs, mitochondria-related differential expressed genes; TFs, transcription factors. and decrease Cav-1, which can spread to adjacent fibroblasts and maintain BC cell growth [31]; (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 [32]; (4) Melanocyte inducing transcription factor (MITF), a differential gene identified in our research, is able to enhance mitochondrial oxidative phosphorylation [33]. It has been reported that enhanced mitochondrial oxidative phosphorylation may induce the resistance to chemotherapy of BC cells [34]; thus, these genes could be related to drug resistance of BC cells. Overall, the validated mrDEGs mentioned above might be enriched in the pathways of 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 [35]. 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 the development from ductal carcinoma in situ to invasive ductal carcinoma, leading to BC progression [36]. Model 3 included zinc finger E-box binding homeobox 1 (ZEB1), vimentin (VIM) and FN1, participating in tissue development and gene expression. ZEB1 increased the expression of vascular endothelial growth factor (VEGF) via paracrine to stimulate angiogenesis in BC [37]. ZEB1 also promoted epithelial-mesenchymal transformation (EMT), proliferation and migration of BC [38]. 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. Ualcan was utilized to validate and visualize the expression of 9 hub mrDEGs. TF-miRNA-hub mrDEGs network was constructed, which had important instruction significance to explore the potential regulatory mechanism of hub mrDEGs in BC. TF-miRNA-hub mrDEGs network showed that MAZ of TF nodes could interact with 3 hub mrDEGs including EFNA3, NRP1 and RHOQ, which implied its significance in BC. Myc-associated zinc finger protein (MAZ) has been considered as a transcription factor with C2H2 zinc finger motif binding to a GA box [39,40], with important roles in BC progression [40,41]. The study suggested that the transactivation and transcriptional alteration of MAZ could modulate the process of aerobic glycolysis in tumor [42]. NRP1 targeting hub mrDEGs of MAZ might be located in mitochondria regulating mitochondrial function and iron-dependent oxidative stress [43]. A GEO dataset (GSE115118) for miRNA mitochondrial sublocalization indicated that hsa-mir-218-5p, hsa-mir-26a-5p and hsa-mir-335-5p regulated by NRP1 were located in mitochondria. Therefore, MAZ was predicted to impact on mitochondria function by interacting with NRP1 and these miRNAs to regulate BC progression. Further investigation is needed to elucidate the function of these hub mrDEGs in BC.
Survival analysis showed that the up-regulated FN1 and down-regulated DDR2 suggested poor BC prognosis (P<0.05) with the potential to be a significant biomarker. FN1 has been demonstrated to be up-regulated in BC epithelial cells without mitochondria DNA [44]. FN1 was also a core gene of mrDEGs network and its encoded fibronection distributed in BC cell matrix affecting tumor progression [44]. It has been reported that FN1 can increase cisplatin resistance in non-small cell lung cancer by modulating β-catenin signaling via interaction with integrin-β1 [45]. FN1 might also be a potential biomarker for radioresistance in head and neck squamous cell carcinoma [46] and a potential therapeutic target for breast cancer [47]. We initially reported that FN1 was closely related to mitochondrial function. It had the highest degree score and could also interact with hsa-mir-218-5p and hsa-mir-26a-5p sublocated in mitochondria (GSE115118). Mitochondria are well accepted to exert crucial roles for anticancer drug target in tumor. Therefore, we could speculate that FN1 might be an anticancer drug target for BC by regulating mitochondrial dysfunction. DDR2 was activated by fibrillar collagen to regulate the synthesis of extracellular matrix and wound healing affecting microenvironment [48]. DDR2 was involved in hypoxia-induced cancer metastasis by accelerating migration, invasion and EMT of BC cells [49]. Our TF-miRNA-hub mrDEGs network showed that DDR2 could interact with hsa-mir-21-5p, hsa-mir-548a-3p and hsa-mir-129-2-3p, which were also sublocated in mitochondria (GSE115118). Further study for these genes would help to elucidate BC etiology from the perspective of mitochondrial dysfunction, and thus to identify diagnostic and prognostic biomarkers as well as molecular targets for BC targeted therapy.
In summary, bioinformatics analyses were employed to discover mrDEGs in BC. Gene enrichment analyses were carried out and three interaction networks were constructed to serve as classical models for predicting breast carcinogenesis. We also selected 5 coregulated hub mrDEGs by TFs and miRNAs including FN1, DDR2, NRP1, EFNA3 and RHOQ. And TF-miRNA-hub mrDEGs network was constructed to explore the potential pathogenesis of hub mrDEGs in BC.

Data Availability
The data that support the results of this manuscript are available from the corresponding author upon reasonable request.