Gene expression profiles were obtained from TCGA (https://portal.gdc.cancer.gov/repository/), ICGC (https://dcc.icgc.org/releases/release_28/) (14) and Gene Expression Omnibus (GEO) databases (https://www.ncbi.nlm.nih.gov/geo/) (15). For TCGA-LIHC cohort, gene expression profiles produced by the Illumina HiSeq RNA-Seq platform were downloaded from TCGA database and normalized using the variance stabilizing transformation (VST) function in ‘DESeq2’ R package (R Core Team; http://www.r-project.org) (16). Clinical data were extracted from TCGA database and MSI scores were downloaded from cBioPortal database (http://www.cbioportal.org/). For the ICGC LIRI-JP cohort, gene expression profiles and clinical information were downloaded from the ICGC database (17) and normalized using VST. Clinical information about the tumor grade was available for only 212 patients. For the GSE14520 cohort, gene expression profiles were downloaded from the GEO database and normalized using the limma package (18). The gene expression profiles of TCGA Pan-Cancer were downloaded from the UCSC database (http://xena.ucsc.edu/) (19) and Student's t-test was used to analyze EPHX2 expression between tumor and normal tissues in different types of cancer.

‘Deseq2’ R package was used to screen the DEGs between HCC and paired normal samples in TCGA-LIHC and ICGC LIRI-JP cohorts (16). Adjusted P<0.05 and |log2fold change (FC)| >1 were set as the cut-off thresholds. The results were presented in a volcano map using ‘ggplot2’ R package. The TCGA-LIHC cohort was set as the training cohort, and the ICGC LIRI-JP cohort as the validation cohort. The co-expression network of DEGs in TCGA-LIHC and ICGC LIRI-JP cohorts was constructed based on TCGA-LIHC cohort using the R package ‘WGCNA’ (6). The soft-thresholding power with a slope close to 1 and a scale-free R² close to 0.9 was selected to transform the adjacency matrix to a topological overlap matrix. The soft-thresholding power was set as 7 (scale-free R²=0.91, slope=−1.49). Cut height was set as 0.25 and the minimal module size was set as 30 for network construction and module detection. The module with the highest correlation with HCC was considered as the key module. Univariate Cox proportional hazard regression analysis of those genes whose gene significance (GS) >0.2 and module membership (MM) >0.6 was conducted to screen overall survival (OS)-associated genes. Identified genes were further used to produce the prognostic multiple-gene signature using the LASSO Cox regression with the ‘glmnet’ package in R (20,21). The following risk score formula was used: Risk score (mRNA-based classifier) = sum of coefficients × expression levels of mRNA. The median risk score was used as a cut-off value to divide the patients into high- and low-risk groups for the prognostic prediction. The prognostic gene signature was validated in the ICGC LIRI-JP cohort using the aforementioned formula.

The Oncomine database (https://www.oncomine.org) is an integrated online cancer microarray database for DNA or RNA sequence analysis (22). In the present study, transcriptional expression of EPHX2 between cancer and matched normal samples was obtained from the Oncomine database. EPHX2 levels in various cancer types were compared using Student's t-test. The cut-off P-value and FC were as follows: P-value, 0.05; FC, 1.5; gene rank, 10%; data type, mRNA. EPHX2 expression in HCC was compared among the five datasets (15,23–25) by Oncomine meta-analysis.

The Kaplan-Meier plotter database (http://kmplot.com/analysis/) is an online tool containing gene expression and clinical data, which is commonly used to evaluate the prognostic value of different genes among 21 types of cancer (26). The sources of this database are GEO, TCGA and European Genome Phenome Archive. The prognostic value of EPHX2 mRNA expression in pan-cancer (n=7,489) including 21 different types of cancer was evaluated using the Kaplan-Meier plotter database (http://kmplot.com/analysis/index.php?p=service&cancer=pancancer_rnaseq). Patient samples were split into two groups by auto select best cut-off. The log-rank P-value and hazard ratio (HR) with 95% confidence intervals (CIs) were obtained.

Spearman correlation analysis (R software 3.6.3) was carried out to identify 500 genes closely correlated with EPHX2. Functions of EPHX2 and 500 EPHX2-associated genes were investigated using GO and KEGG analyses with the ‘clusterProfiler’ R package (26). Adjusted P<0.05 was considered as significant. GO terms and KEGG pathways were presented using the R package ‘GOplot’ (27). GO analysis was based on three factors, including biological process (BP), cellular component (CC) and molecular function (MF), which could predict the functional roles of EPHX2 and the 500 related genes. KEGG analysis was used to identify the pathways associated with EPHX2 and the 500 related genes.

GSEA v4.0.3 (http://www.broad.mit.edu/gsea/) was used to analyze the association between EPHX2 expression and biological pathways (28). Pre-defined gene sets (c2.cp.kegg.v7.0.symbols.gmt) were obtained from the Molecular Signatures Database (MsigDB; http://software.broadinstitute.org/gsea/msigdb). The patients in the TCGA-LIHC or ICGC LIRI-JP cohort were sorted into high- and low-EPHX2 expression groups using median mRNA expression levels of EPHX2. False discovery rate <0.25 and nominal P<0.05 were set as the cut-off thresholds.

Two independent cohorts of patients with HCC were used in the present study. For cohort one, tissue microarrays (TMAs) containing 90 pairs of HCC and matched adjacent non-cancerous tissue samples with complete clinical and follow-up data were purchased from Shanghai Liao Ding Biotechnology Co., Ltd. For cohort two, a total of 12 paired HCC and matched adjacent non-cancerous tissue samples (distance from tumor margin, ≤3 cm; cirrhosis tissue was excluded) were obtained from patients (age range, 40–69 years; five male patients and seven female patients) who were diagnosed with HCC at the Second Affiliated Hospital of Chongqing Medical University (Chongqing, China) between January 2018 and December 2019. All patients did not undergo chemotherapy or radiotherapy before surgery. After performing surgical resection, tissue samples were immediately frozen and stored in liquid nitrogen until further use. All cases were histologically confirmed. The study protocol was approved by the Ethics Committee of the Second Affiliated Hospital of Chongqing Medical University (approval no. 2020-186). All patients provided written informed consent.

Normal human hepatocytes MIHA and liver cancer cell lines Huh7, HepG2, MHCC-97H and MHCC-97L were purchased from the American Type Culture Collection. All cell lines were authenticated by STR profiling. Cells were cultured using Dulbecco's modified Eagle medium (Gibco; Thermo Fisher Scientific, Inc.) supplemented with 10% heat-inactivated fetal bovine serum (Gibco; Thermo Fisher Scientific, Inc.), penicillin (100 U/ml) and streptomycin (100 µg/ml; Beyotime Institute of Biotechnology) and maintained at 37°C in a humidified incubator containing 5% CO₂.

Total RNA was extracted from specimens in cohort two or cells using TRIzol^® reagent (Invitrogen; Thermo Fisher Scientific, Inc.), and 1 µg RNA was reverse transcribed into cDNA using PrimeScript RT Reagent kit with gDNA Eraser (Takara Bio, Inc.) according to the manufacturer's introductions. Subsequently, qPCR was carried out using SYBR Green PCR Master Mix (Takara Bio, Inc.) according to the manufacturer's protocols. PCR amplification was performed as follows: Initial denaturation at 95°C for 10 min, followed by 35 cycles of a two-step PCR at 95°C for 14 sec and 60°C for 1 min. The reaction volume was 10 µl. Relative gene expression was normalized to GAPDH and calculated using the 2^−ΔΔCq method (29). The following primer pairs were used: EPHX2, forward, 5′-CCTTCATACCAGCAAATCCCAACA-3′ and reverse, 5′-TTCAGCCTCAGCCACTCCT-3′; GAPDH, forward, 5′-GATCATCAGCAATGCCTCCT-3′ and reverse, 5′-GAGTCCTTCCACGATACCAA-3′.

Total protein was extracted from tissue samples in cohort two or cultured cells using ice-cold radioimmunoprecipitation assay buffer (Beyotime Institute of Biotechnology) supplemented with protease and phosphatase inhibitor cocktails (Roche Diagnostics). Cell lysates were boiled at 100°C for 10 min. Protein concentration was determined using the bicinchoninic acid Protein Assay kit (Thermo Fisher Scientific, Inc.). Equal amounts (40 µg) of protein samples were separated by 10% sodium dodecyl sulfate-polyacrylamide gel electrophoresis, and further transferred to polyvinylidene difluoride membranes (MilliporeSigma). Membranes were then blocked with 5% nonfat milk in Tris-buffered saline-0.1% Tween at room temperature for 2 h and incubated with primary antibodies against EPHX2 (1:1,000; cat. no. 10833-1-AP; Proteintech Group, Inc.) and β-actin (1:5,000; cat. no. 66009-1-lg; Proteintech Group, Inc.) at 4°C overnight. Subsequently, the membrane was incubated with a horseradish peroxidase-conjugated goat anti-rabbit (1:2,000; cat. no. SA00001-2; Proteintech Group, Inc.) or goat anti-mouse (1:2,000; cat. no. SA00001-1; Proteintech Group, Inc.) secondary antibodies at room temperature for 1 h. Protein bands were visualized using an enhanced chemiluminescence (ECL) detection kit (EMD Millipore). An ECL Western blot analysis system (Bio-Rad Laboratories, Inc.) was used to evaluate the bands. The intensity of protein bands was semi-quantified using Image Lab software (Version 6.0.1; Bio-Rad Laboratories, Inc.) and normalized to β-actin.

The tissue sections on the TMAs (thickness, 4 µm) were deparaffinized using xylene and rehydrated using alcohol. Samples were then subjected to heat-induced antigen retrieval using citrate buffer (0.01 M; pH 6.4) in a pressure cooker for 5 min, and cooled to room temperature. Subsequently, the sections were treated with 3% hydrogen peroxide for 10 min at room temperature to block endogenous peroxidase activity, and incubated with goat serum (Beijing Dingguo Changsheng Biotechnology Co., Ltd.) for 1 h at room temperature to block nonspecific antibody binding and incubated with an EPHX2 antibody (1:100; cat. no. 10833-1-AP; Proteintech Group, Inc.) overnight at 4°C in a humidified chamber. After incubation with the secondary goat anti-rabbit IgG (horseradish peroxidase) antibody (1:200; cat. no. ab150077; Abcam) for 30 min at room temperature, coloration with 3,3-diaminobenzidin was performed for 10 min at room temperature. Subsequently, the samples were counterstained with hematoxylin for 2 min at room temperature, dehydrated in a gradient series of ethanol, and then mounted with neutral gum. The stained tissue slices were analyzed by two different pathologists blinded to patients' clinical characteristics with a light microscope (BX53; Olympus Corporation). The intensity of IHC staining was semi-quantified using ImageJ software (Version 1.50i; National Institutes of Health). IHC scores were calculated using the following formula: IHC score = intensity score × percentage score of stained cells. Staining intensity was scored from 0 to 3 (0, negative; 1, weak; 2, moderate; 3, strong). The percentage of positively stained cells was scored from 0 to 100. Therefore, the IHC score ranged from 0 to 300. The median IHC score was defined as the cut-off value for low and high expression.

Data were presented as the means ± standard deviation and statistical analyses were performed using R software (Version 3.6.3; http://www.r-project.org). Wilcoxon matched-pairs test or Student's t-test was performed to compare differences between two groups, and One-way ANOVA followed by Bonferroni post-hoc test was used to compare differences between multiple groups. The χ² test or Fisher's exact test was performed to determine the association between EPHX2 expression and the clinical characteristics. Kaplan-Meier survival analysis was carried out to evaluate the prognostic value of gene signature and EPHX2, and the log-rank test was performed to analyze significance. When the survival curves crossed, the two-stage procedure was used for significance analysis (30,31). Univariate and multivariate Cox proportional hazard regression analyses were also carried out to investigate the association between EPHX2 expression and OS. Spearman correlation test was performed to identify EPHX2-associated genes and assess the correlation between EPHX2 expression and microsatellite instability (MSI). Receiver operating characteristic (ROC) analysis was used to determine the sensitivity and specificity of survival prediction using the gene signature risk score. The area under the curve (AUC) served as an indicator of prognostic accuracy. All the experiments were performed at least three times. P<0.05 was used to indicate a statistically significant difference.

Overall, 4,453 and 4,257 DEGs were identified in TCGA-LIHC and ICGC LIRI-JP cohorts as presented in the volcano plots (Fig. S1A and B). A total of 2,589 DEGs in both cohorts were selected for subsequent analyses and presented in a Venn diagram (Fig. S1C). Expression data of these 2,589 DEGs from TCGA-LIHC cohort were used for WGCNA (Fig. 1). Eight clustering modules (Fig. 1C and D) were used to set the soft-threshold power as 7 (scale-free R²=0.91, slope=−1.49; Figs. 1B and S2) and cut height as 0.25. The turquoise module was closely associated with clinical traits, especially tumor grade (correlation coefficient=0.38, P=2×10⁻¹³; Fig. 1E) and was considered as the key module. Moreover, 251 associated genes from the turquoise module were selected as novel candidates for further analyses (GS>0.2 and MM>0.6). Subsequently, these genes were subjected to univariate Cox regression analysis, and 239 genes were closely correlated with the OS (Table SI).

The aforementioned genes were further applied to the LASSO Cox regression analysis and the regression coefficient of each gene was calculated. The prognostic gene signature produced the best performance when seven genes (DNASE1L3, EPHX2, ADH4, G6PD, CDCA8, KPNA2 and KIAA1841) were included (Fig. 2A). Regression coefficients of these genes were presented in Fig. 2B. Three genes (DNASE1L3, EPHX2 and ADH4) with a HR<1 were considered as protective genes, whereas four genes (G6PD, CDCA8, KPNA2 and KIAA1841) with a HR>1 were considered as risk factors (Fig. 2C). Risk score was calculated for each patient using the following formula: Risk score = (−0.0011 × DNASE1L3 expression) + (−0.0029 × EPHX2 expression) + (−0.0281 × ADH4 expression) + (0.1459 × G6PD expression) + (0.0509 × CDCA8 expression) + (0.1594 × KPNA2 expression) + (0.0623 × KIAA1841 expression).

The prognostic prediction value of this seven-gene signature was further validated. Kaplan-Meier survival analysis indicated that patients with low-risk scores exhibited better survival in the training (P<0.0001) and validation cohorts (P<0.0001) (Fig. 3A and B). ROC curve analysis was also conducted to evaluate the predictive ability of the signature. Time-dependent ROC for OS at different time points indicated that AUC values for 1-, 3- and 5-year OS were 0.778, 0.705 and 0.657 in the training cohort (Fig. 3C), and the values were 0.730, 0.774 and 0.644 in the validation cohort (Fig. 3D).

Further ROC curve analysis for OS within different groups revealed that a combination of TNM stage and risk score could improve the prognostic prediction ability in the training cohort (Fig. 3E). However, the predictive ability of risk score alone was much higher in the validation cohort (Fig. 3F). In addition, the distribution of risk scores, survival time and therapeutic outcomes in the training and validation cohorts are presented in Fig. 4. The heat map revealed upregulation of the four risk factors and downregulation of the three protective genes in the high-risk group. In conclusion, combination of this novel prognostic model with conventional TNM staging might increase prognostic prediction ability in HCC. Among these seven genes, EPHX2 was rarely reported in HCC. Therefore, EPHX2 was selected for further validation.

The mRNA expression levels of EPHX2 in tumor and normal tissues were compared in different types of cancer using Oncomine and TCGA databases. In the Oncomine database, EPHX2 expression was decreased in numerous types of cancer, including leukemia, melanoma, sarcoma, and esophageal, breast, kidney, liver, colorectal, bladder, ovarian, cervical. brain and central nervous system, head and neck, and pancreatic cancer (Fig. 5A). Moreover, EPHX2 expression in cancer and normal specimens was compared using TCGA database (Fig. 5B). Similarly, the expression levels of EPHX2 were significantly reduced in tumor samples, such as urothelial bladder carcinoma (BLCA), invasive breast carcinoma, cholangiocarcinoma, colon adenocarcinoma, head and neck squamous carcinoma (HNSC), kidney chromophobe, kidney renal clear cell carcinoma (KIRC), kidney renal papillary cell carcinoma (KIRP), LIHC, lung adenocarcinoma (LUAD), lung squamous cell carcinoma (LUSC), pheochromocytoma and paraganglioma, prostate adenocarcinoma, rectum adenocarcinoma, sarcoma, stomach adenocarcinoma (STAD), thyroid carcinoma (THCA) and uterine corpus endometrial carcinoma (UCEC). These results revealed downregulation of EPHX2 in various types of cancer, which could be associated with tumor progression.

A comparison of EPHX2 expression levels in HCC from various studies in the Oncomine database revealed downregulation of EPHX2 expression in HCC tissues compared with those in adjacent normal samples (Fig. 5C). Furthermore, EPHX2 expression was reduced in HCC specimens in TCGA-LIHC, ICGC LIRI-JP and GSE14520 cohorts (Fig. 5D-F). Furthermore, RT-qPCR and western blotting were performed to evaluate the mRNA and protein expression levels of EPHX2 in 12 paired HCC and normal samples. The expression levels of EPHX2 were significantly decreased in HCC tissues compared with those in adjacent non-cancerous tissues (Fig. 6A). As presented in Fig. 6C, the protein expression levels of EPHX2 were also markedly reduced in HCC samples. Furthermore, EPHX2 expression was detected in normal hepatocytes and HCC cells. The mRNA and protein expression levels of EPHX2 were markedly decreased in HCC cells compared with those in normal hepatocytes (Fig. 6B and D). These findings revealed the downregulation of EPHX2 in HCC.

The association between EPHX2 expression and clinicopathological traits in patients with HCC was further determined. The mRNA expression levels of EPHX2 were evaluated in different subgroups in TCGA-LIHC cohort, including with or without TP53 mutations, various tumor grades and pathological/TNM stages (Fig. 7A-D). EPHX2 expression was also examined in subgroups in the ICGC LIRI-JP cohort, including various tumor grades and TNM stages (Fig. 7E and F). In TCGA LIHC cohort, the lowest mRNA expression levels of EPHX2 were detected in patients with TP53 mutations, grade 4 and TNM stage III, whereas in the ICGC LIRI-JP cohort, the lowest expression levels were found in grade 3 and TNM stage IV (Fig. 7). Due to the limited number of stage IV patients (only four HCC patients were at stage IV), no significant difference was observed between EPHX2 expression in stage IV and any other stages in TCGA-LIHC cohort. However, the mRNA expression of EPHX2 in pathological stage III was significantly lower than that in stage I and II. These results indicated that patients with HCC with TP53 mutations and at advanced tumor grade and pathological/TNM stage exhibited lower EPHX2 expression. The present study also assessed the correlation between EPHX2 expression and microsatellite instability (MSI); the results demonstrated that increased EPHX2 expression was associated with decreased MSI in LIHC (R=−0.12; P=0.027; Fig. S3). Furthermore, clinicopathological analysis suggested that downregulation of EPHX2 was associated with advanced tumor grade and female sex in TCGA-LIHC cohort, while in the ICGC LIRI-JP cohort, this was associated with higher TNM stage and tumor grade (Table SII). These results indicated that higher EPHX2 expression was detected in patients with early-stage HCC, whereas patients at the advanced stages exhibited lower EPHX2 levels.

Kaplan-Meier plotter database was used to identify the cancer types whose prognostic values were related to EPHX2 expression in OS and recurrence-free survival (RFS). The results indicated that patients with higher EPHX2 levels exhibited better OS prognosis in eight types of cancer: LIHC (HR=0.57; 95% CI=0.4–0.81; P=0.0013), LUAD (HR=0.59; 95% CI=0.41–0.84; P=0.0034), cervical squamous cell carcinoma (CESC; HR=0.42; 95% CI=0.26–0.68; P=2.5×10⁻⁴), HNSC (HR=0.61; 95% CI=0436-0.86; P=0.0048), KIRC (HR=0.44; 95% CI=0.31–0.63; P=4.5×10-6), KIRP (HR=0.42; 95% CI=0.23–0.76; P=0.003), pancreatic ductal adenocarcinoma (PDAC; HR=0.53; 95% CI= 0.34–0.82; P=0.0041) and UCEC (HR=1.71; 95% CI=1.03–2.84; P=0.036) (Fig. S4). Moreover, eight types of cancer exhibited improved RFS prognosis when EPHX2 was highly expressed; the results indicated that patients with higher EPHX2 levels exhibited lower recurrence rates in BLCA, KIRP, LUAD, LUSC, PDAC, STAD and THCA (Fig. S5). These findings indicated that EPHX2 expression could be associated with the prognosis of different types of cancer, and it may be a putative survival predictor for patients with HCC.

The prognostic value of EPHX2 in HCC was further evaluated. Kaplan-Meier survival curves indicated that patients with higher EPHX2 levels exhibited better prognosis in TCGA-LIHC (P=1.5×10⁻⁴) and ICGC LIRI-JP (P=0.017) cohorts (Fig. 8A and B). Furthermore, Kaplan-Meier analysis of OS was performed in patients with HCC according to age, tumor grade, TNM stage and pathologic stage in the TCGA-LIHC cohort, and age and TNM stage in the ICGC LIRI-JP cohort, respectively. In TCGA-LIHC cohort, EPHX2 expression affected OS rates in the different age, tumor grade, pathological stage and TNM stage III/IV subgroups (P<0.05), but not in the TNM stage I/II subgroup (P=0.06) (Fig. 8C-F). In the ICGC LIRI-JP cohort, patients with higher EPHX2 expression exhibited better OS in the age >60 years (P=0.027) and TNM stage III/IV (P=0.0058) subgroups, but not in the age <60 years (P=0.32) and TNM stage I/II (P=0.31) subgroups (Fig. 8G and H).

The independent prognostic value of EPHX2 expression with regard to OS was determined in TCGA-LIHC and ICGC LIRI-JP cohorts. In the univariate analysis, patients with HCC with higher pathological stage, TNM stage and lower EPHX2 levels exhibited poorer OS in TCGA-LIHC cohort, and patients with lower EPHX2 expression also exhibited poorer OS in the ICGC LIRI-JP cohort (Table SIII). In the multivariate analysis, improved OS was detected in patients with HCC with higher EPHX2 levels in TCGA-LIHC cohort (Table SII; HR=0.8433, 95% CI=0.7442–0.9555, P=0.0075); whereas, no significant value was identified in ICGC LIRI-JP cohort. Taken together, the mRNA expression levels of EPHX2 were associated with prognosis in HCC, and EPHX2 expression may be a promising survival predictor for patients with HCC.

To explore the detailed functions of EPHX2, GO analysis, KEGG analysis and GSEA were carried out. Firstly, A total of 500 genes significantly correlated with EPHX2 were screened by Spearman correlation test. Then, functions of EPHX2 and the 500 associated genes were analyzed using GO and KEGG analyses in TCGA-LIHC and ICGC LIRI-JP cohorts. As presented in Fig. 9A and C, commonly enriched BPs in both cohorts were involved in metabolic reprogramming, such as ‘organic acid catabolic process’, ‘small molecule catabolic process’, ‘fatty acid metabolic process’, ‘carboxylic acid catabolic process’, ‘alpha-amino acid metabolic process’ and ‘cellular amino acid metabolic process’; CCs, such as ‘peroxisome’, ‘mitochondrial matrix’, ‘microbody part’, ‘microbody’, ‘peroxisomal matrix’ and ‘peroxisomal part’. MFs, such as ‘oxidoreductase activity, acting on CH-OH group of donors’, ‘coenzyme binding’, ‘iron ion binding’, ‘vitamin binding’ and ‘oxidoreductase activity and acting on paired donors’. In KEGG analysis, enriched pathways, including ‘peroxisome’, ‘carbon metabolism’, ‘complement and coagulation cascades’, ‘valine, leucine and isoleucine degradation’ and ‘propanoate metabolism’, were associated with the functions of EPHX2 and the 500 associated genes in HCC (Fig. 9B and D).

To identify the biological pathways associated with EPHX2, GSEA was further performed to compare the high- and low-EPHX2 expression groups in TCGA-LIHC and ICGC LIRI-JP cohorts. In general, 33 and 33 KEGG pathways were enriched in the high-EPHX2 group in TCGA-LIHC and ICGC LIRI-JP cohorts, respectively; whereas, only one significantly enriched pathway was identified in the low-EPHX2 group in TCGA-LIHC cohort (Table SIV). Elevated EPHX2 levels were significantly associated with numerous metabolic pathways in both cohorts, such as ‘PPAR signaling pathway’, ‘tyrosine metabolism, ‘propanoate metabolism’, ‘histidine metabolism’, ‘valine, leucine and isoleucine degradation’, ‘retinal metabolism’, ‘primary bile acid biosynthesis’, and ‘fatty acid metabolism’ (Fig. 9E). Furthermore, low EPHX2 expression was negatively associated with ‘olfactory transduction pathway’. These results indicated that EPHX2 was involved in catabolic processes and peroxisome metabolism in HCC, and it might be associated with metabolic reprogramming in HCC.

Finally, clinical and prognostic value of EPHX2 was validated in TMAs with complete clinical and follow-up data. Those eight patients whose missing area of tissue section was >50% were excluded. Protein expression levels of EPHX2 were determined in TMAs using IHC, and the results revealed downregulation of EPHX2 in HCC tissues (Figs. S6 and 10A). Kaplan-Meier survival analysis suggested that the high EPHX2 group exhibited better OS (P=0.0046; Fig. 10B) and lower recurrence (P=0.025; Fig. 10C). Further clinicopathological analysis in the TMAs cohort indicated that EPHX2 expression was associated with tumor encapsulation, tumor multiplicity, vascular invasion and TNM stage (Table I). In addition, univariate Cox regression analysis revealed that lower EPHX2 levels were associated with poorer OS (HR=1.367; 95% CI=1.046–2.185; P<0.01; Fig. 10D) and higher recurrence (HR=1.101; 95% CI=1.042–1.365; P=0.01; Fig. S7). In addition, multivariate Cox regression analysis suggested that EPHX2 was an independent prediction indicator for OS (HR=1.415; 95% CI =1.325–2.049; P<0.01; Fig. 10D). These results revealed the clinical and prognostic value of EPHX2 in HCC and suggested that EPHX2 could be an independent prognostic biomarker for HCC.