The effect of m6A RNA methylation regulators on tumor microenvironment in clear cell renal cell carcinoma

Background: m6A RNA methylation and tumor microenvironment (TME) have been reported to play important roles in the progression and prognosis of clear cell renal cell carcinoma (ccRCC). However, whether m6A RNA methylation regulators affect TME in ccRCC remains unknown. Thus, the current study is designed to comprehensively evaluate the effect of m6A RNA methylation regulators on TME in ccRCC. Methods: Transcriptome data of ccRCC was obtained from The Cancer Genome Atlas (TCGA) database. Consensus clustering analysis was conducted based on the expressions of m6A RNA methylation regulators. Survival differences were evaluated by Kaplan-Meier (K-M) analysis between the clusters. DESeq2 package was used to analyze the differentially expressed genes (DEGs) between the clusters. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis were analyzed by ClusterProler R package. Immune, stromal and ESTIMATE scores were assessed by ESTIMATE algorithm. CIBERSORT algorithm was applied to evaluate immune inltration. The expressions of human leukocyte antigen (HLA), immune checkpoint molecules, and Th1/IFNγ gene signature associated with TME were also compared between the clusters. TIDE algorithm and subclass mapping were used to analyze the clinical response of different clusters to PD-1 and CTLA-4 blockade. Results: The expressions of fteen m6A regulators were signicantly different between ccRCC and normal kidney tissues. Based on the expressions of those fteen m6A regulators, two clusters were identied by consensus clustering, in which cluster 1 had better overall survival (OS). A total of 4,429 DEGs were found between the two clusters, and were enriched into immune-related biological processes. Further analysis of the two clusters’ TME showed that cluster 1 had lower immune and ESTIMATE scores, higher expressions of HLA and lower expressions of immune checkpoint molecules. Besides, immune inltration and the expressions of Th1/IFNγ gene signature also have signicant differences between two clusters. Conclusions: Our study that m6A regulators were important participants in the development of ccRCC, with a close relationship with TME.


Background
Renal cell carcinoma (RCC) is ranked sixth in male and tenth in female among the most common cancers [1]. Clear cell renal cell carcinoma (ccRCC), which accounts for 70-80% of all RCC, is the most frequently primary RCC [2]. In recent years, the treatment of RCC has made great progress with the development of surgical resection, chemotherapy and radiotherapy. However, the overall survival (OS) are still dismal [3].
It is reported that one third of the patients relapsed after an average of 1.9 years [4], and the patients who were local or distant metastasis further led to high mortality [5]. Unfortunately, effective RCC prognostic markers were still lacked due to the high heterogeneity and complex disease process. Therefore, it is of great signi cance to analyze the pathogenesis from molecular level, and to predict the prognosis of ccRCC patients accurately.
The signi cance of reversible RNA modi cations has always been underestimated. N6-methyladenosine (m6A), which was rstly found in 1970s and has gained renewed interests as a new layer of control for gene expression, participates in more than 60% of RNA methylation and has become the most common post transcriptional modi cation on mRNA and [6][7][8]. Studies have shown that the regulation of m6A was in a reversible dynamic equilibrium [9]. M6A could be produced by m6A methyltransferases ("writers") and removed by m6A demethylases ("erasers"). Additionally, m6A-binding proteins ("readers") including YT521-B homology (YTH) domain family, insulin-like growth factor 2 mRNA-binding proteins (IGF2BPs), and heterogeneous nuclear ribonucleoproteins (HNRNPs), were responsible for identi cation of m6Amodi ed targeted RNA [10]. Studies have been demonstrated that m6A RNA methylation regulators were involved in tumorigenesis and impacted the patient's survival in various cancer, including ccRCC [11][12][13][14].
It is well established that cancer progression is inevitably related to abnormal tumor microenvironment (TME). TME is the primary battle eld where host immune system and tumor cells interact. Deterioration of TME promotes cancer cells escape from immune surveillance, and accelerates the occurrence, development and metastasis of RCC [15]. In addition, the interactions of chemokines and its receptors can recruit immune cell subsets into TME and so as to further modulation of tumor progression and metastasis [16]. In order to characterize the TME of ccRCC, Yasin et al employed a computational method based on gene expression to pro le the in ltration levels of 24 immune cell populations in 19 cancer types, and found the following results: (1) the ccRCC had got the highest immune in ltration scores when compared against 18 other tumor types, (2) the immunogenicity of ccRCC tumors was highly correlated with MHC class I antigen presenting machinery expression, but not neo-antigen load or mutation load, (3) Th17 cells and CD8+ T/Treg ratio were associated with improved survival, whereas Th2 cells and Tregs were associated with negative outcomes [17].
Abnormal m6A regulators were associated with TME and immunotherapy responses in melanoma [18], colorectal [19], gastric cancer [20] and Head and Neck Squamous Cell Carcinoma [21] , however, the interaction of m6A regulators and microenvironment in ccRCC remains unclear. Thus, the purpose of this study is to comprehensively evaluate the correlation of m6A RNA methylation regulators with prognosis, immune checkpoint molecules, and TME in ccRCC, so as to better understand the etiology and provide a new insight for immunotherapy of ccRCC.
Consensus clustering analysis "ConsensusClusterPlus" R package was applied in consensus clustering analysis to categorize ccRCC patients into clusters. The difference of OS between different clusters was assessed by the Kaplan-Meier method and log-rank test.

GO and KEGG pathway enrichment analysis of DEGs
DEGs between different clusters were identi ed using DEseq2 package of R. P-value < 0.05 and |log 2 (Fold change)| > 1 were used to judge the statistical signi cance of gene expression differences.
"ClusterPro ler R package" was used to analyze the GO and KEGG pathway enrichment. The GO terms include biological process, cellular component and molecular function. Adjusted P < 0.05 was considered statistically signi cant.

Evaluation of tumor microenvironment (TME)
Immune score and ESTIMATE score were calculated by ESTIMATE algorithm using "estimate" R package. Estimation of immunotherapeutic response TIDE algorithm and subclass mapping were used to analyze the difference of response to PD-1 and CTLA-4 blockade between different clusters.

Statistical analysis
All data were analyzed by R (version 4.0.0), and Mann-Whitney-Wilcoxon test was used to compare the data from different clusters, and P < 0.05 was considered statistically signi cant.

Results
The expressions of m6A regulators were related to the development of ccRCC According to the previous studies, we selected 21 m6A regulators for further investigation. As shown in Figure 1A-B, the expressions of 15 m6A regulators were signi cantly different between ccRCC and normal kidney tissues, among which the expressions of METTL3, WTAP, RBM15, FTO, ALKBH5, YTHDC2 and IGF2BP3 were higher in ccRCC samples than that of normal samples. Compared with ccRCC samples, the expressions of METTL14, ZC3H13, RBM15B, YTHDF2, YTHDF3, IGF2BP2, HNRNPC and HNRNPA2B1 signi cantly increased in normal kidney samples. The result of Kaplan-Meier analysis showed that nine m6A regulators including METTL3, WTAP, YTHDC2, IGF2BP3, METTL14, ZC3H13, YTHDF3, IGF2BP2, and HNRNPA2B1 were associated with the prognosis of ccRCC patients. ccRCC patients with higher expressions of METTL3, IGF2BP3, WTAP, IGF2BP2 and HNRNPA2B1 had a more favorable prognosis, while patients with higher expressions of METTL14, ZC3H13, YTHDF3 and YTHDC2 had a less survival time ( Figure 1C). These results indicated that the expressions of m6A regulators were related to the development of ccRCC.
Identi cation of two subtypes in ccRCC based on the expressions of m6A regulators Consensus clustering was then applied to classify ccRCC patients into clusters based on the expressions of these 15 m6A regulators. Clearly, K=2 was identi ed with optimal clustering (Figure 2A-B). A total of 479 ccRCC patients were clustered into two subgroups with 337 patients in cluster 1 and 142 patients in cluster 2. Besides, the results of PCA clearly showed two distinct distribution of the two clusters ( Figure   2C). Kaplan-Meier analysis showed that the patients of cluster 1 had a better OS compared to those of cluster 2 ( Figure 2D).

Biological function analysis of DEGs between the two clusters
To better characterize the difference of biological process involved in the development of ccRCC between the two clusters, we compared the expressions of genes in the two clusters and found a total of 4,429 DEGs ( Figure 3A). To get further insights into the biological characteristics of the DEGs, we performed GO and KEGG pathway analysis. GO analysis showed that the DEGs were enriched in immune-related functions, such as regulation of immune effector process, production of molecular mediator of immune response, immune response-activating cell surface receptor signaling pathway, lymphocyte mediated immunity, B cell mediated immunity, immunoglobulin mediated immune response and regulation of humoral immune response ( Figure 3B). Moreover, the result of KEGG pathway analysis revealed several immune-related pathways, such as cytokine-cytokine receptor interaction, calcium signaling pathway, complement and coagulation cascades and bile secretion ( Figure 3C). These results demonstrated that m6A regulators were involved in the immunity of ccRCC patients.
Identi cation of different TME between the two clusters It has been reported that TME is closely associated with immune response in ccRCC progression [17]. Therefore, we studied the TME characteristics of the two clusters to further explore the role of m6A regulators in ccRCC. The expression pattern of m6A regulators and TME features in cluster 1 and cluster 2 were shown by the heatmap ( Figure 4A). Speci cally, the immune and ESTIMATE scores were remarkably higher in cluster 2 (P = 0.0027 and P = 0.018, respectively), and the stromal score of cluster 2 was relatively higher with no statistical signi cance compared to cluster 1 (P = 0.8) ( Figure 4B-D), suggesting that m6A regulators were correlated with immune in ltration and tumor purity.
Identi cation of different immune landscape between the two clusters Given the signi cant difference of TME between the two clusters, we performed the CIBERSORT to estimate the immune in ltration of the two clusters. The distribution of immune cells were compared between cluster 1 and cluster 2 ( Figure 5A-C). The results showed that the abundance of 11 immune cell types were signi cantly different between two clusters, with higher proportions of naive B cells, CD4 + naive T cells, regulatory T cells, activated CD4 + memory T cells and M0 macrophages and lower plasma cells, M1 macrophages, resting NK cells, monocytes, resting rest cells and neutrophils in cluster 2. These results indicated that clustering subgroups based on the expressions of m6A regulators were closely related to the immune microenvironment in ccRCC.
Identi cation of different expressions of HLA, immune checkpoint molecules and Th1/IFNγ pathway signature between the two clusters HLA, immune checkpoints and Th1/IFNγ pathway related genes can interact with TME and play vital roles in the immunoregulation of multiple tumors. Therefore, we further investigated whether m6A regulators in uence on the expressions of HLA and immune checkpoint molecules. Interestingly, we found 12 HLA gene expressions signi cantly differed between two clusters, among which expressions of 11 HLA genes (P < 0.05), including HLA-W, HLA-T, HLA-P, HLA-K, HLA-J, HLA-H, HLA-G, HLA-F, HLA-E, HLA-DQB1, and HLA-DPA2, were signi cantly up-regulated in cluster 1 ( Figure 6A). The 19 immune checkpoints molecules were found to be remarkably differentially expressed between two clusters, of which cluster 2 had higher expressions of 17 immune checkpoints (P < 0.05), including BTLA, CD27, CD276, CD28, CD40LG, CD47, CD80, CD86, ICOS, LAG3, LGALS9, PDCD1, PVR, TIGIT, TNFRSF18, TNFSF14 and TNFSF4 ( Figure 6B). As for the genes involved in Th1/IFNγ pathway, the expressions of IFNGR2, IFNγ and STAT1 were up-regulated in cluster 2 (P < 0.05), whereas the expressions of JAK1 and JAK2 were upregulated in cluster 1 (P < 0.05) ( Figure 6C). These results indicated that m6A regulators may affect the immune response and TME by regulating the expressions of HLA, immune checkpoint molecules and Th1/IFNγ pathway gene signatures.
It has been reported that both TME and immune checkpoints are associated with therapeutic response [22]. Therefore, in consideration of different TME and expressions of immune checkpoints in the two clusters based on the expressions of m6A regulators, we investigated the response of two clusters to PD-1 and CTLA-4 blockade, which has been used in the clinical treatment of ccRCC patients [23]. As shown in Figure 6D, cluster 2 presented a very poor response to anti-CTLA-4 therapy (Bonferroni corrected P = 0.048). These results indicated that m6A regulators have a close relationship with therapeutic e cacy.

Discussion
Increasing evidences demonstrated that m6A modi cation acts as an indispensable role in in ammation, innate immunity as well as anti-tumor effect through interaction with different m6A regulators. In the eld of tumor research, it has been found that m6A regulator-mediated methylation was involved in tumorigenesis and angiogenesis, of which the process was partially through altering TME. However, the effect of m6A RNA methylation regulators on the TME of ccRCC remained unclear. The present study is to evaluate the relationship between m6A RNA methylation regulators, TME and TME-related issues, such as the expressions of immune checkpoint molecules in ccRCC.
Firstly, we compared the expressions of 21 m6A regulators between ccRCC and normal kidney tissues, and found 15 m6A regulators had remarkably different expressions, the expression abundance of which between ccRCC and normal kidney tissues was consistent with the previous reports of ccRCC [11,24]. Among these dysregulated m6A regulators, we noticed that patients with higher expressions of METTL3, IGF2BP3, WTAP, IGF2BP2 and HNRNPA2B1 had a better prognosis, whereas patients with higher expressions of METTL14, ZC3H13, YHDF3 and YTHDC2 had worse prognosis ( Figure.1C). Lobo et al indicated that "writers" were regarded as the main m6A regulators and related with patients' OS in ccRCC [25]. Our present study also con rmed the vital roles of "writers" of METTL3, METTL14, WTAP and ZC3H13 in predication of ccRCC patients' OS. As members of the "reader" regulators, we observed that YTHDC2 and YTHDF3 were negatively with ccRCC patients' OS, while IGF2BP2, IGF2BP3 and HNRNPA2B1 were positive, and these ndings were similarly to previous reports [26 -29]. Additionally, in other studies, erasers of FTO (Niu et al in breast cancer [30]) and ALKBH5(Chen et al in hepatocellular carcinoma [31]) and readers of YTHDF2(Li et al in Osteosarcoma [29], Zhang et al in Lung Adenocarcinoma [32]) and HNRNPC(Zhang et al in gastric cancer [33], Zhao et al in head and neck squamous cell carcinoma [34]) have been reported to be related to OS of patients, however, in the current study, we found these regulators exhibited no signi cant relationships with OS in ccRCC. Therefore, we speculated that the contributions of these m6A regulators to OS differed in different cancers, which may be due to different tumorigenesis mechanisms. More solid experiments were needed to clarify the phenomenon.
Based on the expressions of these 15 m6A methylation regulators, consensus clustering was applied to classify ccRCC patients into two subgroups, cluster 1 and cluster 2 ( Figure.2A). We found cluster 1 had better OS than cluster 2 ( Figure.2D). Furthermore, when comparing the expression pro les of the two clusters, we found a total of 4429 DEGs with different expressions. Both GO and KEGG pathway enrichment analysis suggested that the DEGs participate immune related biological process and signalings ( Figure.3B and C), which indicated that m6A regulators might be associated with immunity. Therefore, we next compared the TME and its related indicators (immune cell in ltration, HLA, immune checkpoint and Th1/IFNγ pathway gene expression), which have been reported to be closely related to immunity [35,36], between the two clusters. Interestingly, we found that TME and its related indicators did differ between the two clusters.
In the current study, we found cluster 2 with shorter OS had higher scores of immune, stromal and ESTIMATE. Our results were consistent with the study of Xu et al [37], who divided their gliomas samples into two subgroups and found the cluster with higher stroma (P < 0.05), immune (P < 0.05), and ESTIMATE scores (P < 0.05) were accomplished with shorter OS. However, Liu et al found that patients with higher immune scores were signi cantly positively associated with OS time in endometrial carcinoma [38]. The con icted results may be caused by the diverse immune microenvironment and tumorigenesis mechanism in different tumors.
HLA molecules are mandatory for the immune recognition and subsequent killing of neoplastic cells by the immune system [39]. Our study showed that 11 of 12 HLA genes were high-expressed in cluster 1. It is well known that tumor antigens must be presented in a way restricted by HLA in order to be recognized by T-cell receptors. However, impaired HLA-I or HLA-II expression could inhibit the activation of cytotoxic immune mechanisms or the antigen-presenting capability of antigen presenting cells [39]. Combined with our research, we speculated low-HLA levels in ccRCC contributed to the inactivation of cytotoxic immune mechanisms, which in turn performed a worse OS in cluster 2 than that of cluster 1. Moreover, the expressions of HLA-I, including HLA-A, B, C, E, F, G, H, J, K and L, were reported to in uence the pattern of the immune cell in ltration and stromal cell reaction in the TME [40]. High-grade immune in ltration usually accomplished with low levels of HLA-I positive cancer cells. In consideration of our results, it is inferred that lower HLA-I levels in cluster 2 may contribute to higher immune scores (immune cell in ltration level of tumor tissue) than that in cluster 1.
Then we compared the expressions of immune checkpoint molecules, the activation of which enabled tumor cells to evade immune suppression, and found cluster 2 had relatively higher expressions of immune checkpoints. Liu et al found that the crucial immune checkpoints were mostly positive correlation with "CD8 + T cells", "activated memory CD4 + T cells" and "M1 Macrophages", while having negative correlation with "resting memory CD4 + T cells", "M0 Macrophages" and "activated Dendritic cells" [38]. Similarly, we found higher expressed immune checkpoints of cluster 2 are also accompanied by abundant activated CD4 + memory T cells, M0 Macrophages and lower M1 macrophages. Additionally, due to the differences in immune checkpoint molecules, we inferred that the immunotherapeutic response of cluster 1 and cluster 2 may be different. Therefore, we compared the responses of the two clusters to two commonly used immunotherapy, PD-1 and CTLA-4 treatments. Our results showed that both clusters had poor response to these two treatment and cluster 2 had no signi cant response to CTLA-4 therapy. These results suggested that m6A regulators may in uence the effect of immunotherapy, at least PD-1 and CTLA-4 treatments in ccRCC.
Finally, we detected the expressions of genes involved in Th1/IFN-γ signaling in the two clusters, because Th1/IFN-γ signaling also plays an important role in tumor immune surveillance. It was well-known that JAK1 and JAK2 were activated by IFN-γ when it binded to the cell surface receptors (IFNGR1 and IFNGR2), resulting in STAT dimerization. The complex was then translocated to the nucleus and combined with gamma-activated sequence to activate the interferon regulatory factor 1, which will further blocked T cell-mediated immune response [41,42]. However, in the current study, we found cluster 2 expressed higher levels of IFNGR2, IFN-γ and STAT1, but lower levels of JAK1 and JAK2. Therefore, we speculated that m6A regulators may be also involved in Th1/IFN-γ signaling and disturbed the activation of JAK1 and JAK2, but their roles needed further in-depth research.
Conclusions m6A regulators are important participants in the development of ccRCC, with a close relationship with TME. This study provides a new understanding of the role of m6A methylation in ccRCC, and suggests that m6A regulators may affect the immune response and TME by regulating the expressions of HLA, immune checkpoint molecules and Th1/IFNγ pathway gene signatures.