m⁶A RNA methylation regulators contribute to malignant progression and have clinical prognostic impact in gliomas

Expression of m⁶A RNA methylation regulators is correlated with clinicopathological features in gliomas

Considering the important biological functions of each m⁶A RNA methylation regulator in tumorigenesis and development, we systematically investigated the relationships between each individual m⁶A RNA methylation regulator and the pathological features of gliomas, including WHO grade, IDH status, and 1p/19q codeletion status. The expression level of each m⁶A RNA methylation regulator and WHO grades are presented as heatmaps (Fig. 1A and B), showing that the expressions of most m⁶A RNA methylation regulators are significantly associated with WHO grades. The significant correlations between WHO grades and expression levels of WTAP, RBM15, YTHDF, ALBKH5, and FTO were also confirmed by quantitative analyses in both the CGGA (Fig. 1C) and TCGA datasets (Fig. 1D). As the WHO grade increased, the expression of WTAP, RBM15, YTHDF and ALBKH5 increased, while the expression of FTO decreased.

We then studied the relationship between IDH status and expression levels of each m⁶A RNA methylation regulator in lower-grade gliomas (LGG, Fig. 1E and 1F) and glioblastomas (GBM, Fig. 1G), respectively. The results showed that expression levels of ALKBH5, YTHDF2, RBM15, METTL3, METTL14, FTO and YTHDC1 are significantly different between LGG with mutant-IDH and LGG with wildtype-IDH in both CGGA (Fig. 1E) and TCGA (Fig. 1F) datasets. FTO, YTHDC1 and METTL3 were also differentially expressed between GBM with and without IDH mutation in the CGGA dataset (Fig. 1G). However, we were unable to compare the expression of m⁶A RNA methylation regulators between GBM with and without IDH mutation in the TCGA dataset because of the dramatically imbalanced numbers of IDH-mutant GBM (n=10) and IDH-wildtype GBM (139). In the LGG with mutant-IDH, we observed that the expressions of YTHDF2 and WTAP were highly correlated with the status of 1p/19q codeletion (1p/19q codel).

We noticed that the expression of both METTL3 and METTL14 were not correlated with the WHO grade of gliomas. Considering METTL16 is also a putative m⁶A methyltransferase and knockdown of METTL16 could result in approximately 20% decrease of m⁶A [29], we also investigated the expression of METTL16 in gliomas with different malignancy status in the CGGA dataset (Figure S1). METTL16 expression was also not correlated with the WHO grade of gliomas. Though METTL16 expression correlated to the IDH status in LGG (P<0.001) and GBM (P<0.05), it did not have prognostic value in any pathological subgroup of gliomas (data not shown). In addition, the expression information was also not included in the TCGA dataset, so METTL16 was not included in the following section of this study.

Moreover, we also observed that the genetic change (mutation or coy number change) frequencies of the thirteen m⁶A RNA methylation regulators were very low (all ≤ 1.1%) in gliomas (Figure S2), indicating that the expression changes of these m⁶A RNA methylation regulators were not caused by the genetic changes of the corresponding genes.

Consensus clustering of m⁶A RNA methylation regulators identified two clusters of gliomas with distinct clinical outcomes and clinicopathological features

Based on the expression similarity of m⁶A RNA methylation regulators, k = 3 seemed to be an adequate selection with clustering stability increasing from k = 2 to 10 in the CGGA datasets (Fig. 2A and B). However, we noticed that 293 out of 309 gliomas clustered into one of the two subgroups in the CGGA dataset (Fig. 2C and Figure S3). Thus, we compared the clinicopathological features of these two subgroups clustered by k=2, namely, RM1 and RM2 (Fig. 2C) datasets. The RM1 subgroup is significantly correlated with younger age at diagnosis (P < 0.0001), lower grade (P < 0.0001), proneural or neural subtypes (P < 0.0001), IDH-mutational status (P < 0.0001) and 1p/19q codel status (P < 0.0001). The RM2 subgroup mainly contains gliomas with an older age at diagnosis, glioblastoma (GBM) phenotype, classic or mesenchymal subtypes, IDH-wildtype status and 1p/19q noncodel status (Table S2). Moreover, we observed a significantly shorter overall survival (OS) in the RM2 subgroup than the RM1 subgroup (Fig. 2D).

Categories identified by consensus clustering are closely correlated to the malignancy of gliomas

The above findings suggested that the clustering result was closely correlated to the malignancy of the gliomas. To better understand the interactions among the thirteen m⁶A RNA methylation regulators, we also analyzed the interaction (Fig. 3A) and correlation (Fig. 3B) among these regulators. WTAP seems to be the hub gene of the ‘writers’, and its interactions or coexpressions with METTL3, METTL14, KIAA1429 and ZC3H13 are supported both by experimental data and by text mining in the String database (Fig. 3A). The expression of WTAP was also significantly correlated with the ‘writers’ of METTL14, KIAA1429 and RBM15 in gliomas (Fig. 3B). There are several independent interaction groups in ‘readers’, suggesting the diverse functions of different ‘readers’ (Fig. 3A), but the expressions of YTHDF2, YTHDC2 and YTHDF1 were significantly correlated with each other in gliomas (Fig. 3B). Only text mining evidence supported the interaction of ALKBH5 and FTO in the String database (Fig. 3A), and the expressions of ALKBH5 and FTO were negatively correlated with each other in gliomas (Fig. 3B). Moreover, the expressions of WTAP, RBM15, YTHDF2, YTHDF1 and ALKBH5 were highly correlated with each other, and all of their expressions were negatively correlated with FTO in gliomas (Fig. 3B). These findings were consistent with that the expression levels of WTAP, RBM15, YTHDF2, YTHDF1 and ALKBH5 being positively correlated with the increasing malignancy of gliomas, while the expression levels of FTO are negatively correlated with the increasing malignancy of gliomas (Figs. 1 and 2).

We further used principal component analysis (PCA) to compare the transcriptional profile between RM1 and RM2 subgroups. The results showed a clear distinction between them (Fig. 3C). We identified genes that were significantly upregulated (Score(d) for SAM > 8, fold change >2, and normalized P < 0.01) or downregulated (Score(d) for SAM < -8, fold change < 0.5, and normalized P < 0.01) in the RM2 subgroup, and then annotated their function using gene ontology (GO) pathway analysis for biological processes (Fig. 3D). The results indicated that upregulated genes are enriched in malignancy-related biological processes, including cell proliferation, extracellular matrix organization, angiogenesis, migration, and immune response, among others. Similar changes in corresponding signaling pathways were also observed in the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis (Fig. 3E). Furthermore, gene set enrichment analysis (GSEA) revealed that the malignant hallmarks of tumors, including epithelial-mesenchymal transition (NES=1.88, normalized P < 0.001), TNFα signaling via NF-κB (NES=1.79, normalized P < 0.001), inflammatory response (NES=1.74, normalized P = 0.041) and IL6/JAK/ STAT3 signaling (NES=1.79, normalized P < 0.001), were significantly associated with the RM2 subgroup (Fig. 3F). All of these findings indicated that the two categories identified by consensus clustering are closely correlated with the malignancy of gliomas.

Prognostic value of m⁶A RNA methylation regulators, and a risk signature built using seven selected m⁶A RNA methylation regulators

We next sought to investigate the prognostic role of m⁶A RNA methylation regulators in gliomas. We performed a univariate Cox regression analysis on the expression levels in the CGGA dataset (Fig. 4A). The results indicated that eleven out of thirteen tested genes are significantly correlated with OS (P < 0.05). Among these eleven genes, ALKBH5, YTHDF1, YTHDF2, HNRNPC, RBM15, KIAA1429, and WTAP are risky genes with HR > 1, while FTO, YTHDC1, ZC3H13, and METTL3 are protective genes with HR < 1. We also investigated the prognostic value for each m⁶A RNA methylation regulator in stratified LGG (Table S3) and GBM (Table S4). For the genes with prognostic value in LGG, the expression levels of YTHDF2, WTAP, ALKBH5, RBM15, KIAA1429, HNRNPC, YTHDF1, and FTO are significantly correlated with the overall survival (OS) of patients with IDH-mutant and 1p/19q noncodel LGG. Among these genes, YTHDF2, KIAA1429, HNRNPC, and YTHDF1 also have prognostic value in IDH-wildtype LGG. In the GBM, the expressions of FTO, YTHDF2, and RBM15 have prognostic value in IDH-wildtype GBM; and FTO also has prognostic value in IDH-mutant GBM.

To better predict the clinical outcomes of gliomas with m6A RNA methylation regulators, we applied the least absolute shrinkage and selection operator (LASSO) Cox regression algorithm to the eleven prognosis-associated genes in the CGGA dataset, which we used as a training set (Fig. 4A). Seven genes were selected to build the risk signature based on the minimum criteria, and the coefficients obtained from the LASSO algorithm were used to calculate the risk score for both the training dataset (CGGA) and the validation dataset (TCGA). To investigate the prognostic role of the seven-gene risk signature, we separated the glioma patients in the CGGA (n = 309) and TCGA (n = 595) datasets into low- and high-risk groups based on the median risk score and observed significant differences in OS between the two categories (both P < 0.0001; Fig. 4B, C).

Prognostic risk scores showed strong associations with clinicopathological features in gliomas

The heatmap shows the expression of the seven selected m⁶A RNA methylation regulators in high- and low-risk patients in the CGGA dataset (Fig. 5A). We observed significant differences between the high- and low-risk groups with respect to WHO grade (P < 0.001), age (P < 0.001), IDH status (P < 0.001), 1p/19q codel status (P < 0.001), TCGA subtypes (P < 0.001), and RM1/2 subgroups (P < 0.001). We also examined the association between the risk scores and each clinicopathological feature. We observed that risk scores are significantly different between patients stratified by WHO grade, TCGA subtype, age, IDH status, 1p/19q codel status, and RM1/2 subgroups, but not by gender, in both the CGGA (Fig. 5B–H) and TCGA datasets (Figure S4). The ROC curve shows that the risk score can perfectly predict three-year survival rates for glioma patients (AUC = 90.3%), RM1/2 subgroups (AUC = 98.6%), glioma IDH-mutant status (AUC = 86.3%) and glioma 1p/19q codel status (AUC = 82.3%). Furthermore, the predictive efficiency is better than that of the WHO grade and age (Fig. 5I-L). These results indicate that the risk scores calculated with the signature can accurately predict glioma patient outcomes and clinicopathological features, especially for the RM1/2 subgroups.

We then performed univariate and multivariate Cox regression analyses for the CGGA dataset to determine whether the risk signature is an independent prognostic indicator. By univariate analysis, the risk score, 1p/19q codel status, IDH status, age and WHO grade were all correlated with the OS. When including these factors into the multivariate Cox regression, the risk score and WHO grade remained significantly associated with the OS (both P < 0.001, Fig. 5M). Similar results were found in the validation of the TCGA dataset; the risk score (P = 0.027), IDH status (P < 0.001), age (P < 0.001) and WHO grade (P = 0.005) remained significantly associated with the OS in the multivariate regression (Fig. 5N). These results confirmed that the risk score derived from m⁶A RNA methylation regulators can independently predict prognosis in glioma patients.

We also determined the prognostic value of the risk signature for different WHO grades. We found that patients with high risk scores had significantly shorter OS than those with low scores in WHO grade II and III gliomas (Fig. 6A, B) and GBM (Fig. 6C). Meanwhile, we also observed that the risk score had prognostic value in gliomas stratified by the integrated diagnosis of WHO grade 2016 (Figure S5). Moreover, the risk score could also predict the mesenchymal subtype in GBM for both CGGA (Fig. 6D) and TCGA datasets (Fig. 6E). Consistent with this finding, the patients with a high-risk score were also more sensitive to TMZ treatment in GBM for both CGGA and TCGA datasets (Fig. 6F-I).

Datasets

The RNA-seq transcriptome data and corresponding clinicopathological information were obtained for 309 glioma patients from CGGA (www.cgga.org.cn) and 595 glioma patients from TCGA (http://cancergenome.nih.gov/). For the RNA-seq data, reads per kilobase of transcript per million mapped reads values were obtained for the CGGA samples [27], and values normalized by expectation-maximization were obtained for the TCGA samples. The study protocol was approved by the ethics committee of Beijing Tiantan Hospital. Clinicopathological information for the CGGA and TCGA datasets is summarized in Supplementary Table S1.

Selection of m⁶A RNA methylation regulators

We first collated a list of sixteen m⁶A RNA methylation regulators from published literature [6,15,21], and then we restricted the list to genes with available RNA expression data in the CGGA and TCGA datasets. This yielded a total of thirteen m⁶A RNA methylation regulators. Then, we systematically compared the expression of these m⁶A RNA methylation regulators in gliomas with different clinicopathological features.

Bioinformatic analysis

To investigate the function of m⁶A RNA methylation regulators in gliomas, we clustered the gliomas into different groups with “ConsensusClusterPlus” (50 iterations, resample rate of 80%, and Pearson correlation, http://www.bioconductor.org/). We used PCA with the R package for R v3.4.1 to study the gene expression patterns in different glioma groups. GO and KEGG pathway enrichment analyses were performed with the Database for Annotation, Visualization, and Integrated Discovery (http://david.abcc.ncifcrf.gov/home.jsp) to functionally annotate genes that are differentially expressed in different groups. Interactions among m⁶A RNA methylation regulators were analyzed using the STRING database (http://www.string-db.org/). GSEA was performed to investigate the functions correlated with different subgroups of gliomas.

To determine the prognostic value of m⁶A RNA methylation regulators, we performed univariate Cox regression analyses of their expression in the CGGA dataset. From this, we identified eleven genes significantly associated with survival (P < 0.05), which we selected for further functional analysis and development of a potential risk signature with the LASSO Cox regression algorithm [43–46]. Finally, seven genes and their coefficients were determined by the minimum criteria, selecting the best penalty parameter λ associated with the smallest 10-fold cross validation within the training set. The risk score for the signature was calculated using the formula:

$\text{Risk score = }{\displaystyle {\sum }_{i=1}^n{{\displaystyle \text{Coef}}}_{\text{i}}}\times {{\displaystyle x}}_{\text{i}}$,

where Coef_i is the coefficient, and ${\mathrm{x}}_{\mathrm{i}}$ is the z-score-transformed relative expression value of each selected gene. This formula was used to calculate a risk score for each patient in both the training (CGGA) and validation (TCGA) datasets.

Statistical analysis

One-way ANOVA was used to compare the expression level of m⁶A RNA methylation regulators in gliomas with different WHO grades and TCGA subtypes, and t-tests were used to compare the expression levels in gliomas for age, gender, IDH status, and 1p/19q codel status.

Patients were clustered into two groups by consensus expression of m⁶A RNA methylation regulators or were divided into high- and low-risk groups using the median risk score (derived from the risk signature) as the cutoff value. Chi-square tests were used to compare the distribution of gender, WHO grade, TCGA subtype, IDH status, and 1p/19q codeletion status between the two risk groups.

To compare the risk scores of the signature for gliomas with different clinicopathologies, a one-way ANOVA or t-test was performed to compare the risk scores in patients grouped by clinical or molecular-pathological characteristics. Univariate and multivariate Cox regression analyses were performed to determine the prognostic value of the risk score and various clinical and molecular-pathological characteristics.

The prediction efficiency of the risk signature, WHO grade and age for 3-year survival, RM1/2 groups, IDH-mutant status and 1p/19q codeletion status were tested with operating characteristic (ROC) curves.

The Kaplan–Meier method with a two-sided log-rank test was used to compare the OS of the patients in the RM1/2 groups or in the high- and low-risk groups. All statistical analyses were conducted using R v3.4.1 (https://www.r-project.org/), SPSS 16.0 (SPSS Inc., Chicago, IL) and Prism 7 (GraphPad Software Inc., La Jolla, CA).