N 6 -Methyladenosine Regulator-Mediated Methylation Modification Patterns with Distinct Prognosis, Oxidative Stress, and Tumor Microenvironment in Renal Cell Carcinoma

Objective : Emerging evidence suggests the biological implications of N 6 -methyladenosine (m 6 A) in carcinogenesis. Herein, we systematically analyzed the role of m 6 A modification in renal cell carcinoma (RCC) progression. Methods : Based on 23 m 6 A regulators, unsupervised clustering analyses were conducted to determine m 6 A modification subtypes across 893 RCC specimens in the Cancer Genome Atlas (TCGA) cohort. By performing principal component analysis (PCA) analysis, m 6 A scoring system was developed for evaluating m 6 A modification patterns of individual RCC patients. The activity of signaling pathways was assessed by gene-set variation analysis (GSVA) algorithm. The single-sample gene set enrichment analysis (ssGSEA) algorithm was applied for quantifying the infiltration levels of immune cells and the activity of cancer immunity cycle. Drug responses were estimated by genomics of drug sensitivity in cancer (GDSC), the Cancer Therapeutics Response Portal (CTRP) and Preservice Research Institute for Science and Mathematics (PRISM) database. databases. Results : Five m 6 A modification subtypes were characterized by different survival outcomes, oxidative stress, cancer stemness, infiltrations of immune cells, activity of cancer immunity cycle, programmed cell death 1 (PD-1)/programmed cell death ligand 1 (PD-L1) expression and microsatellite instability (MSI) levels. According to m 6 A score, RCC patients were categorized into high and low m 6 A score groups. Patients with high m 6 A score displayed a prominent survival advantage, and the prognostic value of m 6 A score was confirmed in two anti-PD-1/PD-L1 immunotherapy cohorts. m 6 A score was significantly linked to oxidative stress-related genes, and high m 6 A score indicated the higher sensitivity to axitinib, pazopanib and sorafenib and the lower sensitivity to sunitinib. Conclusion : This study analyzed the extensive regulatory mechanisms of m 6 A modification on oxidative stress, the tumor microenvironment, and immunity. Quantifying m 6 A scores may enhance immunotherapeutic effects and assist in developing more effective agents.

Among the renal cancer-related deaths, over 90% are in relation to RCC metastasis [6].Within the tumor microenvironment (TME), there are complex interactions between tumor cells, immune cells, and stromal cells [7,8].Immune checkpoint inhibitors targeting programmed cell death 1 (PD-1)/programmed cell death ligand 1 (PD-L1) combined with a tyrosine kinase inhibitor have remarkably altered the clinical management of metastatic RCC, and have been recommended as the first-line treatment option [9].Nevertheless, long-term response is low due to high risk of resistance.Thus, novel treatment options are urgently required to enhance the therapeutic effects as well as to determine biomarkers for predicting the responses and better stratifying RCC patients, thereby reducing costs and prolonging the survival duration.
N 6 -methyladenosine (m 6 A) represents the most abundant RNA modification pattern in eukaryotic cells, which typically occurs at the nitrogen-6 position of adenosine [10].M 6 A modification, a dynamically reversible process, is mediated by methyltransferases ("writers"), demethylases ("erasers"), and binding proteins ("readers") [11].This methylation modification in mRNAs exerts key functions in diverse aspects of RNA metabolism (mRNA splicing, translation efficiency, etc.), thereby participating in critical biological processes such as tumorigenesis, immunomodulatory process, and cancer metastasis [12,13].Thus, it is of importance to extensively characterize "writers", "erasers" and "readers" that alter the modification levels and recognize the chemical marks [14].Recently, fat mass and obesity associated (FTO) m 6 A demethylase inhibits ccRCC via FTO-peroxisome proliferator-activated receptor-gamma co-activator-1α (PGC-1α) signaling axis [15].M 6 A demethylase AlkB homolog 5 (ALKBH5) accelerates RCC progression through modulating aurora kinase B (AURKB) expression with an m 6 A-dependent manner [16].Nevertheless, the expression of m 6 A regulators in RCC patients with various clinicopathological features, their roles in the TME, and their prognostic implications are mostly unclear.Oxidative stress is an imbalance between oxidants and antioxidants.Overactivation of oncogenes results in the increased generation of reactive oxygen species (ROS) in tumor cells, accompanied by an increased antioxidant ability to maintain redox homeostasis at increased levels in tumor cells.Accumulated evidence suggests that m 6 A modification modulates cellular ROS levels via distinct mechanisms [17].The effects of m 6 A modification in oxidative stress remain unclear in RCC.Here, this study constructed distinct m 6 A modification subtypes characterized by distinct biological functions, oxidative stress, TME features, tumor immunity and survival outcomes in RCC.Moreover, this study developed a m 6 A scoring system for quantifying the m 6 A modification patterns for individual patients, which could be applied for elucidating immune phenotypes and predicting the prognosis and immunotherapy responses in clinical practice.

Interactions between M 6 A Regulators
The m 6 A regulators were inputted into the STRING online database (http://string-db.org/)[21].The interactions were visualized into a protein-protein interaction (PPI) network via Cytoscape software (version 3.8.2,https://cytoscape.org/)[22].Pearson correlation analysis was performed to evaluate the correlation of the mRNA expression of the m 6 A regulators across RCC samples.

Immunotherapy Response
Through Tumor Immune Dysfunction and Exclusion (TIDE) approach, the response to immunotherapy was predicted following the tumor immune evasion mechanisms [21].The expression similarity between subtypes and the patients who might respond to anti-PD-1 and anti-CTLA4 therapy was detected utilizing Subclass Mapping (SubMap) algorithm [23].

Unsupervised Clustering Analyses of 23 M 6 A Regulators
Unsupervised clustering analyses were applied to estimate the number of unsupervised classes across 893 RCC samples based on the expression profile of 23 m 6 A regulators by ConsensuClusterPlus package (version 1.60.0,https://bioconductor.org/packages/release/bioc/htm l/ConsensusClusterPlus.html) [24].The classification accuracy was assessed via t-distributed stochastic neighbor embedding (t-SNE) method.Survival differences among clusters were compared by Kaplan-Meier curves and logrank tests.

Evaluation of Stemness Features
Two stemness indicators: gene expression-based stemness index (mRNAsi) and DNA methylation-based stemness index (mDNAsi) were calculated in RCC samples by building the stemness index models with one-class logistic regression (OCLR) machine-learning method [30].The stemness value was ranging from 0 (no gene expression) to 1 (complete gene expression).

Assessment of Immunological Characteristics
Cancer immunity cycle contains release of cancer cell antigens (step 1), cancer antigen presentation (step 2), priming and activation (step 3), trafficking of immune cells to tumors (step 4), infiltration of immune cells into tumors (step 5), recognition of cancer cells by T cells (step 6), and killing of cancer cells (step 7) [31].The activation of each step was quantified in RCC samples by ssGSEA function.

Tumor Immune Landscape
Immune score and stromal score were obtained using estimation of stromal and immune cells in malignant tumours using expression data (ESTIMATE) algorithm based on gene expression data, which represented the fractions of immune and stromal cells in RCC samples [32].The enrichment levels of immune cells were estimated in RCC samples through the ssGSEA algorithm according to the expression matrix of gene symbols of tumor-infiltrating immune cells (TIICs) [33].The mRNA expression of immune checkpoints was calculated in each RCC specimen.

Identification of M 6 A-Related DEGs
Differentially expressed genes (DEGs) were screened between m 6 A clusters rough limma package [34].Genes with adjusted p < 0.05 were considered as DEGs.DEGs shared by any two clusters were identified as m 6 A-related DEGs.

Generation of m 6 A Scoring System
RCC patients were clustered into distinct gene clusters through applying unsupervised clustering analyses based on the extracted m 6 A-related DEGs.By using ConsensuClus-terPlus package, the number of gene clusters and their stability were determined across RCC samples.The accuracy of the gene clusters was validated with t-SNE method.Survival differences among gene clusters were observed with Kaplan-Meier curves and log-rank tests.Univariate-cox regression analyses were carried out for extracting prognostic m 6 A-related DEGs with p < 0.05 across RCC patients.These extracted DEGs were used for feature selection through recursive feature elimination (RFE) with random forest, followed by 10-fold cross-verification utilizing caret package.The m 6 A scoring system was quantified in individual tumors based on the curated expression profiles of the finally identified DEGs by conducting a principal component analysis (PCA) utilizing the Boruta algorithm, named the m 6 A score.Both principal component 1 and 2 were chosen as signature scores.The m 6 A score was quantified according to the following formula: m 6 A score = ∑ (PC1 i + PC2 i ), where i indicated the expression of m 6 Arelated DEGs.

Collection of Genomic and Clinical Information of Immunotherapy Cohorts
The mRNA expression profiles and follow-up data of patients with anti-PD-1 therapy were downloaded from the GSE78220 dataset (N = 27) in the Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/)repository [36] and the literature reported by Liu et al. [37] (N = 121).The prognostic value of m 6 A score was externally validated in the two immunotherapy cohorts.

Gene Set Enrichment Analysis (GSEA)
Pathways were probed by GSEA package [28].The "c5.bp.v6.2.symbols.gm"gene set was set as the reference set, which was obtained from the molecular signatures database.Pathways with adjusted p < 0.05 were significantly associated with m 6 A score.

Assessment of Sensitivity of Chemotherapeutic Agents
By employing pRRophetic algorithm [38], the half-maximal inhibitory concentration (IC50) values of chemotherapy agents were estimated for RCCs on the basis of the genomics of drug sensitivity in cancer (GDSC; www.cancerrxgene.org/)cell line expression spectrum [39] and mRNA expression profiling.

Prediction of Drug Response
Drug sensitivity data of human cancer cell lines (CCLs) were retrieved from the the Cancer Therapeutics Response Portal (CTRP) (https://portals.broadinstitute.org/ctrp) and Preservice Research Institute for Science and Mathematics (PRISM) database (https://depmap.org/portal/prism/) datasets.Both two datasets offer the area under the dose-response curve as an evaluation indicator for the responses to hundreds of compounds.Reduced area under the curve (AUC) indicates higher sensitivity to a specific drug.By K-nearest neighbor method, the missing AUCs were imputed.Since the CCLs in both datasets were retrieved from the cancer cell line encyclopedia (CCLE) project (https: //portals.broadinstitute.org/ccle/)[40], expression profiling data in CCLE project were employed for CTRP and PRISM analyses.

Statistical Analyses
Statistical analyses were conducted with R packages (version 4.0.2).Kaplan-Meier curves of overall survival (OS), disease free survival (DFS), disease-specific survival (DSS), progression free survival (PFS) and progression free interval (PFI) were depicted for RCC patients and survival differences between groups were compared with log-rank tests.Comparison between two groups was presented via student's t test or Wilcoxon test.One-way ANOVA or Kruskal-Wallis test was used for conducting comparison between three or more groups.Two-sided p < 0.05 indicated statistical significance.

Systematic Analyses of Genetic Variation, Expression, Prognostic Implication and Interactions of m 6 A Regulators in RCC
This study collected three types of RCC ccRCC, pRCC and chRCC from TCGA datasets and removed batch effects after integration (Fig. 1A,B).Finally, 893 RCC samples (chRCC (N = 65), ccRCC (N = 539) and pRCC (N = 289) and 128 normal samples were included for further analyses.In total, 23 m 6 A regulators were collected and the prevalence of CNAs among these m 6 A regulators was analyzed in RCC samples.In Fig. 1C, YTHDC2 had the prevalent gain and IGFBP2, YTHDF2 and RBM15B had the widespread loss across RCC tissues.The mRNA expression of these m 6 A regulators was compared in RCC and normal specimens.We found that most m 6 A regulators were abnormally expressed in RCC compared to normal tissues (Fig. 1D), including writers (CBLL1, METTL14, ZC3H13), eraser (FTO) and readers (FMR1, HNRNPA2B1, HNRNPC, IGF2BP1/2/3, LRP-PRC, YTHDC1/2, YTHDF2/3).We also compared the expression of m 6 A regulators in stage I&II and III&IV RCC tissues.As depicted in Fig. 1E, writers (CBLL1, METTL14, METTL3, RBM15B, ZC3H13), and readers (FMR1, IGF2BP1, IGF2BP3, LRPPRC, YTHDC1, YTHDF1) displayed remarkable differences between stage I&II and III&IV.Additionally, we noted that there were significant differences in the expression of writers (RBM15, WTAP), erasers (ALKBH5, FTO) and readers (HN-RNPA2B1, HNRNPC, IGF2BP2/3, YTHDC1, YTHDF2) among three RCC types ccRCC, pRCC and chRCC (Supplementary Fig. 1).Through Submap algorithm, we predicted the responses to anti-PD-1 and anti-CTLA4 therapy.As a result, among three RCC types, ccRCC patients were more likely to respond to anti-PD-1 therapy.There was no significant difference in the response to anti-CTLA4 therapy among three RCC types (Supplementary Fig. 1).By univariate-cox regression analyses, prognostic implication of each m 6 A regulator was evaluated across RCC patients.METTL14, ZC3H13, FTO, LRPPRC and YTHDC1 were protective factors of RCC, while HNRNPA2B1 and IGF2BP1/2/3 acted as risk factors of RCC (Fig. 1F).Further multivariate cox regression analysis showed that METTL14 and LRPPRC were protective factors, while IGF2BP3 and HNRNPA2B1 were identified as risk factors, which was consistent with the results of univariate-cox regression analysis (Supplementary Fig. 2).Additionally, we evaluated the relationships between m 6 A regulators and DFS, DSS, and PFS outcomes.Our results demonstrated that the m 6 A regulators IGF2BP1/3 were significantly associated with DFS, DSS, and PFS, indicating that they might contribute to RCC progression (Fig. 1G-I).The PPI network revealed the tight interactions among the m 6 A regulators (Fig. 1J).Also, there were significant correlations between the m 6 A regulators at the mRNA levels across RCC samples (Fig. 1K).Above data highlighted the important implications of the interactions of m 6 A regulators in RCC progression.

Protein Expression of Prognostic m 6 A Regulators and Interactions of M 6 A Regulators with Oxidative Stress in RCC
Using the human protein altas database, we verified the expression of prognostic m 6 A regulators in RCC and normal kidneys.As depicted in Fig. 2A, METTL14, ZC3H13, HNRNPA2B1 were lowly expressed in RCC than normal kidneys.Meanwhile, FTO, LRPPRC, YTHDC1, IGF2BP1/2/3 were highly expressed in RCC compared with normal kidneys.We also evaluated the interactions of m 6 A regulators with 32 oxidative stress-related genes in RCC.In Fig. 2B, m 6 A regulators was significantly linked to most oxidative stress-related genes, indicating that the activity of oxidative stress might be modulated by m 6 A modification during RCC progression.

The M 6 A Modification Subtypes with Different Prognosis, Activation of Pathways, Cancer Stemness and Somatic Copy Number Alterations (SCNA)
By performing consensus clustering analysis, we classified 893 RCC samples into five m 6 A modification subtypes based on the expression profile of the m 6 A regulators (Fig. 3A), named m 6 A modification subtype A (N =   lowest expression of the m 6 A regulators, followed by subtype D (Fig. 3D).Most m 6 A regulators were up-regulated in subtype A, B and E. To validate the reliability of this m 6 A modification classification, the E-MTAB-1980 dataset was adopted.Consistently, ccRCC patients were clearly classified into five m 6 A modification subtypes (Supplementary Fig. 3).The notable survival difference among the subtypes was also proven (Supplementary Fig. 3).The activation of the major biological processes was quantified in each RCC specimen vis ssGSEA.Almost all pathways were inactivated in subtype C and D (Fig. 3E).However, most pathways were activated in subtype A, B and E. Cancer stemness was quantified by two indices: mDNAsi and mRNAsi.In Fig. 3F, subtype B had the lowest mDNAsi value while subtype E had the highest mDNAsi value.But mRNAsi value was the lowest in subtype E (Fig. 3G).The somatic copy number alteration (SCNA) was observed in each RCC specimen.Among subtypes, subtype A exhibited the highest SCNA level (Fig. 3H).

The M 6 A Modification Subtypes with Distinct Oxidative Stress and Tumor Immunity
By ssGSEA method, we quantified the infiltration levels of immune cells in RCC tissues.As a result, subtype A had the highest infiltration levels of most tumor-infiltrating immune cells, but subtype B displayed the lowest infiltration levels of most immune cells (Fig. 4A,B).With the same approach, the activity of oxidative stress was quantified.The significant heterogeneity in oxidative stress was found among five subtypes (Fig. 4C).Among them, subtype B had the lowest activity of oxidative stress.The mRNA expression of immune checkpoints PD-1 and PD-L1 was compared among subtypes.We found that patients in subtype A exhibited the highest mRNA expression of PD-1 (Fig. 4D) and PD-L1 (Fig. 4E).The enrichment levels of each step in the cancer immunity cycle were evaluated across RCC samples.Our results showed that most steps exhibited the highest activation in subtype A, while most steps exhibited the lowest activation in subtype B (Fig. 4F).The overall infiltrations of immune cells and stromal cells were calculated by ESTIMATE method.As expected, the highest immune score and stromal score were found in subtype A (Fig. 4G,H).Meanwhile, we found that there were the lowest immune score and stromal score in subtype B. Furthermore, subtype A displayed the lowest tumor purity while the highest tumor purity was detected in subtype B (Fig. 4I).Microsatellite instability (MSI) was quantified and compared among m 6 A modification subtypes.In Fig. 4J, our data showed that there was the lowest MSI level in subtype E but there was the highest MSI level in subtype B.

Exploration of Biological Implications of m 6 A-Related DEGs
By comparing the DEGs between m 6 A modification subtypes, we finally identified 733 m 6 A-related DEGs (Fig. 5A; Supplementary Table 1).GO enrichment analyses indicated that these m 6 A-related DEGs were significantly correlated to methylation modification processes such as histone modification, covalent chromatin modification and histone H3-H4 methylation (Fig. 5B).Meanwhile, we found that the m 6 A-related DEGs were markedly enriched in mRNA modification pathways (such as mRNA surveillance pathway) and carcinogenic pathways (such as renal cell carcinoma, TGF-beta signaling pathway and central carbon metabolism in cancer; Fig. 5C), which confirmed that m 6 A modification played a nonnegligible role in the tumor progression.Among 733 m 6 Arelated DEGs, 512 were significantly associated with RCC prognosis according to univariate-cox regression analyses (Supplementary Table 2).

Construction of the m 6 A Genomic Subtypes Characterized by Different Survival Outcomes
To further validated the regulation mechanism of m 6 A modification, we conducted unsupervised clustering analyses based on the obtained five m 6 A cluster-related genes.On the basis of the expression profiling of the 512 prognostic m 6 A-related DEGs, 893 RCC samples were clustered into four m 6 A genomic subtypes by unsupervised clustering analyses, named m 6 A genomic subtype A (N = 226), B (N = 318), C (N = 140) and D (N = 195; Fig. 5D).The t-SNE plots confirmed that there was a distinct difference on transcriptome profiles among the four m 6 A genomic subtypes (Fig. 5E).Prognostic analyses showed that there was a significant difference in survival outcomes among m 6 A genomic subtypes (Fig. 5F), where m 6 A genomic subtype B had the most prominent survival advantage and m 6 A genomic subtype A exhibited the poorest survival outcomes.

Oxidative Stress, TME Cell Infiltration and Transcriptome Features in Distinct m 6 A Genomic Subtypes
To explore the biological processes of the four m 6 A genomic subtypes, we used GSVA to analyze the activity of the top 50 signaling pathways.In Fig. 6A, m 6 A genomic subtype B and D displayed the activation of most pathways, while m 6 A genomic subtype A and C had the inactivation of most pathways.Among four m 6 A genomic subtypes, subtype D displayed the highest activity of oxidative stress (Fig. 6B).The roles of the m 6 A genomic subtypes in regulating the TME immune infiltration were further analyzed.

We found that several antitumor immune cells such as activated CD8 T cells, activated CD4 T cells, central memory CD4 T cells, central memory CD8 T cells, effector memory CD4 T cells, effector memory CD8 T cells, activated dendritic cells and natural killer T cells displayed the relatively high infiltration levels in m 6
A genomic subtype A (Fig. 6C).Meanwhile, the most protumor immune cells such as immature dendritic cells, and plasmacytoid dendritic cells exhibited the low infiltration levels in m 6 A genomic subtype A and C as well as the high infiltration levels in m 6 A genomic subtype B and D. The activity of cancer  immunity cycle was also analyzed.As a result, priming and activation of the immune system, recognition of cancer cells by T cells and killing of cancer cells had the highest activity in m 6 A genomic subtype A (Fig. 6D).Above steps exhibited the relatively low activity in m 6 A genomic subtype B and C.This indicated that m 6 A genomic subtype A might be an inflammatory phenotype, while m 6 A genomic subtype B and C might be non-inflammatory phenotypes.Also, we quantified the mRNA expression of immune checkpoints.As expected, PD-1 expression was the highest in m 6 A genomic subtype A, while its expression was the lowest in m 6 A genomic subtype C (Fig. 6E).Furthermore, m 6 A genomic subtype B had the highest PD-L1 expression, while m 6 A genomic subtype C had the lowest PD-L1 expression (Fig. 6F).These data indicated that m 6 A genomic subtype A and B might be more sensitive to anti-PD-1/PD-L1 therapy.

Generation of m 6 A Scoring System and Assessment of the Prognostic Value of m 6 A Score in RCC
Due to the heterogeneity and complexity of m 6 A methylation modification, a m 6 A scoring system was gen-erated for quantifying m 6 A modification subtypes of each RCC patient based on the expression profiling of the prognostic m 6 A-related DEGs that were identified by PCA method using the Boruta algorithm, named m 6 A score (Fig. 7A).For better illustrating the clinical value of m 6 A score, this study analyzed the distribution of m 6 A score in different clinicopathological characteristics (Fig. 7B).Our data showed that patients with age ≥65 had the significantly lower m 6 A score compared to those with age <65.There was a significantly decreased m 6 A score in male patients than female patients.With the increase of grade, stage and T, m 6 A score was gradually reduced across RCC specimens.Metastatic patients had significantly lower m 6 A score than non-metastatic patients.Patients with N1-2 displayed markedly reduced m 6 A score in comparison to those with N0.In addition, we also explored the relationship between clinical factors and m 6 A score.Low m 6 A score was corresponding to more dead cases (p < 0.001), advance stage and grade cases (Supplementary Fig. 4).Based on the median m 6 A score, RCC patients in the TCGA cohort were separated into high and low m 6 A score group.Survival analysis were then presented for comparing the prog- nostic divergence between groups.As a result, patients with high m 6 A score displayed the prominent advantage in OS, DSS and PFI in comparison to those with low m 6 A score (Fig. 7C-E).Consistently, better OS outcomes were ob-served in high m 6 A score group relative to low m 6 A score group in the E-MTAB-1980 dataset (Fig. 7F).The prognostic value of m 6 A score was also externally verified in the two anti-PD-1/PD-L1 immunotherapy cohorts.Both in the GSE78220 dataset and the dataset reported by Liu et al. [37], high m 6 A score indicated undesirable outcomes than low m 6 A score among patients who received anti-PD-1/PD-L1 therapy (Fig. 7G,H).Finally, we estimated the correlation between m 6 A cluster, genomic cluster and m 6 A score using sankey plot (Supplementary Fig. 5).Most of cases from m 6 A cluster E and genomic cluster B were attribute to high m 6 A score, which can partly explain why m 6 A cluster E and genomic cluster B associated with a better survival outcome.

Correlation between m 6 A Score and Signaling Pathways, Tumor Immunity and Oxidative Stress in RCC
The activity of the major signaling pathways was compared in high-and low-m 6 A score RCC patients.We found that most pathways such as carcinogenic pathways and inflammation-related pathways were distinctly activated in high-m 6 A score group compared to low-m 6 A score group.Meanwhile, in comparison to high-m 6 A score group, KRAS signaling, coagulation and myogenesis were distinctly activated in low-m 6 A score group (Fig. 8A).Consistently, GSEA results revealed that apoptosis, focal adhesion, insulin signaling pathway, MAPK signaling pathway, oocyte meiosis, pathway in cancer and WNT signaling pathway were significantly up-regulated in high-m 6 A score samples than low-m 6 A score samples (Fig. 8B).No pathways were significantly enriched in low-m 6 A score specimens.Furthermore, we evaluated the differences in tumor immunity between high-and low-m 6 A score groups.Most immune checkpoints such as LAIR1, HAVCR2 and ICOS displayed the higher mRNA expression in high-m 6 A score group compared to low-m 6 A score group (Fig. 8C).TNFRSF18, TNFSF14, TMIGD2, LAG3 and CD70 were significantly up-regulated in low-m 6 A score group than high-m 6 A score group.There were increased infiltration levels of central memory CD8 T cell, eosinophil, effector memory CD4 T cell, memory B cell, type 2 T helper cell, mast cell, neutrophil, natural killer cell, plasmacytoid dendritic cell and immature dendritic cell in high-m 6 A score group than low-m 6 A score group (Fig. 8D).Meanwhile, activated CD8 T cell, CD56dim natural killer cell, monocyte, MDSC and activated B cell exhibited the higher infiltration levels in low-m 6 A score group in comparison to highm 6 A score group.Most oxidative stress-related genes had higher expression levels in high-than low-m 6 A score group (Fig. 8E).

Drug Sensitivity and Potential Druggable Targets in
High-and Low-m 6 A Score RCC Patients IC50 values of four approved targeted drugs including axitinib, pazopanib, sorafenib and sunitinib were estimated in each RCC sample.The differences in IC50 values were compared between high-and low-m 6 A score RCC patients.In comparison to low-m 6 A score group, there were significantly reduced IC50 values of axitinib, pazopanib and sorafenib in high-m 6 A score group (Fig. 9A-C).This indicated that RCC patients with high-m 6 A score were more sensitive to axitinib, pazopanib and sorafenib.Conversely, lower IC50 value of sunitinib was found in low-m 6 A score group compared to high-m 6 A score group (Fig. 9D), demonstrating that low-m 6 A score patients were more likely to benefit from sunitinib.Potential therapeutic agents were predicted for low-m 6 A score patients.Two drug dataset (CTRP and PRISM) were adopted to identifying candidate CTRP and PRISM-derived agents with higher drug sensitivity in low-m 6 A score patients.Differential drug response analyses between high-and low-m 6 A score groups were carried out for predicting drugs with reduced AUCs in low-m 6 A score patients.As a result, six CTRPderived compounds (CR-1-31B, daporinad, leptomycin B, methotrexate, oligomycin A and ouabain; Fig. 9E) and ten PRISM-derived compounds (bosutinib, gemcitabine, GSK1070916, halobetasol-propionate, mesna, OTX015, PD-0325901, Ro-4987655, romidepsin and vincristine; Fig. 9F) had reduced AUCs in low-m 6 A score group compared to high-m 6 A score group, indicating that they could become potential therapeutic agents against low-m 6 A score patients.Additionally, potential druggable targets were predicted for low-m 6 A score patients.As a result, we found that m 6 A score was positively associated with CERES scores of four druggable targets (NMU, MOCOS, UCHL1, and RARRES1), with negative relationships with their protein expression (Fig. 9G,H).Therefore, NMU, MOCOS, UCHL1, and RARRES1 might become potential druggable targets for low-m 6 A score patients.

Discussion
RCC contains distinct malignancies with different pathologic characteristics and different molecular pathways [41].Due to no obvious symptoms, 30% RCC patients have metastasis at diagnosis.Typically, most patients experience undesirable survival outcomes [41].RCC initiation and development involve diverse and complex mechanisms.Hence, based on the precise mechanisms underlying RCC, development of novel therapeutic strategies is of importance in clinical practice.M 6 A methylation modifications mediated by methyltransferases, demethylases and binding proteins act as critical determinants in mRNA metabolism [42].Emerging evidence suggests that m 6 A methylation modifications occupy 80% of RNA methylation modifications, which may mediate diverse malignancy-associated processes like tumorigenesis, metastasis, and immune escape [43].For example, m 6 A methyltransferase METTL14 mediates tumor immune and development of ccRCC [44].In this study, we found that most m 6 A regulators displayed abnormal expression in RCC compared with normal kidneys.Furthermore, there were remarkable differences in the expression of most m 6 A regulators between the early and late stage of RCC, demonstrating the dynamic process of m 6 A modification.The close interactions between the m 6 A regulators were found according to PPI and pearson correlation analyses, indicating that their interactions might participate in mediating RCC progression [45].
Here, based on the expression matrix of 23 m 6 A regulators, we conducted five distinct m 6 A modification subtypes with different biological functions, survival outcomes, oxidative stress, TME features, cancer stemness and  tumor immunity across RCC patients.Furthermore, a m 6 A scoring system was developed for quantifying the m 6 A modification patterns for each RCC patient, which might be applied for predicting prognosis, oxidative stress, immunotherapy responses, sensitivity to chemotherapy drugs and designing novel agents.This scoring system possessed the advantages of focusing the scores on the set with the largest block of well associated (or anti-associated) genes in the set; meanwhile, down-weighting contributions from genes which do not track other set members.High m 6 A score indicated favorable survival outcomes of RCC patients.Selecting RCC subjects who will respond specifically to ICIs is still a challenge.Distinctive biomarkers are urgently needed to predict clinical outcomes and responses to ICIs.The prognostic value of m 6 A score was confirmed in two anti-PD-1/PD-L1 therapy cohorts.Most immune checkpoints exhibited the up-regulation in RCC patients with high m 6 A score, indicating that the group of patients might benefit from ICIs.
Molecularly targeted treatment is prone to drug resistance [41].Here, we found that high m 6 A score indicated higher sensitivity to axitinib, pazopanib and sorafenib, while low m 6 A score were more sensitive to sunitinib.This indicated that m 6 A score might be applied for estimating the drug sensitivity for RCC patients.Due to poor survival outcomes of patients with low m 6 A score, we predicted six CTRP-derived compounds (CR-1-31B, daporinad, leptomycin B, methotrexate, oligomycin A and ouabain) and ten PRISM-derived compounds (bosutinib, gemcitabine, GSK1070916, halobetasol-propionate, mesna, OTX015, PD-0325901, Ro-4987655, romidepsin and vincristine) that could become potential therapeutic agents against low-m 6 A score patients.Nevertheless, more experiments will be required for verifying the therapeutic efficacy of above agents in RCC.In addition, the immunotherapy cohorts of RCC patients are lack.Therefore, we only validated the significance of m 6 A score in predicting the response of anti-PD-1/PD-L1 immunotherapy in the GSE78220 dataset and the dataset reported by Liu et al. [37].Yet, the immunotherapy results were not in line with our prior results, likely due to the tumor heterogeneity and limited sample size of these datasets.Thus, in our future studies, the m 6 A modification subtypes and score will be validated in a prospective and large number cohort.

Conclusions
Collectively, m 6 A regulator-based methylation patterns and quantification of m 6 A score exerted key functions in prognosis stratification, oxidative stress, TME, tumor immunity and immunotherapy responses, which may assist clinicians to achieve individualized therapy for RCC patients.

Fig. 1 .
Fig. 1.Systematic analyses of genetic variation, expression, prognostic implication and interactions of N6-methyladenosine (m6A) regulators in renal cell carcinoma (RCC).(A,B) Before and after removing batch effects by integrating three types of RCC: clear cell RCC (ccRCC), papillary RCC (pRCC) and chromophobe RCC (chRCC) from the Cancer Genome Atlas (TCGA) datasets.PCA, Principal component analysis.(C) The frequency of copy number variation (CNV) (red: gain and blue: loss) of m 6 A regulators across RCC samples.(D) An overview of the mRNA expression of m 6 A regulators in RCC and normal specimens.p values were calculated by unpaired student's t test, *p < 0.05 and ****p < 0.0001.Red: up-regulation and green: down-regulation.(E) The mRNA expression of m 6 A regulators in stage I&II and stage III&IV RCC cases.*p < 0.05, ***p < 0.001 and ****p < 0.0001.(F-I) Univariate-cox regression analyses of the correlation between m 6 A regulators and RCC patients' overall survival (OS), disease free survival (DFS), disease-specific survival (DSS), and progression free survival (PFS).(J) The protein-protein interaction (PPI) network of writers (green), erasers (blue) and readers (red).(K) The Pearson correlations between m 6 A regulators at the mRNA levels across RCC samples.Red: positive correlation and green: negative correlation, *p < 0.05, **p < 0.01 and ***p < 0.001.

Fig. 2 .
Fig. 2. Protein expression of prognostic m 6 A regulators and interactions of m 6 A regulators with oxidative stress in RCC.(A) Validation of the expression of prognostic m 6 A regulators in RCC and normal kidneys from the Human Protein Altas database.Scale bar, 200 µm.(B) Heatmap of the interactions of m 6 A regulators and oxidative stress-related genes in RCC.Red: positive correlation and green: negative correlation, *p < 0.05; **p < 0.01.

Fig. 3 .
Fig. 3. Construction of the m 6 A modification subtypes with different prognosis, activation of pathways, cancer stemness and somatic copy number alteration (SCNA) across RCC samples.(A) Consensus clustering analyses for identifying the best clustering number (k = 5) across RCC samples based on the expression matrix of m 6 A regulators. (B) The t-SNE for assessing the clustering accuracy across RCC samples.(C) Survival analyses of RCC patients in different m 6 A modification subtypes with log-rank tests.(D) An overview of the mRNA expression of m 6 A regulators in each m 6 A modification subtype.Red: up-regulation and green: downregulation.(E) An overview of the activation levels of the 50 main signaling pathways in RCC specimens from five m 6 A modification subtypes.(F) The mDNAsi score, (G) mRNAsi score and (H) SCNA level among different m 6 A modification subtypes.p values were determined with Kruskal-Wallis test.

Fig. 4 .
Fig. 4. Assessment of the m 6 A modification subtypes with distinct oxidative stress and tumor immunity across RCC samples.(A,B) The infiltration levels of immune cells among five m 6 A modification subtypes.(C) Activity of oxidative stress among five m 6 A modification subtypes.(D,E) The mRNA expression of immune checkpoints programmed cell death 1 (PD-1) and programmed cell death ligand 1 (PD-L1) among five m 6 A modification subtypes.(F) The activation of each step in the cancer immunity cycle among five m 6 A modification subtypes.(G-I) Immune score, stromal score, and tumor purity among five m 6 A modification subtypes.(J) The microsatellite instability (MSI) level among five m 6 A modification subtypes.p values were calculated with Kruskal-Wallis test, ns: not significant; *p < 0.05; **p < 0.01; ***p < 0.001.

Fig. 5 .
Fig. 5. Establishment of the m 6 A genomic subtypes with different survival outcomes across RCC samples.(A) Venn diagram for the m 6 A-related differentially expressed genes (DEGs) by comparing the DEGs between m 6 A modification subtypes.(B,C) Gene ontology (GO) and kyoto encyclopedia of genes and genomes (KEGG) enrichment results of the m 6 A-related DEGs.The length of the rectangle indicated the number of enriched genes in a specific term.The more the color tended to red, the smaller the adjusted p value.(D) The four m 6 A genomic subtypes across RCC samples based on the m 6 A-related DEGs via unsupervised clustering analyses.(E) The t-SNE for the transcriptome profiles among the four m 6 A genomic subtypes.(F) Kaplan-Meier survival curves for the four m 6 A genomic subtypes across RCC samples.p values were determined with log-rank tests.

Fig. 6 .
Fig. 6.The m 6 A genomic subtypes characterized by oxidative stress, tumor microenvironment (TME) cell infiltration and transcriptome features across RCC samples.(A) An overview of the activity of the major 50 signaling pathways in RCC specimens from the four m 6 A genomic subtypes.Red: activation and green: inactivation.(B) Activity of oxidative stress in the four m 6 A genomic subtypes.(C) The enrichment scores of tumor-infiltrating immune cells in the four m 6 A genomic subtypes.(D) The activity of each step in the cancer immunity cycle in the four m 6 A genomic subtypes.(E,F) The mRNA expression of immune checkpoints PD-1 and PD-L1 in the four m 6 A genomic subtypes.p values were calculated with Kruskal-Wallis test, ns: not significant; *p < 0.05; **p < 0.01; ***p < 0.001.

Fig. 7 .
Fig. 7. Generation of m 6 A scoring system and assessment and external validation of the prognostic value of m 6 A score for RCC patients.(A) Identification of the final m 6 A-related DEGs that were used for construction of m 6 A score by principal component analysis (PCA) method utilizing Boruta algorithm.(B) The distribution of m 6 A score in different clinicopathological characteristics across RCC patients, including age, sex, stage, grade and T, N, M. p values were determined with Wilcoxon or Kruskal-Wallis tests.(C-E) Kaplan-Meier curves of OS, DSS and PFI for RCC patients with high and low m 6 A score in the TCGA cohort.(F) Kaplan-Meier curves of OS between high and low m 6 A score ccRCC patients in the E-MTAB-1980 cohort.(G,H) Kaplan-Meier curves of OS between high and low m 6A score groups for patients treated with anti-PD-1/PD-L1 therapy in the GSE78220 dataset and the dataset reported by Liu et al[37].

Fig. 8 .
Fig. 8. Correlation between m 6 A score and signaling pathways, tumor immunity and oxidative stress in RCC samples.(A) An overview of the correlation between m 6 A score and the enrichment score of the major 50 signaling pathways in RCC samples.The m 6 A scores from left to right gradually increased.Red indicated high activation of a specific pathway and green indicated low activation of a specific pathway.(B) Gene set enrichment analysis (GSEA) for the activated signaling pathways in high-m 6 A score samples.(C) An overview of the correlation between m 6 A score and the mRNA expression of immune checkpoints in the RCC cohort.Red represented high expression and green represented low expression for a specific immune checkpoint.(D) An overview of the correlation between m 6 A score and the infiltration levels of immune cells in RCC specimens.Red indicated high infiltration level and green represented low infiltration level for a specific immune cell.(E) The mRNA expression of oxidative stress-related genes in high and low m 6 A score groups.Comparisons between two groups were analyzed with Wilcoxon test, *p < 0.05; **p < 0.01; ***p < 0.001.

Fig. 9 .
Fig. 9. Comparison of the response to drugs in high-and low-m 6 A score RCC patients.(A-D) Violin plots of the estimated IC50 values of chemotherapy drugs including axitinib, pazopanib, sorafenib and sunitinib in high-and low-m 6 A score RCC patients.(E) Differential drug response analyses of six Cancer Therapeutics Response Portal (CTRP)-derived compounds in high-vs low-m 6 A score RCC groups.Comparisons between two groups were analyzed with Wilcoxon test, ***p < 0.001.(F) Differential drug response analyses of ten PRISM-derived compounds in high-vs low-m 6 A score RCC groups.Comparisons between two groups were analyzed with Wilcoxon test, ***p < 0.001.(G,H) Associations between m 6 A score and CRISPR-Cas9 essentiality screens (CERES) score and protein expression of druggable targets.