This article is mentioned in:

Introduction

Osteosarcoma (OS) is the most common non-haematological primary bone tumour that affects children and young adolescents, and pulmonary metastases are the major cause of mortality (1–3). In total, ~20% of patients present with pulmonary metastasis at initial diagnosis, and in ~40% of patients, metastases occur in the advanced stage (4). Furthermore, OS has a low survival rate in patients due to chemoresistance and the high rates of pulmonary metastasis (5). For patients with non-metastatic OS, 5-year OS has increased to 60–70%, but this rate is reduced to 20% when metastases occur (6). Therefore, strategies to prevent tumour metastasis are urgently needed. However, the molecular mechanisms of pulmonary metastasis in patients with OS remain poorly understood. Thus, identifying key molecules associated with pulmonary metastasis of OS is important for the development of more effective metastasis-suppressive therapies.

During the last decade, gene expression profiling has identified critical genes and cellular signaling pathways related to the occurrence, development and change of human OS (7). Microarray techniques combined with bioinformatics analysis can determine the differential expression levels of genes accurately and provide an effective tool for large-scale gene expression studies (8). In the patients with OS, the lung is the most common metastatic site, and metastases can lead to a high rate of mortality (9). Therefore, examining differentially expressed genes (DEGs) between non-metastatic and metastatic OS samples is helpful for identifying the key genes and pathways leading to pulmonary metastasis.

In the current study, by comparing the gene expression of non-metastatic and metastatic OS samples in the GEO database, the authors identified DEGs that may be related to pulmonary metastasis. Subsequently, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment pathway analyses were performed. In combination with protein-protein interaction (PPI) information, not only did the authors identify relevant genes and pathways, but also revealed existing molecular mechanisms. In conclusion, the authors' finding can improve the understanding of OS and identify the genes associated with pulmonary metastasis, thereby providing potential therapeutic targets for further studies.

Materials and methods

Gene expression microarray data

In the present study, the gene expression profiles of GSE85537 were downloaded from Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo/). GSE85537 was based on the Affymetrix GPL570 platform (Affymetrix Human Genome U133 Plus 2.0 Array; Affymetrix; Thermo Fisher Scientific, Inc., Waltham, MA, USA). The GSE85537 dataset contained six samples, including three non-metastatic OS samples and three metastatic OS samples.

Identification of DEGs

The raw data files used for the analysis included TXT files. The analysis was carried out using GEO2R (https://www.ncbi.nlm.nih.gov/geo/geo2r/), which can perform comparisons on original submitter-supplied processed data tables using the GEOquery and limma R packages from the Bioconductor project. The DEGs between the non-metastatic and metastatic OS samples were selected (P<0.05), and overlapped genes were identified.

GO enrichment and KEGG pathway analysis of the DEGs

After obtaining the DEGs, the authors submitted the DEG list to the online software Database for Annotation, Visualization and Integrated Discovery (DAVID, https://david.ncifcrf.gov/) to identify overrepresented GO categories and pathway categories. GO analysis can determine the biological meaning in a large list of genes and categorize gene product functions, including biological process (BP), molecular function (MF) and cellular component (CC) (10,11). KEGG (http://www.genome.jp/) is a knowledge base for systematic analysis of gene functions, linking genomic information with higher-level systemic functions (12,13). Finally, the enriched functions of DEGs were selected via GO and KEGG pathway analysis, and P<0.05 was considered to indicate a statistically significant difference.

Construction of the PPI Network of DEGs

To further investigate the molecular mechanism of pulmonary metastasis in patients with OS, the Search Tool for the Retrieval of Interacting Genes (STRING) database (http://www.string-db.org/) was used to evaluate the interactive relationships among DEGs. Initially, the DEG list was first submitted to STRING, and then, the experimentally validated interactions were selected with a combined score >0.4. Subsequently, the PPI networks were analyzed using Cytoscape software (version 3.5.1; www.cytoscape.org). Then, the plug-in Molecular Complex Detection (MCODE) was used to screen the modules of the PPI network in Cytoscape. Furthermore, the enrichment analyses were performed for DEGs in the corresponding modules. P<0.05 was considered to indicate a statistically significant difference.

Results

Identification of DEGs

The gene expression profile GSE85537 was downloaded from the GEO database, and the GEO2R method was used to identify DEGs in metastatic OS samples compared with non-metastatic OS samples. P<0.05, logFC (fold control) >2.0 or logFC <-2.0 was used as the criteria, and 2,493 genes were identified as DEGs. Among these, 485 genes (19.45%) were upregulated, and the remaining 2,008 genes (80.55%) were downregulated. Subsequently, the top 50 upregulated and downregulated DEGs were selected to generate the heatmap (Fig. 1).

Figure 1. The heat map of the differentially expressed genes (top 50 upregulated and downregulated genes). Red, upregulation; blue, downregulation.

GO term enrichment analysis

To functionally categorize these 2,493 significant genes, the DEGs were analyzed using the online software DAVID. GO analysis revealed that the upregulated DEGs were significantly enriched in BP, including transmembrane transport, ion transport, and release of sequestered calcium ion into cytosol (Table I); the downregulated DEGs were significantly enriched in BP, including cell surface receptor signaling pathway, actin filament-based process and movement of cell or subcellular component (Table I). For MF, the upregulated DEGs were enriched in ATPase-coupled ion transmembrane transporter activity and cation-transporting ATPase activity, and the downregulated DEGs were enriched in Ras guanyl-nucleotide exchange factor activity, guanyl-nucleotide exchange factor activity, and growth factor binding (Table I). In addition, GO CC analysis also indicated that the upregulated DEGs were significantly enriched in the synapse, integral component of plasma membrane and intrinsic component of plasma membrane, and downregulated DEGs were enriched in the proteinaceous extracellular matrix, extracellular matrix and plasma membrane region (Table I).

Table I. Gene ontology analysis of DEGs.

Table I.

Gene ontology analysis of DEGs.

Expression	Category	GO-ID	Term	Gene count	%	P-value
Upregulated	BP	GO:0055085	Transmembrane transport	43	12.6	5.10×10−6
	BP	GO:0006811	Ion transport	45	13.2	1.40×10−5
	BP	GO:0051209	Release of sequestered calcium ion into cytosol	10	2.9	3.20×10−5
	BP	GO:0051283	Negative regulation of sequestering of calcium ion	10	2.9	3.20×10−5
	BP	GO:0097553	Calcium ion transmembrane import into cytosol	10	2.9	3.50×10−5
	MF	GO:0042625	ATPase coupled ion transmembrane transporter activity	7	2.0	1.10×10−3
	MF	GO:0019829	Cation-transporting ATPase activity	6	1.8	2.90×10−3
	MF	GO:0022853	Active ion transmembrane transporter activity	7	2.0	6.10×10−3
	MF	GO:0042626	ATPase activity, coupled to transmembrane movement of substances	7	2.0	7.20×10−3
	MF	GO:0070330	Aromatase activity	4	1.2	7.50×10−3
	CC	GO:0045202	Synapse	30	8.8	5.30×10−6
	CC	GO:0005887	Integral component of plasma membrane	46	13.5	5.10X10−5
	CC	GO:0031226	Intrinsic component of plasma membrane	46	13.5	1.30×10−4
	CC	GO:0098794	Postsynapse	17	5.0	1.50×10−4
	CC	GO:0044456	Synapse part	23	6.7	1.50×10−4
Downregulated	BP	GO:0007166	Cell surface receptor signaling pathway	226	17.2	2.00×10−9
	BP	GO:0030029	Actin filament-based process	79	6.0	4.30×10−9
	BP	GO:0006928	Movement of cell or subcellular component	165	12.5	4.30×10−9
	BP	GO:2000145	Regulation of cell motility	85	6.5	1.10×10−8
	BP	GO:0040012	Regulation of locomotion	87	6.6	1.70×10−8
	MF	GO:0005088	Ras guanyl-nucleotide exchange factor activity	34	2.6	1.10×10−6
	MF	GO:0005085	Guanyl-nucleotide exchange factor activity	41	3.1	1.30×10−6
	MF	GO:0019838	Growth factor binding	23	1.7	4.00×10−6
	MF	GO:0005102	Receptor binding	126	9.6	4.90×10−6
	MF	GO:0004713	Protein tyrosine kinase activity	26	2.0	4.20×10−5
	CC	GO:0005578	Proteinaceous extracellular matrix	55	4.2	3.90×10−10
	CC	GO:0031012	Extracellular matrix	67	5.1	1.10×10−8
	CC	GO:0098590	Plasma membrane region	101	7.7	1.20×10−8
	CC	GO:0044421	Extracellular region part	295	22.4	2.00×10−7
	CC	GO:0005576	Extracellular region	333	25.3	4.00×10−6

[i] DEGs, differentially expressed genes; BP, biological process; MF, molecular function; CC, cellular component.

KEGG pathway analysis

KEGG analysis revealed that the upregulated DEGs were enriched in focal adhesion, ovarian steroidogenesis, Rap1 signaling pathway, chemokine signaling pathway and PI3K-Akt signaling pathway, while the downregulated DEGs were enriched in the Rap1 signaling pathway, pathways in cancer, platelet activation, rheumatoid arthritis and cAMP signaling pathway (Table II).

Table II. KEGG pathway analysis of DEGs.

Table II.

KEGG pathway analysis of DEGs.

Expression	Pathway-ID	Name	Gene count	%	P-value	Genes
Upregulated	4510	Focal adhesion	12	3.5	9.20×10−4	RASGRF1, ROCK1, SHC4, COL4A3, COL4A4, COL6A5, EGFR, FLT1, HGF, IGF1R, SPP1, VAV2
	4913	Ovarian steroidogenesis	6	1.8	1.60×10−3	CYP1A1, CYP19A1, CGA, HSD3B1, IGF1R, INS
	4015	Rap1 signaling pathway	11	3.2	3.60×10−3	GNAO1, RASGRP2, EGFR, FGF22, FGFR4, FLT1, HGF, IGF1R, INS, LPAR4, PLCB2
	4062	Chemokine signaling pathway	9	2.6	1.60×10−2	CCL16, CCL25, GRK4, RASGRP2, ROCK1, SHC4, PLCB2, PPBP, VAV2
	4151	PI3K-Akt signaling pathway	13	3.8	1.70×10−2	COL4A3, COL4A4, COL6A5, CSF3R, EGFR, FGF22, FGFR4, FLT1, HGF, IGF1R, INS, LPAR4, SPP1
Downregulated	4015	Rap1 signaling pathway	31	2.4	2.50×10−6	GNAO1, KITLG, RAPGEF5, ADCY1, APBB1IP, CDH1, CNR1, F2R, CSF1, DRD2, EGFR, FGF14, FGF20, FGF7, FGFR2, IGF1R, IGF1, ITGB1, ITGB3, KDR, LCP2, MAGI1, MAP2K6, PIK3CD, PLCB1, PDGFRB, PFN4, PRKD2, RAC2, SIPA1L2, TLN1
	5200	Pathways in cancer	46	3.5	4.50×10−6	BAD, BCR, CXCL12, CXCL8, CXCR4, FAS, GNG11, KITLG, RUNX1T1, ARHGEF12, SPI1, WNT5B, ADCY1, AGTR1, CDH1, F2R, COL4A1, COL4A2, CDK6, EPAS1, EGFR, ERBB2, FGF14, FGF20, FGF7, FGFR2, FZD4, FZD5, IGF1R, IGF1, ITGB1, IL6, LAMA2, LAMB1, MMP2, MAPK10, NFKB2, PTCH1, PIK3CD, PLCB1, PDGFRB, PML, PTGER3, PIAS2, RAC2, TPM3
	4611	Platelet activation	19	1.4	3.90×10−4	FYN, ARHGEF12, ADCY1, APBB1IP, F2R, COL1A2, COL3A1, COL5A3, GUCY1A3, GUCY1B3, ITGB1, ITGB3, LCP2, PIK3CD, PLA2G4C, PLCB1, PTGS1, TLN1, TBXAS1
	5323	Rheumatoid arthritis	15	1.1	4.20×10−4	ATP6V0D2, ATP6V0E1, ATP6V1B1, CCL2, CXCL12, CXCL8, CD86, CTSK, CSF1, IL1B, IL11, IL6, HLA-DOA, HLA-DRB4, TNFSF13
	4024	cAMP signaling pathway	24	1.8	8.60×10−4	HTR1D, HTR1F, ABCC4, ATP1B4, ATP1B2, BAD, ADCY1, CAMK2A, CHRM2, F2R, DRD2, GIPR, GRIN3A, MC2R, MAPK10, PTCH1, PPARA, PIK3CD, PDE3A, PDE4D, PLN, PLD1, PTGER3, RAC2

[i] KEGG, Kyoto Encyclopedia of Genes and Genomes; DEGs, differentially expressed genes.

PPI network of the DEGs and core genes in the PPI network

Based on the information in the STRING database, the PPI network contained 1,123 nodes and 4,986 edges. The nodes indicate the DEGs, and the edges indicate the interactions between the DEGs. The top 10 high-degree hub nodes included albumin (ALB), epidermal growth factor receptor (EGFR), insulin (INS), interleukin 6 (IL-6), cadherin 1 (CDH1), FYN proto-oncogene (FYN), erb-b2 receptor tyrosine kinase 2 (ERBB2), interleukin 8 (IL8), C-X-C motif chemokine ligand 12 (CXCL12) and Ras-related C3 botulinum toxin substrate 2 (RAC2). Among these genes, ALB presented the highest node degree, which was 132. The core genes and their corresponding degree are shown in Table III. Then, the authors used MCODE to screen the modules of the PPI network (Fig. 2), and performed an enrichment analysis of the genes involved in the top three significant modules. The results demonstrated that the DEGs in modules 1–3 were principally related to neuroactive ligand-receptor interaction, the Rap1 signaling pathway, and protein digestion and absorption (Table IV).

Figure 2. Top three modules from the protein-protein interaction network. (A) Module 1, (B) module 2, (C) module 3.

Table III. The core genes and their corresponding degree.

Table III.

The core genes and their corresponding degree.

Gene	Degree
ALB	132
EGFR	113
INS	107
IL6	77
CDH1	76
FYN	72
ERBB2	71
IL8	67
CXCL12	62
RAC2	57
GNAO1	55
MMP2	54
TOP2A	54
EDN1	53
MYH14	52
ITGB1	50
ACTA2	50
KDR	49
IGF1	48
PPARA	48

Table IV. The enriched pathways of modules.

Table IV.

The enriched pathways of modules.

Modules	Enriched pathways	P-value	False discovery rate	Nodes
1	Neuroactive ligand-receptor interaction	9.60×10−10	2.68×10−14	ADORA3, CHRM2, CHRM4, CNR1, DRD2, GRM3, GRM6, HTR1D, HTR1F, PTGER3, S1PR1, S1PR5
	Chemokine signaling	6.10×10−6	2.03×10−7	CCL16, CCL25, CXCL12, CXCL2, CXCR4, IL8, PPBP
	Cytokine-cytokine	5.00×10−4	1.87×10−6	CCL16, CCL25, CXCL12, CXCL2, CXCR4, IL8, PPBP
2	Rap1 signaling pathway	2.40×10-16	5.46×10−21	CDH1, EGFR, F2R, FGF7, FLT1, GNAO1, HGF, IGF1, IGF1R, INS, KDR, LPAR4, PDGFRB, PLCB1, PLCB2
	PI3K-Akt signaling pathway	3.20×10−9	1.09×10−12	EGFR, F2R, FGF7, FLT1, HGF, IGF1, IGF1R, IL6, INS, KDR, LPAR4, PDGFRB
	Calcium signaling pathway	6.10×10−8	1.09×10−10	AGTR1, EGFR, ERBB2, F2R, GRPR, PDGFRB, PLCB1, PLCB2, PTGFR
3	Protein digestion and absorption	4.30×10−14	1.48×10−16	COL1A2, COL21A1, COL3A1, COL4A1, COL4A2, COL4A3, COL4A4, COL5A3
	ECM-receptor interaction	2.30×10−11	8.52×10−14	COL1A2, COL3A1, COL4A1, COL4A2, COL4A3, COL4A4, COL5A3
	Amoebiasis	7.80×10−11	2.44×10−13	COL1A2, COL3A1, COL4A1, COL4A2, COL4A3, COL4A4, COL5A3

Discussion

Osteosarcoma is the most common primary malignant bone tumour, is commonly observed in children and adolescents and shows a strong tendency for pulmonary metastasis (14). Pulmonary metastasis can cause medical therapy failure and high mortality rate in osteosarcoma patients (15,16). However, the present knowledge of the molecular mechanism of pulmonary metastasis in patients with OS remains insufficient. Therefore, the elucidation of molecular mechanisms to inhibit the metastasis of OS is imperative. High-throughput profiling technologies, such as microarrays, have been widely used and regarded as invaluable tools to identify potential therapeutic targets (17,18). The present study performed a comprehensive analysis and built a gene interaction network based on the gene expression profiles (GSE85537), comprising three non-metastatic OS samples and three metastatic OS samples. The results of the analysis demonstrated that there were 2,493 DEGs in the metastatic OS tumour samples compared with non-metastatic OS tumour samples, and among these DEGs, 485 were upregulated and 2,008 were downregulated. Moreover, GO and KEGG pathway analyses were performed to obtain a better understanding of the interactions of DEGs. Combining with the PPI network, key potential genes and pathways that may be associated with pulmonary metastasis of OS were identified.

The results of GO analyses indicated that the significant ontology categories included transmembrane transport, ion transport and aromatase activity. The system of transport is essential to every living cell and it had multiple functions, including allowing the entry to all essential molecules, exporting macromolecules such as proteins and DNA, providing cellular concentrations of ions (19,20). In a recent study on papillary thyroid carcinoma progression, downregulated DEGs were also predominantly enriched in transmembrane transport process (21). Additionally, ion transport system also stimulates the progression of pithelial-mesenchymal transition (EMT) in cervical carcinoma cells (22). Therefore, transmembrane transport and ion transport are fundamental processes and may serve roles in pulmonary metastasis of OS. Moreover, aromatase is an enzyme that converts testosterone into estradiol, and the control of aromatase activity mainly depends on Ca2+ transients (23). It has been reported that aromatase activity is associated with tumour-node-metastasis staging in human breast cancer (24). Cell motility is essential in the regulation of cancer progression, and aberrant cell motility occurs in malignant cancers and results in tumour metastasis (25). The activation of ERBB3-dependent signaling serves a key role in the regulation of cell motility, and regulation of cell motility contributes to intrahepatic metastasis and early recurrence of hepatocellular carcinoma (HCC) (26). Cell motility is also considered as an important process in the regulation of metastasis of oral cancer (27). Cell motility may therefore be important to pulmonary metastasis of OS. In addition, many studies have demonstrated that receptor binding is closely correlated with tumour metastasis, and these receptor bindings include fibroblast growth factor and vascular endothelial growth factor receptor binding, duffy antigen receptor for chemokines and growth factor receptor bound protein 2 Src homology 2 domain binding (28–30). Therefore, GO analyses can help identify the possible biological processes, molecular functions and cellular components involved in tumour metastasis.

Moreover, the KEGG pathways revealed that DEGs were enriched in focal adhesion, platelet activation, the Rap1 signaling pathway and the PI3K-Akt signaling pathway. Focal adhesion is a dynamic multi-protein complex that connects intracellular actin fibres and extracellular substrates (31,32). Focal adhesion kinase (FAK) is an effective therapeutic target of cell migration in various tumour cells (33–35). Furthermore, Rap1 is considered as a protein that can revert Ras transformation, and it is a key mediator in the control of integrin activation (36). A recent study revealed that the cAMP/Epac/Rap1 signaling cascade had a crucial role in blood-tumour barrier hyperpermeability (37). Cancer cell-induced platelet activation is identified as a critical event responsible for prometastatic activity of platelets, and blocking platelet aggregation can inhibit the progression of skeletal metastases, yet the mechanism underlying this effect is still unknown (38). Moreover, the PI3 K/Akt signaling pathway is involved in various cellular processes, such as inflammation, autophagy and cancer progression (39), and many studies have demonstrated that activation of the PI3K-Akt signaling pathway can increase the metastatic potential (40–42). Therefore, further understanding of these signaling pathways can help us to elucidate the crucial mechanism of tumour metastasis.

Furthermore, the authors analyzed the PPI network and found that ALB, EGFR, INS, IL6, CDH1, FYN, ERBB2, IL8, CXCL12 and RAC2 were the top 10 core genes, which may be potential therapeutic targets for pulmonary metastasis of OS. Node degree refers to the number of interacting partners per protein. Among these genes, ALB showed the highest node degree. Albumin is one of the major plasma proteins and is synthesized and secreted primarily by hepatocytes, and it is often used for diagnosis of HCC (43). Albumin is a major indicator of a favourable outcome of HCC, and it can suppress cell proliferation via reduction of phosphorylation of Rb (retinoblastoma) proteins and increase of p21 and p57 (44). A recent study indicated that the expression of albumin increased in human hepatoma BEL-7402 cells after bufalin was used, and ALB may be associated with bufalin-mediated inhibition of the invasion and metastasis of HCC cells (45). Albumin mRNA was also found in lung by means of reverse-transcription polymerase chain reaction, so albumin may serve a role in pulmonary metastasis (46). The present study indicated that ALB was a key gene with the highest node degree in OS, therefore ALB may be a potential target for inhibiting OS tumour metastasis, but further experimental verification is necessary.

The EGFR is a key component in the mitogen-activated protein kinase (MAPK) pathway, and anti-EGFR therapy can lead to inhibition of tumour growth, invasion and metastasis via inhibition of the MAPK pathway (47). Moreover, nasopharyngeal carcinoma (NPC), which has the highest rate of metastasis, has been extensively studied, and inhibition of EGFR can inhibit NPC cell migration and invasion (48). Additionally, EGFR has been shown to promote survival of prostate tumour-initiating cells and circulating tumour cells that metastasize to bone (49). These results corresponded with the function analysis in the present study, and EGFR is likely to be an important regulated target in pulmonary metastasis of OS. Therefore, EGFR may be a potential therapeutic target and prognostic indicator in OS patients, and a thorough study of EGFR in OS is needed. Furthermore, a recent study demonstrated that inflammatory cytokines, such as IL-6, promoted proliferation, migration, invasion and EMT of gallbladder cancer both in vitro and in vivo (50). Overexpression of human epidermal growth factor receptor 2 (ErbB2 or HER2) occurs in ~30% of breast cancer patients, and overexpression of activated ErbB2 in the mammary epithelium can lead to the rapid induction of metastatic multifocal mammary tumours (51,52). Therefore, these core genes in the current study have been suggested to be associated with tumour metastasis in the previous studies, so the analysis of these core genes is useful for understanding the molecular mechanisms and identifying therapeutic targets of OS with pulmonary metastasis. However, subsequent prospective studies will be required to further confirm the function of these core genes in pulmonary metastasis of OS.

Module analysis of the PPI network demonstrated that the pulmonary metastasis of OS was associated with neuroactive ligand-receptor interaction, the Rap1 signaling pathway, protein digestion and absorption and other processes. It has been demonstrated that neuroactive ligand-receptor interaction is important in HCC progression and it is present in the early-, middle- and late-stages of HCC (53,54). Moreover, the neuroactive ligand-receptor interaction pathway was enriched in prostate tumours from African-American patients (55). Therefore, neuroactive ligand-receptor interaction appears to be important in pulmonary metastasis of OS. Moreover, the protein digestion and absorption pathway has been reported to be associated with pancreatic neuroendocrine tumours and breast cancer (56,57). Thus, the roles of these pathways in pulmonary metastasis of OS need to be confirmed by the future studies, and these results may lead to additional therapeutic alternatives in the patients with OS.

In conclusion, the present study identified 2,493 DEGs, which may be involved in the progress of pulmonary metastasis in OS patients, via a comprehensive bioinformatics analysis. GO term, KEGG pathway and PPI network analyses provided a set of related genes and pathways to help elucidate the molecular mechanisms of pulmonary metastasis. Further experimental studies are needed to confirm these results and should help determine potential targets for inhibiting pulmonary metastasis in OS patients.

Acknowledgements

The present study was supported by the State Key Program of National Natural Science Foundation of China (grant no. 81330042), International Cooperation Program of National Natural Science Foundation of China (grant no. 81620108018), Special Program for Sino-Russian Joint Research Sponsored by the Ministry of Science and Technology, China (grant no. 2014DFR31210) and Key Program Sponsored by the Tianjin Science and Technology Committee, China (grant nos. 13RCGFSY19000 and 14ZCZDSY00044).

References