Screening for Metastatic Osteosarcoma Biomarkers with a DNA Microarray

Osteosarcoma (OS) is a malignant bone tumor commonly developing in children, adolescents, and young adults (Hu et al., 2013). Besides, OS has high metastatic potential (about 15%-20%), especially for pulmonary metastasis, which is a major cause of death in OS patients (Salinas-Souza et al., 2013; Zhou et al., 2013). It has been reported that the 5-year survival rate for non-metastatic patients is 60-70% while for pulmonary metastasis patients, far worse, only less than 30% (Ferrari et al., 2005). Therefore, it is of great importance to investigate novel therapeutic regimens to prevent pulmonary metastases during the early stage, and further to improve the prognosis of OS patients. As reported, many oncogenes and suppressor genes are abnormally expressed in OS tissue and many genes are associated with the metastases of OS (Ragland et al., 2002). Densmore et al. have indicated that p53 can effectively reduce the number and the size of lung metastases from OS tissue, and also alleviate the systemic toxicity by aerosol gene therapy with PEI-p53 complexes (Densmore et al., 2001). Further, the expression of Fas inversely correlates with the metastatic potential of OS to


Introduction
Osteosarcoma (OS) is a malignant bone tumor commonly developing in children, adolescents, and young adults (Hu et al., 2013).Besides, OS has high metastatic potential (about 15%-20%), especially for pulmonary metastasis, which is a major cause of death in OS patients (Salinas-Souza et al., 2013;Zhou et al., 2013).It has been reported that the 5-year survival rate for non-metastatic patients is 60-70% while for pulmonary metastasis patients, far worse, only less than 30% (Ferrari et al., 2005).Therefore, it is of great importance to investigate novel therapeutic regimens to prevent pulmonary metastases during the early stage, and further to improve the prognosis of OS patients.
As reported, many oncogenes and suppressor genes are abnormally expressed in OS tissue and many genes are associated with the metastases of OS (Ragland et al., 2002).Densmore et al. have indicated that p53 can effectively reduce the number and the size of lung metastases from OS tissue, and also alleviate the systemic toxicity by aerosol gene therapy with PEI-p53 complexes (Densmore et al., 2001).Further, the expression of Fas inversely correlates with the metastatic potential of OS to

Screening for Metastatic Osteosarcoma Biomarkers with a DNA Microarray
Chun-Yu Diao 1 , Hong-Bing Guo 1 , Yu-Rong Ouyang 1 , Han-Cong Zhang 1 , Li-Hong Liu 1 , Jie Bu 1 , Zhi-Hua Wang 2 , Tao Xiao 1 * the lung and the inhibition of c-FLIP (an inhibitor of the Fas pathway), which may be a possible therapeutic target for the treatment of lung metastasis (Rao-Bindal et al., 2012).Ezrin, encoded by the EZR gene, has been reported to be associated with the progression and metastasis of OS based on its role in physically and functionally connecting the actin cytoskeleton (F-actin) to the cell membrane (Ren et al., 2008;Di Cristofano et al., 2010).The expression of genes will determine the occurrence, development, metastasis and prognosis of OS, therefore, gene therapy may supply new sights for OS patients to effectively prevent and resist the metastasis.
Besides, microRNAs (miRNAs) are pivotal components of gene-regulatory networks which modulate a great deal of important pathophysiological processes, including the initiation and progression of cancers (ZhouShi et al., 2013).In OS patients, many miRNAs are also in abnormal status which then lead to the abnormal expression of related genes, for example, the down-regulation of miRNA-143 promotes celullar apoptosis and suppresses the development of OS by targeting Bcl-2 (Zhang et al., 2010).Moreover, miRNA-20a can down-regulate the expression of Fas in OS, and further contributes to the metastasis of OS cells by altering the phenotype and allowing survival in the FasL (+) lung microenvironment (Huang et al., 2012).MiR-206 may has therapeutic effect for the metastasis of OS cells mainly through reducing the cellular viability, promoting apoptosis, and inhibiting cell invasion and migration (Bao et al., 2013).
With a rapid expansion of our knowledge about gene therapy, emerging evidences suggest that the metastasis of OS generally is caused by the variation of involved genes and miRNAs.In order to discover the key factors influencing the metastasis of OS, we employed the high throughput method to identify new targets that can be used as the biomarkers and therapeutic approaches of metastatic OS.

Microarray data
The complete data set was deposited in the National Center of Biotechnology Information (NCBI) Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo,GEO series accession number GSE49003) which was based on GPL6947 platform data (Illumina Human HT-12 V3.0 expression beadchip).Total 12 chips were divided into 2 groups: metastatic OS group (M-group, n=6) and non-metastatic OS group (N-group, n=6).We downloaded the raw data and the probe annotation files for further analysis.

Data preprocessing and differentially expressed genes (DEGs) screening
Expression data in TXT document were first matched with the corresponding genes according to their names and then normalized by quantile method (Smythand Speed, 2003).To identify DEGs between M-group and N-group, the limma package (Linear Models for Microarray Data) (Smyth, 2005) in R language was applied in our analysis.Genes that were consistently different from non-metastatic OS tissue were selected as DEGs.P-value of the paired t-test < 0.05, the False Discovery Rate of P-value (FDR) < 0.05 and |log (fold change)| > 1 were chosen as the cutoff criterion, which were calculated using Bonferroni method in multtest package (Benjaminiand Hochberg, 1995;Benjamini, 2010).

Comparison between DEGs of metastatic and nonmetastatic OS
To further investigate the differences between metastatic and non-metastatic OS, the expression values of all selected DEGs in the two groups were compared using t-test and the p-value less than 0.05 was regarded as the threshold.

Functional enrichment analysis
Functional enrichment analysis of genes was widely used in microarray data analysis to uncover the functional changes of large-scale genes (Hwang et al., 2006).DAVID (the Database for Annotation, Visualization and Integrated Discovery), a high throughput and integrated data-mining environment, was the most popular tool to systematically extract biological meaning from large gene or protein lists (Tanay et al., 2004;Huang da et al., 2009).In the current study, we grouped the DEGs as up-regulated group and down-regulated group, and then uploaded the gene groups for the functional enrichment analysis using DAVID system.A p-value less than 0.05 was chosen as the cutoff value.

Protein-Protein interaction (PPI) network construction
Osprey was a widely used network visualization system which could not only represent interactions in a flexible and rapidly expandable graphical format, but also provide options for functional comparisons between datasets (Breitkreutz et al., 2003b).The Osprey software integrated the General Repository for Interaction Datasets (GRID) (Breitkreutz et al., 2003a) and Biomolecular Interaction Network Database (BIND) (Willisand Hogue, 2006) to investigate protein-protein interaction networks and protein complex conveniently.Therefore, we applied Osprey in the current study and collected more than 50, 000 pairs of PPI to construct the interaction network.

Construction of integrated miRNA-DEG regulatory network
MiRNAs have the capability of affecting cellular transformation, acting either as oncogenes or tumor suppressors, and also exhibit differential expression levels in cancer (Zhang et al., 2007).Therefore, we constructed miRNA-DEG regulatory network to investigate the modulatory functions of the significant DEGs-related miRNAs on screened DEGs.First, the miRNAs related with the up-regulated and down-regulated DEGs were separately selected using the Web-based Gene Set  Analysis Toolkit (WebGestalt; http://bioinfo.vanderbilt.edu/webgestalt/) (Zhang et al., 2005;Duncan et al., 2010).Then the integrated regulatory network was constructed based on the selected significantly related miRNAs and PPI network.A p-value less than 0.05 was chosen as the cutoff value for the significantly related miRNAs screening.

DEGs screening
According to the statistical analysis of microarray data of metastatic and non-metastatic OS, a total of 323 DEGs (|log FC| > 1) including 134 up-regulated DEGs and 189 down-regulated DEGs were screened out (Table 1).The expression profile data after standardization were shown in Figure 1A.

Comparison between DEGs of metastatic and nonmetastatic OS
As shown in Figure 1B, the expression values of DEGs in M-group was significantly higher than those in N-group (P = 0.01388).The expression of DEGs of metastatic OS tended to over-express.

Functional enrichment analysis
The functional enrichment analyses were performed for both up-and down-regulated DEGs.The up-regulated DEGs were mainly enriched in 14 subcategories and most significantly in cytoskeleton organization, while the down-regulated DEGs were in 13 subcategories and most significantly connected with wound healing.The number of DEGs participating in the two enriched functions was 16 (19%, 16/134) and 10 (5.29%, 10/189), respectively.

The significant DEGs-related miRNAs
The screening of significant DEGs-related miRNAs for both up-and down-regulated DEGs was performed by the WebGestalt software.As shown in Table 2, totally 10 miRNAs were selected as the significant DEGsrelated miRNAs for the up-regulated DEGs and the most significant one was miR-9, while for the down-regulated DEGs, also 10 miRNAs were screened and the most significant one was miR-202.

Construction of integrated miRNA-DEG regulatory network
The integrated miRNA-DEG regulatory network was constructed by integrating a large number of DEGs and DEG-related miRNAs.The integrated regulatory network also constructed based on the up-and down-regulated DEGs.The regulatory network was shown in Figure 2.

Discussion
In the current study, we aimed to extend our understanding of the differences of metastatic and nonmetastatic OS, further to investigate candidate biomarkers for the prediction of metastatic OS.Based on the above experiment, totally 323 DEGs including 134 up-regulated DEGs and 189 down-regulated DEGs were screened out.In addition, the functional enrichment showed that down-regulated DEGs may be highly associated with cytoskeleton organization, cytoskeleton and proteinaceous extracellular matrix; while upregulated DEGs with wound healing, cell motion and enzyme binding.MiRNAs can regulate the expression of corresponding genes generally through down-regulated translation of target mRNA or up-regulated degradation of target mRNA (Hoffman et al., 2013).To understand the regulatory network of these DEGs, we stretched back to the related miRNAs, and integrated their possible relationships using high throughput methods.Based on bioinformatic algorithms, we identified 10 miRNAs that were significantly related with the down-regulated DEGs, including miR-202, miR-381, miR-98, miR-29, miR-506, miR-495, miR-484, miR-23, miR-182 andmiR-18; and also 10 miRNAs that were significantly related with the down-regulated DEGs, including miR-9, miR-101, miR-429, miR-520D, miR-142-5P, miR-330, miR-525, miR-15, miR-524 and miR-22.Of these miRNAs, miR-29, miR-142-5P and miR-15 are verified to be highly related to OS (Jones et al., 2012), while miR-202 and miR-9, which are both significantly related miRNAs, may be also associated with OS according to our results.
MiR-202 is a short RNA molecular which locates within a chromosomal fragile site in 10q26 and many studies have clarified the association between miR-202 and various tumors, such as breast cancer, colorectal cancer, follicular lymphoma and endometrial and brain tumors (Zhao et al., 2013).In our study, miR-202 was also high-related with the metastasis of OS, which keeps in line with the former surveys.Further, CALD1, regulated by miR-202 according to the Figure 2A, plays a key role in the regulation of microfilament organization, consequently controlling cell shape, adhesion, cytokinesis and motility (Cuomo et al., 2005).Also, it has been reported to participate in splicing events in cancer (Thorsen et al., 2008).Besides, CALD1 is a cytoskeletonassociated protein and CALD1 is significantly enriched in cytoskeleton organization and cytoskeleton functions and other 4 functions in our experiment, which further indicates the validity of our experiment.Cytoskeleton, the fibrin mesh structure present in cytoplasm and nuclei of all cells, not only forms a continuous dynamic connection between cellular structure for keeping the cell morphology and structure, but also participates in the material exchange and transport out of and into the cells, and even a series of important life activities and the pathological and physiological process like tumor metastasis, cell movement (Amsellem et al., 2005).Therefore, we infer that the metastasis of OS may be attributed to the downregulation of miR-202, which consequently leads to the low expression of CALD1, and further influences the changes in cytoskeleton and cellular adhesion.MiR-202 may be a potential key regulator resulting in the metastasis of OS, and CALD1 may be the possible therapeutic target.MiR-9 is also a widely studied tumor suppressor miRNA in acute myelocytic leukemia, nasopharyngeal carcinoma and ovarian cancer, even an essential oncogenic miRNA specifically over-expressed in mixed lineage leukemia-rearranged leukemia (Emmrich et al., 2013;Lu et al., 2013;Sun et al., 2013).Furthermore, miR-9 is reported to be up-regulated in breast cancer cells, influencing the expression of E-cadherin, and further leading to increased cell motility and invasiveness and even cancer metastasis (Khew-Goodalland Goodall, 2010;Ma et al., 2010).Based on these conclusions, we speculate miR-9 may play a critical role in the development and metastasis of OS, and consistent with the speculation, miR-9 was also identified as a significantly related miRNA contributing to the metastasis of OS in the current study.In addition, STX1A is a new screened DEG as a new target of the metastatic OS compared with non-metastatic OS, and it is directly regulated by miR-9 as shown in Figure 2B.STX1A can be useful for prognosing lung cancer survival and early-stage non-small-cell lung cancer (NSCLC) (Lau et al., 2007;Tsao et al., 2011).Moreover, STX1A expressed at a low level in adenocarcinoma, but a high level in squamous cell lung carcinomas (Fan et al., 2012).However, there is no direct or indirect evidence to verify whether STX1A affects cancer metastasis, and even the functions in cell are still not clearly described.Based on the results of functional enrichment in the current study, STX1A was enriched in several functions involving protein dimerization activity, cytoplasmic vesicle, vesicle, endomembrane system and plasma membrane.All the functions are associated with membrane stability, thus, STX1A may influence cancer metastasis through changing the membrane stability in metastatic OS, and the activity may be directly regulated by miR-9.Thus, miR-9 and STX1A gene may be the key factors influencing OS In conclusion, we have identified a total of 323 DEGs including 134 up-regulated DEGs and 189 down-regulated DEGs in metastatic OS compared with non-metastatic OS.This set of DEGs might be used as biomarkers for early diagnosis or therapy targets for metastatic OS treatment.In addition, we identified two important miRNAs (miR-202 and miR-9) pivotal for the metastasis of OS, and their relevant genes, CALD1 and STX1A.However, the whole study was carried out based the bioinformatics methods and the conclusions have not been confirmed by corresponding experiments yet.Therefore, further experiments are in great need to verify the effect and mechanism of miR-202, miR-9, CALD1and STX1A in metastatic OS.

Figure 1
Figure 1.A: Box figure of DEGs Expression Values in Osteosarcoma Tissue Samples after Standardization.B:The comparison between DEGs of metastatic and non-metastatic OS.The white boxes represented the non-metastatic OS samples and the grey ones were metastatic ones.The X axis was sample origin while the Y axis was expression level of DEGs.The black line in the center was the median of expression value and the consistent distribution indicated a good standardization.*P < 0.05, compared with non-metastatic OS group