A TIMM17A Regulatory Network Contributing to Breast Cancer

Background Translocase of inner mitochondrial membrane 17A (TIMM17A) is overexpressed in breast cancer (BRCA), and upregulation can increase the aggressiveness of BRCA cells. This study examined the influence of the TIMM17A gene network on BRCA outcome. Methods Expression levels of TIMM17A were compared between normal and tumor tissues from the OncomineTM database, and the association with patient survival was analyzed using Kaplan–Meier Plotter. Clinical factors influencing TIMM17A expression were studied by UALCAN. cBioPotal was then used to identify genes interacting with TIMM17A, and network relationships were assessed using the R clusterProfiler package. The association between TIMM17A mutation and mRNA expression in BRCA was examined using the LinkFinder application in LinkedOmics, and coexpressed genes were assessed for functional enrichment using the LinkInterpreter application. Furthermore, TIMM17A expression correlation with cell cycle phase distribution was performed by flow cytometry. Finally, the target networks of kinases, microRNAs (miRNAs), and transcription factors were identified using GeneMANIA. The expression and correlation of potential miRNAs and targets were further validated in BRCA cell lines by qRT-PCR. Results Expression of TIMM17A was significantly elevated in BRCA compared with normal tissue (p < 0.05), and overexpression was associated with both poor overall survival (OS) and shorter distant metastasis-free survival (DMFS) (p < 0.05). Expression of TIMM17A was not associated with age, sex, BRCA subclass, clinical stage, or patient ethnicity. The coexpressed TIMM17A network was enriched in genes targeted by cell cycle regulators such as CDK1, miR-331, and E2F family transcription factors (FDR < 0.001). Furthermore, flow cytometry revealed a strong association between higher TIMM17A expression and faster cell cycle progression in these BRCA cell lines. In addition, expression of TIMM17A protein was correlated with CDK1 protein expression in BRCA cell lines as measured by western blotting. Conclusion Elevated TIMM17A expression accelerates the progression of BRCA, thereby reducing OS and DMFS. The TIMM17A-associated networks identified here provide clues to the molecular pathogenesis of BRCA and potential targets for BRCA treatment.


INTRODUCTION
According to recent global estimates, breast cancer (BRCA) is now the most frequent cause of cancer-related mortality among adult females (Bray et al., 2018). However, BRCA is a highly heterogeneous clinical entity, with multiple subtypes according to distinct histological features and treatment sensitivity profiles (Harbeck and Gnant, 2017;Yeo and Guan, 2017;Liang et al., 2020). This heterogeneity has also been confirmed at the gene expression level by high-throughput molecular profiling (Ahn et al., 2020;Krug et al., 2020). More precise delineation of these BRCA subtypes may eventually yield more efficacious individualized treatments (Zardavas et al., 2015). Currently, however, therapeutic options for many cancers are limited by our relatively poor understanding of the underlying molecular and cellular mechanisms (Garrido-Castro and Winer, 2018). Malignancy results from the disruption of multiple coordinated gene networks (Garrido-Castro and Winer, 2018), so elucidation of the gene network changes characterizing various BRCA subtypes may facilitate the design of patient-specific interventions (Garrido-Castro and Winer, 2018;Radvanyi, 2018;Zacharakis et al., 2018).
Translocase of inner mitochondrial membrane 17A (TIMM17A) is an essential component of the highly conserved TIM23 translocase complex involved in mitochondrial protein import and metabolism (Bomer et al., 1996). Expression of TIMM17A is elevated in BRCA (Kannangai et al., 2007;Yang et al., 2016). Furthermore, TIMM17A upregulation can increase the aggressiveness of BRCA cells (Yang et al., 2016) and is associated with poor clinical outcome in human BRCA patients (Salhab et al., 2012). Indeed, Xu at el. identified TIMM17A as a biomarker for poor outcome from BRCA (Xu et al., 2010;Salhab et al., 2012).
All genes operate within multiple networks of coregulated and functionally related genes. However, TIMM17A networks are not well described. In this study, we examined the associations between clinical outcome measures and TIMM17A expression levels among BRCA cases from The Cancer Genome Atlas (TCGA). Furthermore, we constructed a network of genes coregulated with TIMM17A in BRCA patients to reveal new potential molecular targets for diagnosis and treatment.

Statistical Analysis
Statistical analyses were performed using SPSS Software version 26.0 (IBM SPSS Statistics, Chicago, IL, United States) and GraphPad Prism 7 (GraphPad Software, San Diego, CA, United States). Overall survival (OS) was compared between BRCA patients with high or low TIMM17A expression by the Kaplan-Meier method and log-rank test (Mishra et al., 2019;Okano et al., 2019). The associations between TIMM17A expression levels and various clinicodemographic features were analyzed by Spearman's correlations. Group means were compared by independent samples t-test (Chandrashekar et al., 2017). The association between TIMM17A mutation and mRNA expression in the TCGA BRCA cohort was assessed by t-test with correction for multiple comparisons (Vasaikar et al., 2018). p-Value less than 0.05 (two-tailed) was considered significant for all tests.

Oncomine Analysis
Differences in TIMM17A expression and copy number between normal tissue and BRCA were evaluated for TCGA and Curtis datasets in the Oncomine database 1 .

Kaplan-Meier Plotter Analysis
The association of tumor TIMM17A expression level with BRCA outcome was examined in total BRCA datasets, GSE45255, and GSE7390 using Kaplan-Meier Plotter (publicly available at https://kmplot.com/). For this analysis, we chose TIMM17A gene chip results (Affymetrix ID 215171_s_at) for estimation of mRNA expression level, both OS and distant metastasis-free survival (DMFS) as clinical outcome metrics, and a follow-up duration of 120 months.

UALCAN Analysis
Clinical factors influencing TIMM17A expression were examined in the TCGA dataset using the "breast invasive carcinoma" section of UALCAN 2 . The influencing factors assessed were age, sex, tumor stage, BRCA subclass, and metastasis occurrence based on a previous factor analysis (Chandrashekar et al., 2017).

c-BioPortal Analysis
Genes interacting with TIMM17A were identified in the "breast invasive carcinoma" dataset of Firehose Legacy using c-BioPortal 3 . Neighboring genes with alteration frequencies greater than 5% were selected for network analysis.

clusterProfiler Analysis
The R 3.6.0 package cluster Profiler was used to visualize networks of genes related to functional changes identified by Disease Ontology (DO), Gene Ontology (GO), Kyoto   Boxplot showing relative expression of TIMM17A in normal individuals and separately for BRCA patients with regional lymph node metastasis N0, N1, N2, or N3 tumors. (F) Boxplot showing relative expression of TIMM17A in normal individuals and separately for BRCA patients with regional lymph node metastasis N0, N1, N2 or N3 tumors. Data are mean ± SE. *P < 0.05; **P < 0.01; ***P < 0.001.
Encyclopedia of Genes and Genomes (KEGG), and Reactome (Hayes, 2010). DO was used to identify genes associated with BRCA, GO to describe gene functions using standardized terminology, and KEGG to place genes into the following molecular interaction and reaction networks: metabolism, genetic information processing, environmental information processing, cellular processes, biological system, human diseases, and drug development. Reactome was used to supplement the results of KEGG analysis by visualizing various human reactions and biological channels. The Benjamini-Hochberg (BH) adjusted p method was used to determine significant enrichment at a cutoff of 0.05.

LinkedOmics Analysis
The RNAseq dataset of the "breast invasive carcinoma" cohort was obtained from LinkedOmics 4 (Vasaikar et al., 2018). The dataset including TIMM17A expression comprised 1,093 patients. LinkFinder was used to identify the association between TIMM17A mutation and mRNA expression in the TCGA BRCA cohort and Pearson's correlation analysis to assess the strengths of the associations (positive and negative). Kinase and KEGG pathway analyses were then performed using the gene set enrichment analysis (GSEA) tool in LinkInterpreter. 4 http://www.linkedomics.org/

KEGG Pathway Analysis
KEGG is composed of 15 comprehensive database resources, a manually planned database, and four system information category databases generated by calculation. The KEGG pathway system can link gene expression to multiple "cells and biological functions" domains, including "metabolism, " "cellular processes, " "biological functions, " and "human diseases" (Kanehisa et al., 2017). We identified the pathways for 26 enriched neighboring genes by KEGG 5 .

Cell Culture
The human BRCA cell lines BT-549 and SK-BR-3 were obtained from the Shanghai Cell Bank Type Culture Collection Committee (Shanghai, China). Cells were cultured in complete growth medium as recommended by the distributor in a humidified incubator at 37 • C under a 5% CO 2 atmosphere.

Cell Transfection
Expression of TIMM17A in cancer cell lines was manipulated by transfection with small interfering RNAs (siRNAs). A negative control siRNA (si-NC) and TIMM17A-targeting siRNA (si-TIMM17A) were designed and synthesized by LAIDEMENG (Guangzhou, China). Cell lines in logarithmic growth phase were seeded at 6 × 10 5 cells/well in six-well plates, allowed to adhere overnight, and then incubated in transfection medium containing the indicated siRNA and Lipofectamine 3000 (Invitrogen, Waltham, MA, United States) according to the standard protocol of the manufacturer. After 6 h, the transfection medium was exchanged for complete medium and cells incubated for another 48 h to allow expression changes before subsequent experiments.

Quantitative RT-PCR
Total RNA was extracted from cells using Trizol reagent (Invitrogen, CA, United States) according to the instructions of the manufacturer, and cDNAs generated using EvoM-MLV RT Premix. Gene expression levels were then measured by quantitative PCR on a Bio-Rad CFX RT-qPCR detection system using SYBR Premix Ex TaqII (Takara, Dalian, China). The 18S rRNA gene was used as an internal control. All samples were run in triplicate. Fold changes in gene expression were calculated using the Ct method. Details of the primer sequences used for qPCR are listed in Supplementary  Table 8.

Western Blotting
Whole cell extracts were prepared by incubation of harvested cells in ice-cold lysis buffer containing 1:100 PMSF and 1:100 protease inhibitor cocktail for 30 min. Lysates were then centrifuged at 14,000 rpm for 5 min at 4 • C, and the supernatant was collected and stored at −80 • C. Lysates were mixed with 5 × THE buffer at 4:1 (vol/vol), boiled for 10 min, cooled slowly to room temperature, centrifuged briefly, and stored at 20 • C. Proteins were separated by sodium dodecyl sulfate polyacrylamide gel electrophoresis and electro-transferred (100 V at low temperature, 1 min/kDa) to polyvinylidene fluoride (PVDF) membranes pretreated with methanol for 5∼30 s, and transmembrane buffer for 20 min. Blotted membranes were then rinsed with Tris-buffered saline containing Triton-X (TBST) for 5 min, incubated in TBST with 5% skimmed milk powder rinsed in TBST for 8 min, and incubated in anti-TIMM17A, anti-CDK1, or anti-GAPDH (gel-loading control) at 4 • C overnight. Immunoblotted membranes were then incubated with the indicated secondary antibody at 37 • C for 50 min to 3 h, and washed three times in TBST (8 min/wash). Protein bands were visualized with chemiluminescence substrate (1∼5 min) and captured using MiniChemi imaging system (Bio-OI, Beijing, China).

Flow Cytometry
The cell cycle progression of cancer cell lines was evaluated using a PI cell cycle Kit. Cells were harvested by trypsin digestion, washed twice in PBS, and fixed in 75% precooled ethanol for more than 2 h. Fixed cells were then treated with kit reagents according to the manufacturer's protocol and analyzed using FACSCantoll flow cytometer (BD FACSCanto II, China).

TIMM17A Expression in BRCA
We obtained multiple TIMM17A expression datasets from Oncomine 4.5 databases (Rhodes et al., 2004). Both mRNA expression level and DNA copy number variation (CNV) were significantly higher in BRCA tissues than healthy breast tissues (p < 0.05) (Figures 1A-F). Although the fold changes were less than 2, TIMM17A ranked within the top 15% of differentially expressed genes based on mRNA abundance and within the top 1% based on DNA CNV (Figures 1A-F). Kaplan-Meier survival analysis using Kaplan-Meier Plotter revealed that higher TIMM17A expression was associated with poorer OS and shorter distant metastases-free survival (DMFS) among 1,402 BRCA cases from the Kaplan-Meier Plotter database (p < 0.05), and similar results were observed for GSE45255 and GSE7390 cohorts (Figures 2A-F). Among these patients, 1,402 died due to BRCA within 120 months, and patients with high expression of TIMM17A demonstrated worse clinical outcome. We also evaluated TIMM17A transcription levels in the Oncomine database, which again revealed that TIMM17A expression was higher in BRCA patients than healthy controls. In contrast to cancer status, TIMM17A expression was not associated with sex, age, ethnicity, tumor stage, tumor grade, or BRCA subtype in the UALCAN database (Figures 3A-F). Thus, TIMM17A expression level may serve as a prognostic indicator of BRCA.    The TCGA Firehose Legacy database from cBioPortal was used to assess the types and frequencies of TIMM17A alterations in BRCA (Gao et al., 2013). Among the 1,093 samples, TIMM17A was altered in 153 cases (14%) (Figure 4A). The most common alteration was amplification (111 cases, 10.16%), followed by multiple alterations (13 cases, 1.19%), and mutation only (two cases, 0.18%), while there were no samples showing mRNA upregulation (Table 1).

Biological Interaction Networks Associated With TIMM17A Alterations in BRCA
To construct a biological interaction network for TIMM17A in BRCA, we obtained neighbor genes with expression alterations at frequencies >5% using the Network function of cBioPortal ( Figure 4B and Table 1), which retrieves networks of queried genes from public pathway databases, such as the Human Reference Protein Database, Reactome, and Memorial Sloan-Kettering Cancer Center Cancer Cell Map (Tung et al., 2016). Neighbor genes with the highest alteration frequencies were PAM16 (8.78%), ATP5MC1 (7.5%), and TAZ (7.04%). The 50 most frequently altered neighbor genes were then subjected to GO analysis, which indicated that this group was enriched in genes related to "mitochondrial transport" and "mitochondrial protein import" (Figures 5A-D).

GO and KEGG Pathway Analysis of mRNAs Co-regulated With Mutant TIMM17A in BRCA
We then identified mRNAs demonstrating correlated expression with mutant TIMM17A in BRCA from TCGA using the Function module of LinkedOmics (Vasaikar et al., 2018). The volcano plot ( Figure 6A) indicated that 7,811 genes (dark red dots) demonstrated mRNA expression levels that were positively correlated with mutant TIMM17A expression whereas 12,344 genes (dark green dots) demonstrated mRNA expression levels that were negatively correlated with mutant TIMM17A expression [false discovery rate (FDR) < 0.01]. The 50 significant gene sets exhibited both positive and negative correlations with TIMM17A, as shown by a heat map (Figures 6B,C). , and RABIF (r = 0.623, p = 9.86E−120). GO analysis indicated that these networks are significantly enriched in genes with functions related to "mitochondrial inner membrane, " "mitochondrial protein complex, " and "chromosomal region." Furthermore, these genes contribute to "non-coding (nc)RNA processing, " "chromosome segregation, " and "rRNA metabolic process, " and act on "RNA and structural constituents of ribosome" (Figures 7A-C and Supplementary Tables 1-3). KEGG pathway analysis indicated that these genes are related to "oxidative phosphorylation, " "spliceosome, " and "ribosome" pathways ( Figures 7D,E and Supplementary Table 4). Thus, genes coregulated with mutant TIMM17A in BRCA participate extensively in "oxidative phosphorylation" and "cell cycle regulation." To further validate this network-level relationship between TIMM17A and cell cycle regulator, we examined the association of TIMM17A expression with cell cycle phase distribution in BRCA cell lines by flow cytometry. The proportions of cells in S-phases were reduced by TIMM17A knockdown (Figures 8A-E), suggesting that TIMM17A serves to enhance cancer cell proliferation.
The associated transcription factors were mainly of the E2F transcription factor (E2F) family, including E2F_Q6, E2F1_Q6, E2F1DP2_01, E2F_02, and GGAANCGGAANY ( Table 2 and LeadingEdgeNum, the number of leading edge genes; FDR, false discovery rate from Benjamini and Hochberg gene set enrichment analysis (GSEA) vs. the annotation found in Molecular Signatures Database (MSigDB) for transcription factors (TF).
FIGURE 9 | Protein-protein interaction network of ATR kinase targets (GeneMANIA). Protein-protein interaction (PPI) network and functional analysis of the gene set enriched in targets of CDK1. Different colors of the network edge indicate the bioinformatics methods applied: co-expression, co-localization, physical interactions, pathway, predicted, genetic interactions, and shared protein domains. The different colors for the network nodes indicate the biological functions of the set of enrichment genes. 5-7). The protein-protein interaction network constructed by GeneMANIA revealed many correlations with CDK1, miRNA-331, and TF E2F_Q6.

Supplementary Tables
The gene set enriched in CDK1 targets is responsible mainly for regulating cell cycle G2/M-phase transition, cell division, organelle fission, and the cell cycle checkpoint (Figure 9). To validate the association between CDK1 and TIMM17A, we analyzed the correlation between TIMM17A expression and CDK1 expression both in silico and in BRCA cell lines. Expression of CDK1 was positively correlated with TIMM17A expression in theUALCAN database (Pearson's r = 0.32) ( Figure 10A). Furthermore, western blotting revealed that TIMM17A knockdown reduced CDK1expression in BRCA cell lines ( Figure 10B). The gene set enriched in miRNA-331 targets is involved mainly in regulation of cell-cycle adhesion ( Table 2 and Supplementary Figure 2), so we also examined the correlation between TIMM17A and miRNA expression in BRCA cell lines. Indeed, TIMM17A knockdown altered the expression level of miRNA331 as well as expression of miRNA219 and miRNA 326 in BRCA cell lines (Figures 10C,D). Similarly, the gene set enriched in transcription factor E2F_Q6 targets is involved mainly in regulation of DNA replication, cell cycle checkpoint, and mitotic cell cycle (Supplementary Figure 3).

DISCUSSION
Most BRCA cases are caused by a germline mutation in a single gene or multiple genes, of which the most frequent are in BRCA1, BRCA2, and p53 (Tung et al., 2016). Our analysis revealed an inverse association between TIMM17A expression level and both OS and DFMS. About 90% of all tumor recurrences and metastases occur within 10 years of first diagnosis, so we chose 10 years as the limit for survival analyses. Surprisingly, survival analyses of multiple public datasets for this period revealed that patients with high expression of TIMM17A have worse clinical outcome. Further subgroup analysis of multiple clinicopathological features revealed that TIMM17A mRNA expression was higher in BRCA patients than healthy individuals but did not differ according to sex, age, ethnicity, tumor stage, tumor grade, or BRCA subtype. The TIMM17A copy number was also significantly higher in BRCA tissues than normal tissues, and amplification was the most common genomic alteration with a change in copy number >2.
Expression of mutant TIMM17A was strongly correlated with expression levels of UBE2T, TMEM183A, CACYBP, SNRPE, and RABIF (top 5). The functional network of TIMM17A in BRCA participates extensively in chromosome segregation, oxidative phosphorylation, and cell cycle regulation. Indeed, flow cytometry demonstrated that TIMM17A expression was strongly associated with cell cycle phase distribution in BRCA cell lines. Specifically, TIMM17A knockdown reduced the proportions of cells in S-phases, suggesting that TIMM17A enhances the rate of proliferation. Collectively, these results demonstrate that high TIMM17A expression enhances BRCA aggression, potentially by promoting fasting cell cycle progression.
LinkedOmics analysis identified numerous additional network components and other gene sets with high weight coexpression and protein interactions with TIMM17A. Genes coregulated with TIMM17A were enriched in cell cycle modulators and targets of CDK1, MIR-331, and E2F family transcription factors. Activation of CDK1 induces phosphorylation of lamin, disintegration of the nuclear fiber layer and nuclear membrane, phosphorylation of histone H1, and chromosome condensation (Haneke et al., 2020), processes required for cell division. The evolution of kinase signaling is promoted by gain or loss of phosphorylation sites in rapidly evolving regions (Holt et al., 2009). Ravindran Menon et al. (2018 found that tumor initiation is associated with regulation of transcription factor Sox2 activity by CDK1, suggesting the CDK1-Sox2 pathway as a potential therapeutic target. Mitotic chromosome segregation also depends on H2B serine 6 phosphorylation by CDK1 (Seibert et al., 2019). Unregulated CDK1-mediated phosphorylation allows the cell cycle to proceed unchecked, leading to tumor formation (Cheng et al., 2020). The current analysis identifies multiple additional CDK1-targeted pathways for potential therapeutic intervention. Moreover, we confirmed a strong correlation between TIMM17A and CDK1 in BRCA cells at both transcriptional and protein expression levels as evidenced by the marked CDK1 downregulation following TIMM17A knockdown. Xuefang et al. (2020) reported that miR-331-3p promoted apoptosis of nasopharyngeal carcinoma cells by targeting the elf4B-PI3K-AKT pathway, and Chen et al. (2018) recently reported that miR-331 inhibits the proliferation and invasion of melanoma cells. At the beginning of S-phase, CDK2 combines with cyclin A to inactivate E2F transcription factors which is a precondition for the completion of S-phase, while persistent E2F activity leads to apoptosis (Krek et al., 1995;Fueyo et al., 1998;Lukas et al., 1999). Therefore, selective inhibition of CDK2/cyclin A may increase E2F expression, which in turn could lead to S-phase arrest or apoptosis. Transcription factor E2F1 (Schuldt, 2011) plays an important role in cell cycle control and tumor suppressor gene function, and it is also the target of transforming protein of small DNA oncoviruses (Krek et al., 1995;Chen et al., 2009). These E2F family members contain multiple highly conserved domains (Osorio, 2015;Kent and Leone, 2019), including a DNA-binding domain, dimerization domain interacting with transcription factor protein (DP) regulated by differentiation, an acid amino acid-rich transactivation domain, and a tumor suppressor protein-related domain within the transactivation domain. In addition, E2F proteins E2F2 and E2F3 have cyclin-binding domains. These proteins preferentially bind to retinoblastoma protein pRb in a cell cycle-dependent manner, allowing both cell proliferation and p53-dependent or p53-independent apoptosis (Osorio, 2015). These results provide potential directions for future research on the molecular mechanisms of BRCA.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.

AUTHOR CONTRIBUTIONS
JJC and CXW contributed to data analysis. JJC, JYC, and LH drafted the manuscript. LH and JJC performed the experiments. ZHS and QZ designed the experiments and revised the manuscript. All of the authors critically revised the manuscript, gave final approval of the version to be published, and agreed to be accountable for all aspects of the work.

FUNDING
This work was financially supported by a scholarship from the National Natural Science Foundation of China (31901035) and the Guangdong Province Foundation for Basic and Applied Basic Project (2020A1515010951).

ACKNOWLEDGMENTS
We would like to acknowledge the help from the staff at the General Hospital of Southern Theatre Command of PLA.