In the present study, HCC genetic mutation profile, transcriptomic and clinical data were all downloaded by searching the Project ID ‘TCGA-LIHC’ using The Cancer Genome Atlas (TCGA; http://portal.gdc.cancer.gov/). In the section of Cases and File Counts by Experimental Strategy, 371 cases containing both clinical information and RNA sequence data were selected and saved. Mutation data were visualized using the ‘maftools’ package (v2.12.0, bioconductor.org/packages/release/bioc/html/maftools.html) in R software (v4.0.3, http://mran.microsoft.com/snapshot/2020-12-04/bin/windows/base/). MutsigCV represents mutation significance covariants.

The ‘Limma’ package (v3.52.4, bioconductor.org/packages/release/bioc/html/limma.html) of R software was used to screen for DEGs, whereas the ‘ClusterProfiler’ package (v3.0.4, http://rdocumentation.org/packages/clusterProfiler/versions/3.0.4) was used to analyze the Kyoto Encyclopedia of Genes and Genomes pathways and the Gene Ontology (GO) functions. P<0.05 and −1>log(fold-change)>1 were considered to be the thresholds.

Univariate and multivariate Cox regression analyses were performed and the ‘forestplot’ package (v1.10.1, http://rdocumentation.org/packages/forestplot/versions/1.10.1) was used to show the P-values, hazard ratios and 95% confidence intervals of each variable. P<0.05 was considered to be the threshold. A nomogram based on the results of the multivariate Cox analysis, and a list of proposed risk factors found using the ‘rms’ package (v6.1-0, http://rdocumentation.org/packages/rms/versions/6.1-0), was then constructed to predict the 1-, 3- and 5-year OS rates. If Kaplan-Meier curves crossed over, a two-stage procedure was performed using the ‘TSHRC’ package (v0.1-6, http://cran.r-project.org/web/packages/TSHRC/index.html) of R software. Analyses of risk score, survival status and heatmaps were implemented using the ‘ggrisk’ package (v1.3, http://cran.r-project.org/web/packages/ggrisk/)of R software. The time receiver operating characteristic analysis was used to compare the predictive accuracy of alcohol dehydrogenase 4 (ADH4). All analytical methods were performed using R software (version 4.0.3) with the ‘ggplot2’ package (v3.3.3, http://rdocumentation.org/packages/ggplot2/versions/3.3.3). P<0.05 was considered to indicate a statistically significant difference.

The ‘immunedeconv’ package (v2.0.3; rdocumentation.org/packages/immunedeconv/versions/2.0.3) of R software was used to integrate the Estimating the Proportions of Immune and Cancer cells algorithm to estimate the immune cell score and immune infiltration. CD274, cytotoxic T-lymphocyte-associated protein 4 (CTLA4), hepatitis A virus cellular receptor 2 (HAVCR2), lymphocyte activating 3 (LAG3), programmed cell death protein 1 (PDCD1), programmed cell death protein 1 ligand 2 (PDCD1LG2), T cell immunoreceptor with Ig and immunoreceptor tyrosine-based inhibitory motif domains, and sialic acid binding Ig-like lectin 15 were selected to be the immune checkpoint-associated relevant transcripts, before their expression values were extracted. Analyses were implemented using R software (version 4.0.3) with ‘ggplot2’ and ‘pheatmap’ packages (v1.0.12, http://rdocumentation.org/packages/pheatmap/versions/1.0.12).

Surgically resected specimens were obtained in October 2021 from 3 patients with HCC at Renji Hospital Affiliated with Shanghai Jiao Tong University School of Medicine (Shanghai, China). The inclusion criteria included: i) Pathological diagnosis of HCC; ii) the patient had received no prior treatment; and iii) the patient accepted surgical resection without distant disease. The exclusion criteria included: i) Other types of primary liver cancer; and ii) the patient refused the surgical resection. Tissue within 2 cm from the tumor boundary is called adjacent tissue. All specimens were conserved in liquid nitrogen. The present study was approved by the Ethics Committee of Renji Hospital (Shanghai, China) and written informed consent was provided by the enrolled patients.

ADH4 protein expression was examined using an immunohistochemistry assay. Fresh specimens were fixed in 4% phosphate-buffered formalin for at least 24 h at room temperature. All HCC tissues were paraffin-embedded and cut into 4-µm-thick sections. The sections were then dewaxed with xylene and rehydrated with descending graded alcohol. Activity of endogenous catalase would be inhibited by 0.3% H₂O₂. Heat-mediated antigen retrieval of tissue sections was performed at 98°C and maintained for 20 min before the sections were allowed to cool. After blocking non-specific sites with 5% goat serum at 37°C for 30 min, the sections were incubated with anti-ADH4 primary antibody (cat. no. Ab137077; 1:500 dilution; Abcam) overnight at 4°C and then with HRP-conjugated goat anti-rabbit secondary antibody (cat. no. Ab6721; 1:1,000 dilution; Abcam) for 1 h at room temperature. Antibodies were visualized using 3,3′-diaminobenzidine and counterstained using Meyer's Hematoxylin for 2 min at room temperature. Sections were then dehydrated using an increasing gradient of alcohols before being cleared in xylene. Images were observed using a light Olympus Corporation BX51 instrument. The difference of immunostaining results were evaluated using Image J software (v1.52a, USA).

ADH4 protein expression of the resected tissues was examined using western blotting. The tissues were first treated using a lysis solution (RIPA, Radio-Immunoprecipitation Assay buffer; Thermo Fisher Scientific, USA), which contained 1% (v/v) protease inhibitors (cat. no. P8340; Merck KGaA) to extract protein samples. The protein concentration was measured and normalized using Protein Quantification kit (BCA Assay; Beyotime, China). Western blot analysis was then performed as previously described (18). The membranes were blocked in Tris-Buffer Saline Tween 20 (1×TBST) containing 5% non-fat milk at room temperature for 45 min. Then, membranes would be incubated with anti-ADH4 primary antibody (cat. no. Ab137077; 1:1,000 dilution; Abcam) overnight at 4°C and HRP-conjugated goat anti-rabbit secondary antibody (cat. no. Ab6721; 1:1,000 dilution; Abcam) for 1 h at room temperature. Enhanced chemiluminescence reagent (SuperSignal™ West Atto Ultimate Sensitivity Substrate; Thermo Fisher Scientific, Inc.) was used for immunodetection. The level of actin protein expression was measured as an internal standard (anti-β actin antibody; cat. no. Ab8227; 1:1,000 dilution; Abcam). The density of protein band was quantitatively measured using Image J software (v1.52a, USA).

Potential ICB responses were predicted using the Tumor Immune Dysfunction and Exclusion algorithm (19) and software packages ‘ggplot2’ and ‘ggpubr’ (v0.4.0, http://rdocumentation.org/packages/ggpubr/versions/0.4.0). All analytical methods were performed using R software (version 4.0.3). P<0.05 was considered to indicate a statistically significant difference.

All statistical analyses were performed using R software (version 4.0.3). Unpaired Student's t-test was used to compare any differences between the subgroups in DEGs analysis. Data from two groups were compared in cox regression and K-M curve analysis using the Wilcoxon rank-sum test, whereas more than three groups were compared in immune analysis using the Kruskal-Wallis test. Post hoc test was performed using Least Significance Difference method. Data are presented by mean and SD. χ-square test was used to compare significant differences between groups. P<0.05 was considered to indicate a statistically significant difference. Test was performed no less than three independent repeats.

The experimental plan of the present study is presented in Fig. 1. A total of 270 mutant genes in HCC were screened from TCGA. Mutant genes revealed by the screen, including TP53 (28%), Titin (25%), Catenin Beta 1 (24%), mucin 16 (16%) and piccolo presynaptic cytomatrix protein (11%), had higher mutation frequencies than the rest of screened mutant genes, which is demonstrated by the horizontal histogram (Fig. 2A). In particular, missense mutations appeared to be the most common type of TP53 mutation (Fig. 2Ba). Single-nucleotide polymorphisms showed a predominant position compared with deletions or insertions (Fig. 2Bb). C>T was found to be the most distinct mutation type (Fig. 2Bc). The number of mutations per sample is shown in Fig. 2Bd. Each color in the box diagram represents one type of mutation (Fig. 2Be). Fig. 2Bf shows the top 10 mutant genes and TP53 is ranked second. Fig. 2C presents a lollipop plot, which shows that TP53 a highly mutated gene in HCC and missense mutations is the most common type of TP53 mutation.

DEGs between the TP53 mutant and wild-type TP53 groups were next screened. In total, 81 upregulated genes and 27 downregulated genes were identified (Fig. 3A and B). Notably, ‘tyrosine metabolism’ and ‘pyruvate metabolism’ were markedly downregulated, in addition to ‘alcohol metabolic process’ and ‘regulation of immune effector process’ (Fig. 3C). By contrast, mutant TP53 was particularly enriched for pathways associated with the ‘cell cycle’ and cell metabolism, ‘cytochrome metabolism and drug metabolism (Fig. 3Ca and Cb). Furthermore, genes associated with TP53 mutation tended to be more enriched for terms associated with the metabolic process, including nuclear division and humoral immune response (Fig. 3Cc and Cd). These results suggested that genes associated with TP53 mutations are likely to serve an important role in the immunometabolism of the TME in HCC.

A total of 17 genes and corresponding clinical information, including age, pTNM_stage and grade were identified through univariate (Fig. 4A) and multivariate (Fig. 4B) Cox regression analyses, to build a predictive nomogram. The predictors included ADH4, 4-hydroxyphenylpyruvate dioxygenase, patient age and pathological Tumor-Node-Metastasis (pTNM) stage (20) (Fig. 4C), all of which met the P<0.05 criteria during risk assessment. The P-value of ADH4 was found to be <0.001, which suggested that ADH4 can be exploited as a distinct prognostic signature for HCC. The calibration plots for the 1-, 3- and 5-year OS rates were consistent compared with those predicted by the ideal model in the entire cohort (Fig. 4C).

The mutant TP53 and wild-type TP53 groups included 105 and 253 patients, respectively. Patients in the mutant TP53 group displayed significantly worse OS time compared with those in the wild-type group (Fig. 5A). To investigate if ADH4 was independent of the TP53 mutation status, patients with HCC were divided into high- and low-ADH4 groups based on the TP53 mutation status. The high- and low-ADH4 groups in the mutant TP53 group included 52 and 53 patients, respectively. By contrast, the high- and low-ADH4 groups in the wild-type TP53 group included 126 and 127 patients, respectively. The results revealed that the low-ADH4 group in both the mutant and wild-type TP53 groups displayed significantly worse OS time compared with the high-ADH4 group (Fig. 5B and C). The crossed curves in Fig. 5 were determined to be statistically different (Fig. 5A, P=0.0294; Fig. 5B, P=0.0147; Fig. 5C, P=0.0063;) by performing a two-stage procedure, using the ‘TSHRC’ package of R software. In addition, wild-type ADH4 was observed to occupy a particularly higher proportion of patients with HCC than mutant type ADH4 (χ² test; Fig. 5D).

According to the results in Fig. 5, high risk group in Fig. 6 was defined as low ADH4 expression while low risk group in Fig. 6 was defined as high ADH4 expression. To explore the distribution of immune scores (Fig. 6A) and immune infiltration (Fig. 6B) in the groups of different risks, 38 types of immune cells were compared in total. Among them, 23 types of immune cells were found to be statistically different from other rest 15 types of immune cells. It was noted that there was significant immune infiltration by cells expressing markers of B cells in the high-ADH4 group than other immune cells in the Fig. 6B, whilst the expression of markers of T cells was decreased, suggesting a differential immune TME. For the expression distribution of immune checkpoint markers in the HCC tissues, CD274, CTLA4, HAVCR2, LAG3, PDCD1, PDCD1LG2, TIGIT and SIGLEC15 were discovered to be the immune checkpoint-associated transcripts showing a statistical difference, which was different from other transcripts (Fig. 6C).

To validate the identified gene signature, ADH4 expression levels were examined using immunohistochemistry and western blotting on pairs of HCC tissues and adjacent normal tissues from 3 patients. The pathological diagnosis of the 3 male patients was HCC without distant metastasis, and none had received any prior treatment. The clinicopathological information of the 3 patients is shown in Table I. The TP53 status of these patients was mutant type and the molecular hallmark was a mutation at codon 249 in TP53. The results demonstrated that ADH4 protein expression was markedly increased in the normal tissues compared with the tumor tissues (Fig. 7A and D). The risk score of every patient whose RNA Seq files was downloaded from TCGA was then calculated to obtain the median cut-off point to divide them into the high-ADH4 group (n=185) and low-ADH4 group (n=185) (Fig. 7B, top). The survival status of all patients with HCC is shown (Fig. 7B, middle) and the prognostic ADH4 gene expression profile is shown in the heatmap (Fig. 7B, bottom). The Kaplan-Meier survival curves demonstrated that patients in the low-ADH4 group had worse OS time compared with those in the high-ADH4 group (Fig. 7C, P=0.0002). In addition, a time-dependent receiver operating characteristic analysis suggested that ADH4 expression had accurate predictive capabilities for 1-, 3- and 5-year OS (Fig. 7E). A Sankey diagram of ADH4 expression revealed the distribution of the same sample in different characteristic variables including ages, pTNM_stages, grade, ADH4 expression type and status of patients (Fig. 7F). Potential ICB responses indicated that the distribution of the TIDE scores in the low-ADH4 group showed a worse response compared with that in the high-ADH4 group (Fig. 7G).