We systematically identified 24 m6A regulator genes in our study from previous studies (16, 21). These m6A genes included eight writers (METTL3, METTL14, RBM15, RBM15B, WTAP, KIAA1429, CBLL1, and ZC3H13), two erasers (ALKBH5 and FTO), and 14 readers (YTHDC1, YTHDC2, YTHDF1, YTHDF2, YTHDF3, IGF2BP1, IGF2BP2, IGF2BP3, HNRNPA2B1, HNRNPC, FMR1, LRPPRC, ELAVL1, and EIF3A). Unsupervised clustering analysis was then conducted to comprehensively identify differential m6A modification patterns using the ConsensuClusterPlus package (22). Finally, the TCGA-KIRC cohort was classified into several clusters with different biological functions using a consensus clustering algorithm.

First, we downloaded 50 hallmark pathways from the MSigDB database (23). These 50 pathways systematically reflect the majority of the biological functions of humans. The GSVA algorithm was applied to calculate the enrichment scores of these pathways using the “GSVA” R package (24). Then, we analyzed difference in these pathways between different m6A clusters using the LIMMA algorithm (25). An adjusted P < 0.05 was considered statistically significant. Second, the limma R package's empirical Bayesian approach was applied to determine differentially expressed genes (DEGs) between different m6A clusters. The significance criteria for determining DEGs were set as an adjusted P < 0.05 and |logFC|>1. Finally, we performed Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses using the ClusterProfiler R package based on these DEGs.

The anticancer immune response, also called the cancer immunity cycle, is composed of seven key steps in the TME: the release and presentation of cancer cell antigens (Steps 1 and 2), the priming and activation of the immune system (Step 3), trafficking and infiltration of immune cells into tumors (Steps 4 and 5), and recognition and killing of cancer cells by T cells (Steps 6 and 7) (26). The activities of these seven steps were downloaded from http://biocc.hrbmu.edu.cn/TIP/ (27). Then, the single-sample gene-set enrichment analysis (ssGSEA) algorithm was used to quantify the relative abundance of tumor-infiltrating immune cells (TIICs) in the TME based on specific immune cell gene sets obtained from the study of Charoentong (Supplementary Table 2) (28). In addition, to avoid calculation errors caused by different algorithms and mark gene sets, we validated the infiltration level of TIICs using Cibersort-ABS, xCell and TIMER algorithm (29–31).

Mariathasan et al. revealed a set of gene signatures related to immune processes and stromal pathways, such as the CD8 T-effector signature, epithelial-mesenchymal transition (EMT) markers, and the panfibroblast TGF-b response signature (Pan-FTBRS) (20). We also collected 19 gene signatures related to the clinical response to the anti-PD-L1 agent atezolizumab (Supplementary Table 3). The ssGSEA algorithm was used to calculate the enrichment score of individuals.

The R package “WGCNA” was used to develop the gene co-expression network and to identify the m6A cluster-related module (32). First, TPM data from the TCGA-KIRC dataset were tested to determine whether they were good genes or samples. Then, the filtered genes were used to calculate the connection strength and to develop a scale-free network. The gradient method was used to test the scale independence and modules' average connectivity degree. The degree of independence was set as 0.85, and then we chose a suitable power value when the connectivity degree was relatively higher (33). Next, scale-free gene co-expression networks were generated using the selected power value. A heatmap was drawn to describe the interactions between different modules and clinical characteristics, and we chose the module that had the strongest relationship with the m6A cluster.

An m6A score was developed to quantify the m6A modification pattern in an individual patient with KIRC. First, we conducted univariate Cox analysis on genes of the module that had the strongest relationship with the m6A cluster and subsequently identified those genes with prognostic value. Similar to previous studies, we then performed principal component analysis (PCA) on these prognostic genes to calculate principal component 1, which was used for m6A score calculation (16, 34).

where i is the selected gene.

To confirm the robustness of this m6A score, we validated the prognostic value and the association between the m6A score and immunological characteristics of the TME in an independent KIRC cohort (GSE22541).

The functions significantly differed among m6A clusters. We further compared the drug sensitivities between different m6A clusters. First, we collected 184 common anticancer drugs and their target genes from the DrugBank database (www.drugbank.ca). In addition, we validated the predictive value of the m6A score for the response to ICB in three external immunotherapy cohorts.

Correlations between m6A regulators, m6A score and cancer immunity cycle and m6A score and pathways related to the ICB response were explored by Spearman coefficients and distance correlation analyses. Continuous variables fitting a normal distribution between binary groups were compared using a t-test and presented as mean ± standard deviation (SD). Otherwise, the Mann-Whitney U test was applied. Chi-square or Fisher exact tests were used to compare differences between categorical variables. The “survcutpoint” function for the maximum rank statistic was applied to determine the optimal cutoff value of the m6A score. The survival curves for prognostic analyses of categorical variables were generated using the Kaplan-Meier method, while the log-rank test was applied to estimate the statistical significance. The hazard ratio (HR) for m6A regulators was calculated using univariate Cox regression model. The independent prognostic factor of m6A score was conducted using multivariate Cox regression model and the forestplot R package was used to visualize the results. The receiver operating characteristic (ROC) curve and area under the curve (AUC) were conducted to assess the specificity and sensitivity of m6A score using time ROC R package. The mutations of m6A regulators and mutation profiles between high and low m6A score groups were visualized using maftools R package. The level of significance was set at P < 0.05, and all statistical tests were two-sided. Finally, all statistical data analyses were implemented using R software, version 3.6.3 (http://www.r-project.org).

We first analyzed the expression patterns of 24 m6A genes in KIRC and normal tissues. Interestingly, the majority of m6A writers and readers, such as METTL14, EIF3A, YTHDC1, YTHDF1, and YTHDF2, were significantly downregulated in KIRC compared to normal tissues. In contrast, expression of two m6A eraser genes (FTO and ALKBH5) was significantly higher in KIRC (Figure 2A). This expression imbalance between m6A writer and eraser genes may lead to abnormal m6A modification patterns and consequently promote the development of KIRC. Similarly, most of the m6A genes were prognostic factors. METTL14, RBM15, KIAA1429, CBLL1, YTHDC2, ZC3H13, FMR1, RBM15B, YTHDC1, FTO, LRPPRC, YTHDF2, YTHDF3, and EIF3A were favorable prognostic factors. On the other hand, METTL3, IGF2BP1, IGF2BP2, IGF2BP3, and HNRNPA2B1 were adverse prognostic factors (Figure 2B). Based on the expression of these 24 m6A genes, we could completely distinguish KIRC samples from normal samples (Figure 2C). These results suggested that m6A genes are potential diagnostic and prognostic predictors in KIRC.

Next, we assessed the CNV and mutation profiles of 24 m6A genes. Analysis of CNV data revealed prevalent CNV alterations in 24 m6A genes, and most were focused on amplification of YTHDC2, while RBM15 and RBM15B had the highest frequency of CNV deletion (Figure 2D). However, mutations of m6A genes were not frequent. Among 417 KIRC samples, only 66 (15.83%) exhibited mutations in m6A genes. ZC3H13 exhibited the highest mutation frequency at 4%, followed by YTHDC2 (2%) (Figure 2E). Finally, the close connections between the majority of m6A genes laid the foundation for the subsequent m6A clustering analysis (Figure 2F, Supplementary Table 4).

Figure 3A shows the comprehensive landscapes of 24 m6A genes concerning their prognostic value, correlations, and groups. Most of them were prognostic factors and were significantly correlated with each other, which prompted us to perform a comprehensive unsupervised clustering analysis based on these 24 m6A gene expression profiles. The results were robust when the TCGA-KIRC cohort was divided into two independent clusters. One hundred six patients were classified into m6A cluster 1, whereas the remaining 423 patients were classified into m6A cluster 2. m6A cluster 1 exhibited a significantly poorer prognosis (P = 0.00057) (Figure 3B). The DEGs between m6A clusters are displayed in a heatmap and volcano plot (Figures 3C,D, Supplementary Table 5). The results of GO analysis suggested that these DEGs were enriched in several biological processes, including organic anion transport, metal ion transmembrane transporter activity, collagen-containing extracellular matrix, and cellular divalent inorganic cation homeostasis (Supplementary Figures 2A–C, Supplementary Table 6). The results of KEGG analysis indicated that these DEGs were enriched in pathways such as neuroactive ligand-receptor interaction, bile secretion, vascular smooth muscle contraction, mineral absorption, complement and coagulation cascades, serotonergic synapse, protein digestion and absorption, and leukocyte transendothelial migration (Supplementary Figure 2D, Supplementary Table 7). Finally, the enrichment scores of many hallmark signatures significantly differed between the two clusters. As shown in Figure 3E, TGF-beta signaling, Wnt-beta catenin signaling, protein secretion, PI3K-Akt-Mtor signaling, androgen response, heme metabolism, mitotic spindle, and Notch signaling were enriched in m6A cluster 2. In contrast, spermatogenesis, estrogen response late, and KRAS signaling DN were enriched in m6A cluster 1 (Figure 3E, Supplementary Table 8).

We next comprehensively correlated the m6A clusters with immune phenotypes. First, we focused on the activities of anticancer immunity cycles. The activity of priming and activation of the immune system of m6A cluster 1 was significantly higher than that of m6A cluster 2, while the activities of releasing and presenting cancer cell antigens were lower (Figure 4A). In addition, the activities of T cell recruiting, B cell recruiting, and dendritic cell recruiting were consistently higher in m6A cluster 1 (Figure 4A). Finally, activities of recognition of cancer cells by T cells were higher in m6A cluster 1. To confirm these findings, we directly compared the infiltration level of tumor-infiltrating immune cells between m6A clusters. As expected, the abundance of several antitumor immune cells, such as activated CD8 T cells, activated CD4 T cells, CD56bright natural killer cells and type 17 T helper cell, was significantly higher in m6A cluster 1 than in m6A cluster 2 (Figure 4B). However, the abundance of the most recognized protumor immune cells, including regulatory T cells, immature dendritic cells, and plasmacytoid dendritic cells, was significantly downregulated in m6A cluster 1 (Figure 4B). Based on these results, we proposed that m6A cluster 1 may be an inflammatory immune phenotype, while m6A cluster 2 may be a non-inflammatory phenotype. Previous research demonstrated that stroma-associated pathways, such as EMT and Pan-FTBRS signatures, inhibited the anticancer immunity in TME (20). Here, EMT1, EMT3, and Pan-F-TBRS enrichment scores were significantly downregulated in m6A cluster 1 (Figure 4C).

Inflammatory tumor phenotypes are more sensitive to ICB (35, 36). Consistently, pathways that were positively related to the ICB response, such as RNA degradation, the cell cycle, and DNA replication, were enriched in m6A cluster 1 (inflammatory phenotype). In contrast, the pathway cytokine-cytokine receptor interaction negatively related to the ICB response was enriched in m6A cluster 2 (non-inflammatory phenotype) (Figure 4D). Therefore, we confirmed that m6A cluster 1 might represent an inflamed phenotype from the aspect of immunotherapy response.

All tumor data from the TCGA-KIRC dataset were used to develop the gene co-expression network and to identify m6A cluster-related modules. All KIRC samples with full clinical characteristics were included in the co-expression analysis (Figure 5A). The “WGCNA” package was used to allocate genes with similar expression patterns into different modules. In this study, we chose the soft threshold as 7 (scale-free R² = 0.85) to develop a scale-free network. As shown in Figure 5B, a total of 29 modules were recognized. The modules with the most significant association with clinical characteristics had the greatest biological meanings. The turquoise module was found to have the highest association with the m6A cluster (r = 0.64, p = 4e-64; Figure 5C). We chose the turquoise module to be analyzed in the subsequent steps, and the turquoise module was also related to tumor grade and stage. The genes in the turquoise modules were significantly co-expressed (cor = 0.81, P < 1e-200; Figure 5D). Among these genes, 2,214 were significantly related to prognosis (Supplementary Table 9). Then, the m6A score was calculated for individuals using the PCA algorithm.

m6A score was lower in m6A cluster 1 (Figure 6A). Similar to the performance of m6A cluster 1, patients in the low m6A score group exhibited poorer prognosis than patients in the high m6A score group (Figure 6B). Also, m6A score still remained an independent prognosis factor in multivariate Cox regression analysis (p = 0.01, Supplementary Figure 3A). The Q-Q plot of the model showed that the residuals are approximately normally distributed (Supplementary Figure 3B) and the AUC at 5 years showed that the predictive accuracy of m6a score was comparative to tumor stage (Supplementary Figure 3C). There were consistent correlations between the m6A score and the immune phenotype. The CD8 T effector signatures were enriched in the low m6A score group (Figure 6C). The abundance of antitumor immune cells, including activated CD8 T cells, activated CD4 T cells, activated dendritic cells, CD56bright natural killer cells, central memory CD4 T cells, natural killer T cells, type 1 T helper cells, and type 17 T helper cells was significantly upregulated in the low m6A score group (Figure 6D). However, the abundance of protumor immune cells, including immature dendritic cells and plasmacytoid dendritic cells, was downregulated in the low m6A score group (Figure 6D). We validated the infiltration level of TIICs using Cibersort-ABS, xCell, and TIMER algorithm (Supplementary Figures 4–6). Generally, most of the algorithms showed that m6A score was negatively correlated with anti-tumor immune cells, including CD8 T cells, CD4 T cells, and natural killer T cell. Except TIMER algorithm showed that CD8 T cells was positively correlated with m6A score. This could be the calculation errors caused by different algorithms and mark gene sets. In addition, the EMT1 and EMT3 pathways were enriched in the high m6A score group (Figure 6E). Meanwhile, the m6A score was negatively related to the activities of several critical anticancer immunity cycles, such as priming and activation, T cell recruiting, CD8 T cell recruiting, CD4 T cell recruiting, dendritic cell recruiting, Th17 cell recruiting, and infiltration of immune cells into tumors (Figure 6F, Supplementary Table 10). These findings suggested that the low m6A score group may have an inflammatory phenotype.

As expected, m6A scores were negatively correlated with pathways that were positively related to the ICB response, such as RNA degradation, cell cycle, and DNA replication. In contrast, the m6A score was positively related to the cytokine-cytokine receptor interaction pathway, which was negatively related to the ICB response (Figure 6F, Supplementary Table 11). Finally, several common immune checkpoints, such as CTLA-4, PD-1, LAG-3, LAALS3, and TIGIT, were highly expressed in the low m6A score group (Figure 6G).

In summary, the m6A score predicts the immune phenotype and clinical response to ICB.

Genomic mutations are a prominent factor in initiating malignancy. Here, we analyzed distribution differences in the top 20 somatic mutations between m6A score groups using the maftools R package. The most common mutations in KIRC were VHL and PBRM1. There was no difference in the VHL mutation between the m6A score groups (Figure 7A). The mutation frequencies of TTN (32 vs. 23%), SETD2 (19 vs. 9%), BAP1 (16 vs. 7%), and MUC16 (15 vs. 7%) were markedly higher in the low m6A score group suggesting that these mutations may be m6A score-specific mutations in KIRC. In general, a more extensive tumor mutation burden was presented in the low m6A score group than in the high m6A score group (97.4 vs. 90.67%) (Figure 7A). Consequently, the TMB quantification analysis revealed that the low m6A score group was markedly correlated with a higher TMB (Figure 7B). However, there was no difference in MSI status between the two m6A score groups (Figure 7C).

Similar to the performance of the m6A score in the TCGA-KIRC cohort, we found that the low m6A score group had a poorer prognosis in the GSE22541 cohort as well (Figure 8A). Meanwhile, the m6A score was negatively correlated with the activities of many anticancer immunity cycles, such as the recognition of cancer cells by T cells (Figure 8B, Supplementary Table 12). Furthermore, the infiltration levels of activated CD8 T cells, activated CD4 T cells, activated dendritic cells, central memory CD8 T cells, natural killer T cells, type 1 T helper cells, and type 17 T helper cells were significantly higher in the low m6A score group (Figure 8C). Finally, the m6A score was negatively related to most pathways that predicted higher ICB response rates (Figure 8B, Supplementary Table 13). These results confirmed that the m6A score might be a robust predictor of immune phenotype, prognosis, and ICB response.

We further explored the role of the m6A score in guiding clinical decision making in KIRC. First, we found that the sensitivities of many anticancer drugs were significantly different between m6A score groups (Supplementary Table 14). Targeted therapy was the first-line treatment option for advanced KIRC. Here, we collected the targeted therapy drugs used in KIRC and their targeted genes from the DrugBank database: sorafenib with its targeted genes including BRAF, FLT1, FLT3, FLT4, KDR, KIT, and RAF1; sunitinib with its targeted genes including CSF1R, FLT1, FLT3, FLT4, KDR, and RET; pazopanib with its targeted gene SH2B3; and bevacizumab with its targeted gene VEGFA. Interestingly, all targeted therapy drug sensitivities were significantly lower in the low m6A score group (Figure 9A). These results indicate that the m6A score may identify suitable candidates to receive targeted therapy.

Although findings from TCGA-KIRC and GSE22541 cohorts suggested that the m6A score predicts ICB response, it would be more convincing to validate these results in cohorts that received ICB. First, in a KIRC cohort that received anti-PD-1 therapy (nivolumab), we demonstrated that the clinical benefit rate was higher in the low m6A score group than in the high m6A score group (p = 0.26; Figure 9B). Regrettably, because of the small sample size, we didn't find significantly differences. The prognosis of the low m6A score group was better than in the high m6A score group (p = 0.039; Figure 9C). It is worth noting that this survival outcome was contrary to the results showing that the prognosis of the low m6A score group was worse in the TCGA-KIRC and GSE22541 cohorts. These differences in outcome were due to the response rate of immunotherapy being more likely to determine the prognosis of an immunotherapy cohort when compared to other prognostic risk factors, such as the m6A score. Additionally, we successfully validated the role of the m6A score in predicting the response to ICB in two other cancer cohorts, including the IMvigor210 cohort (bladder cancer) and GSE78220 cohort (melanoma) (Figures 9D–G). These findings revealed that this m6A score may represent a generalized predictor for response to ICB in other cancer types as well.