Identifying Prognostic Biomarkers Related to m6A Modification and Immune Infiltration in Renal Cell Carcinoma

Background: Renal cell carcinoma (RCC) is the largest category of kidney tumors and usually does not have a good prognosis. N6-methyladenosine(m6A) and immune infiltration have received increased attention because of their great influence on the clinical outcome and prognosis of cancer patients. Methods: We identified hub genes through multi-dimensional screening, including DEGs, PPI analysis, LASSO regression, and random forest. Meanwhile, GO/KEGG enrichment, cMAP analysis, prognostic analysis, m6A prediction, and immune infiltration analysis were performed to understand the potential mechanism and screen therapeutic drugs. Results: We screened 275 downregulated and 185 upregulated genes using three GEO datasets and the TCGA dataset. In total, 82 candidate hub genes were selected using STRING and Cytoscape. Enrichment analysis illustrated that the top 3 biological process terms and top 1 KEGG term were related to immunity. cMAP analysis showed some antagonistic molecules can be candidate drugs for the treatment of RCC. Then, six hub genes (ERBB2, CASR, P2RY8, CAT, PLAUR, and TIMP1) with strong predictive values for prognosis and clinicopathological features were selected. Meanwhile, P2RY8, ERBB2, CAT, and TIMP1 may obtain m6A modification by binding METTL3 or METTL14. On the other hand, differential expression of CAT, ERBB2, P2RY8, PLAUR, and TIMP1 affects the infiltration of the majority of immune cells. Conclusions: We identified six hub genes through multi-dimensional screening. They all possess strong predictive value for prognosis and clinicopathological features. Meanwhile, hub genes may regulate the progression of RCC via an m6A- and immunity-dependent mechanism.


Introduction
Renal cell carcinoma (RCC) is a common subtype of kidney tumor, with approximately 73,000 new cases and 14,000 deaths each year [1]. Considering its association with an unhealthy lifestyle, genetic alteration, and cell injury, RCC is a complex, multifactorial illness. Despite huge developments in diagnosis, surgery, and medication, RCC still has an unsatisfactory outcome [2]. The prognosis is poor, especially for advanced RCC patients whose 5-year survival rate is less than 5% because of insensitivity to chemotherapy and other treatments [3]. Therefore, more efficient targets for diagnosis and treatment are urgently needed. m6A (N6-methyladenosine) is the most common and reversible RNA modification in eukaryotic cells. It is a dynamic modification regulated by methyltransferases, demethylases, and RNA-binding proteins [4]. Increasing evidence has revealed that m6A modification is associated with several biological processes, such as cell proliferation, apoptosis, tumor progression, damage repair, and immunomodulatory disorders [5]. m6A-related enzymes can alter the m6A levels of cancer inhibitors or promoters affecting tumor progression, such as BRD4, MYC, SOCS2, and EGFR [6].
On the other hand, the tumor immune microenvironment (TIME) has received much attention because of its great influence on the clinical outcome and prognosis of cancer  [7]. Thus, the infiltration of immune cells is considered to be an appealing target for diagnosis and treatment [8]. Meanwhile, some genes can promote or inhibit the infiltration of immune cells. For example, the TBC1D3 family can function as a prognostic biomarker of RCC, and the expression of TBC1D3 has a negative correlation with the infiltrated level of CD4+ T cells [9]. Additionally, TREM2 is regarded as a prognostic marker in several cancers, and the expression of TREM2 in these tumors is positively related to the infiltration level of most immune cells [10].
In this study, we explored six promising biomarkers that have a good correlation with the survival prognosis and clinicopathological features of RCC. We also analyzed their relevance to immune infiltration and m6A modification. Moreover, we presented a preliminary landscape of immune cell infiltration and m6A-related genes in RCC.

Data Collection
Transcription profile data were downloaded from the GEO database with the following criteria. First, the dataset should contain patients and controls. Second, patients and controls of each group should include more than 3 tissue samples from homo sapiens. Third, the data type is expression profiling by array. RNA-seq data of ccRCC and related clinical information were downloaded from the TCGA website. In this study, three GEO datasets and the TCGA dataset were applied to DEG (differentially expressed gene) analysis:

DEG Identification
The DEGs between RCC and normal tissues were identified using the R package "limma" with the threshold set at FoldChange > 2 and adjusted p < 0.05. The results were visualized by the R package "ggplot2". Then, the intersection of the three GEO datasets and the TCGA dataset was calculated and displayed in the form of a Venn diagram.

PPI (Protein-Protein Interaction) Network Construction and Candidate Hub Gene Identification
The PPI network, predicting the functional interactions between genes, could provide insights into the molecular mechanism. The STRING online database (http://string-db.org accessed on 26 September 2022) [11] was utilized to construct the PPI network. Cytoscape is a useful software capable of visualizing molecular interaction networks [12]. MCODE, a plugin of Cytoscape, was used to extract the densely connected subnetwork with the threshold set at MCODE default parameters (node score cutoff ≥ 0.2, degree cutoff ≥ 2, K-core ≥ 2, and max depth = 100). CytoHubb, another plugin of Cytoscape, was applied to identify the top 100 proteins ranked by the MCC algorithm. Then, the intersected molecules were recognized as the candidate hub genes.

GO Enrichment and KEGG Pathway Analysis of Candidate Hub Genes
GO enrichment is a useful bioinformatics method to analyze biological functions, including biological processes, molecular functions, and cellular components. KEGG pathway analysis is also widely utilized to provide functional annotations for a gene set including biological pathways, drugs and chemical substances, and diseases. In this study, we performed GO and KEGG enrichment analyses by the R packages "ggplot2" and "clusterprofiler" with the threshold set at adjusted p < 0.05.

cMAP Analysis of Candidate Hub Genes
The cMap database is the largest perturbation-driven gene expression dataset and is a resource using cellular responses to perturbation to find relationships between diseases, genes, and therapeutics (https://clue.io/, accessed on 26 September 2022) [13]. In this study, antagonistic molecules were screened for RCC by cMAP analysis of 82 candidate hub genes. The query parameters selected were "Gene expression (L1000)" and "Latest". The top 10 ranked molecules and relevant information are listed in Table 1.

Hub Gene Identification
LASSO is a frequently used algorithm to screen the most valuable genes from highdimensional data. Random forest is a supervised learning method. An important property of a Random forest is that it can give an importance measure for every feature. So it was widely used for the extraction of key features from a large gene set [14]. In our research, LASSO and random forest were performed together to select the most powerful variables for survival prognosis. LASSO was performed by "glmnet" package in R. The lambda (λ) value of LASSO was set at lambda.min after ten rounds of cross-validation. Random forest was performed by "randomForestSRC" package in R. "ntree" set at 1000, "nsplit" set at 1, "method" set at "md".

Survival and Clinicopathological Characteristic Analysis of Hub Genes
The overall survival and disease-free survival analyses of the hub genes were performed using the GEPIA website (http://gepia.cancer-pku.cn/, accessed on 26 September 2022) [15]. The correlation analysis between the expression of hub genes and clinicopathological features was carried out by UALCAN (http://ualcan.path.uab.edu/, accessed on 26 September 2022) [16].

Immune Infiltration Analysis
CIBERSORTX is a powerful method for characterizing cell subsets by gene expression profiles [21]. In our study, the CIBERSORTx website (https://cibersortx.stanford.edu/, accessed on 26 September 2022) was used to estimate the abundances of 22 immune cells between cancer tissue and normal tissue, high-and low-expression groups of hub genes, and different clinicopathological feature groups (t1-t2 vs. t3-t4, n0 vs. n1-n3, m0 vs. m1, grade1-grade2 vs. grade3-grade4). Gene expression profiles and Reference gene expression were needed to input into the website. Gene expression profiles came from the ccRCC and control cohort of TCGA. Reference gene expression signatures were LM22 (22 immune cell types) provided by the website. The permutation for significance analysis chose 100. The median value of hub genes was defined as the threshold of high and low expression. Then, R packages ("ggplot2", "vioplot", "corrplot", "survial") were used to visualize the results. p < 0.05 was considered to be significantly different.

DEG Identification in RCC
The workflow of this research is shown in Figure 1. There were 708, 2037, 996, and 1183 upregulated genes and 927, 2442, 1308, and 1034 downregulated genes in the GSE68417, GSE71963, GSE168845, and TCGA datasets, respectively. The DEGs of the three GEO datasets and TCGA dataset are illustrated by volcano maps (Figure 2A

PPI Interaction Network and Identification of Candidate Hub Genes
STRING and Cytoscape were applied together to identify candidate hub genes. Fourteen clusters with 144 genes (gene list 1) were selected using the plugin MCODE of Cytoscape. The MCODE scores ranged from 3.0 to 13.742. The largest cluster contains 32 genes with an MCODE score of 13.74. The smallest cluster has 3 genes with an MCODE score of 3.0. The clusters with MCODE scores > 5.0 are shown in Figure 3. The top 100 genes (gene list 2) ranked by the MCC algorithm were selected using the plugin cytoHubb of Cytoscape. Eighty-two genes overlapping in gene list 1 and gene list 2 were identified as candidate hub genes. Gene list1, gene list2, and candidate hub gene list were displayed in Supplementary Table S1.

GO and KEGG Enrichment Analysis
To understand the potential biological functions of DEGs, GO and KEGG enrichment analyses were performed by the R package "clusterprofiler". In total, 35 KEGG pathways, 4 molecular function terms, the top 10 biological process terms, and the top 10 cellular component terms are shown in Figure 4A,B. The molecular function terms included integrin binding, pattern recognition receptor activity, amyloid−beta binding, and lipoteichoic acid binding. The biological process terms were mainly enriched in lipoteichoic acid binding, T-cell activation, and phagocytosis. The cellular component terms mainly included synapse pruning, platelet alpha granule lumen, and platelet alpha granule lumen. The enriched KEGG pathways were mainly involved in complement and coagulation cascades, Staphylococcus aureus infection, leishmaniasis, and tuberculosis. Meanwhile, the top 3 biological process terms and the top 1 KEGG term were related to immunity. This enrichment result verifies previous research results that immunity is one of the major mechanisms of candidate hub genes [22]. Three GEO datasets and ccRCC cohort of TCGA were included in present study. Through differential gene analysis, PPI, LASSO regression, and random-forest analysis, 6 hub genes were identified. Then, survival analysis, correlation analysis with m6A-related genes, immune infiltration analysis, and association with clinicopathological features were performed to understand 6 hub genes further.

Figure 1.
Schematic sketch of the research design. Three GEO datasets and ccRCC cohort of TCGA were included in present study. Through differential gene analysis, PPI, LASSO regression, and random-forest analysis, 6 hub genes were identified. Then, survival analysis, correlation analysis with m6A-related genes, immune infiltration analysis, and association with clinicopathological features were performed to understand 6 hub genes further.

PPI Interaction Network and Identification of Candidate Hub Genes
STRING and Cytoscape were applied together to identify candidate hub genes. Fourteen clusters with 144 genes (gene list 1) were selected using the plugin MCODE of Cytoscape. The MCODE scores ranged from 3.0 to 13.742. The largest cluster contains 32 genes with an MCODE score of 13.74. The smallest cluster has 3 genes with an MCODE score of 3.0. The clusters with MCODE scores > 5.0 are shown in Figure 3. The top 100 genes (gene list 2) ranked by the MCC algorithm were selected using the plugin cytoHubb of Cytoscape. Eighty-two genes overlapping in gene list 1 and gene list 2 were identified as candidate hub genes. Gene list1, gene list2, and candidate hub gene list were displayed in Supplementary Table S1.

cMAP Analysis of Candidate Hub Genes
After cMAP analysis of 82 candidate hub genes, we found their expression treated by some molecular drugs was basically contrary to that induced by ccRCC. Thus, they can be candidate drugs for the treatment of RCC. The top 10 molecules ranked by the raw score are shown in Table 1. For example, the top molecule is miglustat. It can target UGCG, GAA, and SLC5A4 and act as a glycosylation inhibitor.

Hub Gene Identification
Eighty-two candidate hub genes are too many to guide clinical work. Therefore, we further screened variables using the RandomForest algorithm and LASSO regression. Eighteen genes (gene list3) were selected by the random forest method, and six genes (gene list4) were selected using LASSO regression. The LASSO screening process is visualized in Figure 5. Finally, we integrated the results of LASSO regression and RandomForest to identify six hub genes: ERBB2, CASR, P2RY8, CAT, PLAUR, and TIMP1. Gene list3 and gene list4 were displayed in Supplementary Table S1. To understand the characterized functions of six hub genes, a preliminary summary performed using Genecard (https://www.genecards.org/, accessed on 26 September 2022) was supplied in Supplementary Table S2.

GO and KEGG Enrichment Analysis
To understand the potential biological functions of DEGs, GO and KEGG enrichment analyses were performed by the R package "clusterprofiler". In total, 35 KEGG pathways, 4 molecular function terms, the top 10 biological process terms, and the top 10 cellular component terms are shown in Figure 4A,B. The molecular function terms included integrin binding, pattern recognition receptor activity, amyloid−beta binding, and lipoteichoic acid binding. The biological process terms were mainly enriched in lipoteichoic acid binding, T-cell activation, and phagocytosis. The cellular component terms mainly included synapse pruning, platelet alpha granule lumen, and platelet alpha granule lumen. The enriched KEGG pathways were mainly involved in complement and coagulation cascades, Staphylococcus aureus infection, leishmaniasis, and tuberculosis. Meanwhile, the top 3 biological process terms and the top 1 KEGG term were related to immunity. This enrichment result verifies previous research results that immunity is one of the major mechanisms of candidate hub genes [22].

Survival and Clinicopathological Characteristic Analysis of Hub Genes
Overall survival and disease-free survival analyses were conducted using the GEPIA web ( Figure 6). The results demonstrated that low expression of TIMP1 and PLAUR and high expression of CASR, CAT, ERBB2, and P2RY8 were associated with positive overall survival and disease-free survival rates. Clinicopathological characteristic analysis of the hub genes was carried out by the ULUCAN web ( Figure 7). CASR was highly expressed in normal tissues and expressed at extremely low levels in each clinical subgroup in cancer tissues ( Figure 7A). TIMP1 and PLAUR were expressed at low levels in normal tissues. With increasing clinical stage, histological grade, and the number of metastatic lymph nodes, the expression increased ( Figure 7D,F). CAT and ERBB2 were highly expressed in normal tissues. With increasing clinical stage, histological grade, and the number of metastatic lymph nodes, the expression decreased ( Figure 7B,C). P2RY8 was expressed at low levels in normal tissues. However, we were surprised to find that with increasing

cMAP Analysis of Candidate Hub Genes
After cMAP analysis of 82 candidate hub genes, we found their expression treated by some molecular drugs was basically contrary to that induced by ccRCC. Thus, they can be candidate drugs for the treatment of RCC. The top 10 molecules ranked by the raw score are shown in Table 1. For example, the top molecule is miglustat. It can target UGCG, GAA, and SLC5A4 and act as a glycosylation inhibitor.

Hub Gene Identification
Eighty-two candidate hub genes are too many to guide clinical work. Therefore, we further screened variables using the RandomForest algorithm and LASSO regression. Eighteen genes (gene list3) were selected by the random forest method, and six genes (gene list4) were selected using LASSO regression. The LASSO screening process is visualized in Figure 5. Finally, we integrated the results of LASSO regression and Random-Forest to identify six hub genes: ERBB2, CASR, P2RY8, CAT, PLAUR, and TIMP1. Gene list3 and gene list4 were displayed in Supplementary Table S1. To understand the characterized functions of six hub genes, a preliminary summary performed using Genecard (https://www.genecards.org/) was supplied in Supplementary Table S2.

Immune Infiltration Analysis
To further determine the effect of immune infiltration in RCC and the relationship with hub genes, we performed a further analysis in three steps. First, we aimed to analyze immune infiltration in RCC and normal tissue. Stacked bar charts exhibit the proportions of 22 kinds of immune cells in normal tissue ( Figure 8A) and cancer tissue ( Figure 8B). A large proportion of M1 macrophages, M2 macrophages, CD8 T cells, and resting memory CD4 T cells were observed in the tumor group. Resting memory CD4 T cells, M2 macrophages, monocytes, and resting dendritic cells accounted for a large proportion of cells in normal tissue. B naive cells, resting memory CD4 T cells, resting dendritic cells, and resting mast cells were more highly infiltrated in normal samples ( Figure 8C). CD8 T cells, follicular helper T cells, regulatory T cells, gamma delta T cells, M0 macrophages, M1 macrophages, M2 macrophages, and activated mast cells showed decreased infiltration in tumor tissues ( Figure 8C). A hierarchical cluster plot also showed the state of immune cell infiltration in the tumor and normal groups ( Figure 8D).
Second, the immune infiltration between different groups of clinicopathological features was compared ( Figure 9A-E). Meanwhile, p values were shown together in a heatmap, and red indicated that the immune infiltration was significantly different. The more significant the difference, the darker the color ( Figure 9I). Resting mast cells, regulatory T cells, and M2 macrophages had a greatly different infiltration between different situations of the primary-tumor stage, pathological grade, tumor metastasis, and clinical stage ( Figure 9I). T-cell regulation was a positive predictor for the prognosis of RCC ( Figure 9F). Resting mast cells and resting dendritic cells indicated a poor prognosis ( Figure 9G,H).
Finally, we analyzed the relationship between the expression of hub genes and the infiltration of immune cells (Figure 10). The analysis indicated that differential expression of CAT, ERBB2, P2RY8, PLAUR, and TIMP1 affected the infiltration of the majority of immune cells (Figure 10B-F). p values were also shown together in a heatmap ( Figure 10G).

Survival and Clinicopathological Characteristic Analysis of Hub Genes
Overall survival and disease-free survival analyses were conducted using the GEPIA web ( Figure 6). The results demonstrated that low expression of TIMP1 and PLAUR and high expression of CASR, CAT, ERBB2, and P2RY8 were associated with positive overall survival and disease-free survival rates. Clinicopathological characteristic analysis of the hub genes was carried out by the ULUCAN web (Figure 7). CASR was highly expressed in normal tissues and expressed at extremely low levels in each clinical subgroup in cancer tissues ( Figure 7A). TIMP1 and PLAUR were expressed at low levels in normal tissues. With increasing clinical stage, histological grade, and the number of metastatic lymph nodes, the expression increased ( Figure 7D,F). CAT and ERBB2 were highly expressed in normal tissues. With increasing clinical stage, histological grade, and the number of metastatic lymph nodes, the expression decreased ( Figure 7B,C). P2RY8 was expressed at low levels in normal tissues. However, we were surprised to find that with increasing clinical stage, histological grade, and the number of metastatic lymph nodes, the expression decreased ( Figure 7E).

Correlation Analysis between Hub Genes and m6A-Related Genes
To understand the relationship between hub genes and m6A-related genes in RCC, we first analyzed the expression of m6A-related genes in tumor and normal tissues. The analysis indicated that FTO, ALKBH5, YTHDC2, WTAP, and RBM15 were highly expressed in the tumor group, and METTL14, METTL3, LRPPRC, HNRNPA2B1, FMR1, ZC3H13, YTHDF3, YTHDF2, and YTHDC1 were more highly expressed in normal tissues ( Figures 11B and 12). Kaplan-Meier plots of overall survival demonstrated that higher expression of YTHDC1, YTHDC2, RBM15, METTL14, LRPPRC, FTO, FMR1, ZC3H13, YTHDF3, YTHDF2, WTAP, and YTHDF1 was associated with increased OS (Figure 13). Then, the correlation between hub genes and m6A-related genes was determined.
Detailed information is shown in Table 2. The numbers represent the strength of the correlation, and red, blue, and gray colors indicate positive, negative, and nonsignificant correlations, respectively. Circos visualized the correlation, which showed that P2RY8, ERBB2, CAT, and TIMP1 had a close correlation with m6A-related genes ( Figure 11A). Additionally, several high and very confidently predicted m6A sites were found on the mRNAs of the hub genes, and the specific sequence and location information are shown in Figure 11C. Furthermore, to understand the potential interactive pattern, we drew a predicted interactive diagram using data from the ENCORI and catPARID databases ( Figure 11D). The results indicated that P2RY8, ERBB2, CAT, and TIMP1 may be methylated by METTL3 or METTL14 and then recognized by multiple RNA-binding proteins.  There was a close correlation between the expression of P2RY8, ERBB2, CAT, TIMP1, and m6A-related genes. (B) The heatmap shows the expression of m6A-related genes between tumor and normal tissues. Among the 16 m6A-related genes, they all had differential expression, except KIAA1429 and YTHDF1. * means p < 0.05, ** means p < 0.01, *** means p < 0.001. (C) The distribution of the predicted m6A sites using the SRAMP web server based on the sequence data from Nucleotide-NCBI. Red indicates that the confidence of the predicted m6A sites is very high, and yellow indicates that the confidence is high. The numbers below the line depict the position of m6A sites. "//" represent "Start codon" and "Stop codon". (D) The predicted diagram of the interactive pattern. The interactive diagram was made using data from the ENCORI and catPARID databases. First, we queried the interactive possibility of writers/erasers/readers with hub genes in ENCORI. Then, the interactive score was confirmed and obtained in the catPARID database. The thickness of the line represents the interactive possibility.
tide-NCBI. Red indicates that the confidence of the predicted m6A sites is very high, and yellow indicates that the confidence is high. The numbers below the line depict the position of m6A sites. "//" represent "Start codon" and "Stop codon". (D) The predicted diagram of the interactive pattern. The interactive diagram was made using data from the ENCORI and catPARID databases. First, we queried the interactive possibility of writers/erasers/readers with hub genes in ENCORI. Then, the interactive score was confirmed and obtained in the catPARID database. The thickness of the line represents the interactive possibility.   Detailed information is shown in Table 2. The numbers represent the strength of the correlation, and red, blue, and gray colors indicate positive, negative, and nonsignificant correlations, respectively. Circos visualized the correlation, which showed that P2RY8, ERBB2, CAT, and TIMP1 had a close correlation with m6A-related genes ( Figure 11A). Additionally, several high and very confidently predicted m6A sites were found on the mRNAs of the hub genes, and the specific sequence and location information are shown in Figure 11C. Furthermore, to understand the potential interactive pattern, we drew a predicted interactive diagram using data from the ENCORI and catPARID databases (Figure 11D). The results indicated that P2RY8, ERBB2, CAT, and TIMP1 may be methylated by METTL3 or METTL14 and then recognized by multiple RNA-binding proteins.

Discussion
Some RCC patients, especially advanced cases, usually do not have a good prognosis [23]. The lack of efficient molecular medication is one of the reasons hindering RCC treatment. Our study identified the potential biomarkers and therapeutic targets. The possible molecular mechanisms were exploited at the same time.
In our study, we screened 275 downregulated and 185 upregulated genes using 3 GEO datasets and the TCGA dataset. To choose robust biomarkers, protein-protein interaction (PPI) analysis was performed using the STRING web server. Additionally, the PPI network was further processed through the plugin MCODE and cytoHubba of Cytoscape. Eightytwo genes were screened out. To understand the potential mechanism, GO and KEGG enrichment analyses were carried out. The results illustrated that the top 3 biological process terms and the top 1 KEGG term were related to immunity. RCC is an immunogenic tumor [24]. Meanwhile, through cMAP analysis of 82 candidate hub genes, we found some antagonistic molecules that can be candidate drugs for the treatment of RCC. LASSO and random forest were employed together to select the most powerful variables for prognosis. Finally, six hub genes (ERBB2, CASR, P2RY8, CAT, PLAUR, and TIMP1) were identified. To further determine the predictive power and possible mechanisms of hub genes, multi-dimensional analysis was performed, including survival analysis, association with clinicopathological features, correlation analysis with m6A-related genes, and immune infiltration analysis.
Survival analysis demonstrated that low expression of TIMP1 and PLAUR and high expression of CASR, CAT, ERBB2, and P2RY8 were associated with a good prognosis, including OS and DFS. TIMP1 and PLAUR had a higher expression trend in the better subtype groups of clinicopathological features. CASR, CAT, ERBB2, and P2RY8 had a higher expression trend in the poor subtype groups of clinicopathological features. These analysis results suggest that the six hub genes play an important role in the progression of RCC. This consequence is partially consistent with previous studies, such as those involving ERBB2, which is rarely expressed in renal cell cancer and is involved in RCC oncogenesis [25]. CASR plays a vital role in the mechanism of bone metastasis in RCC [26]. TIMP1 indicates a poor prognosis of RCC and accelerates tumorigenesis [27]. However, P2RY8, PLAUR, and CAT have rarely been researched. In the present study, we found the expression of P2RY8 was interesting. It is expressed higher in tumors. However, with the clinical stage, histological grade, and the number of metastatic lymph nodes increasing, its expression decreased. This seems to be a paradox. In fact, a similar phenomenon does exist in other studies. For example, obesity is a clear risk factor for RCC; however, obesity can also improve the prognosis of RCC [28]. m6A occurs in the sixth-N-position of adenosine, which is the most abundant, prevalent, and conserved RNA modification in higher eukaryotes. Accumulating studies have indicated that m6A modification regulates tumor pathogenesis and progression in ccRCC. For example, METTL3 had a higher expression in ccRCC tissues and enhanced the cell viability, migration, and invasion abilities via regulating m6A modification of HHLA2 [29]. FTO inhibited autophagy of ccRCC cells and promoted tumor progression through an m6A-IGF2BP2-dependent mechanism, and SIK2 was their m6A target [30]. On the other hand, YTHDF2 recognized METTL14-mediated m6A markers on NEAT1_1 to accelerate the degradation of NEAT1_1, thus inhibiting the proliferation and migration ability of RCC cells [31].
To determine whether the mechanism of the six hub genes in RCC is related to m6A modification, a detailed analysis was performed. There are many m6A-related genes, and we chose 16 representative genes based on the literature [17]. According to our research, five m6A-related genes were highly expressed in the tumor group, and nine m6A-related enzymes had higher expression in normal tissues. High expression of twelve m6A-related genes is associated with increased OS. These results are consistent with a previous study [32]: Jiawu Wang indicated most of the m6A regulators have different expressions in RCC tissue, and METTL14 and METTL3 were identified as two powerful independent prognostic regulators after further analysis. Their results illustrate that m6A modification does play an important role in RCC. In our study, we found the expression of P2RY8, ERBB2, CAT, and TIMP1 had a close correlation with the expression of most m6A regulators. Additionally, several high and very confidently predicted m6A sites were found on the mRNAs of the hub genes. Finally, we drew a possible interactive diagram according to the data from ENCORI and catRAPID. These results all strongly support that ERBB2, CAT, P2RY8, and TIMP1 can regulate the development of RCC in an m6A-dependent manner. By searching the PubMed website, we found some similar m6A studies on these four hub genes. In esophageal squamous cell carcinoma, ERBB2 was identified as the target of FTO using m6A-seq and RNA-seq assays, and YTHDF1 then binds and stabilizes ERBB2 mRNA. Finally, the expression of ERBB2 increases and promotes cancer [33]. TIMP1 is more highly expressed in high-grade glioma and highly correlates with m6A markers in glioma, and TIMP1 serves as a potential biomarker for prognosis [34]. However, there are few RCCrelated studies. Regarding P2RY8 and CAT, we found few m6A-relevant studies. However, we must explain as follows: 1. These results only reveal underlying m6A mechanisms with hub genes in ccRCC. In our research, one hub gene can simultaneously promote/inhibit the expression of methylase and demethylase. However, we cannot simply conclude they increase/decrease the m6A level of ccRCC or promote/inhibit tumor progression because one methylase can modify oncogenes and anti-oncogenes to play a synergistic role. Additionally, one m6A regulator can promote proliferation and inhibit invasion and migration, which seem to be a reverse phenotype in a tumor [35]. In addition, the m6A gene can act as an independent prognostic factor rather than an m6A regulator [36]. 2. Although m6Asites localize primarily around 5 and 3 UTR [37], the study also indicates about 35% m6A site in coding regions (CDS) [38]. These results are consistent with our findings. It is worth mentioning that CDS m6A also play an important regulatory role. For example, it can disrupt tRNA selection and slow down translation elongation [39]. Additionally, research showed it induces ribosome pausing in a codon-specific way; removing CDS m6A leads to a further decrease in translation unexpectedly [38].
Tumor-infiltrating immune cells (TIICs) in the tumor microenvironment greatly determine tumor progression; they are promising targets for drug therapy and have demonstrated potential therapeutic value [8]. Our study indicated that TIICs are differently distributed in different subgroups of clinicopathological features. The data showed that M2 macrophages and CD8 T cells obviously infiltrated the tumor group, and this phenomenon is consistent with a previous study that showed that RCC is mainly infiltrated by T cells and myeloid cells [24]. Interestingly, recent research shows that M2 macrophages are heavily enriched in advanced RCC [22]. However, our study indicated that a positive clinicopathological feature group had a higher infiltration of M2 macrophages. This inconformity may be related to the different methods of measurement. The study above applied single-cell RNA sequencing, and our data used ordinary RNA sequencing. On the other hand, although the enrichment of M2 macrophages decreased in the advanced groups, we did not determine their functional state. Finally, studies found that CD8+ T cells and macrophages can interact to induce immune dysfunction [22], and we also did not obtain any data about whether they interacted and the intensity of this interaction. These reasons above may lead to different clinical phenotypes. In addition, specific markers of immune cells can be expressed in a proportion of cancer cells. This immune-like transcriptional reprogramming is considered to be an important hallmark of cancer [40]. Additionally, this immune mimicry may have an impact on the results of immune infiltration.
Our study also showed that regulatory T cells, resting mast cells, and resting dendritic cells are predictors for prognosis in RCC. Meanwhile, the expression of PLAUR, CAT, P2RY8, ERBB2, and TIMP1 is close to the infiltration of TIICs. This finding reveals that the five hub genes above may regulate tumor progression in an immunity-dependent pattern and can be potential biomarkers for immunotherapy. This research has some limitations. The main issue is a bioinformatics design without experimental verification. Additionally, the analysis of potential mechanisms is superficial, and further precise exploration should be carried out, such as the expression of hub genes at the post-transcriptional level, precising m6A sites of hub genes. Moreover, we did not obtain the data on which patients underwent immune therapy and did not explore the relationship between hub gene expression and immune response.

Conclusions
We identified six hub genes through multi-dimensional screening. They all possess strong predictive value for prognosis and clinicopathological features. Meanwhile, we found that m6A modification and immune infiltration may be involved in the development of RCC. Furthermore, hub genes may regulate the progression of RCC in an m6A-and immunity-dependent manner. Overall, our study provides new biomarkers for treatment and diagnosis.