Landscape of genetic and transcriptional changes of methylation modification regulators in ccRCC
Drawing from the published studies, the 9 most common types of methylation modification in DNA, RNA and histone dimensions containing a total of 104 regulators were taken into account (Additional 2: Table S1). 20 regulators of 5mc modification, 23 regulators of m6A modification, 26 regulators of another 6 types of RNA modification and 35 regulators of histone methylation modification were included in the current study [8, 10, 13, 31, 32].
To explore the genetic alterations of different methylation modifications, we first focused on the prevalence of somatic mutations in ccRCC patients. The results were shown in Additional 1: Figure S1A-D. Among 336 samples in TCGA-KIRC cohort, 7.14% (24 samples) experienced mutation of DNA methylation regulators. The mutation frequencies of m6A modification regulators and 6 other types of RNA modification regulators accounted for 7.74% (26 samples) and 4.17% (14 samples) respectively. 17.86% (60 samples) experienced mutation of histone methylation modification regulators. Regulators with higher mutation frequencies were summarized in Fig. 1a. The overall somatic mutation frequencies of regulators were lower and only SETD2 (12%; lysine methylation writer) exhibited a relatively higher mutation frequency. Compared to the wide type group, the expression of SETD2 was downregulated in the mutated group (Fig. 1b). Meanwhile, survival analysis showed that the low expression of SETD2 presented poor overall survival (Fig. 1c).
We further found that copy number variations (CNVs) of these regulators (Additional 1: Figure S1E-H) were prevalent among them, especially in 18 regulators (Fig. 1d). Of these regulators, only YTHDC2, ZBTB38, ADAR and IGF2BP2 showed higher frequencies of CNV gains while others experienced widespread losses of copy number. The positions of copy number alterations on chromosomes of these regulators were shown in Fig. 1e. To inspect whether the copy number alterations had influences on the regulators’ mRNA level in ccRCC patients, we then performed differential expression analysis by using pan-cancer data from TCGA database (Additional 1: Figure S2A-D). As shown in Fig. 1f, only 16 regulators were significantly differentially expressed in ccRCC patients (|logFC| ≥ 0.5; P< 0.05). Of these regulators, 7 were upregulated and 9 were downregulated. Comparing gene expressions with copy number alterations (Fig. 1d and 1g), we found that downregulated genes always showed CNV losses (e.g., KDM1A and SETD2) while upregulated genes experienced CNV gains (e.g., YTHDC2). The results suggested that copy number alterations may play an important role in dysregulation of regulators. At the same time, we noticed that the upregulation of WTAP was accompanied by a high frequency of CNV loss. To explore the reason for these contradictions, we divided the samples into 4 groups according to their CNV values: normal, CNV_Loss, none_ CNV and CNV _Gain. In fact, the CNV loss group had lower expression of WTAP than the normal group (Additional 1: Figure S2E). Nevertheless, contrary to expectations, copy number alterations of IGF2BP2 and TET2 presented the opposite trends with their expression.
Overall, the biogenesis and development of tumors is an intricate process. CNV and somatic mutations couldn’t completely describe the dysregulation of regulators and just offered one possible explanation. Those results showed high heterogeneity, which happened in a landscape of genetic and transcriptional changes of methylation modification regulators in ccRCC. It indicated that the dysregulation of methylation modification regulators had potential values for further study.
Prognostic value and crosstalk patterns of different dimensional methylation modification regulators
In order to improve the accuracy and specificity of analyses. 1033 ccRCC samples from 5 datasets (TCGA-KIRC, E-MTAB-1980, GSE46699, GSE53757, GSE73731) were selected for a more in-depth study. We divided 104 regulators into 4 categories: DNA methylation modification, m6A modification, other RNA modification and histone methylation modification. Univariate Cox regressions were performed to determine the prognostic value. As shown in Fig. 2a-2b, 60 of 104 modification regulators were significantly associated with prognosis.
As the important part of epigenetics, modifications mediated by these regulators enriched the central dogma. We also wondered if there was any mutual influence between them. We performed correlation analyses to demonstrate the interactions of 104 regulators (Additional 1: Figure S2F-I and S3A-F). In DNA methylation modification (5mc), all regulators had shown significantly positive correlations (Additional 1: Figure S2F). Unexpectedly, our results showed that high expression of writers were not associated with low expression of erasers, and even exhibited very strong positive correlations. When it comes to RNA modification, in m6A modification, negative correlations had much lower frequencies than positive correlations (Additional 1: Figure S2G). Only m6A reader IGF2BP2 and IGF2BP3 showed negative correlation with many other regulators, particularly m6A erasers FTO and ALKBH5. In other common types of RNA modification, Additional 1: Figure S2H revealed a critical negative regulatory role of NSUN7 (m5C writer). Only NSUN7 presented negative correlations with other regulators including NOP2 (m5C writer), METTL1 (m7G writer) and PUS1 (Pseudouracil writer). We also calculated correlations between other common types of RNA modification and m6A modification (Additional 1: Figure S3D). The results indicated that the positive correlation was still the predominant form, especially between all writers and erasers which was quite beyond our expectations. At the same time, the expression of regulators was not only closely connected in the same category, but also showed a significant correlation in different types of modifications. For example, the m1A modification writer TRMB61B presented very strong correlations with METTL14 (m6A writer) and METTL2A (m3C writer). Even the methylation site or base was completely different, and there was no direct evidence showing that the neighbor sites or the same bases could interfere with each other, but their regulators still presented a high frequency of synergistic effects. In histone modification, PRMT8 (Arginine methylation writer) attracted our attention. It negatively correlated with the vast majority of histone regulators, especially with some Lysine methylation regulators like EHMT2, SETDB1 and KDM4C (Additional 1: Figure S2I). We summarized the results of univariate Cox regressions and correlations and constructed a huge prognosis-correlation network among all DNA, RNA and histone modification regulators (Fig. 2c). P<0.0001 was considered as significant for co-expression relationships.
Recent studies had reported that crosstalk about different dimensions of methylation modification exerted important biological functions in the occurrence and development of tumors [13]. The existing research was mainly focused on the modifications between DNA methylation and histone methylation. The network suggested that many regulators had shown very strong crosstalk patterns in different dimensions or types of methylation regulators. Positive correlations were dominant and modification writers were not necessarily negatively correlated with modification erasers. Focused on the overall analysis of the network, some regulators (such as PRMT8, IGF2BP2 and especially IGF2BP3) might act as important nodes to negatively regulate the network. The above analyses indicated that the cross-talk among regulators is very important in the process of tumorigenesis and has huge values for further study.
Distinct clusters and immune infiltration characteristics of modification regulators
We used the R package “ConsensusClusterPlus” to perform unsupervised clustering with the expression data of 104 methylation modification regulators. 1033 ccRCC samples were divided into 4 clusters, which were named as methylation-modification-regulator-cluster A-D (MMRC A-D). Group composition was recorded in Additional 2: Table S5. We drew a heatmap to display the expression distribution of regulators in different subtypes (Fig. 3a). The results showed that most of regulators were significantly upregulated in MMRC A and downregulated in MMRC C. Only a few regulators exhibited upregulated expression in MMRC C, such as IGF2BP2, IGF2BP3 and UTY. We performed survival analyses for overall survival and used log-rank test to compare the significance (Fig. 3b). The results showed that the ccRCC patients in MMRC A and D exhibited significant survival advantages than MMRC B and C.
To further investigate the biological functions of different clusters, we performed gene set variation analysis (GSVA) using KEGG gene sets and hallmark gene sets. We drew heatmaps to exhibit the top 20 differences in enrichment between each cluster (Fig3c-3d, Additional 1: Figure S4A-F). MMRC B and C with worse survival outcomes were markedly enriched in stromal activation and tumor proliferation, metastasis pathway including EMT signaling pathway, ECM receptor interaction, inflammatory response, IL6-JAK-STAT3 signaling, KRAS and P53 signaling pathway (Fig. 3c, Additional 1: Figure S4C and S4F). The enriched pathways of MMRC A and D were significantly associated with immune activation and tumor-killing pathways such as interferon alpha response, interferon gamma response and B cell receptor signaling pathway (Fig. 3d, Additional 1: Figure S4D-E).
Many articles have reported that methylation modification was closely associated with immune infiltration [20, 21]. And based on the results of GSVA, we also found that some immune-related pathways were activated in different clusters. Therefore, we wanted to inspect their roles in immune cell infiltration. We performed ssGSEA analyses to calculate the abundance of different immune cells (Fig. 3e). Surprisingly, MMRC B and C, but especially MMRC B, were significantly enriched in immune cell infiltration containing natural killer cells, activated CD4+ T cells, dendritic cells, activated CD8+ T cells, macrophages and activated B cells et al. However, MMRC B and C failed to present the matching survival advantages (Fig. 3b). Chen DS and I Mellman had reported that tumors also could reflect immunosuppressive status, even though they were accompanied by abundant immune cells. These immune cells were blocked outside of tumors by the surrounding stromal cells instead of traversing the parenchyma. It was the activation of stroma that prevented T cells from exerting their beneficial functions [33]. We offered our hypothesis that the same might apply to MMRC B. Subsequently, we calculated the stromal score and immune score using R package “Estimate”. As shown in Fig. 3f, the stromal and immune score of MMRC B was remarkably higher than other clusters. What’s more, from the results of GSVA, some stromal activation pathways such as EMT signaling pathway, angiogenesis, ECM receptor interaction and cytokine receptor interaction were remarkably enriched in MMRC B (Fig. 3c, Additional 1: Figure S4A). These results provided further support for our hypothesis. Next, we further analyzed the correlations between methylation modification regulators and immune cells (Additional 1: Figure S5A-D). Our previous results indicated that IGF2BP3 was more meaningful for prognosis and might act as important nodes to negatively regulate the network (Fig. 2c). We further observed that IGF2BP3 was positively related to numerous immune cells, such as activated CD4+ T cells, activated dendritic cells, CD56Bright NK cells and macrophages (Additional 1: Figure S5B). Moreover, the “Estimate” package was used to quantify the difference of immune cell infiltration abundance between high and low expression of IGF2BP3. As shown in Fig. 3a and 3g, the expression of IGF2BP3 was much higher in MMRC B and the upregulated expression of IGF2BP3 presented a higher immune score. This implied that ccRCC patients with high expression of IGF2BP3 such as MMRC B presented a significantly higher level of immune cell infiltration. Meanwhile, all the other regulators presented negative correlations but only IGF2BP3 showed significantly opposite correlations with CD56Bright NK cells (P<0.001). Recent studies paid more attention to the immunoregulatory functions of CD56Bright NK cells and the mechanism of different methylation modifications regulating the distribution and function of CD56Bright NK cells [34-36]. The enrichment of CD56Bright NK cells depended on the expression of chemokine receptors, such as CCR7, CD62L and CXCR3 [37]. At the same time, CD56Bright NK cells and dendritic cells (DCs) can form a positive feedback regulation to enhance their immune function by promoting IL12, IL15, IL18 and IFN-γ expression [38]. In addition, IL-2 was secreted from T cells to regulate the function of CD56Bright NK cells. In turn, CD56Bright NK cells can inhibit the activation of T cells [39]. We analyzed the expression of these inducing factors (Fig. 3h). Our results suggested that the high expression of IGF2BP3 resulted in the comprehensively increasing expression of chemokine receptors and enhanced the abundance of CD56Bright NK cells, consistent with previous analyses. The upregulation of other interacting factors in high IGF2BP3 group indicated that the high expression of IGF2BP3 might regulate immune activity and promote immune-cells interactions. Meanwhile, high expression of IGF2BP3 was associated closely with PD-L1 and might affect the response to immunotherapy. From these, we speculated the IGF2BP3 mediated methylation modification could strengthen anti-tumor immune responses by activating immune cells and immune-related pathways.
Construction of methylation modification regulator signature
To further explore the potential functional and biological characteristics of the different clusters, we performed differential expression analyses between different clusters respectively and drew Venn diagrams to identify 2856 differential expression genes (Fig. 4a). Subsequently, we performed GO enrichment analyses using the differential expression genes. Top 10 significant enrichments were presented in Fig. 4b. We found that these genes related to pathways like antigen processing and presentation, type I interferon production, RNA splicing and epigenetic modification. This further confirmed that methylation modification regulators played an important role in immune infiltration of the tumor microenvironment. We also performed univariate Cox regression analyses among the 2856 genes, and there were 2170 genes that had prognostic value. To further verify the accuracy of our classification method, we used the R package “ConsensusClusterPlus” to perform unsupervised clustering based on the expression of 2170 DEGs. Consistent with regulator clusters, this analysis also divided 1033 ccRCC patients into 4 subtypes, which were named as methylation-modification-gene-subtype A-D (MMGS A-D). Related to MMRC C, nearly all regulators were remarkably downregulated in MMGS C. The expression of most regulators was elevated in MMGS D, except for a few upregulated genes in MMGS A (Fig. 4c). Interestingly, in line with regulator clusters, gene subtypes B and C also presented high abundances of immune infiltration but demonstrated much poorer prognosis (Fig. 4d, Additional 1: Figure S4G). The concordance of function and prognosis in different clusters revealed the effectiveness and accuracy of our classification method. Due to the heterogeneity of different modifications and different individuals, these clusters based on a large patient population failed to accurately predict the specific situation of individual patients. To calculate the quantitative indicators of the landscape of methylation modification regulators, we used machine learning algorithm to perform PCA analyses. Based on the expression of 2170 DEGs, we obtained a score for each ccRCC patient which was named as methylation modification regulators score (MMR Score). The alluvial diagram presented the transformation of each individual patients’ attribute (Fig. 4e). Its prognostic value was identified by Kaplan-Meier analysis (log-rank, P<0.001). As shown in Fig. 4f, ccRCC patients with high MMR Score had a poor prognosis. In line with this, the increasing MMR Score was accompanied by worse clinical outcomes, including higher degrees of tumor grade and tumor stage, especially pathological T stage (Fig. 4g-4i).
Characteristics of MMR Score in immune checkpoints and somatic mutation
Immune checkpoint blockade treatment provided a huge breakthrough in cancer therapy. Before the initiation of immunotherapy, the results of genetic testing were necessary to guide clinical medication better, including the expression of checkpoints, tumor mutation burden (TMB) and microsatellite instability (MSI). We firstly evaluated the condition of immune activation in ccRCC patients. CD8A, TNF, TBX2, PRF1, GZMA, GZMB, CXCL9 and CXCL10 were selected as immune activity relevant signatures [40]. As shown in Fig. 5a-5b, most of these immune activity related genes were significantly downregulated in high MMR Score group and exhibited strong negative correlations with MMR Score. It might indicate a condition of immunosuppression and reasons for poor prognosis in high MMR Score group. Meanwhile, we also evaluated the expression of common immune checkpoint blockades such as PD-1 and CTLA4 [41]. We found that PD-1 was remarkably upregulated in high MMR Score group (Fig. 5a-5b). Even though the efficacy of immunotherapy did not directly relate to the expression of checkpoints, but it still indirectly implied that ccRCC patients with high MMR Score might be candidates for anti-PD-1 immunotherapy. Survival analyses showed that high expression of PD-1 presented a worse prognosis (Fig. 5c). Combined with MMR Score, the low-PD-1-low-MMRScore group showed a significant survival advantage and provided a more excellent predictive value (Fig. 5d). More and more evidence has shown that the high frequencies of TMB and MSI could increase immune cell infiltration and induce them to specifically recognize and eliminate cancer cells [42, 43]. In the clinical trial (KEYNOTE 012), high frequencies of TMB were confirmed to be able to increase the response to PD-L1 blockades and improve patient outcomes [44]. As shown in Fig. 5e and 5g, we found that ccRCC patients in high TMB status were accompanied by high MMR Score and poor prognosis. The frequencies of TMB have shown a positive correlation (R=0.18, P<0.05) with MMR Score (Fig. 5f). And the combination of high TMB and high MMR Score showed the worst survival outcome (Fig. 5h). The results provided more evidence to imply that ccRCC patients with high MMR Score might be more suitable for immunotherapy. The condition of MSI also drew our attention. Unlike what we expected, the alteration was very slight and failed to reach statistical significance in the low and high MMR Score groups (Additional 1: Figure S5E-F). There was also no statistical significance for survival analysis (Additional 1: Figure S5G). In general, most results indirectly suggested that MMR Score might act as a potential predictive biomarker and patients with a high MMR Score were more likely to derive benefits from immunotherapy, especially anti-PD-1 therapy. To further explore the distribution of somatic mutation, the “Maftools” package was used to present the distribution differences of altered genes. As shown in Fig. 5i-5j, the high MMR Score group had a higher overall mutational frequency. The alteration frequencies of VHL (45%), TTN (17%), SETD2 (15%) and BAP1 (10%) were remarkably higher than in the low MMR Score group. These alterations provided new perspective on molecular mechanisms and choices for treatment. It also helped us better understand the mechanism of somatic mutation in different dimensional methylation modifications and the potential predictive value of the MMR Score in immunotherapy.
The role of MMR Score in response to immunotherapy
ICB therapies represented by anti-PD-L1 and anti-PD-1, which play a key role in blocking the suppression of T cells, has shown better outcomes in cancer treatment, but the efficacy is not available for each individual. The specific criteria of selection and application were still unknown, but necessary. To further verify the accuracy and predictive value of methylation modification model in response to immunotherapy, two immunotherapy cohorts were included in analyses. In the ccRCC immunotherapy cohort (PMID: 32472114), 172 patients who had received Nivolumab (anti-PD-1) treatment were assigned low and high MMR Score and their objective responses were divided into CR/PR (CR: complete response; PR: partial response) and SD/PD (SD: stable disease; PD: progressive disease) groups [45]. It was a pleasant surprise that patients with high MMR Score showed higher rates of therapeutic advantage and clinical benefits (Fig. 6a). At the same time, patients who developed metastases exhibited higher MMR Score (Fig. 6b). We also further validated our analyses in the IMvigor210 immunotherapy cohort (urothelial cancer patients received anti-PD-L1treatment). The results suggested that patients with high MMR Score exhibited similar clinical benefits (Fig. 6c-6d). Especially in distant metastatic lesions, such as renal pelvis and ureter, MMR Score presented more excellent value to guide immunotherapy (Fig. 6e-6g). In consideration of the relationship between MMR Score and tumor microenvironment, we used the “Cibersort” package to quantify the abundance of immune cells. Compared with SD/PD group, dendritic cells, mast cells and macrophages presented significant differences in immune infiltration, and macrophages decreased in CR/PR group (Fig. 6h). At the same time, rather than in low MMR Score patients, the abundance of macrophages downregulated in high MMR Score patients (Fig. 6i). And the M1 and M2 marker genes were also remarkable downregulated in high MMR Score group (Fig. 6j). This gave us new inspiration. Recent studies have shown that the PD-1 monoclonal antibodies could be captured by macrophages and failed to exert its anticancer effects [46]. It implied that the dysregulation of macrophages might play a crucial role in the tolerance in immune checkpoint blocking therapies. To further investigate the molecular mechanism between macrophages and anti-PD-1/PD-L1 immunotherapy would provide a new approach to improve the efficacy of immunotherapy.