To obtain lncRNAs related to m6A regulation, we collected 23 m6A regulators from previous studies, including eight m6A writers, 13 readers and two erasers ( Table S1 ) (18, 20). Afterwards, by using a correlation analysis, 566 m6A-related lncRNAs (|cor value|>0.4, p<0.001) were screened from 14,086 lncRNAs. As shown in Figure 1A , hundreds of lncRNAs were suggested to be involved in m6A regulation, thus implying a potential role for m6A-related lncRNAs in glioma. Furthermore, consensus clustering was performed to clarify the expression characteristics of m6A-related lncRNAs in gliomas. The cumulative distribution function (CDF) of the consensus cluster from k=2 to 9 and the increment in the AUC were analysed. When the consensus matrix k was 2, there existed the least crossover between the glioma samples and the maximum increment of the area under the CDF curve; thus, k=2 was determined ( Figures 1B, C ). To visualize the expression pattern of m6A-related lncRNAs between the two clusters, we plotted a heatmap of 50 randomly selected m6A-related lncRNAs ( Figure 1D ), which exhibited a significant difference in the m6A-related lncRNAs between the two clusters. Moreover, as shown in Figure 1E , the prognosis in Cluster 2 was significantly better than that in Cluster 1. These results suggested the heterogeneous and prognostic values of m6A-related lncRNAs in gliomas.

Due to the fact that the importance of immune regulation in tumors has been proven (21–25), we profiled the immune landscape of m6A-related lncRNAs in gliomas by using the CIBERSORT and ESTIMATE algorithms. Fifteen of the 22 tumor-infiltrating immune cells differed in these two clusters ( Figure S1 ). Among them, immune cells with a proportion higher than 5% showed significant differences in the groups. As shown in Figure 2A , levels of monocytes, resting memory CD4 T cells and activated mast cells were higher in Cluster 1 (with a better OS) than in Cluster 2 (with a worse OS). In contrast, TAMs were significantly higher in Cluster 2 than in Cluster 1. In addition, the immune scores and stromal scores were significantly different between Cluster 1 and Cluster 2, as calculated by the ESTIMATE algorithm ( Figure 2A ). Cluster 2 had higher immune scores and higher stromal scores than Cluster 1. These results demonstrated glioma immune heterogeneity and suggested a correlation between m6A-related lncRNAs and the tumor immune microenvironment.

To investigate the potential biological processes involved in m6A-related lncRNAs, a gene set enrichment analysis (GSEA) was performed between Cluster 1 and Cluster 2. It revealed that 21 gene sets were significantly enriched in Cluster 2 compared with Cluster 1 (p<0.01) ( Supplementary 2 ). Among them, four hallmarks were significantly enriched in Cluster 2 (FDR<0.05), including the focal adhesion pathway, cytokine and cytokine receptor interaction pathway, apoptosis pathway and glycolysis pathway ( Figure 2B ). This result indicated that m6A-related lncRNAs may affect glioma progression by regulating specific pathways and further illustrated the significance of m6A-related lncRNAs in gliomas.

The above results demonstrated the vital role of m6A-related lncRNAs in glioma. Herein, we tried to develop a m6A-related lncRNA prognostic signature (MLPS) to identify key factors from the large number of m6A-related lncRNAs and to assess glioma patient prognosis accordingly. As shown in Table 1 , the TCGA dataset was randomly divided into the following two subsets: the training set (n=297) and the test set (n=283), with no significant difference in age, sex, WHO grade, IDH status or 1p/19q status. Next, 347 DE m6A-related lncRNAs (|log2FC|>1 and p value<0.01) were obtained by comparing gliomas in the training set to normal brain tissues from the GTEx dataset ( Figures 3 and 4A, B ). Furthermore, we performed a correlation analysis between the prognosis of patients and the expression of m6A-related lncRNAs in the training set and found 349 m6A-related prognostic lncRNAs. Subsequently, 73 m6A-related lncRNAs with ln(HR)>0 and log2FC>1 were determined to be promoters, whereas 71 m6A-related lncRNAs with ln(HR)<0 and log2FC<-1 were considered to be suppressors ( Figure 4C ). Afterwards, via the LASSO regression analysis, nine m6A-related lncRNAs were identified to construct the prognostic signature from 144 key m6A-related lncRNAs, including AL390755.1, AL445524.1, AL359643.3, LINC00641, AL117332.1, LNCTAM34A, CRNDE, AP001486.2 and CARD8. AS1 ( Figures 4D, E ). The coefficients of these lncRNAs were used to calculate the risk score.

A positive coefficient indicated that it was a risk factor, whereas a negative coefficient indicated a protective factor in glioma.

According to this risk score, we divided the samples of the training set into low-, mid- and high-risk groups, and the expression profiles of the lncRNAs constituting the signature were found to differ markedly ( Figure 5A ). Figures 5A, B shows that among these three groups, a higher risk score indicated a worsened prognosis of the patients (p<0.0001), including the survival rate and survival time. The univariate and multivariate Cox regression analyses demonstrated the independent prognostic value of the MLPS in the training set ( Figures 5C, D ). The ROC curves showed that the MLPS had a robust predictive ability, with AUCs of 0.881, 0.918 and 0.887 for 1-, 3- and 5-year survival, respectively ( Figure 5E ). Its predictive performance was better (AUC=0.928) than other classic indicators, including age, sex, WHO grade, IDH mutation status and 1p/19q codeletion status ( Figure 5F ). These results indicated that the prognostic signature generated by m6A-related lncRNAs may serve as an efficient and accurate indicator for evaluating prognoses and suggested the potential role of these lncRNAs in gliomas.

To confirm the reliability of the MLPS, we used the test set for validation. As shown in Figures 6A, B , the high-risk group had the worst survival rate and time, and the low-risk group had the best survival rate and time in the test set (p<0.0001), which was consistent with the results in the training set. The univariate and multivariate Cox regression analyses demonstrated the independent prognostic value of the MLPS in the test set ( Figures 6C, D ). The ROC curves showed that the AUCs were 0.856, 0.916 and 0.909 for 1-, 3- and 5-year survival, respectively, and the AUC of the risk score was higher than that of the other indicators ( Figures 6E, F ). These results showed that this MLPS did have an excellent predictive effect on the prognosis of glioma patients.

To further verify the reliability of MLPS, we collected other independent datasets (CGGA-seq-1 and CGGA-seq-2) for further validation. Two lncRNAs of the signature (CRNDE and LINC00641) were found in these independent datasets. The survival curve showed that the survival of patients with high CRNDE expression was significantly worse than that of patients with low CRNDE expression, which was observed both in primary glioma patients and recurrent glioma patients ( Figure S3A, B ). In contrast, LINC00641 exhibited the opposite results of CRNDE in glioma. The survival curve showed that patients with high LINC00641 expression had significantly better survival than patients with low CRNDE expression ( Figure S3C, D ). In addition, the expression of CRNDE was significantly higher in glioma patients with high grade, IDH wildtype, 1p/19q noncodeletion or older age. In comparison, the expression of LINC00641 was higher in glioma patients with lower grades, IDH mutation, 1p/19q codeletion and younger age ( Figure S2 ). These results suggested the clinical value and potential role of lncRNAs composing the signature in gliomas.

To explore the relationship between the m6A-related lncRNA signature and clinicopathological features of glioma, we analysed the risk score in different glioma subgroups ( Figure 7A ). The results showed that there was no significant difference in the score among White, Asian or Black and African patients ( Figure S3E ). In addition, the results showed that high-risk glioma patients had a significantly worse OS in both the younger subgroup (age ≤ 60) and the older subgroup (age>60). Likewise, high-risk patients had a significantly worse OS than low-risk patients in the male or female subgroups, IDH mutant or wild-type subgroups, lower-grade (WHO I-II) or higher-grade (WHO III-IV) subgroups and 1p/19q codeletion or noncodeletion subgroups ( Figure S3F ). These results indicated that this MLPS was effective in predicting the prognosis of glioma patients with different clinicopathological features.

In addition, we found that the high-risk group had significantly higher stromal scores and immune scores than the low-risk group ( Figure 7B ), which implied a potential correlation between the MLPS and the tumor microenvironment of glioma. We screened the immune cell types with more than a 5 percent proportion in all immune cells and analysed their correlations with the risk scores. Two types of immune cells had significant correlations with the risk scores (|R|>0.4) ( Figure 7C ). The correlation analyses showed that monocytes were negatively correlated with the risk score, whereas TAMs were positively correlated with the risk score. These results suggested a potential correlation between the m6A-related lncRNA signature and the distribution of specific immune cells. In recognition of the favourable predictive effect of this signature, there may be a regulation underlying m6A-related lncRNAs and immune cells, which could jointly affect the malignancy of glioma.

Among the nine factors of the MLPS, four lncRNAs (CRNDE, LINC00641, LNCTAM34a and CARD8-AS1) with defined sequences were selected for further in vitro experiments. We detected their expression levels in glioma cell lines (U87, LN229 and U343) and a normal human astrocyte line (SVGp12) by using RT–PCR. We found that the expression of LINC00641 in glioma cells was significantly lower than that in normal astrocytes, and the expression of CRNDE and LNCTAM34a was significantly higher in glioma cells ( Figure 8A ). These results supported the above analysis and verified the validity of the prognostic signature. However, the expression of CARD8-AS1 was inconsistent with the results from the analysis. Specifically, it was expressed at higher levels in normal cells than in glioma cells, which may be attributed to the small number of cell lines included in this study (results not shown).

Among the four markers mentioned above, previous studies have identified that CRNDE and CARD8-AS1 are risk factors, and LINC00641 is a protective factor for glioma, which is consistent with our study results. However, the role of LNCTAM34A in glioma is still unclear. Therefore, we applied a series of experiments to identify the function of LNCTAM34A in glioma. We knocked down LNCTAM34A by using specific siRNAs in U87-MG and LN229 cells and detected LNCTAM34A expression by using RT–PCR ( Figure 8B ). Subsequently, we tested cell proliferation, migration and epithelial-mesenchymal transition (EMT) in normal control glioma cells and LNCTAM34A knockdown glioma cells. The CCK8 assay showed that the proliferation of U87 and LN229 glioma cells was significantly reduced after LNCTAM34A knockdown ( Figures 8C, E ), and the Transwell assay showed that the migration of glioma cells was also significantly reduced after LNCTAM34A knockdown ( Figures 8D, F ), thus demonstrating that LNCTAM34A could promote glioma tumor proliferation and migration. A Western blot analysis showed that E-cadherin was upregulated, but Vimentin was downregulated, in LNCTAM34A knockdown cells, thus proving that LNCTAM34A could regulate the EMT process in glioma cells ( Figures 8G, H ). All of the results showed that LNCTAM34A was a risk factor for glioma, which was consistent with MLPS.

For the TCGA data, the RNA-Seq data of 698 gliomas, including glioblastoma (GBM) and lower grade glioma (LGG), as well as the corresponding clinicopathological data, were obtained from the TCGA website (https://portal.gdc.cancer.gov/). The TCGA dataset was randomly divided into the following two subsets: a training set (N=297) and a test set (N=283). This was done to construct and validate the prognostic signature in glioma, and the cases without clinical data matching were deleted. For the CGGA data, two subsets with RNA-Seq data and clinicopathological data, known as CGGA-seq-1 (N=280) and CGGA-seq-2 (N=657), were downloaded from the CGGA website (http://www.cgga.org.cn/). The RNA-Seq data of the normal tissue (N=423) were downloaded from the Genotype-Tissue Expression (GTEx) Project (http://gtexportal.org/home/) ( Table 2 ).

Based on previous studies, we extracted the expression profile of 23 m6A regulators from the TCGA database, including m6A writers (METTL3, METTL14, METTL16, WTAP, VIRMA, ZC3H13, RBM15 and RBM15B), readers (YTHDC1, YTHDC2, YTHDF1, YTHDF2, YTHDF3, HNRNPC, FMR1, LRPPRC, HNRNPA2B1, IGFBP1, IGFBP2, IGFBP3 and RBMX) and erasers (FTO and ALKBH5) (18, 20). Subsequently, a total of 14,086 lncRNAs were extracted from the TCGA dataset annotated by GENCODE (GRCh37). Moreover, we screened lncRNAs related to the 23 m6A regulators by using a Pearson correlation analysis and defined them as m6A-related lncRNAs. The correlation criteria were |cor value|>0.4 and p value<0.001. All of the analyses were performed and visualized by using R software.

All of the gliomas from the TCGA database were clustered by using the Consensus Cluster Plus package according to the expression of m6A-related lncRNAs. The cumulative distribution function (CDF) of the consensus cluster was from k=2 to 9. The Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis was performed by using GSEA software. The immune cell infiltration characteristics were analysed by using the CIBERSORT package, and we defined the total of M0, M1 and M2 as tumor-associated macrophages (TAMs). The immune scores, stromal scores and estimate scores of the samples were calculated by using the estimate package in R software.

We screened the m6A-related prognostic lncRNAs by merging the survival data with the m6A-related lncRNA expression data and by performing correlation analyses in the training set. Next, we screened the differentially expressed (DE) m6A-related lncRNAs by comparing gliomas in the training set with normal tissues from the GTEx dataset. Significantly differential expression was defined as |log₂FC|>1 and p<0.05. The m6A-related lncRNA with a high hazard ratio and positive fold change was defined as the tumor promotor. Conversely, the m6A-related lncRNA with a low hazard ratio and negative fold change was defined as a tumor suppressor. Seventy-three tumor promotors and 71 tumor suppressors were identified in glioma patients, and nine of them were used to construct the prognostic signature. These nine lncRNAs were selected via a LASSO regression analysis with a cut-off value of p<0.05. The following formula (based on a combination of the Cox coefficient and gene expression) was used to calculate the risk score.

βi represents the coefficients, and Si is the lncRNA expression level.

To evaluate the predictive ability of the MLPS for survival status, samples were classified into high-, mid- and low-risk groups by using the trisection points of all of the patients in the training set, and the corresponding risk scores were used to cut off the patients in the test set.

To test the correlation between the nine prognostic factor lncRNAs and the m6A regulators, we extracted the expression of these 9 lncRNAs and 23 m6A regulators in the TCGA and CGGA databases and detected their correlations by using the Pearson correlation analysis.

To stratify and validate the prognostic signature, samples were classified into high- or low-risk groups by using the medians of all of the risk scores. We further stratified patients into different subgroups according to age (age ≤ 60 or age>60), sex (male or female), WHO grade (lower-grade, including WHO I-II or higher-grade, including WHO III-IV), IDH status (mutant or wild type) and 1p/19q status (codeletion or noncodeletion). We classified patients aged ≤ 60 years as being younger and patients aged >60 years as being older in this study. The risk plot, survival curve and ROC curve were generated by using several packages in R software.

The human glioma cell lines U87-MG, LN229 and U343 were purchased from ATCC (Manassas, VA). Human normal SVGP12 astrocytes were purchased from Shanghai Institutes for Biological Sciences (Shanghai, China). Cells were maintained in DMEM (Invitrogen, Carlsbad, CA) supplemented with 10% foetal calf serum (Gibco BRL) and 1% penicillin plus streptomycin (HyClone, Logan, UT) and incubated in a humidified incubator (37°C, 5% CO2). LNCTAM34A siRNAs were purchased from JTS Scientific and transfected by using Lipofectamine 2000 (Invitrogen, USA) according to the manufacturer’s instructions.

Total RNA was extracted from cells by using TRIzol reagent (Invitrogen, Waltham, Massachusetts). Real-time polymerase chain reaction (RT–PCR) was performed by using the SYBR Green (Applied Biosystems, Foster City, CA) method with a CFX96 Real-Time PCR System (Bio–Rad, Hercules, California); GAPDH was used as the internal control. The settings for amplification were 95°C/120 s, followed by 39 cycles of 95°C/5 s and 60°C/30 s. GAPDH was used as an endogenous control, and the relative RNA expression was calculated by using the 2^−ΔΔCt method. Primers were generated by Sangon Biotechnology (Shanghai, China) (Table 3).

For this assay, 5,000 U87 and LN229 cells were seeded in a 96-well plate. After the cells had been cultured for 0, 12, 36 or 48 hours, cell proliferation was measured via the Cell Counting Kit-8 (CCK-8; Dojindo, Tokyo, Japan) assay, according to the manufacturer’s instructions.

The cell invasion ability was measured by using Transwell chambers (8-μm pore; BD Biosciences). U87 and LN229 cells (3× 10⁴ cells) suspended in 200 μl of serum-free culture medium were added to the upper chamber. The lower chamber was supplemented with DMEM supplemented with 10% FBS. After 24 h, the noninvading cells on the upper surface were separated, and the cells that had invaded to the bottom of the membrane were fixed with methanol and stained with 0.1% crystal violet. Digital image acquisition was performed after air drying. The number of invasive cells was counted by using a microscope.

U87-MG and LN229 cells were lysed by using RIPA lysis buffer (Thermo Scientific, Waltham, MA, USA). Total proteins were loaded into 10% SDS–PAGE gels and probed with E-cadherin (1:500, sc-7870, Santa), Vimentin (1:500, sc-32322, Santa) and β-Actin (1:5000, 3700s, CST), after which they were incubated with horseradish peroxidase-conjugated secondary anti-mouse IgG antibodies (1:2000; W4021, Beijing). The binding antibody was detected by using a hypersensitive ECL chemiluminescence kit (NCMBiotech), and images were collected by using a chemiluminescence imager (Image 800). Relative quantitative analysis was conducted based on the image band density ratio with ImageJ software (NIH, Bethesda, MD, USA).

A student’s t test, Pearson’s chi-square test and one-way ANOVA were used to compare different groups of data. Kaplan–Meier curves were used to evaluate the statistical significance of survival rates between different risk groups. The prediction accuracy of the risk characteristics was determined by using the ROC curve. Univariate and multivariate Cox regression analyses were performed to evaluate the significant prognostic factors. All of the statistical analyses were conducted by using R (ver. 5.0) and SPSS statistics (ver. 21) programs. P values<0.05 were considered to be statistically significant.