Ubiquitin-related genes were acquired from the integrated annotations for Ubiquitin and Ubiquitin-like Conjugation Database (iUUCD, http://iuucd.biocuckoo.org/browse.php), the Molecular Signatures Database (MSigDB, https://www.gsea-msigdb.org/gsea/msigdb/) and the human E3 ubiquitin ligase database (https://esbl.nhlbi.nih.gov/Databases/KSBP2/Targets/Lists/E3-ligases/) [16]. The identified genes were presented in Supplementary Table 1. RNA sequencing data and corresponding clinical information from 374 liver tumor tissues and 50 adjacent tissues were downloaded from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/). The “limma” package in R software was used to identify differentially expressed genes (DEGs) in tumors and adjacent tissues among the ubiquitin-related genes. The screening criteria were false discovery rate (FDR) <0.05 and |log (fold change)| ≥1.0.

The “caret” package in R software was used to randomly divide 370 HCC samples with complete clinical information into two groups, a training cohort (n=224) and a testing cohort (n=146). In the training cohort, the “survival”, “glmnet”, and “survminer” packages in R software were used to construct a ubiquitin-related gene prognostic signature by successively performing univariate Cox, least absolute shrinkage and selection operator (LASSO), and multivariate Cox regression analyses. The risk score of each patient with HCC in the training and testing cohorts was calculated as the sum of the expression values of each ubiquitin-related gene in the model multiplied by its regression coefficient. The median was used as the threshold to divide the patients into high- and low-risk groups. The “survival” and “timeROC” packages in R software were used to analyze Kaplan–Meier survival curve and receiver operating characteristic (ROC) curve for high- and low-risk groups of patients with HCC, and the area under the curve (AUC) was used to evaluate the accuracy and specificity of the risk prediction model.

All data were analyzed using GraphPad Prism software (version 8.0) or R software (version 4.1.0) and presented as the mean±standard deviation. Cox regression analyses were performed to evaluate the effects of ubiquitin-related genes on the survival of patients with HCC. Student’s ttests were used to test the differences between two groups. One-way ANOVA was used to calculate differences among multiple groups. The log-rank test was used to assess the statistical significance of the Kaplan–Meier survival curve. Two-way ANOVA was used to calculate the statistical significance of dynamic curves of tumor volume in mice treated with different compounds. P<0.05 was considered statistically significant.

Based on the iUUCD (including E1s, E2s, E3s, and DUBs), MSigDB, and human E3 ubiquitin ligase databases [16], we collected a ubiquitin-related gene set consisting of 1,423 genes (Fig. 1A). By differential expression analysis of this gene set in the TCGA database, a total of 92 DEGs were identified in HCC and adjacent tissues, of which 91 were upregulated and one was downregulated (Fig. 1B). To explore the clinical significance of these DEGs, 370 patients with HCC that had complete clinical information in the TCGA database were randomly divided into training (n=224) and testing (n=146) cohorts. Next, we conducted a univariate Cox regression analysis in the training cohort and found that 28 DEGs were significantly positively correlated with poor prognosis in patients with HCC (Fig. 1C). Furthermore, we performed LASSO regression and multivariate Cox regression analyses and identified four core ubiquitin-related genes (SSR3, CYHR1, UBE2S, and SKP2) that could be used to construct a risk model for HCC survival prediction (Fig. 1D, E). The regression coefficients of the four genes in the prognostic signature were 0.0248, 0.0783, 0.0613, and 0.0592, respectively. The formula for calculating the risk score for each patient with HCC is the sum of the expression values of these four genes multiplied by their respective regression coefficients. Thus, a risk prediction model for HCC based on ubiquitinrelated genes was successfully constructed.

To evaluate the clinical significance of this risk-prediction model, we calculated the risk score of each patient with HCC in the training cohort and divided them into high- and low-risk groups based on the median score (Fig. 1F, G). In particular, the survival status of patients gradually decreased as the risk score increased (Fig. 1H). Moreover, Kaplan–Meier survival analysis showed that the overall survival of patients with HCC in the high-risk group was significantly lower than that of patients in the low-risk group (Fig. 1I). ROC curve analysis showed that the 1-year, 2-year, and 3-year survival AUC values of the risk prediction model were 0.76, 0.72, and 0.68, respectively (Fig. 1J). In addition, we calculated the risk score of patients with HCC in the testing cohort using the same formula and found a clinical significance consistent with that of the training cohort (Supplementary Fig. 1A–E).

To further examine the potential of the risk score derived from the model as an independent prognostic factor in patients with HCC, we performed univariate and multivariate Cox regression analyses on five variables (age, sex, tumor grade, tumor stage, and risk score) for the entire TCGA cohort. We found a significant association between tumor stage and risk score with the overall survival of patients with HCC (Table 1). The results indicate that the risk score can be used as an independent prognostic factor to evaluate the overall survival of patients.

Next, we examined the expression characteristics and clinical significance of each of the four genes in the risk model. The results showed that only CYHR1 and UBE2S, but not SSR3 and SKP2, had significantly higher RNA and protein levels in HCC tissues than in adjacent normal tissues (Supplementary Fig. 2A). In addition, we evaluated the correlation between the expression levels of the four genes in the model and the survival of patients with HCC. We found a significant correlation of only the expression level of UBE2S with the overall survival and disease-free survival of patients with HCC, wherein the survival of patients with high expression of UBE2S was significantly reduced compared to that of the low expression group (Supplementary Fig. 2B). Therefore, in subsequent experiments, we focused on exploring the biological functions and mechanisms of UBE2S in HCC.

Since the expression level of UBE2S in HCC tissues is significantly higher than that in normal tissues (Supplementary Fig. 2A), we therefore explored the factors regulating UBE2S expression in HCC. Based on JASPAR database analysis, we identified four transcription factors, BATF3, FOXD3, E2F2, and POU2F2, which might bind to the promoter region of UBE2S (Fig. 2A). Further, we found that only the expression level of E2F2, but not the other three transcription factors, in HCC tissues was significantly higher than that in normal tissues (Supplementary Fig. 3A). In particular, the correlation coefficient between the RNA expression levels of UBE2S and E2F2 in HCC tissues was significantly higher than that of the other three transcription factors (Supplementary Fig. 3B). Moreover, we found a significant positive correlation between the protein expression levels of E2F2 and UBE2S by immunohistochemistry in a microarray containing 90 HCC tissue samples (Fig. 2B). qPCR and Western blotting demonstrated that E2F2 overexpression significantly upregulated the RNA and protein levels of UBE2S (Fig. 2C, D). Furthermore, we constructed a luciferase reporter containing the UBE2S promoter to explore whether the transcription of UBE2S is regulated by E2F2. As shown in Figure 2E, E2F2 overexpression increased UBE2S promoter activity. Next, luciferase gene reporters containing truncated domains of the UBE2S promoter were constructed. When the UBE2S promoter domain was truncated to –250 bp/+199 bp, the increase in UBE2S promoter activity mediated by E2F2 overexpression was abolished, indicating that E2F2 could bind to the domain (–500 bp/–250 bp) of the UBE2S promoter (Fig. 2F). To further confirm the binding sites of E2F2 in the UBE2S promoter, we predicted four potential binding sites using the JASPAR database and constructed luciferase gene reporters containing the UBE2S promoter with E2F2 binding site mutations. Mutations in the first binding site of E2F2 did not impede the increase in UBE2S promoter activity mediated by E2F2 overexpression, whereas the other three mutations significantly alleviated the increase in activity, suggesting that E2F2 could bind to binding sites 2, 3, and 4 of the UBE2S promoter (Fig. 2G). In addition, the binding of E2F2 to the UBE2S promoter was confirmed by chromatin immunoprecipitation (ChIP) and qPCR assays (Fig. 2H). Taken together, these results indicate that E2F2 directly binds to the UBE2S promoter and transcriptionally upregulates UBE2S expression.

To systematically elucidate the downstream functions of UBE2S, we constructed Huh-7 cells with UBE2S-knockout using CRISPR-Cas9 and performed RNA sequencing (Supplementary Fig. 4A). A total of 2,458 genes were upregulated, and 2,558 genes were down-regulated after the knockout of UBE2S (Fig. 3A). Further, KEGG analysis revealed a significant relationship between the DEGs and the metabolic pathways (Fig. 3B). Next, a non-target metabolomic analysis was performed to identify specific metabolic pathways regulated by UBE2S (Supplementary Fig. 4B). Differential abundance score analysis revealed a downregulation of all identified metabolites in five pathways, namely glycerolipid metabolism, pentose and glucuronate interconversion, glycolysis/gluconeogenesis, pentose phosphate pathway, and phenylalanine, tyrosine and tryptophan biosynthesis, in Huh-7 cells with UBE2S knockdown (Fig. 3C). Previous studies have shown that glycolysis interacts with several metabolic pathways such as the pentose phosphate pathway, lipid metabolism, and amino acid metabolism and thus plays a central role in the metabolic network [17]. Therefore, we focused on the regulatory role of UBE2S in glycolysis in subsequent experiments.

Cluster analysis showed that the levels of glucose-6-phosphate (G-6-P), fructose-1,6-diphosphate (F-1,6-BP), dihydroxyacetone phosphate (DHAP), 3-phosphoglycerate (3-PG), 2-phosphoglycerate (2-PG), and phosphoenolpyruvate (PEP) in the glycolysis pathway were significantly downregulated following UBE2S knockdown (Fig. 3D, E). Further, we sequentially treated HCC cells with glucose, oligomycin (an inhibitor of ATP synthase), and 2-deoxy-Dglucose (2-DG, an inhibitor of glycolysis) and performed the Seahorse XF glycolysis stress test. The dynamic curve of extracellular acidification rate (ECAR) is shown in Figure 3F and Supplementary Figure 4C. UBE2S knockdown suppressed the glycolysis and glycolytic capacity (Fig. 3F), whereas UBE2S overexpression enhanced these abilities (Supplementary Fig. 4C). Moreover, glucose uptake and lactate production were reduced in Huh-7 and HCCLM3 cells by UBE2S silencing (Fig. 3G, H) and enhanced in MHCC-97H cells following UBE2S overexpression (Supplementary Fig. 4D, E). Crucially, CCK8 assays revealed that the cell proliferation enhanced by UEB2S overexpression in different HCC cell lines was significantly inhibited by 2-DG (Fig. 3I). Taken together, these results suggest that UBE2S promotes the proliferation of HCC cells by positively regulating glycolysis

To elucidate the molecular mechanism by which UBE2S promotes glycolysis in HCC cells, we focused on changes in the metabolic enzymes involved in this process. The cluster heat map analysis revealed significant down-regulation of the RNA levels of a series of enzymes involved in glycolysis, such as GLUT1, HK2, PFKL, and LDHA, in Huh-7 cells with UBE2S knockout (Fig. 4A). These results were further verified using qPCR and Western blotting (Fig. 4B, C, Supplementary Fig. 5A). To investigate the mechanisms underlying the UBE2S regulation of these glycolytic enzymes, we performed gene set enrichment analysis (GSEA) and found the inhibition of several hypoxia-related signaling pathways in Huh-7 cells with UBE2S knockout (Fig. 4D). Moreover, KEGG pathway analysis also revealed a significant relationship between the DEGs and the HIF-1 signaling pathway (Fig. 3B). Furthermore, we found that UBE2S silencing clearly decreased the protein level of HIF-1α in Huh-7 cells, whereas UBE2S overexpression remarkably enhanced HIF-1α protein expression under both normal and hypoxic conditions (Fig. 4E). Since HIF-1α is a key transcription factor regulating cellular glycolysis [18], these results suggest that UBE2S may promote the glycolysis of HCC cells by upregulating the expression of HIF-1α-mediated metabolic enzymes.

To further clarify whether the promotion of glycolysis by UBE2S is dependent on HIF-1α, the Seahorse XF glycolysis stress test was performed. Although UBE2S overexpression in HCCLM3 and MHCC-97H cells increased glycolysis by converting glucose into pyruvate, these functions were significantly inhibited after HIF-1α knockdown (Fig. 4F). Moreover, the levels of key enzymes involved in glycolysis, such as GLUT1, HK2, PFKL, and LDHA, were significantly increased when UBE2S was upregulated in HCCLM3 and MHCC-97H cells, which was decreased with HIF-1α deletion (Fig. 4G). In addition, the increase of glucose uptake and lactate production caused by UBE2S overexpression in HCCLM3 and MHCC-97H cells was also significantly inhibited by PX-478, a HIF-1α inhibitor (Supplementary Fig. 5B, C). Furthermore, CCK8 assays revealed that cell proliferation enhanced by UBE2S overexpression could be significantly inhibited in different HCC cell lines after HIF-1α knockdown or PX-478 treatment (Fig. 4H). Taken together, these results suggest that UBE2S promotes HCC cell proliferation by upregulating HIF-1α-pathway dependent glycolysis.

Given that UBE2S knockout did not alter the RNA level of HIF-1α (Supplementary Fig. 6A), we therefore focused on exploring the regulatory role of UBE2S on HIF-1α at the protein level. Von Hippel-Lindau tumor suppressor (VHL) is a classical ubiquitin ligase that mediates the ubiquitination of HIF-1α, which can destroy the protein stability of HIF-1α by causing its degradation [19]. Interestingly, we found that UBE2S overexpression in different HCC cell lines significantly reduced the expression of VHL and increased the stability of the HIF-1α protein (Supplementary Fig. 6B). Crucially, rescue experiments showed that VHL overexpression could significantly inhibit the increased protein stability of HIF-1α caused by UBE2S overexpression (Supplementary Fig. 6C). These results suggest that UBE2S improves the protein stability of HIF-1α in a VHL-dependent manner. Next, we found that the negative regulatory effect of UBE2S on VHL was significantly inhibited when the cells were treated with MG132, the proteasome inhibitor (Supplementary Fig. 6B). In addition, we also constructed a plasmid encoding catalytically inactive UBE2S in which cysteine was mutated to alanine at site 95 [20]. Mutant UBE2S alleviated the inhibitory effect of wild-type UBE2S on VHL protein expression and could not upregulate the protein expression level of HIF-1α (Fig. 5A, Supplementary Fig. 6D). In particular, UBE2S overexpression had no effect on the RNA level of VHL (Supplementary Fig. 6E), indicating that UBE2S negatively regulates VHL in a ubiquitinproteasome-dependent manner.

Further, we performed coimmunoprecipitation and immunofluorescence assays and found an interaction between endogenous UBE2S and VHL in several HCC cell lines (Fig. 5B, C, Supplementary Fig. 7A, B). Consistently, surface plasmon resonance (SPR) demonstrated that UBE2S (1-156 aa, the conserved “core” ubiquitin conjugating domain [21,22]) showed rapid association with the sensor chipimmobilized VHL, producing a concentration-dependent resonance signal, with a calculated dissociation constant (K_D) of 3.56×10^-5 M (Fig. 5D). These results indicate that UBE2S binds to VHL in a direct manner. In addition, we cotransfected VHL plasmid with the HA tag and ubiquitin plasmid with the Myc-Tag into Huh-7 or HCCLM3 cells and found that UBE2S overexpression significantly enhanced VHL ubiquitination (Fig. 5E). Screening of the mutated forms of ubiquitin for potential lysine ubiquitination types revealed that K11O (Ub with an intact Lys11 residue alone), but not K6O, K27O, K29O, K33O, K48O, or K63O, could be linked to VHL by UBE2S (Fig. 5F). Moreover, UBE2S failed to link K11R (the Ub with only the Lys11 residue was mutated) to VHL in HCCLM3 cells (Fig. 5G). Since substrate proteins are labeled with ubiquitin molecules at lysine residues [23], three lysine sites of VHL were identified and were mutated separately (namely K159R, K171R, and K196R) to determine the potential UBE2S-catalyzed lysine sites. The IP assay demonstrated that mutations in K171R and K196R abolished UBE2S-mediated ubiquitination (Fig. 5H), indicating that the K171 and K196 sites were specifically responsible for the UBE2S-catalyzed ubiquitination of VHL. Next, we wanted to test whether UBE2S directly ubiquitinates VHL in an E3 enzyme-independent manner. In vitro ubiquitination assay showed that coincubation of Ub, E1, UBE2S, and VHL in the absence of additional E3 enzymes leads to the appearance of high-molecular-mass bands corresponding to ubiquitinated VHL (Fig. 5I). Taken together, these results indicate that UBE2S enhances the ubiquitination of VHL and promotes its degradation in an E3 enzyme-independent manner, thus improving the protein stability of HIF-1α.

To further explore the effect of UBE2S on the promotion of HCC in vivo, hepatocyte-specific UBE2S knockout (UBE2S^−/−) mice were established by the crossing of UBE2S^loxp/1oxp mice with the Alb-Cre mice (Supplementary Fig. 8A). The mice were then treated with diethylnitrosamine (DEN)/CCl₄ to induce liver tumors (Fig. 6A). The number of liver tumors, the maximum tumor diameter, and the ratio of liver weight to body weight were significantly reduced in UBE2S^−/− mice compared with that of UBE2S^+/+ mice in the same litter (Fig. 6B, C). However, there was no significant difference in the body weight between the two groups of mice (Fig. 6C). Moreover, the serum levels of alanine aminotransferase (ALT) and aspartate aminotransferase (AST) in UBE2S^−/− mice were significantly lower than those in UBE2S^+/+ mice, indicating attenuated liver damage by UBE2S knockout (Fig. 6D). Kaplan–Meier survival analysis also showed significantly prolonged overall survival of UBE2S^−/− mice compared with UBE2S^+/+ mice (Fig. 6E). Further, hematoxylin-eosin (HE) staining was used to confirm the morphology of tumor cells (Fig. 6F). Immunohistochemical assays showed that the expression levels of UBE2S, PCNA, and HIF-1α in the tumor tissues of UBE2S^−/− mice were lower than those of the control group, whereas the expression of VHL was significantly increased (Fig. 6F). The differences in expression of UBE2S, HIF-1α, and VHL in tumor tissue were also confirmed in subsequent Western blot assays (Fig. 6G). In addition, the protein and RNA levels of HIF-1α-regulated metabolic enzymes in tumor tissues of UBE2S^−/− mice were also downregulated compared with those of UBE2S^+/+ mice (Fig. 6G, H). Since the RNA level of HIF-1α and VHL in tumor tissues between the two groups of mice showed no significant difference (Supplementary Fig. 8B), these in vivo results confirmed that the upregulation of glycolytic enzymes induced by UBE2S is dependent on its post-transcriptional regulation of VHL/HIF-1α signaling pathway. Particularly, central carbon-targeted metabolomic analysis was performed to detect differences in metabolites within tumor tissues between the two groups of mice. The results demonstrated that G-6-P, F-6-P, F-1,6-BP, 3-PG, 2-PG, and PEP in tumor tissues of UBE2S^−/− mice were significantly down-regulated compared with those of UBE2S^+/+ mice (Fig. 6I). Taken together, these results suggest that UBE2S acts as an oncogenic gene that promotes HCC progression by regulating glycolysis in vivo.

Our previous studies have shown that cephalomannine inhibits the expression of UBE2S [24]. To further evaluate its pharmacodynamic effects, we performed cytotoxicity assays and found that cephalomannine significantly inhibited cell viability in Huh-7 and HCCLM3 cell lines, with IC50 values of 0.11 and 0.06 μM, respectively (Fig. 7A, Supplementary Fig. 9A). Although the viability of different HCC cells gradually decreased with increasing concentration of cephalomannine, it eventually stabilized, indicating that the use of cephalomannine as a therapeutic agent for HCC needs to be improved. Since UBE2S promotes the proliferation of HCC cells by regulating HIF-1α-dependent glycolysis, the effect of cephalomannine combined with PX-478 in the treatment of HCC was further evaluated. Cytotoxicity assays revealed that PX-478 significantly inhibited cell viability in Huh-7 and HCCLM3 cells, with IC50 values of 46.51 and 21.85 μM, respectively (Fig. 7B, Supplementary Fig. 9B). Moreover, compared with single-compound therapy, the combination of the two compounds significantly increased the inhibitory effect on the viability of different HCC cells in vitro (Fig. 7C, Supplementary Fig. 9C). In addition, we constructed subcutaneous inoculation models in nude mice using Huh-7 cells that stably expressed firefly luciferase and found that luciferase intensity, tumor size, and tumor weight in the treatment group were significantly reduced compared to those in the negative control group (Fig. 7D–I). In particular, the above indices in the two-compound combination treatment group were further decreased compared to those in the single-compound treatment group (Fig. 7D–I). Taken together, these results indicate that the combination of cephalomannine and PX-478 significantly improves anti-tumor effects, thus providing a potential therapeutic strategy for patients with HCC.