Identification of genes associated with methotrexate resistance in methotrexate-resistant osteosarcoma cell lines

Background This study aimed to better understand the mechanisms underlying methotrexate (MTX)—resistance in osteosarcoma. Methods The raw transcription microarray data GSE16089 collected from three MTX-sensitive osteosarcoma (Saos-2) cell samples and three MTX-resistant osteosarcoma (Saos-2) cell samples were downloaded from Gene Expression Omnibus. After data processing, the differentially expressed genes (DEGs) were identified. Next, DEGs were submitted to DAVID for functional annotation based on the GO (Gene Ontology) database, as well as pathway enrichment analysis based on the KEGG (Kyoto Encyclopedia of Genes and Genomes) database. Transcription factors (TFs) and tumor-associated genes (TAGs) were identified with reference to TRANSFAC and TAG, and TSGene databases, respectively. The protein-protein interaction (PPI) network of the gene-encoded products was constructed, and the subnetwork with the highest score was also detected using Search Tool for the Retrieval of Interacting Genes and BioNet package. Results A total of 690 up-regulated genes and down-regulated 626 genes were identified. Up-regulated DEGs (including AARS and PARS2) were associated to transfer RNA (tRNA) aminoacylation while down-regulated DEGs (including AURKA, CCNB1, CCNE2, CDK1, and CENPA) were correlated with mitotic cell cycle. Totally, 13 TFs (including HMGB2), 13 oncogenes (including CCNA2 and AURKA), and 19 tumor suppressor genes (TSGs) (including CDKN2C) were identified from the down-regulated DEGs. Ten DEGs, including nine down-regulated genes (such as AURKA, CDK1, CCNE2, and CENPA) and one up-regulated gene (GADD45A), were involved in the highest score subnetwork. Conclusion AARS, AURKA, AURKB, CENPA, CCNB1, CCNE2, and CDK may contribute to MTX resistance via aminoacyl-tRNA biosynthesis pathway, cell cycle pathway, or p53 signaling pathway.


Background
Methotrexate (MTX) was first introduced to replace aminopterin to treat acute lymphocytic leukemia, which works via inhibiting dihydrofolate reductase (DHFR), a key enzyme required in intracellular folate metabolism, leading to decreased tretrahydrofolate coenzyme level, accordingly achieving the inhibition of thymidylate and the biosynthesis of DNA and purine. So far, different mechanisms have been presented to address the intrinsic and acquired MTX resistance: (1) decreased MTX transport, (2) impaired MTX polyglutamylation, (3) increased DHFR enzyme activity, (4) altered affinity of MTX for DHFR, and (5) increased MTX efflux due to elevated levels of the multidrug resistance protein (MRP) [1].
Osteosarcoma is a primary malignant, highly vascularized bone tumor, mainly occurring in adolescents and children [2,3]. The unclear understanding of the underlying molecular mechanism greatly hinders the therapy of osteosarcoma [4]. Currently, multiagent chemotherapy, usually using doxorubicin, cisplatin, and high-dose MTX, has improved the survival of osteosarcoma patients from 11 to 70 % [5]. However, MTX resistance has become an issue of growing interest, as little information is available in this disease up to now [6]. TP53, a tumor suppressor gene (TSG), encodes a transcriptional regulator that responds to DNA damage or cellular stress and controls the progression and apoptosis of cell cycle. As previously reported, the accumulation of p53 protein is probably a predictor of response to methotrexate (MTX) [7]. p53 alterations increase the risk for the development of drug resistance by altering MTX transport [8].
Using the transcription profiles of three MTXsensitive osteosarcoma (Saos-2) cell lines and three MTX-resistant Saos-2 cell lines and analyzing the network of the differentially expressed genes (DEGs), Selga et al. have found the alteration in the expression of a number of genes, such as eukaryotic translation elongation factor 1 alpha 1 (EEF1A1) in the MTX-resistant osteosarcomas (Saos-2) cell lines, pancreatic cancer, and erythroblastic leukemia cell lines [9]. However, further systematic analyses, including GO (Gene Ontology) functional and REACTOME pathway enrichment analysis, were not performed for DEGs concerning osteosarcoma cells in their study. REACTOME is a knowledgebase of human reactions and pathways, which provides an integrated view of the molecular details of human biological processes ranging from metabolism to DNA replication and repair to signaling cascades [10], and has been used in various studies [11]. The transcription microarray data GSE16089 deposited in Gene Expression Omnibus (GEO), which includes three chips from MTX-sensitive and three from MTX-resistant osteosarcoma (Saos-2) cell lines [9], were downloaded and analyzed in this study so as to better understand the genetic etiology of osteosarcoma. The DEGs were identified, and the functional and pathway enrichment analysis was performed for them. The protein-protein interaction (PPI) network of the gene productions and its subnetwork were analyzed. These findings in this study will encourage us to investigate the anti-cancer effects of the DEGs or the pathways as well as the MTX resistance in osteosarcoma.

Material and methods
As the study did not involve any human or animal study, the ethical approval was not required.

Microarray data
Gene expression microarray dataset deposited in the National Center of Biotechnology Information (NCBI) GEO (http://www.ncbi.nlm.nih.gov/geo/) with the accession number of GSE16089 was downloaded [9]. The annotation platform was GPL570 [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array. According to contributors, three samples of either MTX-sensitive cells or MTX-resistant cells of the Saos-2 osteosarcoma cell line were used for gene expression analyses [9]. Saos-2 cell line was sensitive to MTX, and its MTX-resistant cells were obtained in the laboratory via incubation with stepwise concentrations of MTX (Lederle) as described previously [12].

Data processing
The raw probe profile data was downloaded from GEO. The processing of the raw microarray data was performed by robust multiarray average (RMA) using R/ Bioconductor package Affy [13]. The preprocessing consisted of background correction, quantile normalization, and probe summarization of expression value. The gene expression matrixes were obtained for further analysis.

Identification of DEGs
Transcriptional sets were mapped to NCBI entrez genes using Gene ID converter [14]. The averaged value was calculated for further analysis if there were multiple probe sets corresponding to the same gene. Probes were filtered if they corresponded to multiple genes. The classical t test was performed among the samples to identify the genes specifically differentially expressed between MTX-sensitive and MTX-resistant Saos-2 cell lines. The cut-off criteria for the DEGs were set at p value <0.05 and |log2FC (fold change)| > 1.

Functional and pathway enrichment analysis
GO and KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway analysis provides prediction of gene function and informs people of how molecules or genes work [15,16]. The DEGs were submitted to DAVID (Database for Annotation, Visualization, and Integrated Discovery) (http://david.abcc.ncifcrf.gov/) to find the significantly enriched biological process (BP) terms, molecular function (MF) terms, and cellular component (CC) terms based on the GO (Gene Ontology) database, as well as pathways based on the KEGG (Kyoto Encyclopedia of Genes and Genomes) database. For the identification of the significantly enriched biological processes in detail, the significantly altered DEGs were subjected to the REAC-TOME knowledgebase (http://www.reactome.org). The thresholds for the significant associated GO functional category and pathways were set at p value < 0.01.

Identification of transcription factors and tumorassociated genes
TRANSFAC (http://www.gene-regulation.com/index2) is a database on transcription factors, their genomic binding sites, and DNA-binding profiles [17]. To identify the DEGs that also act as transcription factors, transcription factor (TF) prediction was performed using the TRANSFAC database. TAG (tumor-associated gene) database (http://www.binfo.ncku.edu.tw/TAG/) is a semiautomatic information retrieving engine which collects specific information about genes from various resources. TSGene database (http://bioinfo.mc.vanderbilt.edu/TSGene/ ) is a resource of tumor suppressor genes (TSGs) that provides a comprehensive TSG catalog for advanced systems biology-based analysis for the cancer research community [18]. TAGs including oncogenes and TSGs were also identified from the DEGs using TAG and TSGene databases, respectively.
Construction of PPI network and the subnetwork analysis STRING (Search Tool for the Retrieval of Interacting Genes) is a web server to retrieve and display the repeatedly occurring neighborhood of a gene which generalizes access to protein interaction data, by integrating known and predicted interactions from a variety of sources [19]. To describe the interactive network of DEGs, the STRING database was used to build the PPI network of encoding products of all of the DEGs. A STRING score of 0.4 was set as the reliability threshold. Cytoscape, a standard tool for integrated analysis and visualization of biological networks, was used to visualize the PPI network [20]. The connectivity degree analysis was performed, and hub nodes were obtained using the scale-free properties of PPI networks. The BioNet package is an R-Package for the functional analysis of biological networks and is used for the mining of the sub-networks in the PPI network [21]. The highest scoring subnetwork was obtained. The threshold of the given false discovery rate (FDR) value was 0.0001.

Functional and pathway enrichment analysis
According to the GO annotation, the up-regulated DEGs were functionally involved in BP terms such as response to endoplasmic reticulum stress (including AARS and CCND1) and transfer RNA (tRNA) aminoacylation for protein translation (including AARS, TARS, YARS, and PARS2), and MF terms such as aminoacyl-tRNA ligase activity (including AARS and PARS2), as well as CC terms such as intracellular membrane-bounded organelle (including CCND1 and GADD45A) ( Table 1). Moreover, the up-regulated DEGs were mainly associated with the KEGG pathways such as aminoacyl-tRNA biosynthesis (including AARS and PARS2). The mainly related REAC-TOME pathways were cytosolic tRNA aminoacylation   (including AARS) and tRNA aminoacylation (including AARS and PARS2) ( Table 2). The down-regulated DEGs were significantly enriched in BP terms such as microtubule cytoskeleton organization and mitotic cell cycle (including AURKA, AURKB, CCNB1, CDK1, and CENPA) and in MF terms such as protein binding and nucleotide binding (including AURKA, AURKB, CCNB1, CDK1, and CENPA), as well as in CC terms related to chromosome, centromeric region (including AURKB, CCNB1, CENPA, and CENPH), and kinetochore (including CCNB1, CENPA, and CENPH) ( Table 3). Also, the down-regulated DEGs were significantly involved in KEGG pathways of cell cycle, oocyte meiosis, and p53 signaling pathway (including CCNB1, CCNE2, and CDK1). The relevant REACTOME pathways were resolution of sister chromatid cohesion and mitotic M-M/G1 phases (including AURKB, CCNB1, CDK1, CENPA, and CENPH) ( Table 4).

Discussion
In this study, a total of 690 up-regulated DEGs and 626 down-regulated DEGs were identified in MTX-resistant osteosarcoma cells. According to the functional and pathway enrichment analysis, the up-regulated DEGs such as AARS, TARS, YARS, and PARS2 were mainly associated with tRNA aminoacylation, and the downregulated DEGs such as AURKA, AURKB, CCNB1, CDK1, and CENPA were mainly correlated with mitotic cell cycle. The up-regulated DEGs, including AARS, TARS, YARS, and PARS2, were biologically related to tRNA aminoacylation. AARS, TARS, and YARS encode aminoacyl-tRNA synthetases alanyl-tRNA synthetase, threonyl-tRNA synthetase, and tyrosyl-tRNA synthetase, respectively, which catalyze the aminoacylation of tRNA by their cognate amino acid, and thus are necessary for protein synthesis. Aminoacyl-tRNA synthetases usually take distinct roles in inflammation and transcriptional regulation [22]. The aberrant expression and cellular localization of aminoacyl-tRNA synthetases disturb normal cell regulatory networks and cause abnormalities through multiple routes [23]. For example, inhibition of osteosarcoma cell migration might be related to the extracellular functions of TARS [24]. MTX has been used to treat antisynthetase syndrome, a type of heterogeneous autoimmune disorder, in which autoantibodies target anti-aminoacyl-transfer RNA synthetase for specific amino acid. The up-regulation of aminoacyl-tRNA synthetases in MTX-resistant cells may indicate that this type of cells has developed an ability to improve the expression level of aminoacyl-tRNA synthetases in the presence of MTX that offset the effect of MTX, leading to MTX resistance. The down-regulated genes were dominantly related to mitotic cell cycle or nuclear division, for example, AURKA, AURKB, CENPA, CDK1, CCNB1, and CCNB2. AURKA is known as an oncogene. AURKA is an important regulator to G2/M transition [25]. AURKA protein is mainly located at the microtubule organizing center at the metaphase I (M1) of oocytes [26]. Further investigations have shown that the localization of AURKA in the area of aligned chromosomes is consistent with the AURKA-dependent phosphorylation of kinetochore component centromere protein A (CENPA). CENPA phosphorylation requires the enrichment of AURKB to maintain the phosphorylation on Ser7 at inner centromeres and for kinetochore function [26][27][28]. Jiang et al. [29] showed that silencing of AURKA expression in osteosarcoma cells significantly decreased both colony formation ability in vitro and tumorigenesis ability in vivo as well as induced cell apoptosis and G2/M cell cycle arrest in osteosarcoma cells. In addition, the median survival time was significantly longer in patients with low-CENPA expression osteosarcomas than in those with high-CENPA expression osteosarcomas. AURKA and CENPA had, respectively, been identified as a susceptibility gene and an independent poor prognostic factor for osteosarcoma [25,30]. The down-regulation of mitosis-related genes seems contrary to the commonly known fact that mitosisrelated genes are usually up-regulated in tumors. However, up-regulation of mitosis-related genes is just recognized in tumors, not in drug-resistant tumor cells. Actually, we know little about the anti-drug mechanisms in tumors. It is possible that down-regulation of mitosis-related genes is a protection mechanism in MTX-resistant cells, or just at certain time point. Apparently, this finding has to be validated experimentally, since the sample number is so small. Furthermore, CDK1, CCNB1, and CCNB2 which were enriched in mitotic cell cycle, oocyte meiosis, and p53 signaling pathways were inactivated by MTX treatment in osteosarcoma cell lines. As we know, deregulated cell proliferation and tumor-associated cell cycle always propel the complexity and idiopathy of cancer [31,32]. Tumor-associated cell cycle is often mediated by the alterations of cyclin-dependent kinase (CDK) activities Fig. 2 The highest score subnetwork for the differentially expressed genes (DEGs). Red stands for the productions of up-regulated DEGs and green for productions of down-regulated DEGs. Round dot indicates the highly associated DEGs with osteosarcoma and square indicates less relevant DEGs [33]. Cyclin B1 (CCNB1) to which p53 is directly bonded mediates G2/M progression and inhibits cell division. As reported, the inhibition of cyclin B and CDK1 led to the arrest of osteosarcoma cell division [34]. Inactivation of CDK1 and CDK2 triggered the apoptosis of osteosarcoma cells [35]. Moreover, MTX prevents tumor cells from proliferating by inhibiting dihydrofolate reductase (DHFR) [36]. The inhibition of CDK reduces the expression of both DHFR mRNA and protein thus enhancing sensitivity of human osteosarcoma cell lines to MTX [36,37]. These studies concluded that the usage of combination of cyclin-CDK inhibitors and MTX which regulated mitotic cell cycle and p53 signaling pathway might overcome MTX resistance in osteosarcoma cells.
However, the DEGs identified in the present study were not the same with those identified to be associated with MTX resistance by Selga et al. using seven cell lines of different types of cancer [9]. This discrepancy may be attributed to the fact that genes identified by bioinformatics methods often vary with the criteria you adopt for analysis, and that Selga et al. focused on seeking genes commonly expressed in different MTX-resistant tumors, whereas we only paid attention to those specifically related to MTX resistance developed in osteosarcoma cell line. However, since there are only three samples for either MTX-sensitive or MTX-resistant cells, a very small sample, the universality and applicability of our findings is impaired, and further experimental proofs are needed to validate the findings.

Conclusions
In conclusion, this study identified several potential molecular targets that might contribute to the MTX resistance in osteosarcoma cells, such as GADD45A, AARS, AURKA, AURKB, CENPA, CCNB1, CCNE2, and CDK1, which may function via aminoacyl-tRNA biosynthesis pathway, cell cycle pathway, or p53 signaling pathway. However, the finding here should be taken prudently.