Landscape of genetic variations in m5C regulators in HCC
Ultimately, 11 writers, 1 eraser and 1 reader were identified. First, the incidence of CNVs and somatic mutations in regulators in HCC were summarized. In 364 samples, 41 showed mutations in m5C regulators, the occurrence of which was 11.26%. DNMT1 was found to be exposed to a higher frequency of mutations, followed by DNMT3A, while ALYREF, NSUN2, NSUN3, and NSUN5 were not (Fig. 1a). CNVs were also detected in 13 other regulators upon exploration of their modification frequencies. Most of the modifications involved a copy number expansion, but TET2, NOP2, and NSUN4 had a broad occurrence of deletions (Fig. 1b). The chromosome sites of the m5C regulators are shown in Fig. 1c. Based on the expression of 13 m5C regulators in HCC patients, HCC samples could be thoroughly differentiated from normal samples (Fig. 1d). To determine whether the expression of m5C regulators was influenced by the genetic mutations mentioned above, the mRNA expression of regulators was explored. We found that a change in m5C was an important factor leading to perturbations in the expression of m5C regulators. Compared with normal hepatic tissues, the expression of m5C regulators with a CNV expansion was significantly higher than that in HCC tissues (e.g., ALYREF and NSUN2) (Fig. 1b and e). The analyses above showed that the genetic and expression alteration landscape of m5C regulators in normal tissues and HCC tissues is highly heterogeneous, suggesting that the expression imbalance of m5C regulators plays an important role in HCC occurrence and progression.
m5C methylation alteration patterns mediated by 13 regulators
Univariate Cox regression analysis showed that 13 m5C modulators have prognostic significance in HCC patients (Fig. 2a). The m5C regulator network revealed m5C modulator interactions, modulator connections and their prognostic significance for patients (Fig. 2b). The R package ConsensusClusterPlus was applied to classify patients with qualitatively different m5C alteration patterns according to the expression of 13 m5C regulators, and unsupervised clustering analysis was performed to identify a total of 3 different modification patterns (120 cases in modification pattern 1, 178 cases in modification pattern 2, and 73 cases in modification pattern 3; referred to as m5C Clusters 1–3, respectively) (Fig. 2c and Table S1). The prognostic analysis of the three major m5C modification subtypes showed that Cluster-2 had a clear survival advantage over the others (Fig. 2d). The above results indicate that the crosstalk between erasers, readers and writers may play an important role in m5C alteration patterns and TME cell infiltration characteristics between individual tumours.
TME cell infiltration characteristics in different m5C modification patterns
To investigate the biological actions associated with m5C modification patterns, GSVA was conducted. As shown in Fig. 2e and Table S2, m5C Cluster-2 was remarkably enriched in the activation of stroma and carcinogenesis pathways, such as the ERBB signalling pathway, cell cycle signalling pathway, and cell adhesion pathway. Cluster-1 was associated with enrichment pathways (Fig. 2e). Cluster-3 was highly associated with carcinogenic activation (Fig. 2f). Unexpectedly, further analysis of TME cell infiltration showed that Cluster-1 was significantly enriched in the infiltration of innate immune cells, including eosinophils, NK cells, macrophages, CD8 T cells, mast cells, and active B cells (Fig. 3a). Compared with patients in Cluster-3, patients in Cluster-2 showed a corresponding survival advantage (Fig. 2d). Prior research has shown that tumours with an immune rejection phenotype exhibit large amounts of immune cells, and these immune cells are in the matrix around the tumour cell nest instead of inside the tissue [30]. GSVA showed that the modification of Cluster-1 was significantly related to matrix activation. Therefore, it was speculated that the Cluster-1 matrix serves as an activation inhibitor of the anti-tumour effect of immune cells. Further analysis showed that matrix activity was greatly upgraded, activating the angiogenesis pathway. These results supported our hypothesis (Fig. 3b). Based on the above analysis, we found that the 3 m5C modification patterns had obvious TME cell infiltration characteristics. Cluster-1 had an immune rejection phenotype, featuring natural infiltrating immune cells and an activated matrix; Cluster-2 had an immunoinflammatory phenotype, featuring adaptively infiltrated immune cells and immune activation; and Cluster-3 had an immune desert phenotype, featuring immune suppression.
The m5C regulator DNMT1 has a strong relationship with infiltrating immune cells
To further explore the role of each m5C regulator in the TME, Spearman correlation analysis was applied to examine the correlation between each TME-infiltrating cell type and m5C regulators (Fig. 4a). An emphasis was placed on the regulator DNMT1, an m5C methyltransferase, and we revealed its positive relationship with the infiltration of many TME immune cells. An estimation method was applied to determine the expression of DNMT1 and the infiltration of immune cells. The results showed that higher DNMT1 expression was related to a higher immune score, which means that a TME with high DNMT1 expression has significantly high immune cell infiltration (Fig. 4b). Based on these results, the specific differences in 23 TME-infiltrating immune cells were explored between patients with high and low DNMT1 expression. We found that tumours exhibiting high DNMT1 expression had markedly more infiltration of 13 TME immune cells than those exhibiting low expression (Fig. 4c). Recently, attention was drawn to the regulatory mechanisms of m5C modification on the activation of DCs, which are the bridge connecting innate immunity with adaptive immunity, the activation of which depends on upregulating the expression of MHC molecules, adhesion molecules and costimulatory molecules (Fig. 4d). As expected, subsequent enrichment analysis showed that tumours with high DNMT1 expression showed remarkable enrichment in immune activation pathways (Fig. 4e). Therefore, it was speculated that m5C methylation modification mediated by DNMT1 may contribute to activated DCs in the TME, thus promoting the anti-tumour immune response in HCC.
High expression of the m5C regulator DNMT1 in tumour tissues is related to a poor prognosis in patients with HCC
IHC staining was used to determine the expression pattern of DNMT1 on a TMA consisting of 90 pairs of HCC tissues and adjacent tissues. Representative micrographs illustrate the various degrees of DNMT1 expression (Fig. 5a and b). The expression of DNMT1 was higher in tumour tissues than in control tissues (Fig. 5c), which was consistent with the findings in the TCGA-LIHC cohort (Fig. 1e). The correlation of DNMT1 expression with the clinicopathological characteristics of patients with HCC is shown in Table S3. In addition, Kaplan-Meier curve analysis showed that patients with high DNMT1 expression had shorter overall survival (OS) than those with low DNMT1 expression (Fig. 5d). Univariable and multivariable Cox regression analyses were used to determine whether the expression of DNMT1 was an independent risk factor. The univariable analysis revealed that DNMT1 expression was associated with tumour size and TB, AFP, and PD-L1 levels (P < 0.05, Table S4). Further analysis demonstrated that DNMT1 might serve as a prognostic predictor for HCC.
Generation of the m5C gene signature and functional annotation
For subsequent exploration of the biological behaviour of each m5C modification pattern, we ascertained 307 m5C phenotype-related DEGs with the limma package (Fig. 6a). clusterProfiler was employed to implement enrichment analysis on the DEGs. Table S5 summarizes the significantly enriched pathways. As expected, we detected enrichment in biological processes that are notably related to m5C modification and immunity, which verified the important role that m5C modification plays in immune regulation in the TME (Fig. 6b).
To further explain the association, we performed unsupervised clustering analysis to classify 307 m5C phenotype-related genes and extracting 27 immune-related genes: VIPR2, CCL7, RBP2, SLC10A2, FGF5, DEFA5, HTR3A, TRH, LCN15, AMBN, ADIPOQ, FGF3, CCK, NTF4, NDP, FGF9, PF4, CMA1, SFTPA2, CGB8, DEFA6, PF4V1, IL25, GH2, FGF8, SST and IAPP. Furthermore, we performed unsupervised clustering analysis based on these genes to categorize patients into different subtypes (Figure S1a-d). In line with the clustering analysis of m5C modification patterns, unsupervised clustering analysis revealed three different m5C-modified phenotypes termed Immu-clusters 1–3, respectively. Thus, there are 3 different distinct immune-related m5C methylation patterns. We observed that tumours in Immu-clusters 2 and 3 were associated with poor differentiation and enriched in diffuse histological subtypes. The opposite pattern was observed in Immu-Cluster 1. Patients whose survival status was known were mainly concentrated in Immu-cluster 1, while patients in clinical stage IV or with a high TNM grade were mainly concentrated in Immu-cluster 2 (Fig. 6c). The analysis also showed that three different gene clusters had different feature genes (Fig. 6c). In total, 114 of the 317 HCC patients clustered in Immu-cluster 1, which was associated with a better prognosis. The prognosis of patients in Immu-cluster 1 (110 patients) and immu-cluster 3 (93 patients) was poor (Fig. 6d). In the three immune clusters, a significant distinction in the expression of m5C regulatory factors emerged. This result was consistent with the m5C methylation modification patterns (Fig. 6e).
Clinical and transcriptional features of the m5C-related phenotypes
To further explain the role that m5C-related phenotypes play in TME immune regulation, the levels of immune cells and expression of chemokines and cytokines in the three Immu-clusters were examined. The chosen cytokines and chemokines were taken from previously existing studies. Our analysis showed that activated CT4 T cells, immature B cells, regulatory T cells, NK cells, macrophages, mast cells, myeloid-derived suppressor cells (MDSCs), monocytes, n neutrophils and plasmacytoid DCs were significantly different among the Immu-clusters (Fig. 7a). Tumour necrosis factor, interferon, CD8A, CXCL9, CXCL10, GZMA, GZMB, PRF1 and TBX2 were associated with immune activation transcription (Fig. 7b) [23, 31]. PD-L1, CD80, CD86, CTLA-4, HAVCR2, etc., were thought to be related to the transcription of immune checkpoints. We compared the transcription of these immune checkpoint genes in the three Immu-clusters (Fig. 7c). ACTA2, CLDN3, VIM, COL4A1, SMAD9, TWIST1, TGFBR2, TGRB1 and ZEB1 are related to the transcription of growth factor β/EMT pathway transformation and exhibited significant differences between the three Immu-clusters (Fig. 7d). We found that mRNAs related to the TGF-β/EMT pathway were significantly upregulated in Immu-cluster 2 and Immu-cluster 3, indicating that this cluster is the matrix-activated group. Immu-cluster 1 showed elevated expression of mRNAs related to activated immune transcription. The above results indicated that Immu-cluster 1 could be categorized into an immune-activated group. The above results also showed that m5C modification has a non-negligible regulatory effect on the formation of different TMEs.