DNA methylation regulators-related molecular patterns and tumor immune landscape in hepatocellular carcinoma.

Increasing evidence showed that the dysregulation of DNA methylation regulators is a decisive feature of almost all cancer types and affects tumor progressions. However, few studies focused on the underlying influences of DNA methylation regulators-related genes (DMRegs) in immune cell-infiltration characteristics, tumor microenvironment (TME) and immunotherapy in HCC patients. In our study, the alterations of DNA methylation regulators modification patterns (DMRPs) were clustered from hepatocellular carcinoma (HCC) samples based on the expression of DNA methylation regulators as well as genetic and transcriptional features. In addition, based on molecular identification of three distinct molecular subtypes, we found that different DMRPs alterations were related to different clinicopathological characteristics, prognosis, and immune cells infiltration features. Moreover, we constructed and validated a DNA methylation regulators-related genes score (DMRegs_score) to predict the survival of HCC patients. A high DMRegs _score, which was characterized by more TP53 wild mutation, high expression of PD-1, CTLA-4, and remarkable immunity activation, was indicative of poor prognosis. Furthermore, we validated the expression of eight genes which were used for the prognostic signature in this risk score by RT-qPCR using tissues from our center. More importantly, DMRegs_score was highly correlated with targeted drug sensitivity. Additionally, we developed a highly accurate scoring system that could be used to improve the clinical applicability of DMRegs _score. In conclusion, these findings may contribute to a better understanding of DNA methylation regulators and provide new strategies for evaluating prognosis and developing more effective combination therapy for HCC patients.

Introduction

Liver cancer, more specifically hepatocellular carcinoma (HCC), is the main leading cause of cancer-related death worldwide in 2020 (1, 2). At present, surgery is still the most effective treatment for HCC. However, due to the occult onset and rapid progress of HCC, patients often have lost the best opportunity for surgical treatment at the time of diagnosis. What’s more, the patients with HCC have poor prognosis because of high metastasis and recurrence rate (3). Therefore, exploring the molecular mechanism of HCC development and finding new early diagnosis and treatment targets are the focus of HCC research.

Epigenetic modifications, such as DNA methylation, play a crucial role in altering gene expression and contributing to disease development in mammals (4). According to present reports, methylation of the fifth carbon of the DNA cytosine within CpG dinucleotides is the most mechanistically understood form of DNA methylation (5). DNA methylation modification is a dynamic and vary process which is modulated by DNA methylation regulators, including DNA methyltransferases, DNA demethylases and DNA binding proteins (6–8). In addition, a growing body of evidence had demonstrated that dysregulation of DNA methylation regulators is a hallmark of almost all cancer types and affects tumor microenvironment (TME) or immunotherapy (7, 9, 10).

Co-inhibitory receptors Cytotoxic T-lymphocyte antigen 4 (CTLA4), programmed cell death protein 1 (PD-1), and programmed cell death ligand 1 (PD-L1) is expressed in the tumor microenvironment. Immunotherapy such as immune-checkpoint inhibitors (ICIs) that target these biomarkers activated the properties of effector T cells which can be able to kill cancer cells. Importantly, ICIs have radically reversed cancer therapy (11). Cancer immunotherapy targeting CTLA4, PD-L1, or PD-1 has become a widely used method of treating various types of cancer (12–14). Recently, anti-CTLA-4 was reported a survival benefit of HCC patients with sorafenib resistance (15). However, these immunotherapies were responding differently with patient to patient, and less than 20% of immune checkpoint blockade therapy was effective (16–18). It has been reported that the expression of PD-L1 and the status of tumor mutation burden (TMB) may be used as biomarkers to assess the effectiveness of immunotherapy (19–21). Interestingly, a recent study has revealed the characteristics of DNA methylation modification patterns of gastric cancer and explored the link between TME and DNA methylation modification, which indicated that DNA methylation may be a new predictor for immunotherapy (8). Moreover, DNA methylation regulators distinguish early HCC stages from chronic liver hepatitis B and C as well as healthy controls, intensify as the disease progresses, and is highly enriched in immune function-related genes such as PD-1 (22). These all suggest that DNA methylation regulators are closely related to immunotherapy and maybe predict the response to immunotherapy. However, it is still unclear how DNA methylation regulators affect tumors, especially cancer immunotherapy in HCC. Therefore, further elucidation of DNA methylation regulators could provide an attractive perspective on cancer immunotherapy.

In this study, we integrated patients from TCGA-LIHC cohort and ICGC LIRI-JP cohort to comprehensively evaluate the correlation between the DNA methylation modification patterns (DMRPs) and tumor immune landscape. First, we explored the expression of 20 DNA methylation regulators between normal and carcinoma tissues, and then identified 3 distinct DMRPs which were tightly correlated with immune cells infiltration and prognosis. In addition, we investigated the functional annotation to distinguish cancer associated signaling pathways to the three patterns. Moreover, we continued to quantify the DMRPs of individual HCC patients and assessed the clinical responses to immunotherapy based on DNA methylation regulators-related genes score (DMRegs_score). In conclusion, our novel DMRegs_score provides a reliable insight by which to identify and feature immune landscape of HCC, and the results suggest the DMRegs_score may be a biomarker for survival and precision treatment.

Materials and method

Data collection and preprocessing of public database

HCC patients with RNA-seq, genetic mutations (VarScan) and clinical information (included age, sex, TNM stage, follow-up time, and survival status) were obtained from the Cancer Genome Atlas (TCGA) data portal (TCGA-LIHC cohort, n=374) (http://portal.gdc.cancer.gov/). The copy number variant profiles (CNV) were downloaded from the UCSC xena (http://xenabrowser.net). The normalized data from another HCC cohort were downloaded from the International Cancer Genome Consortium (ICGC) database (ICGC LIRI-JP, n=231) (http://daco.icgc.org). GSE76427 (n=167) array was downloaded from Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/). The 20 DNA methylation regulators were extracted from previous study (8). The TCGA and ICGC RNA sequence data (fragments per kilobase million, FPKM value) were transformed into TPM (transcripts per kilobase million) format. We excluded patients without complete clinical information and the survival time of 0, thus, a total of 685 HCC patients were further analyzed in this study. These detailed clinical information about 685 patients with HCC was presented in .

Tissue samples and real-time PCR and immunohistochemical staining

Forty-one pairs HCC and nearby non-tumor tissues were collected from HCC patients who underwent hepatic resection in Sun Yat-Sen Memorial hospital between Nov 2020 and Mar 2021. Liquid nitrogen was used to store these samples until further analysis could be completed. The patients’ clinical data were also collected ( ). The study protocol was approved by the Ethics Committee of Sun Yat-Sen Memorial hospital and informed consent was obtained from each patient. We extracted the RNA from the tissues with Trizol (Takara, China), and performed reverse transcription using Prime Script RTase (Takara, China), according to the manufacturer’s protocol, respectively. According to the manufacturer’s instructions, real-time PCR was used to measure mRNA expression levels using SYBR green (Takara, China). A list of the primers used for real-time PCR is provided for . Immunohistochemical (IHC) staining was performed as described previously (23) using the following antibodies: Anti- CDCA3, Anti-CDC20, Anti-YWHAQ, Anti-ADH4, Anti-TRNP1, Anti-CYP2C9, Anti-CALU, Anti-APOC1. All antibodies used in the study are shown in . Quantitative evaluation of protein expression of IHC tissues was measured by ImageJ software. The number of stained cells was identified by trainable Weka segmentation.

Interaction among DNA methylation regulators, copy number variant (CNV) analysis and gene mutation analysis

The crosstalk network diagram of multiple DNA methylation regulators was constructed by using “igraph” package, and presented the categories “Writers”, “Erasers” and “Readers” of these genes. The “RCircos” R package was used to visualize the location of 20 DNA methylation regulators in human chromosomes and the gain or loss status of copy number. The “maftools” R package was applied to evaluate the mutation status of 20 DNA methylation regulators and drawn the waterfall plots in HCC.

Molecular subgroups-based clustering analysis for DNA methylation regulators

We performed the consensus clustering with Euclidean squared distance metric and the K-means clustering algorithm to identify distinct DMRPs based on the expression of 20 DNA methylation regulators by using the “ConsensusClusterPlus” R package. HCC samples were classified into k clusters with k=2 to k=9. Based on the consistent cumulative distribution function (CDF) and delta region graphs, an optimal number of clusters was determined (24). What’s more, we compared the relationships between molecular subgroups, clinical characteristics, and prognosis using the “survival” and “survminer” R packages. The clinicopathologies including age, gender, TNM stage.

Function annotation based on gene set variant analysis (GSVA)

To identify the difference between the biological process of DMRPs, the “GSVA” package in R was utilized to performe GSVA enrichment analysis. GSVA, a nonparametric and unsupervised algorithm, can quantify the gene enrichment results in the sample of a gene expression dataset (25). In addition, we employed the “limma” R package to screen the significant variance in KEGG pathways and Hallmark pathways. The well-defined gene sets of “h.all.v7.4.symbols”, and “c2.cp.kegg.v7.2.symbols” were downloaded from MSigDB database.

Identification of differentially expressed genes (DEGs) between DNA methylation regulators modification patterns

The previous consensus clustering analysis had classified patients into three distinct DMRPs based on 20 DNA methylation regulators, and we identified DNA methylation regulators modification-related differentially expressed genes (DMRegs) among different DMRPs. The Bayesian method of “limma” package was used to statistical analysis, and “venndiagram” R package was applied to visualize the DMRegs. The DMRegs with adjusted p<0.05 and |logFC|=0.5 were considered as screening criterion. To explore the potential functions of DMRegs, the “clusterprofler” package in R software was utilized for Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis with adjusted P values < 0.05. What’s more, we further explored the gene cluster based on the expression profiles of DMRegs using unsupervised clustering methods.

Construction of the DNA methylation regulators-related gene signature

We constructed a set of scoring system to quantify the DMRPs of individual patient with HCC by using the method of LASSO cox regression, and we termed the score as DMRegs_score. The DMRegs_score was developed as follows. Univariate cox regression analysis was performed to identify overlapping DMRegs related to survival with P-values <0.05. Then “glmnet” R package was employed to establish the DMRegs_score based on the expression of significant prognosis of DMRegs among gene clusters. Finally, the DMRegs_score was defined using a formula method like previous study: DMRegs_score = Σ (Expi * coefi), where Coefi and Expi represented the risk coefficient and expression of each gene, respectively. The Kaplan-Meier survival curve, the area under the curve (AUC) of the time-dependent receiver operating characteristics (ROC) curve were implemented to evaluate the predictive ability of the risk model. Combining clinical data with univariate and multivariate cox analysis were done to determine if the risk score was an independent feature.

Estimation the relationship between DMRegs_score and TME, PD-1, PD-L1, and CTLA4

The ESTIMATE algorithm was used to calculate the tumor microenvironment (TME) scores, including stromal scores, immune scores, and estimate scores, which represented the infiltration of immune cells and stromal cells in TME (26). Moreover, based on the transcriptome profiles, we used “GSVA” R package to perform single sample gene set enrichment analysis (ssGSEA) to quantify the relative abundance of 23 immune cell types in the TME among different DMRPs (27). The marker genes of 23 immune cell types were acquired from a previous study, including activated B cell, MDSC, macrophage, regulatory T cell and so on ( ). Furthermore, we analyzed the relationship between the DMRegs_score and the expression of PD-1, PD-L1, CTLA4, and antigen presentation (HLA family).

Characteristics of mutation burden and drug sensitivity analysis

By using “maftools” R package, the mutation annotation format (MAF) from the TCGA database was used to explore whether the mutations status of HCC patients was associated with high- or low-risk group. The tumor mutation burden score (TMB) was also calculated for each patient with HCC in both groups. Using the “pRRophetic” package, we analyzed the IC50 of several chemotherapeutic drugs which were commonly used to treat HCC in both groups.

Validation the DMRegs_score and establishment of a nomogram assessing system

The reliability and predictive ability of this DMRegs_score was validated based on data from GEO dataset using same methods above mentioned. Furthermore, to expand the role of DMRegs_score in clinical practice, we used “rms” package to develop a nomogram predicting the prognosis of HCC patients, which combined the clinical features and DMRegs_score. The time-dependent ROC curves and calibration were performed to describe the predictive value of 1-, 2- and 3-year, respectively.

Statistics analysis

The statistical analysis tools-R software (version 4.0.3, R Foundation for Statistical Computing, Vienna, Austria) was used in this study. Kruskal-Wallis tests or one-way ANOVA were used as nonparametric or parametric methods for comparisons of three groups, respectively using GraphPad Prism 8. And the results of RT-qPCR and IHC were conducted statistical analysis using pair t test. The forest plot and partial violin plots were generated by Sanger Box online tool. The hazard ratio (HR) and 95% confidence intervals (CI) were calculated. All statistical results with a P-value of <0.05 were considered significant.

Results

The landscape of DNA methylation regulators in HCC

A total of 20 DNA methylation regulators were identified in this study, including three writers (DNMT1, DNMT3A, and DNMT3B), three erasers (TET1, TET2 and TET3) and fourteen readers (MBD1, MBD2, MBD3, MBD4, ZBTB33, ZBTB38, ZBTB4, UHRF1, UHRF2, MECP2, UNG, TDG, NTHL1 and SMUG1). Based on a summary analysis of the incidence of somatic mutations in these 20 DNA methylation regulators, there was a low mutation rate for the patients with HCC from TCGA cohort ( ). Thirty-four patients can be found genetic mutation in available samples, and the mutation frequency of 20 DNA methylation regulators range from 1%-2%. The TET1 had the highest number of mutations of all the DNA methylation regulators (2%). The exploration of copy number variation (CNV) alteration frequency indicated common CNV alteration in the 20 DNA methylation regulators. DNMT1, UHRF1, TET2, MBD1/2/3 and ZBTB4 were focused on copy number deletion, while MECP2, DNMT3A, ZBTB33/38, NTHL1, UHRF2 and UNG had widespread frequency of CNV amplification ( ). In addition, the locations of the CNV alterations in the 20 DNA methylation regulators on chromosomes were presented in . The status of CNV alterations indicated that CNV might regulate the mRNA expression of DNA methylation regulators. Further analysis revealed that 19 out of 20 DNA methylation regulators were upregulated in tumor samples except TET2, although the expression of TET2 in tumor tissues is higher than in normal liver tissues ( ). These data suggest the 20 DNA methylation regulators may play important roles in HCC development.

Figure 1

Landscape of genetic variations and transcriptional expression of DNA methylation regulators in HCC from TCGA cohort. (A) 34 of the 364 patients occurred genetic variations of 20 DNA methylation regulators with 9.34% mutation frequencies. (B) The frequency of CNV gain or loss among DNA methylation regulators. (C) The locations of CNV variations in DMRs on 23 chromosomes. (D) The expression level of 20 DNA methylation regulators between normal and HCC tissues. HCC, hepatocellular carcinoma; TCGA, The Cancer Genome Atlas; CNV, copy number variant.

Identification of DNA methylation regulators-related modification patterns in HCC

From above results, we speculated that DNA methylation regulators regulated deep-seated regulatory mechanism, which promotes us to further to investigate their potential functional in HCC. First, we gathered 570 patients from two HCC cohorts (TCGA-HCC and ICGC-LIRC) to explore the expression patterns of DNA methylation regulators involved in tumorigenesis. Spearman correlation analysis was utilized to assess mutual regulation among these DNA methylation regulators ( ). The results revealed MBD2, ZBTB33 and TET2 had a significant positive correlation with other DNA methylation regulators. Next, Cox regression and Kaplan-Meier analysis were performed to classify the prognostic relationship of these regulators with the HCC patients. Forest plot revealed that DNMT1/3A//3B, TET1/3, MBD1/2/3, TDG, UHRF1, SMUG1 and UNG were significantly associated with shorter overall survival and were considered as risk factors in HCC patients ( ). The crosstalk network showed the interaction and the prognostic value among 20 DNA methylation regulators in patients with HCC ( ). In addition, to better understand the role of DNA methylation regulators in tumor immunity, we explored the correlation between the 20 DNA methylation regulators and TME-infiltration immune cells using Spearman correlation analysis, and we found a significant negative relationship between most of these regulators and immune cells interaction ( ). Among them, DNMT1, ZBTB4 and MBD2 presented a strong positive correlation with most types of immune cells, such as activated CD4+ T cells, immature dendritic cells, and regulatory T cells. These results revealed that DNA methylation regulators regulated tumor microenvironment, which might provide strategies for immunotherapy. To understand the heterogeneity of DNA methylation regulators in HCC patients, we then performed unsupervised clustering analysis to classify patients based on the expression profiles of 20 DNA methylation regulators ( ). These results indicated that k=3 could achieve the best cluster efficacy. Therefore, the patients were categorized into three different DMRPs, including pattern A (n=199), pattern B (n=206) and pattern C (n=165) ( ). The Kaplan-Meier curves revealed that the pattern C had the poorer prognosis than pattern A and pattern B ( ). PCA analysis suggested three clusters were apparently discernible dimensions in the 20 DNA methylation regulators transcription profiles ( ).

Figure 2

Prognosis of DNA methylation regulators and patterns of DNA methylation regulators modification. (A) The correlation of 20 DNA methylation regulators in HCC patients. (B) The prognosis of 20 DNA methylation regulators in HCC patients. (C) The interaction among 20 DNA methylation regulators in HCC. The pink and blue line represents positive and negative correlation. The size of the circle represents the p value of the log-rank test. Green points represent favorable factors for OS. Purple points represent risk factors for OS. (D) The correlation between 20 DNA methylation regulators and 23 types of immune cells. (E) The consensus cluster matrix for patients with HCC. (F) The survival analysis for different patterns of patients. (G) PCA analysis indicated significant separation among three patterns. * p < 0.05.

Distinct DMRPs and function pathways analysis

What’s more, we found most of these DNA methylation regulators were highly expressed in pattern B. Pattern A was presented high expression levels of ZBTB4/33/38. Only NTHL1 was relatively highly expressed in pattern C ( ). These data indicated that three DMRPs had distinct characteristics in the DNA methylation regulators modification. Furthermore, as shown in the heatmap, we explored the association between the various clinicopathological characteristics and three patterns based on 20 DNA methylation regulators expression of the metadata set ( ). And we found pattern B was related to female patients and patients younger than 60 years old (p<0.05), and the number of deaths were higher than the other two patterns (p<0.001). Hence, the patients in pattern B had poorer prognosis than other two patterns. The comprehensive comparisons of the clinical features of the three DMRPs suggested most of these DNA methylation regulators played potential roles on oncogenesis. To investigate the underlying molecular mechanism and signal pathways to each DMRP, the GSVA enrichment analysis based on KEGG and Hallmark gene sets were conducted. As presented in , the results indicated significant difference between three patterns. The pattern B was significantly enriched in cell cycle and cancer-associated pathways, including DNA replication, PI3K/AKT/mTOR signaling. The pattern A was highly enriched in processes of metabolism and some carcinogenic activation pathways, such as retinol metabolism, Wnt pathway, mTOR pathway, and TGF-β signaling pathway ( and , ). However, the pattern C mainly presented enrichment metabolism in tyrosine and drug ( ). Thus, our results identified each DMRP is associated with its specific clinicopathological features and signaling pathways. Some previous studies had reported that PI3K/mTOR signaling pathway played a critical regulatory role in the tumor microenvironment. In immunology, mTOR was becoming as a key regulator of immune responses, which played an essential regulatory role in the differentiation and function of both innate and adaptive immune cells (28).What’s more, the previous results revealed that most of DNA methylation regulators were associated with multiple immune cells. Consequently, it’s essential to explore the correlation between DMRPs and tumor immune cells infiltration.

Figure 3

The clinical features of different DNA methylation regulator-related patterns and relevant function mechanism. (A) The expression of 20 DNA methylation regulators among three DNA methylation regulator-related patterns. (B) The heatmap of the differences between the expression of 20 DNA methylation regulators and clinicopathological factors. (C, D) GSVA showed the results of KEGG and Hallmark pathways in distinct DNA methylation regulators related modification patterns, respectively. * P < 0.05, ** P < 0.01, *** P < 0.001, **** P < 0.0001.

Characteristics of the TME immune cell infiltration in distinct DMRPs

Previous studies had reported DNA methylation played a crucial role on tumor immune microenvironment (29, 30). Hence, we evaluated the relationship between three DMRPs and 23 types of immune cell subsets of every HCC sample using ssGSEA. The heatmap displayed significant differences in 23 immune cells infiltrations among these DMRPs with various clinicopathological features ( ). We found that natural killer T cells, eosinophils, gamma T cells and type 1T helper cells had a higher proportion in pattern A than pattern B and C. The infiltration level of activated B cells, activated CD8 T cells, CD56bright natural killer cells, CD56dim natural killer cells, MDSCs, macrophages, monocytes, neutrophils, and type 17 T helper cells were higher in pattern C, while activated CD4 T cells and type 2 T helper cells had significantly higher infiltration in pattern B ( ). The immune landscape stated clearly that the significant differences of the relative expression of multiple immune infiltration cells among three DMRPs. To explore the influence of DNA methylation regulators on the TME of HCC, further analysis of TME scores (immune score, stromal score and estimate score) were evaluated by using the ESTIMATE algorithm. These results showed that the stromal score was the highest in pattern A than other two patterns, but there was no statistical significance between pattern B and pattern C ( ), and the highest immune score was found in pattern C ( ). However, there weren’t any significant differences in estimate scores among the three patterns ( ). In addition, blocking therapy against immune checkpoints was believed to increase the aggressiveness of the host immune system against tumor cells. Hence, we further assessed the expression levels of PD-1, PD-L1 and CTLA4 among three DMRPs. The analysis of immune checkpoints suggested that pattern B exited the highest expression of PD-1compared to pattern A and pattern C ( ). Similarly, the pattern B had a higher expression level of CTLA4 than pattern A and pattern C ( ). We also compared the PD-L1 expression levels in different DMRPs and found a significant upregulation in pattern A ( ). Based on these results, we identified that HCC patients with specific DMRPs were associated with different immune infiltration characteristics, which might influence the development and progression of HCC. What’s more, potential immunotherapy could be selected according to the expression of immune checkpoints in patients with different DNA methylation regulators modification.

Figure 4

Characteristics of the TME immune cells infiltration in distinct DNA methylation regulators related modification patterns. (A) The heatmap of immune cell infiltration in three patterns. (B) The distribution of immune cells among three patterns. (C–E) The differences of TME score (stromal score, immune score and estimate score) among three patterns. (F–H) The expression levels of three important immune checkpoint in the three DMRPs. DMRPs: DNA methylation regulators related modification patterns; TME: tumor microenvironment. * P < 0.05, ** P < 0.01, *** P < 0.001, **** P < 0.0001. ns, no significance.

DMRPs-related DEGs and gene clusters in HCC

To investigate the potential genetic alterations and expression perturbations affected by the three DMRPs in HCC, we screened a total of 151 DMRegs from three DMRPs using “limma” R package based on the metadata set ( ). Function annotation for these genes showed that some DMRegs were significantly correlated with metabolism in biological processes ( ), while material metabolism, glycolysis, and cell cycle were mainly pathways in KEGG analysis ( ). When we used univariate cox regression analysis to explore their relationship with the OS status of the HCC patients, and 112 DMRegs with significant prognostic value were selected to further identify ( ). Based on the expression profiles of these significant genes, we performed consensus clustering analysis and obtained three genomic clusters, namely gene Cluster A-C ( and ). The heatmap displayed the distinct characteristics of three phenotypes on the expression of prognostic DMRegs, and clinical analysis showed geneCluster C tended to relate to the advance TNM stage ( ). Additionally, the survival analysis demonstrated geneCluster C had a poorer survival rate ( ). Obviously, the expressions of DNA methylation regulators were significantly different among three gene clusters in the metadata set ( ). Most of these DNA methylation regulators (15/20) were presented higher expression levels in geneCluster C, such as DNMT1/3A/3B, TET1/3, ZBTB4, MBD1/2/3, TDG, UHRF1/2, SMUG1, UNG and MECP2. However, there was no significant difference between the expression of TET2, MBD4 and three gene clusters. Above all, these results illustrated the existence of specific clusters of genes in different DMRPs, which further supported the important roles of the three DMRPs in HCC.

Figure 5

DNA modification pattern-related DEGs and gene clusters in HCC. (A, B) GO and KEGG enrichment analysis for DEGs of DMRPs. (C) Consensus cluster matrix of 570 patients for k = 3. (D) The relationship of clinical characteristics and unsupervised clustering of DEGs. (E) Kaplan-Meier curves of three gene clusters. (F) The expression levels of 20 DNA methylation regulators among three gene clusters. * P < 0.05, ** P < 0.01, *** P < 0.001, **** P < 0.0001.

Generation of DMRegs_score in HCC

In order to more understand the impact of these DNA methylation regulators on patients on HCC patients and better apply the research results to clinical practice, we constructed a DMRegs_score based on 112 prognostic DMRegs. As displayed in , we performed LASSO cox analysis to build prognostic risk score based on optimal λ. According to the results, we obtained 8 genes (CDCA3, CDC20, YWHAQ, ADH4, TRNP1, CYP2C9, CALU and APOC1) in the signature, including five high-risk genes (CDCA3, CDC20, YWHAQ, TRNP1, and CALU) and three low-risk genes (ADH4, CYP2C9, and APOC1) ( ). We therefore chose these 8 genes to establish the DMRegs_score using following: DMRegs_score = (0.1507*expression of CDCA3) + (0.1259*expression of CDC20) + (0.0408*expression of YWHAQ) + (0.0220*expression of TRNP1) + (0.0443 * expression of CALU) + (-0.0189*expression of CYP2C9) + (-0.0096 * expression of APOC1) + (-0.0055 * expression of ADH4). With an optimal survival cut-point value of 3.75, we divided the patients into high-risk group (n=139) and low-risk group (n=431) ( ). A significant worse prognosis was observed for the patients in the high-risk group compared to the low-risk group ( ). The differential expressions of eight genes between high- and low-risk group were presented in .

Figure 6

Generation of DMRegs_score. (A, B) The screen of candidate prognostic genes through LASSO cox regression analysis. (C) Survival analysis of the OS between different risk groups. (D) The predictive value of DMRs score. (E, F) The differences of DMRegs_score among DMRPs and gene clusters, respectively. (G) Sankey diagram of the DNA patterns, gene clusters, DMRegs_score and survival status. (H) The correlation between DMRegs and DMRegs_score. * P < 0.05, ** P < 0.01, *** P < 0.001, **** P < 0.0001.

Three genes (ADH4, APOC1 and CYP2C9) were highly expressed in low-risk group, and the expression levels of five genes (CDCA3, CDC20, YWHAQ, TRNP1, and CALU) were higher in high-risk group. We also explored the relationship between eight genes expression and clinicopathological factors ( ). These results revealed advanced stage and death patients had higher expression level of high-risk genes. The mutated frequency of these genes was drawn using “maftool” package, and we found the few alterations happened in eight genes based on all HCC patients from TCGA ( ). The loop graph presented the chromosome locations and the gain or loss status of CNV among these risk genes ( ). The mutated frequency and CNV results indicated that these eight genes were epigenetically regulated by DNA methylation rather than DNA mutation or genomic alteration. What’s more, the 1-, 2-, and 3-year survival rate of DMRegs_score were illustrated by AUC values of ROC curves, 0.751, 0.724 and 0.710, respectively ( ). In addition, we discovered significant differences in DMRegs_score between three DMRPs and three gene clusters ( ). The DMRegs_score was highest in pattern B, while that of pattern A was lowest. Differently, the gene cluster C had the highest DMRegs_score than the other two phenotypes. The Sankey diagram showed the distribution of patients in three DMRPs, three gene clusters, DMRegs_score and survival status ( ). In addition, the correlation strength among these genes in DMRegs_score was presented in .

Assessment of TME, immune-infiltrating cells and immune checkpoint between different risk group

In order to clearly understand the relationships between DMRegs_score, TME and immune infiltrating cells, firstly, we evaluated the differences of TME scores between high- and low-risk groups. illustrated the low-risk group had higher stromal score than high-risk group. However, no significant differences of immune score and estimate score were observed between two risk groups. As represented in , the correlation of DMRegs_score and 23 immune-infiltrating cells illustrated DMRegs_score was positively correlated with actively CD4 T cells, MDSC, immune dendritic cells, natural killer T cells and type 2 T helper cells, while was negatively correlated with eosinophil, monocytes, neutrophil, regulatory T cells, and type 1T helper cells. And the heatmap displayed the abundance of 23 types of immune cell infiltration in patients with different clinical features ( ). Notably, we found the DMRegs_score was associated with T cells, so we further explored the expression of human leukocyte antigen (HLA) related genes in different risk groups ( ). Most of the HLA-related genes presented higher expression level in high-risk group. Similarly, analysis of three important immune checkpoints revealed higher expression of PD-1 and CTLA4 in high -risk group, while the expression of PD-L1(CD274) between different risk groups was no significant difference ( ). These results suggested a strong correlation between DMRegs and TME of HCC patients. We guested that the DMRegs played crucial roles on the development of HCC through influenced the immune status of HCC patients.

Figure 7

The immune cell infiltration characteristics of DMRegs_score in HCC patients. (A) The differences of TME score in the two risk groups. (B) The relationship of DMRegs_score and 23 types of immune cells infiltration. (C) The landscape of immune cells infiltration of DMRegs_score in patients with different clinicopathological features. (D) The expression differences of HLA-related genes between two risk group. (E–G) The expression level of CD274, CTLA4 and PDCD1 in two risk groups. * P < 0.05, ** P < 0.01, *** P < 0.001, **** P < 0.0001. ns, no significance.

Evaluation of the relationship between DMRegs_score and clinical characteristics

To investigate the effect of the DMRegs_score on clinical characteristics, we performed univariate and multivariate cox regression to identify whether the DMRegs_score can be an independent predicator to predict the prognosis of HCC patients. The forest plot showed the DMRegs_score could function as an independent prognostic indicator for overall survival in the multivariate analysis ( ). The clinical heat map showed

Figure 8

The clinical features of DMRegs_score in HCC patients. (A) The forest plot of univariate and multivariate cox analysis for independent prognostic factor. (B) The correlation heatmap of two risk groups with different clinical features. (C–F). The correlation between DMRegs_score and survival status, age, gender and TNM stage. (G–J) The percentile of patients with different clinical features in two risk group. * P < 0.05, ** P < 0.01, *** P < 0.001, **** P < 0.0001.

more advanced stage patients gathered in high-risk group ( ). The DMRegs_score was significantly higher in death, female, and advanced stage patients, while no statistical difference was observed in age (<60 or >=60) ( ). What’s more, we also explored the proportion of patients with different clinical characteristics in two risk groups ( ). The Kaplan-Meier curves suggested that, whatever the age, gender and TNM stage, the patients in high-risk group had poorer survival rate than patients in low-risk group ( ).

Function enrichment, tumor mutation burden and drug sensitivity analysis

To evaluate the potential molecular mechanism of this signature, we applied GSEA based on KEGG and Hallmark gene set. The function enrichment analysis demonstrated that cell cycle, glycolysis and cancer-associated pathways were mainly enriched in high-risk group ( ). These results were consistent with previous DMRPs and gene-related phenotypes. In previous studies reported TMB played a crucial role in cancer progress and immunotherapy. Therefore, we also assessed the mutational feature between two risk groups based on TCGA-LIHC cohort. The high-risk group had a higher mutation frequency (91.46%) in 82 patients, while the alteration frequency in the low-risk group was 82.46% ( ). The top ten mutated genes in the high- risk group were TP53, CTNNB1, TTN, MUC16, LRP1B, PCLO, APOB, MUC4, RYR2, and FAT3, and the most common type of mutation was missense mutation. Moreover, the Kaplan-Meier curves showed that patients with lower level of TMB had a more favorable survival rate than high-TMB ( ). Combined TMB with the DMRegs_score to assess the prognosis of HCC patients, we found patients with both high-TMB and high DMRegs_score had the poorest prognosis ( ). We further screened target therapeutic drugs for the treatment of HCC patients to assess the sensitivity in different risk groups. Importantly, we observed that patients in high-risk group had lower IC50 value for sorafenib (VEGFR inhibitor), tipifarnib (farnesyltransferase inhibitor), A.443654 (AKT inhibitor), veliparib (PARP inhibitor), olaparib (PARP inhibitor), IPA-3 (PAK inhibitor), GSK-650394 (SGK inhibitor) and CCT018159 (Hsp90 inhibitor) ( , ), while the IC50 values of axitinib (VEGFR inhibitor), motesanib (VEGFR), CCT007093 (PPM1D inhibitor), and lesteurtinib (JAK inhibitor) were higher in patients with high DMRegs_score ( , ). In a word, DMRegs_score was significantly correlated with TMB and patients’ clinical response to targeted therapy.

Figure 9

The TMB and drug sensitivity analysis for patients with different risk score. (A, B) The GSEA analysis in different risk groups. (C, D) The characteristics of tumor genetic alterations in two risk groups. (E) The survival analysis of patients with different TMB. (F) The survival analysis among four patient groups stratified by both TMB and DMRegs_score. (G–I) The drug sensitivity analysis of patients with different risk score. TMB: tumor mutation burden.

Verification the DMRegs_score and development a nomogram to predict prognosis

To validate the reliability of the DMRegs_score, we used GSE76427 as external validation group. Patients were categorized into high- and low-risk groups. The multivariate cox analysis revealed the DMRegs_score could be an independent prognostic factor ( ). Survival analysis indicated the high DMRegs_score had bad survival rate ( ). And the ROC curves showed the DMRegs_score still had accurate AUC values, 1-year for 0.733, 2-year for 0.789 and 3-year for 0.823( ). We further established a nomogram based on the risk data of DMRegs_score and the patients’ clinical features from metadata set. The nomogram composed of DMRegs_score, gender and TNM stage ( ). The AUC value of nomogram for 1-year, 2-year, and 3-year were 0.768,0.728, and 0.757, respectively ( ). In addition, we used calibration curves to confirm this nomogram prediction model ( ). However, compared with the AUC value of DMRegs_score, the nomogram scoring system had a slightly weaker predictive ability.

Figure 10

Validation of DMRegs_score in GSE76427 set and construction of nomogram. (A) The Kaplan- Meier curves of the OS between the two groups. (B) The predictive value of DMRegs_score was presented by ROC curves. (C) Nomogram to predict the 1-, 2- and 3-year OS of HCC patients in training set. (D, E) ROC and calibration curves of Nomogram for predicting of 1-, 2-, and 3-year OS in training set. ROC: receiver operating characteristic.

Validation of the expression levels of the eight risk genes which are used for the prognostic signature

Forty-one HCC tissues and adjacent normal tissues were used to detect the mRNA and protein expression of eight genes in this risk score by qRT-PCR and IHC. As presented in , the mRNA expression level of ADH4, APOC1, and CYP2C9 were downregulated while those of CALU, CDC20, CDC3A, TRNP1, and YWHAQ were elevated in HCC tissues compared to the levels in the paired normal tissues. The results of IHC staining showed the same result as qRT-PCR, and almost genes expressed in cytoplasmic except TRNP1 expressed in nucleus ( ).

Figure 11

The mRNA expression levels of 8 DMRegs of prognostic signature in hepatocellular carcinoma tissues and corresponding normal tissues by RT-qPCR.

Figure 12

The protein expression levels of 8 DMRegs of prognostic signature in hepatocellular carcinoma tissues and corresponding normal tissues by IHC.

Discussion

DNA methylation is closely related to carcinogenesis, tumor progression and metastasis (31, 32). In addition, by impacting multiple oncogenic pathways and tumor suppressor genes, DNA methylation regulators contribute to carcinogenesis in a broad range of tissue histologist (33, 34). In our study, based on expression levels of 20 DNA methylation regulators, we identified three DMRPs, and each DMRP correlated with different prognosis and signal pathway. Furthermore, Pattern B had the poorer prognosis and enriched in cell cycle and cancer-associated pathways, while the other two patterns mainly enriched in processes of metabolism. Besides, in order to confirm the efficacy of this regulatory mechanism, we applied consensus clustering analysis and found that genes cluster C was closely associated with more advanced TNM stage and poorer prognosis. Therefore, our study elucidates that the involvement of DNA methylation regulators in tumor development from a horizontal perspective, and provides new insights into the molecular networks involved in the regulation of DNA methylation regulators.

DNA methylation regulators can also impact the activation, differentiation, and functional fate of immune cells, which serve as a surveillance system against cancer (9, 35). For example, DNA demethylase TET2 promotes melanoma progression by maintaining the immunosuppressive function of myeloid cells and enhances anti-tumor immunity by governing G-MDSCs and CD8 + T-cell numbers (36, 37). The vast majority of studies, however, focused on a single DNA methylation regulator and its effect on altering TME (38–40). The immune cells infiltration characteristics, which are mediated by multiple synergistic DNA methylation regulator, have remained poorly understood. In this study, the DMRPs and TME immune infiltrating cells were closely related to each other. In addition, the immune scores and immune infiltrating cell types were significantly different in three DMRPs. More importantly, Based on DNA methylation regulators features found in individual tumors, we developed an effective DMRegs_score model and demonstrated its predictive ability. The clinicopathological features, prognosis and stromal score of high- and low-risk score were significantly distinct. Moreover, DMRegs_score was positively correlated with actively CD4+ T cells, MDSC, immune dendritic cells, natural killer T cells and type 2 T helper cells. Previous studies have indicated that more CD4+ T cells suggest a better prognosis (41, 42). However, our study found the tumor stage was more advanced in the population with high-risk group and more CD4+ T cells, and this group had better efficacy with drugs such as sorafenib and higher expression of PD-1 and CTLA-4, suggesting that the group may have a better outcome with combination therapy despite their late staging. In addition, other cancer immunotherapies in clinical trials, including dendritic cell vaccines and oncolytic viruses are also associated with TME immune infiltrating cells, such as T cells and dendritic cells (43, 44). During tumor growth, the stromal component is powerfully constricting to the immune cells, which can both be present in the tumor capsule and throughout the tumor tissue in order to prevent them from exerting anti-tumor effects (45). This is also supported by the evidence of strong stromal activation in DMRP C, where the activation of immune cells was inhibited by the AKT/mTOR pathway.

ICIs have been found to be effective in combinatorial strategies in advanced HCC patients (46). However, there is no effective biomarker for assessing the response to ICIs therapy and the prognosis of patients with HCC (47). Our detailed analyses indicated that the DMRegs_score signature probably is a robust and reliable biomarker to assess HCC patients’ responses to ICIs and TKIs. Patients with high DMRegs_score displayed higher expression of PD-1 and CTLA-4 compared with patients with lower DMRegs_score. Previous studies indicated that higher expression of PD-1, CTLA-4 and TMB might be inclined to respond to ICIs (48, 49). Thus, we concluded that patients with high DMRegs_score, which have high TMB, high expression of PD-1, and CTLA-4 might be more suitable to ICIs. Furthermore, for advanced hepatocellular carcinoma (HCC), tyrosine kinase inhibitors (TKIs) are effective therapeutic strategies. High DMRegs_score group had lower IC50 value for some types of TKIs, such as sorafenib, tipifarnib, veliparib, olaparib, which suggested that the DMRegs_score may be predictive of TKIs and ICIs combination therapy for HCC.

In comparison to existing studies of prognostic signatures of HCC, this study has some notable advantages and limitations. First, the global DNA methylation regulators landscape was modeled in order to systematically examine the effects of DNA methylation regulators on TME in HCC patients, which have not been clarified before. Furthermore, we examined the possible role of DNA methylation regulators-related status in predicting the clinical response to immunotherapy in HCC. Our data give rules about how DNA methylation regulators influenced the multiplicity of TME. Secondly, all analyses and samples were obtained primarily based on bioinformatics analysis, and although we did some clinical validation, it is indispensable to conduct prospective studies to further validate the efficacy of DMRegs_score. Besides, a few important clinical variables like surgery and chemoradiotherapy were missing, which could have affected the prognosis of DNA methylation regulators status and the immune response.

In summary, in our comprehensive research on DNA methylation regulators, we uncovered a broad range of regulatory mechanisms in HCC through which they affected clinicopathological features, TME, and prognosis. Additionally, we investigated the therapeutic effects of DNA methylation regulators in targeted therapy and immunotherapy in HCC. Our study emphasized the important clinical implications of DNA methylation regulators and provide new insights into how to personalize immunotherapy for patients with HCC.

Data availability statement

The original contributions presented in the study are included in the article/ . Further inquiries can be directed to the corresponding authors.

Ethics statement

The studies involving human participants were reviewed and approved by the institutional review committee of the First Affiliated Hospital of Xi’an Jiaotong University, Shaanxi Province, Xi’an, China. The patients/participants provided their written informed consent to participate in this study. Written informed consent was obtained from the individual(s) for the publication of any potentially identifiable images or data included in this article.

Author contributions

All authors listed have made a substantial, direct, and intellectual contribution to the work, and approved it for publication.

Funding

This work was supported by the National Natural Science Foundation of China under Granted (number 82102801& number 81902413); and Key Research and Development Program of Shaanxi under Granted number 2019SF-044 & number.2019SF-129.

Acknowledgments

We would like to thank all teammates for contributing this work.

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.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.