The Prognostic Value and Immune Landscaps of m6A/m5C-related LncRNAs signature in the Low Grade Glioma

doi:10.21203/rs.3.rs-2747371/v1

Download PDF

Research Article

The Prognostic Value and Immune Landscaps of m6A/m5C-related LncRNAs signature in the Low Grade Glioma

https://doi.org/10.21203/rs.3.rs-2747371/v1

This work is licensed under a CC BY 4.0 License

Journal Publication

published 04 Jul, 2023

Read the published version in BMC Bioinformatics →

You are reading this latest preprint version

Background: N6-methyladenosine (m6A) and 5-methylcytosine (m5C) are the main RNA methylation modifications involved in the oncogenesis of cancer. However, it remains obscure whether m6A/m5C-related long non-coding RNAs (lncRNAs) affect the development and progression of low grade gliomas (LGG).
Methods: We summarized 926 LGG tumor samples with RNA-seq data and clinical information from the Cancer Genome Atlas and Chinese Glioma Genome Atlas. 105 normal brain samples with RNA-seq data from the Genotype Tissue Expression project were collected for control. 39 m6A/m5C-related genes were used to draw co-expressed lncRNAs network. We obtained a molecular classification cluster from the expression pattern of sreened lncRNAs. The least absolute shrinkage and selection operator (LASSO) Cox regression was employed to construct a m6A/m5C-related lncRNAs prognostic signature of LGG. In vitroexperiments were employed to validate the biological functions of lncRNAs in our risk model.

Results: The expression pattern of 14 sreenedhighly correlated lncRNAs could cluster samples into two groups, in which various clinicopathological features and the tumor immune microenvironment were significantly distinct. The survival time of cluster 1 was significantly reduced compared with cluster 2. This prognostic signature is based on 8 m6A/m5C-related lncRNAs (GDNF-AS1, HOXA-AS3, LINC00346, LINC00664, LINC00665, MIR155HG, NEAT1, RHPN1-AS1). Patients in the high-risk group harbored shorter survival times. Time-ROC curve analysis was confirmed that our model had promising predictive power in survival prediction for both the TCGA (AUC at 1/3/5years respectively: 0.86, 0.84 and 0.77) and CGGA cohorts (AUC at 1/3/5years respectively: 0.73, 0.76 and 0.76). Notably, WHO class II, IDH mutation, 1p/19q co-deletion and MGMT promoter methylation were decreased in the high-risk group. Immunity microenvironment analysis showed B cells, CD4+ T cells, macrophages, and myeloid-derived DC cells were significantly increased in the high-risk group. In addition, patients in high-risk group possessed more CNVs, and higher G-score. Patients in high-risk group had the worse overall survival time regardless of followed TMZ therapy or radiotherapy. All observed results from the TCGA-LGG cohort could be validated in CGGA cohort. Afterwards, we investigated the biological functions of LINC00664 in glioma cells in vitro. LINC00664 was found to promote cell viability, invasion and migration ability of glioma cells.

Conclusion: Our study elucidated a prognostic prediction model of LGG by 8 m6A/m5C methylated lncRNAs and a critical lncRNA regulation function involved in LGG progression. High-risk patients have shorter survival times and a pro-tumor immune microenvironment.

RNA methylation

long non-coding RNA

prognostic signature

Immune landscape

low grade glioma

Glioma is the most lethal primary intracranial neoplasm, with a reported overall survival of approximately 14–16 months [1]. The incidence rate of glioma is surmised at 3–4 per 100,000 population worldwide and continues to rise as human life expectancy prolongs [2]. In the case of a new diagnosis, the standard treatment consists of surgery followed by concurrent radiotherapy and temozolomide therapy with additional adjuvant temozolomide treatment. More than half of low grade gliomas(LGG) could progressively evolve into high grade gliomas. Despite recent advances in the treatment of glioma, overall survival remains poor and thus requires new treatment strategies [3]. The update of the WHO classification suggests that the role of molecular phenotyping is of enormous importance for the diagnosis and treatment of glioma [4]. Therefore, the development of molecular phenotyping prognostic risk models provides guidance for screening and discrimination of high-risk patients and might contribute to improving clinical outcomes.

RNA modifications have emerged as pivotal post-transcriptional regulators of gene expression programs, boasting over 170 different types of modifications [5]. N6-methyladenosine (m6A) is the most prevalent chemical epigenetic modification in mRNA post-transcriptional modifications, which is a dynamic and reversible methylation modification at the adenosine N6 position, comprising recognition, methylation, and demethylation [6]. Dynamic m6A modifications are involved in the formation and maintenance of cancer stem cells thereby manipulating cancer progression and treatment resistance [7]. The m6A demethylase ALKBH5 promotes the proliferation and invasion of glioblastoma stem cells by maintaining the expression level of FOXM1 [8]. 5-Methylcytosine (m5C) is a wide-ranging mRNA modification that acts on the untranslated regions (UTRs) of mRNA transcripts to exert regulatory effects including RNA export, ribosome assembly, and translation [9]. Loss of m5C methyltransferase NSUN5 in gliomas correlated with long-term survival of glioma patients [10]. Epigenetic RNA methylation modifications are inextricably linked to tumor proliferation and invasion

Long non-coding RNAs (lncRNAs) are multifunctional RNAs with structures similar to mRNAs, including poly-A tails and promoter structures, that regulate gene expression through either the editing of RNA, splicing of pre-mRNA, association with chromatin modifiers, or the abrogation of miRNA-induced repression [11]. A variety of lncRNAs have been studied for their engagement in glioma initiation and progression, including H19, Malat1, CRNDE, XIST, HOTAIR, and SOX2OT [12]. m6A is not only the most prevalent modification in mRNAs but also found in long non-coding RNAs [13]. The m6A reader YTHDF2 plays a crucial role in maintaining the stemness of glioma stem cells by stabilizing RNA expression in glioma stem cells [14]. LncRNA NKILA is regulated by the m5c methyltransferase NSUN2 to increase m5C levels, which promotes NKILA interaction with YBX1 for regulating the progression of cholangiocarcinoma [15]. Therefore, m6A and m5C methylated lncRNAs are robustly correlated with tumor development.

However, the exploration of the biological functions of m6A/m5C-related lncRNAs in low-grade gliomas remains scarce due to the limitations of our understanding of lncRNAs. In this article, we will focus on the biological significance of m6A and m5C modifications of lncRNAs in low-grade gliomas. By analyzing 926 glioma samples, we found that the expression of m6A/m5C-related lncRNAs signatures correlated significantly with glioma survival and immunity.

2.1 Data Collection and Processing

Access the TCGA (The Cancer Genome Atlas) and CGGA (Chinese Glioma Genome Atlas) data platforms to download and arrange the transcriptome profiles of LGG (https://portal.gdc.cancer.gov and http://www.cgga.org.cn/). The corresponding clinical information of glioma patients was also downloaded to match the RNA-seq data. Collate the data by ‘Annoprobe’ R package and filter the lncRNA expression data profile from it. The RNA-Seq data of 105 brain normal control tissues from the Genotype Tissue Expression (GTEx) Project were downloaded and collated from the UCSC XENA database (https://xenabrowser.net/). LGG somatic mutation data in TCGA were analyzed and visualized using the R package ‘maftools’[16].

2.2 mRNAs–lncRNAs co-expression network and differential expression level of m6A/m5C-related lncRNAs

Based on the previous report studies [5, 17], a list of 39 m6A/m5C-related genes was yielded (Table S1). Then, the finding of m6A/m5C genes correlated lncRNAs in TCGA and CGGA cohorts separately by spearman correlation analysis with a setting like absolute value of correlation coefficient > 0.5 and p value < 0.001. The univariate cox regression analysis was used to obtain m6A/m5C-related lncRNA harboring prognostic value from TCGA and CGGA datasets. Subsequently, the intersection of m6A/m5C methylation-related lncRNAs in TCGA and CGGA cohorts to obtain highly correlated lncRNAs. Then we generate gene co-expression network maps using Cytoscape software (Version: v3.9.0) to analyze the results of m6A/m5C spearman correlation analysis. The results of differential analysis of m6A/m5C-related lncRNA gene expression between the expression data of GTEx normal brain samples and LGG expression data in TCGA were presented in a heatmap.

2.3 Consensus clustering glioma cases by m6A/m5C -Related lncRNAs

The expression data of correlated lncRNAs were analyzed using the R package "ConsensusClusterPlus", in which optimal clusterNum was set to 2, resulting in 2 clusters. Subsequently, the distribution of the 2 clusters was visually displayed using the Principal Component Analysis (PCA) method.

2.4 Survival analysis

Survival analysis for the 2 categorical variables was performed by K-M curve analysis for comparison, with statistical differences defined as < 0.05. In another way, patients' survival differences were analyzed by survival status in a 2-category distribution.

2.5 Immune Landscapes Assessment

ESTIMATE algorithm, Timer algorithm, and Cibersort algorithm were used to evaluate immune scores and immune cell infiltration in 2 classified sample species [18, 19]. Frequently detected 29 immune checkpoint expression data were extracted from LGG gene expression profiles for differential expression comparison across clusters or risk groups [20].

2.6 Functional Enrichment Analysis

In order to explore differentially expressed genes (DEGs) of different clusters or risk groups, the R package “limma” based on the TCGA dataset was used with the criteria of | log2(FC)| > 1, p < 0.05, and FDR < 0.05. Gene enrichment analysis (GSEA) was performed to fathom the differentially functional enrichment pathways of samples in distinct clusters. DEGs in the different risk group were opted as candidates for Gene Ontology (GO) analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis based on the R package ‘clusterProfiler’ and ‘enrichplot’ for visualization.

2.7 Tumor-associated somatic mutations and copy number variations analysis

Tumor-associated somatic mutations analysis of 2-classified samples was performed utilizing the ‘maftools’ package. Copy number variation (CNV) evaluations were conducted across clusters or risk groups.

2.8 Prognostic Signature Constructed by m6A/m5C-Related lncRNAs

The matched clinical data were first collated. And then the 14 relevant lncRNAs screened were set as independent variables using least absolute shrinkage and selection operator (LASSO) regression analysis, where the method was set as Cox survival analysis to reduce the variables to fabricate a prognostic model. The respective risk score in each tumor sample was calculated by combining correlation coefficients and gene expression level of obtained m6A/m5C-related lncRNAs with such formula (Risk score =\({\sum }_{i=1}^{n}\text{C}\text{o}\text{e}\text{f}\text{i}\text{*}\text{E}\text{x}\text{p}\text{i}\), where Coefi means the correlation coefficients, Expi is the FPKM value of each m6A/m5C-related lncRNAs).

2.9 The predictive value and accuracy assessment of the risk model

Receiver operating characteristic (ROC) curve analysis, calibration curve, and time-ROC analysis was used to assess the predictive power and precision of the models. Univariate and multivariate Cox regression analyses were used to assess the value of risk models for clinical application.

2.10 In vitro validation assays

2.10.1 Cell transfection and siRNA plasmids.

The human glioma cell lines U87 and U251 were purchased from ATCC (Manassas, VA). Cells were maintained in DMEM (Gibco, USA) supplemented with 10% fetal bovine serum (Gibco, USA) and 1% penicillin plus streptomycin (Gibco, USA) and incubated in a humidified incubator (37°C, 5% CO2).

The sequences of siRNA were designed based on the NCBI Reference Sequence: NR_037194.1 using the siRNA Wizard software 3.1 (InvivoGen). Sequences of siRNAs were as follows. LINC00664 siRNA #1 GGTGATGACAGAATTGTAACA. siRNA #2 GTGACTGTTCTATTCATCATA. siRNA #3 GTTCAAATGGGACATATCTTA. Si-NC GAGGCGAATCGAATTAGATAT. LINC00664 siRNAs were transfected by using Lipofectamine 2000 (Invitrogen, USA) according to the manufacturer’s instructions.

2.10.2 Quantitative Real-Time PCR

Total RNA from U87 and U251 cells was extracted using a Trizol kit (Servicebio, China) and assessed for purity and concentration. The RNA was reverse transcribed into complementary DNA (cDNA). Quantitative RT-PCR was performed on a PCR machine (ABI Q1, CA, USA) with SYBR Green Master Mix (Yeasen, shanghai, China). Primers were adopted as follows: LINC00664, forward, 5′-TGCCTGTTCTCAGGGAAGAT-3′, reverse, 5′-CAGGCAGAGGACTCACATCA-3′. GAPDH, forward, 5′-AATGGGCAGCCGTTAGGAAA-3′, reverse, 5′-GCGCCCAATACGACCAAATC-3′.

2.10.3 CCK8 assay to determine cell proliferation

Firstly, glioma cells were evenly inoculated in 96-well cell plates. After glioma cells were attached to the wall, small interfering RNA to silence LINC00664 was transfected except for the control group. Subsequently, the cells were incubated with CCK8 reagent for 2 hours at 22h, 46h, and 70h after the transfection. Finally, the cell proliferation was determined by measuring the absorbance at 450 nm under the microplate reader (BD Biosciences, USA).

2.10.4 Trans-well assay to verify the migration ability of glioma cells

The culture apparatus used in the transfer well experiments was a transfer well migration champer (8um size; Corning, USA). After Si-LINC00664 for 24h, an equal number of each group of cells 5*10^4 were seeded in the upper chamber and DMEM medium with 30% serum concentration was placed in the lower chamber as an attractant. After 24h of incubation, the upper chamber was fixed, stained, and photographed under a microscope (OLYMPUS BX53).

2.10.5 Wound healing assay to assess migration capacity

Glioma cells transfected with GFP were uniformly grown in a 6-well plate and left to form a dense cell layer at the bottom. Then, vertical scribing with a 200ul microinjection tip to form a scratch. After taking pictures under the microscope, the cells were cultured with DMEM serum-free medium for 24h and then photographed again. Finally, the distance from the edge of the scratch to the growing edge of the cells was measured to assess the migration capacity of the cells under different conditions.

2.11 Statistical Analysis

All statistical analysis was performed on R Studio (version 4.1.1) and GraphPad Prism (version 9.0). Student’s t-test was performed for two-group comparisons. For comparisons among more than two groups, the Wilcoxon test and one-way ANOVA were utilised respectively for non-parametric and parametric data. P ≤ 0.05 was considered statistically significant.

3.1 Analysis of m6A/m5C highly correlated lncRNAs

The flow chart of this study was present in Figure S1. By Spearman analysis, adjusting the correlation coefficient to 0.5, we found 190 correlated lncRNAs with potential prognostic value in the TCGA cohort and 95 correlated lncRNAs in the CGGA cohort (Figure S2A). The TCGA and CGGA m6A/m5C-related lncRNA data sets were intersected to obtain 14 highly correlated lncRNAs that are representative of the majority of cases (Figure S2A, Table S2). Then, we incorporated the intersecting lncRNAs and related m6A/m5C genes to produce co-expression network maps by Cytoscape software (Fig. 1A). By joining the normal sample data in the GTEx database, we found that these 14 highly correlated lncRNAs were differentially expressed in healthy individuals and LGG cases (Fig. 1B).

3.2 Identification of two clusters of LGG cases by consensus clustering and survival analysis

The LGG cases in the TCGA cohort were classified into two clusters by clustering analysis with the best K = 2 consensus matrix (Fig. 1C). The PCA computation method was performed to display the characteristics of the two clusters from the TCGA cohort (Fig. 1D) or CGGA cohort (Figure S2B), which were found to be visually distinguishable between these samples. Subsequently, we performed survival analysis on the two clusters of cases and found that the survival time of cluster 1 was significantly reduced compared with cluster 2, both in the TCGA (Fig. 1E) and in the CGGA cohort (Figure S2C).

3.3 Analysis of oncological differences between the two clusters of LGG cases

Through heat map visualization, we showed the differences in expression levels of 14 highly relevant lncRNAs and the clinicopathological features between the two clustered samples (Fig. 2A). Interestingly, we found that cases in cluster 1 had fewer IDH mutations, 1p19q co-deletions, and MGMT promoter methylation both in the TCGA (Fig. 2A, Table S3) and in the CGGA cohort (Figure S2D, Table S5). In addition, it was found that the ESTIMATE score, immune score, and stromal score of cluster 1 were higher compared with cluster 2 (Fig. 2A). Then, the clusters were then analyzed for differences in immune checkpoint expression levels, where we found that immune checkpoints, including CD276, CD40, IDO1, CTLA4, LAG3, PDCD1, TNFSF14 and et al., were significantly upregulated in cluster 1 (Fig. 2B). Furthermore, the two clusters were remarkably different in immune infiltrating B cells, CD4⁺ T cells, and macrophages by immune cell infiltration analysis with Cibersort and Timer methods (Fig. 2C, D and Figure S2E, F). By GSEA enrichment analysis, it was found that cluster 2 was highly enriched in the Allograft rejection, asthma, GABAergic synapse, graft-versus-host disease, and intestinal immune network for IgA production pathways, while cluster 1 was enriched in morphine addiction, nicotine addiction, and synaptic vesicle cycle pathways (Fig. 2E). Mutation landscape analysis showed that the highly mutated genes both in clusters 1 and 2 were mainly IDH1, TP53, and ATRX, while cluster 1 had a higher mutated EGFR and PTEN gene compared with cluster 2 in the TCGA cohort (Fig. 2F).

3.4 Constructing and evaluating the risk model by m6A/m5C-related lncRNAs

Firstly, LASSO regression analysis was performed to reduce the variables for risk model construction with 14 lncRNAs narrowed down to 8 lncRNAs (Fig. 3A-C). These 8 lncRNAs were combined to build our risk model: GDNF-AS1, HOXA-AS3, LINC00346, LINC00664, LINC00665, MIR155HG, NEAT1, RHPN1-AS1. Subsequently, we grouped patients according to risk scores and found a significantly lower probability of survival among patients in the high-risk group from the TCGA cohort (Fig. 3D, Table S4) and CGGA cohort (Fig. 3G, Table S5). Next, a K-M curve was applied to analyze the survival differences between the different risk groups, verifying that patients in the high-risk group had shorter survival times (Fig. 3E, H). Furthermore, we used time-ROC curve analysis to confirm that our model had promising predictive power in survival prediction for both the TCGA (AUC at 1/3/5years respectively: 0.86, 0.84 and 0.77, Fig. 3F) and CGGA cohorts (AUC at 1/3/5years respectively: 0.73, 0.76 and 0.76, Fig. 3I). The risk scores were in agreement with previous clustering analyses and provided good discriminatory functions in terms of tumor molecular characteristics, such as fewer IDH1 mutations, 1p19q co-deletion, and MGMT methylation in the high-risk group (Fig. 3J, Figure S4A).

3.5 The value of risk models compared with other clinical factors in applications

Univariate Cox regression analysis revealed that our risk score was as favorable a prognosis judgment as other commonly used clinical prognostic indicators such as grade, IDH status, 1p19q co-deletion, and MGMT promoter methylation (Fig. 4A, B). Then, through multivariate COX regression analysis, we verified that the risk scores took a better prognostic predictive power compared with molecular characteristics in different cohorts (Fig. 4C, D). Using ROC curves to analyze the degree of predictive accuracy, we observed that the risk scores possessed the greatest AUC values among the various clinically relevant prognostic factors in the TCGA cohort as well as CGGA cohorts (Fig. 4E, F). Visualization of our multivariate COX model through Nomogram to visually depict that our risk score model have a greater weight in the evaluation system (Fig. 4G, Figure S3A). Moreover, our time-ROC curve analysis (Fig. 4H, Figure S3B) as well as the calibration curves (Fig. 4I, Figure S3C) demonstrated the high accuracy of our prediction model in predicting survival at different time durations.

3.6 Risk model distinguishes immune landscapes and oncological characteristics of different risk groups

The Sankey diagram showed the correspondence between our subgroups, and the risk scores for prognostic prediction yielded similar results to the consensus clustering analysis (Fig. 5A). Then, analyzing the expression levels of immune checkpoints, CD276, CD274, CTLA4, ICOS, ICOSLG, TNFSF14, and PDCD1 were observed to be elevated in the high-risk group (Fig. 5B). In addition, the stromal score, immune score, and ESTIMATE score were all higher in the high-risk group as measured by ESTIMATE immune rating (Fig. 5C, Figure S4B). With immune infiltration cell analysis, we found that B cells, CD4 + T cells, macrophages, and myeloid-derived DC cells were significantly increased among the high-risk group (Fig. 5D, E, Figure S4C, D). Via a further correlation analysis, we verified that neutrophils, myeloid DC cells, and M2 macrophages were all positively correlated with risk scores (Fig. 5F, Figure S4E). Combining the tumor TMB analysis, we found that patients with high TMB had shorter survival times (Fig. 5G). Additionally, risk scores remained well discriminating among patients in TMB loading groups, suggesting that patients with high-risk scores had shorter survival times regardless of whether they were classified in the high or low TMB group (Fig. 5H). Considering that tumor behavior is highly associated with mutations in genetic background, it was found by CNV analysis that the high-risk group possessed more CNVs, and a higher G-score (Fig. 5I).

In addition, the expression level of each screened lncRNA as our risk model could be an individual indicator to predict the prognosis of LGG patients (Figure S5A). Mutation landscape analysis from different risk groups showed that low-risk group had a relatively higher mutation rate in IDH1(low-risk group vs. high-risk group: 93% vs. 60%), while high- risk group had a higher mutated EGFR, NF1 and PTEN gene compared with low-risk group (Figure S5B). GO analysis revealed that cell functions enriched by heterogeneous genes between the high-risk group and its counterpart were as follows: regionalization, skeletal system development, pattern specification process, collagen fibril organization, anterior/posterior pattern specification (Figure S5C, Table S6). KEGG analysis unveiled that the signal pathways enriched by diverse genes between the high-risk group and its counterpart were: complement and coagulation cascades, ECM-receptor interaction, phagosome, focal adhesion and hematopoietic cell lineage (Figure S5D).

3.7 Predictive value of risk scores in clinical treatment subgroups

The routine treatment of glioma includes TMZ therapy and radiotherapy. In both the TMZ-treated and non-TMZ-treated subgroups, high-risk patients had shorter survival times in the TCGA cohort (Fig. 6A). Similarly, the same findings could be obtained in the CGGA group, and the survival differences between high- and low-risk patients were more pronounced (Fig. 6B). Patients in the high-risk group had shorter survival times than those in the low-risk group when grouped according to the presence or absence of radiotherapy both in TCGA and CGGA cohorts (Fig. 6C, D).

3.8 In vitro validation of biological functions of risk model lncRNAs

LINC00664 was sorted from risk model lncRNAs to validate its biological function in U87 and U251 cell lines. Firstly, we downregulated the expression level of LINC00664 in U87 and U251 cell lines by small interfering RNA and confirmed that it could promote the proliferation of glioma cells with the CCK8 proliferation assay (Fig. 7A, B). And then, it was demonstrated by a trans-well assay that the downregulation of LINC00664 could significantly reduce the invasion ability both of U87 and U251 cells (Fig. 7C). Subsequently, it was reconfirmed by wound healing assay that downregulation of LINC00664 reduced the migration ability of glioma cell lines U87 and U251 (Fig. 7D).

Ambulatory RNA methylation modifications, such as m6A and m5C, are involved in tumor proliferation, invasion, and immunological regulation [21]. As research progresses, there is growing evidence that lncRNAs in tumors could undergo RNA methylation and exert regulatory effects on tumors [22, 23]. A number of studies have shown that m6A and m5C methylated lncRNAs could be used as predictors of tumor survival [24, 25]. In our study, we investigated the two most common types of m6A/m5C methylated lncRNAs in LGG and constructed a predictive model for prognostic and immune characteristics. We identified that patients in the high-risk group had a shorter survival time and presented with a pro-tumor immune microenvironment including higher expression levels of immune checkpoints and suppressive immune cell infiltration. Assessment of m6A/m5C methylation-associated lncRNAs in LGG patients based on our risk model could provide innovative approaches to cluster therapy for glioma patients, while also providing a strategy for regulating the immune microenvironment in LGG.

LGG is primary brain tumors that usually occur in younger patients and have a higher long-term survival rate compared to high-grade gliomas. Therefore, a good prognostic prediction model helps us to predict and profile the survival status and tumor characteristics of patients, potentially facilitating the progress of LGG therapy and thus prolonging the survival time and quality of patients. Previously reported single-gene prediction methods no longer meet current needs, as single-gene models fail to characterize tumors adequately [26]. Zhang et al. reported a prognostic risk model for LGG based on six immune-related genes that could predict survival and immune cell infiltration [27]. Zheng et al. identified three clusters of m6A RNA methylation regulators by consensus clustering to predict the prognosis of low-grade gliomas [28]. Along with the improved understanding of functional tumor modifications and epigenetics, there is a need for a wider range of models to help us realize the relationship between LGG prognosis and functional modifications. We constructed predictive models by in-depth profiling of the two most common classes of functional lncRNAs with m6A/m5C methylation functional modifications to reveal the impact of methylation lncRNA profiles in low-grade gliomas on the prognosis.

lncRNAs are a class of functional RNA molecules that regulate the expression of protein-coding genes by recruiting or sequestrating gene regulatory proteins and microRNAs [29]. It was reported that LINC00689 regulated glioma cell growth, invasion, and glycolysis by competitively binding miR-338-3p to regulate the expression of the PKM2 gene [30]. Tang et al. reported that lncRNA TPTEP1 mediated MAPK signaling pathway through binding miR-106a-5p to regulate glioma stemness and radiotherapy resistance [31]. Therefore, lncRNAs are involved in the regulation of glioma from the perspective of tumor stemness maintenance, proliferation, invasion, treatment resistance, etc. In our in vitro validation, we demonstrated that one of the lncRNA LINC00664 participated in regulating the proliferation and invasion of glioma. Wang et al. reported that LINC00664 promoted the invasion, and proliferation of human oral squamous cell carcinoma via the miR-411-5p/KLF9 pathway, which was similar to our findings.

The immune microenvironment of glioma differs significantly from other tumors outside the central nervous system. The most common immune cells in gliomas are not T cells, but macrophages, which can account for more than 20% of the tumor components [32]. Therefore, macrophages are essential contributors to the regulation of the immune microenvironment in glioma and potentially convert tumor-infiltrating T cells from active to ineffective by modulating immune checkpoints [33]. In the present study, we found a higher level of upregulated immune checkpoint expression and M2 macrophage infiltration in our high-risk group, suggesting that our model allows for the initial screening of such patients to conduct macrophage-related immunotherapy. Catalina reported that myeloid-derived suppressor cells mediate immunosuppression by transferring PD-L1 into B cells, thereby promoting the development of a tumor-suppressive microenvironment [34]. The increase in both myeloid-derived cells and B cells found in our study implied that our model would be capable of identifying glioma patients with such characteristics.

The following aspects of our study still need further improvement: 1. Lack of external real-world cohort validation; 2. Our study validated only a single lncRNA biological function in the risk model, and more comprehensive and in-depth validation is needed; 3. The predictive value of m6A/m5C methylated lncRNAs for clinical applications needs further evaluation.

In this study, we constructed a prognostic prediction model of LGG by 8 m6A/m5C methylated lncRNAs and identified a critical lncRNA regulation function involved in LGG progression. High-risk patients have shorter survival times and a pro-tumor immune microenvironment.

Availability of data and materials

This study analyzed data from The Cancer Genome Atlas (TCGA) (https://portal.gdc.cancer.gov/), Chinese Glioma Genome Atlas (CGGA) ( http://www.cgga.org.cn/) and Genotype Tissue Expression (GTEx). These data are free and publicly available.

Acknowledgements

We thank the reviewers for their comments.

Funding

This research was supported by the Natural Science Foundation of China (82203879) and Natural Science Foundation of Hubei Province (2022CFB883) .

Author Contributions

Conceptualization, RL, CY, HH; Data Collection, RL, HC, CL, YQ; Methodology, RL, HC, YQ; Writing–original draft, RL, YQ, CY, HH; Writing–review and editing, KZ, JW, CY, HH. All authors have read and agreed to the published version of the manuscript.

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing Interests

The authors declare no competing interests.

McKinnon C, Nandhabalan M, Murray SA, Plaha P: Glioblastoma: clinical presentation, diagnosis, and management. Bmj 2021, 374:n1560.
Tan AC, Ashley DM, Lopez GY, Malinzak M, Friedman HS, Khasraw M: Management of glioblastoma: State of the art and future directions. CA: a cancer journal for clinicians 2020, 70(4):299-312.
Wen PY, Weller M, Lee EQ, Alexander BM, Barnholtz-Sloan JS, Barthel FP, Batchelor TT, Bindra RS, Chang SM, Chiocca EA et al: Glioblastoma in adults: a Society for Neuro-Oncology (SNO) and European Society of Neuro-Oncology (EANO) consensus review on current management and future directions. Neuro-oncology 2020, 22(8):1073-1113.
Louis DN, Perry A, Wesseling P, Brat DJ, Cree IA, Figarella-Branger D, Hawkins C, Ng HK, Pfister SM, Reifenberger G et al: The 2021 WHO Classification of Tumors of the Central Nervous System: a summary. Neuro-oncology 2021, 23(8):1231-1251.
Wiener D, Schwartz S: The epitranscriptome beyond m(6)A. Nature reviews Genetics 2021, 22(2):119-131.
Huang H, Weng H, Chen J: The Biogenesis and Precise Control of RNA m(6)A Methylation. Trends in genetics : TIG 2020, 36(1):44-52.
Paris J, Morgan M, Campos J, Spencer GJ, Shmakova A, Ivanova I, Mapperley C, Lawson H, Wotherspoon DA, Sepulveda C et al: Targeting the RNA m(6)A Reader YTHDF2 Selectively Compromises Cancer Stem Cells in Acute Myeloid Leukemia. Cell stem cell 2019, 25(1):137-148 e136.
Zhang S, Zhao BS, Zhou A, Lin K, Zheng S, Lu Z, Chen Y, Sulman EP, Xie K, Bogler O et al: m(6)A Demethylase ALKBH5 Maintains Tumorigenicity of Glioblastoma Stem-like Cells by Sustaining FOXM1 Expression and Cell Proliferation Program. Cancer cell 2017, 31(4):591-606 e596.
Haruehanroengra P, Zheng YY, Zhou Y, Huang Y, Sheng J: RNA modifications and cancer. RNA biology 2020, 17(11):1560-1575.
Janin M, Ortiz-Barahona V, de Moura MC, Martinez-Cardus A, Llinas-Arias P, Soler M, Nachmani D, Pelletier J, Schumann U, Calleja-Cervantes ME et al: Epigenetic loss of RNA-methyltransferase NSUN5 in glioma targets ribosomes to drive a stress adaptive translational program. Acta neuropathologica 2019, 138(6):1053-1074.
Kopp F, Mendell JT: Functional Classification and Experimental Dissection of Long Noncoding RNAs. Cell 2018, 172(3):393-407.
He X, Qi Y, Zhang X, Liu X, Li X, Li S, Wu Y, Zhang Q: Current landscape of tumor-derived exosomal ncRNAs in glioma progression, detection, and drug resistance. Cell death & disease 2021, 12(12):1145.
Huang H, Weng H, Chen J: m(6)A Modification in Coding and Non-coding RNAs: Roles and Therapeutic Implications in Cancer. Cancer cell 2020, 37(3):270-288.
Dixit D, Prager BC, Gimple RC, Poh HX, Wang Y, Wu Q, Qiu Z, Kidwell RL, Kim LJY, Xie Q et al: The RNA m6A Reader YTHDF2 Maintains Oncogene Expression and Is a Targetable Dependency in Glioblastoma Stem Cells. Cancer discovery 2021, 11(2):480-499.
Zheng H, Zhu M, Li W, Zhou Z, Wan X: m(5) C and m(6) A modification of long noncoding NKILA accelerates cholangiocarcinoma progression via the miR-582-3p-YAP1 axis. Liver international : official journal of the International Association for the Study of the Liver 2022, 42(5):1144-1157.
Mayakonda A, Lin DC, Assenov Y, Plass C, Koeffler HP: Maftools: efficient and comprehensive analysis of somatic variants in cancer. Genome research 2018, 28(11):1747-1756.
Li M, Zha X, Wang S: The role of N6-methyladenosine mRNA in the tumor microenvironment. Biochimica et biophysica acta Reviews on cancer 2021, 1875(2):188522.
Li T, Fu J, Zeng Z, Cohen D, Li J, Chen Q, Li B, Liu XS: TIMER2.0 for analysis of tumor-infiltrating immune cells. Nucleic acids research 2020, 48(W1):W509-W514.
Yoshihara K, Shahmoradgoli M, Martinez E, Vegesna R, Kim H, Torres-Garcia W, Trevino V, Shen H, Laird PW, Levine DA et al: Inferring tumour purity and stromal and immune cell admixture from expression data. Nature communications 2013, 4:2612.
Zhang C, Zhang Z, Zhang G, Zhang Z, Luo Y, Wang F, Wang S, Che Y, Zeng Q, Sun N et al: Clinical significance and inflammatory landscapes of a novel recurrence-associated immune signature in early-stage lung adenocarcinoma. Cancer letters 2020, 479:31-41.
Cui L, Ma R, Cai J, Guo C, Chen Z, Yao L, Wang Y, Fan R, Wang X, Shi Y: RNA modifications: importance in immune cell biology and related diseases. Signal transduction and targeted therapy 2022, 7(1):334.
Sun Z, Xue S, Zhang M, Xu H, Hu X, Chen S, Liu Y, Guo M, Cui H: Aberrant NSUN2-mediated m(5)C modification of H19 lncRNA is associated with poor differentiation of hepatocellular carcinoma. Oncogene 2020, 39(45):6906-6919.
Yang X, Zhang S, He C, Xue P, Zhang L, He Z, Zang L, Feng B, Sun J, Zheng M: METTL14 suppresses proliferation and metastasis of colorectal cancer by down-regulating oncogenic long non-coding RNA XIST. Molecular cancer 2020, 19(1):46.
Jin Y, Wang Z, He D, Zhu Y, Hu X, Gong L, Xiao M, Chen X, Cheng Y, Cao K: Analysis of m6A-Related Signatures in the Tumor Immune Microenvironment and Identification of Clinical Prognostic Regulators in Adrenocortical Carcinoma. Frontiers in immunology 2021, 12:637933.
Wang T, Kong S, Tao M, Ju S: The potential role of RNA N6-methyladenosine in Cancer progression. Molecular cancer 2020, 19(1):88.
Patel T, Bander ED, Venn RA, Powell T, Cederquist GY, Schaefer PM, Puchi LA, Akhmerov A, Ogilvie S, Reiner AS et al: The Role of Extent of Resection in IDH1 Wild-Type or Mutant Low-Grade Gliomas. Neurosurgery 2018, 82(6):808-814.
Zhang M, Wang X, Chen X, Zhang Q, Hong J: Novel Immune-Related Gene Signature for Risk Stratification and Prognosis of Survival in Lower-Grade Glioma. Frontiers in genetics 2020, 11:363.
Zheng J, Wang X, Qiu Y, Wang M, Yu H, Zhou Z, Wu Z, Jiang X: Identification of Critical m(6)A RNA Methylation Regulators with Prognostic Value in Lower-Grade Glioma. BioMed research international 2021, 2021:9959212.
McCabe EM, Rasmussen TP: lncRNA involvement in cancer stem cell function and epithelial-mesenchymal transitions. Seminars in cancer biology 2021, 75:38-48.
Liu X, Zhu Q, Guo Y, Xiao Z, Hu L, Xu Q: LncRNA LINC00689 promotes the growth, metastasis and glycolysis of glioma cells by targeting miR-338-3p/PKM2 axis. Biomedicine & pharmacotherapy = Biomedecine & pharmacotherapie 2019, 117:109069.
Tang T, Wang LX, Yang ML, Zhang RM: lncRNA TPTEP1 inhibits stemness and radioresistance of glioma through miR‑106a‑5p‑mediated P38 MAPK signaling. Molecular medicine reports 2020, 22(6):4857-4867.
Gutmann DH, Kettenmann H: Microglia/Brain Macrophages as Central Drivers of Brain Tumor Pathobiology. Neuron 2019, 104(3):442-449.
Buonfiglioli A, Hambardzumyan D: Macrophages and microglia: the cerberus of glioblastoma. Acta neuropathologica communications 2021, 9(1):54.
Lee-Chang C, Rashidi A, Miska J, Zhang P, Pituch KC, Hou D, Xiao T, Fischietti M, Kang SJ, Appin CL et al: Myeloid-Derived Suppressive Cells Promote B cell-Mediated Immunosuppression via Transfer of PD-L1 in Glioblastoma. Cancer immunology research 2019, 7(12):1928-1943.

No competing interests reported.

Download PDF

Journal Publication

published 04 Jul, 2023

Read the published version in BMC Bioinformatics →

Editorial decision: Major revision
23 May, 2023
Reviews received at journal
17 May, 2023
Reviews received at journal
16 May, 2023
Reviewers agreed at journal
08 May, 2023
Reviewers agreed at journal
06 May, 2023
Reviewers invited by journal
06 May, 2023
Editor assigned by journal
05 May, 2023
Editor invited by journal
11 Apr, 2023
Submission checks completed at journal
11 Apr, 2023
First submitted to journal
28 Mar, 2023

You are reading this latest preprint version

The Prognostic Value and Immune Landscaps of m6A/m5C-related LncRNAs signature in the Low Grade Glioma

Status:

Journal Publication

Version 1

Abstract

Figures

1. Introduction

2. Methods

2.1 Data Collection and Processing

2.2 mRNAs–lncRNAs co-expression network and differential expression level of m6A/m5C-related lncRNAs

2.3 Consensus clustering glioma cases by m6A/m5C -Related lncRNAs

2.4 Survival analysis

2.5 Immune Landscapes Assessment

2.6 Functional Enrichment Analysis

2.7 Tumor-associated somatic mutations and copy number variations analysis

2.8 Prognostic Signature Constructed by m6A/m5C-Related lncRNAs

2.9 The predictive value and accuracy assessment of the risk model

2.10 In vitro validation assays

2.10.1 Cell transfection and siRNA plasmids.

2.10.2 Quantitative Real-Time PCR

2.10.3 CCK8 assay to determine cell proliferation

2.10.4 Trans-well assay to verify the migration ability of glioma cells

2.10.5 Wound healing assay to assess migration capacity

3. Results

3.1 Analysis of m6A/m5C highly correlated lncRNAs

3.2 Identification of two clusters of LGG cases by consensus clustering and survival analysis

3.3 Analysis of oncological differences between the two clusters of LGG cases

3.4 Constructing and evaluating the risk model by m6A/m5C-related lncRNAs

3.5 The value of risk models compared with other clinical factors in applications

3.6 Risk model distinguishes immune landscapes and oncological characteristics of different risk groups

3.7 Predictive value of risk scores in clinical treatment subgroups

3.8 In vitro validation of biological functions of risk model lncRNAs

4. Discussion

5. Conclusion

Declarations

References

Additional Declarations

Supplementary Files

Status:

Journal Publication

Version 1