Comprehensive Evaluation of the m6A Regulator Prognostic Risk Score in the Prediction of Immunotherapy Response in Clear Cell Renal Cell Carcinoma

Background Clear cell renal cell carcinoma (ccRCC) is known for its high drug resistance. The tumor-immune crosstalk mediated by the epigenetic regulation of N6-methyladenosine (m6A) modification has been demonstrated in recent studies. Therefore, m6A modification-mediated immune cell infiltration characteristics may be helpful to guide immunotherapy for ccRCC. Methods This study comprehensively analyzed m6A modifications using the clinical parameters, single-cell RNA sequencing data, and bulk RNA sequencing data from the TCGA-ccRC cohort and 13 external validation cohorts. A series of bioinformatic approaches were applied to construct an m6A regulator prognostic risk score (MRPRS) to predict survival and immunotherapy response in ccRCC patients. Immunological characteristics, enriched pathways, and mutation were evaluated in high- and low-MRPRS groups. Results The expressional alteration landscape of m6A regulators was profiled in ccRCC cell clusters and tissue. The 8 regulator genes with minimal lambda were integrated to build an MRPRS, and it was positively correlated with immunotherapeutic response in extent validation cohorts. The clinicopathological features and immune infiltration characteristics could be distinguished by the high- and low-MRPRS. Moreover, the MRPRS-mediated mutation pattern has an enhanced response to immune checkpoint blockade in the ccRCC and pan-cancer cohorts. Conclusions The proposed MRPRS is a promising biomarker to predict clinical outcomes and therapeutic responses in ccRCC patients.


INTRODUCTION
Clear cell renal cell carcinoma (ccRCC) is the most common type of renal cancer and accounts for nearly 3% of adult malignant tumors (1). Approximately 30% of patients already have advanced ccRCC or metastases when they are first diagnosed, and have missed the opportunity for surgical intervention (2). Although targeted therapy and immunotherapy have become the main adjuvant therapy for advanced ccRCC, the complete response rate and partial response rate remained low (3,4). So far, the biomarkerbased therapeutic strategies for advanced ccRCC have been missing.
N6-methyladenosine (m 6 A) modification is an important factor for messenger RNA (mRNA) stability, splicing, and translation (5)(6)(7). Serval m 6 A-sequencing studies have revealed that abnormal m 6 A regulatory enzymes are involved in mutagenesis, proliferation, and tumorigenesis through the dysregulation of the m 6 A pathway (8,9). Recently, m 6 A modifications have been shown to play a role in the regulation of immune cells, such as the following: METTL3-mediated m 6 A modification increased the translation of certain immune transcripts and physiologically promoted the activation of dendritic cells (DCs) and DC-based T-cell responses (10), and ALKBH5 regulated m 6 A modification in the 3'UTR region of PD-L1 mRNA and inhibited the expansion and cytotoxicity of T cells by sustaining tumor cell PD-L1 expression (11). The potential relationship between RNA m 6 A dysregulation and tumor-infiltrating immune cells (TIICs) has motivated us to investigate and find the potential biomarkers for predicting immune checkpoint therapy outcomes. Herein, we systematically evaluated the m 6 A regulator-based risk score and its associated gene mutation with the TIICs and revealed a new predictive method that could be used to predict the immunotherapy response in ccRCC and pan-cancer.

Data Collection and Processing
The RNA sequencing (RNA-seq) transcriptome data of patients with ccRCC and the corresponding clinical data and mutation profiles were downloaded from The Cancer Genome Atlas (TCGA) database. The validation datasets (GSE53757, GSE40435, GSE29609, and E-MTAB-3267) were included for analysis from the Gene Expression Omnibus (GEO) database and the European Molecular Biology Laboratory (12)(13)(14)(15). The relative transcriptomic and clinical data of three immunotherapeutic cohorts of patients with ccRCC were obtained from the online supplementary data (16)(17)(18). RNA-seq data of the immunotherapy cohort of bladder cancer (19) and melanoma (PRJEB23709 and phs000452) were collected for testing (20,21). The annotated response and mutational data of patients from a discovery cohort receiving ICB treatment from 4 studies were collected and consolidated to study the relationship between mutated genes and immunotherapy (17,(22)(23)(24). The single-cell dataset of ccRCC ICB treatment was obtained from PMID33861994 (25). The information for all collected data is presented in Table S1.

Single-Cell RNA Sequencing Analysis
The association with m 6 A regulators was established by analyzing the genes related to the immune response in the scRNA-seq results of ccRCC (25). The CellRanger software (version 5.0.0) and STAR were used for preprocessing. Principal component analysis (PCA) was run using the "RunPCA" function on the variable genes identified, and the k-nearest neighbor graph was constructed by the "FindNeighbors" function. Uniform manifold approximation and projection (UMAP) was used to visualize single-cell transcriptional profiles and clusters. Marker genes were visualized on UMAP plots using log-normalized counts.

Cell-Cell Communication Analysis
CellPhoneDB applies an algorithm that considers only receptors and ligands with broad expression among the tested cell types, followed by calculating the likelihood of cell-type specificity of a given receptor-ligand complex with a sufficient number of permutations (26).

Construction and Validation of the m 6 A Gene Signature
The significant m 6 A RNA methylation regulators were established by the least absolute shrinkage and selection operator (LASSO) Cox regression (with the penalty parameter estimated by 20-fold crossvalidation). Those regulator genes with minimal lambda were integrated to build an MRPRS, and it was developed according to the expression level using univariate Cox. The "glmnet" package was used to perform the LASSO Cox regression model analysis.
The limma R package's empirical Bayesian approach was applied to determine differentially expressed genes (DEGs) between high and low m 6 A scores. The significance criteria for determining DEGs were set as the 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. A protein-protein interaction (PPI) network was constructed by STRING (https:// string-db.org/) and evaluated using the Cytoscope software (31).

Statistical Analysis
Statistical tests were carried out using R version 4.0.4, SPSS 25.0 (IBM, NY, USA) and GraphPad Prism 8.0. The expression levels of the m 6 A RNA regulators were compared with the Mann-Whitney U test in ccRCC versus normal tissues. Survival curves were generated using the Kaplan-Meier method, and the difference was compared with the log-rank test. Pearson correlation coefficient was used to compare the correlation between MRPRS and gene expression values. The "oncoplot" function of the R package "maftools" was used to determine the mutation landscape of the TCGA ccRCC cohort and immunotherapeutic cohort. The high-and low-group was divided based on the optimal cut-off value calculated by the function "surv_cutpoint" in the R package "survminer". All the R package used in this study is listed in Table S2. p < 0.05 indicated statistical significance.
We next used the scRNA-seq data (25) (PMID33861994) to evaluate the expression of m 6 A regulators in different subsets of ccRCC cells. A total of 65,535 cells were divided into 7 cell clusters, a total of 6,539 epithelial cells from ccRCC multiple regions (Near, Far, Center, and Lymph node) and normal tissues of 6 ICB-treated and untreated patients were extracted, and 6 ccRCC cell clusters were identified based on AQP3, GPX3, CCNI, STMN1, VCAM1, and VCAN expression (Figures 1B, C; Figures S1C, D; Table S5). The genes with the most significant differential expression in each cell cluster were described in the heatmap ( Figure 1E; Table S6), and significant functional heterogeneity was found among the 7 cell clusters ( Figure 1F). As shown in Figure 1G, m 6 A regulators also presented the expression heterogeneity between normal cell clusters and tumor cell clusters, as well as between the 6 ccRCC cell clusters. Moreover, the expression of WTAP, YTHDC1, YTHDC2, HNRNPC, and HNRNPA2B1 was significantly different between the ICB-resistant-related GPX3 + epithelial cells and ICB-responserelated VCAM1 + epithelial cells ( Figure 1D), indicating that the differential expression level of m 6 A regulators in ccRCC tissue and cell clusters may be related to the efficacy of immunotherapy.

Construction, Validation, and Immunotherapy Response Evaluation of the m 6 A Regulator Prognostic Risk Score
To systematically evaluate these differences in m 6 A regulators, the MRPRS was established by the LASSO Cox algorithm, the 8 regulator genes with minimal lambda were integrated to build an MRPRS (Figures 2A, B) ( Figure S2A), and it was developed according to the expression level using univariate Cox ( Figure  S2B). The ccRCC patients in the TCGA database were divided into the high-MRPRS group (N = 134) and the low-MRPRS group (N = 395) based on the optimal cut-off value calculated by the function "surv_cutpoint" in the R package "survminer", and the ccRCC patients in the high-MRPRS group had a significantly shorter overall survival time than that in the low-MRPRS group (p < 0.001; Figure 2C). The prognostic value of MRPRS was also validated in an independent cohort (GSE29609, p = 0.037; Figure 2D). We continued to extend the MRPRS signature to 16 other tumor types, such as esophageal carcinoma (ESCA), lung adenocarcinoma (LUAD), and pancreatic adenocarcinoma (PAAD) ( Figure S2C). These results present that MRPRS is negatively associated with survival outcomes.
Next, we investigated the correlation between MRPRS and immunotherapy response in three independent ccRCC cohorts (PMID29301960, PMID32472114, and PMID32895571) (16,18,22), and found that the MRPRS was significantly higher in the response group than in the non-response group; the high-MRPRS group presented a markedly prolonged survival ( Figures 2F-I).
Moreover, the increased MRPRS in the VCAM1 + cell cluster was positively correlated with the patients who experienced complete and mixed responses ( Figure 2E). Similar results were also obtained in the extended dataset of bladder cancer (IMvigor210) and melanoma (PRJEB23709) (Figures S3A-E).
In addition, we analyzed the expression of targeted therapyand chemotherapy-related genes between high-and low-MRPRS groups ( Figure 2J). Interestingly, VEGF and mTOR pathway-related genes were found to be highly expressed in the low-MRPRS group, and MRPRS was negatively correlated with the expression level of angiogenesis-related genes including PECAM, FLT1/4, VWF, and CDH5 ( Figures 2J, K). In the E-MTAB-3267 cohort of ccRCC patients treated with sunitinib, the MRPRS was significantly lower in the response group than in the non-response group, and the low-MRPRS group showed a markedly prolonged survival (Figures 2L, M). Collectively, our data suggest that the patients with high MRPRS may benefit from immunotherapy and those with low MRPRS may benefit from targeted therapy.

The Clinicopathological Features and Immune Infiltration Characteristics in Distinct MRPRS
We examined the correlation between the MRPRS and the clinical parameters. No significant association was found between the MRPRS and gender and age, but significant associations in terms of TNM stages, grade, and survival status were observed ( Figure S4A), and MRPRS is positively correlated with TNM stages and grades (p < 0.001; Figure 3A). Moreover, each of the four different T stages including stage I had significantly higher MRPRS when compared with the control subjects, and MRPRS in the metastasis group was significantly higher than that in the non-metastasis group (p < 0.001, respectively; Figure 3A).
To investigate the effects of MRPRS on the immune infiltration characteristic of ccRCC, we evaluated the expression of immunomodulators and the infiltration levels between high-and low-MRPRS groups in ccRCC, as shown in Figure 3B; 5 immunomodulators (chemokine, receptor, immunostimulator, inhibitory immune checkpoint, and MHC) and the infiltration levels of 4 types of TIICs (CD8+ T cells, DC, macrophages, and Th1 cells) were positively correlated with the high-MRPRS group (p < 0.05). The MCP counter, xCell, CIBERSORT, and ssGSEA algorithm were used to calculate an immune score and to estimate the abundance of various types of immune cells. We found significantly higher estimates of Tregs, CD8 + T cells, NK cells, and B cells in ccRCC with  Table S7). Moreover, the VCAM1 + cell cluster presented upregulation of HLA-A and HLA-B ( Figure 1E), and it was also the cell cluster that communicates most frequently with immune cells ( Figure 3E; Table S8).
We performed volcano plots based on the DEGs from the high-and low-MRPRS groups. The results of the volcano plots showed that 1,780 genes were significantly upregulated in the comparison of the high-and low-MRPRS groups (Table S9). In the PPI network from the STRING database with the Cystoscope software, we constructed a co-expression network consisting of 45 nodes and 169 edges ( Figure S4C). These included immunerelated genes, CD19 and CD79A, and membrane proteins on the surface of B cells, which participate in the proliferation and differentiation of B cells. FOXP3 and IL2RA (CD25) are the characteristic markers of Treg cells. We also found that the expression of many costimulatory factors, such as TNFSF14, TNFRSF18, and a large amount of interleukins such as IL2 and IL6, promotes T-cell proliferation and T-cell-mediated killing ( Figure S4B). GO enrichment analysis and KEGG analysis of these signature genes revealed that these DEGs were enriched in several biological processes and pathways related to immune regulation ( Figures S4D, E; Table S10). Moreover, the number and diversity of T-cell receptors (TCRs) were higher in the high-MRPRS group than in the low-MRPRS group (p < 0.01) ( Figure 3D). These findings suggest that the regulation of TCR gene expression may be influenced by the specific tumor cell cluster with abnormal m 6 A modifications.

The Landscape of Genetic Variation of MRPRS Groups in ccRCC
The somatic mutation profile between the high-and low-MRPRS groups in the TCGA-ccRCC cohort used the maftools package, and the top 10 most frequently mutated genes in each group are shown in Figure 4A. Notably, SETD2, TRIOBP, RYR2, ZFPM2, and ABCC6 occupy the top 5 positions among differently mutated genes between the high-and low-MRPRS group (Figure 4B), and a lollipop plot showed the different mutation spots of these mutated genes between two groups ( Figure 4C). Interestingly, the mutation rate of SETD2 was 23.18% in the high-MRPRS group and 3.87% in the low-MRPRS group, and the remaining four genes were mutated only in the high-MRPRS group. In addition, the distribution of variants according to variant classification, variant type, and singlenucleotide variant (SNV) class was displayed as a cohort summary plot, and among all the genomic alterations, missense mutations were the predominant type, with C>T and C>G representing the most common SNV classes (Figures S6A, B). Somatic mutation  gene interaction networks showed a high correlation between VHL and PBRM1, PBRM1 and SETD2, and TTN and MUC16 in the high MRPRS score group (Figures S6C, D). We also applied the MCP counter and ssGSEA algorithm to estimate the tumor-infiltrating immune cells between the group defined by patients with at least one mutation in these five genes or without mutation. As shown in Figures 4D, E, the infiltrating immune cells of T cells, DC cells, and B cells in the mutation status group were higher than those in the non-mutation status group (p < 0.01).

The Role of the MRPRS-Mediated Mutation Pattern in Predicting the Response to Immunotherapy
We next investigated whether the MRPRS-mediated mutation pattern could predict patients' response to immunotherapy. We constructed a pan-cancer cohort with anti-PD-1/PDL1 immunotherapy consisting of 1,959 cases based on four cohorts (17,(22)(23)(24) (Table S11), and patients with mutation exhibited a significantly clinical response to immunotherapy and markedly prolonged survival in ccRCC ( Figures 5A, B). Immunotherapy represented by PD-L1 and PD-1 blockade is a breakthrough in tumor therapy. We continued to extend the potential role of MRPRS-mediated mutation pattern in predicting responses to immunotherapy in pan-cancer ( Figure 5C) and revealed that the OS and PFS in patients with mutations were significantly higher than in those without mutations (Figures 5D, E). However, the MRPRS-mediated mutation pattern had no significance in OS of either TCGA-ccRCC or TCGA-pan-cancer ( Figures S7A, B); by contrast, the PFS of the mutation group was worse than that of the non-mutation group in the TCGA-ccRCC (p = 0.049, Figure S7C), and the PFS in the TCGA-pan-cancer was not significant ( Figure S7D).
In addition, tumor mutation burden (TMB) may serve as a biomarker for predicting the response to ICB treatment. We next divided the pan-cancer cohort patients into three groups according to TMB and MRPRS-mediated mutations and found that the OS of patients with low TMB and mutations was significantly better than that of the patients with high TMB and the patients with low TMB and non-mutations (P < 0.0001, Figure 5F).

DISCUSSION
The TME of ccRCC is known to be highly immunosuppressive (32). In the TME, T cells are continuously exposed to antigens, which leads to the impairment of T-cell function and ultimately to a dysfunctional state called "exhaustion" (33). The use of monoclonal antibodies or small molecules to reverse T-cell exhaustion is the basic strategy of immunotherapy (34). Since the results of the Checkmate-025 study, the immunotherapy of ccRCC has been the focus of attention, and now, combined targeted and immunotherapy has become a key component of the adjuvant treatment of advanced ccRCC (35). However, the complete or mixed response rate of immunotherapy in ccRCC is still low. Relying on biomarkers to screen patients who benefit from immunotherapy and to avoid overtreatment has long been expected in clinical practice. Increasing evidence has demonstrated that m 6 A modification plays an indispensable role in immunity, inflammation, and therapy resistance through various m 6 A regulators (36). In this study, we systematically evaluated the expression level of m 6 A regulators in the ccRCC tissue and cell clusters and focused on the detailed relationship between m 6 A modification and TME to enhance our understanding of the ccRCC-immune crosstalk. We constructed an MRPRS comprising 8 m 6 A regulators by the LASSO algorithm to provide reliable biomarkers able to predict the prognosis and immunotherapy efficacy. For the first time, we analyzed the MRPRS levels in ccRCC cell clusters and found that the increased MRPRS in the VCAM 1 + cell cluster was positively correlated with patients who experienced complete and mixed responses. This is consistent with our finding that the positive correlation between MRPRS and immunotherapy benefits the bulk tissue datasets. It is interesting to note that the spatial localization of this immunotherapy-related ccRCC cell cluster is worthy of further investigation. We further explored the detailed role of m 6 A modification in modifying immune characteristics in ccRCC. The results of the GO and KEGG pathway analyses revealed a significant enrichment of genes in immune-related pathways. GO enrichment analysis showed that these DEGs were enriched in the humoral immune response, immunoglobulin complex, and antigen binding. The results of the KEGG analysis indicated these enriched pathways such as neuroactive ligand-receptor interaction, cytokine-cytokine receptor interaction, and the calcium signaling pathway. These results indicated that DEGs in the ccRCC are enriched in immune-related genes distinguished by the MRPRS. Among these DEGs, numerous immune-related genes were found, such as CD19, CD79A, FOXP3, CXCL13, IL2, and TNFRSF13B. FOXP3 is a hallmark of regulatory T cells, CXCL13 is related to CD8 T cells, and CD79A, CD19, and TNFRSF13B are markers of neoplastic B cells. This was in accordance with results from the single-cell sequencing analysis of ccRCC (25). Moreover, these immune cells comprise the main part of tertiary lymphoid structures (TLSs), which have recently been associated with effective antitumor immune responses in cancer patients (37,38). These findings suggest that m 6 A modification may influence the formation of tertiary lymphatic structures.
The patient with a high MRPRS has a poor prognosis, and this could be due to the observation that several critical inhibitor immune checkpoints were significantly highly expressed in the high-MRPRS group, which may limit cytotoxic immune cell activities in the TME, such as CD8 T cells, causing cytotoxic cells to be in an exhausted functional state (39). Therefore, patients with high MRPRS may be more sensitive to immunotherapy. Several studies have also demonstrated that inflammatory tumor phenotypes are more sensitive to ICB (40,41). We next compared the prognostic value of the MRPRS based on ccRCC immunotherapeutic cohorts, and the high-MRPRS group presented a prolonged survival. These findings suggest that MRPRS could be used as a new predictive biomarker for immunotherapy response in ccRCC. Furthermore, we identified 5 genes (SETD2, TRIOBP, RYR2, ZFPM2, and ABCC6) that show the most significant differences in the comparison of mutated genes between two MRPRS groups. We found that the patients with mutated genes had worse PFS outcomes than the non-mutated group, and this was consistent with the high-MRPRS group showing worse survival than those with lower MRPRS. The potential association of TMB with sensitivity to ICB is based on the hypothesis that in tumors with high TMB, there is an increased production of surface neoantigens, thus stimulating the anti-tumor immune system response (42). The TMB has been investigated in several tumor settings, mainly in NSCLC and melanoma, as a stratification marker to predict the response to immune agents, showing promising yet inconclusive results (43,44). In contrast, it has also been reported that high TMB fails to predict immune checkpoint blockade response across all cancer types (45). Herein, we applied the prediction of immunotherapeutic efficacy with the MRPRS-mediated mutation pattern and TMB in pan-cancer cohort and found that the MRPRS-mediated mutation pattern was a better predictor of immunotherapy outcome than the TMB. The regulatory relationship between the m6A modification and the gene mutation still needs to be studied. Consequently, we provided a new perspective on the immune characteristics and immunotherapy strategies of ccRCC. However, several limitations should be recognized. Although we analyzed immune cell characteristics in a scRNA-seq dataset, the tumorinfiltrating immune cells were obtained based on algorithms, and thus, further experimental validation in vitro/in vivo is needed. Our study was also limited by the lack of clinical datasets to verify the relationship between the MRPRS and patients receiving targeted treatment or ICB combined targeted treatment. The combination of an MRPRS-based panel with prospective clinical trials is worth carrying out in the future.

CONCLUSION
This study revealed a significant association between MRPRS and TIICs of ccRCC. The proposed MRPRS is a promising biomarker to predict clinical outcomes and therapeutic responses in ccRCC patients.

AUTHOR CONTRIBUTIONS
AZ and ZZ designed this work. MY, XL, HX, SS, FW, DC, and ZW integrated and analyzed the data. AZ, MY, and XL wrote this manuscript. AZ, ZZ, ZW, and GL edited and revised the manuscript. All authors contributed to the article and approved the submitted version.