The analysis of N6-methyladenosine regulators impacting the immune infiltration in clear cell renal cell carcinoma

The tumor immune microenvironment (TIME) and N6-methyladenosine (m6A) are related to the progression of several types of cancer. Nevertheless, the impact of m6A on the TIME of clear cell renal cell carcinoma (ccRCC) remains unclear. This study used an unsupervised clustering algorithm to divide the samples into distinct subgroups. The single sample gene set enrichment analysis (ssGSEA) algorithm to estimate the TIME. The correlation between m6A regulators and immune cells in different subgroups was calculated using Spearman analysis. At last, the relationship between IGF2BP2 and HMGA2 was validated in several datasets, including TCGA-KIRC, GEO, and HPA datasets. We found that m6A regulators were differently expressed in several clinical groups. Based on the expression of m6A regulators, we divided the samples into three subgroups. Then, the survival analysis for these three subgroups showed that the cluster 2 subgroup had poor overall survival (OS). Further, we found that IGF2BP2 and IGF2BP3 were essential components in the cluster 2 subgroup using the principal component analysis (PCA) algorithm. In addition, the expression of these two genes was significantly correlated with survival time. At last, we found that HMGA2 was significantly correlated with IGF2BP2 in several datasets, which indicated that HMGA2 is an essential role in affecting IGF2BP2 regulating the TIME. There is a close correlation between m6A regulators and TIME. Moreover, IGF2BP2 is related to the progression of ccRCC and plays an essential role in affecting the TIME.


Background
Kidney cancer is one of the most commonly diagnosed cancers in the globe [1]. According to the World Health Organization statistics, there are more than 140,000 renal cell carcinoma (RCC)-related deaths [2]. Moreover, clear cell RCC (ccRCC) is the most typical subtype of kidney cancer and contributes to most kidney cancer-related deaths [3,4]. Thus, it is necessary to have a further understanding of ccRCC.
Numerous studies revealed that cancer development needs cancer cells and requires other components, like tumor immune cells and stromal cells [5]. In recent studies, researchers have shown that the tumor immune microenvironment (TIME) strongly impacts the progression and prognosis of tumors. For instance, some studies showed that tumor-infiltrating immune cells are correlated with the prognosis of different tumors [6,7]. Hereby, immunotherapy is a new potential cancer therapy strategy that stimulates and improves the immune system's natural ability to fight cancer cells [8]. Moreover, immunotherapy based on immune checkpoint blockade (ICB, PD1/L1, and CTLA4) has already shown remarkable clinical benefit in a small number of patients with advanced tumors [9,10]. Unfortunately, we have to face the fact that the satisfactory drug efficacy was limited to a small number of patients, and the overall survival (OS) in advanced ccRCC is less than 30% [11,12]. Meanwhile, numerous patients remain in no response or resistance to drugs based on immune checkpoint blockade [13]. Therefore, an in-depth investigation of individual differences in the TIME will contribute to precision immunotherapy.
The post-transcriptional regulation of the mammalian RNA contains large numbers of modified nucleosides, such as 5-methylcytidine (m5C), N6-methyladenosine (m6A), N7-methylguanosine (m7G), a portion of the 5′ terminal cap. m6A has been regarded as the most frequent internal modification in mRNA [14,15]. m6A modification is a kind of dynamic reversible process in mammalian cells, which is regulated by methyltransferases, demethylases, and binding proteins, also known as "writers," "erasers," and "readers" [16]. Recent studies have shown that m6A plays a critical role in regulating RNA transcription, processing, translation, and metabolism [17]. Furthermore, emerging evidence suggests that m6A regulators are also closely associated with oncogenic or tumor-suppressive functions, including proliferation [18], tumorigenesis [19], metastasis [20], and other phenotypes of tumors. Meanwhile, researchers demonstrate that there is a close correlation between m6A modification and TIME [21,22]. Han et al. discovered that the m6A molecular in tumor cells might impact dendritic cells' ability to cross-present tumor neoantigens to T lymphocytes [22]. Yi et al. analyzed the m6A regulators of 1938 gastric cancer samples and revealed the m6A regulations played an essential role in TIME complexity of gastric cancer [23]. Chong et al. analyzed 1307 colon carcinoma samples using multiple bioinformatics methods, finding that m6A modification correlated with the tumor microenvironment of colon carcinoma [24]. However, until now, it has not been comprehensively clarified the prognostic effect of m6A modification in ccRCC and the potential roles of m6A modification in host anti-tumor immune response remain primarily unexplored [17].
This work utilized bioinformatics and high-throughput sequencing data to analyze the correlation between m6A regulators and immune cells in TIME. Further, we analyzed the correlation between m6A regulators and immune modulators, which indicated that m6A was associated with the TIME. In addition, we analyzed the correlation between gene expression and immune cells in TIME and found that IGF2BP2 and IGF2BP3 play an essential role in immune cells in TIME. Meanwhile, HMGA2 was significantly related to IGF2BP2 in the poor prognosis subgroup. These genes could be potential immunotherapeutic targets in ccRCC therapy.

Extracting the data
The Cancer Genome Atlas: the gene expression and the clinical data of TCGA-KIRC (count data) were obtained from the UCSC Xena (https:// xenab rowser. net/ datap ages/), which consists of 534 tumor samples and 73 adjacent normal tissue samples. The mutation of genes and copy number variation (CNV) were obtained from The Cancer Genome Atlas (TCGA) using the GDC-client tool.

Bioinformatics analysis
The immune cells infiltration of samples was calculated via single sample gene set enrichment analysis (ssGSEA) using the R package "GSVA" [29]. Then the relative abundance of each immune cell was normalized to − 1 to 1 via the function "scale" of R. To have further understanding of the biological function of genes related to m6A regulators, weighted correlation network analysis (WGCNA) was conducted to obtain the core module associated with m6A subgroups using R package "WGANA" [30]. Gene Ontology (GO) analysis was undertaken to investigate the possible function of the m6A-associated genes using the R package "clusterProfiler" [31]. Further, gene set enrichment analysis (GSEA) with hall marker gene set was performed to understand the biological process involved in different m6A cluster subgroups using R package "cluster-Profiler." R package "pheatmap" was conducted to visualize the expression of m6A genes in different groups [32]. The STRING database is a powerful function for analyzing protein interaction. We utilized the "STRING" app to calculate and visualize the interaction network of m6A regulators [33]. The gene expression relationship was calculated by R package "corrplot" using Spearman correlation analysis [34]. To estimate the tumor microenvironment, we utilized R package "estimate" to assess the stromal cells and immune cells of each sample [35], by which the immune scores and stromal scores were obtained. In this study, we used R package "survminer" to determine the best cutoff [36] and used R package "survival" to conduct K-M survival analysis [37]. The TMB of each sample was calculated and visualized via R package "maftool" [38].

Machine learning analysis
Unsupervised clustering analysis was conducted to detect distinct m6A modification patterns based on m6A regulators gene expression. In this process, we utilized R package "ConsensusClusterPlus" to classify the samples based on the k-means method [39]. Principal component analysis (PCA) is an algorithm to reduce the dimensionality of multivariate data to two or three that can be visualized graphically with minimal loss of information. Furthermore, PCA also can detect the essential components in multivariate data. We utilized R package "factoextra" [40] and "FactoMineR" [41] to calculate the principle factors in m6A regulators-related subgroups.

Statistical analysis
This study used R software (version 4.0.3; The R Foundation for Statistical Computing) for data analysis and statistics. The expression of m6A genes in different groups was compared via Wilcox test. K-M survival analysis and log-rank tests analyzed the OS, DSS, and PFI in different subgroups. Hazard ratio (HR) was generated via univariate Cox proportional hazards regression when examining the role of immune cells and m6A genes in OS. The correlation between m6A genes and tumor-infiltrating immune cells was evaluated by Spearman correlation. All statistical tests were two-sided. P < 0.05 was regarded as statically.

The expression of m6A regulators in various clinical characteristics
The workflow of this study is shown in Fig. S1. In this work, we enrolled 23 m6A RNA methylation regulatory genes in ccRCC. Figure 1A shows the modification process of three types of m6A regulators, including methylation, demethylation, and recognition of m6A. Among the 503 patients with ccRCC, the expression of the 23 m6A regulators was discordant between tumor tissues and adjacent normal tissues (Fig. 1B, C). The expression of "readers," "writers," and "erasers" were all significantly differently expressed in different clinical characters. FMR1, HNRNPA2B1, IGF2BP2, LRPPRC, YTHDF2, YTHDF3, METTL14, RBM15B, and ZC3H13 were downregulated considerably in tumors tissues, while IGF2BP3, YTHDC2, ALKBH5, KIAA1429, METTL13, RBM15B, and WTAP were significantly upregulated in tumor tissues (Fig. 1C). Meanwhile, the m6A regulators were differently expressed across clinical grade and pathological stage (Fig. 1D, E). The correlation analysis showed that most m6A regulators were related to each other, and FTO was significantly correlated with the most "writers" of m6A regulators (Fig. 1F). The protein-protein internetwork (PPI) demonstrated that these m6A regulators were associated considerably in protein level (Fig. 1G). To sum, these results implied that the m6A regulators of RNA methylation play an essential role in ccRCC.

The distinguished subgroups obtained from cluster analysis
The unsupervised clustering result showed that the 503 patients could be categorized into three clusters with high stability when the cluster method was k-means and k = 3 using R package "ConsensusCluster" ( Fig. 2A-C). Then, the K-M survival analysis for m6A cluster subgroups showed that cluster 2 had a worse prognosis than the other two subgroups (Fig. 2D), indicating that cluster 2 differed from cluster 1/3. Therefore, we analyzed the expression of m6A regulators in these three subgroups. 13 Genes were significantly differentially expressed in subgroups (Fig. 2F). Moreover, IGF2BP1, IGF2BP2, and IGF2BP3 were highly expressed in cluster 2 (Fig. 2E). This result indicated that the subgroups screened via m6A regulators expression were closely associated with the progression and prognosis of ccRCC.
To further analyze cluster 2, we use the PCA algorithm to find the vital components in cluster 2. Figure   tors expressed in different grades, E the boxplot for m6A regulators expressed in different stages, F the correlations among m6A regulators were analyzed by Spearman correlation, and G a PPI network was established to examine the interaction of m6A RNA methylation regulators. ***P < 0.001, **P < 0.01, *P < 0.05, ns not significant showed that the result of each principal component contributes to the distribution, and dimension 1 contributes to the distribution most extensive. Figure 3B showed a two-dimensional image based on the principal component classification. We found that the difference between cluster 2 and cluster 1/3 was mainly in the principal component 2 (PC2), so we primarily analyzed the vital factors in principal component 2 (PC2). We found that IGF2BP2 and IGF2BP3 contributed to the top two in the principal component 2 (PC2) (Fig. 3C, D). These results indicated that IGF2BP2 and IGF2BP3 were important roles in cluster 2. The m6A cluster-related modules and genes Based on the above analysis, the m6A cluster 2 and IGF2BPs attracted our attention. We used the TCGA-KIRC dataset to construct the gene co-expression network to screen m6A cluster-related modules using WGCNA tool. We calculated the correlation coefficient of each gene with IGF2BP2 and IGF2BP3, then these genes with P-value < 0.05 were included in the WGCNA. The scale-free R 2 was set at 0.85, and the soft threshold was set at 7 (Fig. 4A, B). The selected genes were recognized as seven modules (Fig. 4C). As shown in Fig. 4D, E, the result demonstrated that the Brown module was significantly correlated with gene significance in cluster 2 (cor = 0.83, P < 0.001). To further understand the module genes, we utilized GO analysis to obtain an insight into the biological activities of 23 RNA methylation regulators and cluster-related genes. GO analysis demonstrated that the result of enrichment was related to nuclear activation (Fig. 4F). The GSEA showed that these cluster 2-associated genes were significantly associated with CD8+ T cell infiltration (Fig. 4G). These results suggested that m6A methylation regulators, especially IGF2BP2 and IGF2BP3, have a close relationship to immune cells and TIME.

The landscape of immune infiltration cells in ccRCC
The above analysis demonstrated that m6A cluster 2 is related to the immunological response. Consequently, we focused our research on the characteristics of TIME cell infiltration. The immune infiltration of ccRCC patients was obtained from ssGSEA algorism. The correlation between  (Fig. 5A). In addition, the immune cells were differently infiltrated in different clinical grades and pathological stages (Fig. 5B, C). These results demonstrated that tumor immune infiltrating cells were associated with clinical characters of ccRCC, and tumor immune cells can perhaps affect the progression of ccRCC. Fig. 4 Identification of the core modules related to the m6A cluster via WGCNA. A The scale-free topology model fit index was analyzed for various soft threshold powers (β), and the scale-free R 2 was set at 0.85. B Analysis of the mean connectivity for various soft threshold powers and the soft threshold was set at 7. C A cluster dendrogram was established to detect co-expression gene modules with corresponding color assignments in these genes related to IGF2BPs. The gray module refers to no co-expression between the genes. D Heatmap showing the correlation between six modules and three m6A cluster subgroups. E The scatterplot of gene significance for cluster 2 subgroup and module membership in the Brown module. F the result of GO analysis for the genes related to cluster 2. G the result of GSEA for the genes related to cluster 2.

The relationship between m6A regulators subgroups and immune infiltration cells
Based on the above analysis, the association between immune infiltration cells and m6A regulators is necessary to analyze. Our study revealed that cluster 2 was significantly highly exhibited CD4+ T cell, CD8+ T cell, MDSC, T follicular helper cell, Th1, and Th2 cells (P < 0.05, Fig. 6A). This result further showed that the m6A cluster was closely correlated with immune cell infiltration, especially cluster 2 strongly associated with several T cells. To further analyze the correlation between m6A regulators and the immune microenvironment, we used the ESTIMATE algorithm to calculate the ESTI-MATE score, immune score, and stomal score for each sample, assessing the purity of tumor samples. When comparing the purity of tumor samples in three m6A regulators clusters, we found out the cluster 2 exhibited a high ESTIMATE score (P < 0.001, Fig. 6B) and immune scores (P < 0.001, Fig. 6C), while cluster 2 exhibited no significant difference in the stromal score (P = 0.075, Fig. 6D). This result revealed a difference in immune infiltration

The relationship between m6A regulators subgroups and immunomodulator agonists
Antigen-presenting molecular and immunomodulators play an essential role in tumorigenesis. Therefore, we explore the expression of antigen-presenting genes and immunomodulatory genes in m6A cluster subgroups. The result showed that several antigen-presenting genes were significantly differentially expressed in m6A cluster subgroups. Compared with cluster 1/3, HLA-DBM, HLA-DOB, HLA-DPB2, HLA-DRA, and HLA-DRB6 were highly expressed in cluster 2 (Fig. 6E). Moreover, the expression of these antigen-presenting genes is consistent in cluster 1 and cluster 3, suggesting that cluster 2 has different immune characters compared with cluster 1/3. Previous studies have revealed that immune checkpoint was significantly correlated with the progression and prognosis of various cancers. We accessed the expression of several immunomodulator agonist genes between tumor samples and adjacent normal samples. The results showed that these immunomodulator agonist genes were differently expressed in tumor and normal samples (Fig. S2). Further, we analyzed the correlation between m6A clusters and immunomodulator agonist genes. The result showed that the LAG3, TIGIT, CTLA4, PD1, CD27, and PDL2 were significantly highly expressed in m6A cluster 2 (Fig. 6F). In summary, m6A cluster 2 is closely related to immunomodulatory factors in ccRCC, and cluster 2 perhaps has a better immune response than cluster 1/3.

The role of m6A regulators in the progression and prognosis of ccRCC
The above analyses showed that the m6A regulators were highly correlated with the TIME cells, and some tumorinfiltrating immune cells were associated with the prognosis of patients with ccRCC. Therefore, we explored the role of m6A regulators in the progression and prognosis of ccRCC. The univariate Cox analysis for 23 m6A regulators found that 16 genes were significant prognostic factors (Fig. 7A). ELAVL1, HNRNPA2B1, HNRNPC, IGF2BP1, IGF2BP2, IGF2BP3, YTHDF1, METTL3, and MTAP were adverse prognostic factors, while LRPPRC, YTHDC1, YTHDF1, FTO, CBLL1, KIAA1429, METTL14, and ZC3H13 were favorable prognostic factors. Furthermore, the risk factors have relatively consistent expression, and the correlation analysis among these significant prognostic factors was consistent with the univariate Cox analysis (Fig. 7B). Then the K-M survival analysis for IGF2BP2 and IGF2BP3 showed that these two genes were significantly correlated with the OS, DSS, and PFI of patients with ccRCC ( Fig. 7C-H), suggesting that these two m6A regulators play a vital role in the prognosis of ccRCC.

Mutation and copy number variation in m6A cluster-related subgroups
Next, we analyzed the mutation profiles and CNV of 23 m6A regulators in all samples and different m6A cluster subgroups. The result of CNV in all samples revealed that the most amplification is YTHDC2, and the highest frequency of CNV deletion is RBM15B (Fig. S3A). The result of the mutation in all samples indicated that the most amplification of m6A genes is VHL (Fig. S3B). Interestingly, we found that CNV and TMB in m6A cluster subgroups were different from all samples, and the result of these top highest mutation genes has some difference in different m6A cluster subgroups ( Fig. 8A-C). The mutation frequency of VHL gene was 54%, 41%, and 41% in clusters 1, 2, and 3, respectively. The mutation of PBRM1 (37%, 45%, 43%), SETD2 (12%, 29%, < 4%), BAP1 (5%, 20%, 12%), and TTN (12%, 20%, 14%) was significantly higher in the cluster 2, suggesting that the mutation in cluster 2 is different from the cluster 1/3. These different mutations among m6A cluster subgroups may lead to different immune responses.
The CNV analysis of three m6A cluster subgroups showed that the highest CNV frequency of m6A genes was IGF2BP3 in cluster 2 (Fig. 8E), higher than the CNV frequency in cluster 1/3 (Fig. 8D, F). Meanwhile, the CNV frequency of IGF2BP2 was about 40%, which is significantly higher than that in cluster 1/3. The result of CNV analysis indicated that CNV of m6A genes in cluster 2 might lead to different mutations among three groups. Further, these CNV and mutation differences among the m6A cluster subgroup influenced the immune infiltration.

The association between m6A regulators and immune cells
Based on the above analysis, cluster 2 has a close relationship to immune cells and immunomodulators. Specifically, IGF2BP2 and IGF2BP3 were significantly highly expressed in cluster 2 and worked as vital components in cluster 2. Therefore, we analyzed the association of these two genes with immune cells. Figure 9A, E shows that these two genes are closely related to several tumor-infiltrating immune cells in all samples, inferred that these two genes were associated with the immune response. To further explore the mechanism of these two genes in affecting immune cells and their role in different cluster subgroups, we analyzed the relationship between these two genes and immune cells in three m6A cluster subgroups, respectively. In cluster 1, IGF2BP2 and IGF2BP3 were significantly associated with several types of T cells (Fig. 9B, F). Interestingly, we found that IGF2BP2 was only positively correlated with activated CD4+ T cell in cluster 2 (Fig. 9C), while IGF2BP3 was positively correlated with the activated CD4+ T cell and central CD8+ T cell (Fig. 9G). The relationship found in cluster 2 was significantly different from cluster 1/3 (Fig. 9D, H).
To further understand the role of these two genes, we used univariate Cox analysis to explore the role of immune cells in the prognosis of patients with ccRCC. As is shown in Fig. 9I, ten types of immune cells were significantly correlated with the OS of patients with ccRCC. In these immune cells, activated CD4+ was a risk factor in prognosis, which was validated in K-M survival analysis (Fig. 9J). The correlation between IGF2BPs and T cells indicated that IGF2BPs might interfere and influence the infiltration of activated CD4+ T cells and CD8+ T cells in ccRCC, which further influence the progression of cancer cells.

Screening the genes related to IGF2BP2 in cluster 2
Based on the above analysis, IGF2BPs are crucial components of cluster 2 closely related to TIME. To further analyze the mechanism of IGF2BPs in affecting immune infiltration, we analyzed the expression of these two genes in ccRCC in protein level and found that IGF2BP3 was negative in IHC of ccRCC (Fig. S4A). Therefore, we focused on analyzing the IGF2BP2 and related genes. And we found that HMGA2 was most closely related to IGF2BP2 in cluster 2 (R = 0.5901, P < 0.001, Fig. 10A). According to published literature, aberrant alterations of the HMGA2 have a relationship with the occurrence of several types of tumors [42,43], which means HMGA2 potentially is an oncogene. In our study, we found that the expression of HMGA2 was different in tumor tissue and normal tissue (Fig. 10B), and HMGA2 was significantly higher expressed in cluster 2 subgroup than other two cluster subgroups (Fig. 10C), which indicated that HMGA2 affected the tumor progression in cluster 2. To validate the correlation between IGF2BP2 and HMGA2, we conducted the correlation analysis in GSE33689 and GSE53757. The results showed that HMGA2 was significantly related to IGF2BP2 in GSE33689 (R = 0.5695, P = 0.0013) (Fig. 10D) and in GSE53757 (R = 0.2322, 0.0497, Fig. 10E). These external validations confirmed the relationship between HMGA2 and IGF2BP2.
To further analyze the relationship between HMGA2 and IGF2BP2, we downloaded the IHC of HMGA2 and IGF2BP2 in renal cancer from HPA. Moreover, we found that HMGA2 is expressed mainly in the nucleus (Fig. 10F;  Fig. S4B), IGF2BP2 is principally expressed in the cytoplasm and cell membrane ( Fig. 10G; Fig. S4B). Meanwhile, we found that HMGA2 and IGF2BP2 have corresponding expressions in the same patients, confirming a specific positive correlation between HMGA2 and IGF2BP2. Based on the expression location of these two genes and the relationship between these two genes, we speculated that HMGA2 regulates the expression of IGF2BP2 in the nucleus, then IGF2BP2 affects the ccRCC tumor cells proliferation.

Discussion
Kidney cancer is one of the most diagnosed cancer diseases, and ccRCC accounts for the majority of kidney cancer [1,44]. The traditional treatment for kidney cancer is conducting surgery to remove the tumor node. Nowadays, many studies have shown that kidney cancer is associated with genetic alteration [45]. Therefore, understanding kidney cancer at the molecular genetic level will be more beneficial to treating kidney cancer. m6A modification is the most commonly encountered alteration in eukaryotic mRNA. Although numerous studies have shown that m6A plays an indispensable role in the immune response and anti-tumor effect through methylation [21][22][23]46], the interaction between the TIME and m6A regulators in ccRCC has not been comprehensively understood. Therefore, it is essential to identify the role of m6A in TIME of ccRCC, which will be beneficial to develop novel targets for immunotherapy strategy.
In this work, we found that most m6A genes were differently expressed in ccRCC compared with the normal tissue of the kidney, and most m6A genes were significantly differentially expressed in clinical characters groups, suggesting that these m6A regulators played a crucial role in the tumor progression. In several previous studies, researchers found FLG  DNAH9  COL5A3  ARID1A  AKAP9  USH2A  SCAF4  MYH2  MUC4  FBN2  ATM  AHNAK2  BAP1  HMCN1  MUC16  KDM5C  SETD2  TTN  PBRM1  VHL   TP53  NOTCH2  LRP2  KIAA1549L  FAM178A  FAM135A  COL1A2  ATM  XIRP2  TRIOBP  SPEN  RANBP2  CSMD3  DNAH9  MTOR  TTN  BAP1  SETD2  VHL  PBRM1   RYR3  HYDIN  DST  ATM  ALMS1  DAMTS12  MACF1  LRP2  KMT2D  HMCN1  ANK3  ROS1  KMT2C  DNAH9  ARID1A  MTOR  BAP1  TTN  VHL  that the expression level of m6A regulators was associated with the progression of tumors, and these m6A regulators can be oncogene or tumor suppressors [14]. Ma et al. showed that the low expression level of METTL14 has a worse prognosis in liver cancer patients [47], while, as shown in the studies of Chen et al., METTL3 has a positive correlation with progress liver cancer progression, and METTL3 promotes the tumor progression via YTHDF2-dependent posttranscriptional silencing of SOCS2 [48]. Although METTL4 and METTL3 were m6A "readers," the effect of these m6A readers on tumor progression was different, indicating the complex mechanism of the m6A regulating process. Therefore, our study found that m6A regulators were differently expressed in tumor samples and normal samples, which was consistent with the previous studies. Moreover, the role of m6A in ccRCC deserved to be further analyzed.
We divided the samples into three subgroups via unsupervised clustering based on the expression of m6A regulators. Meanwhile, we calculated the abundance ratio of tumorinfiltrating immune cells via ssGSEA algorithm. Then, we found that the relative abundance of immune cells was distinctly distributed among these three subgroups. Especially in cluster 2, T cells were highly infiltrated compared with the other two groups. Interestingly, the cluster 2 subgroup has worse OS than cluster 1/3, indicating that the different distribution of immune cells contributes to the worse prognosis.
The functional enrichment analysis showed that the genes correlated with IGF2BPs in cluster 2 were enriched in some immune response and regulating T cells, suggesting that T cells in TIME played a crucial role in the progression of ccRCC. To further comprehend the relationship between m6A regulators and TIME, we analyzed the correlation between the m6A regulators and immune modulators. This procedure indicated that the immune modulators were differently expressed in m6A cluster groups, indicating that the m6A regulators perhaps affected the immune response in tumor progression. Next, the IGF2BP2 and IGF2BP3 were identified as essential genes in cluster 2, and these two genes were significantly correlated with activated CD4+ T cells and central memory CD8+ T cells compared with the other Fig. 10 The relationship between IGF2BP2 and HMGA2. A The correlation between IGF2BP2 and HMGA2 in TCGA-KIRC dataset. B The box plot for the expression of HMGA2 in tumor samples and adjacent normal tissue in ccRCC. C The box plot for the expression of HMGA2 in three m6A cluster subgroups. D and E The correlation between IGF2BP2 and HMGA2 in GSE36895 dataset and GSE53757 dataset. F The IHC of HMGA2 in renal cancer. G The IHC of IGF2BP2 in renal cancer two clusters. Therefore, the alteration of these two genes may affect the correlated immune cells to influence the progression of ccRCC.
We learn from previous articles that m6A molecular may have a regulatory effect on immune cell infiltration. Han et al. revealed that YTHDF1 could affect the tumor antigen presentation ability and immune recognition ability of DCs to influence the TIME [22]. This effect of m6A regulators was also observed in several studies [23,24]. Zhang et al. found that after knocking down METTL14/YTHDF1 in gastric cell lines, the transcript levels of IFN-α/β were downregulated [19]. IFN-α/β is critical for tumor cell suppression and anti-tumor immune stimulation through promoting and balancing the anti-inflammatory response [49]. Therefore, m6A readers, METTL14, could regulate the immune microenvironment of the tumor cell to affect the progression of gastric cancer. As shown in these studies, the effect of m6A regulators on immune cells infiltration in the tumor microenvironment is a crucial factor that cannot be ignored in the progression of the tumor.
Our study found that high expression of IGF2BP family genes (IGF2BP2 and IGF2BP3) was significantly correlated with the poor OS, DPI, and PFI of ccRCC patients. Meanwhile, the expression of IGF2BPs was positively related to the activated CD4+ T cells and central memory CD8+ T cells in cluster 2. Further, our study found that IGF2BP2 is closely correlated with HMGA2 in cluster 2, and HMGA2 and IGF2BP2 have corresponding expressions in the patients. In previous studies, Jan R. et al. found that IGF2BP2 was inhibited in HMGA2-null mice, suggesting that HMGA2 regulated IGF2BP2 in mice [50]. In our research, the HMGA2 was highly expressed in the nucleus, which indicated that HMGA2 could also regulate the IGF2BP2 in ccRCC. Zhou et al. have revealed that the IGF2BPs enhance the stability and translation of their target mRNAs [51]. Some studies also showed that IGF2BPs proteins post-transcriptionally regulate the stability of target mRNAs to affect the progression of tumor cells, in which the MYC is the essential target [52,53]. Huang et al. have revealed that IGF2BPs regulate the MYC expression by binding to MYC mRNA and affect the cancer cell proliferation [53]. Our research can infer that m6A regulators can impact the ccRCC tumor cell by regulating the HMGA2 expression. Therefore, taking our research sum, IGF2BP2 and IGF2BP3 could affect the tumor microenvironment through multiple functional pathways, and these two genes and their correlated genes could be potential immunotherapy targets in ccRCC patients.
Inevitably, there are several limitations of this study that are supposed to discuss. First, the gene expression data were extracted from the TCGA and GEO databases, in which most of the patients were white people. Therefore, the data will exert some selection of bias of the population. Secondly, we comprehensively analyzed the correlation between m6A regulators and TIME; however, the mechanism of how m6A molecular affected the immune modulators and immune cells did not clarify. The vitro or vivo experiments could be conducted to analyze further how m6A affects the ccRCC.

Conclusion
In conclusion, our work demonstrated a close relationship between m6A regulators and immune infiltrating cells in ccRCC. The difference of m6A modification subgroups has some effect on the regulating tumor microenvironment. Further, IGF2BP2 and IGF2BP3 are essential m6A regulators that cannot be neglected in affecting the TIME of ccRCC.
have obtained ethical approval. Our study is based on open-source data, so there are no ethical issues and other conflicts of interest.

Consent for publication
The UCSC Xena is a public database and patients included in the database have provided informed consent for the use of their data and for publication.