This article is mentioned in:

Introduction

Osteosarcoma (OS) is one of the most common and destructive primary bone malignancies, with an annual incidence of ~5 cases per 1,000,000 individuals worldwide (1,2). The incidence of OS exhibits bimodal distribution, with the primary peak occurring in adolescence, while the incidence of OS in the elderly is relatively low (3,4). OS originates from mesenchymal tissues and is characterized by malignant osteogenesis and osteoblast differentiation, which confers a high degree of malignancy and strong metastatic ability. The incidence of common pulmonary metastasis through the blood is currently >85%, and due to therapeutic limitations in previous years, amputation was the principal method of treatment for OS. However, the survival rate of patients treated by surgery alone is 15–17%, and the appearance and function of the limbs after amputation can negatively impact patient welfare (5). However, with rapid developments in medicine and technology, limb-salvage surgery following OS has become possible. To inhibit the growth and metastasis of tumors, OS treatment has progressed from surgery alone to a combination of surgery and adjuvant chemotherapy (6). However, the overall survival rate of patients with OS remains unsatisfactory, with ~50% developing fatal lung metastases in the advanced stages of disease (7). Due to the limitations of current diagnostic techniques, only 15–20% of patients with metastases can be diagnosed by auxiliary examination, thus metastasis remains the leading mortality-associated factor for patients with OS (8). Therefore, to identify novel biomarkers and therapeutic targets, it is of great significance to study the mechanisms of OS occurrence and metastasis at the molecular level.

Bioinformatics tools are widely used to explore genes involved in the development of OS. Using differential gene expression analysis, Xie et al (9) identified ZNRD1, GPR68, CAT, FUT3, ANPEP and CDK1 as key genes associated with chemotherapy resistance in OS. Analyzing the DNA methylation profiles of OS, Chen et al (10) revealed that the abnormal methylation of NMU and NMUR1 may contribute to the development of OS. Furthermore, based on integrated bioinformatics analysis, Wang et al (11) showed that microRNA-203 may be a suppressor of OS. Through WGCNA, Tian et al (12) showed that IGFBP5, IGFBP6, WISP3, and MYL2 which involved in insulin-like growth factor binding may play key roles in the metastasis of osteosarcoma. Similarly, through WGCNA analysis, Wang et al (13) indicated that MEPE, BPIFB1, HBA2, and SERPINB3 were key genes involved in the metastasis of osteosarcoma. However, the majority of these findings were only obtained using prediction tools, and lacked experimental validation.

Mitogen-activated protein kinases (MAPKs) are a class of serine/threonine protein kinases expressed in almost all cell types, which regulate evolutionarily conserved signal transduction pathways and various cellular functions, including proliferation, migration, apoptosis and differentiation (14). Abnormal MAPK signaling plays a key role in the occurrence and development of different types of cancer (15). MAPK15 is the most recently identified member of the MAPK family (16), and it is widely accepted to be upregulated in a variety of cancer types. Studies have shown that unlike other members of the MAPK family, the phosphorylation of MAPK15 is mainly self-phosphorylation, so its total protein level is positively correlated with the phosphorylation level (17,18). Additionally, the overexpression of MAPK15 promoted gastric cancer cell proliferation by stabilizing the expression of c-Jun (19). MAPK15 is also regarded as an oncogene in male germ cell tumors (20), and was revealed as a core gene involved in the radio-resistance of nasopharyngeal carcinoma cells (21). However, there are few studies focused on the relationship between phosphorylation location and protein activity of MAPK15. Similarly, the role of MAPK15 in OS remains unknown.

In present study, a combination of bioinformatics and related experimental methods were used to identify metastasis-associated genes in OS. The present study revealed that MAPK15 may be a novel biomarker for the diagnosis of OS, as well as an effective target for clinical treatment.

Materials and methods

Data source

The analysis and experimental procedures of the present study are outlined in Fig. 1. To identify key genes associated with the metastasis of OS, human gene expression data from OS tissues (dataset GSE87624) were downloaded from the Gene Expression Omnibus (http://www.ncbi.nlm.nih.gov/geo/), and subsequently used to perform WGCNA and DEG analysis. The GSE87624 dataset, comprising 44 human patient OS samples and 3 normal bone samples, was first published by Scott et al (22).

Figure 1. Flow diagram of current study. The current study was divided into two parts: Bioinformatics methods for exploring hub genes for the metastasis of OS (yellow) and biological experimentation used to detect the effects of hub genes on the metastasis of OS cells (green). In the part of bioinformatics methods, WGCNA was performed to select core genes in the gene module associated with the metastasis of OS and DEG analysis for exploring DEGs between metastatic or non-metastatic OS tissues. In the part of biological experiments, wound healing assays, transwell assays and a lung metastasis model were utilized. Additionally, the molecular mechanism was explored using western blotting and a rescue experiment. WGCNA, weight gene co-expressed network analysis; DEG, differentially expressed gene; OS, osteosarcoma; GEO, Gene Expression Omnibus.

Identification of DEGs

EdgeR is a Bioconductor software package (version number: 3.9; http://www.bioconductor.org/packages/release/bioc/html/edgeR.html) with the capacity to determine the differential expression of genes (23). Therefore, in the present study, the level of differential gene expression between OS and normal bone tissues was assessed using the EdgeR package (version number: 3.9; http://www.bioconductor.org/packages/release/bioc/html/edgeR.html) in R. The DEGs between metastatic and non-metastatic OS tissues were also analyzed. An over-dispersed Poisson model was used to explain biological and technological variability, and the Empirical Bayesian method was used to mitigate the excessive dispersion of transcription and improve the reliability of reasoning. A combination of adjusted P<0.05 and |log2 fold change (FC)|≥1 was used as the stringent cutoff to determine the significance of the differences in gene expression. Then, WGCNA was conducted to determine the DEGs between normal bone and OS tissues.

WGCNA

In the present study, the R package WGCNA (v.1.66; http://cran.r-project.org/web/packages/WGCNA/index.html) was employed to conduct WGCNA for all DEGs between normal bone and OS tissues. First, the DEG expression profiles and their associated clinical information were imported. Pearson's correlation analysis was then performed to cluster the samples and detect outliers; the outlier threshold was set as 10,000. Then, all genes-pairs were analyzed using Pearson's correlation analysis and a matrix of similarity was constructed. In addition, to identify specific modules, WGCNA uses a soft-thresholding procedure to avoid the selection of an arbitrary cut-off. The β value represented a soft-thresholding parameter that could emphasize strong correlations between genes and penalize weak correlations to ensure a scale-free network. The network connectivity (k) of the gene was defined as the sum of its adjacency with all the other genes for network generation. The decision value of the threshold power was determined on the basis of the scale-free topology criterion, which aims to mimic a network structure commonly found in nature. In the present study, β=3 was selected based on the scale-free topology criterion >0.85. The cutreeDynamic function was used for adaptive branch pruning of hierarchical clustering dendrograms and the dynamicTreeCut package (v.1.66; http://cran.r-project.org/web/packages/WGCNA/index.html) was adopted to generate co-expression modules. Subsequently, to further analyze the gene modules, the dissimilarity of the module eigengenes (ME) was calculated using the moduleEigengenes function in the R WGCNA package, which was defined as the first principal component of a given module and considered to be representative of the gene expression profiles in a module. A cut-off line for the module dendrogram was selected and the module was merged. Eventually, the adjacency was converted into a topological overlap matrix (TOM), and modules were subjected to hierarchical cluster analysis according to the TOM-based dissimilarity measure with the mini-size set as 20.

Screening of metastasis-related hub genes

The correlation between the co-expression modules and the metastasis of OS was analyzed. Firstly, gene significance (GS) was defined as the log10 transformation of the corresponding P-value (GS=lgP) of the correlation between gene expression and pathological stage. Secondly, module eigengenes (MEs) were determined by principal component analysis and selected to represent the mean measure of the overall co-expression module. Finally, the degree of correlation between MEs and clinical traits was determined using Pearson's correlation analysis. Significant modules were identified according to the absolute value of moderate intensity correlation (>0.3), with a significance threshold of P<0.05. All genes in the significant modules related to the significance of metastatic traits were screened to identify the core genes, according to GS>0.2 and MM>0.8. To obtain credible metastasis-related hub genes, Bioinformatics and Evolutionary Genomics (http://bioinformatics.psb.ugent.be/webtools/Venn/) was used to identify genes which were not only hub genes in the metastasis-associated gene modules in WGCNA, but also differentially expressed genes in differential expression analysis. The metastasis-related hub genes were then employed for experimental evaluation.

Tissue ethics

All OS tissues used in the current study were obtained from the Zhujiang Hospital affiliated South Medical University (Guangzhou, China) between July 2014 and September 2019. Patients were enrolled using the advanced procedures if they met the following criteria: i) The tissues were obtained from operation and two pathologists diagnosed as osteosarcoma; ii) patients diagnosed and treated for the first time; and iii) patients willing to participate. The exclusion criteria were as follows: i) Patients complicated with other malignancies; ii) patients with other systemic disease; iii) patients received treatment prior to admission; and iv) patients and/or their families refused to participate. A total of 26 pairs of OS and adjacent non-tumor tissues were obtained from extremities of OS patients, of which 11 pairs were obtained from patients without, and 15 pairs were provided from patients with OS metastasis at diagnosis. The procedures of the current study were approved by the Clinical Ethics Management Committee of Zhujiang Hospital, South Medical University, and all patients provided informed consent in writing.

Reverse transcription-quantitative (RT-q) PCR

Total RNA was extracted from OS tissues using TRIzol® reagent (Takara Bio, Inc.), and an ultraviolet spectrophotometer (Bio-Rad Laboratories, Inc.) was used to assess the concentration and quality; samples with OD260/OD280 of 1.8–2.1 were considered to qualify for subsequent experimentation. cDNA synthesis was conducted using the 1st Strand cDNA Synthesis kit (Shanghai Yeasen Biotech Co., Ltd.,) according to the manufacturer's instructions, and HieffTM qPCR SYBR® Green Master Mix (Shanghai Yeasen Biotech Co., Ltd.,) was used to performed RTq-PCR. The thermocycling conditions were as follows: 95°C for 5 min, 40 cycles of 95°C for 30 sec, annealing at 60°C for 30 sec, and a final elongation step at 72°C for 30 sec. The 2−ΔΔCt method (24) was employed to measure the relative expression levels of target genes, and GAPDH was used as the loading control. The primers used in the present study were as follows: Deleted in lung and esophageal cancer protein 1 (DLEC1) forward, 5′-CCAAAACGCGGAGGTCTTTAG-3′ and reverse, 5′-GGGAGGAATACAAGGAGGACT-3′; forkhead box J1 (FOXJ1) forward, 5′-GCCTCCCTACTCGTATGCCA-3′ and reverse, 5′-GCCGACAGGGTGATCTTGG-3′; MAPK15 forward, 5′-GGGCCTATGGCATTGTGTG-3′ and reverse, 5′-TCTCTGGGCATCTGTCTTATCC-3′; and GAPDH forward, 5′-GGAGCGAGATCCCTCCAAAAT-3′ and reverse, 5′-GGCTGTTGTCATACTTCTCATGG-3′.

Immunohistochemical staining (IHC)

The OS tissues were dehydrated and embedded in paraffin (Wuhan Servicebio Technology Co., Ltd.); the tissues were subsequently cut into 4-µm slices for use in the present study, where they were deparaffinized using xylene and rehydrated in a descending alcohol series under room temperature. After restoration with sodium citrate the samples were treated with 3% H2O2 to block endogenous peroxidase activity, and then blocked using 5% bovine serum albumin (BSA; Servicebio, Wuhan, China) for 30 min under room temperature. The specimens were subsequently incubated with a primary anti-MAPK15 antibody (dilution: 1:400; cat. no. ab137619; Abcam, USA) for 12 h at 4°C, followed by a second incubation with horseradish peroxidase-conjugated secondary antibodies (dilution: 1:200; cat. no. G1210-2-A-100; Wuhan Servicebio Technology Co., Ltd.) for 2 h at room temperature. After subsequent development using the Cell and Tissue Staining HRP-DAB kit (Beyond) according to the manufacturer's protocol, images were captured with an orthophoto microscope (magnification, ×400).

Cell culture and transfection

The human OS MG63 and 143B cell lines were purchased from the American Type Culture Collection. Cells were cultured in Dulbecco's modified Eagle's medium (DMEM; Invitrogen; Thermo Fisher Scientific, Inc.) supplemented with 100 U/ml penicillin and 100 µg/ml streptomycin (Wuhan Boster Biological Technology, Ltd.,), and 10% fetal bovine serum (FBS; Gibco; Thermo Fisher Scientific, Inc.,) at 37°C with 5% CO2. Human MAPK15-targeting short hairpin RNA (shRNA) and the corresponding scramble shRNA oligonucleotide sequences were cloned into the pSuper-retro-puro vector to generate pSuper-retro-MAPK15-RNAi(s) and pSuper-retro-scramble-RNAi(s). The shRNA sequences were as follows: scramble shRNA, 5′-CCGGAATTCTCCGAACGTGTCACGTCTCGAGACGTGACACGTTCGGAGAATTTTTTTG-3′; and MAPK15 shRNA, 5′-CCGGGCTTGGAGGCTACTCCCCTCGAGGGGAGTAGCCTCCAAGCTTTTTG-3′. In order to obtain cell lines that stably exhibited low-level MAPK15 expression (MAPK15-knockdown), the cells were treated with puromycin (0.5 µg/ml) for 12 days following transfection. MAPK15 (NM_139021.3) was used as backbone and subcloned into pcDNA3.1 vector to construct MAPK15 overexpression plasmid (Guangzhou GeneCopoeia Co., Ltd.). The transfection processes of MAPK15 overexpression plasmid or pcDNA3.1 empty vector were preformed using Lipo2000 (Invitrogen; Thermo Fisher Scientific, Inc), and cells transfected pcDNA3.1 empty vector were seted as negative control (NC).

Wound-healing assay

After transfection, total 1×106 MG63 and 143B cells were plated into per well of 6-well plates and incubated at 37°C for 24 h. At 100% confluence, the cells were serum starved for 24 h. Then, the cell monolayer was then scraped with a 200-µl sterile pipette tip to form a central linear wound. After washed by PBS for three times, cells were treated with serum free DMEM medium with or without a c-Jun inhibitor SP600125 (MCE) and cultured at 37°C with 5% CO2. Then, the cells were photographed under an optical microscope (magnification, ×40) and the wound closure rate was recorded after 24 h. At 0 h, the wound was regarded as 100% of the average clearance. The relative migration rate was calculated as follows: Relative migration rate=(Streat-0 h -Streat-24 h)/(SNC-0 h-SNC-24 h) ×100%, where Streat-0 h and Streat-24 h was the area of the scratch at 0 h and 24 h in the treatment group, SNC-0 h and SNC-24 h was the area of the scratch at 0 h and 24 h in the NC group.

Transwell invasion assay

A total of 2×104 MG63 and 143B cells were resuspended in serum-free DMEM and placed into the upper chamber of a Transwell insert (Invitrogen; Thermo Fisher Scientific, Inc.) with 8-µm pores, which had been pre-coated with Matrigel at 37°C (Becton, Dickinson and Company); DMEM containing 10% FBS was placed in the lower chamber as a chemical attractant, and the inserts were incubated at 37°C (5% CO2) for 24 h. After removing the non-invasive cells on the top of the membrane, the cells were fixed with 4% paraformaldehyde for 15 min, stained with 0.5% crystal violet for a further 30 min, counted and photographed under an optical microscope equipped with a digital camera (magnification, ×40). The number of invasive cells was counted in five random fields per sample.

In vivo lung metastasis model

A total of 10 4-week-old female nude mice (Beijing HuaFukang Biological Technology Co. Ltd) were housed in a facility at 23–24°C, and the light-dark cycle was set at 12-h intervals. After one week of adaptive feeding, mice were randomly divided into the sh-scramble (n=5) and sh-MAPK groups (n=5). A total of 1×107 MAPK15-knockdown 143B cells and the corresponding control 143B cells were injected into the tail veins of the mice in the corresponding groups. Animal health and behavior after injecting were monitored each day. While the mouse had the features of hard breathe and limitation of motion, the mouse was sacrificed in order to reduce animal suffering. When half of mice in any group (n≥3) was sacrificed due to the above reasons, this experiment should be ended. Cervical dislocation was performed manually in all mice and resulted in euthanasia within ~10 sec in 5 weeks post-injection, and the death of mice was verified by the absence of a heart beat and the onset of rigor mortis. Ensuring the death of mice, the lungs were removed and the tissues were fixed in 4% paraformaldehyde for 30 min under room temperature, and subsequently stained with hematoxylin and eosin prior for 1 min under room temperature to image capture. The animal care and experimental procedures used in the present study were approved by the Institutional Animal Care and Use Committee of the Southern medical university. In addition, all mouse procedures, euthanasia and surgery, including injections of 143B cells, were conducted painlessly or under anesthesia using a combination of hydrochloric acid medetonidine 0.3 mg/kg + midazolam 4 mg/kg + butorphanol tartrate 5 mg/kg through tail vein injection (25). In the present study, the success of anesthesia is evaluated as the decrease of respiratory rate, the increase of respiratory depth and the disappearance of eyelid and corneal reflexes.

Western blot analysis

The samples were homogenized in RIPA (including 1% PMSF) to extract the total protein, according to the manufacturer's instructions. Equal amounts of protein (25 µg) were loaded into each lane of a 10% gel, separated by SDS-PAGE and transferred onto PVDF membranes (Merck KGaA). The membranes were then blocked with 5% non-fat milk (Wuhan Servicebio Technology Co., Ltd.,) at room temperature for 2 h, and incubated overnight at 4°C with primary antibodies (dilution: 1:1000) against GAPDH (cat. no. 60004-1-Ig; ProteinTech Group, Inc.), MAPK15 (cat. no. 13452-2-AP; ProteinTech Group, Inc.), c-Jun (cat. no. 66313-1-Ig; ProteinTech Group, Inc.), p-c-Jun (cat. no. 5464; Cell Signaling Technology, Inc.), MMP-9 (cat. no. 10375-2-AP; ProteinTech Group, Inc.,) and MMP-2 (cat. no. 10373-2-AP; ProteinTech Group, Inc.). The membranes were washed three times with TBST and incubated with a horseradish peroxidase-conjugated goat anti-mouse or goat anti-rabbit IgG secondary antibody (dilution: 1:3000; Wuhan Servicebio Technology Co., Ltd.). The protein bands were visualized with an enhanced chemiluminescent (ECL) kit (Wuhan Servicebio Technology Co., Ltd.) using Bio-rad Chemidoc image software (v.5.2.1; Bio-Rad Laboratories, Inc.). In the present study, GAPDH was used as a reference gene for normalizing the expression of MAPK15, MMP2 and MMP9, while relative expression of p-C-Jun was calculated according to formula (p-c-Jun/GAPDH)/(total c-Jun/GAPDH).

Statistical analysis

All statistical analyses were carried out using GraphPad Prism (v.5.0; GraphPad Software, Inc.). One-way analysis of variance and the LSD-t test were performed to analyze the statistical differences between≥3 groups, while the independent sample t-test was used to detect differences between two groups. P<0.05 was considered to indicate a statistically significant difference.

Results

Co-expression module construction by WGCNA

Prior to the construction of a co-expression module, DEG analysis was performed for normal bone and OS tissue data retrieved from the GSE87624 dataset, and 1,043 DEGs were identified (Table SI). After removing outlier samples using cluster analysis, the DEGs and corresponding clinical information of 32 OS tissue were used to construct a co-expression module by WGCNA (Fig. S1). When the soft power β was set as 3, the scale independence of the topology network reached higher than 0.85 (Fig. 2A) with the mean connectivity lower than 50 (Fig. 2B). So a soft power of β=3 was selected as the soft threshold for performing subsequent analyses. The result showed that when β=3, the topological overlap matrix was able to meet the scale-free topology criterion with R2=0.9 (Fig. 2C and D). The dynamic tree cut function was used to prune the branches in hierarchical clustering dendrograms that determined the generation of the co-expression modules. Then, the MEs were calculated by the moduleEigengenes function to quantify the co-expression similarity of the modules and the clustered modules were merged based on the similarity. A total of 16 co-expressed modules were identified, while unclustered genes were distributed in the grey module (Fig. 3). The 16 co-expression modules were then used for subsequent analysis.

Figure 2. Analysis of network topology for different soft-thresholding powers. (A) Scale independence and (B) mean connectivity of various soft-thresholding values (β). (C) Histogram of connectivity distribution when β=5. (D) Checking the scale-free topology when β=5.

Figure 3. Clustering dendrograms of all DEGs between normal bone tissues and osteosarcoma tissues, with dissimilarity based on topological overlap, together with assigned module colors. Dynamic tree cut algorithm was applied to the dendrogram for module identification, and the closed modules (difference of feature vectors <0.2) were merged into new modules. Different colors represent different gene modules and there are 16 co-expressed modules in the weighted gene co-expression network analysis network, including tan, black, magenta, midnight blue, salmon, pink, green, turquoise, purple, light cyan, brown, cyan, red, green yellow, blue and yellow. DEG, differentially expressed gene.

Identification of metastasis-associated modules and screening of metastasis-related hub genes

Using WGCNA, the relationship between modules and tumor metastasis was investigated, and the MEtan (cor=0.45; P=0.0094) and MEmidnightblue modules (cor=0.41; P=0.018) (Fig. 4A and B), containing 33 and 24 genes, respectively, were found to be associated with OS metastasis. All significantly associated genes were screened according to GS>0.2 and MM>0.8 and ultimately, 16 module core genes were obtained from the metastasis-related modules, including PRR15, PLEKHH1, MIR200A, OVOL1, PCSK6, MIR200B, RHOV, TRPM1, MAPK15, DLEC1, GAGE2A, FOXJ1, FOXA3, PAGE1, C7orf57 and GZMB (Fig. 4C and D). DEG analysis was also performed between metastatic and non-metastatic OS tissue data and 628 DEGs were identified (Table SII). The intersection of the core genes in the metastasis related-modules, and DEGs between metastatic and non-metastatic OS tissues was then determined. The expression change fold of 16 module core genes obtained from the metastasis-related modules between metastatic and non-metastatic OS tissues were also showed in Table I. According to the significance difference of FDR<0.05, the results indicated DLEC1, FOXJ1 and MAPK15 as core genes in the metastasis related-modules, which were also highly expressed in metastatic OS tissues (Fig. 5A).

Figure 4. Identification of significance modules associated with metastasis and screening module core genes. (A) Association of module eigengenes with metastatic traits. Each cell contains the corresponding correlation and P-value (indicated in brackets). The table is color-coded by correlation according to the color legend. (B) Distribution of metastasis-related genes in all modules. Modules are presented on the x-axis, and the enrichment significance is shown on the y-axis. All the modules correspond to the modules in (A). Scatterplots of gene significance for metastasis vs. MM in (C) tan and (D) midnight blue modules. The corresponding correlation coefficient and P-value are listed above the scatterplot. The red lines indicate GS=0.2 and MM=0.8, respectively. MM, module membership.

Figure 5. MAPK15 was most significantly upregulated in metastatic OS tissues. (A) Venn diagram showing intersected genes of core genes in the gene module associated with the metastasis of OS and DEGs between metastatic and non-metastasis OS tissues. (B) RT-qPCR was used to detect the expression of FOXJ1 in metastatic (n=15) and non-metastasis (n=11) OS tissues. *P<0.05. (C) RT-qPCR was used to detect the expression of DLEC1 in metastatic (n=15) and non-metastasis (n=11) OS tissues. (D) qRT-PCR was used to detect the expression of MAPK15 in metastatic (n=15) and non-metastatic (n=11) OS tissues. **P<0.01. (E) Immunohistochemistry was performed to detect the expression of MAPK15 in metastatic (n=15) and non-metastasis (n=11) OS tissues. Scale bar, 100 µm. MAPK15, mitogen-activated protein kinase 15; DEG, differentially expressed gene; OS, osteosarcoma; RT-qPCR, reverse transcription-quantitative PCR.

Table I. Expression change of 16 module core genes obtained from the metastasis-related modules between metastatic and non-metastatic osteosarcoma tissues.

Table I.

Expression change of 16 module core genes obtained from the metastasis-related modules between metastatic and non-metastatic osteosarcoma tissues.

Gene	logFC	FDR
FOXJ1	4.35	<0.01
MAPK15	2.44	<0.01
DLEC1	1.78	0.01
RHOV	1.36	0.07
GZMB	1.74	0.10
OVOL1	1.69	0.12
GAGE2A	3.43	0.17
MIR200B	2.48	0.21
PRR15	1.52	0.21
FOXA3	1.65	0.23
PAGE1	2.21	0.26
TRPM1	1.45	0.28
MIR200A	1.98	0.30
C7orf57	1.08	0.32
PCSK6	1.07	0.34
PLEKHH1	0.82	0.38

[i] LogFC, Log (Fold change); FDR, false discovery rate.

Validating the expression of DLEC1, FOXJ1 and MAPK15 in OS tissues

To verify that the DLEC1, FOXJ1 and MAPK15 hub genes were associated with OS metastasis, RT-qPCR was performed to quantify the expression of these genes in OS tissues from 26 patients. The results revealed that the mRNA expression levels of MAPK15 and FOXJ1, but not DLEC1, were significantly increased in the OS tissues of patients who exhibited tumor metastasis at diagnosis (Fig. 5B-D); the increased expression level was most significant for MAPK15. Therefore, IHC was also performed and the results showed that the MAPK15 protein was also highly expressed in OS tissues from patients with metastasis at diagnosis (Fig. 5E).

MAPK15-knockdown suppresses OS cell migration and invasion in vitro

Targeting MAPK15 shRNAs were used to construct MG63 and 143B cell lines stably expression low levels of MAPK15 (MAPK15-knockdown). The results of the wound-healing assays showed that the migration rate of the sh-MAPK15 group was significantly decreased compared with that of the sh-scramble group (Fig. 6A and B). Furthermore, the results of the Transwell assays demonstrated that inhibiting the expression of MAPK15 decreased the number of invasive MG63 and 143B cells per field (Fig. 6C and D). Collectively, these results suggest that MAPK15-knockdown reduces the migratory and invasive abilities of OS cells.

Figure 6. Inhibition of MAPK15 decreases MG63 and 143B cell migration and invasion in vitro. (A) Wound-healing assay was used to detect the migration ability of the sh-scramble and sh-MAPK15 groups in 143B and MG63 cells. Sale bar, 100 µm. (B) Statistical chart showing the migration ability of the sh-scramble and sh-MAPK15 groups in 143B and MG63 cells. (C) Transwell assays were used to detect the invasion ability of both groups in MG63 and 143B cells. Scale bar, 100 µm. (D) Statistical chart showing the number of invasive cells per field. **P<0.01. MAPK15, mitogen-activated protein kinase 15; sh-, short hairpin (RNA).

MAPK15 inhibition decreases the metastasis of 143B-cell tumors in vivo

Studies have shown that about 15–20% of patients have evidence of metastases at diagnosis, mostly in the lungs (26,27). Therefore, A mouse model of lung metastasis was used to detect the effect of MAPK15 inhibition on the metastasis of OS cells in vivo. 143B cells transfected with sh-MAPK15 and sh-scramble lentivirus were injected into the tail vein of 4-week-old female nude mice. The time of mice in sh-scramble and sh-MAPK15 group free from pulmonary embolism induced breathing and movement difficultly was recorded (Fig. 7A). The results of lung tissues diagram (Fig. 7B and C) and HE stain (Fig. 7D) both indicated that metastasis was significantly decreased in the lungs of mice inoculated with MAPK15-knockdown 143B cells, suggesting that MAPK15 promotes the metastatic potential of OS.

Figure 7. Inhibition of MAPK15 significantly decreases the metastasis of 143B cells tumors in vivo. (A) Rate of mice free from pulmonary embolism mediated breathing and movement disorders before experiment end. (B) Typical examples of the lungs in the sh-scramble and sh-MAPK15 group mice. (C) Statistical chart showing the metastatic lungs of mice in the sh-scramble and sh-MAPK15 groups. (D) Typical images of H&E staining exhibiting the metastasis of 143B cells in lungs of the sh-scramble and sh-MAPK15 group mice. Scale bar, 1,000 µm. **P<0.01. MAPK15, mitogen-activated protein kinase 15; sh-, short hairpin (RNA).

MAPK15 significantly regulated the c-Jun/MMPs pathway in OS cells

Previous studies have indicated that MAPK15 increases the phosphorylation level of c-Jun and activated Jun, and therefore the levels of c-Jun and p-c-Jun were detected after MAPK15 inhibition or overexpression. The results showed that following MAPK15 inhibition, the expression of p-c-Jun was significantly decreased in both MG63 and 143B cells. Then, the expression of MMP2 and MMP9, downstream proteins of c-Jun, were also detected, revealing that both MMP2 and MMP9 expression was significantly decreased following MAPK15-knockdown (Fig. 8A and B). Similarly, we found that p-C-Jun was significantly increased while MAPK15 overexpression, as well as increasing the expression of MMP2 and MMP-9 (Fig. 8C and D).

Figure 8. MAPK15 significantly regulates C-Jun/MMPs signaling pathway. (A) Western blotting was used to detect the expression levels of MAPK15, c-Jun, P-c-Jun, MMP2 and MMP9. MAPK15 expression was inhibited. (B) Statistical analysis for MAPK15, c-Jun, p-c-Jun, MMP2 and MMP9. (C) Western blotting was used to detect the expression of MAPK15, c-Jun, P-c-Jun, MMP2 and MMP9. (D) Statistical analysis for MAPK15, c-Jun, p-c-Jun, MMP2 and MMP9. **P<0.01. MAPK15, mitogen-activated protein kinase 15; MMP, matrix metalloproteinase; NC, normal control; p-, phosphorylated-.

Inhibition of c-Jun reverses the effects of MAPK15-overexpression on the migration and invasion of OS cells

To determine whether the c-Jun/MMP pathways are involved in the migration and invasion induced by MAPK15, a c-Jun inhibitor SP600125 was used. Wound-healing assays showed that the overexpression of MAPK15 increased the migratory ability of OS cells, while SP600125 blocked this effect (Fig. 9A and B). Similarly, Transwell assays showed that the number of invasive cells was significantly increased by MAPK15 overexpression, which was subsequently blocked by SP600125 (Fig. 9C and D).

Figure 9. Inhibition of c-Jun blocks the promoting effects of MAPK15 overexpression on the migration and invasion of osteosarcoma cells. Osteosarcoma cells were divided into three groups: NC group, MAPK15 overexpression (MAPK15-U) group and MAPK15 overexpression plus SP600125 (40 µM) treatment group (MAPK-U+SP600125). (A) Wound-healing assay to detect the migration ability of each group in 143B and MG63 cells. Scale bar, 100 µm. (B) Statistical chart showing the migration ability of each group in 143B and MG63 cells. (C) Transwell assay to detect the invasion ability of each group in MG63 and 143B cells. Scale bar, 100 µm. (D) Statistical chart showing the number of invasion cells per field of each group in 143B and MG63 cells. **P<0.01. MAPK15, mitogen-activated protein kinase 15; NC, negative control.

Discussion

OS is the most common primary malignant tumor in the skeletal system with strong characteristics of invasion, metastasis and recurrence. Current clinical treatments are not effective for patients with OS-associated metastasis and recurrence, thus the prognosis of patients with end-stage OS with metastasis has not improved (28). In addition, epidemiological studies have shown that metastasis is the leading cause of death in patients with OS (29). Therefore, exploring the molecular mechanism of OS and discovering novel therapeutic targets has important clinical significance for improving the survival and prognosis of patients with OS.

In the present study, two bioinformatics methods (WGCNA and DEG analysis) were used to analyze data from the GSE87624 dataset, in order to identify key oncogenes associated with the metastasis of OS. According to WGCNA, 16 genes including DLEC1, FOXJ1 and MAPK15 was predicted as core genes in metastasis-related modules. Furthermore, through DEG analysis for non-metastatic and metastatic OS tissues, we showed that DLEC1, FOXJ1 and MAPK15 were also upregulated DEGs in metastatic OS tissues. Then, the mRNA and protein level of DLEC1, FOXJ1 and MAPK15 in OS tissues was detected. We found that both the mRNA and protein expression levels of FOXJ1 and MAPK15 were significantly upregulated in the OS tissues of patients with tumor metastasis at diagnosis.

MAPK15 is one of the most recently discovered MAPKs, but has been researched in a series of different tumor types (30). Wu et al (31) reported that MAPK15 was highly expressed in several human lung cancer cell lines; this study also demonstrated that siRNA-silencing of MAPK15, or NF-kB inhibition, reduced the sensitivity of lung cancer cells to arsenic trioxide, suggesting that MAPK15 regulates the drug resistance of cancer cells via the NF-kB pathway. Xu et al (32) found that the overexpression of MAPK15 could further increase the trans-activation of activator protein-1 by promoting the phosphorylation of proto-oncogene c-Jun, thus inducing the proliferation and transformation of cancer cells. Similarly, Jin DH (and a number of others) have also found that the MAPK15 copy number and level of expression in gastric cancer tissues were significantly increased, compared with non-gastric cancer tissues; additionally, siRNA knockdown of MAPK15 inhibited the proliferation of gastric cancer cells, leading to cell cycle arrest G1/S phase (19). Consistently, the present study confirmed that MAPK15 also plays an oncogenic role in OS, and that its inhibition significantly decreased OS cell metastasis in vitro and in vivo.

MMPs are involved in matrix degradation and play key roles in tumor growth, invasion and angiogenesis (33,34). c-Jun is a major transcription factor, which regulates the expression of MMPs (35). The role of the c-Jun/MMP pathways in OS has been widely reported; by activating JNK/c-Jun/MMP-2 signaling, the oncogene Astrocyte elevated gene-1 promotes the development of OS (36). Cryptochrome 2, an OS suppressor, inhibited the proliferation of OS cells by decreasing the expression of c-Jun/MMPs (37). Similarly, c-Jun/MMPs were also the target for OS therapy; by inhibiting c-Jun-associated pathways, glabridin effectively inhibited the proliferation of OS cells (38). Likewise, nobiletin inhibited the proliferation, migration and invasion capacities of OS cells by blocking the c-Jun/ MMP signaling pathways (39). In colorectal cancer, MAPK15 binds to and increases the level of c-Jun phosphorylation. Consistently, the present study illustrated that MAPK15 inhibition significantly decreased the level of c-Jun phosphorylation, as well as that of its target genes, including MMP2 and 9; c-Jun inhibition blocked the effects of MAPK15 overexpression on the migration and invasion of OS cells.

In the present study, WGCNA and DEG analysis of the GSE87624 dataset revealed that DLEC1, FOXJ1 and MAPK15 were predicted hub genes for the metastasis of OS, and that MAPK15 was most significantly upregulated in OS tissues. A series of in vitro and in vivo experiments subsequently confirmed that MAPK15 promoted the proliferation and metastasis of OS cells and tumors, respectively, by activating the c-Jun/MMP signaling pathways. Therefore, MAPK15 may be a novel biomarker for the diagnosis of OS, as well as an effective therapeutic target.

Supplementary Material

Supporting Data

Supporting Data

Supporting Data

Acknowledgements

Not applicable.

Funding

The present study was funded by the Science and Technology Program of Guangzhou, China (grant no. 201704020129).

Availability of data and materials

The datasets used and/or analyzed during the present study are available from the corresponding author on reasonable request.

Authors' contributions

ZS and LL conceived and designed the study. ZS, BY and ZZ collected the data, performed the associated related assays and wrote the manuscript. SZ, SL and CW performed data analysis and validation. YJ contributed to the study design and proofread the manuscript, and all authors read and approved the final manuscript.

Ethics approval and consent to participate

The OS tissue related study was approved by the Clinical Ethics Management Committee of Zhujiang Hospital, South Medical University, and all patients provided informed consent in writing. The animal care and experimental procedures used in the present study were approved by the Institutional Animal Care and Use Committee of the Southern Medical University.

Patient consent for publication

Not applicable

Competing interests

The authors declare that they have no competing interests.

References