The prognostic model of low-grade glioma based on m6A-associated immune genes and functional study of FBXO4 in the tumor microenvironment

Data acquisition and differential analysis

Clinical information, mutation data, and agene expression dataset (HTSeq-FPKM) for LGG patients were obtained from the TCGA. To obtain an external validation set, we also downloaded datasets and clinical data for mRNAseq_693 and mRNAseq_325 from the CGGA website, and m6A-regulated genes were obtained from previous publications (Xu et al., 2020b). Immune-related genes were extracted from the gene set enrichment analysis (Subramanian et al., 2005) (GSEA: www.gsea-msigdb.org) website. All the samples were obtained from tissues obtained from primary glioma patients via clinical surgical excision. All analyses were restricted to WHO grade I–II cases, and high-grade tumors (grades III–IV) were excluded during the preprocessing stage. All the data were preprocessed via the sva and Limma packages (Ritchie et al., 2015).

To control for multiple comparisons in differential gene expression analysis, we applied the Benjamini-Hochberg false discovery rate (FDR) correction method using the p.adjust() function in R (method = ‘fdr’). This approach was chosen to balance the control of false positives while minimizing the loss of true positives, particularly given the large-scale nature of genomic data. Genes with an FDR-adjusted P-value < 0.05 and |log2 fold change| > 0.5 were considered statistically significant and included in subsequent analyses.

To obtain m6A-associated immune genes, we calculated correlations between immune genes and m6A-regulated genes and further selected the resulting results, which calculated an absolute correlation coefficient greater than 0.4 (P-value < 0.05). The Kaplan-Meier method and univariate Cox analysis were then used to perform survival analysis of the m6A-regulated genes, which in turn led to the desired hazard ratio (HR) value as well as its P value. A package (the Igraph package) with a correlation coefficient greater than 0.8 was used to plot the co expression network between these two genes. To obtain a heatmap that expresses the differences between the tumor and normal groups, we subsequently used the limma package, as others have done (Luo et al., 2022), to identify m6A-related immune genes (false discovery rate (FDR) < 0.05 and |log fold change (FC)| > 0.5).

Prognostic modeling, mutation analysis and dynamic nomograms

We divided all the LGG samples from the TCGA database into two groups: the training group and the internal validation group, while also seeking the corresponding external validation group in the CGGA database. To obtain statistically significant m6A-related immune genes for prognosis, we analyzed the training group correlation data via multivariate Cox stepwise regression analysis, univariate Cox regression, and the LASSO algorithm. In addition, highly correlated genes were excluded in order to optimize the prognostic model. In order to distinguish between the high-risk and low-risk groups, we determined the median risk scores of the training group and divided them accordingly. Next, we plotted receiver operating characteristic (ROC) curves and survival curves to test the predictive power of the model.

After constructing the prognostic model, we visualized and analyzed the data from the maftools package (Mayakonda et al., 2018) by producing mutation waterfall plots. We next analyzed the effect of mutational load on survival in both groups. To make our prognostic model more clinically useful, we downloaded the DynNom software package (Jalali et al., 2019) from the web and created a corresponding web-based dynamic histogram application, which helped us to determine the prognosis of our patients more accurately.

ssGSEA analysis and cluster typing

The samples from both TCGA and CGGA databases were analyzed via ssGSEA via the GSVA software package to determine the number of immune cells present in each sample. To distinguish between high and low expression of m6A-regulated genes, we performed a cluster analysis of the genes obtained separately via the ConsensusClusterPlus package (Wilkerson & Hayes, 2010). We also plotted the violin, heatmap, and survival curve of the tumor microenvironment for the above three groups to determine the correlation between m6A and immunity.

We then performed a cluster analysis of m6A-associated immune genes affecting prognosis (P < 0.05, top 50), via the same method, to obtain groups with high, medium, and low expression of immune genes. We then also plotted violin plots, heatmaps, and survival curves of the tumor microenvironment capable of expressing correlations between m6A-related immune genes and immunity. Finally, we mapped multigene enrichment curves by using the org.Hs.eg.db R package. Our results demonstrate significant enrichment of five Gene Ontology (GO) terms and pathways from the Kyoto Encyclopedia of Genes and Genomes (KEGG) in the high expression group of immune genes compared to the low expression group (P < 0.05).

m6A score calculation and sankey plot

We used the limma package to locate DEGs in the low-, medium-, and high-expression groups of m6A-regulated genes and then used their intersections via the VennDiagram package to obtain the final list of DEGs. The prognosis-related genes were then screened again via the one-way Cox method (P < 0.05). After obtaining information on the expression of prognosis-related genes, we calculated the m6A score for each sample via PCA. The resulting m6A score is then sorted into high and low groups via the surv_cutpoint function, and the corresponding survival curves are plotted.

We plotted Sankey diagrams to visualize the relationships among immune gene type, immune gene type, m6A score, and risk score. Correlation matrices for m6A score and immune cells were also generated via the corrplot package. Finally, we determined the correlation between them.

Comprehensive analysis of model genes

Furthermore, we also collated correlations between individual model genes and parameters such as stem cell indices, immune subtypes, the tumor microenvironment, and clinical stage. In addition, to validate the prognostic role of the model genes, the Gene Expression Profile Interaction Analysis (GEPIA) database was selected for analysis. Finally, to validate the accuracy of the prognostic model, we used a 5-year ROC curve, which enables the comparison of risk values, tumor inflammatory characteristic (TIS) scores, and TIDE scores and assesses the prognostic predictive efficacy between them.

Cell lines and siRNA transfection

The human astrocyte cell line HA1800 was obtained from the Shanghai Institutes for Biological Sciences (Chinese Academy of Sciences, Shanghai, China). The glioma cell lines U87-MG (HTB-14™), and U118-MG (HTB-15™) were obtained from the American Type Culture Collection (ATCC). U251-MG were purchased from Boster Biological Technology, Ltd. (Wuhan, China). The culture temperature was 37 °C, and the concentration of CO₂ used was 5%. Our culture conditions included Dulbecco’s modified Eagle medium (DMEM) (HyClone, Logan, UT, USA) supplemented with 10% fetal bovine serum (Gibco, Waltham, MA, USA). Our siRNA (100 nM) was purchased from RIBO Biotechnology, and after 24–48 h of incubation, the cells were transfected with the siRNA via an In vitro RNATMtransfection reagent (InvivoGen Biotechnology, Hong Kong, China).

qRT-PCR

Total RNA was extracted from the HA1800, U87, U118, and U251 cell lines, as well as siRNA-transfected U87 and U251 cell lines (including their respective blank controls), using TRIzol reagent (Invitrogen, Waltham, MA, USA). The RNA concentration and purity were measured using a spectrophotometer. Then, 1 µg o Next, quantitative PCR (qPCR) was performed using gene-specific primers (listed below) and SYBR Green reagent in a 96-well plate. The qPCR conditions were as follows: initial denaturation at 95 °C for 30 s, followed by 40 cycles consisting of denaturation at 95 °C for 10 s and extension at 60 °C for 30 s. β-actin was used as the internal reference gene, and relative expression levels were calculated using the ΔΔCt method. RNA was reverse-transcribed into cDNA. The primers used for manipulation were as follows:

FBXO4-Forward: GTCTTCGGGAAGACCATTGTTGG;

FBXO4-Reverse: GAGTTTCAGCCTCTGTATCCTGG;

β-actin-Forward: GACCACAGTCCATGCCATCA;

β-actin-Reverse: GTCAAAGGTGGAGGAGTGGG.

Cell proliferation

CCK-8 can be used to detect the proliferation of cells. We cultured U87 and U251 cells transfected with the FBXO4 siRNA for 48 h in a 96-well plate. In 96-well plates, 2,000 cells were cultured in each well with 100 μl of DMEM containing 10% FBS. We examined the value-added intensity of the treated cells every 24 h and obtained data for 5 days. In accordance with the instructions, we added equal amounts of CCK8 reagent (Yeasen, Shanghai, China) to each well of the plate and, after 1 h, absorbance values were obtained at 450 nm via a microplate reader (BioTek, Winooski, VT, USA).

Migration assay

Transwell assays were used to examine cell migration. For the transwell assay, 40,000 cells were inoculated into the upper chamber of a Matrigel transwell plate (BD Biosciences, Milpitas, CA, USA). After 24 h of incubation, the inserts were cleaned with 1× phosphate buffer solution, and the cells were fixed in 95% ethanol and finally stained with 1% crystal violet for 20 min at room temperature. Finally, we photographed the perforated cells in their microscopic state and then used ImageJ to calculate the number of cells that crossed the membrane.

Statistical analyses

We used R software (version 4.0.3) for routine statistical analyses. Some of the analysis data were also obtained from publicly available information bases and special software. Graph Pad (version 8.0.1) was also used for some of the graphs. The construction of the prognostic models relied on univariate Cox regression, multifactor Cox stepwise regression, and the LASSO algorithm. P < 0.05 was considered statistically significant.

Prognostic and immune co-expression of m6A genes

The clinical information we needed, derived from the TCGA database (N = 512), CGGA693 database (N = 443) and CGGA325 database (N = 186), has been collated in Table 1, and the results of survival analyses with respect to the m6A regulatory genes have been uploaded in File S1. Among the 22 m6A regulatory genes, except for METTL3, HNRNPC, LRPPRC, IGFBP1, and RBMX, all m6A regulatory genes were prognostic factors for LGG. The m6A prognostic network map clearly revealed a positive correlation between most m6A-regulated genes (Fig. 1A), with a few negative correlations indicated by blue lines. It also presents the P-values from the univariate Cox regression analysis, represented by the size of the dots, as well as whether each gene is a risk factor or a favorable factor, indicated by the color of the right half of the dots. In the co-expression network of immune genes and m6A-regulated genes, when the absolute value of the correlation coefficient is greater than 0.8, YTHDF2 exhibits the strongest correlation with most immune genes (Fig. 1B). Then, a correlation analysis was performed between m6A genes and immune genes, followed by batch correction of TCGA and CGGA samples, resulting in 7,757 m6A-associated immune genes. Subsequently, differential analysis was conducted between the normal group (N = 5) and the tumor group (N = 529), identifying 7,654 differentially expressed m6A-associated immune genes (FDR < 0.05 and |logFC| > 0.5). The differential expression results of these m6A-associated immune genes are visualized in a heatmap (Fig. S5).

Identification and validation of prognostic m6A-related immune genes: construction of a predictive model

We combined the differential expression profiles of m6A-related immune genes with survival time and survival status, and used the coxph function to evaluate the impact of these genes on survival time, identifying a total of 3,229 m6A-related immune genes associated with prognosis (P < 0.05). We then performed LASSO regression on the 3,229 genes, using penalty terms to shrink some coefficients to zero, narrowing the candidate genes down to nine (Figs. 2A, 2B). These nine genes were subsequently incorporated into a multivariate Cox regression model to assess their independent effects on survival. Ultimately, we identified four genes with a significant impact on prognosis. These four genes form the primary features of our prognostic model. The risk score formula for each sample is as follows: risk score = 0.157CRNDE + 0.160CUL7 + 0.505FBXO4 + 0.265TNFAIP6 (Fig. 2C, Table 2). Specific information on the training group, the internal validation group, and the external validation group is included in Files S2–S5, respectively. We calculated the median risk score from the training group and used this value to classify the samples in the training, internal validation, and two external validation groups into low-risk and high-risk groups. Survival analysis showed that, although survival varied across groups, the survival rate in the low-risk group was significantly higher than that in the high-risk group. We also performed AUC testing under the ROC curve in all four groups to validate the predictive efficacy of the model (Figs. 2D–2G). The results revealed an AUC value greater than 0.9 in the training group, an AUC ranging from 0.831 to 0.894 in the internal validation group, and an AUC between 0.623 and 0.813 in the external validation groups. These results suggest that our prognostic model has high predictive power in LGG.

The prognostic value of m6A-related immune genes and mutational load

We systematically evaluated whether riskScore has independent prognostic predictive ability for patients through univariate and multivariate Cox regression analyses. In the multivariate analysis, we further validated whether riskScore could independently predict the prognosis of patients after controlling for clinical variables such as age, sex, and grade. The results showed that riskScore had a significant P-value and reasonable HR in both the univariate Cox regression and multivariate analysis (Figs. 3A, 3B), indicating that riskScore can serve as an independent prognostic factor and help predict the survival of LGG patients. Mutation waterfall plots revealed the genes with the highest mutation frequencies, which were IDH1, TPP3, and ATRX (Figs. 3C, 3D). Further analysis revealed that mutational load had a significant effect on survival in both groups (Fig. 3E). To make it easier to calculate patient survival, we have designed a web-based dynamic histogram application (available at https://u20131050.shinyapps.io/LGG-m6A_ImmRNA-Dynamic_nomogram/) in which the expression levels of the model genes can be entered directly to obtain the desired results.

Sample clustering based on m6A regulatory genes and tumor TME scoring

In the clustering analysis, we determined the optimal number of clusters to be three based on the cumulative distribution function (CDF) plot, as the curve for K = 3 was tighter and smoother compared to K = 2. Based on the expression levels of m6A-regulating genes in TCGA and CGGA samples, the samples were divided into high, medium, and low expression groups. The specific classification and expression levels of each gene are shown in the heatmap (Fig. 4A). Survival analysis revealed a negative correlation between the expression of m6A-regulating genes and survival rate (Fig. 4B). Although no statistical difference was observed between the medium and low expression groups, significant statistical differences were found between the medium-high and high-low expression groups (P < 0.05). Through the analysis of heatmaps and violin plots of the tumor microenvironment, we observed that the high-expression group exhibited higher immune cell and stromal cell content, leading to enhanced immune activity and lower tumor purity. In contrast, the medium- and low-expression groups showed the opposite trend (Figs. 4C–4G).

Sample clustering based on prognostic m6A-related immune genes and tumor TME scoring

We grouped the TCGA samples into high, medium, and low expression categories based on the expression levels of prognostic m6A-related immune genes. The specific classification of samples can be visualized in the heatmap (Fig. 5A). Survival analysis revealed that the low immune gene expression group had a significantly higher survival rate than the other two groups (Fig. 5B). Further analysis showed that in the high-expression group, the levels of immune and stromal cells were higher, indicating increased immune activity and lower tumor purity. In contrast, the medium and low-expression groups showed lower levels of immune and stromal cell content, with higher tumor purity (Figs. 5C–5G).

In addition, GSEA enrichment analysis showed that the high-expression group of m6A-related immune genes was significantly enriched in pathways such as the negative regulation of alpha-beta T cell activation, negative regulation of CD4-positive alpha-beta T cell activation, negative regulation of macrophage migration, negative regulation of T cell receptor signaling pathway, and T cell tolerance induction in GO enrichment analysis (Fig. 6A, only the five most significant pathways are listed). In the KEGG enrichment analysis, this group was significantly enriched in the IL-17 signaling pathway, JAK-STAT signaling pathway, NF-κB signaling pathway, TNF signaling pathway, and Toll-like receptor signaling pathway (Fig. 6B, only the five most significant pathways are listed). Taken together, these results reveal an interesting phenomenon: in patients with high expression of prognostic m6A-related immune genes, multiple immune-suppressive pathways are activated, which may be one of the reasons for the lower survival rate observed in this group.

Prognostic significance of m6A scores and their association with immune cell infiltration in LGG

All information about the sample m6A scores can be found in File S6. As shown by the survival curves, the survival of the high m6A score group was significantly worse that of than the low m6A score group (Fig. 7A). The Sankey plots and box line plots revealed important correlations: the high m6A score group corresponded to the moderate immune gene expression group, the high immune gene expression group, and the high m6A regulatory gene expression group (Figs. 7B–7D). In addition, the correlation matrix clearly revealed significant positive correlations between m6A score and most immune cell types, although negative correlations were observed with CD56dim natural killer cells and monocytes. The above results revealed that the m6A score was strongly positively correlated with immune infiltration (Fig. 7E).

Association of model genes with clinical stage, immune subtype, tumor microenvironment, and stem cell index

The associations between the selected model genes and clinical stage, immune subtype, tumor microenvironment, and stem cell index were compared with the help of our constructed prognostic model. The results of the present study revealed that the immune subtype and clinical stage were significantly correlated with the selected model genes (Fig. S1). Additionally, stratification of both the tumor microenvironment and stem cell indices yielded significant correlations with model genes (Fig. S2). From the data extracted from the GEPIA database, we can conclude that all model genes are high-risk prognostic factors for LGG patients. The area under the 5-year ROC curve (AUC) values revealed that the efficacy of our prognostic model was significantly greater than that of the TIS and TIDE scores (Fig. S3).

Downregulation of FBXO4 inhibits the proliferation and migration of glioma cell lines

The biological role of model genes in low-grade gliomas needs to be further investigated. Using the Gene Expression Profile Interaction Analysis (GEPIA) server (Fig. S3C), we identified a gene, FBXO4, that was independently associated with overall survival in low-grade glioma patients. Our results revealed that, at the RNA level glioma cell lines expressed higher levels of FBXO4 than did normal human astrocytic brain cells, especially the U87 and U251 cell lines (Fig. 8A). Next, we silenced the FBXO4 target gene in the U87 and U251 cell lines (Figs. 8B, 8C). In addition, the siRNA target sequences we used (Si1, Si2, and Si3) have been uploaded (File S7). Next, we uniformly used the identified relatively efficient gene silencing sequences, i.e., si1, for the corresponding experimental studies. The CCK8 assay revealed that the proliferative capacity of glioma cells in the FBXO4-silenced group was significantly lower than that of the normal group (Figs. 8D, 8E). Next, in transwell assays, FBXO4 knockdown significantly inhibited the migration of U251 and U87 cells (Figs. 8F–8H). These results suggest that FBXO4 gene knockdown can inhibit the proliferation and migration ability of glioma cells.