Zhongguang Wu
Xiaobo Zhang2†
Dongjie Chen
Xin Wu
Youwen Deng- 1Department of Laboratory Medicine, The Third Xiangya Hospital, Central South University, Changsha, China
- 2Department of Spine Surgery, Third Xiangya Hospital, Central South University, Changsha, China
- 3Department of Hepatopancreatobiliary Surgery, The Third Xiangya Hospital, Central South University, Changsha, China
- 4Department of Clinical Laboratory, Qinghai Provincial People’s Hospital, Xining, China
N6-Adenosine methylation, yielding N6-methyladenosine (m6A), is a reversible epigenetic modification found in messenger RNAs and long non-coding RNAs (lncRNAs), which affects the fate of modified RNA molecules and is essential for the development and differentiation of immune cells in the tumor microenvironment (TME). Osteosarcoma (OS) is the most common primary bone tumor in children and adolescents, and is characterized by high mortality. Currently, the possible role of m6A modifications in the prognosis of OS is unclear. In the present study, we investigated the correlation between m6A-related lncRNA expression and the clinical outcomes of OS patients via a comprehensive analysis. Clinical and workflow-type data were obtained from the Genotype-Tissue Expression Program and The Cancer Genome Atlas. We examined the relationship between m6A modifications and lncRNA expression, conducted Kyoto Encyclopedia of Genes analysis and also gene set enrichment analysis (GSEA), implemented survival analysis to investigate the association of clinical survival data with the expression of m6A-related lncRNAs, and utilized Lasso regression to model the prognosis of OS. Furthermore, we performed immune correlation analysis and TME differential analysis to investigate the infiltration levels of immune cells and their relationship with clinical prognosis. LncRNA expression and m6A levels were closely associated in co-expression analysis. The expression of m6A-related lncRNAs was quite low in tumor tissues; this appeared to be a predicting factor of OS in a prognostic model, independent of other clinical features. The NOD-like receptor signaling pathway was the most significantly enriched pathway in GSEA. In tumor tissues, SPAG4 was overexpressed while ZBTB32 and DEPTOR were downregulated. Tissues in cluster 2 were highly infiltrated by plasma cells. Cluster 2 presented higher ESTIMATE scores and stromal scores, showing a lower tumor cell purity in the TME. In conclusion, m6A-related lncRNA expression is strongly associated with the occurrence and development of OS, and can be used to as a prognostic factor of OS. Moreover, m6A-related lncRNAs and infiltrating immune cells in the TME could serve as new therapeutic targets and prognostic biomarkers for OS.
Introduction
Epigenetic silencing, for instance via histone modifications or RNA methylation, is a crucial tumorigenic mechanism. In fact, epigenomic alterations affect tumor immunogenicity, and immune cells are involved in antitumor responses (1). Furthermore, epigenetic reprogramming via inhibition of methylation enhances plasma cell polarization and decreases the response efficiency of antitumorigenic T cells (2). Cell differentiation, immunity, inflammation, and carcinogenesis are all regulated by non-coding RNAs (ncRNAs) such as microRNAs (miRNAs), long non-coding RNAs (lncRNAs), and circular RNAs (circRNAs) (3). Recent epigenetic studies have focused on N6-methyladenosine (m6A). To develop effective antitumor treatments and predict individual survival and recurrence risks, it is essential to unravel the role of m6A-related lncRNAs in the development and differentiation of immune cells.
Among the most common invasive bone tumors, human osteosarcoma (OS) is characterized by a high mortality rate in children and adolescents (4). Specifically, approximately 5% of all childhood malignancies and approximately 9% of tumor-related deaths depend on OS (5). This disease often affects the distal femur, proximal tibia, and humerus (6). The prognosis of OS patients is uncertain even when they are subjected to a combination of treatments, including chemotherapy and amputation (5). Despite significant advances in treatment strategies over the past few decades, the overall survival rates of OS patients remain low, especially for metastatic OS (7). Therefore, it is crucial to identify novel biomarkers and ensure effective treatment for patients with OS.
The most common internal chemical modification of mRNA is N6-adenosine methylation, a post-transcriptional modification of RNA yielding N6-methyladenosine (m6A) (8). m6A network components are distinguished into three subtypes: writers, readers, and erasers (9). Writers facilitated the N6-Adenosine methylation and erasers removed; moreover, m6A recognition involves specific reader proteins (10). These m6A-related regulatory factors play crucial roles in many physiological and pathological processes by modulating RNA stability, mRNA translation and splicing, and miRNA processing (11–13). In addition, N6-adenosine methylation contributes to the initiation and development of several cancers, namely, breast cancer, liver cancer, colorectal cancer, glioma, hepatoblastoma, and cervical cancer (14–19). Nevertheless, the role and clinical value of m6A-related regulatory factors remain unclear in OS.
m6A-Methylated lncRNAs regulate a number of pathological and biological processes. For example, by inducing mRNA splicing to mediate immunotherapy and regulate gene expression, m6A residues prolong patient survival and affect lactate levels in the TME (20). Additionally, overexpression of METTL14 significantly inhibited the proliferation, invasion, migration, and apoptosis of OS cells, whereas knockdown of METTL14 exerted the opposite effects (21). Nevertheless, the causes of abnormal expression and methylation of lncRNAs in OS are unknown. Thus, it is vital to construct transcriptome maps of lncRNA expression profiles and m6A modification profiles to understand how these processes affect the prognosis of OC.
Patients with OS can benefit from immune cell score analysis. Indeed, immune risk scores aid in assessing prognosis, facilitate adjuvant treatment of OS, and help explain the link between chemotherapy and immune response, thereby demonstrating the importance of immunotherapy (22). As previously reported in solid tumors, treatment outcomes can be influenced by the extent of immune cell infiltration. According to the type of immune cells that are present in the TME, cancers can be distinguished into tumors that are immunologically active or inflammatory (hot) and tumors that are inactive or non-inflammatory (cold). Cold tumors are associated with poor prognosis because they are insensitive to chemotherapy and immunotherapy. A characteristic feature of febrile neoplasms is the infiltration of leukocytes, including CD8+ T cells (23). The immune score is a cytotoxic immune response index that can improve the prognosis of collision tumors and guide therapeutic strategies (24). However, few studies have examined the relationship between m6A-associated lncRNAs and immune cell infiltration in OS. Therefore, it is necessary to study immune cell infiltration within the TME to clarify the relationship between the latter and OS clinicopathological parameters.
This study focused on the epigenetic mechanisms by which m6A-related lncRNAs mediate the development and differentiation of immune cells. We first evaluated the prognostic value and expression status of m6A-related lncRNAs by analyzing public genome databases of OS patients. Additionally, we developed clustering subtypes based on the expression levels of prognostic m6A-related lncRNAs to determine their relationships with N6-adenosine methylation, TME scores, target gene expression, prognosis, and infiltration of immune cells. Pathway and functional gene set enrichment analyses (GSEA) allowed to further pinpoint the pathways involved in the process. Finally, Lasso regression was used to test the constructed m6A-related lncRNA prognostic model.
Materials and Methods
Sample Data Acquisition and Collation
The TARGET-OS gene expression (RNA-Seq) dataset (https://ocg.cancer.gov/programs/target) provided clinical data of 88 patients with OS. The GTEx datasets were obtained from the Genotype-Tissue Expression database (http://commonfund.nih.gov/GTEx). Using PERL software (https://www.perl.org/), we created an mRNA matrix and processed transcriptomic data with the corresponding script for gene ID transformation, which enabled to process the clinical data. The software was installed in a path named in English with no spaces.
Identification of m6A-Related lncRNAs
To identify mRNAs and lncRNAs, we developed a human configuration file and a gene expression matrix using PERL software, including the expression levels of m6A-related lncRNAs derived from the collated transcriptomic data. After running the script file, we were able to determine the nature of the genes and also the expression levels of mRNAs and lncRNAs. Gene IDs were converted into gene names using the Ensembl database (http://asia.ensembl.org/info/data/index.html). Using the limma package (http://bioconductor.riken.jp/packages/3.0/bioc/html/limma.html) in R software (https://www.r-project.org/), we extracted data on m6A-related gene expression, and also the corresponding gene names, based on gene types (i.e., readers, writers, and erasers). We then used co-expression analysis to examine the relationship between the expression of m6A-related genes and that of lncRNAs. The resulting correlation coefficients indicated regulatory hypotaxis between these elements. Using the limma package, we obtained expression data of m6A-related lncRNAs. Correlations were visualized in a network plot generated using the igraph package. Using the limma package, clinical survival data were cross-referenced with m6A-related lncRNA expression data. This analysis revealed lncRNAs related to prognosis, which were selected and utilized to calculate hazard ratios and confidence intervals using the survival package. A forest plot was generated to visualize univariate Cox regression coefficients. Finally, we computed and analyzed differences in m6A-related lncRNA expression between tumor tissues and normal tissues using the reshape2, pheatmap, limma, and ggpubr packages in R software. Differences were considered statistically significant at P <0.05. Heatmaps and boxplots were constructed to visualize differences in gene expression.
Exploring the Role and Function of m6A-Related lncRNAs
Based on lncRNA expression data (https://bioconductor.org/install/), we classified m6A-associated lncRNAs that were related to prognosis into clusters 1 and 2 using the Consensus Cluster-Plus and limma packages. In our calculations, the following parameters were used: clusterAlg = “km” and clusterNum = “2”. Furthermore, we used the survminer and survival packages to assess the prognostic value of m6A-related lncRNAs, and the pheatmap package to generate a heatmap enabling to identify differentially expressed lncRNAs related to prognosis in each cluster, and analyze their correlation with clinicopathological parameters. Using the limma package, we identified different types of tissues and differential target gene expression in the related OS subtypes. The common names of the target genes were obtained from the NCBI (https://www.ncbi.nlm.nih.gov). Finally, using the limma package, we conducted correlation analyses to clarify the relationship between the expression of prognostic m6A-associated lncRNAs and that of target genes in OS. Differences were considered statistically significant at P <0.05.
Analysis of the Tumor Microenvironment and Immune Cell Infiltration
We applied the preprocessor, limma, and e1071 packages to explore the levels of immune cell infiltration in the samples and calculated the number of infiltrating immune cells. We analyzed the TME using the limma package, and calculated its immune score, ESTIMATE score, and stromal score. These scores predicted tumor purity. Moreover, we estimated differences in the abundance of various immune cell subtypes between the clusters using the limma package, and visualized them using the vioplot package. The degree of infiltration of each immune cell type was also analyzed using the limma package, and the results were plotted on a boxplot. In parallel, we used the limma package to perform differential analysis of the TME based on the immune score, ESTIMATE score, and stromal score, and produced boxplots to determine the tumor cell purity of different OS subtypes. Next, gene set enrichment analysis (GSEA) was performed (http://www.gsea-msigdb.org/gsea/index.jsp) to identify pathways and related functions showing differential enrichment among samples, and PERL software was used to analyze the results. To understand whether the gene set enrichment of different clusters was dynamic or static, using the calculated scores and the generated plots, we classified each sample as “H” or “L” based on its assignation to a high- or low-risk lncRNA cluster, respectively. This analysis was performed by setting the ordering mode of the gene list as “descending”, the ranking gene metric as “Signal2Noise”, the sorting mode of the gene list as “real”, the normalization model as “meandiv”, and the permutation type and number to “no collapse” and “phenotype”, and “1000”, respectively (25). FDR <0.05 was set as threshold for statistical significance.
Prognostic Value of m6A−Related lncRNAs
We performed Lasso regression to generate a prognostic model. Considering the median risk score of prognostic m6A-associated lncRNAs, the samples were categorized as high- or low-risk. Lasso regression was carried out using 50% training samples and 50% test samples, and the corresponding plots were obtained. We then compared the high- and low-risk survival curves. To assess the accuracy of our model in predicting the survival rate of patients with OS, we generated a receiver operating characteristic (ROC) curve using the time ROC package. Data on the risk score, the risk associated with m6A-related lncRNAs, and survival status were combined to create a risk curve. Through independent prognostic factor analysis, we evaluated the independence of our model from other prognostic clinical factors. Multivariate and univariate analyses were performed, and the hazard ratios were calculated. We then validated our model using clinical data from various sources to assess its applicability to different clinical contexts, and conducted correlation analyses of clinical characteristics and risk scores to assess the relationship between these factors and m6A-associated lncRNAs. Using the limma and pheatmap packages, we generated a heatmap. Boxplots were constructed to assess the relationship between clinical data and risk scores. Using differential gene expression analysis, we evaluated differences in target gene expression between the risk groups, performed correlation analyses to determine the relationship between risk scores and levels of infiltrating immune cells, and visualized this nexus on a scatter plot to determine the beneficial role of immune cells.
Immunofluorescence
Primary antibodies against ZBTB32, DEPTOR, and SPAG4 were purchased from Abcam, Cell Signaling Technology, and Santa Cruz Biotechnology, respectively. Tissue samples embedded in paraffin were prepared for sectioning and antigen retrieval. According to the standard procedure, slides were incubated with primary antibodies overnight at 4°C. The slides were then washed three times with PBS, and incubated with secondary antibody at room temperature for 50 min in the dark. Subsequently, DAPI was used to stain nuclei. Finally, images were captured using a fluorescence microscope.
Results
Identification of m6A-Related lncRNAs
As shown in Figure 1, we created a detailed workflow for our entire study. First, the expression data of m6A-related genes were derived from the collated transcriptomic data by discriminating mRNAs from lncRNAs. Co-expression analysis showed a correlation between the expression of lncRNAs and that of m6A-related genes. To visualize this correlation, we generated a network plot (Figure 2A). Also, a forest plot was generated to visualize the results of univariate Cox regression analysis (Figure 2B). m6A-Associated lncRNAs were deemed to be prognostic at P <0.05, and DGCR5 (DCGR10) presented the highest hazard ratio. The heatmap and box plot in Figures 2C, D show differences in the expression levels of prognostic m6A-related lncRNAs between normal and tumor tissues. Among these, a total of 25 lncRNAs displayed differential expression between tumor tissues and normal tissues. The expression of some of these lncRNAs was lower in tumors, whereas that of others was higher (P <0.05).
Figure 1 Flow diagram of this study.
Figure 2 Differential expression of m6A-associated lncRNAs in relation to the prognosis of OS. (A) Network plot showing the relationship between the expression of lncRNAs and that of m6A-associated genes. (B) Forest plot showing the results of Cox regression analysis using univariate data. We extracted lncRNA-associated prognostic data and calculated hazard ratios and confidence intervals. Red indicates high risk, while green indicates low risk. (C) Boxplot showing differences in the expression levels of m6A-related lncRNAs associated with the prognosis of OS between normal tissues and tumor tissues. ***P < 0.001; **P < 0.01; *P < 0.05. (D) Heatmap showing differential expression of m6A-related lncRNAs in tumor tissues and normal tissues. ***P < 0.001; **P < 0.01; *P < 0.05. Red indicates high expression levels, while blue indicates low expression levels. Samples are reported on the abscissa and lncRNAs with prognostic value are reported on the ordinate.
The Role of m6A-Associated lncRNAs
The minimum overlap and lowest cumulative distribution function (CDF) value, in terms of lncRNA expression, occurred at K = 2. Therefore, lncRNAs were grouped into two groups (Figure 3 and Figures S1, S2), named clusters 1 and 2. The prognostic value of m6A-related lncRNAs was assessed based on survival analysis stratified by lncRNA subtype; this analysis revealed that lncRNAs from cluster 1 were associated with a lower survival rate than that of lncRNAs from cluster 2 (Figure 4; P = 0.007). Moreover, differences in the expression of three target plasma cell genes between the different types of tissue specimens, and also among the related subtypes, were investigated (Figure 5). The expression levels of ZBTB32 and DEPTOR were higher in normal tissues than in tumor tissues (Figure 5B, E; P <0.001), while the opposite was true for SPAG4 expression (Figure 5F, P <0.001). In addition, we performed gene correlation analysis to determine the relationship between OS-prognostic m6A-associated lncRNAs and target genes; this showed a correlation between the expression of DEPTOR and that of the prognostic m6A-related lncRNA AC036214.2; and between the expression of ZBTB32 and SPAG4 and that of the prognostic m6A-related lncRNA DGCR5 (DCGR10) (Figures 5C, F, I; P <0.05).
Figure 3 Classification of m6A-associated lncRNAs with prognostic value. The minimum overlap based on expression levels occurred at K = 2, together with the lowest cumulative distribution function (CDF) value; therefore, lncRNAs were classified into two clusters: cluster 1 and cluster 2. (A) Empirical CDF graph for K = 2–9. (B) Relative changes in the area under the CDF curve for K = 2–9. (C) Trace plot providing an overview of projected cluster members for different K values and of the cluster history relative to earlier clusters. (D) Consensus clustering matrix at K = 2.
Figure 4 Survival analysis and correlation of m6A-associated lncRNAs with clinical parameters. (A) Survival analysis revealed that, between the lncRNA subtypes, cluster 2 was associated with a higher 5-year survival rate (P = 0.007). (B) Heatmap showing differential expression of lncRNAs associated with prognosis and their correlation with clinicopathological parameters in the different clusters. Red indicates high expression levels, while blue indicates low expression levels. LncRNAs related to prognosis are reported on the ordinate and samples are reported on the abscissa.
Figure 5 Differential expression of target genes in different OS subtypes, and correlation analysis between m6A-associated lncRNA levels and target gene expression. Differences in the expression of ZBTB32, DEPTOR, and SPAG4 in various subtypes (A, D, G) and different tissue types (B, E, H). Correlation analysis of the relationship between the expression of the target genes ZBTB32, DEPTOR, and SPAG4 and that of prognostic m6A-related lncRNAs in OS tissues (C, F, I). Red indicates positive relationships, while blue indicates negative relationships. *P < 0.05, **P < 0.01, ***P < 0.001.
Correlation Between the Expression of ZBTB32, DEPTOR, and SPAG4 and That of Prognostic m6A-Related lncRNAs in OS
Immunofluorescence assays were carried out to verify the role of the three above-mentioned genes in OS (Figure 6). Compared to those of normal bone tissue (left), ZBTB32 (upper panel) and DEPTOR (middle panel) protein levels were consistently lower in OS tissues (right), while SPAG4 (bottom panel) protein levels were higher.
Figure 6 The expression of ZBTB32 and DEPTOR decreased while that of SPAG4 increased in OS tissues, as revealed by immunofluorescence staining of ZBTB32, DEPTOR, and SPAG4 in OS tissues and paired normal bone tissues.
Role of the TME and Infiltration of Immune Cells
TME analysis was performed to examine the degree of immune cell infiltration in the samples. The results of differential infiltration analysis in the two clusters were reported in a violin plot (Figure 7A). We also plotted the degree of immune cell infiltration within each OS cluster (Figure 7B and Figures S4–S6). These analyses revealed that cluster 2 exhibited a higher concentration of plasma cells than that of cluster 1 (P <0.05). We then examined the correlation between risk scores and immune cell infiltration levels in the TME of different OS subtypes, and generated boxplots to illustrate tumor cell purity in each of these subtypes (Figures 8A–C). Cluster 2 had a significantly higher stromal score than that of cluster 1, indicating a lower purity in the TME (P <0.05). We then conducted GSEA to identify pathways and related functions that showed differential enrichment between the samples (Figure 9 and Figure S3). A partial list of the GSEA results is reported in Figure 9. The gene set with the highest enrichment was NOD-like receptor signaling pathway. Cluster 2 was positively associated with all gene sets (P <0.05).
Figure 7 Comparison of immune cell infiltration in different clusters. (A) Violin plot. (B) Boxplot. Cluster 2 exhibited higher levels of immune cells, in particular of plasma cells (P <0.05).
Figure 8 Comparison of the tumor microenvironment (TME) in the two subtypes. Analysis of differences in the TME of different OS subtypes. The ESTIMATE score (A) and the immune score (B) did not differ significantly between the two types, whereas the stroma score (C) was higher in cluster 2 (P <0.05).
Figure 9 Differences in the enrichment of pathways and related functions between the two clusters based on gene set enrichment analysis (GSEA). Upper panel, partial results of GSEA. (A–I) Multiple pathways and functions related to cancer were found to be enriched in low-risk m 6 A-related lncRNAs, belonging to cluster 2. The latter is positively correlated to all gene sets (P <0.05).
Prognostic Value of m6A-Associated lncRNAs
To develop a prognostic model, we carried out Lasso regression and categorized all samples into a high-risk group and a low-risk group based on the median value of the risk score. Lasso regression was performed on a test group (50%) and a training group (50%) (Figures 10A, B). Patients were then divided into two subgroups, that is, a high-risk group and a low-risk group, based on their median risk score. Subsequent survival curve analysis revealed that patients exhibiting higher risk scores, from both the test and the training groups, had poorer prognoses. Conversely, patients in the low-risk group showed a higher survival rate than that of patients in the high-risk group (P <0.05; Figures 10C, D). Moreover, the analysis of ROC curves showed that a prognostic signature based on m6A-associated lncRNAs was a reliable predictor of test subtype [area under the curve (AUC), 0.929, 0.787, and 0.822 for 1, 3, and 5 years, respectively] and training subtype (AUC, 0.740, 0.787, and 0.855 for 1, 3, and 5 years, respectively; Figures 10E, F). Figure 11 illustrates the expression patterns (heatmap, upper panel), risk score (middle panel), and survival time (bottom panel) associated to m6A-related lncRNA prognostic regulators. Notably, the death rate and high-risk ratio increased with increasing risk scores. To examine the independence of our model from other clinical prognostic factors affecting patient outcomes, we performed independent prognostic factor analysis (Figure 12), which demonstrated that the risk score acts as an independent prognostic factor for OS (P <0.05). As illustrated in Figure 13, clinical validation of the model confirmed its suitability under different clinical parameters, including sex and age (P <0.05), suggesting that our model is broadly applicable. Furthermore, correlation analysis of clinical characteristics and risk scores was performed to evaluate how differences in the expression of m6A-related lncRNAs between low- and high-risk groups correlated with clinical characteristics (Figure 14A), and also with risk prognosis (Figures 14B–D). Notably, AC110015.1, SNHG6, and AC036214.2 showed high expression in the high-risk group, whereas AC073073.2, AC0086213.5, and BX322562.1 were poorly expressed in the low-risk group. Moreover, the clusters were closely related to the risk score (P <0.05). Finally, differences in the expression of target genes between the different risk groups in our OS model were observed (Figures 15A–C). In particular, ZBTB32 and DEPTOR expression levels were lower in the high-risk group than in the low-risk group (P <0.001), while the opposite was true for SPAG4 expression. A scatter plot was constructed to visualize this nexus and assign a beneficial or detrimental role to immune cells (Figure 15D). Plasma cells correlated negatively with risk scores (R <0 and P <0.01).
Figure 10 Prognostic model. (A, B) Results of Lasso regression for the construction of the prognostic model. (C, D) Survival curves related to the test group (C) and the training group (D) at P <0.05. (E, F) Receiver operating characteristic (ROC) curve indicating the accuracy of the model for survival prediction (left, test group; right, training group; area under the curve [AUC] >0.5).
Figure 11 Heatmap, risk-related curve, and spot plot of the test group (A, C, E) and the training group (B, D, F). The death rate and high-risk ratio increased with increasing risk scores.
Figure 12 Multivariate (A, C) and univariate (B, D) independent prognostic factor analyses in the test group (A, B) and the training group (C, D). The risk score was found to be an independent risk factor for OC prognosis (P <0.05).
Figure 13 Survival curves for model validation. The model resulted valid for various clinical groups, differing by age (A) or gender (B, C), at P <0.05.
Figure 14 Correlation between risk scores and clinical outcome. (A) Heatmap showing the correlation between risk scores and clinical outcome. (B–D) Boxplots showing the relationship between risk scores and clusters, gender, and immune scores. Clusters and risk score were closely correlated (P <0.05, ***P < 0.001).
Figure 15 Differential expression analysis of target genes and scatter plots showing the correlation of immune cell levels with risk scores. The expression of ZBTB32 (A) and DEPTOR (B) was lower in the high-risk group, whereas that of SPAG4 (C) was higher in the high-risk group (P <0.05). (D) Plasma cell levels were negatively correlated with the risk score (R <0, P <0.05).
Discussion
Previous studies have demonstrated that the downregulation or upregulation of N6-adenosine methylation regulators is related to the onset and development of various tumors, such as non-small cell lung cancer, acute myeloid leukemia, glioblastoma, breast cancer, and liver cancer (14, 26–30). Moreover, the same methylation regulators may function differently in different tumor tissues. Epigenetic modifications of proteins or DNA are crucial for cancer progression and prognosis (31, 32). Previous studies have focused primarily on m6A-related genes and pathways associated with cancer metastasis and prognosis; however, m6A-related lncRNAs have not been systematically analyzed in OS. Thus, the identification and analysis of lncRNAs associated with m6A are of considerable importance, as they can guide future OS research. Additionally, in this study N6-adenosine methylation of RNA molecules was also found to correlate with prognosis and levels of immune cell infiltration in the TME of OS tissues.
In the current study, we identified mRNAs and lncRNAs associated with m6A gene expression. Furthermore, co-expression analysis was performed to determine whether the expression of m6A-related genes was correlated with lncRNA expression. The resulting co-expression network plot revealed that the expression of several lncRNAs is correlated with that of m6A-related genes in OS. This finding sparked our interest in m6A-related lncRNAs and their role in OS; therefore, we identified lncRNAs associated with the prognosis of OS, calculated their hazard ratios and confidence intervals, and observed a close correlation between the expression of m6A-related lncRNAs and prognosis of OS, as revealed by Cox regression analysis of univariate data. Cancer development and m6A modifications are closely linked. Indeed, N6-adenosine methylation influences lncRNA splicing, thereby possibly altering cancer progression (33, 34). In our study, we found 25 lncRNAs associated with the prognosis of OS and m6A modifications, whose expression varied between tumor and normal tissues. A number of these lncRNAs were highly expressed in tumor tissues, while others were highly expressed in normal tissues (P <0.05).
Interestingly, PVT1, a gene significantly associated with OS survival, is upregulated in cancers through ALKBH5-mediated N6-adenosine methylation (35, 36). Recent studies have suggested that PVT1 may promote glycolysis, metastasis, and doxorubicin (DXR) resistance in OS cells (37, 38). LncRNA modifications can have a profound effect on the protein interaction network and can thus affect cancer progression (39). Altogether, the above findings revealed that a number of m6A-associated lncRNAs are highly expressed in tumor tissues, while others are highly expressed in normal tissues. Additionally, m6A-related lncRNAs can function as either tumor suppressors or oncogenes.
The role of m6A-associated lncRNAs was further analyzed. To evaluate the prognostic value of m6A-related lncRNAs, survival analysis based on lncRNA expression in different OS subtypes was conducted. The expression of lncRNAs associated with the low-risk group was found to be beneficial for the prognosis of OS. Furthermore, lncRNA expression was associated with OS survival rates. Similar results were reported by Chen et al., who demonstrated that the m6A demethylase ALKBH5 inhibits the binding of the reader protein YTHDF2 to the lncRNA PVT1, thereby suppressing its degradation (36). Notably, PVT1, whose expression is mediated by m6A, serves as a valuable prognostic predictor for patients with OS. Interestingly, m6A-related lncRNAs were expressed at similar levels in all clusters, possibly because of their low expression levels in OS. In addition, the study of the role of modified lncRNAs in OS progression is still limited. Therefore, further research on N6-adenosine methylation of lncRNAs and m6A recognition is urgently needed to validate our findings.
ZBTB32 (zinc finger and BTB domain containing 32), a gene located on chromosome 19, is widely expressed in the testis, lymph nodes, and other tissues. By targeting the Zpo2 to GATA3 promoters, ZBTB32 promotes the deregulation of GATA3 target genes, leading to the progression of aggressive breast cancer (40). For instance, human activated B cells in diffuse large B-cell lymphoma (DLBCL) carry ZBTB32 as the most consistently expressed gene (41). Conversely, CIITA, encoding an MHC class II transactivator, is silenced during the differentiation of B cells into plasma cells (42). However, ZBTB32 was identified as an early repressor of CIITA in an ex vivo plasma cell differentiation model. Therefore, ZBTB32 expression seems to primarily affect the maintenance of terminal plasma cells (43). In our study, ZBTB32 was found to be downregulated in OS samples, pointing at this gene as a potential tumor suppressor. However, higher expression of ZBTB32 was observed in cluster 2, possibly resulting in better prognosis of OC. Moreover, co-expression analysis of target genes revealed that ZBTB32 expression was significantly correlated with that of m6A-associated lncRNAs in OS tissues. In particular, ZBTB32 expression was highly correlated with DGCR5 (DCGR10) expression. As shown in Figure 2B, DGCR5 (DCGR10) expression was decreased in tissues from OS patients, resulting in decreased expression of ZBTB32. This observation further supports our hypothesis that DGCR5 (DCGR10) acts as a tumor suppressor.
In mature myeloma cells, DEPTOR (DEP domain containing MTOR interacting protein) increases sensitivity to therapeutic agents by inhibiting the activity of mTOR kinases (44, 45). Similarly, through the inhibition of mTORC1/2 signals, DEPTOR serves as a tumor suppressor in human prostate cancer (46, 47). In this study we observed low DEPTOR levels in OS patients, indicating that this gene might function as an oncosuppressor. In addition, DEPTOR expression was higher in cluster 2, suggesting its role in ensuring a favorable prognosis for OS. The results of our study concur with those of Hu et al., who found that downregulation of DEPTOR reduces the survival of OS patients via altered regulation of the PI3K/Akt/mTOR pathway (48). Correlation analysis between the expression of target genes and that of prognostic m6A-associated lncRNAs in OS have revealed that DEPTOR expression is closely associated with that of several m6A-related lncRNAs. Additionally, DEPTOR expression was negatively correlated with that of NBR2 and positively correlated with that of AC036214.2 and AL161785.1. Thus, these genes may serve as therapeutic targets for OS treatment.
SPAG4 (sperm associated antigen 4), a SUN family protein, was first discovered in the mammalian sperm tail. This protein interacts with outer dense fiber of sperm tails 1 (ODF1) and plays a crucial role in spermatogenesis and sperm motility (49). SPAG4 has been reported as a clinically relevant cancer marker (50). Recently, SPAG4 knockdown has been shown to reduce the invasion and migration abilities of cancer cells, while its overexpression enhanced these abilities of tumor cells (51, 52). In the present study, we found that the expression of SPAG4 was higher in patients with OS, suggesting that this gene might act as a tumor promoter. In addition, SPAG4 expression was low in cluster 2, suggesting that it may be detrimental to OS survival. Correlation analysis between the expression of SPAG4 and that of prognostic m6A-associated lncRNAs in OS revealed that SPAG4 expression was most negatively correlated with AC010609.1 expression.
Additionally, we measured the levels of infiltrating immune cells in order to determine the significance of immune cell infiltration in the TME of OS. Differential immune cell infiltration analysis revealed that cluster 2 was enriched in immune plasma cells. As previously stated, cluster 2 represents lower-risk OS cases characterized by good prognosis. Therefore, plasma cells infiltrating the TME may benefit the prognosis of the patient. This conclusion is consistent with the findings of Weiner et al., who showed that plasma cell enrichment correlated with prognosis and a fully functional regulation of immunity in human prostate cancer (53). Furthermore, it is becoming increasingly evident that plasma cells and intratumoral antibody production can significantly impact the antitumor response. In particular, various studies have shown that plasma cells are very strongly associated with long-term survival in many cancers (54–57). In this study, we also performed differential analysis of TME for various subtypes of cancer to examine the purity of tumor cells within different clusters. Notably, the scores of cluster 2 were consistently higher than those of cluster 1, suggesting a lower tumor cell purity and a greater number of immune and stromal cells within the TME. These results support the hypothesis that OS tissues belonging to cluster 2 are characterized by low risk scores and good prognosis. Similarly, multiple studies have demonstrated that the levels of tumor-infiltrating immune cells are closely associated with the prognosis of OS and are reliable predictors of OS outcome (22, 58–60). Altogether, our findings and previous studies confirm our hypothesis that immune cell infiltration affects patient prognosis. In particular, higher stromal scores result in less pure tumors and better prognosis than lower immune scores.
Recently, N6-adenosine methylation has been found to play a critical role in the tumor immune microenvironment and to be associated with cancer prognosis (61–65). In fact, N6-adenosine methylation regulates not only the development or function of immune cells and other stromal cells within the TME, but also the responses of the TME to various stimuli, such as hypoxia, internal or external cellular stress, metabolic dysregulation, chronic inflammation, and TGF-β (66, 67). Meanwhile, the TME also plays a key role in the complex regulatory network of N6-adenosine methylation, thereby affecting tumor initiation, progression, and response to therapy (67). Nevertheless, such regulatory relationship between the TME and N6-adenosine methylation and its role in tumor pathology require further study.
Next, we conducted GSEA. The most significantly enriched gene set was the NOD-like receptor signaling pathway. The function of this signaling pathway in OS is relatively unknown, and may thus be explored in future research to further improve the diagnosis and treatment of OS. We speculate that the expression of m6A-related lncRNAs may impact OS antitumor immunity by modulating the NOD-like receptor signaling pathway. Finally, we established a prognostic model based on m6A-related lncRNAs using Lasso regression. In both the test and the training group, the low-risk subtype displayed a higher survival rate than that of the high-risk subtype. Therefore, m6A-related lncRNA expression allows to accurately predict the outcome of patients with OS. In addition, this model for predicting disease survival has considerable accuracy. In fact, the number of deaths and the proportion of high-risk individuals increased as the risk score increased. Moreover, the model was not affected by other clinical prognostic factors that may affect the outcome of OS, and was found to be applicable to a variety of clinical groups. N6-adenosine methylation of lncRNAs may modulate their interactions with proteins and modify their structure, thereby repressing gene expression (60, 68). Moreover, such modification of lncRNAs may change the subcellular distribution of these molecules, thereby altering their stability and promoting metastasis and tumorigenesis (68).
Conclusion
In summary, data from the literature and the present study indicate that m6A-associated lncRNAs function as epigenetic regulators of plasma cell development. Meanwhile, the expression of m6A-related lncRNAs can serve as a clinical indicator for OC prognosis. In the developed OS prognostic model, differential gene expression analysis showed a higher expression of SPAG4 in the high-risk group, supporting the hypothesis that SPAG4 acts as a tumor promoter in OS. Although this function of SPAG4 has been shown in other cancers, only few studies have investigated the role of SPAG4 in OS cells; therefore, further studies are necessary to confirm these findings. Furthermore, the three examined target genes were found to regulate plasma cell differentiation. Correlation analysis between risk scores and levels of immune cells revealed a negative correlation between the risk score and plasma cell levels. Similarly, Weiner et al. proposed that the presence of plasma cells may improve tumor outcomes and thus represents a positive factor during cancer therapy (53). This hypothesis is consistent with our results of immune cell infiltration analysis stratified by cluster. In conclusion, our results suggest that epigenetic regulation of the development and differentiation of immune cells by m6A-associated lncRNAs yields promising predictive markers, and may provide valuable insights into OS treatment by guiding effective immunotherapy.
Data Availability Statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.
Ethics Statement
The studies involving human participants were reviewed and approved by The Third Xiangya Hospital. The patients/participants provided their written informed consent to participate in this study.
Author Contributions
ZW and XZ: Acquisition and analysis of data. YD, JW, and XW: Conceptual design. ZW, XZ, and XW: Wrote the manuscript. DC and ZL: Data analysis. All authors contributed to the article and approved the submitted version.
Funding
This work was supported by The Natural Science Foundation of China (81472058); Hunan Science and Technology Innovation Plan (2018SK2105; 422000008); The Fundamental Research Funds for the Central Universities of Central South University (under Grant No. 2020zzts896); Qinghai Key Laboratory of Laboratory Medicine, Qingke Fazheng [2019] No. 105.
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s Note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Acknowledgments
XW would like to thank his wife Ethna for her unwavering support of this study.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2021.806189/full#supplementary-material
Abbreviations
OS, Osteosarcoma; lncRNA, Long non-coding RNA; m6A, N6-methyladenosine; TME, Tumor microenvironment; GSEA, Gene set enrichment analysis; ROC, receiver operating characteristic; AUC, area under the curve; DEPTOR, DEP domain containing MTOR interacting protein; ZBTB32, zinc finger and BTB domain containing 32; SPAG4, sperm associated antigen 4.
References
1. Hogg SJ, Beavis PA, Dawson MA, Johnstone RW. Targeting the Epigenetic Regulation of Antitumour Immunity. Nat Rev Drug Discov (2020) 11:776–800. doi: 10.1038/s41573-020-0077-5
2. Timp W, Feinberg AP. Cancer as a Dysregulated Epigenome Allowing Cellular Growth Advantage at the Expense of the Host. Nat Rev Cancer (2013) 7:497–510. doi: 10.1038/nrc3486
3. Yi YC, Chen XY, Zhang J, Zhu JS. Novel Insights Into the Interplay Between M(6)A Modification and Noncoding RNAs in Cancer. Mol Cancer (2020) 1:121. doi: 10.1186/s12943-020-01233-2
4. Liu W, Zhao Z, Wang Y, Li W, Su Q, Jia Q, et al. Dioscin Inhibits Stem-Cell-Like Properties and Tumor Growth of Osteosarcoma Through Akt/Gsk3/β-Catenin Signaling Pathway. Cell Death Dis (2018) 3:343. doi: 10.1038/s41419-018-0363-x
5. Han G, Guo Q, Ma N, Bi W, Xu M, Jia J, et al. LncRNA BCRT1 Facilitates Osteosarcoma Progression via Regulating miR-1303/FGF7 Axis. Aging (Albany NY) (2021) 11:15501–10. doi: 10.18632/aging.203106
6. Buddingh EP, Kuijjer ML, Duim RA, Bürger H, Agelopoulos K, Myklebost O, et al. Tumor-Infiltrating Macrophages are Associated With Metastasis Suppression in High-Grade Osteosarcoma: A Rationale for Treatment With Macrophage Activating Agents. Clin Cancer Res (2011) 8:2110–9. doi: 10.1158/1078-0432.Ccr-10-2047
7. Siegel RL, Miller KD, Fuchs HE, Jemal A. Cancer Statistics, 2021. CA Cancer J Clin (2021) 1:7–33. doi: 10.3322/caac.21654
8. Wu Y, Zhou C, Yuan Q. Role of DNA and RNA N6-Adenine Methylation in Regulating Stem Cell Fate. Curr Stem Cell Res Ther (2018) 1:31–8. doi: 10.2174/1574888x12666170621125457
9. Yang Y, Hsu PJ, Chen YS, Yang YG. Dynamic Transcriptomic M(6)A Decoration: Writers, Erasers, Readers and Functions in RNA Metabolism. Cell Res (2018) 6:616–24. doi: 10.1038/s41422-018-0040-8
10. Heyn H, Esteller M. An Adenine Code for DNA: A Second Life for N6-Methyladenine. Cell (2015) 4:710–3. doi: 10.1016/j.cell.2015.04.021
11. Zhao BS, Roundtree IA, He C. Post-Transcriptional Gene Regulation by mRNA Modifications. Nat Rev Mol Cell Biol (2017) 1:31–42. doi: 10.1038/nrm.2016.132
12. Lence T, Akhtar J, Bayer M, Schmid K, Spindler L, Ho CH, et al. M(6)A Modulates Neuronal Functions and Sex Determination in Drosophila. Nature (2016) 7632:242–7. doi: 10.1038/nature20568
13. Du Y, Hou G, Zhang H, Dou J, He J, Guo Y, et al. SUMOylation of the M6a-RNA Methyltransferase METTL3 Modulates its Function. Nucleic Acids Res (2018) 10:5195–208. doi: 10.1093/nar/gky156
14. Visvanathan A, Patil V, Arora A, Hegde AS, Arivazhagan A, Santosh V, et al. Essential Role of METTL3-Mediated M(6)A Modification in Glioma Stem-Like Cells Maintenance and Radioresistance. Oncogene (2018) 4:522–33. doi: 10.1038/onc.2017.351
15. Sadikovic B, Yoshimoto M, Chilton-MacNeill S, Thorner P, Squire JA, Zielenska M. Identification of Interactive Networks of Gene Expression Associated With Osteosarcoma Oncogenesis by Integrated Molecular Profiling. Hum Mol Genet (2009) 11:1962–75. doi: 10.1093/hmg/ddp117
16. Liu L, Liu X, Dong Z, Li J, Yu Y, Chen X, et al. N6-Methyladenosine-Related Genomic Targets are Altered in Breast Cancer Tissue and Associated With Poor Survival. J Cancer (2019) 22:5447–59. doi: 10.7150/jca.35053
17. Chen M, Wei L, Law CT, Tsang FH, Shen J, Cheng CL, et al. RNA N6-Methyladenosine Methyltransferase-Like 3 Promotes Liver Cancer Progression Through YTHDF2-Dependent Posttranscriptional Silencing of SOCS2. Hepatology (2018) 6:2254–70. doi: 10.1002/hep.29683
18. Cui X, Wang Z, Li J, Zhu J, Ren Z, Zhang D, et al. Cross Talk Between RNA N6-Methyladenosine Methyltransferase-Like 3 and miR-186 Regulates Hepatoblastoma Progression Through Wnt/β-Catenin Signalling Pathway. Cell Prolif (2020) 3:e12768. doi: 10.1111/cpr.12768
19. Liu X, Liu L, Dong Z, Li J, Yu Y, Chen X, et al. Expression Patterns and Prognostic Value of M(6)A-Related Genes in Colorectal Cancer. Am J Transl Res (2019) 7:3972–91.
20. Li N, Kang Y, Wang L, Huff S, Tang R, Hui H, et al. ALKBH5 Regulates Anti-PD-1 Therapy Response by Modulating Lactate and Suppressive Immune Cell Accumulation in Tumor Microenvironment. Proc Natl Acad Sci USA (2020) 33:20159–70. doi: 10.1073/pnas.1918986117
21. Liu Z, Liu N, Huang Z, Wang W. METTL14 Overexpression Promotes Osteosarcoma Cell Apoptosis and Slows Tumor Progression via Caspase 3 Activation. Cancer Manag Res (2020) 12:12759–67. doi: 10.2147/cmar.S284273
22. Chen Y, Zhao B, Wang X. Tumor Infiltrating Immune Cells (TIICs) as a Biomarker for Prognosis Benefits in Patients With Osteosarcoma. BMC Cancer (2020) 1:1022. doi: 10.1186/s12885-020-07536-3
23. Sharma P, Allison JP. The Future of Immune Checkpoint Therapy. Science (2015) 6230:56–61. doi: 10.1126/science.aaa8172
24. Pagès F, Mlecnik B, Marliot F, Bindea G, Ou FS, Bifulco C, et al. International Validation of the Consensus Immunoscore for the Classification of Colon Cancer: A Prognostic and Accuracy Study. Lancet (2018) 10135:2128–39. doi: 10.1016/s0140-6736(18)30789-x
25. Yu ZL, Zhu ZM. N6-Methyladenosine Related Long Non-Coding RNAs and Immune Cell Infiltration in the Tumor Microenvironment of Gastric Cancer. Biol Proced Online (2021) 1:15. doi: 10.1186/s12575-021-00152-w
26. Du M, Zhang Y, Mao Y, Mou J, Zhao J, Xue Q, et al. MiR-33a Suppresses Proliferation of NSCLC Cells via Targeting METTL3 mRNA. Biochem Biophys Res Commun (2017) 4:582–9. doi: 10.1016/j.bbrc.2016.11.077
27. Yang X, Wang J. Precision Therapy for Acute Myeloid Leukemia. J Hematol Oncol (2018) 1:3. doi: 10.1186/s13045-017-0543-7
28. Cui Q, Shi H, Ye P, Li L, Qu Q, Sun G, et al. M(6)A RNA Methylation Regulates the Self-Renewal and Tumorigenesis of Glioblastoma Stem Cells. Cell Rep (2017) 11:2622–34. doi: 10.1016/j.celrep.2017.02.059
29. Zhang C, Samanta D, Lu H, Bullen JW, Zhang H, Chen I, et al. Hypoxia Induces the Breast Cancer Stem Cell Phenotype by HIF-Dependent and ALKBH5-Mediated M⁶A-Demethylation of NANOG mRNA. Proc Natl Acad Sci USA (2016) 14:E2047–56. doi: 10.1073/pnas.1602883113
30. Ma JZ, Yang F, Zhou CC, Liu F, Yuan JH, Wang F, et al. METTL14 Suppresses the Metastatic Potential of Hepatocellular Carcinoma by Modulating N(6) -Methyladenosine-Dependent Primary MicroRNA Processing. Hepatology (2017) 2:529–43. doi: 10.1002/hep.28885
31. Morrow JJ, Khanna C. Osteosarcoma Genetics and Epigenetics: Emerging Biology and Candidate Therapies. Crit Rev Oncog (2015) 3-4:173–97. doi: 10.1615/critrevoncog.2015013713
32. Lu B, He Y, He J, Wang L, Liu Z, Yang J, et al. Epigenetic Profiling Identifies LIF as a Super-Enhancer-Controlled Regulator of Stem Cell-Like Properties in Osteosarcoma. Mol Cancer Res (2020) 1:57–67. doi: 10.1158/1541-7786.Mcr-19-0470
33. Alarcón CR, Goodarzi H, Lee H, Liu X, Tavazoie S, Tavazoie SF. HNRNPA2B1 Is a Mediator of M(6)A-Dependent Nuclear RNA Processing Events. Cell (2015) 6:1299–308. doi: 10.1016/j.cell.2015.08.011
34. Huarte M. The Emerging Role of lncRNAs in Cancer. Nat Med (2015) 11:1253–61. doi: 10.1038/nm.3981
35. Yan M, Gao H, Lv Z, Liu Y, Zhao S, Gong W, et al. Circular RNA PVT1 Promotes Metastasis via Regulating of miR-526b/FOXC2 Signals in OS Cells. J Cell Mol Med (2020) 10:5593–604. doi: 10.1111/jcmm.15215
36. Chen S, Zhou L, Wang Y. ALKBH5-Mediated M(6)A Demethylation of lncRNA PVT1 Plays an Oncogenic Role in Osteosarcoma. Cancer Cell Int (2020) 20:34. doi: 10.1186/s12935-020-1105-6
37. Wan J, Liu Y, Long F, Tian J, Zhang C. Circpvt1 Promotes Osteosarcoma Glycolysis and Metastasis by Sponging miR-423-5p to Activate Wnt5a/Ror2 Signaling. Cancer Sci (2021) 5:1707–22. doi: 10.1111/cas.14787
38. Li D, Huang Y, Wang G. Circular RNA Circpvt1 Contributes to Doxorubicin (DXR) Resistance of Osteosarcoma Cells by Regulating TRIAP1 via miR-137. BioMed Res Int (2021) 2021. doi: 10.1155/2021/7463867
39. He RZ, Jiang J, Luo DX. The Functions of N6-Methyladenosine Modification in lncRNAs. Genes Dis (2020) 4:598–605. doi: 10.1016/j.gendis.2020.03.005
40. Shahi P, Wang CY, Lawson DA, Slorach EM, Lu A, Yu Y, et al. ZNF503/Zpo2 Drives Aggressive Breast Cancer Progression by Down-Regulation of GATA3 Expression. Proc Natl Acad Sci USA (2017) 12:3169–74. doi: 10.1073/pnas.1701690114
41. Care MA, Barrans S, Worrillow L, Jack A, Westhead DR, Tooze RM. A Microarray Platform-Independent Classification Tool for Cell of Origin Class Allows Comparative Analysis of Gene Expression in Diffuse Large B-Cell Lymphoma. PloS One (2013) 2:e55895. doi: 10.1371/journal.pone.0055895
42. Yoon HS, Scharer CD, Majumder P, Davis CW, Butler R, Zinzow-Kramer W, et al. ZBTB32 is an Early Repressor of the CIITA and MHC Class II Gene Expression During B Cell Differentiation to Plasma Cells. J Immunol (2012) 5:2393–403. doi: 10.4049/jimmunol.1103371
43. Zhu C, Chen G, Zhao Y, Gao XM, Wang J. Regulation of the Development and Function of B Cells by ZBTB Transcription Factors. Front Immunol (2018) 580:580. doi: 10.3389/fimmu.2018.00580
44. Peterson TR, Laplante M, Thoreen CC, Sancak Y, Kang SA, Kuehl WM, et al. DEPTOR is an mTOR Inhibitor Frequently Overexpressed in Multiple Myeloma Cells and Required for Their Survival. Cell (2009) 5:873–86. doi: 10.1016/j.cell.2009.03.046
45. Quwaider D, Corchete LA, Misiewicz-Krzeminska I, Sarasquete ME, Pérez JJ, Krzeminski P, et al. DEPTOR Maintains Plasma Cell Differentiation and Favorably Affects Prognosis in Multiple Myeloma. J Hematol Oncol (2017) 1:92. doi: 10.1186/s13045-017-0461-8
46. Chen X, Xiong X, Cui D, Yang F, Wei D, Li H, et al. DEPTOR is an In Vivo Tumor Suppressor That Inhibits Prostate Tumorigenesis via the Inactivation of Mtorc1/2 Signals. Oncogene (2020) 7:1557–71. doi: 10.1038/s41388-019-1085-y
47. Yao H, Tang H, Zhang Y, Zhang QF, Liu XY, Liu YT, et al. DEPTOR Inhibits Cell Proliferation and Confers Sensitivity to Dopamine Agonist in Pituitary Adenoma. Cancer Lett (2019) 459:135–44. doi: 10.1016/j.canlet.2019.05.043
48. Hu B, Lv X, Gao F, Chen S, Wang S, Qing X, et al. Downregulation of DEPTOR Inhibits the Proliferation, Migration, and Survival of Osteosarcoma Through PI3K/Akt/mTOR Pathway. Onco Targets Ther (2017) 10:4379–91. doi: 10.2147/ott.S143518
49. Tarnasky H, Gill D, Murthy S, Shao X, Demetrick DJ, van der Hoorn FA. A Novel Testis-Specific Gene, SPAG4, Whose Product Interacts Specifically With Outer Dense Fiber Protein ODF27, Maps to Human Chromosome 20q11.2. Cytogenet Cell Genet (1998) 1:65–7. doi: 10.1159/000014990
50. Kennedy C, Sebire K, de Kretser DM, O’Bryan MK. Human Sperm Associated Antigen 4 (SPAG4) Is a Potential Cancer Marker. Cell Tissue Res (2004) 2:279–83. doi: 10.1007/s00441-003-0821-2
51. Knaup KX, Monti J, Hackenbeck T, Jobst-Schwan T, Klanke B, Schietke RE, et al. Hypoxia Regulates the Sperm Associated Antigen 4 (SPAG4) via HIF, Which is Expressed in Renal Clear Cell Carcinoma and Promotes Migration and Invasion In Vitro. Mol Carcinog (2014) 12:970–8. doi: 10.1002/mc.22065
52. Ji Y, Jiang J, Huang L, Feng W, Zhang Z, Jin L, et al. Sperm−associated Antigen 4 (SPAG4) as a New Cancer Marker Interacts With Nesprin3 to Regulate Cell Migration in Lung Carcinoma. Oncol Rep (2018) 2:783–92. doi: 10.3892/or.2018.6473
53. Weiner AB, Vidotto T, Liu Y, Mendes AA, Salles DC, Faisal FA, et al. Plasma Cells are Enriched in Localized Prostate Cancer in Black Men and are Associated With Improved Outcomes. Nat Commun (2021) 1:935. doi: 10.1038/s41467-021-21245-w
54. Gentles AJ, Newman AM, Liu CL, Bratman SV, Feng W, Kim D, et al. The Prognostic Landscape of Genes and Infiltrating Immune Cells Across Human Cancers. Nat Med (2015) 8:938–45. doi: 10.1038/nm.3909
55. Schmidt M, Hellwig B, Hammad S, Othman A, Lohr M, Chen Z, et al. A Comprehensive Analysis of Human Gene Expression Profiles Identifies Stromal Immunoglobulin κ C as a Compatible Prognostic Marker in Human Solid Tumors. Clin Cancer Res (2012) 9:2695–703. doi: 10.1158/1078-0432.Ccr-11-2210
56. Lohr M, Edlund K, Botling J, Hammad S, Hellwig B, Othman A, et al. The Prognostic Relevance of Tumour-Infiltrating Plasma Cells and Immunoglobulin Kappa C Indicates an Important Role of the Humoral Immune Response in Non-Small Cell Lung Cancer. Cancer Lett (2013) 2:222–8. doi: 10.1016/j.canlet.2013.01.036
57. Tsou P, Katayama H, Ostrin EJ, Hanash SM. The Emerging Role of B Cells in Tumor Immunity. Cancer Res (2016) 19:5597–601. doi: 10.1158/0008-5472.Can-16-0431
58. Zhang C, Zheng JH, Lin ZH, Lv HY, Ye ZM, Chen YP, et al. Profiles of Immune Cell Infiltration and Immune-Related Genes in the Tumor Microenvironment of Osteosarcoma. Aging (Albany NY) (2020) 4:3486–501. doi: 10.18632/aging.102824
59. Hong W, Yuan H, Gu Y, Liu M, Ji Y, Huang Z, et al. Immune-Related Prognosis Biomarkers Associated With Osteosarcoma Microenvironment. Cancer Cell Int (2020) 20:83. doi: 10.1186/s12935-020-1165-7
60. Liu N, Zhou KI, Parisien M, Dai Q, Diatchenko L, Pan T. N6-Methyladenosine Alters RNA Structure to Regulate Binding of a Low-Complexity Protein. Nucleic Acids Res (2017) 10:6051–63. doi: 10.1093/nar/gkx141
61. Gu Y, Wu X, Zhang J, Fang Y, Pan Y, Shu Y, et al. The Evolving Landscape of N(6)-Methyladenosine Modification in the Tumor Microenvironment. Mol Ther (2021) 5:1703–15. doi: 10.1016/j.ymthe.2021.04.009
62. Zhao R, Li B, Zhang S, He Z, Pan Z, Guo Q, et al. The N(6)-Methyladenosine-Modified Pseudogene HSPA7 Correlates With the Tumor Microenvironment and Predicts the Response to Immune Checkpoint Therapy in Glioblastoma. Front Immunol (2021) 12:653711. doi: 10.3389/fimmu.2021.653711
63. Sun M, Xie M, Zhang T, Wang Y, Huang W, Xia L. m(6)A Methylation Modification Patterns and Tumor Microenvironment Infiltration Characterization in Pancreatic Cancer. Front Immunol (2021) 12:739768. doi: 10.3389/fimmu.2021.739768
64. Quan C, Belaydi O, Hu J, Li H, Yu A, Liu P, et al. N(6)-Methyladenosine in Cancer Immunotherapy: An Undervalued Therapeutic Target. Front Immunol (2021) 12:697026. doi: 10.3389/fimmu.2021.697026
65. Xu W, Tian X, Liu W, Anwaier A, Su J, Zhu W, et al. m(6)A Regulator-Mediated Methylation Modification Model Predicts Prognosis, Tumor Microenvironment Characterizations and Response to Immunotherapies of Clear Cell Renal Cell Carcinoma. Front Immunol (2021) 11:709579. doi: 10.3389/fonc.2021.709579
66. Li Y, Gu J, Xu F, Zhu Q, Chen Y, Ge D, et al. Molecular Characterization, Biological Function, Tumor Microenvironment Association and Clinical Significance of m6A Regulators in Lung Adenocarcinoma. Brief Bioinform (2021) 22(4):bbaa225. doi: 10.1093/bib/bbaa225
67. Liu Z, Zhong J, Zeng J, Duan X, Lu J, Sun X, et al. Characterization of the m6A-Associated Tumor Immune Microenvironment in Prostate Cancer to Aid Immunotherapy. Front Immunol (2021) 735170. doi: 10.3389/fimmu.2021.735170
Keywords: epigenetic, N6-methyladenosine (m6A), osteosarcoma, plasma cell, long non-coding RNAs (lncRNAs), tumor microenvironment
Citation: Wu Z, Zhang X, Chen D, Li Z, Wu X, Wang J and Deng Y (2022) N6-Methyladenosine-Related LncRNAs Are Potential Remodeling Indicators in the Tumor Microenvironment and Prognostic Markers in Osteosarcoma. Front. Immunol. 12:806189. doi: 10.3389/fimmu.2021.806189
Received: 31 October 2021; Accepted: 06 December 2021;
Published: 12 January 2022.
Edited by:
Guan-Jun Yang, Ningbo University, ChinaReviewed by:
Bin Wu, Wuhan Hospital of Traditional Chinese and Western Medicine, ChinaJessie Song, University of Macau, China
Copyright © 2022 Wu, Zhang, Chen, Li, Wu, Wang and Deng. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Xin Wu, NTIyODAxNTg1QHFxLmNvbQ==; Jianlong Wang, MTM3ODcwNzM3ODdAMTM5LmNvbQ==; Youwen Deng, ZHJ5d2RlbmdAMTYzLmNvbQ==
†These authors have contributed equally to this work