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 (m⁶A) modification has been demonstrated in recent studies. Therefore, m⁶A modification-mediated immune cell infiltration characteristics may be helpful to guide immunotherapy for ccRCC.

Methods: This study comprehensively analyzed m⁶A 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 m⁶A 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 m⁶A 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 biomarker-based therapeutic strategies for advanced ccRCC have been missing.

N6-methyladenosine (m⁶A) modification is an important factor for messenger RNA (mRNA) stability, splicing, and translation (5–7). Serval m⁶A-sequencing studies have revealed that abnormal m⁶A regulatory enzymes are involved in mutagenesis, proliferation, and tumorigenesis through the dysregulation of the m⁶A pathway (8, 9). Recently, m⁶A modifications have been shown to play a role in the regulation of immune cells, such as the following: METTL3-mediated m⁶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⁶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⁶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⁶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.

Materials and Methods

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–15). The relative transcriptomic and clinical data of three immunotherapeutic cohorts of patients with ccRCC were obtained from the online supplementary data (16–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–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⁶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).

Selection of m⁶A RNA Methylation Regulators

Based on previous studies (5, 27–30), 23 m⁶A RNA methylation regulators, namely, ALKBH5, CBLL1, FMR1, IGF2BP1/2/3, FTO, YTHDC1/2, YTHDF1/2/3, HNRNPC, LRPPRC, METTL3/14/16, WTAP, KIAA1429, RBM15/15B, ZC3H13, and HNRNPA2B1, were used for our analysis. Immunohistochemistry (IHC) images of m⁶A regulators have been used in the tissue atlas and pathology atlas panels in the Human Protein Atlas. The protein and gene expression of m⁶A regulators in normal individuals and ccRCC patients were analyzed on University of Alabama Cancer Database (UALCAN) portal.

Immune Infiltration Analysis in RCC

CIBERSORT and MCP counter were used to transform the RNA-seq data into the proportion of TIICs. The MCP counter R package was used to evaluate the expression of nine TIICs types. CIBERSORT (https://cibersort.stanford.edu/) was used to quantify the 22 infiltrated immune cells according to normalized gene expression profiles, which included different types of B cells, T cells, NK cells, DC cells, and mast cells. As a verification method, the single-sample gene-set enrichment analysis (ssGSEA) and xCell algorithm were applied.

Construction and Validation of the m⁶A Gene Signature

The significant m⁶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 cross-validation). 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⁶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⁶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.

Results

The Expressional Alteration Landscape of m⁶A Regulators in ccRCC Tissues and Cell Clusters

On reviewing the literature (5, 27–30), 23 genes were found that mainly regulate m⁶A modification including 9 writers (METTL3, METTL14, METTL16, RBM15, RBM15B, CBLL1, ZC3H13, KIAA1429, and WTAP), 12 readers (FMR1, HNRNPA2B1, HNRNPC, YTHDF1/2/3, YTHDC1/2, IGF2BP1/2/3, and LRPPRC), and 2 erasers (FTO and ALKBH5). We utilized the bulk TCGA-ccRCC data (529 cases of ccRCC and 74 cases of normal tissues) to analyze the expression of these m⁶A regulators, revealing that 15 out of 23 m⁶A regulators were differentially expressed (Figure 1A; Table S3), and this phenomenon was also found in two other GEO datasets (GSE53757 and GSE40435) (Figure S1A). In addition, the protein levels of these regulators were also evaluated from the IHC results (Figure S1B; Table S4).

FIGURE 1

Figure 1 (A) Expression of 23 m⁶A RNA methylation regulators between renal cancer and normal tissues in the TCGA-ccRCC cohort. (B, C) The UMAP plot and overview of epithelial cells by the origin and cell type of the cells. (D) Composition of various epithelial cells in different immunotherapeutic responses. (E) The heatmap of marker gene expression in 7 identified epithelial cell subsets. (F) Dot plot analysis of KEGG pathway enrichment of 7 epithelial cell subsets. (G) Violin plots showing the partial expression of m⁶A regulators for each epithelial cell type ("*p < 0.05; **p < 0.01; ***p < 0.001; ****p < 0.0001; ns, not significant).

We next used the scRNA-seq data (25) (PMID33861994) to evaluate the expression of m⁶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⁶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-response-related VCAM1⁺ epithelial cells (Figure 1D), indicating that the differential expression level of m⁶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⁶A Regulator Prognostic Risk Score

To systematically evaluate these differences in m⁶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.

FIGURE 2

Figure 2 (A) LASSO coefficient profiles of the 23 m⁶A RNA methylation regulators in the TCGA-ccRCC cohort. (B) The prognostic analyses for 23 m⁶A RNA methylation regulators in the TCGA-ccRCC cohort using the univariate Cox regression model. (C) Kaplan–Meier analysis of patients between high- and low-MRPRS groups in the TCGA-ccRCC cohort. (D) Validation cohort of MRPRS from GSE29609. (E) Box plot of different MRPRSs in 3 epithelial cell subsets. (F) The MRPRS between response and non-response groups in PMID32472114. (G–I) Kaplan–Meier analysis of three validation cohorts of immunotherapy in ccRCC (PMID32472114, PMID29301960, and PMID32895571). (J) Heatmap of chemotherapy and targeted drug-related genes between high- and low-MRPRS groups. (K) The correlation of MRPRS and genes associated with angiogenesis. (L) The MRPRS between response and non-response groups in ccRCC with sunitinib (E-MTAB-3267). (M) Kaplan–Meier analysis between high- and low-MRPRS groups in ccRCC with sunitinib (E-MTAB-3267) (*p < 0.05; **p < 0.01; ***p < 0.001).

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 therapy- and 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).

FIGURE 3

Figure 3 (A) Box plot of the relationship between stage T, N, M, grade, and MRPRS. (B) The heatmap of markers on multiple immune infiltrates. (C) The MCP_counter algorithm was used to estimate the abundance of various types of immune cells between high- and low-MRPRS groups. (D) The abundance and diversity of TCR clone in high- and low-MRPRS groups. (E) Crosstalk between immune cells and epithelial cells (*p < 0.05; **p < 0.01; ***p < 0.001; ****p < 0.0001).

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 high MRPRS (Figure 3C and Figures S5A–C; 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 immune-related 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⁶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 single-nucleotide 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).

FIGURE 4

Figure 4 (A) Waterfall plot of the distribution of mutations found in the high- and low-MRPRS groups of the TCGA-ccRCC cohort. (B) The top 5 genes of high- vs. low-MRPRS group mutation status. (C) Lollipop plot of somatic mutations in SETD2, TRIOBP, RYR2, ZFPM2, and ABCC6. (D, E) MCP_counter and ssGSEA algorithm were used to estimate the abundance of various types of immune cells in high- and low-MRPRS groups (*p < 0.05; **p < 0.01; ***p < 0.001; ****p < 0.0001; ns, not significant).

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–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).

FIGURE 5

Figure 5 (A) Different distribution ratio of response and non-response in the immunotherapeutic cohort of ccRCC. (B) Kaplan–Meier analysis of patients in the mutated and non-mutated groups in the immunotherapeutic cohort of ccRCC (Van_2018, Morris_2019, PMID29337640 and PMID29301960). (C) The composition of major cancer types in the immunotherapeutic cohort contains mutations of pan-cancer. (D, E) Kaplan–Meier analysis (OS and PFS) of patients in the mutation and non-mutation groups in the immunotherapeutic cohort of pan-cancer. (F) Kaplan–Meier analysis of patients in the low TMB of mutated, low TMB of non-mutated, and high TMB groups.

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⁶A modification plays an indispensable role in immunity, inflammation, and therapy resistance through various m⁶A regulators (36). In this study, we systematically evaluated the expression level of m⁶A regulators in the ccRCC tissue and cell clusters and focused on the detailed relationship between m⁶A modification and TME to enhance our understanding of the ccRCC-immune crosstalk. We constructed an MRPRS comprising 8 m⁶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¹+ 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⁶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⁶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 tumor-infiltrating 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.

Data Availability Statement

All data used in this work can be acquired from the Gene-Expression Omnibus (GEO; https://www.ncbi.nlm.nih.gov/geo/) under the accession numbers GSE53757, GSE40435 and GSE29609, The Cancer Genome Atlas (TCGA) data portal (https://portal.gdc.cancer.gov/), The European Molecular Biology Laboratory (https://www.embl.org/), UALCAN (http://ualcan.path.uab.edu/), HumanProtein Atlas (HPA, https://www.proteinatlas.org/) and https://www.ebi.ac.uk/arrayexpress/experiments/E-MTAB-3267/. The authors would like to thank the above database for the data provided.

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.

Funding

This work was supported by the National Natural Science Foundation of China (Reference number: 81402117) and the Qianjiang Talent Project of Zhejiang Province (Reference number: QJD1602025)

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2022.818120/full#supplementary-material

Abbreviations

ccRCC, clear cell renal cell carcinoma; m⁶A, N6-methyladenosine; MRPRS, m⁶A regulator prognostic risk score; TMB, tumor mutational burden; mRNAs, messenger RNAs; DCs, dendritic cells; TIICs, tumor-infiltrating immune cells; ccRCC, clear cell renal cell carcinoma; TCGA, The Cancer Genome Atlas; GEO, Gene Expression Omnibus; ICB, immune checkpoint blockade; WES, whole-exome sequencing; IHC, immunohistochemistry; scRNA-Seq, single-cell RNA sequencing; PCA, principal component analysis; UMAP, uniform manifold approximation and projection; MsigDB, molecular signature database; TME, tumor microenvironment; ssGSEA, single-sample gene-set enrichment analysis; LASSO, least absolute shrinkage and selection operator; DEGs, differentially expressed genes; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; PPI, protein–protein interaction; MHC, major histocompatibility complex; TCR, T-cell receptor; SNV, single-nucleotide variant; OS, overall survival; PFS, progression-free survival; TLSs, tertiary lymphoid structures; VEGF, vascular endothelial growth factor; mTOR, mammalian target of rapamycin; ESCA, esophageal carcinoma; LUAD, lung adenocarcinoma; PAAD, pancreatic adenocarcinoma; ACC, adrenocortical carcinoma; BLCA, bladder cancer; COAD, colon adenocarcinoma; GBM, glioblastoma multiforme; NSCLC, non-small-cell lung cancer; KIRP, kidney renal papillary cell carcinoma; LAML, acute myeloid leukemia; LGG, brain low-grade glioma; LIHC, liver hepatocellular carcinoma; MESO, mesothelioma; READ, rectum adenocarcinoma; SARC, sarcoma; STAD, stomach adenocarcinoma.

References

1. Siegel RL, Miller KD, Jemal A. Cancer Statistics, 2020. CA Cancer J Clin (2020) 70(1):7–30. doi: 10.3322/caac.21590

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Li QK, Pavlovich CP, Zhang H, Kinsinger CR, Chan DW. Challenges and Opportunities in the Proteomic Characterization of Clear Cell Renal Cell Carcinoma (ccRCC): A Critical Step Towards the Personalized Care of Renal Cancers. Semin Cancer Biol (2019) 55:8–15. doi: 10.1016/j.semcancer.2018.06.004

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Rini BI, Plimack ER, Stus V, Gafanov R, Hawkins R, Nosov D, et al. Pembrolizumab Plus Axitinib Versus Sunitinib for Advanced Renal-Cell Carcinoma. N Engl J Med (2019) 380(12):1116–27. doi: 10.1056/NEJMoa1816714

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Rini BI, Powles T, Atkins MB, Escudier B, McDermott DF, Suarez C, et al. Atezolizumab Plus Bevacizumab Versus Sunitinib in Patients With Previously Untreated Metastatic Renal Cell Carcinoma (IMmotion151): A Multicentre, Open-Label, Phase 3, Randomised Controlled Trial. Lancet (2019) 393(10189):2404–15. doi: 10.1016/S0140-6736(19)30723-8

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Zaccara S, Ries RJ, Jaffrey SR. Reading, Writing and Erasing mRNA Methylation. Nat Rev Mol Cell Biol (2019) 20(10):608–24. doi: 10.1038/s41580-019-0168-5

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Huang H, Weng H, Sun W, Qin X, Shi H, Wu H, et al. Recognition of RNA N(6)-Methyladenosine by IGF2BP Proteins Enhances mRNA Stability and Translation. Nat Cell Biol (2018) 20(3):285–95. doi: 10.1038/s41556-018-0045-z

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Zhao X, Yang Y, Sun BF, Shi Y, Yang X, Xiao W, et al. FTO-Dependent Demethylation of N6-Methyladenosine Regulates mRNA Splicing and is Required for Adipogenesis. Cell Res (2014) 24(12):1403–19. doi: 10.1038/cr.2014.151

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Lin S, Choe J, Du P, Triboulet R, Gregory RI. The M(6)A Methyltransferase METTL3 Promotes Translation in Human Cancer Cells. Mol Cell (2016) 62(3):335–45. doi: 10.1016/j.molcel.2016.03.021

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Liu J, Eckert MA, Harada BT, Liu SM, Lu Z, Yu K, et al. M(6)A mRNA Methylation Regulates AKT Activity to Promote the Proliferation and Tumorigenicity of Endometrial Cancer. Nat Cell Biol (2018) 20(9):1074–83. doi: 10.1038/s41556-018-0174-4

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Wang H, Hu X, Huang M, Liu J, Gu Y, Ma L, et al. Mettl3-Mediated mRNA M(6)A Methylation Promotes Dendritic Cell Activation. Nat Commun (2019) 10(1):1898. doi: 10.1038/s41467-019-09903-6

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Qiu X, Yang S, Wang S, Wu J, Zheng B, Wang K, et al. M6A Demethylase ALKBH5 Regulates PD-L1 Expression and Tumor Immunoenvironment in Intrahepatic Cholangiocarcinoma. Cancer Res (2021) 81(18):4778–93. doi: 10.1158/0008-5472.CAN-21-0468

PubMed Abstract | CrossRef Full Text | Google Scholar

12. von Roemeling CA, Radisky DC, Marlow LA, Cooper SJ, Grebe SK, Anastasiadis PZ, et al. Neuronal Pentraxin 2 Supports Clear Cell Renal Cell Carcinoma by Activating the AMPA-Selective Glutamate Receptor-4. Cancer Res (2014) 74(17):4796–810. doi: 10.1158/0008-5472.CAN-14-0210

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Wozniak MB, Le Calvez-Kelm F, Abedi-Ardekani B, Byrnes G, Durand G, Carreira C, et al. Integrative Genome-Wide Gene Expression Profiling of Clear Cell Renal Cell Carcinoma in Czech Republic and in the United States. PLoS One (2013) 8(3):e57886. doi: 10.1371/journal.pone.0057886

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Edeline J, Mottier S, Vigneau C, Jouan F, Perrin C, Zerrouki S, et al. Description of 2 Angiogenic Phenotypes in Clear Cell Renal Cell Carcinoma. Hum Pathol (2012) 43(11):1982–90. doi: 10.1016/j.humpath.2012.01.023

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Beuselinck B, Job S, Becht E, Karadimou A, Verkarre V, Couchy G, et al. Molecular Subtypes of Clear Cell Renal Cell Carcinoma are Associated With Sunitinib Response in the Metastatic Setting. Clin Cancer Res (2015) 21(6):1329–39. doi: 10.1158/1078-0432.CCR-14-1128

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Braun DA, Hou Y, Bakouny Z, Ficial M, Sant' Angelo M, Forman J, et al. Interplay of Somatic Alterations and Immune Infiltration Modulates Response to PD-1 Blockade in Advanced Clear Cell Renal Cell Carcinoma. Nat Med (2020) 26(6):909–18. doi: 10.1038/s41591-020-0839-y

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Miao D, Margolis CA, Gao W, Voss MH, Li W, Martini DJ, et al. Genomic Correlates of Response to Immune Checkpoint Therapies in Clear Cell Renal Cell Carcinoma. Science (2018) 359(6377):801–6. doi: 10.1126/science.aan5951

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Motzer RJ, Robbins PB, Powles T, Albiges L, Haanen JB, Larkin J, et al. Avelumab Plus Axitinib Versus Sunitinib in Advanced Renal Cell Carcinoma: Biomarker Analysis of the Phase 3 JAVELIN Renal 101 Trial. Nat Med (2020) 26(11):1733–41. doi: 10.1038/s41591-020-1044-8

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Mariathasan S, Turley SJ, Nickles D, Castiglioni A, Yuen K, Wang Y, et al. TGFbeta Attenuates Tumour Response to PD-L1 Blockade by Contributing to Exclusion of T Cells. Nature (2018) 554(7693):544–8. doi: 10.1038/nature25501

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Gide TN, Quek C, Menzies AM, Tasker AT, Shang P, Holst J, et al. Distinct Immune Cell Populations Define Response to Anti-PD-1 Monotherapy and Anti-PD-1/Anti-CTLA-4 Combined Therapy. Cancer Cell (2019) 35(2):238–55.e6. doi: 10.1016/j.ccell.2019.01.003

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Van Allen EM, Miao D, Schilling B, Shukla SA, Blank C, Zimmer L, et al. Genomic Correlates of Response to CTLA-4 Blockade in Metastatic Melanoma. Science (2015) 350(6257):207–11. doi: 10.1126/science.aad0095

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Miao D, Margolis CA, Vokes NI, Liu D, Taylor-Weiner A, Wankowicz SM, et al. Genomic Correlates of Response to Immune Checkpoint Blockade in Microsatellite-Stable Solid Tumors. Nat Genet (2018) 50(9):1271–81. doi: 10.1038/s41588-018-0200-2

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Samstein RM, Lee CH, Shoushtari AN, Hellmann MD, Shen R, Janjigian YY, et al. Tumor Mutational Load Predicts Survival After Immunotherapy Across Multiple Cancer Types. Nat Genet (2019) 51(2):202–6. doi: 10.1038/s41588-018-0312-8

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Rizvi H, Sanchez-Vega F, La K, Chatila W, Jonsson P, Halpenny D, et al. Molecular Determinants of Response to Anti-Programmed Cell Death (PD)-1 and Anti-Programmed Death-Ligand 1 (PD-L1) Blockade in Patients With Non-Small-Cell Lung Cancer Profiled With Targeted Next-Generation Sequencing. J Clin Oncol (2018) 36(7):633–41. doi: 10.1200/JCO.2017.75.3384

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Krishna C, DiNatale RG, Kuo F, Srivastava RM, Vuong L, Chowell D, et al. Single-Cell Sequencing Links Multiregional Immune Landscapes and Tissue-Resident T Cells in ccRCC to Tumor Topology and Therapy Efficacy. Cancer Cell (2021) 39(5):662–77.e6. doi: 10.1016/j.ccell.2021.03.007

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Efremova M, Vento-Tormo M, Teichmann SA, Vento-Tormo R. CellPhoneDB: Inferring Cell-Cell Communication From Combined Expression of Multi-Subunit Ligand-Receptor Complexes. Nat Protoc (2020) 15(4):1484–506. doi: 10.1038/s41596-020-0292-x

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Li H, Hu J, Yu A, Othmane B, Guo T, Liu J, et al. RNA Modification of N6-Methyladenosine Predicts Immune Phenotypes and Therapeutic Opportunities in Kidney Renal Clear Cell Carcinoma. Front Oncol (2021) 11:642159. doi: 10.3389/fonc.2021.642159

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Chen L, Hu B, Song X, Wang L, Ju M, Li Z, et al. M(6)A RNA Methylation Regulators Impact Prognosis and Tumor Microenvironment in Renal Papillary Cell Carcinoma. Front Oncol (2021) 11:598017. doi: 10.3389/fonc.2021.598017

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Huang H, Weng H, Chen J. M(6)A Modification in Coding and Non-Coding RNAs: Roles and Therapeutic Implications in Cancer. Cancer Cell (2020) 37(3):270–88. doi: 10.1016/j.ccell.2020.02.004

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Chen XY, Zhang J, Zhu JS. The Role of M(6)A RNA Methylation in Human Cancer. Mol Cancer (2019) 18(1):103. doi: 10.1186/s12943-019-1033-z

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Szklarczyk D, Franceschini A, Wyder S, Forslund K, Heller D, Huerta-Cepas J, et al. STRING V10: Protein-Protein Interaction Networks, Integrated Over the Tree of Life. Nucleic Acids Res (2015) 43:D447–52. doi: 10.1093/nar/gku1003

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Diaz-Montero CM, Rini BI, Finke JH. The Immunology of Renal Cell Carcinoma. Nat Rev Nephrol (2020) 16(12):721–35. doi: 10.1038/s41581-020-0316-3

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Thommen DS, Schumacher TN. T Cell Dysfunction in Cancer. Cancer Cell (2018) 33(4):547–62. doi: 10.1016/j.ccell.2018.03.012

PubMed Abstract | CrossRef Full Text | Google Scholar

34. He QF, Xu Y, Li J, Huang ZM, Li XH, Wang X. CD8+ T-Cell Exhaustion in Cancer: Mechanisms and New Area for Cancer Immunotherapy. Brief Funct Genomics (2019) 18(2):99–106. doi: 10.1093/bfgp/ely006

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Motzer RJ, Escudier B, McDermott DF, George S, Hammers HJ, Srinivas S, et al. Nivolumab Versus Everolimus in Advanced Renal-Cell Carcinoma. N Engl J Med (2015) 373(19):1803–13. doi: 10.1056/NEJMoa1510665

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Shulman Z, Stern-Ginossar N. The RNA Modification N(6)-Methyladenosine as a Novel Regulator of the Immune System. Nat Immunol (2020) 21(5):501–12. doi: 10.1038/s41590-020-0650-4

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Helmink BA, Reddy SM, Gao J, Zhang S, Basar R, Thakur R, et al. B Cells and Tertiary Lymphoid Structures Promote Immunotherapy Response. Nature (2020) 577(7791):549–55. doi: 10.1038/s41586-019-1922-8

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Cabrita R, Lauss M, Sanna A, Donia M, Skaarup Larsen M, Mitra S, et al. Tertiary Lymphoid Structures Improve Immunotherapy and Survival in Melanoma. Nature (2020) 577(7791):561–5. doi: 10.1038/s41586-019-1914-8

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Sanmamed MF, Chen L. A Paradigm Shift in Cancer Immunotherapy: From Enhancement to Normalization. Cell (2018) 175(2):313–26. doi: 10.1016/j.cell.2018.09.035

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Ji RR, Chasalow SD, Wang L, Hamid O, Schmidt H, Cogswell J, et al. An Immune-Active Tumor Microenvironment Favors Clinical Response to Ipilimumab. Cancer Immunol Immunother (2012) 61(7):1019–31. doi: 10.1007/s00262-011-1172-6

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Gajewski TF, Corrales L, Williams J, Horton B, Sivan A, Spranger S. Cancer Immunotherapy Targets Based on Understanding the T Cell-Inflamed Versus Non-T Cell-Inflamed Tumor Microenvironment. Adv Exp Med Biol (2017) 1036:19–31. doi: 10.1007/978-3-319-67577-0_2

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Conway JR, Kofman E, Mo SS, Elmarakeby H, Van Allen E. Genomics of Response to Immune Checkpoint Therapies for Cancer: Implications for Precision Medicine. Genome Med (2018) 10(1):93. doi: 10.1186/s13073-018-0605-7

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Hellmann MD, Ciuleanu TE, Pluzanski A, Lee JS, Otterson GA, Audigier-Valette C, et al. Nivolumab Plus Ipilimumab in Lung Cancer With a High Tumor Mutational Burden. N Engl J Med (2018) 378(22):2093–104. doi: 10.1056/NEJMoa1801946

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Snyder A, Makarov V, Merghoub T, Yuan J, Zaretsky JM, Desrichard A, et al. Genetic Basis for Clinical Response to CTLA-4 Blockade in Melanoma. N Engl J Med (2014) 371(23):2189–99. doi: 10.1056/NEJMoa1406498

PubMed Abstract | CrossRef Full Text | Google Scholar

45. McGrail DJ, Pilie PG, Rashid NU, Voorwerk L, Slagter M, Kok M, et al. High Tumor Mutation Burden Fails to Predict Immune Checkpoint Blockade Response Across All Cancer Types. Ann Oncol (2021) 32(5):661–72. doi: 10.1016/j.annonc.2021.02.006

PubMed Abstract | CrossRef Full Text | Google Scholar