Comprehensive Analysis of a ceRNA Network Identifies lncR-C3orf35 Associated with Poor Prognosis in Osteosarcoma

Shi, Yi; Ren, Jianhua; Zhuang, Ze; Zhang, Wenhui; Wang, Zhe; Liu, Yuangao; Li, Jinze; Liang, Tangzhao; He, Ronghan; Wang, Kun

doi:https://doi.org/10.1155/2020/3178037

BioMed Research International

On this page

Abstract Introduction Materials and Methods Results Discussion Data Availability Conflicts of Interest Authors’ Contributions Acknowledgments References Copyright Related Articles

Special Issue

Bioinformatics, Experimental and Computational Biology for Personalized Medicine in Chronic Diseases

View this Special Issue

Research Article | Open Access

Volume 2020 | Article ID 3178037 | https://doi.org/10.1155/2020/3178037

Comprehensive Analysis of a ceRNA Network Identifies lncR-C3orf35 Associated with Poor Prognosis in Osteosarcoma

Yi Shi,^1,2Jianhua Ren,¹Ze Zhuang,¹Wenhui Zhang,¹Zhe Wang,¹Yuangao Liu,¹Jinze Li,¹Tangzhao Liang,¹Ronghan He,¹and Kun Wang¹

Academic Editor: Bing Niu

Received06 Jun 2020

Revised27 Aug 2020

Accepted11 Sept 2020

Published21 Sept 2020

Abstract

Osteosarcoma is a highly malignant bone cancer which primarily occurs in children and young adults. Increasing evidence indicates that long noncoding RNAs (lncRNAs), which function as competing endogenous RNAs (ceRNAs) that sponge microRNAs (miRNAs) and messenger RNAs (mRNAs), play a pivotal role in the pathogenesis and progression of cancers. The regulatory mechanisms of lncRNA-mediated ceRNAs in osteosarcoma have not been fully elucidated. In this study, we identified differentially expressed lncRNAs, miRNAs, and mRNAs in osteosarcoma based on RNA microarray profiles in the Gene Expression Omnibus (GEO) database. A ceRNA network was constructed utilizing bioinformatic tools. Kaplan-Meier survival analysis showed that lncR-C3orf35 and HMGB1 were associated with poor prognosis of osteosarcoma patients. Furthermore, results of Gene Set Enrichment Analysis (GSEA) suggested that lncR-C3orf35 may be involved in cellular invasion, the Toll-like receptor signaling pathway, and immune cell infiltration in the tumor microenvironment. Further analysis showed that patients with osteosarcoma metastasis expressed higher levels of lncR-C3orf35 and HMGB1 compared to metastasis-free patients. Moreover, the metastasis-free survival rate of the high lncR-C3orf35/HMGB1 expression group was significantly lower than that of the low expression group. The ESTIMATE algorithm was used to calculate the immune score and stromal scores for each sample. High lncR-C3orf35 and HMGB1 levels were correlated with low immune scores. ImmuCellAI analysis revealed that a low proportion of macrophage infiltration was associated with high lncR-C3orf35 and HMGB1 expression. The differential expression of lncR-C3orf35, miR-142-3p, and HMGB1 was further verified by quantitative real-time PCR. This study indicates that lncR-C3orf35 could be considered as a novel potential biomarker and therapeutic target of osteosarcoma, which may contribute to a better understanding of ceRNA regulatory mechanisms.

1. Introduction

Osteosarcoma, the most common bone malignancy, is a major cause of cancer-associated mortality in children and adolescents [1]. Osteosarcoma primarily affects the terminus of long bones, such as proximal tibias and distal femurs. Approximately 5 new osteosarcoma cases occur per million people under the age of 20 each year in America [2]. The incidence of osteosarcoma has increased by 0.3% per year over the last decades [3]. Although surgery and chemotherapeutic regimens have recently achieved notable progress in the treatment of osteosarcoma, the prognosis of osteosarcoma is still unsatisfactory. For patients exhibiting osteosarcoma metastasis at diagnosis, the 5-year survival rate is less than 30% [4]. The mechanisms underlying osteosarcoma carcinogenesis and progression have not been fully elucidated, and novel biomarkers such as therapeutic targets and predictors of prognosis are needed to be investigated.

Long noncoding RNAs (lncRNAs), of more than 200 nucleotides in length, are a major class of noncoding RNAs. Initially, lncRNAs were regarded as “junk genes” during transcription due to the absence of protein-coding capacity [5]. Recently, studies have shown that lncRNAs are involved in various biological processes, including gene transcription, chromatin modification, and epigenetic regulation [6]. Salmena et al. [7] determined that lncRNAs functioned as competing endogenous RNAs (ceRNAs). The ceRNAs act as sponges for target microRNAs (miRNAs) and regulate the miRNA-induced gene silencing. miRNAs are a class of short noncoding RNAs which contain 18 to 25 nucleotides. miRNAs can bind to base-complementary mRNAs, leading to inhibition of translation or degradation of mRNAs [8]. Moreover, aberrant expression of lncRNAs has been reported to be involved in the occurrence and development of osteosarcoma [9]. Increased lncRNAs including MALAT1 [10], TUG1 [11], and 91H [12] in osteosarcoma act as oncogenes to promote tumor proliferation and invasion. Conversely, decreased expressions of loc285194 [13] and MEG3 [14] have been reported to act as tumor suppressors in osteosarcoma. However, the majority of lncRNAs have not yet been functionally characterized in the pathogenesis of osteosarcoma and warrant further investigation.

In this study, osteosarcoma data were collected from the Gene Expression Omnibus (GEO) (http://www.ncbi.nlm.nih.gov/geo/), which includes 41 osteosarcoma samples and 8 normal controls. A ceRNA network with 14 lncRNAs, 9 miRNAs, and 58 mRNAs was obtained based on bioinformatic prediction analysis. Kaplan-Meier survival analysis of an independent dataset identified lncR-C3orf35 and HMGB1 as genes associated with survival. In addition, the potential biological activity and clinical features of lncR-C3orf35/HMGB1 were analyzed. The expression of the lncR-C3orf35/miR142-3p/HMGB1 axis was validated by quantitative real-time PCR (qRT-PCR) in osteosarcoma cells. Our study was aimed at analyzing the epigenetic mechanisms of osteosarcoma prognosis and providing support for lncRNAs as potential biomarkers.

2. Materials and Methods

2.1. Collection of RNA Expression Profiles

RNA expression profiles of osteosarcoma were collected from GEO. lncRNA and mRNA microarray datasets were deposited by Sadikovic et al. [15] and Fritsche-Guenther et al. [16] with accession numbers GSE12865 and GSE14593, respectively (Table 1). The GSE12865 dataset comprised 12 osteosarcoma tissue samples and 2 normal tissue samples, while the GSE14359 dataset included 10 osteosarcoma and 2 normal tissue samples. The GSE28423 dataset was deposited by Namløs et al. [17] and included miRNA expression profiles of 19 osteosarcoma cell lines and 4 normal osteoblast cell lines. The 3 datasets were used to screen differentially expressed RNAs in osteosarcoma. We downloaded the GSE21257 [18] dataset which includes 53 osteosarcoma patients with mRNA and lncRNA expression profiles and clinical features for further analysis.

2.2. Identification of Differentially Expressed RNAs

We directly downloaded RNA expression series-matrix profiles, which had been processed by the original authors. The names of lncRNAs, miRNAs, and mRNAs in the microarrays were annotated using the corresponding GEO platforms. Expression measurements of multiple probe sets mapping to the same genes were averaged, and the probe sets mapping to multiple genes were removed. The and logfold change were used as the cut-off criteria. Differentially expressed RNAs with statistical significance between osteosarcoma and normal controls in each dataset were identified using the limma R package. Differentially expressed (DE) RNAs in the GSE12865 and GSE14359 datasets were merged. The miRcode [19] lncRNA annotation was used to screen lncRNA genes. Hence, using this approach, we could identify DE mRNAs (DEmRNAs), lncRNAs (DElncRNAs), and miRNAs (DEmiRNAs).

2.3. Establishment of the ceRNA Regulatory Network

lncRNA-miRNA interactions for DElncRNAs and DEmiRNAs were predicted using miRcode [19]. miRNA-mRNA regulatory pairs among DEmiRNAs and DEmRNAs were obtained using miRTarBase [20], TargetScan [21], and miRDB [22]. Overlapping regulatory pairs were identified by Venn diagrams, which were then used to construct the lncRNA-miRNA-mRNA network based on the union of DElncRNA-DEmiRNA pairs and DEmiRNA-DEmRNA pairs. Cytoscape v3.7.1 [23] software was used to visualize the ceRNA work.

2.4. Survival Analysis

Fifty-three samples with RNA expression profiles and survival information in GSE21257 were used for the Kaplan-Meier (KM) survival analysis [24]. All lncRNAs and mRNAs in the ceRNA network were included. A value < 0.05 was set as the cut-off criterion. Then, the survival-associated lncRNA-mRNA pairs were selected and visualized with KM survival curves using the survival R package.

2.5. Regression Analysis of the lncR-C3orf35/HMGB1 Pair

Regression analysis of the expression of survival-associated lncRNA and mRNA was completed by Spearman’s correlation analysis and was visualized by the ggpubr R package. and were used as the selection criteria.

2.6. Gene Set Enrichment Analysis for lncR-C3orf35

To identify the key biological process and signaling pathway affected by lncRNA, we classified samples into high and low expression lncR-C3orf35 groups. Gene Ontology (GO) biological process and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were conducted using Gene Set Enrichment Analysis (GSEA) v4.0.2 (https://software.broadinstitute.org/gsea/downloads.jsp) to analyze the mRNA expression profile of GSE21257.

2.7. Clinical Feature Analysis for lncR-C3orf35/HMGB1

The GSE21257 dataset contains demographic and metastasis information for 53 osteosarcoma patients. We compared the lncR-C3orf35 and HMGB1 expression levels of the metastasis and nonmetastasis groups using unpaired -test. The samples were divided into the high expression group and low expression group based on median lncR-C3orf35 or HMGB1 levels. The high expression group consisted of 26 samples, and the low expression group contained 27 samples. Demographic data of patients were compared between the high and low expression groups by a two-way ANOVA test. Then, KM metastasis-free survival analysis was conducted using the survival R package. Huvos grades consequent to chemotherapy were recorded in the dataset. Huvos grades 1 and 2 were regarded as poor responses. Conversely, Huvos grades 3 and 4 were regarded as good responses. We compared lncR-C3orf35 and HMGB1 expression between the good and poor response groups.

2.8. Tumor Immune Infiltration Analysis

ESTIMATE (Estimation of STromal and Immune cells in MAlignant Tumor tissues using Expression data) is an algorithm used to evaluate the tumor composition by calculating the immune score and stromal score [25]. Using the ESTIMATE R package, the immune and stromal scores were obtained based on the RNA expression data of osteosarcoma samples in the GSE21257 dataset. Then, the scores were compared between the high and low lncR-C3orf35/HMGB1 expression groups. ImmuCellAI is a website tool that estimated the abundance of 24 types of immune cell infiltration in tumor tissue from a gene expression dataset [26]. Normalized expression data of GSE21257 was uploaded to the web portal (http://bioinfo.life.hust.edu.cn/web/ImmuCellAI/). Then, the estimated proportion of immune cell types can was obtained for each tumor sample. For each cell subset, differences in high lncR-C3orf35/HMGB1 and low lncR-C3orf35/HMGB1 expressions were compared.

2.9. Cell Culture and qRT-PCR Validation

The normal human osteoblast cell line hFOB 1.19 and osteosarcoma cell line SASJ-2 were purchased from ATCC (USA). The cells were cultured in Dulbecco’s Modified Eagle’s Medium (DMEM, Gibco, Invitrogen, USA) supplemented with 10% FBS (Gibco, Invitrogen, USA) and grown at 37°C in the presence of 5% CO₂ air atmosphere.

Total RNA was isolated from cells using TRIzol (Invitrogen) according to the manufacturer’s instructions. For the detection of lncR-C3orf35 and HMGB1, the total RNA of each sample was reverse transcribed to cDNA using the PrimeScript RT reagent kit (Takara, Tokyo, Japan). GAPDH served as the internal control. All-in-One™ miRNA real-time quantitative polymerase chain reaction (qRT-PCR) was employed to detect miR-142-3p expression levels using U6 snRNA as the internal control. qRT-PCR was performed with SYBR Premix Ex Taq II (TaKaRa, Tokyo, Japan) on the ABI 7500 thermocycler (Applied Biosystems; Thermo Fisher Scientific, USA). Relative expression was calculated by the 2^-ΔΔCt method, and each experiment was performed in triplicate. . The experiments were repeated three times. The sequences of real-time PCR primers are listed in Table 2.

2.10. Statistical Analysis

R v3.5.3 and GraphPad Prism v8.0 were used for the statistical analysis. Student’s -test was used to determine differences between two groups. Any value < 0.05 was considered statistically significant.

3. Results

3.1. Differentially Expressed lncRNA, mRNA, and miRNA

With the differentially expressed criteria of absolute and , a total of 824 DEmRNAs and 4 DElncRNAs were identified in the GSE143599 dataset (Figure 1). Meanwhile, 1794 DEmRNAs and 32 DElncRNAs were extracted in the GSE12865 dataset. With the union of the two sets, a total of 36 DElncRNAs and 2357 DEmRNAs were obtained. Of these, 21 lncRNAs and 1142 mRNAs were upregulated, while 15 lncRNAs and 1215 mRNAs were downregulated. Additionally, 205 DEmiRNAs were identified in the GSE28423 dataset, of which 101 miRNAs were upregulated and 104 miRNAs were downregulated.

(a)

(b)

(c)

(d)

Figure 1

Differentially expressed (DE) lncRNAs and mRNAs in osteosarcoma and normal samples from the GSE14359 and GSE1865 datasets. (a) Heat map of DElncRNAs in GSE12865. (b) Volcano plot of DEGs in GSE12865. (c) Heat map of DElncRNAs in GSE14359. (d) Volcano plot of DEGs in GSE14359. In the heat maps, orange blocks indicate upregulated lncRNAs. Blue blocks indicate downregulated lncRNAs. In the volcano plots, red and green points represent significantly upregulated and downregulated genes, respectively, with and . DEGs: differentially expressed genes.

3.2. Construction of ceRNA Network

lncRNAs and mRNAs targeted by miRNAs were screened based on the interactions among the DElncRNAs, DEmRNAs, and DEmiRNAs described above. Nine of 205 DEmiRNAs were predicted to interact with 14 of 36 DElncRNAs according to miRcode. Next, 58 mRNAs of 2357 DEmRNAs were screened that targeted 9 DEmiRNAs by searching TargetScan, miRTarBase, and miRDB. Finally, the ceRNA network consisting of 14 lncRNAs, 9 miRNAs, and 58 mRNAs was constructed by Cytoscape software and visualized using the networkD3 R package (Figure 2).

3.3. Survival Analysis

Overall survival is considered the ultimate detection standard for the prediction of patient prognosis; hence, screening of survival-associated lncRNAs and mRNAs is important. KM analysis was carried out to explore associations between survival of osteosarcoma patients and DElncRNAs/DEmRNAs included in the ceRNA network. High expressions of both lncR-C3orf35 and HMGB1 were significantly associated with poor overall survival of osteosarcoma patients (Figure 3).

(a)

(b)

(c)

(d)

3.4. Regression Analysis

Since the data of lncR-C3orf35 and HMGB1 expression was not normally distributed, we performed Spearman’s correlation analysis to compare lncR-C3orf35 and HMGB1. The results showed and , which indicated that HMGB1 expression was significantly positively correlated with lncR-C3orf35.

3.5. GSEA for lncR-C3orf35

Samples were divided into high and low lncR-C3orf35 expression groups, and analyses of GO biological processes and KEGG GSEA were conducted for the two groups. Forty GO terms were obtained which mainly included terms associated with leukocyte chemotaxis, lymphocyte migration, chemokine production, and response to interferon gamma. Furthermore, KEGG pathways include antigen processing and presentation, B cell receptor signaling, chemokine signaling, Nod-like receptor signaling, and Toll-like receptor (TLR) signaling, as well as cell adhesion molecules (cams), cell cycle, and mismatch repair processes (Figure 4). These biological processes and pathways were associated with cell proliferation, invasion, and immune cell infiltration in the tumor microenvironment.

3.6. Clinical Features with lncR-C3orf35/HMGB1

For further analysis, we investigated correlations between lncR-C3orf35/HMGB1 and clinical features including metastasis and chemotherapy response. The demographic data of patients are displayed in Table 3. There was no difference in the sex and age distribution between the high and low expression groups. Thirty-four of 53 osteosarcoma patients in the GSE21257 dataset presented metastasis at diagnosis or during the follow-up period. lncR-C3orf35 and HMGB1 expression in the metastasis group was significantly higher than that in the nonmetastasis group (lncR-C3orf35: ; HMGB1: ) (Figure 5). Next, samples were classified into the high and low expression groups based on the median expression of either lncR-C3orf35 or HMGB1. KM analysis for metastasis-free survival between high and low lncR-C3orf35/HMGB1 expression groups was performed. The metastasis-free survival rate of the high lncR-C3orf35/HMGB1 expression group was lower than that of the low lncR-C3orf35/HMGB1 expression group (lncR-C3orf35: ; HMGB1: ), which indicated that lncR-C3orf35/HMGB1 was correlated with osteosarcoma metastasis. The GSE21257 dataset recorded the Huvos grades of 47 osteosarcoma patients. There are 18 samples in the good response group (Huvos grades 3 and 4) and 29 samples in the poor response group (Huvos grades 1 and 2). lncR-C3orf35 expression in the good response group () was similar to that in the poor response group (, ). Likewise, there was no significant difference in HMGB1 expression between the good () and poor (, ) response groups.

(a)

(b)

(c)

(d)

(e)

(f)

3.7. Immune Infiltration Analysis

The ESTIMATE algorithm assessing gene signatures was utilized to calculate a stromal score and immune score which represented the proportion of stromal cells and infiltrating immune cells, separately. The stromal score and immune score of each sample in the GSE21257 dataset were obtained using the ESTIMATE R package. The high lncR-C3orf35 expression group presented a significantly lower immune score ( vs. , ) than the low lncR-C3orf35 expression group (Figure 6). Similarly, the high HMGB1 expression group presented a lower immune score ( vs. , ) than the low lncR-C3orf35 expression group, while differences in stromal scores between the high and low lncR-C3orf35/HMGB1 expression groups were not significant. Hence, lncR-C3orf35/HMGB1 expression was associated with immune cell infiltration in osteosarcoma. ImmuCellAI was used to evaluate the proportion of 24 immune cell subtypes in tumor tissue basing on gene expression data. Twenty-four subtypes of immune cell proportion of each sample were calculated. The macrophage proportion was lower in the high lncR-C3orf35 expression group () than in the low lncR-C3orf35 expression group (, ) and was also lower in the high HMGB1 expression group () than in the low HMGB1 expression group (, ) (Figure 7). The results indicated that macrophage infiltration in osteosarcoma tissue was associated with lncR-C3orf35/HMGB1 expression.

(a)

(b)

(c)

(d)

(a)

(b)

3.8. Validation by qRT-PCR

To validate the bioinformatic analysis results, we utilized qRT-PCR to detect the expression of lncR-C3orf35, miR-142-3p, and HMGB1 in cell lines. lncR-C3orf35 () and HMGB1 () were overexpressed in the osteosarcoma cell line SASJ-2 compared to the human osteoblast cell line hFOB 1.19. In addition, miR-142-3p was pressed in low levels in SASJ-2 cells () (Figure 8). The results were consistent with the results described above.

(a)

(b)

(c)

4. Discussion

Clinicians face a great challenge when confirming early diagnosis and attempting to predict the prognosis of osteosarcoma patients. Discovering predictive biomarkers and exploring the molecular mechanism of osteosarcoma may provide novel insights for the diagnosis and treatment of osteosarcoma patients. Growing evidence has indicated that dysregulated expression of lncRNAs is associated with carcinogenesis [27]. According to the ceRNA theory, lncRNAs could act as sponges by the binding of base pairs with miRNAs, through which they inhibit miRNA-induced gene silencing. Construction of a ceRNA network contributes to reveal the occurrence and progression of cancer. Systematic analyses of ceRNA networks have been reported in lung cancer [28], gastric cancer [29], and colon cancer [30], but there have been few reports in osteosarcoma.

In this study, we identified 36 lncRNAs, 2357 mRNAs, and 205 miRNAs that were differentially expressed comparing osteosarcoma and normal controls based on RNA microarray data from the GEO database. We constructed a ceRNA regulatory network by utilizing several predictive bioinformatic tools. Through KM analysis, the lncR-C3orf35/HMGB1 pair was found to be associated with the overall survival of osteosarcoma patients. GSEA analysis was performed and revealed that lncR-C3orf35 was involved in cell adhesion, regulation of the cell cycle, and immune cell infiltration. Next, we found that lncR-C3orf35 and HMGB1 expression was correlated with osteosarcoma metastasis. Patients with high expression of lncR-C3orf35/HMGB1 presented lower metastasis-free survival. Using the ESTIMATE algorithm, immune scores of the high lncR-C3orf35/HMGB1 expression group were significantly lower than those of the low expression group, which suggests that immune cell infiltration differed between these two groups. We calculated the proportion of immune cell subtypes for each sample utilizing the ImmuCellAI online tool and found that the macrophage proportion in the high lncR-C3orf35/HMGB1 expression group was significantly lower than that in the low expression group. Moreover, the expression of lncR-C3orf35, miR-142-3p, and HMGB1 was validated in cell lines by qRT-PCR.

LncRbase (http://bicresources.jcbose.ac.in/zhumur/lncrbase/index.html) showed that lncR-C3orf35 is located in chr3:3744065-37476988 and has a length of 574 bases. There have been no reports on the association between lncR-C3orf35 and cancer. Our study showed that lncR-C3orf35 was overexpressed in osteosarcoma tissue and cells. In addition, upregulation of lncR-C3orf35 was found to be associated with poor overall survival and metastasis-free survival, which indicated that lncR-C3orf35 was a potential biomarker of osteosarcoma development.

The GSEA KEGG analysis identified an enriched TLR pathway, which could bind HMGB1 and activate a series of inflammatory responses [31]. Our results showed that lncR-C3orf35/HMGB1 expression was correlated with a low immune score and low macrophage infiltration. Previous studies have reported that a high macrophage proportion in osteosarcoma tissue was associated with better prognosis and inhibition of metastasis [32, 33]. Macrophages are a crucial member of the innate immune system and are also the most abundant immune cell in the tumor microenvironment [34]. Tumor-associated macrophages differentiate into M1-type and M2-type macrophages under various cytokine stimulation. M1 macrophages mediate inflammation and kill tumor cells, while M2 macrophages promote tumor growth and induce immune suppression. The balance of M1/M2 is closely associated with tumor metastasis [35]. HMGB1 has been reported to induce tumor-associated macrophage polarizing to M2 macrophages through the NF-κB signaling pathway [36]. Hence, the lncR-C3orf35/HMGB1 axis may affect the polarization of macrophages, which promotes osteosarcoma metastasis.

HMGB1 is a highly conserved protein with 215 amino acids and 3 structural domains. It was initially regarded as a nuclear protein regulating transcription. HMGB1 has been shown to bind TLRs and RAGEs, through which it induces the secretion of various proinflammatory cytokines and triggers a series of inflammatory responses [37]. Under HMGB1 stimulation, angiogenesis and immune inhibition promote tumor metastasis [38]. Conversely, high levels of HMGB1 in osteosarcoma may induce autophagy of osteosarcoma cells, which contributes to resistance to cisplatin, doxorubicin, and methotrexate treatment [39]. Nonetheless, our study did not identify any correlation between HMGB1 expression and the chemotherapy response due to the small sample size. Extracellular HMGB1 protein results in higher osteosarcoma cell proliferation, migration, and osteogenic differentiation [40]. Our study showed that high HMGB1 expression was associated with metastasis and immune cell infiltration of osteosarcoma, which was consistent with previous studies. In addition, we propose that HMGB1 may mediate the influence of lncR-C3orf35 on osteosarcoma prognosis through a ceRNA mechanism.

Xu et al. [41] reported that miR-142-3p was downregulated in osteosarcoma tissue and cell lines. Several target mRNAs of miR-142-3p have been proposed including HMGA1, HMGB1, and Rac1 [42]. In hepatocellular cancer, miR-142-3p was shown to bind the 3-UTR region of HMGB1 and resulted in the inhibition of proliferation and invasion of tumor cells [43]. Moreover, miR-142-3p has the potential to predict prognosis of colorectal cancer [44] and renal carcinoma [45] patients. In our study, miR-142-3p could interact with lncR-C3orf35 and HMGB1 based on predictive bioinformatic tools. We further confirmed differences in gene expression by microarray data and qRT-PCR, which supported the expression of the lncR-C3orf35/miR142-3p/HMGB1 axis in osteosarcoma. Further studies are needed to explore the biological functions of the lncR-C3orf35/miR142-3p/HMGB1 axis in osteosarcoma.

In conclusion, we identified dysregulated lncRNAs, miRNAs, and mRNAs in osteosarcoma and constructed a ceRNA regulatory network, which contribute to a better understanding of osteosarcoma progression. Several genes in the network were shown to be dysregulated and to function as oncogenes or tumor suppressors in osteosarcoma, such as HMGB1 and miR-142-3p. lncR-C3orf35 and HMGB1 were associated with poor prognosis of osteosarcoma patients. We propose that lncR-C3orf35 may upregulate HMGB1 by sponging miR-142-3p. Moreover, lncR-C3orf35 and HMGB1 may affect immune cell infiltration and tumor metastasis. Thus, our study supports lncR-C3orf35 as a new biomarker of osteosarcoma carcinogenesis and as a potential therapeutic target.

Data Availability

All the datasets used for this study can be found in GEO: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE14359, https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE12865, https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE28423, and https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE21257.

Conflicts of Interest

The authors declare that they have no competing interests.

Authors’ Contributions

Yi Shi was responsible for the conceptualization, methodology, software, and writing of original draft. Jianhua Ren was responsible for the formal analysis, investigation, data curation, and project administration. Ze Zhuang was responsible for the validation and data visualization. Wenhui Zhang was responsible for the formal analysis and data curation. Zhe Wang was responsible for the investigation and resources. Yuangao Liu was responsible for the software. Jinze Li was responsible for the software and data visualization. Tangzhao Liang was responsible for the project administration. Ronghan He was responsible for the writing (review and editing) and supervision. Kun Wang was responsible for the project administration and funding acquisition. All authors read and approved the submitted version. Yi Shi and Jianhua Ren contributed equally to this work.

Acknowledgments

This work was supported by grants from the National Natural Science Foundation of China (81772368 and 81802184). This work was also funded by the Natural Science Foundation of Guangdong Province (2017A030313496 and 2017A030310226).

References

B. A. Lindsey, J. E. Markel, and E. S. Kleinerman, “Osteosarcoma overview,” Rheumatology and Therapy, vol. 4, no. 1, pp. 25–43, 2017.
View at: Publisher Site | Google Scholar
G. Ottaviani and N. Jaffe, “The epidemiology of osteosarcoma,” Cancer Treatment and Research, vol. 152, pp. 3–13, 2009.
View at: Publisher Site | Google Scholar
K. Song, J. Song, K. Lin et al., “Survival analysis of patients with metastatic osteosarcoma: a surveillance, epidemiology, and end results population-based study,” International Orthopaedics, vol. 43, no. 8, pp. 1983–1991, 2019.
View at: Publisher Site | Google Scholar
A. Misaghi, A. Goldin, M. Awad, and A. A. Kulidjian, “Osteosarcoma: a comprehensive review,” SICOT-J, vol. 4, p. 12, 2018.
View at: Publisher Site | Google Scholar
J. X. Yang, R. H. Rastetter, and D. Wilhelm, “Non-coding RNAs: an introduction,” Advances in Experimental Medicine and Biology, vol. 886, pp. 13–32, 2016.
View at: Publisher Site | Google Scholar
X. Li, Z. Wu, X. Fu, and W. Han, “Long noncoding RNAs: insights from biological features and functions to diseases,” Medicinal Research Reviews, vol. 33, no. 3, pp. 517–553, 2013.
View at: Publisher Site | Google Scholar
L. Salmena, L. Poliseno, Y. Tay, L. Kats, and P. P. Pandolfi, “A ceRNA hypothesis: the Rosetta Stone of a hidden RNA language?” Cell, vol. 146, no. 3, pp. 353–358, 2011.
View at: Publisher Site | Google Scholar
J. Li, S. Tan, R. Kooger, C. Zhang, and Y. Zhang, “MicroRNAs as novel biological targets for detection and regulation,” Chemical Society Reviews, vol. 43, no. 2, pp. 506–517, 2014.
View at: Publisher Site | Google Scholar
Z. Yang, X. Li, Y. Yang, Z. He, X. Qu, and Y. Zhang, “Long noncoding RNAs in the progression, metastasis, and prognosis of osteosarcoma,” Cell Death & Disease, vol. 7, no. 9, article e2389, 2016.
View at: Publisher Site | Google Scholar
K. T. Gao and D. Lian, “Long non-coding RNA MALAT1 is an independent prognostic factor of osteosarcoma,” European Review for Medical and Pharmacological Sciences, vol. 20, no. 17, pp. 3561–3565, 2016.
View at: Google Scholar
Y. Li, T. Zhang, Y. Zhang, X. Zhao, and W. Wang, “Targeting the FOXM1-regulated long noncoding RNA TUG1 in osteosarcoma,” Cancer Science, vol. 109, no. 10, pp. 3093–3104, 2018.
View at: Publisher Site | Google Scholar
J. Zhao and S. T. Ma, “Downregulation of lncRNA H19 inhibits migration and invasion of human osteosarcoma through the NF-κB pathway,” Molecular Medicine Reports, vol. 17, no. 5, pp. 7388–7394, 2018.
View at: Publisher Site | Google Scholar
I. Pasic, A. Shlien, A. D. Durbin et al., “Recurrent focal copy-number changes and loss of heterozygosity implicate two noncoding RNAs and one tumor suppressor gene at chromosome 3q13.31 in osteosarcoma,” Cancer Research, vol. 70, no. 1, pp. 160–171, 2010.
View at: Publisher Site | Google Scholar
B. Shen, N. Zhou, T. Hu, W. Zhao, D. Wu, and S. Wang, “LncRNA MEG3 negatively modified osteosarcoma development through regulation of miR-361-5p and FoxM1,” Journal of Cellular Physiology, vol. 234, no. 8, pp. 13464–13480, 2019.
View at: Publisher Site | Google Scholar
B. Sadikovic, M. Yoshimoto, S. Chilton-MacNeill, P. Thorner, J. A. Squire, and M. Zielenska, “Identification of interactive networks of gene expression associated with osteosarcoma oncogenesis by integrated molecular profiling,” Human Molecular Genetics, vol. 18, no. 11, pp. 1962–1975, 2009.
View at: Publisher Site | Google Scholar
R. Fritsche-Guenther, A. Noske, U. Ungethüm et al., “De novo expression of EphA2 in osteosarcoma modulates activation of the mitogenic signalling pathway,” Histopathology, vol. 57, no. 6, pp. 836–850, 2010.
View at: Publisher Site | Google Scholar
H. M. Namløs, L. A. Meza-Zepeda, T. Barøy et al., “Modulation of the osteosarcoma expression phenotype by microRNAs,” PLoS One, vol. 7, no. 10, article e48086, 2012.
View at: Publisher Site | Google Scholar
E. P. Buddingh, M. L. Kuijjer, R. A. Duim et al., “Tumor-infiltrating macrophages are associated with metastasis suppression in high-grade osteosarcoma: a rationale for treatment with macrophage activating agents,” Clinical Cancer Research: An Official Journal of the American Association for Cancer Research, vol. 17, no. 8, pp. 2110–2119, 2011.
View at: Publisher Site | Google Scholar
A. Jeggari, D. S. Marks, and E. Larsson, “miRcode: a map of putative microRNA target sites in the long non-coding transcriptome,” Bioinformatics, vol. 28, no. 15, pp. 2062-2063, 2012.
View at: Publisher Site | Google Scholar
C. H. Chou, N. W. Chang, S. Shrestha et al., “miRTarBase 2016: updates to the experimentally validated miRNA-target interactions database,” Nucleic Acids Research, vol. 44, no. D1, pp. D239–D247, 2016.
View at: Publisher Site | Google Scholar
D. M. Garcia, D. Baek, C. Shin, G. W. Bell, A. Grimson, and D. P. Bartel, “Weak seed-pairing stability and high target-site abundance decrease the proficiency of lsy-6 and other microRNAs,” Nature Structural & Molecular Biology, vol. 18, no. 10, pp. 1139–1146, 2011.
View at: Publisher Site | Google Scholar
W. Liu and X. Wang, “Prediction of functional microRNA targets by integrative modeling of microRNA binding and target expression data,” Genome Biology, vol. 20, no. 1, p. 18, 2019.
View at: Publisher Site | Google Scholar
P. Shannon, A. Markiel, O. Ozier et al., “Cytoscape: a software environment for integrated models of biomolecular interaction networks,” Genome Research, vol. 13, no. 11, pp. 2498–2504, 2003.
View at: Publisher Site | Google Scholar
V. S. Stel, F. W. Dekker, G. Tripepi, C. Zoccali, and K. J. Jager, “Survival analysis I: the Kaplan-Meier method,” Nephron. Clinical Practice, vol. 119, no. 1, pp. c83–c88, 2011.
View at: Publisher Site | Google Scholar
K. Yoshihara, M. Shahmoradgoli, E. Martínez et al., “Inferring tumour purity and stromal and immune cell admixture from expression data,” Nature Communications, vol. 4, no. 1, p. 2612, 2013.
View at: Publisher Site | Google Scholar
Y. R. Miao, Q. Zhang, Q. Lei et al., “ImmuCellAI: a unique method for comprehensive T-cell subsets abundance prediction and its application in cancer immunotherapy,” Advanced science, vol. 7, no. 7, article 1902880, 2020.
View at: Publisher Site | Google Scholar
A. Fatica and I. Bozzoni, “Long non-coding RNAs: new players in cell differentiation and development,” Nature Reviews Genetics, vol. 15, no. 1, pp. 7–21, 2014.
View at: Publisher Site | Google Scholar
R. Li, Y. E. Yang, J. Jin et al., “Identification of lncRNA biomarkers in lung squamous cell carcinoma using comprehensive analysis of lncRNA mediated ceRNA network,” Artificial Cells, Nanomedicine, and Biotechnology, vol. 47, no. 1, pp. 3246–3258, 2019.
View at: Publisher Site | Google Scholar
W. Gu, J. H. Ren, X. Zheng, X. Y. Hu, and M. J. Hu, “Comprehensive analysis of expression profiles of long non-coding RNAs with associated ceRNA network involved in gastric cancer progression,” Molecular Medicine Reports, vol. 20, no. 3, pp. 2209–2218, 2019.
View at: Publisher Site | Google Scholar
M. Wu, W. Li, F. Huang et al., “Comprehensive analysis of the expression profiles of long non-coding RNAs with associated ceRNA network involved in the colon cancer staging and progression,” Scientific Reports, vol. 9, no. 1, article 16910, 2019.
View at: Publisher Site | Google Scholar
Y. Zhang, R. Karki, and O. J. Igwe, “Toll-like receptor 4 signaling: a common pathway for interactions between prooxidants and extracellular disulfide high mobility group box 1 (HMGB1) protein-coupled activation,” Biochemical Pharmacology, vol. 98, no. 1, pp. 132–143, 2015.
View at: Publisher Site | Google Scholar
C. Dumars, J. M. Ngyuen, A. Gaultier et al., “Dysregulation of macrophage polarization is associated with the metastatic process in osteosarcoma,” Oncotarget, vol. 7, no. 48, pp. 78343–78354, 2016.
View at: Publisher Site | Google Scholar
A. Gomez-Brouchet, C. Illac, J. Gilhodes et al., “CD163-positive tumor-associated macrophages and CD8-positive cytotoxic lymphocytes are powerful diagnostic markers for the therapeutic stratification of osteosarcoma patients: an immunohistochemical analysis of the biopsies from the French OS2006 phase 3 trial,” OncoImmunology, vol. 6, no. 9, article e1331193, 2017.
View at: Publisher Site | Google Scholar
M. F. Heymann, F. Lézot, and D. Heymann, “The contribution of immune infiltrates and the local microenvironment in the pathogenesis of osteosarcoma,” Cellular Immunology, vol. 343, article 103711, 2019.
View at: Publisher Site | Google Scholar
A. Sica, M. Erreni, P. Allavena, and C. Porta, “Macrophage polarization in pathology,” Cellular and Molecular Life Sciences: CMLS, vol. 72, no. 21, pp. 4111–4126, 2015.
View at: Publisher Site | Google Scholar
A. Rojas, P. Araya, J. Romero et al., “Skewed signaling through the receptor for advanced glycation end-products alters the proinflammatory profile of tumor-associated macrophages,” Cancer Microenvironment, vol. 11, no. 2-3, pp. 97–105, 2018.
View at: Publisher Site | Google Scholar
G. P. Sims, D. C. Rowe, S. T. Rietdijk, R. Herbst, and A. J. Coyle, “HMGB1 and RAGE in inflammation and cancer,” Annual Review of Immunology, vol. 28, no. 1, pp. 367–388, 2010.
View at: Publisher Site | Google Scholar
X. Wang, L. Xiang, H. Li et al., “The role of HMGB1 signaling pathway in the development and progression of hepatocellular carcinoma: a review,” International Journal of Molecular Sciences, vol. 16, no. 9, pp. 22527–22540, 2015.
View at: Publisher Site | Google Scholar
J. Huang, K. Liu, Y. Yu et al., “Targeting HMGB1-mediated autophagy as a novel therapeutic strategy for osteosarcoma,” Autophagy, vol. 8, no. 2, pp. 275–277, 2012.
View at: Publisher Site | Google Scholar
S. Martinotti, M. Patrone, M. Manfredi et al., “HMGB1 osteo-modulatory action on osteosarcoma SaOS-2 cell line: an integrated study from biochemical and -omics approaches,” Journal of Cellular Biochemistry, vol. 117, no. 11, pp. 2559–2569, 2016.
View at: Publisher Site | Google Scholar
G. Xu, J. Wang, Y. Jia, F. Shen, W. Han, and Y. Kang, “MiR-142-3p functions as a potential tumor suppressor in human osteosarcoma by targeting HMGA1,” Cellular Physiology and Biochemistry: International Journal of Experimental Cellular Physiology, Biochemistry, and Pharmacology, vol. 33, no. 5, pp. 1329–1339, 2014.
View at: Publisher Site | Google Scholar
P. Shabani, S. Izadpanah, A. Aghebati-Maleki et al., “Role of miR-142 in the pathogenesis of osteosarcoma and its potential as therapeutic approach,” Journal of Cellular Biochemistry, vol. 120, no. 4, pp. 4783–4793, 2019.
View at: Publisher Site | Google Scholar
Y. Fu, L. Q. Sun, Y. Huang et al., “miR-142-3p inhibits the metastasis of hepatocellular carcinoma cells by regulating HMGB1 gene expression,” Current Molecular Medicine, vol. 18, no. 3, pp. 135–141, 2018.
View at: Publisher Site | Google Scholar
W. Gao, D. Pang, and S. Yu, “Serum level of miR-142-3p predicts prognostic outcome for colorectal cancer following curative resection,” The Journal of International Medical Research, vol. 47, no. 5, pp. 2116–2125, 2019.
View at: Publisher Site | Google Scholar
X. Peng, X. Pan, K. Liu et al., “miR-142-3p as a novel biomarker for predicting poor prognosis in renal cell carcinoma patients after surgery,” The International Journal of Biological Markers, vol. 34, no. 3, pp. 302–308, 2019.
View at: Publisher Site | Google Scholar

Copyright

Copyright © 2020 Yi Shi et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

PDF Download Citation

Download other formats

Order printed copies

Views

293

Downloads

1126

Citations