The RNA-sequencing (RNA-Seq) transcriptome profile (TCGA- LIHC) of transcripts per kilobase of exon model per million mapped reads format, including 365 HCC tissues and 50 normal hepatic tissues, were acquired from TCGA database (TCGA-LIHC, http://portal.gdc.cancer.gov) using the TCGAbiolinks package ((https://bioconductor.org/p-ackages/release/bioc/html/TCGAbiolinks.html, version 2.28.4) (20) in R (version 4.3). The accompanying clinical information of 365 patients with HCC was also obtained by the TCGAbiolinks package in R (version 4.3). For duplicate samples, the mean was calculated and used as the final sample gene expression value. Samples clearly identified as HCC were included and samples with missing or incomplete clinical information [age, sex, tumor, node, metastasis (TNM) staging, survival time and status] were removed in subsequent prognostic and modeling analyses. The patients were then divided into high- and low-FOXO3 expression groups according to X-tile software (version 3.6.1, /x-tile.software.informer.com/), with a cut-off value of 4.206.

The RNA-seq transcriptome data and corresponding clinical information of 206 patients with HCC (ICGC-LIRI-JP), including 206 HCC tissues and 177 normal hepatic tissues, were downloaded from the ICGC database (ICGC-LIRI-JP, http://dcc.icgc.org/). The data from ICGC was processed in the same way as the aforementioned database.

An mRNA co-expression network of the TCGA-LIHC dataset was constructed using the WGCNA R package (version1 73) (21). Briefly, the Pearson correlation between each pair of genes was calculated and the similarity matrix was acquired. The power function was used to transform the similarity matrix to an adjacency matrix using the WGCNA R package and the β was determined through the scale-free topological fit test for constructing scale-free weighted network. A threshold of R²>0.9 (soft threshold=32) was selected to acquire a high-confidence scale free network. Co-expression modules were acquired through the pairwise topological overlap between genes and highly correlated modules were further merged. The hub module was designated as that with the highest Pearson correlation coefficient and P<0.05. Furthermore, module membership (MM) and gene significance (GS) were calculated and genes with MM >0.8 and GS >0.2 in the hub module were selected as the hub genes.

LASSO regression analysis was performed using the glmnet (version 4.1-4) (22) and survival (version 3.8-3) (23) R packages. LASSO was then used to identify the final gene signature for constructing the prognostic model according to the minimum λ value. The PPI network was constructed using the STRING database (version 12.0; string-db.org/) and visualized by Cytoscape (version 3.9.1, http://cytoscape.org/).

Univariate Cox analysis was used to identify the hub genes and genes with P<0.05 were retained for subsequent analysis. LASSO analysis was performed and five genes were identified for constructing the model according to the minimum λ value (24). The risk score (RS) of each HCC sample was calculated using the following formula: RS= Σⁿ_i=₁coef (geneⁱ) × expr (geneⁱ). The patients were then divided into high- and low-risk groups according to X-tile software.

Time-dependent receiver operating characteristic (ROC) curves were generated using the pROC R package (version 1.18.5) (25). Univariate and multivariate Cox analyses were performed and visualized using the survival R package, SPSS (version 23; IBM Corp.) and the ggplot2 R package (version 3.5.1) (26). Nomogram and calibration curves were generated using the rms R package (version 7.0) (27).

GSEA of the training (TCGA cohort) and validation sets (ICGC cohort) was performed using GSEA software (version 4.3.2, gsea-msigdb.org/gsea/index.jsp). Hallmark gene sets from the molecular signature database (MSigDB; http://www.gseamsigd-b.org/gsea/msigdb/index.jsp) were defined as the reference gene set database. Pathways with P<0.05 and false discovery rate (FDR) <0.25 were considered statistically different between the high- and low-risk groups. Infiltration of 22 immune cell types in the tumor microenvironment was analyzed using the CIBERSORTx website (cibersortx.stanford.edu/). Results were visualized using the ggplot2 package.

A total of 10 paired tumor and tumor-adjacent tissues from patients with HCC were collected by hepatectomy between May 2024 and August 2024 from the Second Affiliated Hospital of Guangdong Medical University (Zhanjiang, China). The use of samples was approved by Research Ethics Committee of the Second Affiliated Hospital of Guangdong Medical University (approval no. PJKT-2024-042).

A total of 10 HCC samples obtained were included in accordance with the following inclusion criteria: i) The patient received hepatectomy; ii) all samples were confirmed to have a clinicopathological diagnosis of HCC through pathology reports; and iii) patients had no severe infection and stable vital signs. The exclusion criteria were: i) Patients with secondary and recurrent liver cancer; ii) patients with liver cancer who had received any medication prior to surgery; and iii) patients with multiple primary tumors. The clinical information of 10 patients with HCC is presented the Table SI. Total RNA from HCC and paired paracancerous tissues was extracted using the Total RNA Isolation Kit (cat. no. RE-03011; Foregene Co., Ltd.) according to the manufacturer's instructions. RNA was reverse transcribed to cDNA using a cDNA synthesis kit (cat. no. 1708891; Bio-Rad Laboratories, Inc.) and the following conditions: 5 min at 25°C, 20 min at 46°C, 1 min at 95°C and holding at 4°C. cDNA was quantified using Universal SYBR Green Supermix (cat. no. 1708891; Bio-Rad Laboratories, Inc.) according to the manufacturer's instructions (denaturation for 5 sec at 95°C, extension: 30 sec at 60°C, 40 cycles). Relative mRNA expression was obtained using the 2^−∆∆Cq method (28) with GAPDH as the internal reference. The following primers were used: GAPDH forward, 5′-GGAGCGAGATCCCTCCAAAAT-3′; GAPDH reverse, 5′-GGCTGTTGTCATACTTCTCATGG-3′; DDX55 forward, 5′-AGCTGGGCTTCCCGTACAT-3′; DDX55 reverse, 5′-CAGCGACATCTTTGTTTCGCA-3′; RAB10 forward, 5′-CTGCTCCTGATCGGGGATTC-3′, RAB10 reverse, 5′-TGATGGTGTGAAATCGCTCCT-3′; RAB7A forward, 5′-GTGTTGCTGAAGGTTATCATCCT-3′; RAB7A reverse, 5′-GCTCCTATTGTGGCTTTGTACTG-3′; TAF1B forward, 5′-AAAGAACGCTGTACTCAGTGTG-3′; TAF1B reverse, 5′-CCCCGGTTGAGGGCTTTTA-3′; TAF3 forward, 5′-ATGTGCGAGAGTTACTCCAGG-3′; and TAF3 reverse, 5′-GGGTCTGTTCGGCCATAGAG-3′.

Protein from HCC and paired paracancerous tissues were extracted using radioimmunoprecipitation assay buffer with protease and phosphatase inhibitors (cat. no. P0013C; Beyotime Institute of Biotechnology) according to the manufacturer's instructions. Protein samples were quantification by the BCA kit (cat. no. P0012S, Beyotime Institute of Biotechnology). Protein samples were heated at 70°C for 10 min after mixing with sodium dodecyl sulfate-polyacrylamide gel electrophoresis (SDS-PAGE) buffer (cat. no. P0015A; Beyotime Institute of Biotechnology) and then added into a pre-prepared 10% gel for electrophoresis (15 µg per lane). The protein was transferred to polyvinylidene fluoride membranes. Subsequently, the membranes were blocked with 5% skim milk at room temperature for 2 h and then incubated with primary antibody at 4°C overnight. After incubation, the membranes were placed on a shaker and washed with TBST buffer (1 l TBS buffer with 1 ml Tween-20) for 7 min. After washing three times, the membranes were incubated with the secondary antibody at room temperature for 1 h. Finally, the membranes were washed and developed by chemiluminescence. The gray values were obtained by ImageJ (National Institutes of Health, version 1.54h) and analyzed using GraphPad Prism (Dotmatics, version 10.2.3). The following antibodies were used: β-actin (1:1,000; cat. no. ET1702-52; HUABIO), DDX55 (1:1,000; cat. no. ER63225; HUABIO), RAB10 (1:1,000; cat. no. 11808-1-AP; Proteintech Group, Inc.), RAB7A (1:1,500; cat. no. 55469-1-AP; Proteintech Group, Inc.), TAF1B (1:500; cat. no. 12818-1-AP; Proteintech Group, Inc.), FOXO3 (1:1,000; cat. no. A9270; ABclonal Biotech Co., Ltd.) and TAF3 (1:500; cat. no. 18901-1-AP; Proteintech Group, Inc.), HRP conjugated goat anti-rabbit IgG polyclonal antibody (1:30000, HA1001, HUABIO), HRP conjugated goat anti-mouse IgG polyclonal antibody (1:30,000, HA1006, HUABIO).

The Huh7 cell was obtained from Cell Bank of the Chinese Academy of Sciences (Shanghai, China) and cells (10,000/well) were seeded into 6-well plates. Cells were cultured using DMEM (Gibco, USA) supplemented with 10% FBS) (ScienCell, USA) and 1% penicillin and streptomycin (HyClone). After the cell density reached 20%, the cells were transfected with small interfering RNAs (siRNA; Shanghai GenePharma Co., Ltd.). Briefly, the siRNA, serum-free medium and siRNA-Mate transfer agent (cat. no. G04003; Shanghai GenePharma Co., Ltd.) were mixed, added to the cells at a final concentration of 100 nM and shaken well according to the manufacturer's instructions. After 72 h transfection, the protein was extracted to detect the knockdown efficiency via western blot analysis. The following siRNA sequences were used: RAB10 sense, 5′-CCAUAGGAAUAGACUUCAAGA-3′; antisense, 5′-UUGAAGUCUAUUCCUAUGGUG-3′; RAB7A sense, 5′-GGAAGAAAGUGUUGCUGAAGG-3′; antisense, 5′-UUCAGCAACACUUUCUUCCUA-3′; TAF3 sense, 5′-GCGGGAUGUGCGAGAGUUACU-3′; antisense, 5′-UAACUCUCGCACAUCCCGCUG-3′; FOXO3 sense, 5′-AACUAAACCCUUUAGUGACAU-3′; antisense, 5′-GUCACUAAAGGGUUUAGUUUU-3′; negative control sense, 5′-CAUAAAUCUACAGGAUGAUTT-3′; antisense, 5′-AUCAUCCUGUAGAUUUAUGTT-3′.

The cells (5,000/well) were seeded in 96-well plates with three replicates per group. According to the manufacturer's instructions, the absorbance at 450 nm was detected using a CCK-8 kit after incubation for 2 h (cat. no. ZP328-3; Beijing Zoman Biotechnology Co., Ltd.) at 6, 24, 48 and 72 h after seeding the cells. The absorbance value detected 6 h after seeding the cells was considered as 0 h.

Data with normal distribution and skewed distribution was analyzed by t test and Wilcoxon rank-sum test, respectively. The results of the bioinformatics analysis were analyzed using R software (version 4.2.0) and the Wilcoxon rank-sum test. Kaplan-Meier survival curves were compared using the log-rank test or two-stage method. The survival curves were plotted and cut-off was selected by the survival package in R. The results of RT-qPCR, western blot analysis and CCK-8 assay were analyzed using Student's t-test. All statistical analyses were performed using R, SPSS and GraphPad Prism (version 9.3.1; Dotmatics). P<0.05 was considered to indicate a statistically significant difference.

To evaluate the mRNA expression of FOXO3, the transcriptome and corresponding clinical profiles of 365 patients with HCC and 206 patients with HCC were obtained from TCGA and ICGC-LIRI-JP, respectively. The mRNA expression of FOXO3 in HCC tissues was higher compared with that in normal tissues in the TCGA dataset (Fig. 1A). Moreover, the patients with pathology grade G3/4 had higher expression of FOXO3 compared with grade G1/2 (Fig. 1B), but the mRNA expression of FOXO3 was not associated with TNM staging and vascular invasion (Fig. S1A and B). The ICGC database was used to further validate these findings. FOXO3 expression was significantly higher in HCC tissues compared with normal tissues (Fig. 1C). Moreover, FOXO3 was significantly higher in tissues with expressed in bile duct invasion HCC compared with HCC (Fig. 1D), but the mRNA expression of FOXO3 did not differ with higher TNM staging or with vascular invasion (Fig. S1C and D). Both databases indicate that high expression of FOXO3 was strongly associated with poor prognosis in HCC (Fig. 1E and F).

FOXO3, as a transcription factor, may serve an important role in HCC; therefore, WGCNA was used to identify a FOXO3-associated gene module. A scale-free network was constructed using WGCNA in 365 patients with HCC from TCGA (Fig. 2A). The average-linkage hierarchical clustering method was used to cluster genes and modules with >80% similarity were merged (Fig. 2B), resulting in six modules. The correlations between each module and the mRNA expression of FOXO3 in HCC were analyzed (Fig. 2C), which indicated that the red module was most closely related to FOXO3 (Pearson coefficient of 0.72). The gene distribution in the red module was further analyzed, which demonstrated that GS and MM were significantly correlated, suggesting that genes in the red module were strongly associated with FOXO3 (Fig. 2D).

Based on the criteria of MM >0.8 and GS >0.2, hub genes in the red module were identified, resulting in 875 hub genes. Univariate Cox regression analysis was then used to screen 294 prognosis-associated genes from 875 genes with a P-value of <0.05 (Fig. 3A). To further narrow the range of variables, LASSO regression analysis was performed, which identified five genes (DDX55, TAF1B, TAF3, RAB10 and RAB7A) for constructing the model according to the minimum λ value (Fig. 3B and C). Fig. 3D shows the PPI network between FOXO3 and the 294 proteins obtained from univariate Cox regression analysis.

According to the results of LASSO regression analysis, the expression values and regression coefficients of these five genes for each sample were used to calculate the RS of each sample. Subsequently, the cut-off values of the RS were obtained using X-tile software and the patients were divided into high-risk and low-risk groups (Fig. 4A). All five genes were more highly expressed in the high-risk group and the patients with HCC in the high-risk group had a worse prognosis compared with the low-risk group (Fig. 4B). The time-dependent ROC curve showed that the area under the curve (AUC) values for 1, 3 and 5 years were 0.73, 0.69 and 0.71, respectively (Fig. 4C). ICGC data was used to evaluate the robustness of the model and similar results were obtained (Fig. 4D-F). Moreover, the present model was compared with other prognostic models which found the AUC values (1-year) of the present model were higher compared with the AUC values (1-year) of other models in training and validation sets (Table SII) (29–33).

To further evaluate the clinical value of the model, univariate and multivariate Cox regression analyses were used to analyze the RS and other clinical factors (age, sex, TNM staging and vascular invasion) in both datasets. Both univariate Cox analysis (Fig. 5A and D) and multivariate Cox analysis (Fig. 5B and E) showed that RS and TNM staging were risk factors in both datasets. Considering that TNM staging has been accepted as a standard prognostic method, these results suggested that RS has potential for clinical prognostic evaluation. In addition, RS and TNM staging in TCGA and ICGC Cox analysis showed similar statistical significance, although TNM staging has a higher hazard ratio. However, age, sex and vascular invasion were statistically significant only in partial Cox analyses. TNM staging and RS were used to construct nomograms in both datasets to quantify prognostic scores for patients with HCC. For TCGA dataset, the survival probability of patients with HCC with total points >100 was <0.6 (Fig. 5C). Similarly, the 4-year survival probability of patients with HCC with total points >100 was <0.6 in the ICGC dataset (Fig. 5F). The calibration curves of TCGA and ICGC datasets are shown in Fig. S2.

Activation of carcinogenic pathways is an integral part of tumor progression (32). In HCC, tumor development is often accompanied by activation of cell cycle-associated pathways and phosphorylation of the AKT pathway and the activity of these pathways suggests poor prognosis of patients (34). In the present study, GSEA was performed for patients in the high- and low-RS groups, which demonstrated that the cell cycle-associated pathways (G₂/M checkpoint and E2F targets) and AKT pathway were more active in the high-RS group (Fig. 6A-C), partly explaining why patients in the high-risk group had worse outcomes. Similar results were obtained from the ICGC dataset (Fig. S3A-C).

The tumor microenvironment, in which several immune cells play a major role, is another important factor mediating tumor prognosis (35). Thus, the infiltration of 22 immune cell types was quantified using transcriptome data from patients in the high- and low-risk groups. In TCGA and ICGC datasets, M2 macrophages and CD4⁺ T cells were the main immune cells enriched in the tumor microenvironment (Figs. 6D and S3D). In addition, there was more infiltration of M2 macrophages and resting CD4⁺ memory T cells in the high-risk group compared with the low-risk group in TCGA dataset (Fig. 6E). Compared with the high-risk groups, a higher infiltration of activated natural killer cells was found in the low-risk groups in TCGA and ICGC datasets (Figs. 6E and S3E).

The expression of the five genes (DDX55, RAB10, RAB7A, TAF1B and TAF3) used for constructing the model in HCC was further investigated. TCGA and ICGA datasets both showed that these five genes were highly expressed in HCC tissues compared with normal tissues (Figs. 7A and S4A). Moreover, the high expression of these five genes in TCGA dataset was significantly associated with poor prognosis of HCC, while only the expression of the RAB10, RAB7A and TAF3 genes was associated with poor prognosis of patients with HCC in the ICGC dataset (Figs. 7B-F and S4B-F). To validate these results, tumor and matched non-tumor tissues were collected from ten patients with HCC. The RT-qPCR analysis indicated that the mRNA expression of the five genes in tumor tissues was higher compared with that in non-tumor tissues (Fig. 7G). Immunohistochemical results from the Human Protein Atlas (HPA) database suggested that DDX55, RAB10 and RAB7A have higher expression in tumor tissue compared with non-tumor tissue (Fig. 7H) (36). However, it was not possible to obtain representative immunohistochemical staining results for TAF1B and TAF3 from the HPA database. Western blot analysis demonstrated that there was higher protein expression levels of all five proteins in tumor tissues compared with non-tumor tissues (Fig. 7I).

To further explore the influence of the genes of interest on HCC progression, three genes (RAB10, RAB7A and TAF3) that were significantly associated with prognosis in both TCGA and ICGC datasets were selected for functional experiments. The CCK-8 assay demonstrated that knockdown of RAB10, RAB7A and TAF3 inhibited the proliferation of Huh7 cells (Figs. 7J and S4G). Moreover, the AUC values (1-year) among five genes, RS and TNM in TCGA and ICGC datasets were compared which demonstrated the five genes and RS were not inferior to TNM in predicting first-year survival status of patients with HCC (Fig. S5A and B). In addition, the protein expression of FOXO3 was knocked down using siRNA (Fig. S5C). The protein expression levels of the other five proteins (DDX55, RAB10, RAB7A, TAF3 and TAF1B) were also downregulated after knocking down FOXO3 (Fig. S5D). Additionally, mRNA levels of FOXO3 were positively correlated with mRNA levels of the 5 genes in both TCGA and ICGA datasets (Fig. S6).