Skip to main content
Download PDF

Development and analysis of long non-coding RNA-associated competing endogenous RNA network for osteosarcoma metastasis

Hereditas volume 158, Article number: 9 (2021) Cite this article

Abstract

Background

Osteosarcoma is the primary bone malignant neoplasm that often develops metastasis. Increasing evidences have shown that non-coding RNAs (ncRNAs) relate to the progression of osteosarcoma. However, the ncRNAs’ roles in osteosarcoma metastasis are still unknown.

Methods

Differentially expressed (DE) RNAs were identified from Gene Expression Omnibus (GEO) database. Protein-protein interaction (PPI) of DE messenger RNAs (DEmRNAs) was built through STRING database. The target mRNAs and long ncRNAs (lncRNAs) of microRNAs (miRNA) were predicted through miRDB, Targetscan and Genecode databases, which then cross-checked with previously obtained DERNAs to construct competing endogenous RNA (ceRNA) network. All networks were visualized via Cytoscape and the hub RNAs were screened out through Cytoscape plug-in Cytohubba. The gene functional and pathway analyses were performed through DAVID and MirPath databases. The survival analyses of hub RNAs were obtained through Kaplan-Meier (KM) survival curves.

Results

Five hundred sixty-four DEmRNAs, 16 DElncRNAs and 22 DEmiRNAs were screened out. GO functional and KEGG pathway analyses showed that DERNAs were significantly associated with tumor metastasis. The ceRNA network including 6 lncRNAs, 55 mRNAs and 20 miRNAs were constructed and the top 10 hub RNAs were obtained. Above all, PI3K/AKT signaling pathway was identified as the most important osteosarcoma metastasis-associated pathway and its hub ceRNA module was constructed. The survival analyses showed that the RNAs in hub ceRNA module closely related to osteosarcoma patients’ prognosis.

Conclusions

The current study provided a new perspective on osteosarcoma metastasis. More importantly, the RNAs in hub ceRNA module might act as the novel therapeutic targets and prognostic factors for osteosarcoma patients.

Introduction

Osteosarcoma, one of the most common primary musculoskeletal system malignant tumors, often affects children and young adults [1]. With the development of multidiscipline therapy, especially complete resection of tumor with multi-agent chemotherapy, the 5-year survival rate of patient has significantly improved to over 70% [2]. However, a large number of patients who receive standard or aggressive treatment still develop metastasis and their 5-year survival rates dramatically decrease to less than 20% [3, 4]. Therefore, the promising strategy to improve the prognosis of metastasis patients has been widely investigated. Unfortunately, little progression has been acquired since the underlying mechanisms of invasion and metastasis are still unclear [5].

Non-coding RNAs (ncRNAs) are a group of RNAs that have no protein-coding functions and play important roles in cellular physiology and pathology. Based on the length of nucleotide, ncRNAs are divided into small ncRNAs (< 200 bp) and long ncRNAs (200 to 100 kb). MicroRNAs (miRNAs), one of the most widely studied small ncRNAs, are associated with cell proliferation, differentiation, apoptosis [6]. After gene transcription, miRNAs could bind to 3’UTR (untranslated regions) of targeted messenger RNAs (mRNAs) and prevent their translation, which play vital roles in tumor occurrence and metastasis [7]. Unlike miRNAs, long ncRNAs (lncRNAs) take part in epigenetic, transcriptional, post-transcriptional and translational processes in various cells [8]. A large number of studies have revealed that lncRNAs correlate with tumor progression and act as the biomarkers for diagnosis and prognosis. Some researches even find lncRNAs would be optimistic therapeutic targets of cancers [9,10,11,12,13].

The competing endogenous RNA (ceRNA) hypothesis is proposed as a novel regulatory mechanism between coding RNAs and ncRNAs. Both lncRNAs and miRNAs are important components in this network. LncRNAs can bind to miRNAs through miRNA-binding sites (MREs) and act as miRNA sponges to regulate mRNA expression [14]. Plenty of researchers have constructed different ceRNA networks in various cancers. Lei et al analyzed the data from The Cancer Genome Atlas (TCGA) database and built up a lncRNA-miRNA-mRNA regulatory network in muscle-invasive bladder cancer [15]. Fang et al constructed a ceRNA network of head and neck squamous cell carcinoma from 546 patients’ profiles [16]. However, studies of ceRNA networks in osteosarcoma patients are limited and the role of lncRNA-miRNA-mRNA regulatory network between metastatic and non-metastatic osteosarcoma is still unclear.

In the present study, we obtained the differentially expressed (DE) miRNAs, mRNAs and lncRNAs between metastatic and non-metastatic patients from GSE39040 and GSE39055 datasets. The online databases was used to predict target mRNAs and lncRNAs of miRNAs. Then the ceRNA regulatory network was constructed. We also performed Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway and protein-protein interaction (PPI) network analyses of DEmiRNAs and DEmRNAs. Moreover, we identified the osteosarcoma metastasis-associated hub ceRNA module based on the KEGG pathway analysis. Kaplan-Meier (KM) survival analysis was conducted to verify the relationships between hub ceRNAs and prognosis. We hope that the present study would extend our knowledge of osteosarcoma metastasis-associated ceRNA network and provide potential therapeutic targets and prognostic markers for osteosarcoma patients.

Materials and methods

Microarray dataset selection

The expression profiles of osteosarcoma were searched from the Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo). The microarray datasets involved in this study should meet the following criteria: (1) the tumors were confirmed as osteosarcoma by histology; (2) the datasets contained metastatic and non-metastatic patients; (3) the datasets had complete prognosis-associated information; (4) the sample size in the dataset was more 10. Finally, two GEO datasets were selected for further study. The miRNA data were downloaded from GSE39040(platform GPL15762, Illumina Human v2 MicroRNA Expression BeadChip, Illumina, Inc., San Diego, CA, USA) which contained 65 osteosarcoma patients. The mRNA and lncRNA data were downloaded from GSE39055(platform GPL14951, Illumina HumanHT-12 WG-DASL V4.0 R2 Expression BeadChip, Illumina, Inc., San Diego, CA, USA) which had 37 osteosarcoma patients. The clinical characteristics of patients in GSE39040 and GSE39055 were shown in Supplementary Table 1.

Differential expression analysis

The probe names in the microarray datasets were converted into RNA names and multiple values of the same gene were averaged by R software (version 3.61). Patients from selected datasets were divided into metastatic and non-metastatic groups. The Linear Models for Microarray Analysis (Limma) package was applied to perform DERNAs analysis between two groups. The cutoff value of DERNAs was P<0.05 and |fold change (FC)|>1.5. The heat map was drawn through ‘pheatmap’ package in R software. According to the Gene database information (https://www.ncbi.nlm.nih.gov/gene) and classification of RefSeq accession numbers and molecule types (https://www.ncbi.nlm.nih.gov/books/NBK21091/table/ch18.T.refseq_accession_numbers_and_mole/?report=objectonly), the DERNAs in GSE39055 were divided into DEmRNA (‘NM_’ or ‘XM_’) and DElncRNA (‘NR_’ or ‘XR_’) groups. The flow chart of this study was shown in Fig. 1.

Fig. 1
figure 1

Flow chart of the study. GO:gene ontology; DEmRNA:differentially expressed messenger RNA; DEmiRNA:differentially expressed microRNA; DElncRNA:differentially expressed long non-coding RNA; KEGG:Kyoto Encyclopedia of Genes and Genomes; ceRNA:competing endogenous RNA; PPI:Protein-protein interaction

Full size image

Prediction of miRNA targets

The target mRNAs of DEmiRNA were screened from Targetscan (http://www.targetscan.org/vert_72) and miRDB (http://mirdb.org/) respectively [17, 18]. Then the intersectional mRNAs in two databases were obtained for further study. The target lncRNAs of DEmiRNA were selected from Genecode (https://www.gencodegenes.org/) by overlapping miRNA’s seed region to lncRNAs.

Functional and pathway enrichment analyses

GO functional and KEGG pathway enrichment analyses of all DEmRNAs and mRNAs in ceRNA network were performed by Database for Annotation, Visualization and Integrated Discovery (DAVID, https://david.ncifcrf.gov/) [19]. GO functional enrichment analysis included cellular component (CC), molecular function (MF), and biological process (BP) [20]. MirPath (http://snf-515788.vm.okeanos.grnet.gr) was used to elucidate the GO functional and KEGG pathway enrichment of DEmiRNAs [21]. P-value<0.05 was defined as significantly enriched.

Protein-protein interaction (PPI) network construction and analysis

To study the interaction of all DEmRNAs, the online tool Search Tool for the Retrieval of Interacting Genes database (STRING, https://string-db.org/) was used to construct PPI network [22]. The cutoff confidence score was 0.4 and the disconnected nodes were hid in the network. Then the network was visualized in Cytoscape(v3.7.1) and the top 10 hub mRNAs were screened out by Cytoscape plug-in CytoHubba. KEGG pathway analysis of hub mRNAs was also conducted by DAVID and p-value<0.05 was defined as significantly enriched.

ceRNA network construction and analysis

The overlapping mRNAs were obtained from the intersection of DEmRNAs from GSE39055 and target mRNAs of DEmiRNA. The same method was used to identify overlapping lncRNAs. Then the inverse expression pairs of miRNA-mRNA and miRNA-lncRNA were selected to construct ceRNA network via Cytoscape. Top 10 hub RNAs in ceRNA network were also screened out by CytoHubba. The KEGG pathway analysis of miRNAs in hub RNAs was conducted through their targeted mRNAs in ceRNA network by DAVID. P-value<0.05 was defined as significantly enriched. Moreover, the most important osteosarcoma metastasis-associated hub ceRNA module was identified and constructed through KEGG pathway analysis results. The hub ceRNA module was visualized by Cytoscape.

Survival analysis

The prognostic information of osteosarcoma patients was obtained from the selected dataset. All Patients were divided into high- and low-expression groups based on the expression level of RNAs in hub ceRNA module. The Kaplan-Meier survival analysis of RNAs in hub ceRNA module was conducted by ‘survival’ package in R software. P < 0.05 was considered as statistically significant.

Results

Identification of DEmiRNAs, DEmRNAs and DElncRNAs

The miRNA expression data were obtained from the GSE39040 dataset, which contained 24 metastatic and 41 non-metastatic osteosarcoma patients. A total of 41 DEmiRNAs were found and the distribution characteristics of these genes were shown in Fig. 2a. 22 of 41 DEmiRNAs IDs could convert into MiRBase IDs for further study and the information of these DEmiRNAs were shown in Table 1.

Fig. 2
figure 2

The heat map of DEmiRNAs in GES39040(a) and DERNAs in GES39055(b). Red and greed indicated up and down-regulated genes between metastasis and non-metastasis groups. M:metastasis; NM: non-metastasis; FC:fold change; DEmiRNA:differentially expressed microRNA

Full size image
Table 1 The information of all selected DEmiRNAs in GSE39040
Full size table

The GSE39055 dataset, which contained 19 metastatic and 18 non-metastatic osteosarcoma patients, was used to analyze DEmRNAs and DElncRNAs. A total of 600 DERNAs were identified and the heat map was used to visualize the different expression values (Fig. 2b). Five hundred sixty-four demrnas and 16 DElncRNAs in these DERNAs were identified for further study. The top 10 up and down-regulated DEmRNAs and all DElncRNAs were listed in Tables 2 and 3.

Table 2 The information of top 10 up and down regulation DEmRNAs in GES39055
Full size table
Table 3 The information of all DElncRNAs inGES39055
Full size table

GO functional and KEGG pathway analyses of DEmiRNAs and DEmRNAs

In order to understand the potential biological functions of these DEmiRNAs, GO functional and KEGG pathway analyses were conducted by MirPath. In BP ontology, the DEmiRNAs were involved in several tumor metastatic processes such as biosynthetic process, cellular protein modification process, gene expression and catabolic process (Fig. 3a). As for CC analysis, the DEmiRNAs were mainly distributed in organelle, nucleoplasm, protein complex, cytosol, microtubule organizing center (Fig. 3b). In MF category, DEmiRNAs were encriched in ion binding, nucleic acid binding transcription factor activity, enzyme binding, protein binding transcription factor activity, enzyme regulator activity, cytoskeletal protein binding, RNA binding and transcription corepressor activity (Fig. 3c). Most of these functions took part in tumor metastasis processes through regulating gene expression, protein binding or enzyme activity [23,24,25].

Fig. 3
figure 3

GO functional and KEGG pathway analyses of DEmiRNA. a top 10 significantly enriched GO biology process; b7 significantly enriched GO cellular components; c top 10 significantly enriched molecular function; d top 10 significantly enriched KEGG pathways. GO:gene ontology;KEGG:Kyoto Encyclopedia of Genes and Genomes;DEmRNA: differentially expressed messenger RNA

Full size image

In addition, the KEGG pathway analysis indicated that DEmiRNAs enriched in 53 different pathways (Supplementary Table 2) and the top 10 pathways were shown in Fig. 3d. A large amount of common tumor metastasis-related pathways such as TGF-beta signaling pathway, Hippo signaling pathway, Ras signaling pathway, Wnt signaling pathway and PI3K/Akt signaling pathway were involved.

The online tool DAVID was used to perform GO functional and KEGG pathway analysis for DEmRNAs. The significant result of GO functional enrichment was shown in Fig. 4a. Plenty of items, such as transcription corepressor activity, cAMP response element binding, protein serine/threonine/tyrosine kinase activity were closely related to the process of tumor metastasis [26,27,28].

Fig. 4
figure 4

GO functional and KEGG pathway analyses of DEmRNA in GSE39055. a significantly enriched GO items in CC, BP and MF; b significantly enriched KEGG pathways. GO:gene ontology; KEGG:Kyoto Encyclopedia of Genes and Genomes; FC:fold change; DEmRNA: differentially expressed messenger RNA;BP: biology process;CC:cellular component; MF: molecular function

Full size image

KEGG pathway analysis showed that DEmRNAs significantly enriched in 8 pathways (Fig. 4b). Several tumor metastasis-associated pathways such as PI3K/Akt signaling pathway and Hippo signaling pathway were included.

PPI network construction and analysis

The online tool STRING was used to identify the correlation of DEmRNAs and construct the PPI network. A total of 564 nodes and 900 edges were involved in the PPI network. The disconnected nodes were deleted and the PPI network was input to Cytoscape for further study (Fig. 5). The Cytoscape plug-in CytoHubba was used to screen out the top 10 hub mRNAs, which included IL6, ITGB2, PPP2CA, HNRNPC, CPSF1, YWHAB, RPS6KB1, UBA7, MAP 2 K1 and RPS3A. The KEGG pathway enrichment of these hub mRNAs was performed in DAVID and the result was shown in Table 4. It is worth noting that the tumor metastasis-associated PI3K/Akt signaling pathway was included once again.

Fig. 5
figure 5

Protein-protein interaction network of DEmRNA in GSE39055. Red nodes represent up-regulated mRNAs, blue nodes represent down-regulated mRNAs.DEmRNA: differentially expressed messenger RNA;FC:fold change

Full size image
Table 4 The KEGG pathway enrichment of hub DEmRNAs
Full size table

ceRNA network construction and analysis

Targetscan and miRDB were used to predict the target mRNAs of DEmiRNAs just as previously mentioned. In total, 4030 targeted mRNAs were screened out. Then these mRNAs were cross-checked with DEmRNAs in GSE39055 and the intersectional mRNAs which presented with an opposite expression trend to DEmiRNAs were selected. After that, 61 mRNA-miRNA regulatory relationships including 16 miRNAs and 55 mRNAs were obtained. The same method was used to establish the lncRNA-miRNA network through Genecode and 36 lncRNA-miRNA regulatory relationships including 6 lncRNAs and 15miRNAs were screened out. Finally, the ceRNA network which contained 6 lncRNAs,20 miRNAs and 55 mRNAs were constructed through Cytoscape (Fig. 6a).

Fig. 6
figure 6

ceRNA network of DERNAs. a ceRNA network of all DERNAs; b the osteosarcoma metastasis-associated hub ceRNA module. Red nodes represent upregulated RNAs,blue nodes represent downregulated RNAs.DERNAs: differentially expressed RNA;FC:fold change;mRNA:messenger RNA;miRNA: microRNA; lncRNA: long non-coding RNA; ceRNA: competing endogenous RNA

Full size image

GO functional and KEGG pathway analyses of mRNAs in ceRNA network were conducted through DAVID. The GO functional enrichment result was shown in Fig. 7a and several biological processes such as regulation of growth, positive regulation of transcription from RNA polymerase II promoter, mitotic nuclear envelope reassembly had a significant impact on tumor metastasis [29, 30]. The outcome of KEGG pathway enrichment of mRNAs was shown in Fig. 7b and Supplementary Table 3. The tumor metastasis-associated pathways such as PI3K/Akt signaling pathway, TGF-beta signaling pathway, AMPK signaling pathway and HIF-1 signaling pathway were significantly enriched.

Fig. 7
figure 7

GO functional and KEGG pathway analyses of mRNAs in ceRNA network. a significantly enriched GO items in CC, BP and MF; b top 10 significantly enriched KEGG pathways. GO:gene ontology; KEGG:Kyoto Encyclopedia of Genes and Genomes; FC:fold change; DEmRNA: differentially expressed messenger RNA;BP: biology process;CC:cellular component; MF: molecular function

Full size image

The top 10 hub RNAs in ceRNA network were screened out through Cytoscape plug-in CytoHubba. All hub RNAs were lncRNA and miRNA, including TPT1-AS1, hsa-miR-612, hsa-miR-1304, MIR137HG, MRPL20-AS1, hsa-miR-1183, KMT2E-AS1, hsa-miR-885-3p, hsa-miR-202*:9.1 and hsa-miR-1268. Then the KEGG pathway analysis of these miRNAs was conducted through their targeted mRNAs in ceRNA network (Table 5). The tumor metastasis-associated PI3K/Akt signaling pathway was enriched once again and we identified it as the most important osteosarcoma metastasis-associated pathway because it was the only pathway that significantly enriched in every KEGG analysis result. Based on the regulatory relationships between hub RNAs and mRNAs related to PI3K/Akt signaling pathway, the osteosarcoma metastasis-associated hub ceRNA modules were constructed (Fig. 6b).

Table 5 The KEGG pathway enrichment of hub miRNA targeted mRNAs in ceRNA network
Full size table

Survival analysis of RNAs in osteosarcoma metastasis-associated hub ceRNA modules

The prognosis information of patents was download from the same datasets. The KM survival curve analysis was used to identify the relation between the expression levels of RNAs in hub ceRNA module and the progression-free survival rate (PFS) of patients (Fig. 8). The results showed that different expression levels of these RNAs had a distinct impact on patients’ PFS rate and most of them were statistically or borderline (8/11) significant, which indicated that these RNAs might be the promising prognostic biomarkers for osteosarcoma patients. What’s more, we searched these RNAs in pubmed database and the results also showed that they were correlated with osteosarcoma progression and prognosis [31,32,33,34].

Fig. 8
figure 8

Kaplan-Meier survival curves of RNAs in hub ceRNA module. a-d: mRNA; e-h: lncRNA; i-k: miRNA.mRNA:messenger RNA; miRNA: microRNA; lncRNA: long non-coding RNA

Full size image

Discussion

Osteosarcoma is one of the most common bone malignant tumors and about one-fifth of patients appear clinically detectable metastases at diagnosis. The most common metastatic site of osteosarcoma is pulmonary and the 5-year survival rate of these patients is extremely low [35]. In order to improve the prognosis of patients with metastasis, the mechanism underlying metastasis was studied and first mentioned by Paget in 1889, just known as the “seed and soil” theory [36]. After that, with the development of science and technology, more theories of cancer metastasis had been proposed, such as epithelial-mesenchymal transition (EMT) enhancement, gene fusion and mutation, intravasation, enhancement of cancer stemness, increased angiogenesis and local invasion [37]. Unfortunately, the primary mechanism of osteosarcoma metastasis was still unclear and the survival time of metastatic patients was not significantly prolonged in decades.

In the last few years, the roles of lncRNA in various cancers had been widely studied. One important function of lncRNA was to bind to miRNAs and act as the miRNA sponge to regulate mRNA expression, which was correlated with tumor formation, progression and metastasis. These lncRNA-miRNA-mRNA regulatory relationships, also known as ceRNA network, had been reported in different cancers [15, 16]. However, the ceRNA network in metastatic osteosarcoma had not been revealed. Therefore, we conducted this study and looked forward to identify the key RNAs involved in the process of metastasis.

In the present study, we got 41 DEmiRNAs from GSE39040 and 22 of them were selected for further study. 564 DEmRNAs and 16 DElncRNAs were identified from GSE39055. The GO functional analysis of these aberrantly expressed RNAs showed that they were closely related to the process of cancer metastasis, such as protein binding transcription factor activity, enzyme regulator activity, protein serine/threonine/tyrosine kinase activity and so on. As for the KEGG pathway enrichment analysis, some common cancer metastasis-associated pathways were significantly enriched, including PI3K/Akt signaling pathway, Hippo signaling pathway, Ras signaling pathway and Wnt signaling pathway. The PPI network of DEmRNAs was constructed and the hub mRNAs were screened out. These hub mRNAs were closely related to the metastatic processes in various cancers. For instance, the overexpression of miR-449a directly inhibited the translation of mitogen-activated protein (MAP) kinase kinase MAP 2 K1, which decreased the activity of MEK1/ERK1/2/c-Jun pathway and suppressed the invasion and metastasis of lung cancer cells [38]. IL6, one of the proinflammatory cytokines, promoted breast cancer growth and metastasis through IL6/JAK/STAT3 loop [39]. The transcriptional co-activator Yes-associated protein (YAP) induced the expression of integrin subunit ITGB2, which increased cancer cell invasion [40]. In addition, the KEGG pathway analysis revealed that the hub mRNAs were significantly enriched in PI3K/Akt signaling pathway once again. To sum up, these analyses showed that the DERNAs obtained from two datasets played important roles in the process of cancer metastasis and the results of our study were reliable.

Next, we predicted the target mRNAs and lncRNAs of DEmiRNAs from online databases and cross-checked them with previously obtained DEmRNAs or DElncRNAs. After that, the osteosarcoma metastasis-related ceRNA network was constructed and the KEEG pathway analysis showed that the mRNAs in ceRNA network were enriched in PI3K/Akt signaling pathway once again. At the same time, the hub RNAs in ceRNA network were screened out and the miRNAs among them were enriched in PI3K/Akt signaling pathway again. Based on the results showed above, the PI3K/Akt signaling pathway was treated as the most important osteosarcoma metastasis-associated pathway because apart from it, no other pathways could be significantly enriched every time. Next, based on the PI3K/Akt signaling pathway associated hub RNAs and their targeted mRNAs, the hub osteosarcoma metastasis-associated ceRNA module was constructed. At last, the KM survival analysis showed that the vast majority of RNAs in hub ceRNA module had a significant impact on the prognosis of osteosarcoma patients.

The PI3K/AKT signaling pathway is one of the most crucial intracellular oncogenic pathways in cancers. It could be activated by G-protein-coupled receptor signaling and receptor tyrosine kinases (RTKs) such as vascular endothelial growth factor receptor (VEGFR), platelet-derived growth factor receptor (PDGFR) and insulin-like growth factor 1 receptor (IGF1R). The activation or inhibition of PI3K/AKT signaling pathway would result in a cascade of events that played important roles in pathological processes, such as oncogenesis, invasion, metastasis, angiogenesis and proliferation [41]. What’s more, almost every component of this pathway had amplifications or mutations in cancers. In the present study, we identified PI3K/AKT as the most important signaling pathway in osteosarcoma metastasis. A large number of studies had reported the same results. Zhang et al found that the Fibulin-4 could induce EMT by activating PI3K/AKT signaling pathway, which promoted osteosarcoma cell invasion and metastasis [42]. Jin et al reported that suppressing the Glycoprotein non-metastatic melanoma protein B (GPNMB) would inhibit osteosarcoma cell proliferation and metastasis via blocking PI3K/Akt/mTOR signaling pathway [43]. In the last few years, numerous studies had also revealed that ncRNAs could regulate osteosarcoma metastasis through PI3K/AKT signaling pathway. LINC00628, one of the common lncRNAs, was decreased in osteosarcoma tissues and overexpression of it would inhibit the proliferation, invasion and migration of osteosarcoma cells through PI3K/AKT signaling pathway [44]. To sum up, the PI3K/AKT signaling pathway is involved in the process of osteosarcoma metastasis and alterations in this pathway might improve the prognosis of patients.

With the development of microarray and next generation sequencing technology, numerous aberrantly expressed oncogenes were detected and most of them related to the occurrence and development of cancers. On the basis of these findings, we constructed the hub ceRNA module associated with PI3K/AKT signaling pathway, which contained three miRNAs (hsa-miR-1304, hsa-miR-202*:9.1 and hsa-miR-1183), four mRNAs (PPP2R2A, ATF2, PPP2CA and RPS6KB1) and four lncRNAs (MIR137HG, TPT1-AS1, KMT2E-AS1 and MRPL20-AS1). A number of studies found that the RNAs in our ceRNA network correlated with the proliferation, invasion and metastasis of various cancers. PPP2CA, the catalytic subunit of PP2A, was higher expressed in osteosarcoma tissues. The reduction of PPP2CA suppressed the proliferation and metastasis of osteosarcoma both in vitro and in vivo [34]. RPS6KB1 was a conserved serine/threonine kinase that served as the transcription factor in protein synthesis process. Cai et al found that the overexpression of RPS6KB1 in prostate cancer patients related to aggressive progression and poor prognosis. Meanwhile, suppressing RPS6KB1 activity by miR-195 inhibited prostate cancer invasion and migration, which provided a potential target for prostate cancer treatment [45]. MiR-1304, usually regarded as a tumor suppressor, inhibited non-small cell lung cancer (NSCLC) cell growth through suppressing heme oxygenase-1 (HO-1) and was down-regulated in NSCLC samples [46]. In a word, the genes in hub ceRNA module influenced the metastasis of cancers and might act as the therapeutic targets. However, the roles of these genes in osteosarcoma were still limited and more studies should be conducted in the future.

The KM survival analysis of our study showed that RNAs in hub ceRNA network were correlated with the patients’ survival rate, which indicated that they could serve as the predictors of prognosis. PPP2R2A, the upstream effectors of Akt/mTOR pathway, encoded an alpha isoform of the regulatory subunit B55 subfamily. Our research showed that osteosarcoma patients with higher expression of PPP2R2A had a worse prognosis, which was consistent with the study by Hein et al. [47] Zhou et al reported that the expression level of RPS6KB1 encoded protein p70S6K was negatively correlated with the progression of osteosarcoma patients. The same result was obtained in our KM survival analysis [33].

In conclusion, the present study screened out DElnRNAs, DEmiRNAs and DEmRNAs between metastatic and non-metastatic osteosarcoma patients through microarray datasets in GEO database. Importantly, the ceRNA network had been constructed and the PI3K/AKT signaling pathway was identified as the most important osteosarcoma metastasis-associated signaling pathway. The hub ceRNA module including three miRNAs (hsa-miR-1304, hsa-miR-202*:9.1 and hsa-miR-1183), four mRNAs (PPP2R2A, ATF2, PPP2CA and RPS6KB1) and four clncRNAs (MIR137HG, TPT1-AS1, KMT2E-AS1 and MRPL20-AS1) was constructed. Our KM survival analysis and other researchers’ results showed that the RNAs in hub ceRNA module might act as the novel therapeutic targets and prognostic markers for osteosarcoma patients. Based on the present study, we would continue to conduct researches in the future. It might reveal the potential mechanisms underlying osteosarcoma metastasis and improve the prognosis of osteosarcoma patients.

Availability of data and materials

The authors declare that the data supporting the findings of this study are available within the article.

Abbreviations

ncRNA:

Non-coding RNA

DERNA:

Differentially expressed RNA

GEO:

Gene Expression Omnibus

PPI:

Protein-protein interaction

mRNA:

Messenger RNA

DEmRNA:

Differentially expressed messenger RNA

lncRNA:

Long non-coding RNA

miRNA:

microRNA

ceRNA:

Competing endogenous RNA

KM:

Kaplan-Meier

UTR:

Untranslated region

MRE:

miRNA-binding sites

TCGA:

The Cancer Genome Atlas

GO:

Gene Ontology

KEGG:

Kyoto Encyclopedia of Genes and Genomes

Limma:

Linear Models for Microarray Analysis

DAVID:

Database for Annotation, Visualization and Integrated Discovery

FC:

Fold change

CC:

Cellular component

MF:

Molecular function

BP:

Biological process

STRING:

Search Tool for the Retrieval of Interacting Genes database

EMT:

Epithelial-mesenchymal transition

MAP:

Mitogen-activated protein

YAP:

Yes-associated protein

RTK:

Receptor tyrosine kinase

VEGFR:

Vascular endothelial growth factor receptor

PDGFR:

Platelet-derived growth factor receptor

IGF1R:

Insulin-like growth factor 1 receptor

GPNMB:

Glycoprotein non-metastatic melanoma protein B

NSCLC:

Non-small cell lung cancer

HO-1:

Heme oxygenase-1

References

  1. Picci P. Osteosarcoma (osteogenic sarcoma). Orphanet journal of rare diseases. 2007;2:6.

    Article  PubMed  PubMed Central  Google Scholar 

  2. Anderson ME. Update on Survival in Osteosarcoma. The Orthopedic clinics of North America. 2016;47(1):283–92.

    Article  PubMed  Google Scholar 

  3. Khanna C, Fan TM, Gorlick R, Helman LJ, Kleinerman ES, Adamson PC, et al. Toward a drug development path that targets metastatic progression in osteosarcoma. Clinical cancer research : an official journal of the American Association for Cancer Research. 2014;20(16):4200–9.

    Article  Google Scholar 

  4. Kansara M, Teng MW, Smyth MJ, Thomas DM. Translational biology of osteosarcoma. Nat Rev Cancer. 2014;14(11):722–35.

    Article  CAS  PubMed  Google Scholar 

  5. Cao J, Wang Y, Dong R, Lin G, Zhang N, Wang J, et al. Hypoxia-induced WSB1 promotes the metastatic potential of osteosarcoma cells. Cancer Res. 2015;75(22):4839–51.

    Article  CAS  PubMed  Google Scholar 

  6. Sana J, Faltejskova P, Svoboda M, Slaby O. Novel classes of non-coding RNAs and cancer. J Transl Med. 2012;10:103.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  7. Di Leva G, Garofalo M, Croce CM. MicroRNAs in cancer. Annu Rev Pathol. 2014;9:287–314.

    Article  PubMed  CAS  Google Scholar 

  8. Khorkova O, Hsiao J, Wahlestedt C. Basic biology and therapeutic implications of lncRNA. Adv Drug Deliv Rev. 2015;87:15–24.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  9. Nie FQ, Sun M, Yang JS, Xie M, Xu TP, Xia R, et al. Long noncoding RNA ANRIL promotes non-small cell lung cancer cell proliferation and inhibits apoptosis by silencing KLF2 and P21 expression. Mol Cancer Ther. 2015;14(1):268–77.

    Article  CAS  PubMed  Google Scholar 

  10. Parasramka MA, Maji S, Matsuda A, Yan IK, Patel T. Long non-coding RNAs as novel targets for therapy in hepatocellular carcinoma. Pharmacol Ther. 2016;161:67–78.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  11. Peng F, Shi X, Meng Y, Dong B, Xu G, Hou T, et al. Long non-coding RNA HOTTIP is upregulated in renal cell carcinoma and regulates cell growth and apoptosis by epigenetically silencing of LATS2. Biomedicine & pharmacotherapy = Biomedecine & pharmacotherapie. 2018;105:1133–40.

    Article  CAS  Google Scholar 

  12. Rupaimoole R, Lee J, Haemmerle M, Ling H, Previs RA, Pradeep S, et al. Long noncoding RNA Ceruloplasmin promotes Cancer growth by altering glycolysis. Cell Rep. 2015;13(11):2395–402.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  13. Wu W, Gao H, Li X, Zhu Y, Peng S, Yu J, et al. LncRNA TPT1-AS1 promotes tumorigenesis and metastasis in epithelial ovarian cancer by inducing TPT1 expression. Cancer Sci. 2019;110(5):1587–98.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  14. Salmena L, Poliseno L, Tay Y, Kats L, Pandolfi PP. A ceRNA hypothesis: the Rosetta stone of a hidden RNA language? Cell. 2011;146(3):353–8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  15. Lyu L, Xiang W, Zhu JY, Huang T, Yuan JD, Zhang CH. Integrative analysis of the lncRNA-associated ceRNA network reveals lncRNAs as potential prognostic biomarkers in human muscle-invasive bladder cancer. Cancer Manag Res. 2019;11:6061–77.

    Article  PubMed  PubMed Central  Google Scholar 

  16. Fang XN, Yin M, Li H, Liang C, Xu C, Yang GW, et al. Comprehensive analysis of competitive endogenous RNAs network associated with head and neck squamous cell carcinoma. Sci Rep. 2018;8(1):10544.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  17. Agarwal V, Bell GW, Nam JW, Bartel DP. Predicting effective microRNA target sites in mammalian mRNAs. eLife. 2015;4.

  18. Liu W, Wang X. Prediction of functional microRNA targets by integrative modeling of microRNA binding and target expression data. Genome Biol. 2019;20(1):18.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  19. Huang da W, Sherman BT, Lempicki RA. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protoc 2009;4(1):44–57.

  20. Consortium GO. The gene ontology (GO) project in 2006. Nucleic Acids Res. 2006;34(Database issue):D322–6.

    Article  CAS  Google Scholar 

  21. Vlachos IS, Zagganas K, Paraskevopoulou MD, Georgakilas G, Karagkouni D, Vergoulis T, et al. DIANA-miRPath v3.0: deciphering microRNA function with experimental support. Nucleic Acids Res. 2015;43(W1):W460–6.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  22. Szklarczyk D, Franceschini A, Wyder S, Forslund K, Heller D, Huerta-Cepas J, et al. STRING v10: protein-protein interaction networks, integrated over the tree of life. Nucleic Acids Res. 2015;43(Database issue):D447–52.

    Article  CAS  PubMed  Google Scholar 

  23. Wang Y, Kumar S, Rachagani S, Sajja BR, Xie Y, Hang Y, et al. Polyplex-mediated inhibition of chemokine receptor CXCR4 and chromatin-remodeling enzyme NCOA3 impedes pancreatic cancer progression and metastasis. Biomaterials. 2016;101:108–20.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  24. Heeg S, Das KK, Reichert M, Bakir B, Takano S, Caspers J, et al. ETS-Transcription Factor ETV1 Regulates Stromal Expansion and Metastasis in Pancreatic Cancer. Gastroenterology. 2016;151(3):540–53 e14.

    Article  CAS  PubMed  Google Scholar 

  25. Li Z, Takino T, Endo Y, Sato H. Activation of MMP-9 by membrane type-1 MMP/MMP-2 axis stimulates tumor metastasis. Cancer Sci. 2017;108(3):347–53.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  26. Valencia-Sama I, Zhao Y, Lai D, Janse van Rensburg HJ, Hao Y, Yang X. Hippo component TAZ functions as a co-repressor and negatively regulates DeltaNp63 transcription through TEA domain (TEAD) transcription factor. J Biol Chem. 2015;290(27):16906–17.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  27. Wang YW, Chen X, Gao JW, Zhang H, Ma RR, Gao ZH, et al. High expression of cAMP-responsive element-binding protein 1 (CREB1) is associated with metastasis, tumor stage and poor outcome in gastric cancer. Oncotarget. 2015;6(12):10646–57.

    Article  PubMed  PubMed Central  Google Scholar 

  28. Hu L, Chen HY, Cai J, Zhang Y, Qi CY, Gong H, et al. Serine threonine tyrosine kinase 1 is a potential prognostic marker in colorectal cancer. BMC Cancer. 2015;15:246.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  29. de Las Heras JI, Batrakou DG, Schirmer EC. Cancer biology and the nuclear envelope: a convoluted relationship. Semin Cancer Biol. 2013;23(2):125–37.

    Article  PubMed  CAS  Google Scholar 

  30. Zhang X, Chiang HC, Wang Y, Zhang C, Smith S, Zhao X, et al. Attenuation of RNA polymerase II pausing mitigates BRCA1-associated R-loop accumulation and tumorigenesis. Nat Commun. 2017;8:15908.

    Article  PubMed  PubMed Central  Google Scholar 

  31. Sun Z, Zhang T, Hong H, Liu Q, Zhang H. miR-202 suppresses proliferation and induces apoptosis of osteosarcoma cells by downregulating Gli2. Mol Cell Biochem. 2014;397(1–2):277–83.

    Article  CAS  PubMed  Google Scholar 

  32. Ma T, Liu A, Xu D, Zhang T. Mechanisms underlying the promotion of osteosarcoma cell proliferation and invasion by lncRNA PBB12. Oncol Rep. 2020;43(2):736–46.

    CAS  PubMed  Google Scholar 

  33. Zhou Q, Deng Z, Zhu Y, Long H, Zhang S, Zhao J. mTOR/p70S6K signal transduction pathway contributes to osteosarcoma progression and patients' prognosis. Medical oncology (Northwood, London, England). 2010;27(4):1239–45.

  34. Yang D, Okamura H, Morimoto H, Teramachi J, Haneji T. Protein phosphatase 2A Calpha regulates proliferation, migration, and metastasis of osteosarcoma cells. Laboratory investigation; a journal of technical methods and pathology. 2016;96(10):1050–62.

  35. Isakoff MS, Bielack SS, Meltzer P, Gorlick R. Osteosarcoma: current treatment and a collaborative pathway to success. Journal of clinical oncology : official journal of the American Society of Clinical Oncology. 2015;33(27):3029–35.

    Article  CAS  Google Scholar 

  36. Fokas E, Engenhart-Cabillic R, Daniilidis K, Rose F, An HX. Metastasis: the seed and soil theory gains identity. Cancer Metastasis Rev. 2007;26(3–4):705–15.

    Article  PubMed  Google Scholar 

  37. Han Y, Guo W, Ren T, Huang Y, Wang S, Liu K, et al. Tumor-associated macrophages promote lung metastasis and induce epithelial-mesenchymal transition in osteosarcoma by activating the COX-2/STAT3 axis. Cancer Lett. 2019;440–441:116–25.

  38. You J, Zhang Y, Li Y, Fang N, Liu B, Zu L, et al. MiR-449a suppresses cell invasion by inhibiting MAP 2K1 in non-small cell lung cancer. Am J Cancer Res. 2015;5(9):2730–44.

    CAS  PubMed  PubMed Central  Google Scholar 

  39. Chang Q, Bournazou E, Sansone P, Berishaj M, Gao SP, Daly L, et al. The IL-6/JAK/Stat3 feed-forward loop drives tumorigenesis and metastasis. Neoplasia (New York, NY). 2013;15(7):848–62.

    Article  CAS  Google Scholar 

  40. Liu H, Dai X, Cao X, Yan H, Ji X, Zhang H, et al. PRDM4 mediates YAP-induced cell invasion by activating leukocyte-specific integrin beta2 expression. EMBO reports. 2018;19:6.

    Google Scholar 

  41. Porta C, Paglino C, Mosca A. Targeting PI3K/Akt/mTOR Signaling in Cancer. Front Oncol. 2014;4:64.

    Article  PubMed  PubMed Central  Google Scholar 

  42. Zhang D, Wang S, Chen J, Liu H, Lu J, Jiang H, et al. Fibulin-4 promotes osteosarcoma invasion and metastasis by inducing epithelial to mesenchymal transition via the PI3K/Akt/mTOR pathway. Int J Oncol. 2017;50(5):1513–30.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  43. Jin R, Jin YY, Tang YL, Yang HJ, Zhou XQ, Lei Z. GPNMB silencing suppresses the proliferation and metastasis of osteosarcoma cells by blocking the PI3K/Akt/mTOR signaling pathway. Oncol Rep. 2018;39(6):3034–40.

    CAS  PubMed  Google Scholar 

  44. He R, Wu JX, Zhang Y, Che H, Yang L. LncRNA LINC00628 overexpression inhibits the growth and invasion through regulating PI3K/Akt signaling pathway in osteosarcoma. Eur Rev Med Pharmacol Sci. 2018;22(18):5857–66.

    CAS  PubMed  Google Scholar 

  45. Cai C, Chen QB, Han ZD, Zhang YQ, He HC, Chen JH, et al. miR-195 inhibits tumor progression by targeting RPS6KB1 in human prostate Cancer. Clinical Cancer Res. 2015;21(21):4922–34.

    Article  CAS  Google Scholar 

  46. Li CG, Pu MF, Li CZ, Gao M, Liu MX, Yu CZ, et al. MicroRNA-1304 suppresses human non-small cell lung cancer cell growth in vitro by targeting heme oxygenase-1. Acta Pharmacol Sin. 2017;38(1):110–9.

    Article  PubMed  CAS  Google Scholar 

  47. Hein AL, Seshacharyulu P, Rachagani S, Sheinin YM, Ouellette MM, Ponnusamy MP, et al. PR55alpha subunit of protein phosphatase 2A supports the tumorigenic and metastatic potential of pancreatic Cancer cells by sustaining hyperactive oncogenic Signaling. Cancer Res. 2016;76(8):2243–53.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

Download references

Acknowledgments

Not applicable.

Funding

This work was supported by the National Natural Science Foundation of China (NSFC, NO. 81702661), Youth Foundation of Shanghai Health and Family Planning Commission (NO. 2016Y0143).

Author information

Authors and Affiliations

  1. Department of Orthopedics, Ruijin Hospital, Shanghai Jiaotong University School of Medicine, 197 Ruijin Er Road, Shanghai, 200025, PR China

    Yucheng Fu, Qi Liu, Qiyuan Bao, Junxiang Wen, Zhuochao Liu, Yuehao Hu, Guoyu He, Cheng Peng, Yiqi Xu & Weibin Zhang

Authors
  1. Yucheng Fu
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Qi Liu
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Qiyuan Bao
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Junxiang Wen
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Zhuochao Liu
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Yuehao Hu
    View author publications

    You can also search for this author in PubMed Google Scholar

  7. Guoyu He
    View author publications

    You can also search for this author in PubMed Google Scholar

  8. Cheng Peng
    View author publications

    You can also search for this author in PubMed Google Scholar

  9. Yiqi Xu
    View author publications

    You can also search for this author in PubMed Google Scholar

  10. Weibin Zhang
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

YCF and WBZ designed the study. QL and YHH searched the data from database. YCF, JXW and QYB performed in silico analysis of the data. YCF and ZCL wrote the manuscript. GYH and WBZ revised the manuscript. CP and YQX modified language. All authors have read and approved the final manuscript.

Corresponding author

Correspondence to Weibin Zhang.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1 Supplementary Table S1.

The clinical characteristics of patients in GSE39040 and GSE39055. Supplementary Table S2.The KEGG pathway enrichment of DEmiRNAs. Supplementary Table S3.The KEGG pathway enrichment of mRNAs in ceRNA network.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Fu, Y., Liu, Q., Bao, Q. et al. Development and analysis of long non-coding RNA-associated competing endogenous RNA network for osteosarcoma metastasis. Hereditas 158, 9 (2021). https://doi.org/10.1186/s41065-021-00174-0

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s41065-021-00174-0

Keywords

Hereditas

ISSN: 1601-5223

Contact us