Molecular Characterization of m6A Modifications in Non-Clear Cell Renal Cell Carcinoma and Potential Relationship with Pathological Types

Introduction

Renal cell carcinoma (RCC) originates from the renal tubular epithelium and is currently recognized as one of the 10 most common cancers, killing approximately 180,000 individuals worldwide each year.^1–3 In recent years, the classification of renal cancer has been repeatedly revised based on its histological and genomic molecular features,^4–7 and 16 different disease subtypes have been designated in the latest classification.⁷ Of these RCC subtypes, clear cell RCC is the most common, accounting for approximately 70–80% of the total cases,^8,9 while papillary RCC (pRCC) (10–20%) and chromophobe RCC (chRCC) (> 5%) are predominant among non-ccRCC.^10,11 Other subtypes, with an incidence of less than 1%, are rare and include translocation Xp11.2 or 6p21 RCC, collecting duct carcinoma, renal medullary carcinoma, and oncocytoma, among others.¹²

All RCCs in the lining of the proximal convoluted tubule, and even though they are clinically treated similarly, the complex heterogeneity between or within subtypes leads to different clinical outcomes.^13–17 Several studies have used the epigenomics and gene transcriptomics approach to explore ccRCC, while non-ccRCC has received relatively less attention. However, non-ccRCC has a more complex tissue phenotype, intra-tumor heterogeneity, poor immunotherapy response, and prognosis. Therefore, studies evaluating the mechanisms underlying non-ccRCC are essential and may be beneficial for developing individualized treatment strategies and improving the clinical benefit of RCC patients.

RNA modifications are widely observed in nature. More than 170 types of RNA modifications in different types of RNAs, such as mRNA, tRNA, rRNA, snRNA, miRNA, and lncRNA, have been identified. N6-methyladenosine (m6A) is the most common form of RNA modification, accounting for approximately 60% of all modifications; it is dynamic and reversible.^18–20 The m6A RNA modification involves three types of regulatory proteins—methyltransferases (writers), m6A binding proteins (readers), and demethylases (erasers).²¹

Methyltransferases are responsible for catalyzing the modification of m6A sites, while the demethylases play an opposite role and mediate the reversal of m6A modification. The m6A-binding proteins are mainly involved in the recognition of m6A methylation sites.^22,23 m6A modifications are associated with several biological functions, including translation,²⁴ nuclear export,²⁵ RNA stability,²⁴ and alternative splicing.^26,27 In addition, m6A is also involved in regulating key cellular processes, such as stem cell renewal, cell differentiation, and DNA damage response.²⁸

Several studies have reported that m6A modifications promote the development and progression of various malignancies by regulating gene expression and play a fundamental role in prognosis.^18,29 However, little is known about the function of m6A modification and their associated genes in non-ccRCC. In this study, we systematically investigated the roles of 23 m6A regulators in non-ccRCC based on relevant data from The Cancer Genome Atlas (TCGA) database using bioinformatics scores.

Materials and Methods

Data Acquisition

All non-ccRCC-related data (RNA-seq transcriptome, genetic mutations, and corresponding clinical datasets) were obtained from the public network TCGA³⁷ (https://portal.gdc.cancer.gov/repository), including 291 kidney renal papillary cell carcinoma (KIRP) samples and 113 kidney Chromophobe (KICH) samples. Copy number variation (CNV) data were obtained from the shared database (https://xena.ucsc.edu/). A total of 23 regulators were extracted from previous studies, including eight writers²¹ (METTL3, METTL14, METTL16, WTAP, VIRMA, ZC3H13, RBM15/RBM15B), 13 readers²² (YTHDC1/YTHDC2, YTHDF1/YTHDF2/YTHDF3, HNRNPC, FMR1, LRPPRC, HNRNPA2B1, IGFBP1/IGFBP2/IGFBP3, RBMX), and two erasers³⁰ (FTO and ALKBH5).

Cluster Analysis Based on m6A Regulators

Gene expression levels were quantified using the metric log2 (TPM +1). Clustering analysis based on the expression levels of 23 m6A regulators was conducted using the ‘ConsensusClusterPlus 1.57.0’ R package. The optimal number of clusters was determined based on the Bayesian information criterion.

Estimation of Immune Infiltrating Cells

The single-sample gene-set enrichment analysis (ssGSEA) algorithm was used to quantify the relative abundance of each type of infiltrated immune cell in the non-ccRCC tumor microenvironment. The enrichment scores calculated by ssGSEA analysis were utilized to represent the relative abundance of each TME infiltrating cell in each sample. Based on a similar study by Charoen Tong, we evaluated a total of 22 types of immune cells,^38,39 including activated CD8 T cell, natural killer T cell, regulatory T cell, activated dendritic cell, macrophage and so on.

Gene Set Variation Analysis (GSVA) and Functional Annotation

We used the “GSVA 1.41.3” R package to analyze the differences in pathways and biological processes with the different m6A modification patterns. Gene sets of “c2.cp.kegg.v7.4.symbols.gmt” and “c5.go.v7.4.symbols.gmt” were downloaded from the Molecular Signatures Database (MSigDB) (https://www.gsea-msigdb.org/gsea/index.jsp). P-values less than 0.05, after adjustment, were considered statistically significant.

Identification of Differentially Expressed Genes (DEGs) Among m6A Clusters

To identify m6A-related DEGs between different modification patterns, we used the Limma R package to perform a differential analysis of the two groups of m6A clusters.⁴⁰ The significance criteria for determining DEGs were set as an adjusted P-value < 0.001.

Construction of the m6A Gene Signature

We computed the m6Ascore to quantify the m6A modification pattern in a single non-ccRCC sample. First, we performed a univariate Cox analysis on the DEGs of two m6A clusters (in the case of multiple clusters, considered the overlapping genes of every two clusters first and then the intersections) to clarify the relationship between the DEGs and prognosis. DEGs that significantly correlated with clinical prognosis were used for subsequent analyses (P < 0.001). We then performed principal component analysis (PCA) to construct the m6A correlation signature. Both principal components 1 and 2 were selected to act as signature scores. Finally, similar to the genomic grade index (GGI) approach, an m6Ascore was obtained for each sample.

Correlation Between m6Ascore and Other Relevant Phenotypes

The correlation between m6Ascore and other phenotypes was determined using Spearman correlation analysis. Clinical correlation survival analysis was performed based on the Kaplan-Meier curve using the Survival R package.

Statistical Analysis

The correlation between m6A regulators and immune infiltrating cells was calculated based on Spearman correlation analysis. One-way ANOVA and Kruskal–Wallis tests were used to compare differences among three or more groups.⁴¹ The optimal cut-off points for each data set subgroup were determined using the Survminer R package. The Kaplan-Meier method was used to plot prognosis-related survival curves, and the Log rank test was used to assess the significance of the differences. Univariate and multiple Cox regression models were used to assess the risk ratios between m6A-associated modifier genes and various phenotypes. The visualization of mutation data and CNV was done using the Maftools and RCircos R package, respectively. The aforementioned operations were performed using R software version 4.1.0, the statistical P values were all bilateral, and P < 0.05 was considered statistically significant.

Results

Molecular Characterization of m6A Regulators in Non-ccRCC

Based on previous studies, we identified 23 m6A regulators, including eight writers²¹ (METTL3, METTL14, METTL16, WTAP, VIRMA, ZC3H13, RBM15/RBM15B), 13 readers²² (YTHDC1/YTHDC2, YTHDF1/YTHDF2/YTHDF3, HNRNPC, FMR1, LRPPRC, HNRNPA2B1, IGFBP1/IGFBP2/IGFBP3, RBMX), and two erasers³⁰ (FTO and ALKBH5). We first analyzed the prevalence of somatic mutations and CNV in the 23 m6A regulators in non-ccRCC. We found that the frequency of both, somatic mutations and CNV in non-ccRCC, was low, and only a few genes showed copy number amplification or deletion (Figure 1A and B, Figure S1A). KIRP has a higher mutation frequency than KICH (Figure S1B), however, the total mutation rate in KIRP samples did not exceed 8% (Figure 1C). In addition, we compared the differential expression of m6A regulators in normal and tumor tissues in non-ccRCC. We found that expression levels of 17 out of the 23 regulators were significantly different in non-ccRCC (Figure 1D). The above results suggest significantly heterogenous tumorigenesis in patients with non-ccRCC and that gene mutations and CNV are not always the main causes of non-ccRCC, and instead, aberrant modifications of m6A regulators may be crucial in non-ccRCC. We also explored the correlation between m6A regulators and non-ccRCC prognosis and constructed a prognostic correlation network graph (Figure 1E). We found that 17 out of the 23 m6A regulator expression variants had significant oncogenic properties associated with the poor prognosis of non-ccRCC and a synergistic relationship was observed among m6A readers, erasers, and writers. Of these 17 regulators, abnormal modifications of IGFBP2/IGFBP3, VIRMA, RBM15, and RBMX genes play an important role in non-ccRCC and are strongly associated with poor prognosis.

Figure 1 Molecular characteristics and clinical relevance of m6A regulators in non-ccRCC. (A) Dumbbell plot depicting the copy number variation (CNV) alteration frequency of m6A regulators in kidney renal papillary cell carcinoma. Amplification frequencies are marked with red dots and deletion frequencies are marked with blue dots. (B) CNV alterations of m6A regulators on 23 chromosomes in the TCGA-KIRP cohort. (C) Mutation frequency of m6A regulators in KIRP. (D) Gene expression alterations among m6A regulators. (E) Interaction of m6A regulators in ccRCC. Erasers are indicated in red readers are indicated in yellow, and writers are indicated in blue. The size of each circle represents the survival impact of each m6A regulator. Green dots represented favorable and purple represent risk factors of overall survival. The lines connecting m6A regulators present their interactions, and the thickness of the lines correspond to the correlation strength among the. regulators. Positive correlation is indicated in red and negative in blue. *P < 0.05; **P < 0.01; ***P < 0.001.

Sample Clustering Based on 23 m6A Regulators

Based on the expression levels of the 23 m6A regulators, we clustered the non-ccRCC samples (Figure 2A–C). The 357 samples were split into two clusters (m6A cluster_A and m6A cluster_B). To investigate the clinical significance of the two subtypes, we first performed a survival analysis of the two sample groups. We found that m6Acluster_A had a significant (P < 0.001) survival advantage over m6A cluster_B (Figure 2D). We then investigated the expression of m6A regulators in the two clusters. As shown in Figure 2E, almost all the m6A regulators were highly expressed in m6A cluster_B. To further investigate whether there are differences in relevant functions and pathways between the different subtypes, we next conducted a GSVA. We found that gene transcription, cell differentiation, and other related pathways including the cell cycle, basal transcription factors, spliceosome, RNA degradation, regulation of gene expression, epigenetics, mRNA transport, and Notch signaling pathway (Figure 2F and G), were abnormally enriched in samples that formed the m6A cluster_B. The m6A cluster_A was mainly enriched in some metabolism-related pathways, such as phenylalanine metabolism, tyrosine metabolism, and drug metabolism. This suggests that in non-ccRCC patients, aberrant modifications of RNA have different functions that affect transcriptional activity by regulating upstream pathways on the one hand and that affect metabolism by regulating downstream pathways on the other. Therefore, we postulate that aberrant RNA modifications that regulate upstream pathways have a higher malignant potential and affect prognosis more than those that regulate downstream pathways.

Figure 2 Typing and biological properties of m6A methylation modifications in non-ccRCC. (A) Consensus clustering matrix for k = 2. (B) Relative change in area under the CDF curve for k = 2–9. (C) Principal component analysis of non-ccRCC samples from TCGA dataset. (D) Survival curve for the two clusters (m6A cluster-A, blue; m6A cluster-B, yellow) of non-ccRCC patients. (E) Comparison of heatmap and clinicopathological features between different subtypes. Red indicates upregulated regulators and green indicates downregulated regulators. The intensity of the shade corresponds to the expression. The m6A cluster, fustat, T stage, gender, age, and histology were used as patient annotations. (F and G) Heat map demonstrating the differences in biological processes between different m6A clusters. Red indicates activation and blue indicates inhibition of biological processes. F shows the KEGG and G shows the BP analysis respectively.

Immunological Characteristics in Different Subtypes of m6A

To determine whether infiltrating immune cells between the different m6Aclusters in the non-ccRCC tumor microenvironment (TME) are different, we conducted a ssGSVA. We found that the immune cell infiltration characteristics in the two m6A modification patterns were different (Figure 3). The m6Acluster_A was abundant in natural immune infiltrating cells, including activated CD8+ T cells, monocytes, CD56 dim/bright natural killer cells (NK), and eosinophils. The m6Acluster_B was abundant in acquired immune cells, including Th17 cells, Th2 cells, immature dendritic cells (DCs), neutrophils, and plasmacytoid dendritic cells. This may explain the better survival advantage of cluster A patient. It has previously been reported that the degree of Th2 cell infiltration and IGF2BP3 expression are positively correlated with and impaired survival is observed in anti-PD-1 cohorts.³¹ Our observation in the m6Acluster_B corroborates this finding. In conclusion, different m6A modification patterns lead to different patterns of immune cell infiltration that ultimately alter the prognostic outcome of patients.

Figure 3 Box plot comparing the abundance of immune cell infiltrates between different m6Aclusters. The upper and lower ends of the boxes represented the interquartile range. The lines in the boxes represented median value, and blue(yellow) dots indicate outliers. The asterisks represent the statistical p-value. *P < 0.05; **P < 0.01; ***P < 0.001.

Genomic Typing and Clinical Features of m6A Characteristic Genes

To investigate the potential biological consequences under different m6A modification patterns, we identified 1251 phenotypically related differential genes based on the Limma R package and conducted Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis. Notably, we observed an enrichment of tumor-associated biological behaviors such as cell development, proliferation, differentiation, and migration, including the development of the renal system, ameboidal-type cell migration, PI3K-Akt signaling pathway, and TGF-β signaling pathway (Figure 4A and B, Figure S1C). Additionally, we found considerable activation of the cellular matrix environment, such as cellular matrix organization, cytoarchitectonic organization, skeletal system development, cell-substrate adhesion, and growth factor binding (Figure S1D). We then subjected 1251 identified differential genes to unsupervised clustering analyses and successfully classified the samples into three genomic subtypes (gene Cluster A-C) (Figure 4C). Although the clustering grouping with m6A modification pattern was slightly different, we found that (Figure 4D) geneCluster C had the same m6A modification pattern as m6A cluster_B, and geneCluster A-B had the same modification pattern as m6A cluster A. Notably, all TCGA-KICH samples were distributed in geneCluster_A with the same methylation modification pattern. Whereas TCGA-KIRP samples were distributed in different gene clusters (B-C) with different m6A modification patterns. These results suggest that there are different RNA methylation modification patterns and regulatory mechanisms in the same pathological subtype of RCC (pRCC) and highlight again, the obvious heterogeneity of tumorigenesis. The survival analysis of the three gene clusters also showed significant survival differences, with geneCluster A-B having a significant survival advantage over geneCluster_C (Figure 4E). Significant variability in the expression of m6A regulators was also observed (Figure 4F), which is consistent with the expected results for the m6A methylation modification pattern.

Figure 4 Genomic typing and clinical features of m6A characteristic genes. (A and B) Depict GO and KEGG analysis of m6A characteristic differential genes, respectively. (C) Consensus clustering matrix for k = 3. (D) Heat map of unsupervised clustering of characteristically different genes dividing the sample into geneCluster A–C. The m6Acluster, fustat, stage, gender, age, and histology were used as patient annotations. (E) Survival curve for geneClusters (geneCluster_A, blue; geneCluster_B, yellow; geneCluster_C, red) of non-ccRCC patients. (F) Box plot of differential expression of m6A-related regulators in genecluster A–C. The lines in the boxes represented median value, and blue(red) dots indicate outliers. The asterisks represented the statistical p-value. ***P < 0.001.

m6Ascore Signature and Related Prognostic Value

Considering that the aforementioned findings are based on population studies and cannot accurately predict the m6A methylation modification patterns in individual samples, we computed the m6Ascore to further assess the m6A modification patterns of individual non-ccRCC patients. Univariate Cox regression analysis was first performed on 1251 differential genes associated with m6Acluster typing, and 714 prognosis-associated genes were extracted (P < 0.001) for conducting a PCA to establish the m6Agene signature. Samples were divided into two groups—high-m6Ascore and low-m6Ascore, according to the optimal cutoff value. Survival analysis (Figure 5A) revealed that the low-m6Ascore group showed significantly better survival than the high-m6Ascore group (P < 0.01), with a significant difference in 5-year survival (70% vs 20%), suggesting that the m6Ascore signature has potential prognostic value as a biomarker. Furthermore, including m6Ascore, clinical stage, sex, and age in multivariate Cox analysis, the m6Ascore was able to act as an independent prognostic factor for non-ccRCC (Figure 5B). Visualizing the relationship between m6Acluster, m6Ascore, and survival status of non-ccRCC patients (Figure 5C), we observed that almost all samples in the m6Acluster_B group were distributed in the high-m6Ascore region, while the highest proportion of patients who died originated from the high-m6Ascore group. Subsequently, we analyzed the correlation of m6Ascore with immune cells (Figure 5F). The results showed that the m6Ascore was positively correlated with Th17 and Th1 and negatively correlated with eosinophils, CD56dim NK cells, and neutrophils. In addition, both m6Acluster_B and geneCluster C had a higher m6Ascore than the other subtypes, and the difference was significant (Figure 5D and E). Investigating the correlation between m6Ascore and mutational load we found a significant survival impairment in the H-TMB+H-m6Ascore group (P < 0.01) (Figure 6A). We further analyzed the correlation between m6Ascore and clinical phenotype, and significant differences were also observed (Figure 6B–D). To further explore the response of non-ccRCC to immunotherapy, we compared the differential PD-L1 expression between the high- and low-m6Ascore subgroups. The results did not reveal any significant differences in expression (Figure 6E).

Figure 5 m6Ascore signature and related prognostic value. (A) Survival curve based on m6Ascore (blue, low- m6Ascore group; red, high-m6Ascore group). (B) Forest plot of multivariate Cox regression analysis for m6Ascore in TCGA (non-ccRCC) cohort. (C) Alluvial diagram depicting the relationship between m6Aclusters, geneClusters, m6Ascore, and fustat. (D and E) Differences in m6Ascore between m6Aclusters (geneClusters) in TCGA cohort. P < 0.001, Kruskal–Wallis test. (F) Relationship between m6Ascore and immune cells, red represents positive and blue represents negative correlation.

Figure 6 Clinical phenotypes associated with the m6Ascore. (A) Survival curves based on H/L-TMB and H/L-m6Ascore grouping. (B) Comparison of differences in m6Ascore between different survival states. (C) Comparison of survival rates between the high- and low-m6Ascore groups. (D) Comparison of survival differences between the high- and low-m6Ascore groups in T1-2 staging. (E) Comparison of differential PD-L1 expression between the high-m6A and low-m6Ascore groups.

Discussion

Several studies have reported the fundamental role of m6A modifications in development of various diseases, including cancers. However, most studies have mainly focused on the effect of a single m6A regulon on tissue morphological phenotypes or disease phenotypes. A comprehensive understanding of how multiple regulators synergistically influence disease phenotypes is lacking. Clarifying the interrelationship between different regulators, methylation patterns, and disease phenotypes can help us gain a comprehensive understanding of pathogenesis, as well as, facilitate the formulation of individualized medical strategies for better prognosis.

In this study, we first identified two m6A-RNA modification patterns in non-ccRCC based on 23 m6A regulators. Subsequently, we clustered the samples into three gene clusters based on differentially expressed m6A signature genes. Notably, m6A_geneCluster C had the same methylation modification pattern and poorer prognosis as m6ACluster_B, while m6A_geneCluste A-B had the same methylation modification pattern and better prognosis as m6Acluster_A (Figure 4D). Surprisingly, all pathological types of chRCC (KICH) fit almost perfectly with m6A_geneClusterA. This suggests that the methylation modification pattern of RNA in non-ccRCC cells is consistent with the pathological tissue phenotype. We also observed that pathological tissue typing for pRCC (KIRP) corresponded to two methylation modification patterns, as well as m6A_geneCluste B-C.

According to the 2016 version of the WHO renal tumor classification, pRCC can be divided into pRCC-I and pRCC-II; pRCC-I has a better prognosis than pRCC-II. In this context, we speculate that geneClusterB and geneClusterC may correspond to the clinicopathological types of pRCC-I and pRCC-II, respectively. This further suggests that m6A methylation modifications play a decisive role in regulating the histopathological phenotype as well as the prognostic process of tumors. In addition, we found that the methylation modification pattern and prognosis of gene cluster_B was closer to that of gene cluster_A among the three gene clusters. That is, some samples of pRCC (probably pRCC- I) and chRCC have similar modification patterns and prognoses. This suggests that the use of methylation modification patterns as a typing criterion in histotyping may more accurately reflect the prognosis of patients with non-ccRCC.

The unexpected finding of our study that m6A-based genetic cluster typing is highly compatible with non-ccRCC pathological types may not be just a coincidence but some as yet unknown correspondence between m6A genotypes and extrinsic pathological phenotypes. It is well known that pathological type is inextricably linked to clinical prognosis, and m6A is reversibly modified, whether the survival prognosis of patients can be improved by changing m6A modification is unknown, but from the perspective of scientific research, this may provide a novel idea for future individualized treatment of tumors. In addition, m6A-based genotyping for non-ccRCC may have more prognostic value than clinicopathological classification and may serve as an important reference for new non_ccRCC typing.

Analysis of non-ccRCC TME immune-infiltrated cells revealed that the relative enrichment of natural immune cells in m6Acluster_A was associated with better prognosis, while the relative enrichment of acquired immune cells in m6Acluster_B was associated with poorer prognosis. However, no enrichment of immune-related pathways was observed when GSVA was performed. This is in contrast to ccRCC, where infiltration of a large number of immune cells (immunoinflammatory phenotype),^32–34 as well as considerable enrichment of immune activation-related pathways,³⁵ is observed. Notably, in the subsequent GO/KEGG analysis of m6A characteristic differential genes (Figure 4B), we observed a significant enrichment of the TGF-β signaling pathway, which is thought to be a major contributor to the immune-excluded phenotype. The immune rejection phenotype is characterized by significant immunosuppression. Although rejection phenotype is also infiltrated with abundant immune cells, they are mainly distributed in the stroma around the tumor rather than in the tumor parenchyma and do not exert an effective anti-tumor effect, and the immune-excluded phenotype is also associated with a lack of activated T cells.³² We also observed a significant enrichment of biological processes such as matrix activation in the GO analysis (Figure 4A), and previous studies have confirmed that matrix activation has an immunosuppressive profile. In addition, we observed a significant increase in Th2 cell content in m6Acluster_B. Th2 cell infiltration, a marker of immune tolerance has been reported to be positively associated with methylation modifications mediated by high IGF2BP3 expression and to influence poor prognosis after PD-1 blockade.³⁶ This suggests a fundamental role of m6A modification in the development of an immune-excluded phenotype in non-ccRCC. We also revealed the relationship between different m6A modification patterns and tumor immune infiltration heterogeneity. Whether altering m6A modification or inhibiting matrix activation may improve the sensitivity of immunotherapy to improve survival rates in patients with immune rejection such as non-ccRCC has not yet been established but may be a potential treatment in the future.

We used the m6Ascore to assess methylation patterns in individual samples and to eliminate the effect of heterogeneous factors in the overall sample. The combined analysis revealed that m6Ascore had the ability to independently predict patient survival and show significant survival differences across patients with different T stages (Figure 6D). Although there was no direct correlation between the m6Ascore and tumor burden mutation (Figure S1E), H-TMB combined with H-m6Ascore showed a significant survival disadvantage. We also analyzed the relationship between m6Ascore and PDL1 expression and no significant correlation was found between the high or low m6Ascore and PDL1 expression. These results suggest that ccRCC and non-ccRCC may have different m6A methylation modification patterns, regulating different immune phenotypes. Our study is limited in that it does not involve a larger sample size due to the low prevalence of non-ccRCC, as well as because our data included only pRCC and chRCC under non-ccRCC.

Conclusions

Overall, our study unravels the important role of m6A methylation modifications in non-ccRCC and deepens the understanding of the pattern of immune infiltration in non-ccRCC. In particular, m6A-based gene clusters were found to be highly correlated with clinicopathological types, as well as significant prognostic differences. This provides a new perspective for the development of individualized medical strategies.