Clinicopathological data and mRNA-seq profiles of HCC patients were obtained from three separate databases. A total number of 371 HCC samples with 50 cases of paired non-cancerous liver tissues from TCGA (https://portal.gdc.cancer.gov; ID: LIHC) were included in the exploratory cohort. Clinicopathological parameters including age, gender, grade, AJCC stage, serum alpha fetoprotein (AFP) value, Child–Pugh grade, vascular invasion (VI) type, and mutation information of TP53 are shown in Supplementary Material S1.

Two external validation cohorts GSE76427(n = 115) and LIRI-JP (n = 231) with transcriptomic profiling were extracted from gene expression omnibus (GEO) repository (https://www.ncbi.nlm.nih.gov/geo; Code: GSE76427) and International Cancer Genome Consortium (ICGC) database (https://dcc.icgc.org/projects; Code: LIRI-JP), respectively. Detailed information is exhibited in Supplementary Material S2. As necessary, Log2-transformed data of probe-level expression values were used for further analysis.

Six cell cycle-related gene sets were obtained from the Molecular Signatures Database (MSigDB; Figure 1A). Then, Gene Set Enrichment Analysis (GSEA) was performed to screen out genes that varied at the transcriptome level between cancerous tissue and normal liver epithelium with 1,000 permutations. For each analysis, statistical significance was deemed by false discovery rate (FDR) <0.25, normalized P < 0.05, and |Normalized enrichment score (NES)| >1.65 with 1,000 gene set permutations.

Using TCGA cohort, univariate Cox analysis was conducted on the cell cycle-related genes to identify mRNAs associated with OS. Subsequently, these mRNAs were screened out for risk score model construction using a multivariate Cox analysis. The risk score of each sample was calculated by the following formula (Coef: coefficient β; x: value of gene expression):

Median risk score served as the cutoff point to stratify HCC samples into high- or low-risk groups.

Enrichment of potential pathways in the high-risk group was qualified by GSEA (9). BioCarta (c2.cp.biocarta.v7.2.symbols.gmt), Reactome (c2.cp.reactome.v7.2.symbols.gmt), PID (c2.cp.pid.v7.2.symbols.gmt), and Hallmark (h.all.v7.2.symbols.gmt) gene sets were acquired from the MSigDB database. |NES| >1.65 and normalized P < 0.05 were considered as significant enrichment.

For quantitative evaluation of prognosis, a nomogram combining the CCS with clinicopathologic characteristics of HCC samples was plotted using the “rms” R package. Calibration analysis and the AUC were used to examine the efficiency of the nomogram.

Cancerous samples (n = 23) with paired benign lesions (n = 23) were obtained from patients who underwent hepatectomy at the First Affiliated Hospital of Chongqing Medical University from July 2020 to February 2021. Samples were diagnosed as HCC in the Pathology Department of Chongqing Medical University. This research was approved by the Ethics Committee of Chongqing Medical University (Ethical Approval Number: 2021-556). Relevant clinical parameters of the patients, including age, gender, Edmondson–Steiner grade and AJCC stage, serum AFP value, and microvascular invasion (MVI) information, are provided in Supplementary Material S3.

Human HCC cell lines (Hep3B, HepG2, and Huh-7) procured from the Institute of Chinese Academy of Sciences were routinely maintained in our laboratory. All cell lines were authenticated using STR analysis, and the STR profiling report is shown in Supplementary Material S4. HCC cells were cultured in Dulbecco’s modified Eagle’s medium (HyClone, Logan, UT, USA) with 10% fetal bovine serum (FBS; Gibco, Australia) and 1% penicillin/streptomycin. Small interfering RNA (siRNA) oligonucleotides targeting TICRR were designed and synthetized by GenePharma Biological Technology (China). TICRR-siRNA sequences were as follows: si-TICRR-1, 5′-GGGCCUUCAAGUUCUUUGATT-3′ (sense) and 5′-UCAAAGAACUUGAAGGCCCTT-3′ (anti-sense); si-TICRR-2, 5′-GCCAGCUUCAGGUAUUUCUTT-3′ (sense) and 5′-AGAAAUACCUGAAGCUGGCTT-3′ (anti-sense); si-TICRR-3, 5′-GACCAAAGUUCGAAGAAAUTT-3′ (sense) and 5′-AUUUCUUCGAACUUUGGUCTT′(anti-sense). Cells were transfected with 50 pmol siRNAs by the Lipofectamine 2000 (Invitrogen) for 6 h. Functional tests were performed at 3 days after transfection.

Total RNA was extracted from HCC tissues or cell lines using TRIzol reagent (TAKARA, Japan). Reverse transcription of total RNA samples was performed using RT Master Mix for qPCR kit (MCE; No.: HY-K0511). RT-qPCR was performed using ABI Applied Biosystems Prism 7500 Sequence Detection System. Each reaction was conducted with 10-μl mixture containing SYBR Green qPCR Master Mix (MCE; No.: HY-K0501), and glyceraldehyde-3-phosphate dehydrogenase (GAPDH) was used for normalization. Primer sequences are provided in Supplementary Material S5.

Cell viability was detected by the Cell Counting Kit-8 (CCK-8, Biosharp, China) and the 5-ethynyl-2′-deoxyuridine (EdU) kit (RiboBio, C10310-1, China), according to the manufacturer’s instructions. Here, 2,000 cells and 8,000 cells seeded in 96-well plates were used for CCK-8 and EdU prefiltration assay respectively. Images of EdU assay were acquired by a fluorescence microscope (ZEISS, Germany).

Briefly, HCC cells were fixed in 70% ice-cold ethanol overnight at 4°C and stained with propidium iodide (50 μg/ml). Cell cycle distribution was then detected using a FACSCalibur flow cytometer (BD Biosciences).

All analyses were performed using GraphPad Prism 9 software and R studio 4.0.3. Survival difference was assessed using log-rank test. The variables were presented as the mean ± standard deviation by three separate experiments. For all quantitative analyses, the analyses were blinded, and three observers performed them. P < 0.05 was taken as statistically significant.

To investigate the cell cycle-related genes that are aberrantly expressed in HCC, we performed a GSEA from the exploratory cohort (TCGA). Six gene sets were identified to be significantly enriched in HCC samples based on FDR <0.25, NES > 1.65, and nominal P < 0.02 (Figure 1A). Meanwhile, 907 differentially expressed genes (DEGs) were extracted from these gene sets (CORE ENRICHMENT: YES) for further analysis. Thirty-five genes were discovered to be closely related to OS by univariate Cox analysis (P < 0.001; Figures 1B, C and Supplementary Material S6). Notably, the 35 mRNAs were all upregulated in the HCC samples relative to normal peritumoral tissues, and TICRR (also called Treslin) had the highest hazard ratio (HR; HR = 2.323). Stepwise multivariate Cox analysis was then performed to construct the prognostic signature. Subsequently, 13 cell cycle-related genes (TICRR, SPDL1, SGO1, HJURP, CENPA, GINS1, EZH2, BRSK1, NUF2, PLK1, HMMR, E2F2, and CDCA8) were finally developed as an independent indicator of poor prognosis (Supplementary Material S7).

With the median risk score of the exploratory cohort as the threshold, all samples were dichotomized into two categories (Figure 2A). Kaplan–Meier survival analysis indicated that patients with low risk score possessed significantly longer OS (HR = 2.699, 95% CI 1.965–3.969, P < 0.0001; Figure 2D). The predictive efficacy of the CCS was assessed by ROC, and the area under the curve (AUC) values at 1, 2, 3, 4, and 5 years were 0.785, 0.735, 0.721, 0.654, and 0.682, respectively (Figure 2G).

To further confirm prognostic accuracy, two other cohorts from ICGC and GEO databases were employed for validation. Both validation cohorts were divided into two groups by the same algorithm aforementioned. Significantly, the high risk score also implied worse survival in both ICGC (Figures 2B, E) and GSE76427 (Figures 2C, F) cohorts (high risk vs. low risk: HR = 6.85, 95% CI 2.680–9.008; HR = 6.71, 95% CI 2.562–13.47, respectively), consistent with the exploratory cohort described before. Indeed, the model presented a better performance in the two external validation cohorts, especially in GSE76427 cohort (AUC = 0.835/0.822/0.808/0.821/0.826 at 1/2/3/4/5 years, Figure 2I). In the ICGC cohort, all the AUC values were more than 0.7 except for the “5-year” (AUC = 0.740/0.768/0.779/0.838/0.606 at 1/2/3/4/5 years; Figure 2H). These results above suggested that the predictive power of the CCS had an appropriate sensitivity and specificity.

On the basis of the CCS, the univariate and multivariate Cox analyses were performed after adjusting for clinical factors to assess the independence of the risk score in predicting survival. Due to lack of the sample size in M stage (n = 4) and N stage (n = 4), the two parameters were excluded from our analysis. In univariate analysis, the AJCC stage, T stage, and the risk score were observed to be associated with OS (Figure 3A). The multivariate analysis suggested that the risk score based on the CCS was a prognostic indicator independent of other parameters (P < 0.001; Figure 3C).

To verify the prognostic performance of the risk core with OS, the ROC analysis was further implemented to compare the predictive power of clinicopathological factors with that of the risk score. In the exploratory cohort, the AUC value of the CCS-based risk score reached 0.786, which was higher than that of other clinical factors (Figure 3B). Given that T stage had the highest HR from the univariate analysis and the second highest AUC value (AUC = 0.679) from the ROC analysis, we appraised the prognostic accuracy of this factor. The AUC of the signature was larger than that of T stage (signature-AUC = 0.785/0.735/0.721/0.654/0.682 at 1/2/3/4/5 years vs. T stage-AUC = 0.679/0.628/0.650/0.674/0.672 at 1/2/3/4/5 years; Figure 3D) in TCGA cohort, demonstrating that the prognostic model could provide better clinical guidance than T stage.

Several clinical parameters were stratified to measure the predictive power of the CCS in each subgroup. The 13-CCS had no OS-predictive power for the two subgroups of Child–Pugh grade (P > 0.05; Figure 4G). However, high risk score presented shorter survival in the other subgroups, including age ≤65 years (HR = 2.164), age >65 years (HR = 2.821; Figure 4A), men (HR = 2.703), women (HR = 2.264; Figure 4B), G1–2 (HR = 2.404), G3–4 (HR = 2.203; Figure 4C), stages I–II (HR = 2.624), stages III–IV (HR = 1.921; Figure 4D), T1–2 (HR = 2.422), T3–4 (HR = 2.139; Figure 4E), AFP ≤400 (HR = 1.712), AFP >400 (HR = 3.641; Figure 4F), VI(-) (P = 0.0002), and VI (+) (HR = 3.081; Figure 4H). Stratification analysis of the 13-mRNA model further verified that the CCS predicted the OS with precision in division of these parameters.

The relevance of the CCS-based risk score to each clinical feature was evaluated. Stratification was performed according to the same method, and the results suggested that the risk score had no association with age, gender, grade, Child–Pugh grade, and VI of patients. However, HCC samples in advanced clinicopathological stage (Stages III–IV or T3–4, P < 0.001) and high level of AFP (AFP >400, P = 0.0492) were positively correlated with the risk score (Figures 5A–H), reflecting that the CCS-based risk score was involved in HCC progression.

Besides, analysis of distinct expression of the 13 genes among different AJCC stages in HCC samples showed that the expression of the 13 mRNAs was positively correlated with AJCC stage (Figures 6A–M). Similarly, because the sample size in Stage VI (n = 4) was too small, we did not count it for statistics. Almost all the 13 genes were upregulated in later AJCC stage, such as stage III vs. Stage I (P < 0.05) and Stage II vs. Stage I (except BRSK1, P < 0.05). Interestingly, only TICRR and HJURP expression had a difference between Stage II and Stage III (Figures 6C, M). In addition, through Kaplan–Meier curve analysis, a higher TICRR level resulted in significantly worse prognosis in TCGA datasets (HR = 1.842; Figure 6N) and ICGC datasets (HR = 4.848; Figure 6O), while for the GSE76427, no significant association with survival was found (Figure 6P).

The nomogram integrating clinicopathological parameters and the CCS-based risk score was developed to provide a clinical usability of the prognostic model (Figure 7A). Calibration plots for both 3-year OS and 5-year OS confirmed considerable coherence of ideal prediction with actual observations (Figure 7B). AUC at 3-year OS and 5-year OS reached 0.768 and 0.75, respectively (Figure 7C), showing the robustness and accuracy of the nomogram.

To obtain mechanistic grasp of the cell cycle-related processes in biological functions, we performed a pathway analysis in the high-risk group using BIOCARTA, Hallmark, REACTOME, and PID gene sets by GSEA. As shown in Figures 8A–D, several canonical pathways, such as MCM, MAPK, Glycolysis, Wnt, and Myc pathways, which associated with proliferation/aggressiveness phenotype, were enriched in the high-risk samples. Consistent with this, the high-risk group had higher expression of the proliferation marker gene “KI67” than that in the low-risk group in TCGA cohort (P < 0.001; Figure 8E). Given that patients of the “proliferation class” usually have more frequent TP53 mutations and worse outcomes, we further examined the relationship among these factors. As expected, patients with TP53 mutations have a higher risk score (P < 0.001) and worse prognosis (HR = 1.495) (Figures 8F, G).

The function of the 13 genes in cell cycle regulation and their roles in HCC were analyzed. Ten of the 13 cell cycle genes, including SGO1 (10), HJURP (11), CENPA (12), GINS1 (13), EZH2 (14), NUF2 (15), PLK1 (16), HMMR (17), E2F2 (18, 19), and CDCA8 (20), have been confirmed by previous studies, indicating that their overexpression could increase the malignant phenotype of HCC (Figures 9A, B). Five of these have been identified as activators for G1/S transition in HCC cells, including HJURP, CENPA, EZH2, NUF2, and E2F2. GINS1 and SGO1 regulated S-phase and M-phase duration respectively. HMMR participated in both G1/S and G2/M transitions (17), while elevated E2F2 or PLK1 activated G2/M transitions in HCC cells. However, TICRR, SPDL1, and BRSK1 had no confirmed role in HCC.

All 13 genes were all strikingly overexpressed in the HCCs relative to adjacent non-neoplastic tissues in TCGA cohort (P < 0.001; Figure 10A), and detailed gene expression data are available in Supplementary Material S8. RT-qPCR was performed to further detect the expression of the three function-unknown genes. SPDL1 expression was decreased in HCCs compared with non-neoplastic liver tissues (Figure 10B). The expression of BRSK1 mRNA tended to be higher in HCC relative to paired normal liver tissue, but without reaching statistical significance (P = 0.198; Figure 10C). Importantly, we observed that expression of TICRR mRNA was elevated in HCCs relative to respective controls (P < 0.05; Figure 10D), consistent with TCGA RNA-seq findings. We also noticed that TICRR at mRNA level was negatively correlated with OS in HCC patients (Figures 6N, O), suggesting its potentially important role in tumor progression.

After comparing TICRR expression in three HCC cells (Figure 10E), HepG2 and Hep3B cells were chosen for transfection with three independent siRNAs targeting TICRR. Based on knockdown efficiency, two siRNAs (si-TICRR-1 and si-TICRR-3) were employed in depletion of TICRR (Figure 10F). The effect of TICRR on HCC cell proliferation was studied by CCK-8 (Figure 10G) and EdU assays (Figures 11A, C), and the results demonstrated that depletion of TICRR in Hep3B and HepG2 cells by 2 separate TICRR siRNAs led to decreased cell proliferation. In addition, the cell cycle analysis revealed that silencing TICRR restrained G1/S transition in Hep3B and HepG2 cells by 2 independent si-TICRRs (Figures 11B, D). Taken together, downregulation of TICRR attenuated proliferation by promoting G1/S arrest in HCC cells.