Interactome Analysis of Microtubule-targeting Agents Reveals Cytotoxicity Bases in Normal Cells

Cancer causes millions of deaths annually and microtubule-targeting agents (MTAs) are the most commonly-used anti-cancer drugs. However, the high toxicity of MTAs on normal cells raises great concern. Due to the non-selectivity of MTA targets, we analyzed the interaction network in a non-cancerous human cell. Subnetworks of fourteen MTAs were reconstructed and the merged network was compared against a randomized network to evaluate the functional richness. We found that 71.4% of the MTA interactome nodes are shared, which affects cellular processes such as apoptosis, cell differentiation, cell cycle control, stress response, and regulation of energy metabolism. Additionally, possible secondary targets were identified as client proteins of interphase microtubules. MTAs affect apoptosis signaling pathways by interacting with client proteins of interphase microtubules, suggesting that their primary targets are non-tumor cells. The paclitaxel and doxorubicin networks share essential topological axes, suggesting synergistic effects. This may explain the exacerbated toxicity observed when paclitaxel and doxorubicin are used in combination for cancer treatment.


Introduction
In 2013, 14.9 million cases of cancer, 8.2 million deaths, and 196.3 million disability adjusted life years (DALYs) were reported [1]. The standardized incidence levels by age for all cancer types combined increased by more than 10% in 113 countries from 1990 to 2013 [1], with developing countries among the most affected. For this reason, the development of effective cancer drugs is a global priority.
Microtubules are among the most extensively studied therapeutic targets in cancer treatment [2,3]. Made up of a and b tubulin protofilaments, microtubules can serve as the basis of a cellular network, which are essential for determining cell shape and organization through the movement of organelles, such as the nucleus, mitochondria, endoplasmic reticulum, and Golgi apparatus [4][5][6]. In addition, this network plays a regulatory role in cell migration and adhesion as well [6,7].
Microtubule-targeting agents (MTAs) can affect microtubule stability, leading to disruption of the mitotic spindle and cell death [8], and are therefore one of the most effective classes of drugs used in chemotherapy against cancer. According to their binding property with tubulin, MTAs are classified as stabilizing agents (known as vinca alkaloids) [9,10], which bind the tubulin polymer, and destabilizers (known as taxanes) [11][12][13], which bind tubulin dimers. The stabilizing vinca-type drugs, including vincristine, vinblastine, and vinorelbine, are used for haematological cancers such as lymphomas and leukaemia, whereas destabilizing taxane-type drugs, including paclitaxel and docetaxel, are used to treat solid cancers such as breast, ovary, and oesophageal cancer [14].
However, high in vivo cytotoxicity and increased tumor resistance have been observed with the usage of MTAs [15,16], since they can recognize microtubules of interphase cells [17]. MTAs kill rapidly-dividing cells by arresting them in mitosis [18], but how MTAs kill slowly-dividing tumor cells has not been fully described. It has been proposed that MTAs interfere with microtubule trafficking system in prostate cancer patients [19], enabling the use of MTAs in a combination regimen with DNA-damaging agents (DDAs), such as doxorubicin, for cancer treatment [20].
One method for studying the interactions between drugs and the human proteome is to use network theory, in which proteins and/or chemical compounds are represented by nodes, and the interactions between them are represented by edges. As a result, inference of cellular processes affected by certain drugs can be inferred by studying the connections between nodes. Protein network reconstruction has allowed a better understanding of the physiopathological mechanisms of cancer, identified genes associated with specific pathologies to improve disease classification, and defined specific cancer sub-networks important for the identification of therapeutic targets [21][22][23][24][25].
In this study, we created an integrated interaction network between MTAs and the proteome of a non-cancerous cell to identify the essential topological axes of its functional profile. Potential cellular mechanisms associated with the cytotoxicity of MTAs are also described.

Different biological function detected from MTA subnetworks
The related data for MTA subnetworks (Table S1) showed that the microtubule inhibitor 2-methoxyestradiol triggers the metabolism of flavonoids, lipids, vitamins, and podophyllotoxins. In addition, association of 2-methoxyestradiol with cellular adhesion and adiponectin secretion was also detected. Colchicine and vincristine were involved in the cellular stress, death, and proliferation as well. Besides, combretastatin A4, epothilone, epothilone B, and vinblastine are shown to induce the formation of protein complexes and microtubule-based organelle movements, and might induce apoptosis of lymphocytes. Spongistatin was involved in the regulation of NF-jb, whereas vindesine, vinorelbine, and tasidotin induce cytotoxicity. Finally, nocodazole, noscapine, and paclitaxel were involved in protein phosphorylation, nuclear membrane reorganization, and apoptosis.
Integrated MTA interactome shows a common proteome core In order to establish the common interactome for the MTAs, we merged all the subnetworks for MTAs together. It was determined that 71.4% of the integrated network shared a common core formed by 363 nodes and 2327 connections ( Figure 1).
To create a comparison pattern for the integrated network, a randomized network was generated with similar topological characteristics. We observed that the integrated network had a positive correlation between the number of nodes and the node degree (R = 0.935; R 2 = 0.799) compared with the random network (R = À0.045; R 2 = 0.04; Figure 1A and B). Additionally, the distribution of betweenness centrality (BC) and closeness centrality (CC) indicated that the integrated network was more compact than the random network ( Figure 1CÀF). The mean value for the clustering coefficient for the integrated network (0.583) was significantly higher than that for the random network (0.035) (Wilcoxon Z-value = À12.2321, P = 1EÀ5, significant at P < 0.01). Likewise, the mean k value for the integrated network (26.122) was significantly higher than that for the random network (12.822) (Wilcoxon Z-value = À12. 322, P = 1EÀ5, significant at P < 0.01). These results suggest that all nodes in the integrated network have a distribution that could be attributed to their functional properties, indicating the presence of functionally-essential nodes within the network (topological axes).
To identify the topological axes in the integrated network, the distribution of the networks according to the node degree was analyzed. This distribution showed that the different tubulin conformations were the most connected nodes. Moreover, there are also other topological axes, including TP53, AKT1, CDK1, PCCNA, JUN, VEGFa, CASP3, CCNB1, and RAC1. The encoded proteins regulate cellular processes such as apoptosis, cell differentiation, cell cycle control, replication control in eukaryotes, stress response, and energy metabolism regulation. In addition, we also observed high degree values for paclitaxel and docetaxel, indicating that these MTAs can make many connections with the human proteome ( Table 1).
To determine the role of these topological axes in information flow, we analyzed the network topology in terms of stress and BC. We found a positive correlation between these two parameters in the integrated network (R = 0.897; R 2 = 0.814) ( Figure 2). The nodes with the highest stress and BC coincided with those identified using node degree. These results suggest that paclitaxel, TUBB2A, TP53, and AKT1 are the topological axes that control and regulate the information flow in the integrated network.
We then performed the cluster analysis of the network to identify the functional profile of its main subnetworks. Consequently, two functional domains were revealed. The first functional domain was related to apoptosis regulation, proliferation, locomotion, and cell migration, whereas the second domain was related to drug catabolism ( Table 2). To identify the role of these functional domains in cancer metabolism, we mapped the nodes into the KEGG database. As shown in Table 3, nodes are mapped to pathways related to apoptosis, cell proliferation, and cancer metabolism. Figure 1 Comparative analysis of the topological properties between the random network and the MTA interactome A. Node distribution according to the node degree relative to the number of neighbours for random network. B. Node distribution according to the node degree relative to the number of neighbours for MTA interactome. C. Node distribution for betweenness centrality relative to neighbour distribution for random network. D. Node distribution for betweenness centrality relative to neighbour distribution for MTA interactome. E. Node distribution for closeness centrality relative to the number of neighbours for random network. F. Node distribution for closeness centrality relative to the number of neighbours for MTA interactome. MTA, microtubule-targeting agent. These results suggest the presence of functionally-essential nodes within the network (topological axes) ( Figure 3). The topological axes in the DDA network were obtained based on node degree. We found that TP53, UBC, AKT1, MYC, EGFR, SRC, EP300, JUN, CREBBP, and CCND1 proteins were key components. These proteins regulate cellular processes such as apoptosis, cell differentiation, cell cycle control, replication, stress response, and transcriptional activation ( Table 3).
The paclitaxel subnetwork was the principal hub of the MTA interactome with 166 nodes and 1497 edges and represents 33.8% of the MTAs network. The paclitaxel and doxoru- bicin subnetworks were merged to produce the MTA-DDA integrated network. The MTA-DDA network shows a positive correlation for all the topological parameters evaluated, including node degree (R = 0.734; R 2 = 0.805), BC (R = 0.953; R 2 = 0.821), and CC (R = 0.953; R 2 = 0.821). In addition, the merged MTA-DDA network shared the same essential topological axes for the DDA network (Table S1). The functional profile of the MTA-DDA network revealed only one functional domain associated with the regulation of biological process, apoptosis, cellular cycle, and response to cellular stress.

Discussion
The MTA integrated network was developed to describe the basal drug-protein interactions in a non-cancerous cell. The induction and regulation of apoptosis by MTAs should be considered as the central hub for all the subnetworks. However, a direct interaction between tubulin and any MTA was not observed, which opens two possibilities: (1) the apoptotic induction is mediated by secondary targets and (2) there is a common MTAs interactome core in the cell.
Our network analysis indicates that MTAs have common functional mechanisms, consistent with experimental observations showing that low concentrations of these compounds inhibit mitosis in a synergistic fashion [26]. The network analysis also allows the identification of interactions related to mitosis and microtubule activity during interphase.
It has been proposed that the functional advantages of MTAs are based on their inhibition of mitosis during interphase [27]. Our MTA interactome analysis revealed proteins and signaling pathways that are key points of regulation by these compounds, which include apoptosis, cell proliferation, cell cycle, and drug catabolism. However, these advantages of MTAs are primarily observed in vitro at the cell culture level [28,29], as these compounds are associated with cytotoxicity in vivo [30,31] and have harmful effects on rapidly dividing cells [32].
The cells of solid tumors show a high proliferation rate but a low apoptotic index [27], also known as the proliferation rhythm paradox [33]. The secondary effects of these compounds are observed mostly in fast-growing cells, such as haematopoietic and epithelial cells, leading to the conclusion Figure 2 Essential nodes for information flow in the integrated MTA network A fitted line and its slope were calculated to identify the most significant nodes and MTAs that modulate the information flow in the integrated network, using the betweenness centrality and stress values. Betweenness centrality represents the number of times a node is visited and stress indicates how many times a particular node is part of different shortest paths. The most relevant MTAs in the information flow of the network are underlined. that the main target of MTAs is the microtubules of interphase cells [27]. One of the main functions of interphase microtubules is to serve as ''railways" for protein, vesicle and organelle trafficking [34]. The MTA interactome in this study revealed a large number of potential secondary targets (nodes), which were represented as client proteins of the interphase microtubules. Two of these microtubule accessory proteins, dyneins DYNC1H1 and SYNC2H1 that transport cellular elements, were identified as topological axes for the MTA interactome. Taxanes antagonize the androgen receptor signaling pathway in prostate cancer cells by blocking dynein-mediated protein trafficking along interphase microtubules [35,36], and paclitaxel decreases the endocytic trafficking of epidermal growth factor receptor in lung cancer cells [37]. In normal cells, these effects on dynein function lead to the accumulation of proteins in the cytosol and ultimately to apoptosis. Poruchynsky et al showed that in cells treated with doxorubicin, the proteins ATM, ATR, DNA-PK, RAD50, MRE11, p95/NBS1, p53, 53BP1, and p63, which are involved in DNA repair, were sequestered when vincristine and paclitaxel were added [20].
One of the most significant biological functions of the MTA interactome as represented by microtubule accessory proteins was the regulation of apoptosis. In the MTA interactome, both the catalytic a and the constitutive b subunits of chaperone HSP90 were revealed to be essential topological axes. HSP90 is vital in the folding process of several proteins involved in cell proliferation and apoptosis. Both subunits bind to the acetylated tubulin in a similar fashion as Akt and p53 [38], two other topological axes of the network. In fact, it is known that the cytosolic accumulation of p53 is induced by paclitaxel, vincristine, and nocodazole in lung cancer cells [39].
Furthermore, Pcl-2, which is known to interact with paclitaxel and initiates the pro-apoptosis cascade [40], was also identified as a topological axis for the MTA interactome. In human leukaemia cell lines and in the cells of patients with chronic lymphocyte leukaemia treated with paclitaxel, there is a significant increase in the expression of cytochrome c, caspase 3 and PARP, as well as a decreased JNK activity [41]. It is of note that the proteins JNK, caspase 3, and PARP were also identified as topological axes.
Another essential node that is tightly associated with apoptosis is the protein Bim EL (MCL2L11), which is known as the most important physiological antagonist of survival proteins that are predominant in T lymphocytes. When phosphorylated, Bim EL can leave the microtubule and be cleaved by caspase 3, thus activating mitochondrial-and/or receptor-based apoptotic mechanisms [42]. Our study is consistent with the observation that a cell can enter apoptosis through the interaction of caspase 3 and Bim EL [42].
A second functional domain identified in the MTA interactome is related to xenobiotic metabolism. This type of metabolism consists of two phases, the first involving CYP450 proteins that detoxify or bioactivate chemical compounds through chemical functionalization and the second involving the exposure of molecules to conjugation reactions to make them more hydrophilic and susceptible to degradation [43]. The metabolism of paclitaxel is mediated by CYP1A2, 1B1, 2A6, 2C9, 2E1, and 3A4 [44], with 1A2 and 2C9 identified as the essential topological axes via the interactome analysis.
Using the non-cancerous cell MTA interactome, the protein components responsible for the induction and regulation of apoptosis both at the mitochondria and receptor levels could be revealed. Another functional component we find is the metabolic mechanisms for drug elimination, which was identified through several proteins involved in xenobiotic metabolism. These results suggest that the MTAs disrupt essential regulators of normal cell physiology present in both healthy and cancerous cells, mainly dependent on their effects on interphase microtubules.

Protein-compound interaction network design
A list of agents targeting microtubules was obtained via a literature search, and the classification system developed by Bazan et al [45] was used as a reference. The names of the compounds were searched in the PubChem database to obtain the corresponding International Chemical Identifier (InChI) code. Interaction networks were created for each of the following compounds: (1) MTAs, including 2-methoxyestradiol, colchicine, combretastatin A4, docetaxel, epothilone, epothilone B, estramustine, nocodazole, noscapine, paclitaxel, podophyllotoxin, spongistatin, tasidotin, vinblastine, vincristine, vindesine, and vinorelbine; and (2) DDAs, including doxorubicin and darboplatin. The initial subnetworks for MTAs and DDA were obtained using STICH 4.0 database (http://stitch.embl.de/), with the following criteria: a confidence score of 0.500 with 500 interactions, forcing search saturation, and all prediction methods being active.

Construction of merged MTA networks and topological analysis
The MTA subnetworks were imported into Cytoscape version 3.4.0 [46] and merged using the Merge Networks plugin [http:// www.cytoscape.org/plugins2.php] included in Cytoscape by default. The following integrated network topology parameters were obtained using the NetworkAnalyzer plugin: level (k, degree), indicating how one node is connected with the others; intermediation (BC), indicating the number of shortest paths that pass through a node; CC, indicating which nodes are closer to the centre of the network; stress, indicating how many times a particular node is part of different shortest paths; and the clustering coefficient. As a control, random networks with 363 nodes and 2327 connections were generated with the mean k values of the integrated network using the Randomizer Network plugin [47]. The integrated network was analyzed using the MCODE plugin [48] to detect the main clusters, and GO analysis was performed for both the integrated network and the clusters using the Biological Network Gene Ontology (BiNGO) plugin. Hypergeometric distribution and false discovery rate (FDR) correction were adopted to determine the functional richness level for each network in each category. The FDR correction with a significance level of 0.05 is shown for two descriptive categories of the clusters, the first representing biological function and the second representing the node mapping in KEGG [49].
A network between paclitaxel and doxorubicin was merged and analyzed to test the hypothesis about a synergistic effect between MTAs and DDAs similarly as described above.