Analysis and identification of m6A RNA methylation regulators in metastatic osteosarcoma.

Osteosarcoma (OS) is characterized by rapid growth and early metastasis. However, its mechanism remains unclear. N6-methyladenosine (m6A) modification and its regulatory factors play essential roles in most cancers, including OS. In this study, we screened out 21 m6A modifiers using the Therapeutically Applicable Research to Generate Effective Treatments (TARGET) database, followed by the identification of the critical m6A methylation modifiers. The results revealed that the expression levels of three m6A methylation regulators, namely RBM15, METTL3, and LRPPRC, were associated with the low survival rate of patients with OS. We further studied the independent prognostic factors by performing univariate and multivariate Cox analyses and found that metastasis was an independent prognostic factor for patients with OS. Furthermore, we found for the first time that RBM15 was specific for metastatic OS rather than non-metastatic OS. Moreover, the significant overexpression of RBM15 was validated in metastatic OS cell lines and in actual human clinical specimens. We also revealed that RBM15 promoted the invasion, migration, and metastasis of OS cells through loss-functional and gain-functional experiments and an animal metastatic model. In conclusion, RBM15 has a high correlation with OS metastasis formation and the decreased survival rate of patients with OS, and this may serve as a useful biomarker for predicting metastasis and prognosis of patients with OS.

Introduction

Osteosarcoma (OS), characterized by features of rapid growth and early metastasis, is highly prevalent among children and teenagers. The annual incidence rate of OS has reached three cases per million individuals. Despite extensive research on OS conducted in the past 30 years, no targeted treatment is available for OS. Presently, the primary treatment for OS is chemotherapy, followed by the surgical removal of bone tumors. Although the treatment regimen for OS has improved markedly, and the 5-year survival rate for patients with OS has been maintained at 60%–75%, the prognosis remains poor and may cause the occurrence of other side effects that considerably influence the quality of life of patients. This influence is attributable to the fact that OS easily metastasizes, and its pathogenesis remains unclear. Therefore, a thorough understanding of the mechanism underlying the development of metastatic OS is essential for improving new OS therapeutic regimens.

In various physiological and pathological processes, RNA modification has been found to be a key regulator of transcript expression. As one of the most common types of RNA modification, N6-methyladenosine (m6A) has attracted considerable attention in recent studies. According to recent evidence, m6A modification is an epigenetic change with dynamic and reversible characteristics, which can regulate the splicing, nuclear export, positioning, translation, degradation, and stability of RNA through multiple regulatory proteins, including methyltransferases (METTL3, WTAP, METTL14, RBM15, CBLL1, ZC3H13, RBM15B, and KIAA1429, referred to as “writers”), demethylases (ALKBH5 and FTO, referred to as “erasers”), and critical m6A-binding proteins (HNRNPA2B1, YTHDF1-3, IGF2BP1-3, FMR1, LRPPRC, HNRNPC, ELAVL1, and YTHDC1-2, referred to as “readers”)., Several researchers have suggested that m6A regulators may serve as novel prognostic biomarkers for a wide variety of tumors. A high expression of METTL3 and YTHDF1 is correlated with the prognosis of patients with hepatic cell carcinoma (HCC). The expression of the m6A regulator HNRNPA2B1 reportedly increases markedly in esophageal cancer (ESCA); moreover, it plays an oncogenic role in the development of ESCA cells, thereby indicating that it may be a therapeutic target and a promising prognostic biomarker for human ESCA. HNRNPA2B1 expression can be significantly misregulated, thus serving as a reliable biomarker for endometriosis. These findings suggest that m6A regulators are strongly implicated in the development of cancer and may be promising prognostic biomarkers and therapeutic targets for cancer.

The m6A regulators are also linked with the progression and prognosis of OS. Overexpression of the m6A demethylase ALKBH5 notably restrains the growth, migration, and invasion of OS cells. Overexpression of METTL14, a type of m6A methyltransferase, significantly lowered the proliferation, migration, invasion, and viability of OS cells and promoted apoptosis. Moreover, the m6A methyltransferase WTAP downregulates homeobox-containing 1 (HMBOX1) expression in an m6A-dependent manner and regulates the PI3K/AKT pathway, which promotes the development of OS. However, current research mainly focuses on the m6A regulators at the early and middle stages of OS; these are unsuitable as prognostic biomarkers because of the principal factor of high mortality in patients with OS with their distant metastasis and frequent recurrence. Thus, to improve target therapy and to reduce the number of patients with OS mortality, it is essential to investigate the underlying mechanisms of m6A regulatory factors in metastatic OS.

Our study aimed to identify key m6A methylation modifiers in metastatic OS using a comprehensive bioinformatics analysis of data of OS samples from the Therapeutically Applicable Research to Generate Effective Treatments (TARGET) dataset, followed by the subsequent validation of the 22 metastatic OS samples. Furthermore, to predict the possible regulatory potential of the key m6A regulators on OS metastasis, cell colony formation, migration, Cell Counting Kit-8, wound-healing, and invasion assays were performed using OS cell lines in vitro and a metastatic experiment in vivo. Our study provides novel insights into metastatic OS therapy.

Results

Interaction and correlation of m6A RNA methylation regulation factors

Considering the important role played by m6A RNA methylation regulation factors in the initiation and progression of tumors, we used genome-wide data of 88 OS cases from the TARGET database to comprehensively study the effects and functions of 21 modulators of m6A RNA methylation in OS. We first analyzed the interaction between these 21 m6A RNA methylation regulation factors using a protein-protein interaction (PPI) network. The results of the analysis indicated that among the 21 m6A RNA methylation regulation factors, the 5 “writers” (METTL3, METTL14, WTAP, KIAA1429, and RBM15) exhibited the highest number of interactions with other factors (Figure 1A), with the highest number of related nodes (Figure 1B). Additionally, we examined the relationship between 21 regulators through the Pearson correlation test and found that the “eraser” FTO showed a negative correlation with ALKBH5 (correlation coefficient: −0.15) and RBM15 (correlation coefficient: −0.34), while it showed a positive correlation with METTL3, METTL14, and WTAP. The “writer” WTAP showed a positive correlation with METTL3 (correlation coefficient: 0.48) and METTL14 (correlation coefficient: 0.25). Interestingly, we also observed that the “reader” HNRNPC showed the least correlation with other factors (Figure 1C). Thus, most of the m6A RNA methylation regulation factors showed strong associations with each other. To investigate the underlying mechanisms of these regulation factors, we conducted a Gene Ontology (GO) analysis. As shown in Figure 1D, most regulation factors were enriched in the regulation of mRNA metabolic processes and the regulation of mRNA stability (Figure 1D and S1), which was in keeping with the functions of m6A RNA methylation regulators.

Figure 1

Correlation and interaction networks of 21 N6-methyladenosine modifiers

(A) The creation of the protein-protein interaction (PPI) network for the 21 modulators of N6-methyladenosine RNA methylation. (B) Interaction node counts of the 21 regulators in PPI networks. (C) The correlation among 21 N6-methyladenosine (m6A) modification regulators. (D) Gene Ontology analysis of 21 regulators.

Consensus clustering based on the expression levels of regulation factors

To explore the relationship between the expression level of key m6A RNA methylation regulation factors and the prognosis of patients with OS, all 88 OS cases were subjected to consensus clustering based on the expression level of 21 m6A RNA methylation regulation factors. By comparing the different number of clusters, a K of 3 was selected as an acceptable criterion; this was because the best tumor typing was observed when K = 3 (Figures 2A and 2B). Additionally, consensus clustering showed that the samples could be distinguished within three categories according to the degree of correlation existing between the groups (Figures 2C and S2A–2I). More importantly, we analyzed the overall survival of patients with these 3 clusters and observed a significant difference in overall survival probabilities between patients with sets 1, 2, and 3 (p < 0.05), which indicated a potential association between the m6A RNA methylation regulation factors and the prognosis of OS cases (Figure 2D). Furthermore, we compared the clinical and pathological characteristics (including metastasis, tumor site, age, and sex) of these three subgroups and found no remarkable differences between them. However, they exhibited distinct differences in the expression levels of a few m6A modifiers in different clusters, such as HNRNPA2B1 and ALKBH5. HNRNPA2B1 was expressed at a low level in subgroup 1 but was highly expressed in subgroups 2 and 3, whereas ALKBH5 exhibited the opposite trend. YTHDC2 was highly expressed in subgroup 2 but not in subgroups 1 and 3, while IGF2BP1 exhibited the opposite trend (Figure 2E).

Figure 2

Consensus clustering for m6A RNA methylation regulators and predictive analysis

(A) The consensus clustering of cumulative distribution function (CDF) when the index k ranged from 2–10. (B) The relative variation in the area under the CDF curve for the index k ranging from 2–9. (C) The consensus clustering matrix showing that the samples were divided into three clusters when k = 3. (D) The comparison of survival curves of three clusters (clusters 1, 2, and 3). (E) The heatmap and clinicopathologic features of the three clusters (clusters 1, 2, and 3) defined by the m6A RNA methylation regulators’ consensus expression in osteosarcoma.

Screening of key m6A modulators with poor prognosis in OS

The above-mentioned analysis implies the potential application of these 21 m6A modulator combinations in the prediction of OS overall survival probability. Subsequently, a survival analysis was conducted for each modulator of m6A RNA methylation to perform a screening of the key regulators related to the probability of overall survival of OS. As shown in Figure 3A, 3 modulators (RBM15, METTL3, and LRPPRC) simultaneously exhibited the smallest p value (p < 0.1) and the highest hazard ratio (HR > 1.99) in the survival analysis results. Herein, the increased expression of RBM15, METTL3, and LRPPRC was correlated with poor survival, indicating that they were high-risk markers for poor prognosis of OS. Thus, three candidate regulators, RBM15, METTL3, and LRPPRC, were selected for further analysis. To better forecast the clinical outcome of OS with dysregulated m6A RNA methylation regulation factors, we performed a least absolute shrinkage and selection operator (LASSO)-Cox regression analysis in R that was based on the expression levels of these three candidates in OS cases, according to the minimum criterion (Figures 3B and 3C).

Figure 3

Selection of three key modulators of m6A RNA methylation

(A) Overall-survival-related regulators in patients with osteosarcoma analyzed via univariate Cox regression hazard analysis. (B) The least absolute shrinkage and selection operator (LASSO) regression analysis of the selected three factors: RBM15, METTL3, and LRPPRC. (C) 10-fold cross-validation used for adjusting parameter selection in the LASSO regression model. Vertical solid lines indicate partial likelihood deviance with standard errors. The vertical dashed line represents the best value of the adjustment parameter (λ) in terms of the minimum standard.

Construction of the survival model consisting of three modulators of m6A RNA methylation

Based on the risk score results obtained herein, patients with OS were separated into 2 subgroups, namely high-risk and low-risk groups, with a cut-off point of 5 (risk score = 5). The rank of the risk scores, the status of survival time, and the RBM15, METTL3, and LRPPRC expression profile-based signatures in OS cases are shown in Figure 4A. The heatmap results showed that three m6A RNA methylation regulation factors, RBM15, METTL3, and LRPPRC, were highly expressed in the high-risk group. Moreover, we observed a remarkable difference in the overall survival rate between the two subgroups, in which compared with the patients with OS from the low-risk group, those in the high-risk group presented with a worse prognosis (p < 0.05) (Figure 4B). To test the sensitivity and specificity of the constructed prognostic model, a time-dependent receiver operating characteristic (ROC) curve analysis was conducted for the three factors; then, the area under the curve (AUC) value was determined. As shown in Figure 4C, the AUC value demonstrated that the signatures of the three modulators yielded an ideal model with high accuracy, with an AUC >0.7 (AUC = 0.83 in 1 year, AUC = 0.75 in 3 years, and AUC = 0.77 in 5 years) for predicting the survival rate of patients with OS. This indicates that they may exhibit crucial functions in OS. We then mapped the relationship between risk scores (high- and low-risk groups) and clinicopathological characteristics (including tumor site, metastasis, age, and sex). Understandably, most patients with OS from the high-risk group had an increased expression of RBM15, METTL3, and LRPPRC, but no significant association was established with the tumor site, metastasis, age, and sex (Figure 4D).

Figure 4

Generation of prognostic signature consisting of RBM15, METTL3, and LRPPRC

(A) The top panel shows the risk score distribution; the middle panel shows the survival status of each patient; and the bottom panel shows the heatmap illustrated for the analysis of expression of RBM15, METTL3, and LRPPRC. (B) The survival curves for low- and high-risk groups. (C) The accuracy of prognostic models was analyzed by performing receiver operator characteristic (ROC) curve analysis. The ROC curve and the area under the curve statistics show that the prognosis impact in 1, 3, and 5 years exhibits excellent specificity and sensitivity. (D) Heatmap illustrated for the analysis of expression of RBM15, METTL3, and LRPPRC and the exhibition of their clinicopathological characteristics. (E and F) The univariate (E) and multivariate (F) Cox regression analysis of the risk score and variables of clinical characteristics in patients with OS.

Univariate and multivariate Cox regression analyses were performed to assess whether the risk signature was an independent prognostic indicator for patients with OS. The univariate analysis for all variables exhibited that the poor prognosis of patients with OS was significantly correlated with metastasis (p < 0.001, HR = 4.770), risk score (p < 0.05, HR = 0.401), and cluster (p < 0.05, HR = 1.802), but they were not associated with sex, age, or the tumor site (Figure 4E). Cox multivariate regression analysis for all variables also indicated that metastasis (p < 0.001, HR = 4.038) significantly influenced the overall survival of patients with OS (Figure 4F).

Expression of m6A RNA methylation modulators in metastatic and non-metastatic OS samples

Based on the results of the univariate and multivariate analyses that indicated that metastasis was an independent risk factor for patient prognosis, we then investigated the expression profile of 21 m6A RNA methylation regulation factors. We used the RNA sequencing data of metastatic (n = 66) and non-metastatic (n = 22) OS samples derived from the TARGET database to perform a screening of the most critical modulators of m6A RNA methylation. The differential expression analysis results revealed that among the 21 m6A regulatory factors, the expression of RBM15 and WTAP was significantly upregulated in the metastatic OS samples (p < 0.05) when compared with the non-metastatic OS samples. However, combining these results with those described above, we found that among these three candidate regulators (METTL3, LRPPRC, and RBM15) associated with poor prognosis, only RBM15 was significantly overexpressed in metastatic OS samples (Figure 5). Thus, RBM15 may be one of the key m6A RNA methylation regulators that may have a vital function in the metastatic mechanisms of OS. Therefore, we focused on the functions of RBM15 in further analyses.

Figure 5

Dot plot generated for the analysis of expression of 21 modulators of m6A RNA methylation in the non-metastatic and metastatic osteosarcoma samples

Expression of RBM15 in other types of cancer

Considering the limited sample size of OS from the TARGET database and the different subtypes of OS with additional features, we analyzed the copy number of RBM15 in different types of OS using the Oncomine database. RBM15 exhibited exceptionally high frequencies of copy number gain events in four subtypes of OS (chondroblastic OS, telangiectatic OS, fibroblastic OS, and osteoblastic OS) (Figure S3A). Interestingly, in recurrent OS samples, the copy number of RBM15 was high in five out of seven samples (log2 copy number units > 0) (Figure S3B), further indicating the crucial role of RBM15 in OS. Furthermore, we also investigated the expression of RBM15 in other types of tumors and observed that RBM15 was highly expressed in several tumors, such as cervical cancer (Figure S4). We also observed that the expression level of RBM15 was related to disease-specific survival in a variety of tumors by evaluating the relationship between them in 33 tumors using data obtained from The Cancer Genome Atlas (TCGA) (Figure S5). These data suggest that RBM15 may also play a vital role in the development of other types of cancer and that it exhibits a specific association with the survival time of patients; however, the underlying mechanism warrants further analysis.

Mutation analysis of RBM15 and its relationship with immune-infiltrating cells

To further explore the functions of RBM15 in OS, we analyzed the mutation status of RBM15 using sarcoma (SARC) mutation data. As shown in Figure 6A, we observed that RBM15 showed a somatic mutation rate of 0.84%, in which two missense mutations might occur in the RNA recognition motif (Figure 6A). Furthermore, possible methylation sites of RBM15 in different OS cell lines were predicted using the Cancer Cell Line Encyclopedia (CCLE) database (Figure 6B). Moreover, we further identified proteins that established interactions with RBM15 by conducting the PPI network analysis; we found that RBM15 mostly established interactions with two regulators, namely RRN3 and DIDO1 (Figure S6A). We observed that RBM15 showed a positive correlation with both regulators (Figures S6B and S6C). Independent tumor-infiltrating immune cells may be used to predict overall survival efficiently. Therefore, we examined the possible association between the expression of RBM15 and the level of immune infiltration in OS using the Tumor Immune Estimation Resource (TIMER) database. As shown in Figure 6C, RBM15 expression showed a negative correlation with CD4 T cells, Timer macrophages, and dendritic cells (p < 0.05), while there was no significant association with the other three immune infiltrating cells (B cells, CD8 T cells, and Timer neutrophils). These observations demonstrate that RBM15 may play a key role in the immune infiltration of OS cells.

Figure 6

Analysis of the mutation pattern and the methylation sites of RBM15 and immune-infiltration cells

(A) The analysis of mutation pattern of RBM15 in the sarcoma samples. (B) The distribution of methylation sites of RBM15 in different osteosarcoma cell lines. (C) Correlations between RBM15 expression and immune-infiltration levels of six immune cells in OS.

Verification of RBM15 expression in the metastatic OS cell line and in actual human clinical specimens

Based on the bioinformatics analysis results, which showed that RBM15 was highly expressed in metastatic OS samples, and the findings that its relationship was associated with a poor prognosis in patients with OS, we investigated the expression level of RBM15 in metastatic and non-metastatic OS cell lines by performing quantitative real-time PCR and western blot assays. As shown in Figure 7A, our results showed that the mRNA expression of RBM15 was distinctly increased in metastatic OS cell lines (MNNG and 143B) compared with non-metastatic OS cell lines (HOS and MG63). The western blot analysis also showed the same results (Figure 7B). In addition, we also detected RBM15 expression in actual human clinical specimens using immunohistochemistry (IHC). As shown in Figure 7C, RBM15 protein expression was markedly increased in metastatic OS tissue compared with non-metastatic OS tissue. To the best of our knowledge, there has been no unbiased study conducted on the function of RBM15 in metastatic OS. To further determine the role of RBM15 in metastatic OS, we conducted the loss-of-functional and gain-of-functional analyses. We firstly transfected three siRNAs against RBM15 into MNNG and 143B cells, and then we assessed the inhibitory effect of the RBM15 siRNA by quantitative real-time PCR analysis. RBM15 expression was significantly reduced after RBM15-specific siRNA transfection in 143B cells and MNNG cells (Figure S7A and S7B). Based on the results of siRNA transfection, the most effective target for RBM15 interference was chosen as the lentivirus-mediated sh-RBM15 target to generate stable knockdown of RBM15 in metastatic OS cells (MNNG cell line), while the expression vector of RMB15 was transfected into non-metastatic OS cells (HOS cell line). As shown in Figures 7D and 7E, RBM15 expression was significantly increased after transfection of the expression vector of RMB15 in the HOS cell line. Meanwhile, in the MNNG cell line, RBM15 expression was significantly decreased after RBM15 knockdown (Figures 7F and 7G). These results indicate that we successfully generated OS cell lines with RBM15 overexpression or knockdown.

Figure 7

Validation of RBM15 expression in OS cell lines and in human clinical specimens

(A) Quantitative real-time PCR analysis of RBM15 expression in the osteosarcoma cell lines. (B) Western blot analysis of RBM15 expression in osteosarcoma cell lines. GAPDH was used as a loading control. (C) The representative image of RBM15 protein expression in human primary and metastatic osteosarcoma tissue by immunohistochemistry. (D) Quantitative real-time PCR validation of RBM15 overexpression in the HOS cell line. (E) Western blot validation of RBM15 overexpression in the HOS cell line. GAPDH was used as a loading control. (F) Quantitative real-time PCR validation of RBM15 knockdown in the MNNG cell line. (G) Western blot validation of RBM15 knockdown in the MNNG cell line. GAPDH was used as a loading control. All bar plot data have been presented as mean ± standard deviation of three independent experiments. ∗∗p < 0.01, ∗∗∗p < 0.001.

RBM15 promotes the migration, proliferation, and metastasis of OS cells

To further determine the regulatory potential of RBM15 on OS metastasis, we performed a Transwell migration experiment. As shown in Figure 8A, the migration ability of the HOS cells was significantly enhanced after RBM15 overexpression compared with that of the negative control. In contrast, the migratory ability was prominently obstructed after RBM15 knockdown in MNNG cells. Moreover, the invasion ability in HOS cells with RBM15 overexpression was also augmented, but inhibited, in MNNG cells with RBM15 defects, as evidenced by a Matrigel invasion assay (Figures 8A and S7C). Additionally, the wound-healing assay results further confirmed the effect of RBM15 on the migration of OS cells, in which RBM15-overexpressing HOS cells healed better than normal untreated cells when subjected to scratches, whereas the MNNG cells with RBM15 deletion exhibited worse healing effects on subjection to scratches than did normal non-treated cells (Figures 8B and S8A). Alterations in the proliferation capability were further assessed by a colony-forming assay. As shown in Figure 8C, overexpression of RBM15 significantly enhanced the colony-forming ability of the non-metastatic OS cell line, while the silencing of RBM15 remarkably disrupted the colony-forming ability of the metastatic OS cell line (Figures 8C and S8B). In addition, to further confirm the influence of RBM15 on the proliferative capacity of OS cells, we performed a Cell Counting Kit-8 (CCK-8) assay. As shown in Figure 8D, RBM15 exhibited a critical impact on maintaining a high proliferation rate in OS cells, in which the overexpression and knockdown of RBM15 significantly enhanced and suppressed the cell viability in OS cells, respectively. To further investigate the metastatic potential of RBM15 in vivo, MNNG cells with RBM15 knockdown were injected into the tail vein of nude mice. After 4 weeks, the lungs of nude mice were collected and analyzed. As shown in Figure 8E, the loss of RBM15 remarkably reduced the pulmonary metastasis of OS cells compared with negative control. Altogether, these findings reveal that RBM15 played a role in the motility and metastasis of OS cells.

Figure 8

Validation of the role of RBM15 on the invasion, migration, proliferation, and metastasis of osteosarcoma

(A) The effect of RBM15 on cell migration and invasion ability was investigated in HOS and MNNG cells by Transwell assay. The left panel shows the representative images; scale bars: 100 μm. The quantitative analyses of the cell migration and invasion level have been shown in the right panel. (B) The effect of RBM15 on the cell migration capability was evaluated by a wound-healing assay in HOS and MNNG cells. The left panel shows the representative images acquired under an inverted microscope; scale bars: 100 μm. The right panel shows the percentage of areas exhibiting relevant healing. (C) The proliferation capacity of HOS cells with RBM15 overexpression and MNNG cells with RBM15 knockdown was detected by a colony formation assay. The left panel shows the overall view of the colony formation in the entire dish. The right panel shows the number of cell colonies in each dish. (D) The effects of RBM15 knockdown and overexpression on cell vitality were determined by a CCK-8 assay. All bar plot data have been presented as mean ± standard deviation of three independent experiments. ∗∗p < 0.01, ∗∗∗p < 0.001. (E) Representative images of gross view of lungs and H&E staining of lung sections of nude mice in the model injection with indicated OS cells through the tail vein.

Discussion

The pathogenesis of OS is a complicated multistep and multifactorial process implicated in a broad spectrum of molecular abnormalities and tumor heterogeneity. It has been reported that m6A is the most plentiful RNA inner modification in eukaryotic cells, and abnormal m6A RNA methylation plays an essential role in tumorigenesis through multiple molecular mechanisms. However, the precise functions of these regulators in the modification of m6A RNA methylation in metastatic OS remain unknown.

In our study, a multivariate Cox proportional risk model was employed to identify 21 prognostic regulators of m6A RNA methylation in OS. The results of the analysis demonstrated that the high expression of METTL3, LRPPRC, or RBM15 was validated as a valuable prognostic risk factor for the low overall survival of patients with OS. A LASSO-Cox regression analysis and the construction of survival models to perform screenings of the tumor markers have been widely conducted in studies of several cancers, such as in liver cancer, pancreatic cancer, and bladder cancer. Three m6A RNA methylation regulators (RBM15, METTL3, and LRPPRC) were selected to construct a prognostic risk model. According to the risk characteristics, patients with OS were separated into two subgroups, and a lower survival rate was observed in the high-risk group. The ROC analysis revealed an AUC >0.7, which revealed that the risk profile exhibited a good performance. In conclusion, we constructed a reliable prognostic model and successfully predicted that the highly expressed METTL3, LRPPRC, and RBM15 regulators might be potential biomarkers of OS.

For the three candidate regulators, only RBM15 showed high expression levels in the 22 metastatic OS samples and in metastatic OS cell lines (MNNG and 143B cell lines). Furthermore, previously it has been reported that the increased expression of RBM15 was observed in hepatocellular carcinoma, testicular germ cell tumors, and laryngeal squamous cell cancer. RBM15 may also promote the progression of acute/chronic myeloid leukemia and Kaposi's SARC. By comparing the expression levels before and after RBM15 silencing in the MNNG cell line and RBM15 overexpression in the HOS cell line, our data uncovered that the suppression of RBM15 expression in OS inhibited cell proliferation, migration, and invasion, consistent with the results reported in a previous study. Importantly, we further confirmed the metastatic potential of RBM15-knockdown cell lines in vivo, indicating that RBM15 could promote the metastasis of OS cells. Thus, our findings suggest that RBM15 plays a non-negligible role in the metastasis of OS, which may be a potential predictor for OS survival outcomes. We also demonstrated that this prognostic signal could be used as a useful tool to predict OS survival outcomes.

Furthermore, we demonstrated that RBM15 showed a negative correlation with CD4 T cells, Timer macrophages, and dendritic cells in our bioinformatics analysis. This indicates that RBM15-mediated m6A RNA methylation may be implicated in the immune infiltration of OS. RBM15 is a critical methyltransferase that regulates the fate of cellular mRNA molecules by catalyzing the formation of m6A in mRNA; it has been previously uncovered in cortical development, liver maturation, and Xist-mediated X chromosome inactivation., It is worth noting that we downloaded mutation data from TCGA for mutant processing of SARC and found that RBM15 is mainly mutated in a RNA recognition motif (RRM) with a somatic mutation rate of only 0.84%. It is possible that the mutated forms of RBM15 may not be for OS, especially metastatic OS, due to the mutation data that is derived from SARC, which includes other SARC types. And, to clarify this point, further genome sequencing of the gene is required in future studies. We also found that RBM15 established interactions with two proteins, RRN3 and DIDO1, using a PPI analysis. RRN3 is a transcription initiation factor that acts as a bridge between RNA polymerase I and the pre-initiation complex at the ribosomal DNA promoter to initiate pre-ribosomal RNA transcription. It has been demonstrated to participate in the differentiation and proliferation of different mammalian cells, thereby triggering the initiation of breast cancer by damaging mammary epithelial morphogenetic processes., As another protein that establishes interactions with RBM15, DIDO1 has been revealed to be highly expressed in several tumor types, such as in bladder cancer and melanoma. Furthermore, DIDO1 affects the invasion, proliferation, and migration of tumor cells.29, 30, 31 Thus, it is reasonable to speculate that RBM15 promotes the development and metastasis of OS, possibly by regulating the expression and function of RRN3 and DIDO1 through m6A RNA methylation. These hypotheses should be further verified with future experimentation.

In conclusion, we analyzed the expression and prognostic value of 21 modulators of m6A RNA methylation in metastatic OS and further validated these using cytological experiments in vitro, both before and after RBM15 silencing and overexpression, and a metastatic model in vivo. We identified that the m6A RNA methyltransferase RBM15 played a critical role in the development and metastasis of OS; thus, it might be a useful new candidate biomarker for the prognosis of patients with metastatic OS.

Materials and methods

Data acquisition

The RNA sequencing transcriptome data of 88 patients with OS along with their corresponding clinicopathological information were first downloaded from the TARGET database (https://ocg.cancer.gov/programs/target). For the RNA sequencing data, TARGET samples (n = 88) were normalized based on the expression levels of genes determined by fragments per kilobase of exon model per million mapped reads. Detailed, clinically relevant information on the 88 patients is listed in Table 1.

Table 1

Clinicopathological characteristics of patients with OS from TARGET

Characteristics	OS patients (n = 88)
	Amount	Percentage (%)
Sex
Male	51	57.95
Female	37	42.05
Age
<18	68	77.27
>18	20	22.73
Metastatic status
Metastasis	22	25.00
Non-metastasis	66	75.00
Specific tumor site
Femur	40	45.98
Other	47	54.02
Living state
Alive	57	66.28
Dead	29	33.72

Extraction of data on 21 m6A RNA methylation regulation factors

According to the TARGET data on OS and information obtained from other published studies,32, 33, 34 21 modulators of m6A RNA methylation have been investigated in this study, including 8 methyltransferases, 2 demethylases, and 11 m6A-binding proteins.

Correlation analysis and creation of PPI network

The PPI network comprising 21 selected modulators of m6A RNA methylation was constructed to understand their molecular regulatory mechanisms through the use of the Retrieval of Interacting Genes protein database 11.0 (STRING, http://string-db.org/). Writers, erasers, and readers were organized within the PPI network. The correlation of 21 modulators of m6A RNA methylation was investigated by performing the Pearson correlation test, and the visualization of results was performed utilizing the R “corrplot” package.

GO analysis

To understand the molecular functions of the 21 regulation factors of m6A RNA methylation, we performed a GO analysis using the “clusterProfiler” package in R. Moreover, we used the R “GO-Plot” visualization package to better visualize the interaction of the 21 regulation factors, as well as their functional categories.

Screening of the key regulators

To perform a screening of the three key m6A RNA methylation regulation factors in connection with the probability of overall survival in OS samples, we carried out a univariate Cox regression analysis for the 21 regulation factors. The LASSO-Cox regression analysis was used to create a prognostic model for the moderators of m6A RNA methylation in OS.

Establishment of the prognostic model and illustration using a heatmap

To compare the difference in overall survival of patients with OS in the high-risk and low-risk groups, the Kaplan–Meier analysis method was used. Moreover, a ROC curve was produced, and the area under the ROC curve was determined to assess the sensitivity and specificity of the prognostic signature in OS samples. Clinical information (tumor site, metastasis, age, and sex) based on the risk groups related to m6A RNA methylation regulation factors was illustrated using a heatmap by applying the “pheatmap” package in R. To obtain the correlation between different clinical features (including sex, age, metastasis, site, cluster, and risk) and the overall survival rate for patients with OS, both univariate and multivariate Cox regression analyses were conducted.

Analysis of the mutation pattern and methylation sites of RBM15

The mutation pattern of RBM15 was analyzed using mutation data obtained from the SARC sample from TCGA database. The methylation sites of RBM15 were evaluated in different OS cell lines using the CCLE database.

Expression of m6A RNA methylation modulators in other OS samples

The Oncomine (https://www.oncomine.org/) database was used as a data platform to further examine the expression levels of selected m6A RNA methylation modulators in both the recurrent and different subtypes of OS samples.

Relationship between RBM15 expression and immune-infiltrating cells

The data on the scores of the six immune-infiltrating cells in SARC, acquired from the Tumor Immune Estimation Resource (TIMER) database (https://cistrome.shinyapps.io/timer), were subjected to further analysis. This was performed to assess whether RBM15 was associated with the six immune infiltrating cells.

Cell culture

Human OS cells (HOS, MG63, MNNG, and 143B cell lines) were obtained from the American Type Culture Collection (ATCC, Manassas, VA, USA). The MG63 and HOS were used as non-metastatic OS cells, while 143B and MNNG were used as metastatic OS cells. All cells were incubated in the Dulbecco's modified Eagle medium (Gibco, Shanghai, China) supplemented with 10% (vol/vol) fetal bovine serum (Tianhang, Zhejiang, China), 100 μg/mL streptomycin, and 100 U/mL penicillin (Solarbio, Beijing, China). The OS cells were grown in plates at 37°C in an incubator with a humid atmosphere and 5% CO2.

Transfection

Three small interfering RNAs (siRNAs) against RBM15 (si-RBM15#1 sequence: sense (5′-3′) CCGCAACAAUGAAGGGAAATT, antisense (5′-3′) UUUCCCUUCAUUGUUGCGGTT; si-RBM15#2 sequence: sense (5′-3′) GGGUUCAACUGGAGGCAAATT; antisense (5′-3′) UUUGCCUCCAGUUGAACCCTT; si-RBM15#3 sequence: AGGUGAUAGUUGGGCAUAUAUTT, antisense (5′-3′) AUAUAUGCCCAACUAUCACCUTT) and control non-targeting siRNA (si-Ctrl sequence: sense (5′-3′) UUCUCCGAACGUGUCACGUTT, antisense (5′-3′) ACGUGACACGUUCGGAGAATT) were provided by Sangon Biotech (Shanghai, China). Transient transfection was conducted with Lipofectamine 3000 (Invitrogen, Carlsbad, CA, USA) following the manufacturer’s protocols. Briefly, the 143B and MNNG cells were cultured in a 6-well plate for 24 h, and, subsequently, 30 pmol of siRNAs and 5 μL of the Lipofectamine 3000 transfection reagent were added after they were mixed for 20 min at room temperature (20°C). The transfection efficacy was evaluated by analyzing the expression level of RBM15 using a quantitative real-time PCR assay. A human lentivirus-mediated sh-RBM15 and an expression vector of RMB15 were generated by GeneChem (Shanghai, China). The lentiviruses were ultracentrifuged, concentrated, validated, and added to the cell culture medium. After infection, the cells were selected with puromycin (Sangon Biotech, Shanghai, China) for 10 days, and the surviving cells were continuously cultured as stable mass transfectants.

Quantitative real-time PCR assay

Total RNA was purified from OS cells by utilizing the Magen kit (Guangzhou, China). According to the instructions, 1 μg of RNA was reverse-transcribed into cDNA using the PrimeScriptTM RT reagent Kit (Takara, Beijing, China). For mRNA analyses, SYBR Premix Ex Taq II (Takara, Beijing, China) was used, and PCRs were conducted using a LightCycler@96 PCR instrument (Roche, Switzerland) following the manufacturer’s instructions. The thermocycling conditions were set at 95°C for 600 s, 45 cycles of 95°C for 10 s, and followed by 60°C for 1 min, as per methods described previously. The expression level of the key genes was normalized to that of glyceraldehyde phosphate dehydrogenase (GAPDH) and was analyzed utilizing the comparative method 2−ΔΔCt. The PCR primers used in this study are listed in Table 2.

Table 2

PCR primers used in this study

Gene	Sequence (5′ to 3′)
RBM15-Forward	GAGTTCTCCCAGCAGTTCCT
RBM15-Reverse	TATAACAGGGTCAGCGCCAA
GAPDH-Forward	CCACTCCTCCACCTTTGAC
GAPDH-Reverse	ACCCTGTTGCTGTAGCCA

Wound-healing assay

The metastatic OS cell lines were cultured in a 6-well plate; thereafter, a line was scratched onto the cell monolayer in each well after the cells reached 100% confluence using a sterile disposable p200 pipette tip. The cells were then subjected to washing steps using fresh non-serum medium and cultured for 24 h. An image of the scratched area was acquired at 0 h; then, the same region was photographed again 24 h using a microscope (Leica, Germany) at 10× magnification after scratching was performed. The relative migration rate was evaluated by measuring the healing distance across the wound region, which was normalized to the original scratched part at 0 h. These assays were performed three times, with three replicates each.

Migration assay and invasion analysis

The migration assay was carried out using a Transwell plate (24-well) based on the manufacturer’s protocols. Briefly, the OS cells (1 × 105) were grown in the upper chamber of the Transwell plate, while the lower chambers were immersed in 1 mL of a complete medium. After the completion of 24 h of incubation at 37°C, the OS cells that did not migrate to the upper layer were discarded, while the cells that migrated to the lower layer of the chamber were subjected to staining using a 0.5% crystal violet solution after fixation with 4% paraformaldehyde for 30 min. Representative images were obtained using an inverted fluorescence microscope (Leica, Germany). For the invasion assay, the Transwell chamber was coated with Matrigel (Solarbio, Beijing, China), and a procedure similar to that described above was conducted according to the manufacturer’s protocols. The number of migrated or invaded cells was determined to evaluate the migration and invasion rates based on the three randomly acquired field digital pictures. Based on these triplicate experiments, the average number of migrated or invaded cells and the standard deviation values were obtained.

Cell colony formation experiments

The cells (1,000 cells/well) were placed in 12-well plate and cultured for 5 days. Next, the cells were subjected to washing steps using PBS after the removal of the culture medium. Then, 0.5% crystal violet solution was used to stain the cells after fixation with 4% paraformaldehyde. Images were acquired under an camera (Cano, Japan), and the cell colonies were counted in each well. This experiment was repeated three times.

Western blot

Proteins were extracted from OS cell lines with a radio-immunoprecipitation assay (RIPA) lysis buffer. Lysates were sonicated on ice and centrifuged at 12,000 RPM for 20 min at 4°C. The protein concentration of the supernatant was determined using a BCA protein assay kit. Equal amounts of protein (60 μg) were separated by 10% SDS-PAGE and transferred to PVDF membranes. The membranes were incubated with blocking buffer 5% non-fat milk in TBS containing 0.1% Tween 20 (TBST) for 2 h at room temperature and then probed with the primary antibodies against RBM15 (dilution 1:2,000) and GAPDH (dilution 1:2,000) overnight at 4°C. After washing three times with TBS containing 0.1% Tween 20 for 5 min, the membranes were then incubated with the secondary antibody by Sangon Biotech (Shanghai, China). Finally, High Sensitivity Plus ECL luminescence reagent (Shanghai, China) was added to the PVDF membranes, and the color was developed in Amersham Imager 600 (General Electric, USA).

CCK-8 assay

The transfected cells were seeded in 96-well plates at a concentration of 2 × 103 per well at different time points (0, 24, 48, and 72 h), 10 mL of CCK-8 reagent (Biosharp, Beijing, China) was added to each well after cell attachment, and cells were incubated at 37°C for 2 h. The absorbance of the solution was measured at 450 nm using a Labsystems (Thermo, Waltham, MA, USA).

IHC

Primary OS and lung metastatic OS tissues from the same patient were collected from Guangxi Medical University Cancer Hospital. IHC staining was performed to identify the protein expression of RBM15 in OS tissues. Briefly, the OS tissue sections were incubated with 3% H2O2 for 15 min and blocked by goat serum for another 15 min at room temperature. Subsequently, knee tissue sections were incubated with the primary antibody of RBM15 (1:200 dilution; Proteintech, China) at 37°C for 90 min in the humidified chamber. After being rinsed three times with PBS, OS tissue sections were incubated with the biotin-labeled goat anti-mouse/rabbit immunoglobulin G (IgG) (ZSGb Bio, China) for 30 min. Then, samples were counterstained with hematoxylin (Solarbio, China) after being stained with 3, 3′-diaminobenzidine tetrahydrochloride (DAB; Boster, China). Finally, the tissue sections sealed with neutral resin were imaged using the upright microscope (Olympus BX53, Tokyo, Japan).

Human tumor samples and animal experiment

All procedures involving human tumor samples in this research were approved by the Medical Ethic Committee of Guangxi Medical University. All the animal operations involved in this study were carried out as prescribed by the Guide for the Care and Use of Laboratory Animals published by the National Institutes of Health (Eighth Edition) and were approved by the Animal Ethic Committee of the Guangxi Medical University. The MNNG cells with RBM15-knockdown (sh-RBM15) or control cells (sh-ctrl) (equivalent volume and density) were suspended in sterile PBS and injected into the tail vein of each nude mouse (male, 4 weeks old). After 4 weeks, the lungs of nude mice were taken for observation and analysis. Then, the samples were embedded in paraffin and sliced into 5 μm sagittal tissue sections. The slices were stained with hemotoxylin and eosin (H&E).

Statistical analysis

The normality test of the expression of RNA throughout the OS samples was performed using the Kolmogorov-Smirnov test. A two-sided t test was performed to estimate gene expression levels in metastatic and non-metastatic OS cells. Statistical significance was set at p < 0.05.