METTL14 Acts as a Potential Regulator of Tumor Immune and Progression in Clear Cell Renal Cell Carcinoma

Clear cell renal cell carcinoma (ccRCC) is characterized by its insensitivity to chemoradiotherapy and lacks effective diagnostic and prognostic biomarkers. In this study, we focused on the role of m6A RNA methylation regulators for tumor immunity. Based on the expression of 20 m6A regulators, consensus clustering was performed to divide patients into cluster1/cluster2 and showed that there was a survival difference between the two clusters. Through cox regression analysis, five hub m6A regulators were screened to construct a risk model. Further analysis showed that the risk score was an independent prognostic factor. GSEA, GSVA, and KEGG analysis revealed that immune cell pathways played a critical role between the high risk group and low risk group. Combined with CIBERSORT and survival analysis, five hub tumor-infiltrating immune cells (TIICs) were identified for further study. Meanwhile, correlation analysis indicated that IGF2BP2 was positively associated with activated memory CD4 T cell and METTL14 was negatively correlated to the regulatory T cell. Therefore, IGF2BP2 and METTL14 were regarded as key genes. Further study verified that only METTL14 possessed good diagnostic and prognostic value. Then, GSEA exhibited that METTL14 was mainly enriched in chemokine related pathways. We also found that CCL5 was negatively correlated to METTL14 and might serve as a potential target of METTL14. In conclusion, these findings suggest that the METTL14/CCL5/Tregs axis is a potential signaling pathway for regulating tumor immunity, and might become novel therapeutic targets for ccRCC.

Introduction

Renal cell carcinoma (RCC) is one of the most common solid tumors in the urinary system. In the United States, the latest research infers that there will be 73,750 new cases and 14,830 people will die of this cancer in 2020 (Siegel et al., 2020). According to the survival data from the Surveillance, Epidemiology, and End Results (SEER) database, the 5-year relative survival is 75.2%¹. Clear cell renal cell carcinoma (ccRCC) is the most frequent subtype of RCC, and accounts for 70–80% (Moch et al., 2016) of cases. Due to the insensitivity of chemoradiotherapy and hidden clinical manifestations, the mortality of ccRCC is persistently high and surgery is still the primary treatment (Ljungberg et al., 2015; Capitanio and Montorsi, 2016). Moreover, approximately 30% of patients have distant metastasis at the first time of diagnosis, which leads to the loss of the chance for surgery (Gong et al., 2016). It is, therefore, necessary to discover novel diagnostic and prognostic biomarkers and therapeutic strategies for ccRCC.

At present, many studies have verified that post-translational RNA modifications play a vital role in various diseases (Prieto et al., 2020), especially cancers (Huang et al., 2020). Among these modifications, m6A RNA methylation has been considered as the most pervasive modification in regulating RNA transcription (Barbieri et al., 2017), RNA splicing (Zhao et al., 2014), RNA degradation (Du et al., 2016), etc. It has been reported that m6A methylation was performed via a methyltransferase complex composed of “writer,” “eraser,” and “reader” (Chen et al., 2019; Lan et al., 2019). Methyltransferase like 14 (METTL14), a member of “writer,” are capable of catalyzing the formation of m6A. Previous studies have demonstrated that METTL14 was closely associated with tumorigenesis and progression. In bladder cancer, Gu et al. revealed that METTL14 could inhibit tumor initiating cells (TIC) self-renewal and tumorigenesis (Gu et al., 2019). It has been reported that METTL14 acts as a tumor suppressor for inhibiting tumor metastasis via regulating miRNA in hepatocellular carcinoma (Ma et al., 2016). Similarly, Yang et al. also verified that METTL14 suppressed the proliferation and metastasis of colorectal cancer through mediating lncRNA XIST (Yang et al., 2020). In breast cancer, METTL14 serves as an oncogene and drives tumor cell migration and invasion (Yi et al., 2020). However, the study of METTL14 is lacking and the molecular mechanisms are unclear in ccRCC.

In this study, we focused on the functions of m6A RNA methylation regulators and selected METTL14 as a key gene for in-deep study in ccRCC. Through comprehensive bioinformatics analysis, the correlation and potential mechanisms between METTL14 and tumor immune were further clarified.

Materials and Methods

Data Collection and Study Design

All gene expression data and clinical traits data were obtained from The Cancer Genome Atlas (TCGA) database². The study design is shown in detail in the flow diagram Supplementary Figure 1.

Identification of Differentially Expressed Genes (DEGs) or m6A Regulators (DE-m6A)

The “DESeq2” (Love et al., 2014) and “edgeR” (Robinson et al., 2009) packages were used to screen DEGs or DE-m6A with adjusted p-value <0.05 and | log2(FoldChange)| > 2. The wilcoxon test was utilized to identify differentially expressed m6A regulators according to a cutoff criterion of adjusted p-value <0.05.

Correlation Analysis

The “cor.mtest” algorithm and “corrplot” package³ were used to perform spearman correlation analysis and visualize the results. P-value <0.05 was selected as a cutoff point.

Construction of Protein-Protein Interaction (PPI) Network

The m6A RNA methylation regulators were submitted to the STRING database⁴ for constructing the PPI network. Moreover, cytoscape (Shannon et al., 2003) software was used to visualize and analyze the network.

Consensus Clustering of m6A Regulators

The “ConsensusClusterPlus” package (Wilkerson and Hayes, 2010) was utilized to divide ccRCC patients into different clusters (k = 2–9). Moreover, we performed survival analysis to evaluate the prognostic value of different clusters. P-value <0.05 was considered statistically significant.

Principal Component Analysis (PCA)

In this study, PCA was performed by using “limma” package (Ritchie et al., 2015) and visualized via using “ggplot2” package⁵ in R program.

Construction of Cox Proportional Hazard Regression Model

The “survival” package⁶ was applied to perform a univariate cox regression analysis for preliminarily assessing the prognostic value of m6A regulators. Least Absolute Shrinkage and Selection Operator (LASSO) regression analysis was then used to further screen key prognostic factors and build a risk model via running “glmnet” (Friedman et al., 2010) and “survival” packages in R. The risk score (RS) formula as following:

$RS=\sum _{i=1}^{n}Coef(i)X(i)$

The Coef (i) represents the coefficient and the X(i) represents the expression of selected genes. Moreover, univariate and multivariate cox regression analyses were used to assess the prognostic value of the risk model.

Survival Analysis and Receiver Operator Characteristic (ROC) Curve Analysis

According to the median of gene expression/risk score/TIICs abundance, ccRCC patients were divided into two groups (high group and low group). The “survival” package and GraphPad Prism software were utilized to draw the survival curves. A log rank p-value <0.05 was considered statistically significant. We performed further ROC curves to evaluate the diagnostic value.

Gene Set Enrichment Analysis (GSEA) and Gene Set Variation Analysis (GSVA)

All ccRCC patients were divided into high group and low group. Then, GSEA online tool (Subramanian et al., 2005) was used to find potential molecular mechanisms. In addition, the “GSVA” package (Hänzelmann et al., 2013) was applied to screen differential signaling pathways between the high risk group and the low risk group.

Kyoto Encyclopedia of Genes and Genomes (KEGG) Pathway Enrichment Analysis

To further analyze the biological function of DEGs, KEGG enrichment analysis was performed via running the “clusterProfiler” package (Yu et al., 2012) in the R program. P-value <0.05 was selected as the cutoff point.

Purity Analysis of Stromal Cell and Immune Cell

The “estimate” package (Yoshihara et al., 2013) was applied to identify the purity of stromal cells and immune cells. Then, we performed a survival analysis to evaluate the prognostic value.

Distribution of Tumor-Infiltrating Immune Cells (TIICs)

CIBERSORT algorithm⁷ with LM22 signature was utilized to calculate the relative abundance of 22 TIICs at 1,000 permutations. P-value <0.05 was considered statistically significant.

Genetic Alteration Analysis and Immunohistochemistry (IHC) Analysis

The “cBioPortal” online tool (Cerami et al., 2012; Gao et al., 2013) was used to analyze the genetic alteration of hub m6A regulators. Furthermore, the IHC results of METTL14 were obtained from The Human Protein Atlas (HPA) database.

Statistical Analysis

In this study, SPSS 22.0 and GraphPad Prism 7.0 software were applied to process data. All data were represented as mean ± SD. Student’s t-test and Chi-square test were used to analyze the data. P-value <0.05 was considered statistically significant.

Results

The Expression Level and Interaction of m6A RNA Methylation Regulators in ccRCC

Gene expression data were obtained from the TCGA KIRC dataset. As shown in Figures 1A,B, a total of 15 differentially expressed m6A regulators (DE-m6A) were identified. There were six up-regulated m6A regulators (YTHDF3, IGF2BP2, METTL14, YTHDF2, HNRNPA2B1, and ZC3H13) and nine down-regulated m6A regulators (IGF2BP3, KIAA1429, RBM15, FTO, YTHDC2, ALKBH3, ALKBH5, METTL3, and WTAP). Moreover, correlation analysis indicated that METTL14 was likely to be the hub gene and strongly positively correlated to YTHDF3 (R = 0.60), KIAA1429 (R = 0.53), ZC3H13 (R = 0.52), and YTHDC1 (R = 0.66) (Figure 1C). Further PPI network analysis exhibited that all m6A regulators had close interactions with each other, especially METTL14 (Figure 1D).

FIGURE 1

Figure 1. The expression level and interaction of twenty m6A RNA methylation regulators in ccRCC. (A,B) The expression level of m6A regulators between ccRCC tissues and normal tissues. (C) Spearman correlation among the 20 m6A regulators. (D) PPI network of m6A regulators following STRING database. ccRCC, clear cell renal cell carcinoma; PPI, protein-protein interaction. *P < 0.05; **P < 0.01; ***P < 0.001.

Consensus Clustering of m6A RNA Methylation Regulators

Based on the expression similarity of m6A regulators, k = 2 was selected as the threshold for dividing the patients into two subgroups (cluster1 and cluster2) in the TCGA KIRC dataset (Figures 2A,B). PCA further verified that there was an obvious difference in transcription profile between cluster1 and cluster2 (Figure 2C). Moreover, we found that patients in the cluster2 subgroup had poor overall survival (OS) (Figure 2D). Cox regression analysis showed that the cluster was an independent prognostic factor (Table 1). Combined with clinical characteristics, results revealed that M stage, gender, and survival state had a significant distinction between the above two subgroups (Figure 2E).

FIGURE 2

Figure 2. Consensus clustering of m6A regulators. (A) Cumulative distributive function (CDF) for k = 2–9. (B) ccRCC patients were divided into two clusters for k = 2. (C) PCA for total RNA expression pattern. Red represents cluster1 and green represents cluster2. (D) Patients within cluster1 had better OS. (E) The distinction of clinical traits between cluster1 and cluster2. ccRCC, clear cell renal cell carcinoma; PCA, principal component analysis; OS, overall survival. *P < 0.05; **P < 0.01; ***P < 0.001.

TABLE 1

Table 1. The cox regression analysis of two clusters.

Cox Regression Analysis and Construction of Risk Model

To evaluate the prognostic value of m6A regulators, univariate and LASSO cox regression analyses were performed. As shown in Figure 3A, a total of fourteen m6A regulators were regarded as prognostic factors (p < 0.05). In LASSO cox regression analysis, five hub m6A regulators (IGF2BP3, IGF2BP2, METTL14, METTL3, and HNRNPA2B1) were screened for the construction of a risk model (Figures 3B,C). Furthermore, the associated coefficients were obtained from LASSO analysis. The risk score = 0.147^∗IGF2BP3 expression – 0.247^∗METTL14 expression + 0.023^∗IGF2BP2 expression + 0.006^∗ HNRNPA2B1 expression + 0.030^∗METTL3 expression. According to the risk score, the ccRCC patients were divided into a high risk group and a low risk group. Survival analysis indicated that patients with a high risk score had shorter overall survival times (Figure 3D). In addition, the ROC curve exhibited that the risk score had a good diagnostic value (Figure 3E, AUC = 0.724). Further univariate (Figure 3F, HR = 2.355, p < 0.001) and multivariate (Figure 3G) cox regression analysis demonstrated that risk score was a potential independent prognostic factor in ccRCC. Moreover, we also found that the risk score was significantly correlated to the T stage, N stage, M stage, TNM stage, grade, gender, and survival state (Figure 3H).

FIGURE 3

Figure 3. Construction of risk model and assessment of its diagnostic and prognostic value. (A) Univariate cox regression analysis of twenty m6A regulators. (B,C) Five genes (IGF2BP3, IGF2BP2, METTL14, METTL3, and HNRNPA2B1) were selected as the key variate for the risk model via LASSO cox regression. (D) Patients with high risk have a poor OS. (E) ROC curve showed that the risk model had a good diagnostic value (AUC = 0.724). (F,G) Univariate and multivariate cox regression verified that risk score was an independent prognostic factor. (H) The distinction of pathological characteristics between high risk group and low risk group. LASSO, least absolute shrinkage and selection operator; OS, overall survival; ROC, receiver operating characteristic; AUC, area under the curve. *P < 0.05; **P < 0.01; ***P < 0.001.

The Risk Model Was Closely Associated With Immune-Related Pathways

To investigate potential molecular mechanisms of the risk model, GSEA, GSVA, and KEGG pathway analyses were performed. GSEA indicated that the high risk group was mainly enriched in immune cell-related pathways (Figures 4A,B). Further GSVA also verified that there were many differentially expressed immune cell-related pathways between the high risk group and the low risk group (Figures 4C,D). In addition, we screened differentially expressed genes (DEGs) via DESeq2 and edgeR algorithms. As shown in Figures 4E–G, 162 DEGs were identified. KEGG pathway enrichment analysis proved that the DEGs primarily enriched in tumor-immune related signaling pathways including TNF signaling pathway, JAK-STAT signaling pathway, and IL-17 signaling pathway (Figure 4H). These findings suggested that the risk model consisted of hub m6A regulators and played a vital role through mediating tumor immunity.

FIGURE 4

Figure 4. Risk signature was closely associated with immune-related pathways. (A,B) GSEA indicated that the high risk group was mainly enriched in immune cell related pathways. (C,D) GSVA showed that there were many differentially expressed immune cell-related pathways between the high risk group and the low risk group. (E–G) Identification of DEGs between high risk group and low risk group via “DESeq2” and “edgeR” packages. (H) KEGG pathway enrichment analysis of DEGs between high risk group and low risk group. GSEA, gene set enrichment analysis; GSVA, gene set variation analysis; DEGs, differentially expressed genes; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Prognostic Role of TIICs and Correlation With Hub m6A Regulators

According to the previous analysis, we further focused on the role of TIICs in ccRCC. ESTIMATE was performed to calculate the purity of immune cells and stromal cells. Survival analysis indicated that a high ImmuneScore predicted poor OS (Figures 5A–C). Therefore, CIBERSORT was utilized to analyze the distribution of 22 TIICs. As shown in Figure 5D, there were twelve differential TIICs between ccRCC samples and normal samples. Further survival analysis revealed that high infiltration of resting dendritic cells/resting mast cells and low infiltration of activated CD4 memory T cells/follicular helper T cells/regulatory T cells predicted shorter OS in ccRCC (Figures 5E–I), which were regarded as hub TIICs. Correlation analysis exhibited that IGF2BP2 was significantly positively associated with activated CD4 memory T cells (R = 0.46) and METTL14 was significantly negatively associated with regulatory T cells (Tregs) (R = −0.46) (Figure 5J). Thus, we inferred that IGF2BP2 might promote the infiltration of activated CD4 memory T cells and METTL14 might decrease the infiltration of Tregs in ccRCC tissues. Moreover, we found that the five hub m6A regulators all exhibited genetic alterations in ccRCC (Figure 5K). Taken together, IGF2BP2 and METTL14 were selected for subsequent analysis.

FIGURE 5

Figure 5. The prognostic role of TIICs and the correlation between TIICs and five hub m6A regulators. (A–C) ESTIMATE analysis indicated that only the immune score had good prognostic value and patients with a high immune score had poor OS. (D) Identification of differentially infiltrated TIICs between ccRCC tissues and normal tissues via CIBERSORT algorithm. Further survival analysis revealed that high infiltration of resting dendritic cells (E)/resting mast cells (F) predicted better OS and high infiltration of activated memory CD4 T cells (G)/helper follicular T cells (H)/regulatory T cells (I) predicted poor OS. (J) The correlation between hub m6A regulators and hub TIICs. (K) The genetic alteration of hub m6A regulators. TIICs, tumor-infiltrating immune cells; ESTIMATE, Estimation of STromal and Immune cells in Malignant Tumor tissues using Expression data; OS, overall survival; ccRCC, clear cell renal cell carcinoma.

METTL14 Was Down-Regulated and Closely Associated With Clinical Traits in ccRCC

Survival analysis indicated that the high expression of METTL14 predicted better OS (Figure 6A) and disease-free survival (DFS) (Figure 6B). Patients with high expression of IGF2BP2 had shorter OS (Figure 6C), but DFS analysis of IGF2BP2 had no statistical difference (Figure 6D). Therefore, METTL14 was selected as the key gene for further study. In the TCGA KIRC dataset, the expression of METTL14 decreased in ccRCC tissues (Figures 6E,F). As shown in Figures 6G,H, the expression of METTL14 was negatively correlated with the TNM stage and Grade. The ROC curve demonstrated that the expression of METTL14 could effectively differentiate ccRCC tissues from normal tissues (Figures 6I,J). In the HPA database, we also found that the protein expression level of METTL14 was down-regulated in ccRCC tissues (Figure 6K).

FIGURE 6

Figure 6. The METTL14 expression level was closely associated with various clinical traits in ccRCC. (A,B) High expression of METTL14 predicted better OS and DFS. (C,D) High expression of IGF2BP2 predicted poor OS and did not affect DFS. (E,F) METTL14 was down-regulated in ccRCC samples. (G,H) The expression level of METTL14 was significantly negatively correlated to TNM stage and Grade. (K) IHC result from the HPA database verified that the expression of METTL14 was decreasing in ccRCC. ccRCC, clear cell renal cell carcinoma; OS, overall survival; DFS, disease-free survival; TNM, Tumor Node Metastasis; IHC, immunohistochemistry; HPA, The Human Protein Atlas. (I,J) Expression of METTL14 differentiates the adjacent normal tissues and ccRCC cancer tissues. *P < 0.05; **P < 0.01; ***P < 0.001.

METTL14 Might Inhibit the Infiltration of Tregs via Suppressing CCL5

GSEA was performed to analyze the specific molecular mechanisms of METTL14 in ccRCC. The results indicated that the low expression of METTL14 was mainly enriched in chemokine related pathways (Figure 7A). We then analyzed the correlation between METTL14 and members of the CCL/CXCL superfamily. According to the threshold of correlation, CCL5 (Figure 7B, Spearman = −0.43, Pearson = −0.41) and CCL17 (Figure 7C, Spearman = −0.34, Pearson = −0.36) were regarded as the potential downstream targets of METTL14 (Table 2). Survival analysis showed that only CCL5 had a good prognostic value and high CCL5 expression predicted poor OS (Figures 7D,E). These results suggest that CCL5 might play a more important role than CCL17 in ccRCC. Moreover, CCL5 was up-regulated in ccRCC (Figure 7F) and positively correlated to FOXP3 (Figure 7G). Therefore, we speculate that the METTL14/CCL5/Tregs axis might be a potential pivotal regulated pathway of tumor immunity and progression in ccRCC (Figure 7H), a subject worthy of further in-depth study.

FIGURE 7

Figure 7. CCL5 was a potential downstream targeted molecule of METTL14. (A) GSEA indicated that low expression of METTL14 was enriched in chemokine-related pathways. Further analysis uncovered that METTL14 was significantly negatively associated with CCL5 (B) and CCL17 (C). (D) Survival analysis exhibited that patients with high expression of CCL5 had poor OS. (E) The expression of CCL17 did not affect OS. (F) CCL5 was elevated in ccRCC tissues. (G) The expression of CCL5 was markedly positively correlated to FOXP3 in ccRCC. (H) The schematic representation of the potential molecular mechanism of METTL14/CCL5/Tregs axis in ccRCC. GSEA, gene set enrichment analysis; OS, overall survival; ccRCC, clear cell renal cell carcinoma. *P < 0.05; **P < 0.01; ***P < 0.001.

TABLE 2

Table 2. The correlation between METTL14 and chemokines in ccRCC.

Discussion

ccRCC is a common lethal tumor that is characterized as involving high morbidity and mortality. For advanced and metastatic ccRCC, molecularly targeted therapy (Cowey and Hutson, 2010) was selected as the first-line treatment, which improves the prognosis of patients. In recent years, immunotherapy has become the latest research hotspot. A lot of immunotherapeutic agents that target cytokine or checkpoint pathways are designed to enhance the response of the immune system in attacking tumor cells (Carosella et al., 2015; Deleuze et al., 2020). However, there are still some patients with medicine resistance and poor survival. Therefore, it is necessary to further clarify the molecular mechanisms of tumor immunity and find novel effective therapeutic targets in ccRCC.

This study analyzed the potential diagnostic and prognostic value of m6A regulators. Combined with TIICs analysis, we found that the METTL14/CCL5/Tregs axis was a potential tumor immune regulative pathway in ccRCC.

Recently, increasing studies have reported that m6A modification of RNA, which is mediated by three type regulators (“writer,” “eraser,” and “reader”), and is involved in central events in biology (Meyer and Jaffrey, 2014; Liu and Pan, 2016). METTL14 is a member of “writer” and identified as a catalytic component of the methyltransferase complex responsible for creating m6A modifications in various biological processes (Huang et al., 2019). Animal experiments verified that METTL14 could regulate striatal function in a mouse model (Koranda et al., 2018). In the hematological system, METTL14 suppressed hematopoietic stem/progenitor differentiation and promoted leukemogenesis via mediating MYB/MYC (Weng et al., 2017). Moreover, Chen et al. revealed that METTL14 inhibited colorectal cancer growth and metastasis through targeting miR-375. Our study uncovered that METTL14 was down-regulated in ccRCC and significantly negatively correlated to tumor stage. Survival analysis also demonstrated that the high expression of METTL14 predicted better OS and DFS.

We still lack studies of the correlation between METTL14 and tumor immunity. In this study, we indicated that METTL14 might act as a modulator of Tregs to regulate tumor immunity in ccRCC. Previous studies have reported that immunosuppressive cells including regulatory T (Treg) (Sakaguchi, 2000) cells, myeloid derived suppressor cells (MDSCs) (Gabrilovich, 2017), and tumor-associated macrophages (TAMs) (Mantovani et al., 2017) inhibit effective anti-tumor immune responses. According to molecular markers, Tregs fall into two major categories: resting/naïve Treg (rTreg, FOXP3^loCD45RA⁺CD25^lo) and effector/activated Treg (eTreg, FOXP3^hiCD45RA^–CD25^hi) (Ohue and Nishikawa, 2019). Effector/activated Treg has a strong immunosuppressive function and inhibits anti-tumor immunity through various molecular mechanisms, such as the regulation of the immune checkpoint molecules (Kamada et al., 2019), induction of apoptosis for immune effector cells (Gotot et al., 2018), and production of immunosuppressive cytokines (Collison et al., 2007). It has been reported that inhibition of Nr4a suppresses tumor progression via weakening Treg-mediated immune tolerance (Hibino et al., 2018). Our study also demonstrated that infiltration of Tregs was a potential risk factor in ccRCC.

Previous studies have also revealed that Tregs accumulate in the tumor microenvironment (TME) through chemokine gradients (Bromley et al., 2008). CCR4-CCL17, CCR8-CCL1, and CXCR3-CXCL9/10/11 are the major signaling pathways for Tregs chemotaxis (Ohue and Nishikawa, 2019). Wang et al. verified that tumor-derived FOXP3 recruited FOXP3(+)Treg cells by activating CCL5 in pancreatic ductal adenocarcinoma (Wang et al., 2016). Analogously, You et al. also revealed that ovarian cancer stem cells recruited Tregs via the CCL5-CCR5 axis (You et al., 2017). In this study, GSEA and correlation analysis indicated that CCL5 was negatively associated with METTL14, which inferred that CCL5 was a potential target for METTL14. CCL5 (C-C Motif Chemokine Ligand 5), also named RANTES, can activate several chemokine receptors including CCR1, CCR3, CCR4, and CCR5 for chemotaxis of immune cells. Our research also indicated that the expression of CCL5 was positively correlated to FOXP3 in ccRCC tissue.

Conclusion

In conclusion, METTL14 was identified as a key m6A RNA methylation regulator via integrated bioinformatics analysis. Combining this approach with tumor immune analysis, we found that METTL14 might suppress tumor progression through mediating Tregs. Further study indicates that CCL5 served as a potential downstream molecule of METTL14 to regulate Tregs chemotaxis. Taken together, these findings suggest that the METTL14/CCL5/Tregs axis is expected to be a novel therapeutic target for ccRCC. However, the specific mechanisms require further research.

Data Availability Statement

Publicly available datasets were analyzed in this study. This data can be found here: https://www.cancer.gov/tcga.

Ethics Statement

We certify that this manuscript is original and has not been published and will not be submitted elsewhere for publication. No data have been fabricated or manipulated (including images) to support our conclusions.

Author Contributions

XZ and KC designed this study. TX and SG performed data collection and analysis, wrote the manuscript, and contributed to preparing and making figures and tables. TX, SG, HR, and JL performed the majority of the experiments. YL, DL, and JT collected the clinical samples and managed the clinical data. JT and JS reviewed the relevant literature. HY, XZ, and KC provided conceptual advice and critically reviewed the manuscript. All authors read and approved the final manuscript.

Funding

This study was supported by the Key Research and Development Plan in China (grant no. 2017YFB1303100), the National Natural Science Foundation of China (grant nos. 81672524, 81672528, and 81874090), the Hubei Provincial Natural Science Foundation of China (grant no. 2018CFA038), the Independent Innovation Foundation of Huazhong University of Science and Technology (grant no. 118530309), and the Clinical Research Physician Program of Tongji Medical College, Huazhong University of Science and Technology (grant no. 5001530015).

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.

Acknowledgments

We would like to express our gratitude to all the members who participated in the discussion and assisted in this study.

Supplementary Material

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

Supplementary Figure 1 | Flow chart of study design. The processes of data analysis are shown in the flow chart. TCGA, The Cancer Genome Atlas; KIRC, Kidney Clear Cell Carcinoma; DE-m6A, differentially expressed m6A; ESTIMATE, Estimation of STromal and Immune cells in Malignant Tumor tissues using Expression data; TIICs, tumor-infiltrating immune cells; LASSO, least absolute shrinkage and selection operator; GSEA, gene set enrichment analysis; GSVA, gene set variation analysis; OS, overall survival; DFS, disease-free survival; HPA, The Human Protein Atlas; ROC, receiver operating characteristic.

Footnotes

1. ^ https://seer.cancer.gov/
2. ^ https://www.cancer.gov/tcga
3. ^ https://CRAN.R-project.org/package=corrplot
4. ^ https://string-db.org
5. ^ https://CRAN.R-project.org/package=ggplot2
6. ^ https://CRAN.R-project.org/package=survival
7. ^ https://cibersort.stanford.edu/

References

Barbieri, I., Tzelepis, K., Pandolfini, L., Shi, J., Millán-Zambrano, G., Robson, S. C., et al. (2017). Promoter-bound METTL3 maintains myeloid leukaemia by m(6)a-dependent translation control. Nature 552, 126–131. doi: 10.1038/nature24678

PubMed Abstract | CrossRef Full Text | Google Scholar

Bromley, S. K., Mempel, T. R., and Luster, A. D. (2008). Orchestrating the orchestrators: chemokines in control of T cell traffic. Nat Nat. Immunol. 9, 970–980. doi: 10.1038/ni.f.213

PubMed Abstract | CrossRef Full Text | Google Scholar

Capitanio, U., and Montorsi, F. (2016). Renal cancer. Lancet. 387, 894–906. doi: 10.1016/s0140-6736(15)00046-x

CrossRef Full Text | Google Scholar

Carosella, E. D., Ploussard, G., LeMaoult, J., and Desgrandchamps, F. (2015). A Systematic Review of Immunotherapy in Urologic Cancer: Evolving Roles for Targeting of CTLA-4, PD-1/PD-L1, and HLA-G. Eur. Urol. 68, 267–279. doi: 10.1016/j.eururo.2015.02.032

PubMed Abstract | CrossRef Full Text | Google Scholar

Cerami, E., Gao, J., Dogrusoz, U., Gross, B. E., Sumer, S. O., Aksoy, B. A., et al. (2012). The cBio cancer genomics portal: an open platform for exploring multidimensional cancer genomics data. Cancer Dis.covery. 2, 401–404. doi: 10.1158/2159-8290.Cd-12-0095

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, X. Y., Zhang, J., and Zhu, J. S. (2019). The role of m(6)A RNA methylation in human cancer. Mol. Cancer. 18:103. doi: 10.1186/s12943-019-1033-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Collison, L. W., Workman, C. J., Kuo, T. T., Boyd, K., Wang, Y., Vignali, K. M., et al. (2007). The inhibitory cytokine IL-35 contributes to regulatory T-cell function. Nature. 450, 566–569. doi: 10.1038/nature06306

PubMed Abstract | CrossRef Full Text | Google Scholar

Cowey, C. L., and Hutson, T. E. (2010). Molecularly targeted agents for renal cell carcinoma: the next generation. Clin. Adv. Hematol. Oncol. 8, 357–360,361–374.

Google Scholar

Deleuze, A., Saout, J., Dugay, F., Peyronnet, B., Mathieu, R., Verhoest, G., et al. (2020). Immunotherapy in renal cell carcinoma: the future is now. Int. J. Mol. Sci. 21:2532. doi: 10.3390/ijms21072532

PubMed Abstract | CrossRef Full Text | Google Scholar

Du, H., Zhao, Y., He, J., Zhang, Y., Xi, H., Liu, M., et al. (2016). YTHDF2 destabilizes m(6)A-containing RNA through direct recruitment of the CCR4-NOT deadenylase complex. Nat. Commun. 7:12626. doi: 10.1038/ncomms12626

PubMed Abstract | CrossRef Full Text | Google Scholar

Friedman, J., Hastie, T., and Tibshirani, R. (2010). Regularization paths for generalized linear models via coordinate descent. J. Stat. Softw. 33, 1–22.

Google Scholar

Gabrilovich, D. I. (2017). Myeloid-derived suppressor cells. Cancer Immunol. Res. 5, 3–8. doi: 10.1158/2326-6066.Cir-16-0297

PubMed Abstract | CrossRef Full Text | Google Scholar

Gao, J., Aksoy, B. A., Dogrusoz, U., Dresdner, G., Gross, B., Sumer, S. O., et al. (2013). Integrative analysis of complex cancer genomics and clinical profiles using the cBioPortal. Sci. Signal. 6:l1. doi: 10.1126/scisignal.2004088

PubMed Abstract | CrossRef Full Text | Google Scholar

Gong, J., Maia, M. C., Dizman, N., Govindarajan, A., and Pal, S. K. (2016). Metastasis in renal cell carcinoma: biology and implications for therapy. Asian J. Urol. 3, 286–292. doi: 10.1016/j.ajur.2016.08.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Gotot, J., Dhana, E., Yagita, H., Kaiser, R., Ludwig-Portugall, I., and Kurts, C. (2018). Antigen-specific helios-, neuropilin-1-tregs induce apoptosis of autoreactive B cells via PD-L1. Immunol. Cell. Biol. 96, 852–862. doi: 10.1111/imcb.12053

PubMed Abstract | CrossRef Full Text | Google Scholar

Gu, C., Wang, Z., Zhou, N., Li, G., Kou, Y., Luo, Y., et al. (2019). Mettl14 inhibits bladder TIC self-renewal and bladder tumorigenesis through N6-methyladenosine of Notch1. Mol. Cancer. 18:168. doi: 10.1186/s12943-019-1084-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Hänzelmann, S., Castelo, R., and Guinney, J. (2013). GSVA: gene set variation analysis for microarray and RNA-Seq data. BMC Bioinf.ormatics. 14:7. doi: 10.1186/1471-2105-14-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Hibino, S., Chikuma, S., Kondo, T., Ito, M., Nakatsukasa, H., Omata-Mise, S., et al. (2018). Inhibition of Nr4a receptors enhances anti-tumor immunity by breaking Treg-mediated immune tolerance. Cancer Res. 78, 3027–3040. doi: 10.1158/0008-5472.Can-17-3102

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, A. Z., Delaidelli, A., and Sorensen, P. H. (2020). RNA modifications in brain tumorigenesis. Acta Neuropathol. Commun. 8:64. doi: 10.1186/s40478-020-00941-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, J., Dong, X., Gong, Z., Qin, L.-Y., Yang, S., Zhu, Y.-L., et al. (2019). Solution structure of the RNA recognition domain of METTL3-METTL14 N6-methyladenosine methyltransferase. Protein Cell 10, 272–284. doi: 10.1007/s13238-018-0518-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Kamada, T., Togashi, Y., Tay, C., Ha, D., Sasaki, A., Nakamura, Y., et al. (2019). PD-1(+) regulatory T cells amplified by PD-1 blockade promote hyperprogression of cancer. Proc. Natl. Acad. Sci. U.S.A. 116, 9999–10008. doi: 10.1073/pnas.1822001116

PubMed Abstract | CrossRef Full Text | Google Scholar

Koranda, J. L., Dore, L., Shi, H., Patel, M. J., Vaasjo, L. O., Rao, M. N., et al. (2018). Mettl14 is essential for epitranscriptomic regulation of striatal function and learning. Neuron 99, 283–92e5. doi: 10.1016/j.neuron.2018.06.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Lan, Q., Liu, P. Y., Haase, J., Bell, J. L., Huttelmaier, S., and Liu, T. (2019). The critical role of RNA m(6)a methylation in cancer. Cancer Res. 79, 1285–1292. doi: 10.1158/0008-5472.CAN-18-2965

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, N., and Pan, T. (2016). N6-methyladenosine–encoded epitranscriptomics. Nat. Struct. Mol. Biol. 23, 98–102. doi: 10.1038/nsmb.3162

PubMed Abstract | CrossRef Full Text | Google Scholar

Ljungberg, B., Bensalah, K., Canfield, S., Dabestani, S., Hofmann, F., Hora, M., et al. (2015). EAU guidelines on renal cell carcinoma: 2014 update. Eur. Urol. 67, 913–924. doi: 10.1016/j.eururo.2015.01.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Love, M. I., Huber, W., and Anders, S. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biology 15:550. doi: 10.1186/s13059-014-0550-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Ma, J. Z., Yang, F., Zhou, C. C., Liu, F., Yuan, J. H., Wang, F., et al. (2016). METTL14 suppresses the metastatic potential of HCC by modulating m(6) a-dependent primary miRNA processing. Hepatology 65, 529–543. doi: 10.1002/hep.28885

PubMed Abstract | CrossRef Full Text | Google Scholar

Mantovani, A., Marchesi, F., Malesci, A., Laghi, L., and Allavena, P. (2017). Tumour-associated macrophages as treatment targets in oncology. Nat. Rev. Clin. Oncol. 14, 399–416. doi: 10.1038/nrclinonc.2016.217

PubMed Abstract | CrossRef Full Text | Google Scholar

Meyer, K. D., and Jaffrey, S. R. (2014). The dynamic epitranscriptome: N6-methyladenosine and gene expression control. Nat. Rev. Mol. Cell Biol. 15, 313–326. doi: 10.1038/nrm3785

PubMed Abstract | CrossRef Full Text | Google Scholar

Moch, H., Cubilla, A. L., Humphrey, P. A., Reuter, V. E., and Ulbright, T. M. (2016). The 2016 who classification of tumours of the urinary system and male genital organs-part a: renal. Penile Testicular Tumours. Eur. Urol. 70, 93–105. doi: 10.1016/j.eururo.2016.02.029

PubMed Abstract | CrossRef Full Text | Google Scholar

Ohue, Y., and Nishikawa, H. (2019). Regulatory T (Treg) cells in cancer: can treg cells be a new therapeutic target? Cancer Sci. 110, 2080–2089. doi: 10.1111/cas.14069

PubMed Abstract | CrossRef Full Text | Google Scholar

Prieto, M., Folci, A., and Martin, S. (2020). Post-translational modifications of the fragile X mental retardation protein in neuronal function and dysfunction. Mol. Psychiatry. 25, 1688–1703. doi: 10.1038/s41380-019-0629-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Ritchie, M. E., Phipson, B., Wu, D., Hu, Y., Law, C. W., Shi, W., et al. (2015). limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 43:e47. doi: 10.1093/nar/gkv007

PubMed Abstract | CrossRef Full Text | Google Scholar

Robinson, M. D., McCarthy, D. J., and Smyth, G. K. (2009). edgeR: a bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 26, 139–140. doi: 10.1093/bioinformatics/btp616

PubMed Abstract | CrossRef Full Text | Google Scholar

Sakaguchi, S. (2000). Regulatory T cells: key controllers of immunologic self-tolerance. Cell 101, 455–458. doi: 10.1016/s0092-8674(00)80856-9

CrossRef Full Text | Google Scholar

Shannon, P., Markiel, A., Ozier, O., Baliga, N. S., Wang, J. T., Ramage, D., et al. (2003). Cytoscape: a software environment for integrated models of biomolecular interaction networks. Geno. Res. 13, 2498–2504. doi: 10.1101/gr.1239303

PubMed Abstract | CrossRef Full Text | Google Scholar

Siegel, R. L., Miller, K. D., and Jemal, A. (2020). Cancer statistics, 2020. CA Cancer J. Clin. 70, 7–30. doi: 10.3322/caac.21590

PubMed Abstract | CrossRef Full Text | Google Scholar

Subramanian, A., Tamayo, P., Mootha, V. K., Mukherjee, S., Ebert, B. L., Gillette, M. A., et al. (2005). Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc. Natl. Acad. Sci. U.S.A. 102, 15545–15550. doi: 10.1073/pnas.0506580102

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, X., Lang, M., Zhao, T., Feng, X., Zheng, C., Huang, C., et al. (2016). Cancer-FOXP3 directly activated CCL5 to recruit FOXP3(+)Treg cells in pancreatic ductal adenocarcinoma. Oncogene 36, 3048–3058. doi: 10.1038/onc.2016.458

PubMed Abstract | CrossRef Full Text | Google Scholar

Weng, H., Huang, H., Wu, H., Qin, X., Zhao, B. S., Dong, L., et al. (2017). METTL14 inhibits hematopoietic stem/progenitor differentiation and promotes leukemogenesis via mRNA m6A modification. Cell Stem Cell 22, 191–205.e9. doi: 10.1016/j.stem.2017.11.016

PubMed Abstract | CrossRef Full Text | Google Scholar

Wilkerson, M. D., and Hayes, D. N. (2010). Consensus cluster plus: a class discovery tool with confidence assessments and item tracking. Bioinformatics 26, 1572–1573. doi: 10.1093/bioinformatics/btq170

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, X., Zhang, S., He, C., Xue, P., Zhang, L., He, Z., et al. (2020). METTL14 suppresses proliferation and metastasis of colorectal cancer by down-regulating oncogenic long non-coding RNA XIST. Mol. Cancer. 19:46. doi: 10.1186/s12943-020-1146-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Yi, D., Wang, R., Shi, X., Xu, L., Yilihamu, Y., and Sang, J. (2020). METTL14 promotes the migration and invasion of breast cancer cells by modulating N6-methyladenosine and hsa-miR-146a-5p expression. Oncol. Rep. 43, 1375–1386. doi: 10.3892/or.2020.7515

PubMed Abstract | CrossRef Full Text | Google Scholar

Yoshihara, K., Shahmoradgoli, M., Martínez, E., Vegesna, R., Kim, H., Torres-Garcia, W., et al. (2013). Inferring tumour purity and stromal and immune cell admixture from expression data. Nat. Commun. 4:2612. doi: 10.1038/ncomms3612

PubMed Abstract | CrossRef Full Text | Google Scholar

You, Y., Li, Y., Li, M., Mi, L., Wu, M., Qu, Y., et al. (2017). Ovarian cancer stem cells promote tumour immune privilege and invasion via CCL5 and regulatory T cells. Clin. Exp. Immunol. 191, 60–73. doi: 10.1111/cei.13044

PubMed Abstract | CrossRef Full Text | Google Scholar

Yu, G., Wang, L. G., Han, Y., and He, Q. Y. (2012). Cluster profiler: an R package for comparing biological themes among gene clusters. Omics 16, 284–287. doi: 10.1089/omi.2011.0118

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, X., Yang, Y., Sun, B. F., Shi, Y., Yang, X., Xiao, W., et al. (2014). FTO-dependent demethylation of N6-methyladenosine regulates mRNA splicing and is required for adipogenesis. Cell Res. 24, 1403–1419. doi: 10.1038/cr.2014.151

PubMed Abstract | CrossRef Full Text | Google Scholar