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.
Mitochondria, the only extranuclear organelle carried with genetic material, plays an important role in carcinogenesis through its communication and retrograde regulation of nucleus . The reactive oxygen species in mitochondria were suggested to promote proliferation, migration and apoptosis of tumor cells . 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  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 . 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–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 . 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 , 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.
Materials and methods
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 .
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.
GO enrichment analysis and KEGG mapping
Gene ontology (GO) analysis was performed for validated mrDEGs with Search Tool for the Retrieval of Interacting Genes (STRING, https://string-db.org/), including cellular component, molecular function and biological process . P<0.05 was considered as statistical significance. KEGG mapper (https://www.genome.jp/kegg/mapper.html) was applied to Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway map of mrDEGs.
Protein–protein interaction (PPI) network construction and modeling analysis
PPI network of validated mrDEGs was constructed with STRING. The cut-off value of Interaction score was set as 0.4, and PPI network was visualized. Subsequently, 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 regarded as cut-off criteria.
Screening and tissue identification of hub mrDEGs
Hub mrDEGs were selected by cytohubba plug-in of Cytoscape_3.7.2 according to MCC ≥ 6 . Next, Ualcan (http://ualcan.path.uab.edu/) was utilized to validate the expression levels of mrDEGs in BC tissue samples .
Construction of TF-hub mrDEGs network and miRNA-hub mrDEGs network
NetworkAnalyst (https://www.networkanalyst.ca/), a website for comprehensive gene expression analysis, meta-analysis and network biology, was applied to search TFs and miRNAs targeting hub mrDEGs [22–24]. Hub mrDEGs were uploaded to NetworkAnalyst to acquire the TFs targeting hub mrDEGs from ENCODE database and miRNA-hub mrDEGs pairs from TarBase and miRTarBase. Two interaction lists were downloaded and Cytoscape_3.7.2 was used to visualize TF-hub mrDEGs network and miRNA-hub mrDEGs network respectively.
TF-miRNA-hub mrDEGs network construction and validation
In order to construct the TF-miRNA-hub mrDEGs network, TF-hub mrDEGs network and miRNA-hub mrDEGs network were overlapped by ‘MATCH function’. Coregulated hub mrDEGs of TFs and miRNAs were selected to construct TF-miRNA-hub mrDEGs network and visualized with Cytoscape_3.7.2.
Subsequently, we searched miRNA-hub mrDEGs network in starBase database (http://starbase.sysu.edu.cn/starbase2/mirMrna.php) and TF-hub mrEGs network in ChEA database [25,26] to screen overlapped hub mrDEGs between the two networks and validate TF-miRNA-hub mrDEGs network. The validated TF-miRNA-hub mrDEGs network was also visualized with Cytoscape_3.7.2.
Survival analysis based on coregulated hub mrDEGs by TFs and miRNAs
The correlation between coregulated hub mrDEGs and the overall survival of 1402 BC cases was analyzed by using the Kaplan–Meier plotter database (www.kmplot.com) . The mRNA gene chip of BC was chosen in this database. All BC patients were divided into two groups according to the median of gene expression to perform survival curve analysis. The hazard ratios with 95% confidence intervals and P values of log-rank test were calculated and displayed in the figure. P<0.05 was considered to be statistically significant.
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).
GO function enrichment and KEGG pathway map for mrDEGs
GO enrichment analysis for 94 validated mrDEGs was performed by STRING. The results were presented in Table 2 and Supplementary Table S1. The mrDEGs in BC were shown to function in cancer-related biological processes, such as neural crest cell migration involved in autonomic nervous system development (e.g. FN1, semaphorin 3A (SEMA3A), semaphorin 3A (SEMA3F) and neuropilin 1 (NRP1)), regulation of cell migration (e.g. C-X-C motif chemokine ligand 2 (CXCL2), DDR2 and FN1), cell surface receptor signaling pathway (e.g. Microtubule Affinity Regulating Kinase 1 (MARK1), biglycan (BGN) and CXCL2) and cell differentiation (e.g. Ankyrin 2 (ANK2), Glypican 2 (GPC2) and melanocyte inducing transcription factor (MITF)). However, no significant molecular function or cellular component was observed.
The mrDEGs were mapped in KEGG pathway. They were suggested to participate in the following cancer-related pathways: PI3K-AKT pathway (e.g. phosphoinositide-3-kinase regulatory subunit 1 (PIK3R1)), TGF-beta pathway (e.g. MITF), evading apoptosis (e.g. microsomal glutathione S-transferase 1 (MGST1)) and resistance to chemotherapy (e.g. MITF) (Figure 3).
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 dimerization protein 4 (MXD4), Kruppel-like factor 9 (KLF9), Kruppel-like factor (KLF16), HDGF and AT-rich interaction domain 4B (ARID4B) regulated 3 hub mrDEGs.
Analysis of miRNA-hub mrDEGs network
Subsequently, miRNA-hub mrDEGs pairs of 9 hub mrDEGs were performed with TarBase and miRTarBase. Finally, only 6 hub mrDEGs were mapped except for BGN, SEMA3F and SEMA3A, then 31 links among 25 miRNAs and 6 hub mrDEGs were obtained by Cytoscape software (Figure 7). We found that hsa-mir-21-5p could interact with 3 hub mrDEGs including COL5A2, DDR2 and RHOQ.
Construction and validation of TF-miRNA-hub mrDEGs network analysis
The hub mrDEGs (FN1, EFNA3, NRP1, RHOQ and DDR2) coregulated by TFs and miRNAs were selected and their interactive regulators were extracted. Then, TF-miRNA-hub mrDEGs network was constructed with Cytoscape (Figure 8). A total of 5 hub mrDEGs, 21 miRNAs and 117 TFs were included in the TF-miRNA-hub mrDEGs network. Next, we analyzed the interactive results of TF-mrDEGs and miRNA-hub network, respectively (Table 5). We found that MAZ, HDGF and SP2 could regulate 3 hub mrDEGs. Simultaneously, 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. In addition, FN1, EFNA3 and NRP1 had the highest degree in the network.
Subsequently, validation of TF-miRNA-hub mrDEGs network was performed by starBase and ChEA. It was shown that 3 hub mrDEGs (FN1, EFNA3 and NRP1), 10 miRNAs and 6 TFs were validated and visualized with Cytoscape software (Supplementary Figure S1).
Survival analysis of coregulated mrDEGs by TFs and miRNAs
Prognostic significance of the 5 coregulated hub mrDEGs by TFs and miRNAs in BC was evaluated using Kaplan–Meier plotter. Two hub mrDEGs with statistical significance were identified in survival analysis (P<0.05, n=1402) including FN1 (HR = 1.28 (1.03–1.59), P=0.023) and DDR2 (HR = 0.77 (0.62–0.96), P=0.017). Therefore, the up-regulated FN1 and down-regulated DDR2 might confer to poor BC prognosis (Figure 9).
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 . 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 , 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 ; (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 fibroblasts and maintain BC cell growth ; (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 ; (4) Melanocyte inducing transcription factor (MITF), a differential gene identified in our research, is able to enhance mitochondrial oxidative phosphorylation . It has been reported that enhanced mitochondrial oxidative phosphorylation may induce the resistance to chemotherapy of BC cells ; 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 . 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 . 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 . ZEB1 also promoted epithelial–mesenchymal transformation (EMT), proliferation and migration of BC . 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 . NRP1 targeting hub mrDEGs of MAZ might be located in mitochondria regulating mitochondrial function and iron-dependent oxidative stress . 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 . FN1 was also a core gene of mrDEGs network and its encoded fibronection distributed in BC cell matrix affecting tumor progression . 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 . FN1 might also be a potential biomarker for radioresistance in head and neck squamous cell carcinoma  and a potential therapeutic target for breast cancer . 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 . DDR2 was involved in hypoxia-induced cancer metastasis by accelerating migration, invasion and EMT of BC cells . 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.
The data that support the results of this manuscript are available from the corresponding author upon reasonable request.
The authors declare that there are no competing interests associated with the manuscript.
This work is supported by the National Key R&D Program of China [grant number 2018YFC1311600] and partly by grants from the Natural Science Foundation of Liaoning Province in China [grant number 201705040987].
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.
AT-rich interaction domain 4B
collagen type V alpha 2 chain
C-X-C motif chemokine ligand 2
discoidin domain receptor tyrosine kinase 2
heparin binding growth factor
Kyoto Encyclopedia of Genes and Genomes
Kruppel-like factor 16
Kruppel-like factor 9
microtubule affinity regulating kinase 1
MYC associated zinc finger protein
microsomal glutathione S-transferase 1
melanocyte inducing transcription factor
mitochondria-related differential expressed genes
MAX dimerization protein 4
phosphoinositide-3-kinase regulatory subunit 1
ras homolog family member Q
Sp2 transcription factor
zinc finger E-box binding homeobox 1