N6-Methyladenosine-Related LncRNAs Are Potential Remodeling Indicators in the Tumor Microenvironment and Prognostic Markers in Osteosarcoma

Abstract

N6-Adenosine methylation, yielding N6-methyladenosine (m⁶A), 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 m⁶A modifications in the prognosis of OS is unclear. In the present study, we investigated the correlation between m⁶A-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 m⁶A 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 m⁶A-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 m⁶A levels were closely associated in co-expression analysis. The expression of m⁶A-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, m⁶A-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, m⁶A-related lncRNAs and infiltrating immune cells in the TME could serve as new therapeutic targets and prognostic biomarkers for OS.

Keywords: N6-methyladenosine (m6A); epigenetic; long non-coding RNAs (lncRNAs); osteosarcoma; plasma cell; tumor microenvironment.

Figure 1

Flow diagram of this study.

Figure 2

Differential expression of m⁶A-associated lncRNAs in relation to the prognosis of OS. (A) Network plot showing the relationship between the expression of lncRNAs and that of m⁶A-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 m⁶A-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 m⁶A-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.

Figure 3

Classification of m⁶A-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 m⁶A-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 m⁶A-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 m⁶A-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.

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.

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).

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).