The data with HCC were gathered by TCGA database [17], and all samples with survival information were used. The measured degrees of methylation were reported as beta values ($\beta$) [18]. CRGs were sourced from previous literature [19], and the appropriate details of DNA-methylation site were acquired from the methylation annotation file of the TCGA database. Two-thirds of the HCC patients were picked as training-cohort to find and create prognostic signature based on CRG-located DNA-methylation sites, while the remaining one-third of patients served as validation-cohort for validation of the signature.

To identify sites associated with survival, Univariate-Cox-regression [20] was performed on CRG-located DNA-methylation sites. DNA-methylation sites with statistical significance were picked (p 0.01) for multivariate-Cox-regression in order to form models incorporating all conceivable arrangements of 2–4 sites, aiming to screen for biomarker combinations associated with survival. Cox regression was conducted utilizing the “survival” (v.3.5-0) [21] R package. The methylation levels of CRG-located DNA-methylation sites and accompanying regression-coefficients were used calculated patient’s prognostic risk score. In accordance with the median, all recipients of HCC were broken down into high- and low-risk groups.

The Kaplan-Meier analysis [22] and Log-rank test were carried out to contradistinguish survival diversities between patients in low- and high-risk groups. The Receiver-operating-characteristic (ROC) curve [23] was adopted to cast the sensitivity and specificity of CRG-located DNA-methylation combinations and to judge the capacity of risk prediction models to predict HCC prognosis at 1, 3 and 5 years, which were displayed applying the “survivalROC” (v.1.0.3.1) R package. The nomogram [24] incorporating the age, clinical stage, gender, risk score and tumor grade was constructed to forecast overall survival (OS). The nomograms were created using the “rms” (v.6.5-0) R package.

The Stromal scores, ESTIMATE scores, Immune scores and immune cell infiltration levels were computed for all HCC individuals employing “CIBERSORT” [25] and “ESTIMATE” [26] R packages. The correlation between signature and immune-checkpoint-associated genes (ICGs) were analyzed using spearman correlation and the expression of ICG was contrasted between different groups utilizing the Mann-Whitney-U test. Furthermore, differentially-expressed genes (DEGs) comparing two group patients were obtained ($|$logFC$|$ $\ge$1, p 0.01, FC means the ratio of the gene expression value between two group). Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis of DEGs were carried out implementing the David [27].

The LX-2 cell lines were obtained from Guangzhou Otwobio Biotech Inc (Guangzhou, China). The WRL68 and SNU449 cell lines were obtained from Shanghai Genechem Co., Ltd (Shanghai, China). The MHCC-97L and MHCC-97H cell lines were obtained from Shanghai iCell Bioscience Inc (Shanghai, China). All cell lines have been authenticated by short tandem repeat (STR) profiling. All cell lines were tested for presence of mycoplasma using the polymerase chain reaction (PCR) Mycoplasma Test Kit (Hangzhou HuaAn biotechnology Co., Ltd, Hangzhou, China, catalogue number K0103) and found to be free of contamination. These cell lines were cultured in a high sugar medium. The medium was supplemented with 10% fetal bovine serum (FBS) and 1% P/S. The cell culture was placed in an incubator, which was maintained at 37 °C and 5% CO${}_2$ atmosphere. The cell culture was passaged every 2 days. Genomic DNA of the cell line was extracted using Animal Tissues/Cells Genomic DNA Extraction Kit (Beijing solarbio science & technology Co., Ltd., Beijing, China) and was modified by the EpiTect Fast DNA Bisulfite Kit (QIAGEN, Duesseldorf, Germany). The Light cycler 480 Real Time PCR system (Roche, Basel, Switzerland) was used in conjunction with the primers listed in Supplementary Table 1 to conduct Real-time quantitative methylation-specific PCR (RT-qMSP) to identify gene methylation. The 2${}^{-\mathrm{\Delta }\mathrm{\Delta }\text{CT}}$ method, where $\mathrm{\Delta }$CT = CT${}_{\text{M}}$ – CT${}_{\text{U}}$ (M: methylated, U: unmethylated), was implemented to measure the methylation.

A total of 96 CRGs were meticulously gathered through extensive literature searches, and 1374 DNA-methylation sites located in the CRGs were extracted. 277 CRG-located DNA-methylation sites that were substantially linked with prognosis were sought out employing univariate-Cox-regression, and then executed with multivariate-Cox-regression to create a prognostic signature. The signature comprising four CRG-located DNA-methylation sites (cg01123518, cg05337637, cg05398307, cg05706061) was selected as the best prognostic model for predicting survival through ROC analyses. The sites were linked to patient survival in both univariate- and multivariate-Cox models. These sites corresponded to the genes NUBP1, NDOR1, MPI and SLC31A2. The information of the site displays in Table 1.

The coefficients and methylation levels were used to generate the CRG-located DNA-methylation risk score for patients. Score = (17.48412 $\times$ $\beta$${}_{cg01123518}$) + (–3.71968 $\times$ $\beta$${}_{cg05337637}$) + (34.83459 $\times$ $\beta$${}_{cg05398307}$) + (–1.42065 $\times$ $\beta$${}_{cg05706061}$). In accordance to the median, participants in the training cohort were classified into high- and low-risk groups. The Kaplan-Meier analysis revealed substantial survival disparity between them and ROC analysis demonstrated that the CRG-located DNA-methylation signature had good predictive accuracy (Fig. 1A,B). The Kaplan-Meier and ROC analysis in the validation cohort revealed a considerably decreased survival odds as risk increased (Fig. 1C,D).

Fig. 1.

Kaplan-Meier and Receiver Operating Characteristic (ROC) analyses of signature in The Cancer Genome Atlas (TCGA) training (A,B) and validation cohort (C,D). HR, the estimate of the ratio of the hazard rates between the high-risk and low-risk group.

Age, gender, clinical stage, and tumor grade can also be used as predictors for determining HCC patient prognosis. All individuals were broken down into two groups determined by median, and it ultimately emerged that the proportion of patients in each age, clinical stage, gender, and tumor grade was relatively even in the different risk groups (Fig. 2A). To assess the independence and applicability of the signature, patients were regrouped via different clinical characteristics. The regrouping of patients according to age revealed no significant difference between the high-age and low-age group, demonstrating that the signature had no bearing on patients’ age. Depending Kaplan-Meier analysis, patients in the low-risk group across all age groups had a considerably greater outcome than those in the high-risk group, and the methylation signature had superior predictive performance (Fig. 2B). Grouping patients by gender revealed that female had higher risk scores, and both male and female groups’ survival rates differed drastically with different risk scores (Fig. 2C). Tumor grade is an important prognostic characteristic and categorizing patients into high-grade (G3 and G4) and low-grade (G1 and G2) cohorts found that signature effectively distinguished patients’ risk for any tumor grade group (Fig. 2D). DNA methylation changes pertained with disease stage, and patient survival outcomes may vary widely even at the same stage. Due to the limited number of samples per stage, patients were split up into late-stage (III and IV) and early-stage (I and II) cohorts. Despite the fact that patients’ disease severity varied greatly, there was an equally large survival variation between two risk groups at the same stage (Fig. 2E). All results demonstrate the validity of risk stratification, and the CRG methylation signature may provide a better reference for different regrouped cohorts, suggesting that this signature is an independent prognostic predictor applicable to patient survival.

Fig. 2.

The independence analysis of signature in different clinical and pathological groups. (A) Ratio of clinical characteristics. The low-risk group in outer circle, high-risk group in inner circle. (B–E) The distribution of risk score and Kaplan-Meier curves stratified by clinical characteristics.

The independence of CRG-located DNA-methylation signatures from characteristics such as age, clinical stage, gender and tumor grade was explored using Cox regression. The findings revealed that stage and risk score were connected to OS in univariate-Cox-regression for HCC patients (Fig. 3A) and maintained strongly correlated in multivariate-Cox-regression (Fig. 3B), indicating that it might be considered as a standalone prognostic sign. The clinical characteristics and risk score were thereafter utilized to create the prediction nomogram of OS, which can provide information for clinical purposes (Fig. 3C).

Fig. 3.

The independence of signatures from clinical characteristics (A,B) and nomogram construction for the survival prediction (C).

The distinct variations in immune cell infiltration between various risk groups were evaluated to understand the potential relevance of DNA-methylation signature to the immune status. The findings revealed that Stromal and ESTIMATE scores, which were considerably greater in the low-risk group, were linked adversely with risk score (Fig. 4A,B). Pursuant to the immune cell infiltration research, risk scores were strongly inversely linked with T-cells-CD4-memory-resting, T-cells-gamma-delta and Macrophages-M1, and favorably tied with T-cells-regulatory-Tregs and Macrophages-M0. The immune cell infiltration level of T-cells-regulatory-Tregs and Macrophages-M0 were drastically greater in the high-risk group (Fig. 4A,C). The checkpoint inhibitors are critical in clinical treatment and therefore the differential expression of ICGs between two risk groups was further regarded. The findings demonstrated an elevated relationship between risk scores and PD-1, CTLA4, CD276, CD80, CD86 (Fig. 5A), and there were considerably distinct in CTLA4, CD276, and CD80 between two groups (Fig. 5B). In addition, it was found that these genes participate in the cell adhesion molecules pathway, and 67% of 157 genes in the pathway are expression varied between normal and malignant tissues. These findings imply a potential link between that the CRG-located DNA-methylation signature and tumor immune microenvironment, and thus affect the prognosis of patients.

Fig. 4.

The association of signature with immune (A) and expression differences between different risk hepatocellular carcinoma (HCC) patients (B,C). ns: no significance, *p 0.05, **p 0.01, ***p 0.001.

Fig. 5.

The association of signature with immune-checkpoint-associated genes (ICGs) (A) and expression differences of ICGs between different risk HCC patients (B). ns: no significance, *p 0.05, **p 0.01, ***p 0.001.

DEGs between two risk groups were looked at in order to gain an awareness of the fundamental mechanisms in various risk groups. 1316 up-regulated genes and 35 down-regulated genes were found in the high-risk group. GO enrichment revealed that up-regulated genes were primarily linked to cation channel complex, ion channel activity, passive transmembrane transporter activity and voltage-gated potassium channel activity (Fig. 6A). KEGG pathway enrichment revealed that up-regulated genes were enriched in GABAergic synapse, glutamatergic synapse, neuroactive ligand-receptor interaction and protein digestion and absorption (Fig. 6B).

Fig. 6.

The Gene ontology (GO) (A) and Kyoto Encyclopedia of Genes and Genomes (KEGG) (B) analysis of differentially-expressed genes (DEGs).

DNA methylation alterations can influence the onset, development, and treatment of various tumors by regulating the transcription of corresponding genes. Our analysis found that cg05706061 contained in prognosis signature was located in the promoter region of the cuproptosis-related gene SLC31A2 (Fig. 7A). The DNA methylation level of cg05706061 was significantly different between tumor tissue and normal tissue (Supplementary Fig. 1), and significantly correlated with the expression of SLC31A2 (Supplementary Fig. 2). We further investigated the promoter methylation status of SLC31A2 by qMSP, the result showed that the DNA-methylation level of SLC31A2 in HCC cell lines (SNU449, MHCC-97L and MHCC-97H) was significantly decreased compared to that in normal liver cell lines (WRL68 and LX-2 cell) (p 0.01, Fig. 7B, Supplementary Fig. 3). Furthermore, the DNA-methylation level of SLC31A2 exhibited a gradual decrease with increasing metastasis potential of the cell lines. The DNA-methylation levels in primary liver cancer cell line (SNU449) were significantly higher than in low metastatic cell lines (MHCC-97L) (p 0.01), which were in turn significantly higher than that in high metastatic cell line (MHCC-97H) (p 0.05) (Fig. 7B, Supplementary Fig. 3).

Fig. 7.

The methylation level of SLC31A2 in HCC cell lines. (A) Genomic organization of SLC31A2 and the qMSP assay. The X-axis represents chromosome location (GRCh38). qMSP represents qMSP section location. cg05706061 represents cg05706061 site location in the promoter region of SLC31A2. (B) The DNA methylation level of SLC31A2 was determined by qMSP in HCC cell lines (SNU449, MHCC-97L and MHCC-97H) and control liver cell lines (WRL68 and LX2 cell). ns: no significance, *p 0.05, **p 0.01, ***p 0.001.