This article is mentioned in:

Introduction

Osteosarcoma (OS), which is one of the common primary bone tumors in children and adolescents is highly aggressive and metastatic (1); it comprises ~5% of all pediatric tumors and patients with pulmonary metastasis present a poor prognosis (2). There are definite differences in the mechanism between primary and metastatic OS, and numerous previous studies have been performed to expound them (3–5).

Multiple signaling pathways, including the Notch signaling pathway, are involved in cancer migration initiation through the extracellular matrix (ECM) towards the bloodstream (6). Survival of cancer cells in the bloodstream is mediated by the integrin signaling pathway to activate the phosphatidylinositol 3-kinase/protein kinase B survival pathway (7). Signaling pathways, including those of Src and Wnt, are involved in apoptosis resistance throughout metastasis (8,9). In addition, the immune system is also invaded and modulated for survival advantage (10). Subsequently, cell adhesion is arrested, and the cells extravasate to their new environment and re-establish cell-cell adhesions (11). There are numerous various viewpoints regarding the mechanism underlying metastatic tumor cell arrest to distant organ sites (12,13). Angiogenesis is necessary for cancer cell proliferation and metastasis (12); however, there remain a number of unknown mechanisms, including the precise mechanism underlying how cells alternate between exhibiting primary and metastatic properties.

Therefore, in the present study, gene expression profiles were used to further investigate the metastatic mechanism of OS. The GSE49003 dataset considered in the present study was previously analyzed by Diao et al (14) to identify differentially expressed genes (DEGs) between metastatic and non-metastatic patients with OS, and crucial microRNAs associated with OS metastasis, by merging data from different metastatic or non-metastatic OS cell lines. However, different types of metastatic OS may be regulated by different molecular mechanisms. In addition, transcription factors (TFs) may also serve a vital role in this pathomechanism. Therefore, the dataset was reanalyzed in the present study to emphasize the different mechanisms of different metastatic OS cell lines.

Materials and methods

Gene expression profile data

The raw expression data (dataset GSE49003; http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE49003), as provided by Endo-Munoz et al on July 18, 2013, were used in the present study. The microarray expression profile was obtained from two metastatic OS cell lines and two non-metastatic OS cell lines. The metastatic KHOS and KRIB cell lines and the non-metastatic HOS cell line were used with three-duplicated samples, and gene expression data from each of these cell lines were used. The platform of this dataset was GPL6847 Illumina HumanHT-12 V3.0 expression beadchip (Illumina Inc., San Diego, CA, USA).

Data preprocessing and DEG screening

The Limma package (15) in Bioconductor was used to add probe annotation files for each Illumina chip in order to preprocess the expression profile. Background correction, quantile normalization and probe summarization were performed to generate the gene expression data matrix.

DEGs between KHOS and HOS and between KRIB and HOS were determined using the Limma package. The differential expression of genes were analyzed by Student's t-test, and those with a false discovery rate adjusted P-value of <0.01 and |fold change| ≥2 were screened.

Functional enrichment analysis of DEGs

The Gene Ontology (GO) (16) project was established for gene classifications by molecular function, biological process (BP) and cellular component. DEGs of KHOS vs. HOS and KRIB vs. HOS were functionally enriched by GO. In addition, the Kyoto Encyclopedia of Genes and Genomes (KEGG) (17), which has a become a major database resource for understanding high-level functions of genes, was utilized. The default threshold of P<0.01 was selected for the hypergeometric enrichment test.

Protein-protein interaction (PPI) network construction and subnetwork mining

The Search Tool for the Retrieval of Interacting Genes/Proteins is a biological database that provides known and predicted PPIs (13). The tool was applied in the present study to identify interacting protein pairs between DEGs with a PPI score of 0.9. Subsequently, Cytoscape software version 2.8.0 (18) was used to visualize the constructed PPI network. Subnetworks (modules) with a hypergeometric P-value <0.05 were identified by the ClusterONE plugin (19) from Cytoscape. Furthermore, the Database for Annotation, Visualization and Integrated Discovery (20,21) was utilized to perform the KEGG pathway cluster analyses of DEGs in modules with P<0.05.

Overlapping gene analysis

Overlapping DEGs that were upregulated in KHOS and KRIB cells compared with HOS cells were identified, and downregulated DEGs that were common in KHOS and KRIB cells were identified. Thereafter, KEGG signaling pathways of these two types of overlapping genes were enriched.

Furthermore, based on the regulatory association between TFs and target genes recorded in the University of California Santa Cruz (UCSC) (22) database, regulatory associations between TFs and their target DEGs were identified.

Results

DEGs of various groups

A total of 1,552 (711 downregulated and 841 upregulated) and 1,330 DEGs (570 downregulated and 760 upregulated) were obtained from the KHOS vs. HOS and KRIB vs. HOS comparisons, respectively.

Functional enrichment analyses of DEGs

Significant enriched terms of GO and KEGG pathway enrichment analyses in KHOS and KRIB groups are presented in Tables I and II, respectively. Upregulated genes of KHOS were associated with GO-BP terms of anatomical structure morphogenesis and cellular response to extracellular stimulus, and the downregulated genes were enriched in BP terms of multicellular organismal development and anatomical structure morphogenesis. Using KEGG, focal adhesion, axon guidance and peroxisome proliferator-activated receptor signaling pathways were significantly enriched by upregulated genes of KHOS, and pathways in cancer, tight junctions and cell adhesion molecules were enriched by downregulated DEGs in KHOS cells (Table I). GO and pathway enrichment analyses of KRIB were similar to those of KHOS, and the significantly enriched GO term of KRIB overexpressed genes was anatomical structure morphogenesis, and the significantly enriched pathway was focal adhesion (Table II).

Table I. Gene ontology and pathway enrichment analysis of genes in KHOS cells.

Table I.

Gene ontology and pathway enrichment analysis of genes in KHOS cells.

Category identity	Name	Count	P-value
Upregulated genes
GO:0009653	Anatomical structure morphogenesis	132	8.44×10−7
GO:0031668	Cellular response to extracellular stimulus	20	1.75×10−6
GO:0031669	Cellular response to nutrient levels	18	1.98×10−6
GO:0071496	Cellular response to external stimulus	24	2.72×10−6
GO:0072358	Cardiovascular system development	61	4.90×10−6
Pathway:4510	Focal adhesion	19	5.07×10−4
Pathway:4360	Axon guidance	14	7.61×10−4
Pathway:250	Alanine, aspartate and glutamate metabolism	6	1.66×10−3
Pathway:3320	PPAR signaling pathway	9	2.08×10−3
Pathway:4512	ECM-receptor interaction	10	2.35×10−3
Downregulated genes
GO:0007275	Multicellular organismal development	310	0.00
GO:0009653	Anatomical structure morphogenesis	192	0.00
GO:0032502	Developmental process	338	0.00
GO:0044767	Single-organism developmental process	285	0.00
GO:0048468	Cell development	158	0.00
Pathway:5200	Pathways in cancer	35	2.61×10−5
Pathway:5146	Amoebiasis	14	1.07×10−3
Pathway:5144	Malaria	9	1.08×10−3
Pathway:4530	Tight junction	16	1.24×10−3
Pathway:4514	Cell adhesion molecules	16	1.34×10−3

[i] Count, number of genes enriched in the functional term; GO, Gene Ontology; PPAR, peroxisome proliferator activated receptor; ECM; extracellular matrix.

Table II. Gene ontology and pathway enrichment analysis of genes in KRIB cells.

Table II.

Gene ontology and pathway enrichment analysis of genes in KRIB cells.

Category identity	Name	Count	P-value
Upregulated genes
GO:0009653	Anatomical structure morphogenesis	106	3.88×10−6
GO:0001944	Vasculature development	38	9.85×10−6
GO:0010646	Regulation of cell communication	108	1.24×10−5
GO:0071294	Cellular response to zinc ion	5	1.50×10−5
GO:0023051	Regulation of signaling	107	1.88×10−5
Pathway:4510	Focal adhesion	16	7.36×10−4
Pathway:480	Glutathione metabolism	7	1.05×10−3
Pathway:4360	Axon guidance	11	3.03×10−3
Pathway:3320	PPAR signaling pathway	7	7.33×10−3
Pathway:4520	Adherens junction	7	9.19×10−3
Downregulated genes
GO:0000902	Cell morphogenesis	104	0.00
GO:0000904	Cell morphogenesis involved in differentiation	86	0.00
GO:0006928	Cellular component movement	132	0.00
GO:0007275	Multicellular organismal development	299	0.00
GO:0007399	Nervous system development	153	0.00
Pathway:4514	Cell adhesion molecules	19	1.14×10−5
Pathway:5144	Malaria	11	1.74×10−5
Pathway:5146	Amoebiasis	16	2.83×10−5
Pathway:5200	Pathways in cancer	32	4.49×10−5
Pathway:4512	ECM-receptor interaction	13	1.40×10−4

[i] Count, number of genes enriched in the functional term; GO, Gene Ontology; PPAR, peroxisome proliferator activated receptor; ECM; extracellular matrix.

PPI network and subnetwork

A total of 856 pairs of interacting proteins and 495 nodes were noted in the PPI network of KHOS (Fig. 1A). The top 6 nodes with a greater degree of connectivity compared with the others were jun proto-oncogene (JUN; degree=33), vascular endothelial growth factor A (VEGFA; degree=29), epidermal growth factor receptor (EGFR; degree=23), cyclin-dependent kinase inhibitor 1A (degree=19), interleukin 6 (IL-6; degree=19) and signal transducer and activator of transcription 3 (STAT3; degree=17). Furthermore, two modules with P<0.05 were obtained (Fig. 1B and C). A total of 9 genes, including 7 collagen molecules and 2 integrin family members involved in module 1 were significantly enriched in terms of ECM-receptor interaction and focal adhesion.

Figure 1. PPI network and subnetwork (module) of DEGs in KHOS cells. (A) PPI network. (B) Module 1. (C) Module 2. The diamond nodes denote the downregulated DEGs. The circular nodes represent the upregulated DEGs. Lines between nodes reveal the interactions between them. PPI, protein-protein interaction; DEGs, differentially expressed genes.

There were 719 interacting pairs and 397 nodes of the PPI network of KRIB (Fig. 2). VEGFA (degree=29), Fyn oncogene related to Src, FGR, Yes (degree=25), JUN (degree=23), EGFR (degree=22) and IL-6 (degree=19), were the top 5 hub nodes. In addition, one module was screened (Fig. 2), and the significant pathways enriched by its genes were ECM-receptor interaction and focal adhesion.

Figure 2. PPI network and subnetwork (module) of DEGs in KRIB cells. (A) PPI network. (B) Module one. The diamond nodes denote the downregulated DEGs. The circular nodes represent the upregulated DEGs. Lines between nodes reveal the interaction between them. PPI, protein-protein interaction; DEGs, differentially expressed genes.

KEGG pathway enrichment and transcriptional regulatory network analysis of overlapping DEGs

A total of 421 upregulated and 595 downregulated overlapping genes were obtained. Pathway of focal adhesion, glutathione metabolism and endocytosis were enriched by their upregulated genes, whereas pathways in cancer, focal adhesion and cell adhesion molecules were enriched by downregulated genes. In addition, the mitogen-activated protein kinase (MAPK) and transforming growth factor-β (TGF-β) signaling pathways were enriched by overlapping genes with decreased expression levels (Table III).

Table III. Kyoto Encyclopedia of genes and genomes pathway of overlapping genes.

Table III.

Kyoto Encyclopedia of genes and genomes pathway of overlapping genes.

Category identity	Name	Count	P-value
Upregulated overlapping gene
4510	Focal adhesion	13	2.96×10−3
480	Glutathione metabolism	5	3.11×10−2
4144	Endocytosis	10	3.18×10−2
4360	Axon guidance	8	3.46×10−2
4520	Adherens junction	6	3.67×10−2
Downregulated overlapping gene
5200	Pathways in cancer	27	2.52×10−4
4510	Focal adhesion	17	3.93×10−3
4514	Cell adhesion molecules	13	4.25×10−3
4512	ECM-receptor interaction	10	4.40×10−3
4010	MAPK signaling pathway	19	1.26×10−2
4350	TGF-β signaling pathway	9	1.72×10−2

[i] Count, number of genes enriched in the functional term; ECM, extracellular matrix; MAPK, mitogen-activated protein kinase; transforming growth factor-β.

Fig. 3 demonstrates the transactional regulatory network of overlapping DEGs. A total of 7 TFs, including activating transcription factor 3, forkhead box (FOX)A1, STAT3, FOXA2, hes-related family bHLH transcription factor with YRPW motif 1 (HEY1), SUZ12 polycomb repressive complex 2 subunit and zinc finger ZZ-type containing 3 were identified. The first three of these were downregulated and the others were overexpressed.

Figure 3. Transcription regulatory network of overlapping DEGs. The diamond nodes denote the downregulated DEGs. The circular nodes represent the upregulated DEGs. The triangular nodes reveal the transcription factors. DEGs, differentially expressed genes.

Discussion

Due to the high risk for recurrence and the poor survival of metastatic OS patients (23), two different metastatic OS cell lines and one non-metastatic cell lines were used to investigate the mechanism underlying metastasis. The GO and pathway enrichment analysis results were similar of these two metastatic types of OS. Besides, a large number of overlapping genes, including 7 TFs, were obtained. These genes and biological functions may be important in OS tumor metastasis.

The most prominent node in the PPI networks of KHOS and KRIB was VEGFA, which is a growth factor that is activated in vasculogenesis, angiogenesis and endothelial cell growth (24). It is not only a commonly upregulated angiogenic factor during tumor growth, but also the most frequent lymphangiogenic factor contributing to the lymphatic metastasis (25). A previous study suggested that VEGFA promoted tumor growth and metastasis via the VEGF-VEGFR1 signaling pathway (26). In addition, EGFR, which is a member of the EGF family, was another common hub node with a high degree (27). It was revealed that EGFR participated in the promotion of cell division, migration and angiogenesis, and inhibited cell apoptosis (28). Furthermore, these two genes were enriched in the development and differentiation process, including anatomical structure morphogenesis and vasculature development. Yang et al (29) demonstrated that Twist, which is a master regulator of embryonic morphogenesis, contributed to metastasis by promoting an epithelial-mesenchymal transition (29). Therefore, tumor metastasis is the recapitulation of morphogenetic processes (30). Morphogenetic factors, including hepatocyte growth factor/scatter factor have metastasis-promoting effects on breast carcinoma cells (31). Based on these studies, VEGFA and EGFR may be factors involved in the morphogenetic processes to promote tumor metastasis.

The KEGG pathway enrichment analysis results of module 1 of the KHOS group was similar to that of module 1 of the KRIB group. ECM-receptor interaction and focal adhesion were the observed pathways enriched by multiple collagen genes and integrin genes. These results were in accordance with previous studies (32,33) and demonstrated the accuracy and practicability of them. Additionally, functional enrichment analysis of overlapping genes of these two group would be more direct to exhibit commonality of metastatic OS.

Focal adhesion was the significant pathway enriched by upregulated and downregulated overlapping genes. Furthermore, these downregulated genes remained enriched in signaling pathways that included the MAPK and TGF-β signaling pathways. The former pathway serves an important role in the early metastasis of the tumor (34,35), whereas the latter has positive and negative effects on carcinogenesis (36), and further investigation for is required to fully elucidate the molecular mechanism.

Finally, a total of 7 TFs, which regulated multiple genes, were screened. Among them, FOXA1, STAT3 and HEY1 were more important in the regulatory network. Candy et al (37) indicated that overexpressed HEY1 was a negative prognostic factor for colorectal cancer, and it is closely associated with vascular and perineural invasion and lymph node metastasis. Therefore, the upregulated HEY1 may be associated with OS metastasis. In addition, FOXA1 was revealed to be a participant in embryonic development, mediating tissue-specific gene expression in differentiated tissues (38). This gene encoded protein may be another important factor involved in the similar process of morphogenesis.

In conclusion, recapitulation of morphogenetic processes, a disordered signaling system, and activated cell adhesion and migration were involved in the facilitated tumor metastasis. Furthermore, morphogenesis and vasculature development induced angiogenesis, and tumor growth and adhesion, or even the lymph node metastasis induced the metastatic OS. The present study further illustrated the OS tumorigenesis mechanism and displayed the crucial roles of genes such as VEGFA, EGFR and HEY1. The practical clinical significance has not yet been clearly described and requires further investigation.

References