The role of m6A RNA methylation regulator in meningioma

The expression level of m6A RNA methylation regulator in meningioma and normal tissues

The datasets (GSE55609 and GSE43290) were merged and processed in batches, and the 20 m6A regulatory genes were matched separately. Fifteen m6A regulatory genes were differentially expressed in meningioma and normal samples. Interestingly, the expression levels of ZC3H13, Y1HDC1, IGFBP3, YIHDF1, FTO, YTHDF2, YTHDF3, and HNRNPC in normal samples were substantially higher than those in tumor tissues (Figure 1A); the expression levels of the samples were also notably different between meningioma (Grade I and II) and normal tissues (Figure 1B, 1C). ZC3H13 and YTHDF3 were moderately correlated (COR = 0.59), and there were two negatively correlated genes (RAM15B and METTL3). Therefore, it was speculated that dysregulation of the m6A regulatory gene played a crucial role in the mechanism of meningiomas.

Using m6A RNA methylation regulator as a prognostic factor for meningiomas

To further verify the clinical effects of m6A methylation, a prognostic analysis of 18 m6A regulatory genes was conducted in the GSE16581 dataset. First, univariate and multivariate Cox regressions screened for prognosis-related and independent prognostic m6A regulators (Table 1). We believed that IGFBP2 and METTL3 were survival-related m6A regulatory genes with critical prognostic value in meningiomas.

Second, a risk-scoring system based on the level of gene expression and risk coefficient was constructed to comprehensively explore the relationship between the expression level of the m6A regulator and prognosis of meningioma. To avoid overfitting the model construction, LASSO regression was used, and a risk model for the two prognostic genes was constructed based on partial likelihood deviation, lambda value (lambda = 0.04) (Figure 2A, 2B), and RiskScore = (–0.34 * IGF2BP2) + (–1.36 * METTL3). Moreover, the “cutoff” package was used to identify the best cutoff value to group the prognostic models (high-risk and low-risk groups). Overall, univariate (C-index: 0.81) and multivariate Cox (C-index: 0.8) regression analyses further confirmed that the RiskScore was an independent prognostic factor for the prognosis of meningioma (Figure 2C, 2D). The high-risk group had a shorter survival time than the low-risk group (p < 5.3e-3) (Figure 3E). The receiver operating characteristic (ROC) curve was used to evaluate the prognostic model and prediction accuracy of the two m6A genes (Figure 2F). The prognostic RiskScore models showed 1-, 3-, and 5-year area under the ROC curve (AUC; 0.67, 0.85, and 0.79, respectively). RiskScore proved to be an independent prognostic factor of meningiomas with good sensitivity and specificity.

Construction of m6A RNA methylation gene network based on the WGCNA

The results revealed that RiskScore was an independent prognostic factor for meningiomas. The gene regulation of methylation in meningiomas was clarified by a designed gene regulatory network associated with the RiskScore; this was also used to establish separate groups (high- and low-risk groups) through the “cutoff” package, and WGCNA was performed. The soft threshold power was set to four (R² = 0.89) and modules were merged (threshold value = 0.2) to obtain 12 modules (Figure 3A–3D). The network heatmap indicated close relation between the modules (Figure 3D), and the eigengene adjacency heatmap showed that the modules were independent of each other (Figure 3E). The module-trait relationships heatmap (Figure 4A) helped in determining that the RiskScore was negatively correlated with the red module (COR = –0.3, p = 0.01).

To further understand the features of the characteristic genes in the red module, the correlation between module membership and significance of the gene was first analyzed (COR = 0.46, p = 2.6e-19) (Figure 4B and Supplementary Table 2); this accelerated the understanding of the regulatory mechanism of m6A methylation. Next, Gene Ontology (GO) function and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses were performed for all genes in the red module (Figure 5A, 5B). The red module GO (Biological Progress) analysis mainly comprised I-kappaB kinase/NF-kappaB signaling, and KEGG pathway analysis comprised necroptosis and natural killer cell-mediated cytotoxicity. The RiskScore regulatory network was identified, and the disease gene characteristics and potential biological mediation pathways of RiskScore as independent prognostic factors were preliminarily described.

m6A RNA methylation regulator gene network visualization

The Search Tool for the Retrieval of Interacting Genes/Proteins (STRING; https://cn.string-db.org/) database serves as an online platform that evaluates and integrates protein-protein interactions and transmits mutual information between organisms [10]. Here, all genes in the red module (n = 100) were analyzed. Moreover, Cytoscape software was used for network visualization and processing (Figure 6A). The m6A RNA methylation regulatory network landscape in meningiomas was also described. In other words, the mechanism of meningiomas could involve a gene network regulated by m6A RNA methylation and occur due to changes in the abovementioned genes. The results of this study lay a foundation for future research on the genetic landscape of meningiomas.

Low expression of m6A methylated regulators (METTL3 and IGF2BP2) in meningioma cells

To verify whether m6A methylated hub regulators (METTL3 and IGF2BP2) played an important role in meningioma, their mRNA expression levels in normal meningeal cells (HMC) and meningioma cells (IOMM-Lee cells) were first detected by reverse transcription polymerase chain reaction (Figure 7A, 7B). The results showed that the expression levels of METTL3 and IGF2BP2 in IOMM-Lee cells were lower than those in HMC. In addition, low expression of METTL3 and IGF2BP2 in meningioma cells was verified by western blotting (Figure 7A).

Bioinformatics analysis

The GEO database is a comprehensive and free public gene expression database. Here, three datasets, GSE16581 [26] (clinical information in Table 2), GSE55609, and GSE43290 [27], were included (Supplementary Table 1). Overall, 20 m6A methylation regulators were obtained from previous studies [28, 29], and the expression of GSE16581 was matched separately with the m6A methylation regulators.

Construction of the prognostic characteristics of methylation regulators

An expression matrix of 20 m6A methylation regulators was obtained from the GSE16581 dataset. First, univariate Cox regression analysis was used to identify m6A methylation regulators related to the prognosis of meningioma (p < 0.2 was used as the screening condition). Thereafter, multivariate Cox regression analysis was then used to assess independent prognostic factors associated with overall survival in patients with meningioma (with p < 0.05 as a screening condition). The LASSO Cox regression algorithm was used to construct the best prognostic prediction model for the m6A methylation regulators. The optimal cutoff value was defined by the expression level of the m6A methylation regulators (“cutoff” package, version 1.3). Patients with meningioma were divided into high- and low-risk groups based on the RiskScore value. The glmnet installation package in the R environment was used to perform LASSO regression analysis. A prognostic prediction model was constructed based on the m6A methylation regulators. Univariate (p < 0.05) and multivariate Cox (p < 0.05) regression analysis was used to evaluate the prognostic effect of RiskScore, and Kaplan–Meier curve analysis was performed using the log-rank test. The TimeROC package (version 0.4) was used to compare the predictive accuracy of the RiskScore of genetic characteristics of patients with meningioma, and p < 0.05 was considered to indicate statistical significance.

WGCNA

To better reveal the mechanism of m6A methylation in meningiomas, a gene co-expression network was established using the WGCNA package (version:1.71). This network identified genes closely related to the RiskScore and aimed at constructing a methylation gene regulatory network. First, the optimal cutoff value was determined through the “cutoff” package, and depending on the RiskScore value, high and low expression groups were established. All samples were clustered to identify outliers (height = 100), and PickSoftThreshold was used to determine the soft threshold (β = 4) to construct a scale-free topology model (R² > 0.85). To minimize the influence of noise, the adjacency matrix was converted into a topological overlap matrix (TOM). TOM-based dissimilarity was used to cut the module through a dynamic tree. The number of genes in the smallest module was 30 and cutting height was 0.25. The correlation analysis between the m6A regulator, color module, and genes in the module was visualized. The clusterProfiler installation package (version 4.3.3) was used to perform KEGG pathway enrichment and GO analyses of the genes in the RiskScore module.

Construction of the m6A gene network in meningioma

The gene in the red module was explored using the online PPI website (STRING); the network was visualized using Cytoscape_v3.8.2 and processed using cytoHubba. Finally, an m6A gene regulation network was constructed for meningiomas.

RNA extraction and real-time fluorescent quantitative PCR detection

Both meningioma and normal meningioma cells were purchased from the Cell Resource Center of the Shanghai Institute for Biological Sciences, Chinese Academy of Sciences. Total RNA was extracted from meningioma and normal meningioma cells using TRIzol reagent (Invitrogen, USA) according to the manufacturer’s instructions, and complimentary DNA was synthesized from 1 μg of total RNA using the PrimeScript μRT kit with gDNA eraser (Takara, Japan). The expression of METTL3 and IGF2BP2 mRNA was detected by quantitative reverse transcription PCR, using the SYBR premixed dimer eraser, and GAPDH was used as a normalization control. The following primers were used: METTL3 F: 5'-TCTTCAGCATCGGAACCAGC-3’, R: 5'-TGGGGATTTCCTTTGACACCAA-3’. IGF2BP2 F: 5'-GACAGGTCCTGCTGAAGTCC-3’, R: 5'-TGTTGACTTGTTCCACATTCTCC-3’. GAPDH F: 5′-CCCATCACCATTCCAGGAG-3′, R: 5′-GTTGTCATGGATGACCTTGGC-3′.

Western blot

Protein concentrations were estimated using the bicinchoninic acid protein assay. Total protein was extracted using cold RIPA buffer containing phenylmethanesulfonyl fluoride (Beyotime, China) and phosphatase inhibitor cocktail 2 (Sigma, USA) for 30 min on ice and at 12,000 × g at 4° C. The mixture was centrifuged for 20 min. Antibodies against METTL3 (Cat. No. 67733-1-Ig, 1:10000) and anti-IGF2BP2 (Cat. No. 11601-1-AP, 1:8000) were purchased from Proteintech (China). The GAPDH antibody (Cat. No. 11601-1-AP, 1:5000) was used as a normalization control.

Statistical analysis

R version 4.0.3 software (R Core Team (2020); R: Language and environment for statistical computing; R Foundation for Statistical Computing, Austria: https://www.R-project.org/) was used for data preprocessing, DEG screening, WGCNA, and functional enrichment analysis. One-way ANOVA followed by Dunnett’s multiple comparisons test was performed, using GraphPad Prism software version 8.0.0 for Windows (GraphPad Software, USA; https://www.graphpad.com), to analyze the expression of selected genes in different groups. The “pROC” R package (version 1.18.0) calculated the ROC curve and the AUC to evaluate the diagnostic value of RiskScore.

Data availability statement

The data used in this study were downloaded from the GEO database (https://www.ncbi.nlm.nih.gov/gds/?term), and the data used to support the results of this study were obtained from the corresponding authors.