Raw RNA-seq data (FPKM files) and clinical data on GBM were extracted from the TCGA data portal (https://tcga-data.nci.nih.gov/tcga/: accessed September 2019). These data on GBM tissue samples (n = 169) and normal brain tissue samples (n = 5) were downloaded from the TCGA.

The thirteen m6A RNA methylation regulators (METTL3, METTL14, WTAP, KIAA1429, RBM15, ZC3H13, YTHDC1, YTHDC2, YTHDF1, YTHDF2, HNRNPC, FTO, and ALKBH5) were collated based on a previously published article (7–12). Then, we systematically compared the expression levels of these m6A RNA methylation regulators between GBM and normal brain tissues and used the Wilcoxon signed-rank test to perform a differential expression analysis. Heatmaps and violin charts were generated with R 3.6.1 (https://www.r-project.org/).

Interactions and Gene Ontology (GO) analysis of the m6A RNA methylation regulators were performed using the STRING database (http://www.string-db.org/). GEPIA (16) (http://gepia.cancerpku.cn/) is a web server for analyzing the RNA sequencing expression data of 9,736 tumor samples and 8,587 normal samples from the TCGA and the Genotype-Tissue Expression (GTEx) project using a standard processing pipeline. |log₂FC| is defined as median of the expression and the cut point of |log₂FC| is 1, the cutoff of p-value is 0.01.

The relation among the regulators was calculated by Pearson's correlation based on gene expression. To explore the function of m6A RNA methylation regulators in GBM, we clustered gliomas into different groups with “ConsensusClusterPlus” (50 iterations, resample rate of 80%, and Pearson's correlation http://www.bioconductor.org/). Then, principal component analysis (PCA) was performed with R 3.6.1 to examine gene expression patterns in different glioma groups.

To determine the prognostic value of m6A RNA methylation regulators, we performed univariate Cox regression analyses of their expression in the TCGA. Then, all regulators were selected for the functional analysis and development of a potential risk signature with the least absolute shrinkage and selection operator (LASSO) Cox regression algorithm. Finally, three genes and their coefficients were selected with the minimum criteria, choosing the best penalty parameter λ associated with the smallest 10-fold cross-validation within the training set. The risk score for the signature was evaluated by using the following formula:

where Coef_i and x_i are respectively the coefficient and the z-score-transformed relative expression value of each selected gene. This formula was used to calculate the risk score for each GBM patient in the TCGA dataset.

The 116 primary glioma samples and clinical information were acquired from postoperative patients, who undergone surgical treatment first time without chemotherapy or radiotherapy between 2013 June and 2014 January at Xiangya Hospital of Central South University. Eighteen normal brain tissues from patients with cerebral trauma surgery were collected as controls. The informed consents were provided to all patients or their family members. All the participants agreed the informed consents. All of samples were instantly frozen in liquid nitrogen and stored at −80°C until the extraction of total RNA. Because of multiple parameters like the lost to follow-up, only 61 of the patients experienced a follow-up period lasting 5 years since the surgery. This study was agreed by Central South University Xiang ya Hospital Medical Ethics Committee.

According to manufacturer's instructions, total RNA was extracted from glioma and normal samples using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) and cDNA was synthesized from 1 μg of total RNA using a Primescript™ RT reagent kit with gDNA eraser (TakaRa, Japan). RT-PCR was carried out by using SYBR Premix Dimer Eraser TM (Takara, Dalian, China) to detect the expression of mRNA HNRNPC, with GAPDH as a normalizing control. The following primers were used: HNRNPC F: 5′-CCTTACCATCAAACACGATGGC-3′, R: 5′-ACTTCGAAAAGATTGCCTCCACA-3′. GAPDH:F: 5′-CCCATCACCATCTTCCAGGAG-3′, R: 5′-GTTGTCATGGATGACCTTGGC-3′.

Protein concentrations were estimated using the BCA protein assay. Total proteins were extracted using cold RIPA buffer with Phenylmethanesulfonyl fluoride (Beyotime, China) and phosphatase Inhibitor Cocktail 2 (sigma, USA) for 30 min on ice, and centrifuged at 12,000 × g for 15 min at 4°C. Antibodies against HNRNPC were obtained from Proteintech. Antibody to β-actin was used as a normalizing control.

The expression of HNRNPC was detected by IHC in 25 cases of glioma tissues. The 5-μm-thick glioma tissue sections were dewaxed in xylene, rehydrated through decreasing concentrations of ethanol and washed in distilled water. According to the standard protocols, sections were processed and stained with hematoxylin and eosin (H&E) and diaminobezidine. At last, sections were dehydrated through increasing concentrations of ethanol and xylene to the transparent state, and sealed with neutral gum. Primary antibody was diluted at 1:100 according to the manufacturer's instructions. HNRNPC positive cells in microscopic fields at ×400 and ×100 were observed and the ICH results were assessed by two pathologists, respectively.

Wilcoxon rank-sum and signed-rank tests were texted to evaluate the expression levels of m6A RNA methylation regulators in GBM for different groups. Univariate and multivariate Cox regression analyses were tested to determine the prognostic value of the regulators. Based on the consensus expression of thirteen regulators, 169 patients with GBM were divided into two subgroups by using consensus cluster and divided into low-risk and high-risk groups by using the median risk score (originated from the risk signature) as a break point. Chi-square tests were performed to compare the characteristics of clinical features (n = 158) such as gender, race, age, IDH, and p53 mutant information, and survival time between the two risk groups and two clusters. The receiver operating characteristic (ROC) curves were used to show the prediction efficiency of the risk model for 5-year survival in clusters 1 and 2. The Kaplan–Meier method was used to compare the overall survival (OS) rates of patients in the high- and low-risk groups, while the log-rank test was tested to compare the survival distributions between two groups. The Cut-off point was used the mean of the expression of HNRNPC. All statistical analyses were conducted using R 3.6.1. Experimental data analysis was performed with the Graphpad Prism 8.0 (GraphPad, San Diego, CA, USA). The t-test was performed to detected differential expression of HNRNPC between glioma in different grades and normal brain tissues. The error bars in bar graphs were represented the mean ± SD (Standard Deviation). A p-value of < 0.05 was considered statistically significant.

To study the functions of the m6A RNA methylation regulators in the tumorigenesis and progression of GBM, differences between tumor and normal brain tissues were investigated and presented as a heatmap (Figure 1A, ^*p < 0.05, ^**p < 0.01, and ^***p < 0.001), showing that almost all the regulators were highly related to the oncogenesis and development of GBM. The differences between genes were presented as a violin plot (Figure 1B). The expression levels of METTL3, WTAP, KIAA1429, ZC3H13, YTHDC2, YTHDF1, YTHDF2, HNRNPC, FTO, and ALKBH5 were significantly different between GBM and normal tissues. Considering the small number of normal brain tissues in the TCGA, we investigated the expression of the m6A regulators on the online web server GEPIA, which contains large dataset of 9,736 tumors and 8,587 normal samples from the TCGA and the GTEx project, respectively (Figure 1C, ^*p < 0.05). A comparison of all results showed that the HNRNPC, WTAP, YTHDF2, and YTHDF1 were up-regulated in tumor samples and might influence the tumorigenesis of GBM.

Originated from the STRING online database (https://string-db.org/) and all m6A methylation regulators were filtered into the PPI network complex (Figure 2A). The correlations among the regulators were analyzed (Figure 2B). WTAP and FTO appeared to have a negative correlation, while METTL14 and YTHDC1 showed the most positive correlation. We also identified GO terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways and classified them into three functional categories: biological process (BP), cellular component (CC), and molecular function (MF). The top three enriched BP terms were “RNA modification,” “mRNA processing,” and “regulation of mRNA metabolic process.” For CCs, the top three terms were “RNA N6-methyladenosine methyltransferase complex,” “nuclear speck,” and “nucleoplasm.” The top three enriched MF terms were “N6-methyladenosine-containing RNA binding,” “RNA binding,” and “oxidative RNA demethylase activity.” All the results were shown in Figure 2C. The pathway enrichment analysis showed that the genes were associated with the reversal of alkylation damage by DNA dioxygenases and processing of the capped intron-containing pre-mRNA signaling pathway (Supplementary Table 1). These findings showed that the regulators had a broad connection to the RNA modification and cancer processes.

Considering the dramatic imbalance in the numbers of GBM (n = 169) and normal brain tissue samples (n = 5) and the in-depth knowledge of GBM, we identified two new clusters of GBM based on the expression of all of the genes related to m6A RNA methylation (k = 2, which appeared to fit with the selection based on clustering stability increasing from k = 2 to 9 in the TCGA dataset) (Figures 3A–C). However, 159 of 169 gliomas clustered into one of the two subgroups in the TCGA dataset. Different clinical characteristics and the expression of m6A regulators between the two clusters (Figure 3D) showed that the different clusters are related to the survival status but the age, race, gender, IDH, and p53 mutant and HNRNPC, ALKBH5, WTAP, YTHDF2, YTHDC2, and FTO were markedly different between the two groups (Figure 3F). HNRNPC, WTAP, and YTHDF2 were upregulated in cluster 2, while ALKBH5, YTHDC2, and FTO were down-regulated in cluster 2. Principal component analysis (PCA) was performed to compare the transcriptional profiles between the cluster 1 (n = 105) and cluster 2 (n = 54) subgroups. There was a clear distinction between them (Figure 3E). Then, we compared the clinical survival outcomes of these two subgroups (cluster 1 and 2) clustered by k = 2. Interestingly, cluster 2 seemed to have a good survival trend (Figure 3G). Therefore, we believed that HNRNPC, WTAP, and YTHDF2 might be predictors of a high survival rate.

We identified genes that were significantly changed (|log₂-fold-change| > 1 and normalized p < 0.05) in the cluster 2 subgroup and then annotated their function using GO and pathway analyses. The results indicated that the downregulated genes in cluster2 are enriched in BPs including “chemical homeostasis,” “transmembrane transport,” and “ion transport” (Figure 4A), CCs including “intrinsic component of membrane,” “integral component of membrane,” and “plasma membrane” (Figure 4B) and MFs including “transmembrane signaling receptor activity,” “signaling receptor activity,” and “transmembrane transporter activity” (Figure 4C). Additionally, genes involved in KEGG pathways were enriched in “cGMP-PKG signaling pathway,” “neuroactive ligand-receptor interaction,” and “calcium signaling pathway” (Figure 4D). The upregulated genes were enriched in intracellular protein transport and ribosome (Supplementary Figure 1). All these findings indicate that the cluster groups divided by the expression of regulators are associated with the development of cancer.

Considering the prognostic indices, we sought to investigate the prognostic role of m6A RNA methylation regulators in GBM. We performed univariate and multivariate Cox regression analysis on the expression levels in the TCGA dataset (Figures 5D,E). The results indicated that HNRNPC might be a protective gene that is significantly related to OS, with a hazard ratio (HR) < 1 (p < 0.05). The remaining tested genes were not (p > 0.05). To better predict the clinical outcomes of GBM with m6A RNA methylation regulators, the LASSO Cox regression algorithm was applied to thirteen genes. Three genes (HNRNPC, ZC3H13, and YTHDF2) were used to build a risk signature based on minimum criteria, and the coefficients obtained from the LASSO algorithm were used to calculate the risk score for the TCGA dataset (Figures 5A,B). We separated patients with GBM from the TCGA (n = 159) dataset into low-risk (n = 80) and high-risk (n = 79) groups and observed that GBM patients in the low-risk group had significantly longer OS than those in the high-risk group (p = 1.076e−02) (Figure 5C), showing that the model could predict prognosis well.

To study the prognostic value of the model. Figure 5F showed the distributions of the three gene signature-based risk scores. The distributions of risk scores and OS status were displayed in Figure 5G, and the heatmap showed the expression levels of the m6A regulators in the low- and high-risk groups. Different clinical characteristics and the expression of m6A regulators between the two groups (Figure 5H, n = 158) showed that the different group are related to the survival status but the age, race, gender, IDH, and p53 mutant. As shown in Figure 5I, the ROC curve revealed that the risk score could perfectly predict the five-year survival rates for GBM patients (AUC = 81.9%). We also calculated the OS and disease-free survival (DFS) rates of HNRNPC in the TCGA GBM dataset by GEPIA (Figures 5J,K). We observed that high expression of HNRNPC was associated with a good prognosis. Also, we identified the relationship between the clinical and the prognosis in TCGA database by univariate and multivariate Cox regression. As is shown in the Figures 6A,B, only the risk-score are associated to the prognosis. These results confirmed that the risk score originating from the selected m6A RNA methylation regulators could predict prognosis in GBM patients. And HNRNPC might be the key gene for prognostic value.

To verify whether HNRNPC played an important role in the development and progression of glioma, we firstly evaluated its mRNA expression level in glioma (n = 116) and normal brain (n = 28) tissues by RT-PCR, which indicated that HNRNPC expression was upregulated in glioma tissues compared with normal brain tissues (Figure 7A, p < 0.01). As shown is Additionally, we observed a relationship between the HNRNPC expression level and histological malignancy in tissues with different grades. As shown in Figure 7B, with the increase in HNRNPC expression, the malignancy of glioma showed an increasing tendency (p < 0.01). We examined the relationship of HNRNPC levels with overall survival (OS) rates through Kaplan–Meier analysis and log-rank test in 61 glioma cases which were collected between 2013 and 2014 with a 5-years followed up information, showing that the high expression of HNRNPC seems to be correlated with a good prognosis (Figure 7E, Hazard Ratio (HR) = 1.892, 95%Confidence Interval (CI) = 1.103–3.244, and p = 0.0205). All the results verified our prediction that HNRNPC might be a key gene in progression and development of GBM. The upregulation of HNRNPC in glioma tissues was further evaluated by Western blot (WB) and immunohistochemical (IHC) staining. For WB, we used 5 normal brain tissues, 11 low-grade (grades I-II) malignant glioma tissues, and 16 high-grade (grades III-IV) malignant glioma tissues and found that the protein level of HNRNPC was dramatically higher in high-grade and low-grade malignant glioma tissues than in normal brain tissues (Figure 7C, P < 0.05). HNRNPC expression was further examined by IHC staining in another 25 glioma tissues. We investigated the positive immunoreactivity of HNRNPC in different grades of human glioma tissues, as shown in Figure 7D, which indicated that the expression of HNRNPC was significantly related to the grades.

To determine the clinical significance of HNRNPC in glioma, we examined the correlation of HNRNPC expression with clinicopathological parameters (Table 1), observing that increased HNRNPC expression level was significantly associated with tumor grade of glioma (p = 0.0001), rather than age, gender, and tumor location. All the results implied that HNRNPC might play a vital role in glioma progression.

To further elucidate the correlation of HNRNPC expression with prognosis of glioma patients, univariate analysis demonstrated that the clinical stage (HR = 2.508, 95% CI = 1.450–4.338, and p = 0.001) and HNRNPC expression level (HR = 0.510, 95% CI = 0.282–0.923, and p = 0.026) were associated with prognosis. Multivariable Cox regression analysis confirmed that the clinical stage (HR = 2.688, 95% CI = 1.359–5.314, and p = 0.004) and the highly expressed HNRNPC (HR = 0.520, 95% CI = 0.283–0.955, and p = 0.035) were independent prognostic factors for OS of glioma patients (Table 2). These results suggested that upregulated HNRNPC seemed to be a prognostic marker for glioma patients.