A Novel Predictive Model Associated with Osteosarcoma Metastasis

Introduction

OS is the most common primary malignant bone tumor¹ with a bimodal age distribution. The first peak occurs during adolescence, while the second peak occurs later in life.^2–6 With the advancements of treatments, surgery and intensive multiagent chemotherapy became the mainstream strategies.^7–9 The high incidence of metastasis in OS, particularly lung metastasis, results in a poor prognosis.^10–12 The 5-year overall survival rate for non-metastatic OS is about 75%; while it is only 20% in metastatic patients.^13,14 Thus, developing biomarkers for early detection and treatment is essential.

LncRNA is a non-coding RNA ranging length from 200 nucleotides (nt) to more than 100 kilobases (kb).¹⁵ LncRNA modulates gene expression on both transcriptional and translational levels, including chromosomal rearrangement, histone modification, alternative splicing formation, and RNA stabilization.¹⁶ LncRNA has been found to play a role in tumor differentiation, invasion, and metastasis.^17–21 LncRNA is implicated in the onset and progression of cancer metastasis. MALAT1, for example, primarily regulates lung cancer metastasis through epithelial-mesenchymal transition (EMT).^22,23 It has also been shown to promote metastasis in bladder cancer by activating Wnt signaling.^24,25 HOTAIR promotes pancreatic cancer metastasis by suppressing interferon-related genes.²⁶ MEG3 plays a decisive role in meningioma and gastric cancer by regulating autophagy and DNA repair.^27,28 However, the role of lncRNA in OS remains unclear, and the lncRNA-based prediction model is currently lacking.

Advances in high-throughput technologies allow researchers to identify biomarkers that can be used to construct prediction models for personalized cancer treatment. In this study, we presented a new risk-score model for predicting the metastasis of OS patients based on lncRNAs. In addition, we would detect the critical biological process, immune-related cells, and immune checkpoint between high and low risk group.

Materials and Methods

Target Data and Processing

The clinical information and expression profiles, including mRNA and lncRNA expression profiles, were downloaded from the TARGET datasets (https://ocg.cancer.gov/programs/target). All datasets were subjected to standardized data preprocessing. In this study, a total of 22 metastasis and 66 primary OS tissue-related data were analyzed. Data pre-processing was performed using RStudio and included background correction, normalization, and batch effect adjustment.

Feature Selection and Construction of Prognostic Model

The Limma R package was used to compare differential gene expression of lncRNA and mRNA in 22 metastasis and 66 primary OS patients. P-value<0.05 and |log 2 (fold change [FC]) | >1 were set as the threshold criteria. Differentially expressed lncRNA in this study, were screened by univariate cox regression, least absolute shrinkage, and selection operator (LASSO) regression, and multivariate cox regression analysis in survival and glmnet R packages. The lncRNA-based prognosis risk score was established using a linear formula with the expression level multiplied regression model (β). The median risk score was used to categorize patients into the high-risk group and low-risk group. The survival and survival ROC R packages were used to plot Kaplan–Meier (KM) survival and receiver operating characteristic (ROC) curves.

Scale-Free Co-Expression Network and Enrichment Analysis

Prognosis-related lncRNAs and metastasis-related mRNA were used to reconstruct the scale-free co-expression network. A total of 24 metastasis-related mRNAs associated with 7 metastasis-related lncRNAs (correlation coefficient R>0.3, p<0.001) were selected. Cytoscape software was used for network visualization. Using DAVID, we investigated the interactions between 24 metastasis-related mRNA and transcription factors (TF). GO pathway enrichment analyses were performed using the R package clusterProfiler. P<0.05 was considered as used as the cut-off for statistical significance.

Gene Set Enrichment Analysis and Gene Set Variation Analysis

GSEA2.2.4 software was used to perform gene set enrichment analysis (GSEA) and generate an enrichment score. The gene sets used for the enrichment analysis were downloaded from the Molecular Signatures Database (MsigDB). After 1000 permutations, gene sets with a false discovery rate (FDR) less than 0.05 were reported to be significantly enriched. Gene Set Variation Analysis (GSVA) is a method for estimating the variation of gene set enrichment using samples from an expression data set. GSVA was performed using the Hallmark gene set in the GSVA package (from R Project 3.5.1) of R software with default parameters.

Infiltrating Immune Cells

The CIBERSORT algorithm (http://cibersort.stanford.edu), a deconvolution algorithm based on gene expression profiles, was used to evaluate the 22 different types of tumor-infiltrating immune cells in the metastatic and primary OS group. CIBERSORT can assess the composition of tumor-infiltrating immune cells, by using P-value as a criterion for measuring confidence in findings. The relationship between 7 prognosis-related lncRNAs and tumor-infiltrating immune cells was also investigated.

Cell Lines and qPCR

Highly metastasis human OS cell lines (MNNG) and weakly metastasis human OS cell lines (U2S2, SAOS-2, SJSA-1, and HOS) were purchased from iCell (Shanghai, China). Trizol reagent (Invitrogen Life Technologies) was used to isolate total RNA. Approximately 2 μg RNA was reverse-transcribed into complementary DNA following the manufacturer’s instructions and normalized using glyceraldehyde 3-phosphate dehydrogenase mRNA levels. Reverse-transcription quantitative PCR assay with PrimeScript RT Master Mix (Takara) was performed on a QuantStudio6 Flex real-time PCR system (Applied Biosystems). The following were the PCR amplification conditions: 2 min at 95°C for one cycle, followed by denaturation for 15 sec at 95°C, and extension for 60 sec at 60°C for 38 cycles. Table S1 presents the primers used.

Result

Differential Gene and lncRNA Screening

This study included 22 metastatic and 66 primary OS samples. The Limma R package was used to select differentially expressed genes and lncRNA. As a result, 135 differentially expressed lncRNAs were identified, including 79 upregulated and 56 downregulated lncRNAs (Figure 1A). In addition, there were 171 differentially expressed genes (Figure 1B).

Figure 1 Differential expression profiles. (A) Differential lncRNA between metastatic and non-metastatic OS; (B) Differential mRNA between metastatic and non-metastatic OS.

Construction of Prognosis Assessment Model

The findings from the univariate Cox regression analysis of 135 lncRNAs were used in LASSO regression to identify potential markers. Multivariate cox regression analysis identified seven prognosis-related lncRNAs (Figure 2A–C), including AL512422.1, AL008718.3, C5orf66-AS1, AL360182.2, CEBPA-DT, AC006033.2, and AL357507.1. All of these lncRNAs were used to establish the prognosis model, based on the following formula, y= 2.603*AL512422.1+3.735*AL008718.3+18.047*C5orf66-AS1+1.292*AL360182.2+0.475*CEBPA-DT+0.038*AC006033.2+1.812*AL357507.1. The patients were divided into high-risk and low-risk groups based on the median risk score calculated using the above formula. Figure 2D depicts the expression of seven prognosis-related lncRNAs in high and low-risk groups. Moreover, the model had higher prediction ability, with all area under the curve (AUC) values >0.8 (1 year: AUC = 0.92, 95% Cl = 0.83–0.99; 3 years: AUC = 0.87, 95% Cl = 0.79–0.96; and 5 years: AUC = 0.86, 95% Cl = 0.76–0.96) (Figure 2D). According to the KM survival curves (Figure 2D), the high-risk group had significantly poorer overall survival time than the low-risk group (p<0.0001, HR=1.09, 95% CI = 1.06–1.13). In addition, the KM survival curves indicated that five lncRNAs, including AL512422.1, AL008718.3, C5orf66-AS1, AL360182.2, and AL357507.1, were associated with poor overall survival. CEBPA-DT and AC006033.2, on the other hand, were associated with improved survival (Figure 3A–G).

Figure 2 Identification of the prognosis-related lncRNA. (A and B) LASSO coefficient profiles; (C) multivariate cox analysis of seven model lncRNAs; (D) the risk scores, ROC curve, and Kaplan–Meier curve of patients.

Figure 3 Kaplan–Meier curve of seven prognosis-related lncRNAs. (A) AL512422.1. (B) AL008718.3. (C) AC006033.2. (D) AL357507.1. (E) AL360182.2. (F) C5orf66-AS1. (G) CEBPA-DT.

Co‑Expression Network and Enrichment Analysis

To explore the mRNAs associated with the model-related lncRNAs and their potential TF, we integrated the lncRNA and gene data to construct a co-expression network. The metastasis-related genes that were closely associated with the seven prognosis-related lncRNAs were selected (R>0.3, p<0.001). A total of 6 prognosis-related lncRNAs and 24 metastasis-related genes were used to construct the co-expression network (Figure 4A). AC006033.2 and AL512422.1 were the most important nodes. To further understand the functions of these genes, we assessed the interactions between 24 metastasis-related genes and 13 TFs (Figure 4B, Table 1). GO analysis was then carried out to determine the functional enrichment of the nodes in this co-expression network (Figure 4C). Several biological processes (BP) and cellular component (CC) pathways were found to be significantly enriched. According to BP, the nodes were mainly enriched in “catecholamine metabolic process”, “lymphocyte proliferation”, “mononuclear cell proliferation”, “lymphocyte differentiation” and “leukocyte proliferation”, while for CC they were enriched in, “integral component of synaptic membrane”, “intrinsic component of synaptic membrane”, and “glutamatergic synapse”.

Table 1 Transcription Factors Associated with 24 Intersecting mRNAs

Figure 4 Co-expression network and enrichment analysis. (A) Co-expression network of lncRNA and mRNA; (B) Co-expression network of mRNA and TF; (C) GO analysis of the genes in co-expression network.

GSVA Analysis of Seven Prognosis-Related lncRNAs

We performed GSVA to figure out the dynamics of biological processes and pathways for Hallmark gene sets based on seven prognosis-related lncRNAs. As Figure 5 showed, “G2M_CHECKPOINT”, “DNA_REPAIR” and “PI3K_AKT_MTOR_SIGNALING” were remarkably activated in the high AL512422.1 group, while the low AL512422.1 group was enriched in “WNT_BETA_CATENIN_SIGNALING” and “KRAS_SIGNALING_DN”. Upregulation of AL008718.3 activated “WNT_BETA_CATENIN_SIGNALING” and “TGF_BETA_SIGNALING”, but downregulation of AL008718.3 led to “HEDGEHOG_SIGNALING” and “PI3K_AKT_MTOR_SIGNALING”. In addition, AC006033.2 high level was significantly enriched for “HYPOXIA”, “APOPTOSIS” and “NOTCH_SIGNALING”, whereas low level triggered “DNA_REPAIR” and “G2M_CHECKPOINT” process. Notably, “GLYCOLYSIS” and “NOTCH_SIGNALING” were all enriched in AL357507.1 and AL360182.2 low groups. For the C5orf66-AS1 low group, “TGF_BETA_SIGNALING”, “HYPOXIA” and “NOTCH_SIGNALING” were enriched. The C5orf66-AS1 high group is related to “WNT_BETA_CATENIN_SIGNALING”. At last, we observed enrichment of “G2M_CHECKPOINT”, “DNA_REPAIR” and “PI3K_AKT_MTOR_SIGNALING” in the CEBPA-DT high group and “WNT_BETA_CATENIN_SIGNALING” in the low group. The top high-frequency enriched pathways were “NOTCH_SIGNALING”, “HYPOXIA”, “WNT_BETA_CATENIN_SIGNALING” and “HEDGEHOG_SIGNALING”.

Figure 5 GSVA enrichment analysis of seven prognosis-related lncRNAs. (A) AL512422.1. (B) AL008718.3. (C) AC006033.2. (D) AL357507.1. (E) AL360182.2. (F) C5orf66-AS1. (G) CEBPA-DT.

Correlation Between Prognosis-Related lncRNAs and Infiltrating Immunocyte Fractions

To investigate the differences between the primary and metastatic OS, 22 infiltrating immunocyte fractions were compared. As shown in Figure 6A, the difference between the two groups was insignificant. Then, the relationship between seven prognosis-related lncRNAs and 22 infiltrating immune cells was analyzed (Figure 6B). AC006033.2 was found to have the strongest relationship with resting dendritic cells, macrophages M0, macrophages M2, activated mast cells, activated neutrophils, and memory T cells CD4. AL357507.1 had a close relationship with dendritic cells resting, macrophages M1, neutrophils and plasma cells. CEBPA-DT was associated with monocyte, T cells CD4 naïve, T cells CD8 and T cells regulatory (Tregs). C5orf66-AS1 and AL512422.1 were only related to dendritic cells activated and dendritic cells resting, respectively. On the contrary, for AL360182.2 and AL008718.3, there was no association with immune cells.

Figure 6 LncRNA and immune analysis. (A) Relationship between high/low group and immune cells; (B) Relationship between seven prognosis-related lncRNAs and immune cells. *<0.05, **<0.01.

mRNA Expression of Seven lncRNAs in OS Cell Lines

The mRNA level of seven lncRNAs (AL512422.1, AL008718.3, C5orf66-AS1, AL360182.2, CEBPA-DT, AC006033.2, and AL357507.1) in five OS cell lines were examined. MNNG cell line was found to have a high OS metastasis rate, while the other four had a low OS metastasis rate (Figure 7). AL512422.1 was found to be highly expressed in the U2OS and SJSA-1 cell lines than MNNG (Figure 7A). AL360182.2 and C5orf66-AS1 mRNA levels were found to be higher in four weakly metastatic cell lines (Figure 7E and F). However, AL008718.3 mRNA level was significantly down-regulated in USO2, SAOS-2, SJSA-1, and HOS cell lines (Figure 7B). In addition, the mRNA level of AL357507.1 was only up-regulated in the SAOS-2 cell line (Figure 7D). Interestingly, even though they were all weak metastatic cell lines, the expression pattern differed. The mRNA level of AC006033.2 was higher in U2OS, SAOS-2, and HOS cell lines, and lower in SJSA-1 cell lines (Figure 7C). In the case of CEBPA-DT, the mRNA levels were higher in the HOS cell line and lower in the SJSA-1 cell line (Figure 7G).

Figure 7 The qPCR result of seven lncRNAs in five OS highly metastasis and weakly metastasis cell lines. (A) AL512422.1. (B) AL008718.3. (C) AC006033.2. (D) AL357507.1. (E) AL360182.2. (F) C5orf66-AS1. (G) CEBPA-DT. *<0.05, **<0.01, ***<0.001, ****<0.0001.

Discussion

OS is an aggressive malignancy with poor outcomes, particularly in patients with metastases.²⁹ The development of advanced treatments has resulted in significant improvement in the overall survival of patients with OS. Surgery is still the cornerstone of primary tumor treatment, whereas adjuvant chemotherapy is the most significant predictor of prognosis in patients with metastases.^30,31 Chemotherapies, such as high-dose methotrexate, doxorubicin, cisplatin, and ifosfamide/etoposide are the most effective agents. However, the incidence of metastasis remains high in a significant proportion of OS patients.³¹ Therefore, to predict or inhibit the metastatic process, it is essential to identify targets to provide potential diagnostic and therapeutic strategies.³²

Accumulating evidence suggests that lncRNAs are closely associated with tumorigenesis and progression. lncRNAs play an essential role in a variety of biological processes implicated in carcinogenesis, including cell proliferation, migration, and tumor metastasis.^33–36 Identifying an effective lncRNA biomarker may aid in the prevention of metastasis. For example, mtlncRNA LIPCAR is a novel biomarker for cardiac remodeling that can predict the survival rate of heart failure patients.³⁷ Immunoreactive mice showed that targeted therapy of ASncmtRNAs could prevent lung metastases following melanoma resection.³⁸ The lentiviral construct targeting ASncmtRNAs could inhibit B16F10 primary tumor proliferation.³⁹ In vivo studies in isogenic mouse renal cancer models revealed that ASncmtRNAs-targeting ASOs could completely reverse tumor growth.⁴⁰ In summary, lncRNA targeted therapy may be a potential therapeutic approach for improving prognosis.

In this study, we established a new risk-score model based on prognosis-related lncRNA generated from a series of bioinformatics analyses to predict the metastasis of OS patients. Based on their median risk score, all patients were classified as high or low-risk groups. In comparison to the low-risk group, patients in the high group had a poor survival rate. By calculating the AUC value, this model exhibited a strong prediction capacity. The co-expression network revealed an association between TF and 24 prognosis-related mRNAs, with the seven lncRNAs having the closest relationship. Zhong et al reported that HSF2 was involved in the growth and metastasis of lung cancer.⁴¹ It is important to note that the most common site of OS metastasis in the lungs and lung metastases are associated with poor prognosis.³⁰ NRF2 is a transcription factor that modulates the level of reactive oxygen species (ROS), detoxifying agents, and antioxidants in the body. The increase in NRF2 in the nuclei is closely related to bone metastases in OS patients.⁴² Furthermore, Cindy et al revealed that catecholamines and their receptors could be potential molecular markers for OS progression.⁴³ The catecholamine metabolic process and catechol-containing compound metabolic process were reported as a consequence of 24 prognosis-related mRNA enrichment.

Several overlapping pathways were discovered among the seven lncRNAs enrichment findings, confirming the critical biological function of the seven lncRNAs. The terms “NOTCH_SIGNALING”, “HYPOXIA”, “WNT_BETA_CATENIN_SIGNALING”, and “HEDGEHOG_SIGNALING” appeared more often. Several studies have shown that the Notch pathway regulates OS metastasis and may be useful for therapeutic targeting in the near future.^44,45 Fang et al reported that the Notch pathway could be pharmacologically inhibited, which was a huge step in targeted therapeutic strategies.⁴⁶ CNOT1, in collaboration with LMNA, may exacerbate OS tumorigenesis through the Hedgehog signaling pathway.⁴⁷ In addition, the hypoxia-HIF-1α-CXCR4 pathway plays a critical role during the migration of human OS cells.⁴⁸ Furthermore, NRP2 was shown to be overexpressed in OS cell lines; depleting NRP2 through downregulation of the active Wnt-signaling pathway would significantly reduce the tumor burden and metastasis rate in OS cell lines.⁴⁹ Additionally, we found that three of the seven lncRNAs were associated with resting dendritic cells. According to Jones et al, IL23 expression and monocyte-derived dendritic cells have both activating and suppressive effects on cytotoxic T cells, as well as direct effects on osteoblasts or associated osteoclasts.⁵⁰ It is worth mentioning that under hypoxic conditions, dendritic cells secrete significant amounts of osteopontin, which contributes to tumor metastasis.⁵¹

There were some limitations in the present study. One limitation is the relatively small sample size, even though recruitment is still ongoing. Another limitation is that no additional analysis stratified by age was conducted to account for metastasis. In conclusion, we established a novel model with excellent reliability and accuracy for predicting the metastasis of OS patients.