Gene expression profiles of RNA-sequencing data and corresponding clinical information of HCC patients were extracted from TCGA database (https://portal.gdc.cancer.gov/repository). Another dataset consisting of 231 samples was downloaded from the ICGC database. The TCGA-LIHC dataset was used as the training cohort and ICGC dataset as the validation cohort. As the data from both TCGA and ICGC are publicly available. Therefore, this study was exempted from approval by the local ethics committee. Patients with no follow up data from both TCGA and ICGC were excluded from the analysis. All read count values were normalized.

A total of 475 CSC genes from 109 CSC-related gene sets were obtained from a published Articles (Liang et al., 2020). The 109 CSC-related gene sets are provided in Supplementary Table S1.

The “limma” R package was used to identify the differentially expressed genes (DEGs) between tumor and adjacent normal samples in the TCGA cohort, with false discovery rate (FDR) < 0.05 set as the cut-offs. Next, univariate Cox regression analysis of overall survival (OS) was performed to determine the prognostic value of CSC-related genes. 79 prognostic genes was obtained after examining the intersection between the two gene groups. Furthermore, to minimize overfitting, prognostic CSC-related genes were assessed by least absolute shrinkage and selection operator (LASSO) Cox proportional hazards regression using the “glmnet” R package. The value of penalty parameter (λ) corresponding to the lowest partial likelihood deviance was used to select the best model by 10-fold cross-validation. We got a list of genes with non-zero beta coefficients. Finally, a stepwise multivariate Cox regression analysis were used to establish a score system for calculating the survival risk for each HCC patient based on the expression level of each prognostic gene and its related regression coefficient. The following formula was used: risk score = sum (gene expression level * regression coefficient). Moreover, all patients were classified into high-risk or low-risk group according to the median value of risk score, which was used as the cut-off value. To observe the clustering conditions of the gene signature, principal component analysis (PCA), and t-SNE were performed using the “prcomp” function of the “stats” R package and the “Rtsne” R package, respectively. Time-dependent receiver operating characteristic (ROC) curves were generated using the “survivalROC” R package for evaluating the predictive capacity of the novel gene signature. The survival curves for different groups were analyzed with the Kaplan-Meier method with log-rank test. Finally, univariate and multivariate Cox regression analyses were performed to determine whether the CSC-related gene signature possessed the independent prognostic value. To verify the effectiveness of the model, the β value derived from the TCGA set was applied to the ICGC set, all patients come from ICGC were also classified into high-risk or low-risk group according to the same median value of risk score in TCGA set.

Based on other clinical features (grade, age, TNM stage, and T stage) of TCGA-LIHC patients, univariate and multivariate Cox regression analyses were conducted to explore whether the prognostic model was an independent variable. To confirm the prognostic significance of the predictive model, TCGA-LIHC patients were divided into two groups according to different clinical characteristics. Patients were separately classified into the following subgroups: grade I/II, grade III/IV, stage I/II, stage III/IV, age <65, age ≥65, T1-T2, and T3-T4 subgroups. Survival outcome analysis was then performed to verify the independent prognostic significance of the gene signature in specific subgroups. The ideal cut-off value of the risk score was established using the surv_cutpoint function of “survminer” in R package.

Single-sample gene set enrichment analysis (ssGSEA) was performed using “ssGSEA” R package to calculate the infiltrating score of 16 immune cells and the activity of 13 immune-related pathways of patients in the high-risk and low-risk groups. Supplementary Table S2 shows the representative gene sets of immune cells and related pathways.

Gene set variation analysis (GSVA) is a non-parametric and unsupervised GSE method which can calculate enrichment scores of predefined gene sets representing various biological processes in each sample (Hanzelmann et al., 2013). The GSVA was performed using the “clusterProfiler” R package to convert the gene expression profiles of patients in high-risk and low risk groups from TCGA and ICGC cohorts into the enrichment scores in biological signaling pathways or functions. The predefined pathway gene sets were obtained from the Kyoto Encyclopedia for Genes and Genomes (KEGG) database.

To identify the risk score-related hub genes, WGCNA was applied to HCC samples from the TCGA database using the “WGCNA” R package. This analysis method has two parts. One part classifies genes into different modules according to their expression patterns and the other identifies modules that are highly correlated with traits (Langfelder and Horvath, 2008). The WGCNA analysis was carried out on the identified DEGs in HCC tumor samples. First, a gene expression similarity matrix was constructed through calculating the Pearson correlation coefficient between any two genes. Next, a soft threshold of β = 8 was used to judge the similarity of two genes, followed by converting the similarity matrix into a weighted adjacency matrix. Finally, a topological overlap matrix (TOM) was applied to further measure the connectivity of genes in the co-expression network. Genes were clustered according to the value of 1-TOM. In addition, different modules were divided from the identified DEGs by building the dynamic pruning tree, with each module containing at least 30 genes. Similar modules were merged at the cut-off value of 0.25.

Co-expression network was constructed and genes were divided into modules, after which the relationships among modules and traits were determined, which included risk score and risk groups (low-risk = 0, high-risk = 1). Highly correlated modules were selected for further analyses. Gene ontology (GO) and KEGG analyses were then performed using the “clusterProfiler” R package to evaluate the biological functions of genes in the modules.

A nomogram was constructed based on gender, stage, grade, age, and risk score (Iasonos et al., 2008). The area under the ROC curve (AUC), 1-, 3-, and 5-years calibration curves, and decision curve analyses (DCA) (Vickers and Elkin, 2006) were then used to assess the nomogram’s prediction accuracy and discriminatory capacity.

Human HCC cell lines Huh-7 and human kidney cells lines HEK293T were purchased from American Type Culture Collection. All cells were cultured in Dulbecco’s modified Eagle medium (HyClone, USA) supplemented with 10% fetal bovine serum (Gibco, USA) and 100U penicillin/streptomycin (HyClone, USA) at 37°C in a humidified thermal incubator with 5% CO₂.

Two optimal guide RNAs (CCA​GTC​TAT​TGG​CCG​CAG​AT and AAT​AGT​TCA​GAT​GTC​GAT​AT) targeting come from the different sites of the human SGO2 gene exon 6 were designed on the CRISPOR website (http://crispor.tefor.net/). Next, the gRNAs were cloned into pLentiCRISPRv2 (Addgene plasmid #52961). HEK293T cells (2x10⁵) were first seeded in 6-well plates and then transfected with pLentiCRISPRv2-gRNA plasmid or empty pLentiCRISPRv2 (2 µg) together with two lentiviral packaging plasmid psPAX2 (1 μg, Addgene plasmid #12260) and pMD2.G (1 μg, Addgene plasmid #12259) to generate lentivirus. After transfection for 48 h, the upper media containing lentivirus was collected and filtered using 0.45 µm membrane filters. On the other hand, Huh-7 cells were seeded in 6-well plates at a density of 2x10⁵ cells per well. After incubation for 24 h, the media were replaced with a mixture of two lentiviruses and fresh media (1:1:2) containing polybrene (10 μg/ml). The infected cells were selected by treatment with 1 μg/ml puromycine for 4 days until the death of control Huh-7 cells. Finally, the expression of SGO2 was determined using quantitative real-time polymerase chain reaction (qRT-PCR).

Cell Counting Kit-8 (CCK-8, Dojindo Laboratories, Kumamoto, Japan) was used to explore the effect of SGO2 knockdown on cell proliferation of HCC cells. Briefly, cells were first seeded into 96-well plates at a density of 1x10³ cells per well. Next, the CCK-8 solution (10 µl in 90 µl DMEM) was added to each well and incubated at 37°C for 1 h. Finally, the optical density of the medium was measured at a wavelength of 450 nm.

Colony formation experiment: Cells were placed in 6-well plates (1,000 cells per well). After 14 days of incubation, the colonies were fixed in 4% paraformaldehyde (Cat. 15,700, Electron Microscopy Sciences, USA) for 30 min, stained with 0.1% crystal violet (Beyotime, Beijing, China) for 15 min, and then washed it. Finally, colonies were photographed and counted using ImageJ software.

Total RNA of wildtype and SGO2 knockdown Huh-7 cells were extracted using TRIzol reagent (TAKARA, Japan). qRT-PCR was then performed with MonAmp™ ChemoHS qPCR Mix (Monad) in accordance with the manufacturer’s instructions. Relative expression of SGO2 mRNA was normalized to β-actin mRNA (internal control) and calculated using the 2^-ΔΔCt method. The sequences of primer used were as follows: SGO2: 5'- GCC​CAG​TCT​ATT​GGC​CGC​AG -3' (forward) and 5'- TTC​AAT​CTT​TTC​CCC​AAT​AT -3' (reverse); β-actin: 5'- CAC​CAT​TGG​CAA​TGA​GCG​GTT​C -3' (forward) and 5'- AGG​TCT​TTG​CGG​ATG​TCC​ACG​T-3' (reverse).

The students t test was used to evaluate gene expressions differences between normal and tumor tissues, while the chi-square test was used to compare proportional differences between normal and tumor tissues. The Kaplan Meier and log-rank tests were used for comparisons of OS for various clinical subgroup. Univariate and multivariate Cox regression analysis were used to identify independent predictors of OS. For data analysis, the R program (version 4.1.0) was utilized. Unless otherwise stated, p < 0.05 denotes significance.

Supplementary Figure S1 shows the flow chart of the study. The detailed clinical characteristics of these patients come from TCGA are summarized in Table 1. A total of 160 CSC-related DEGs were identified between HCC cancer samples and adjacent normal samples (Supplementary Figures S2A,B). Among them, 79 DEGs were found to be associated with OS in the univariate Cox regression analysis (Figures 1A,B). Analysis of the heatmap plot showed that most of these prognostic DEGs were upregulated in the cancer samples (Figure 1C). GO analysis revealed that these CSC-related genes were mainly enriched in response to stem cell differentiation, mesenchyme development, mesenchymal cell differentiation, maintenance of stem cell population, and stem cell development (Supplementary Figure S2C). Moreover, KEGG pathway analysis revealed that these genes were correlated with the hippo signaling pathway, Wnt signaling pathway, TGF-beta pathway, and notch signaling pathway (Supplementary Figure S2D).

To further filter the genes and reduce the risk of overfitting, LASSO Cox regression was employed. According to the multivariate Cox regression of OS, A 3-gene signature (RAB10, TCOF1, and PSMD14) was eventually identified based on the minimum value of λ (Figures 1D,E). Next, the risk score of each HCC patient in the TCGA dataset was calculated using the following formula: 0.238*expression level of RAB10 + 0.609*expression level of TCOF1+0.360*expression level of PSMD14. Patients were stratified into two groups based on the risk score; patients with risk scores higher than the median value were classified in high-risk group (n = 182), whereas patients with risk scores lower than the median value were classified in low-risk group (n = 183) (Figure 2A). Results showed that the three genes were significantly upregulated in the high-risk group (Figure 2B). PCA and t-SNE analysis revealed that patients in the two groups were well distributed in different directions (Figures 2D,E). Further analysis showed that patients in the high risk score group had shorter the survival time and mortality risk (Figure 2C). Time-dependent ROC curves were generated to evaluate the prediction power of the risk score. Results demonstrated that the area under the curve (AUC) for 1-, 2-, and 3-years OS was 0.775, 0.686, and 0.675, respectively, indicating that the gene signature could predict the prognosis of HCC (Figure 2F). Analysis of Kaplan-Meier survival curves showed that patients in the high-risk group had poor OS and platinum-free interval (PFI) (p < 0.001; p = 0.025, Figures 2G,H).

The dataset downloaded from the ICGC was used to verify the generality of the novel CSC-related prognostic model. The risk score was used to divide patients from the ICGC dataset into high-risk groups and low-risk groups (Figure 3A). Similarly, the three prognostic risk genes were upregulated in the high-risk group (Figure 3B). PCA and t-SNE analyses indicated that patients in the different risk groups were distributed in two discrete sections (Figures 3C,D). Results of the scatterplot and Kaplan-Meier survival curves showed that patients with high-risk scores had short survival time and high mortality rate (Figures 3E,F), consistent with findings in the TCGA cohort. The AUC of time-dependent ROC ranged from 0.705 to 0.703 for 3 years (Figure 3G), indicating that the gene signature could predict the prognosis of patients in the ICGC cohort. Moreover, the risk score showed better prediction than most the clinical characteristics such as gender and age (Figure 3H).

To further evaluate the prognostic value of the risk score, univariate and multivariate Cox regression analyses were applied among risk score and some clinicopathological characteristics in both TCGA and ICGC cohorts. The univariate Cox regression analysis results showed that the risk score (p < 0.001, HR = 1.623, 95% CI = 1.422-1.854) and tumor stage (p < 0.001, HR = 2.500, 95% CI = 1.721-3.632) were significantly associated with OS in the TCGA cohort (Figure 4A). In the multivariate Cox regression analysis of the two variables, the risk score (p < 0.001, HR = 1.635, 95% CI = 1.418-1.884) and tumor stage (p < 0.001, HR = 2.361, 95% CI = 1.624-3.432) were independent prognostic factor in the TCGA dataset (Figure 4C). The risk score was also an independent prognostic factor in the ICGC cohort (p < 0.001, HR = 3.534, 95% CI = 1.862-6.709) (Figure 4B, Figure 4D). Collectively, these results suggest that the constructed risk score was an independent prognostic predictor of OS.

To assess the prognostic value of the model in various clinic-pathological subgroups, clinical variables and samples were randomized into two subgroups in term of TNM stage, age, grade, and T stage. The obtained results indicated that the 3-CSC-related genes signature was significantly correlated with the survival of patients in different clinic-pathological subgroups in the TCGA-LIHC cohort (Supplementary Figure S3).

HCC patients were classified into high-risk and low-risk groups based on the risk score. It was evident that patients with high-risk scores were more likely to have poor prognosis, which has been confirmed above. However, the correlation between poor prognosis and high-risk score was not clear. We suspected that there might be a difference in immune status between the two risk groups. To this effect, the ssGSEA method was applied to all HCC samples in the high-risk and low-risk groups to evaluate immune cells infiltration and the activities of immune-related pathways. Figure 5A and Figure 5B show the enrichment scores of 16 types of infiltrating immune cells. In TCGA and ICGC cohorts, the scores of aDCs, macrophages, Th2 cells, and Tregs were significantly upregulated in the high-risk score group, whereas the scores of neutrophils and NK cells were lower than those in the low-risk group (p < 0.05). Among the screened immune-related functions or pathways, the expression of MHC class I was increased in the high-risk score group. In contrast, the activities of type I IFN response and type II IFN response were inhibited in the high-risk score group (p < 0.05, Figures 5C,D).

To further explore the differences in the activities of various biological signaling pathways between high-risk and low-risk groups, GSVA enrichment analysis was performed for each sample from the TCGA and ICGC cohorts. Results showed that the common pathways and functions enriched in samples from high-risk groups in both TCGA and ICGC cohorts were homologous recombination, cell cycle, RNA degradation, splicesome, and ubiquitin mediated proteolysis (Figures 6A,B). On the other hand, the activities of various metabolism-related pathways in the low-risk groups were substantially higher than those in the high-risk groups, including glycine-serine and threonine metabolism, fatty acid metabolism, linoleic acid metabolism, and the PPAR signaling pathway. In addition, one immune-related pathway, complement and coagulation cascades, was enriched in the low-risk groups (Figures 6A,B).

To screen out the hub genes responsible for the high-risk of poor prognosis in HCC, WGCNA was performed on the DEGs between tumor and normal tissues identified from TCGA database. A co-expression network of 5,301 genes was built and 14 modules were divided from it through setting the optimal soft-thresholding power at 8 (scale free R² = 0.85, Figures 7A–C). The correlation coefficients between each module and risk score or risk groups (high-risk = 1; low-risk = 0) were calculated and presented in a heatmap. The results suggested that the cyan module was the most closely and positively associated with risk score and the high-risk group, with correlation coefficients of 0.8 and 0.67, respectively (p < 0.001, Figure 7D). Therefore, the cyan module was selected for further analysis. GO and KEGG enrichment analyses were conducted on the 2,815 genes in the cyan module to reveal the functional link between key genes. The top 10 GO terms and KEGG pathways enriched from the genes of cyan module are shown in Supplementary Figure S4.

Genes in the module were ranked by the scores and p values of gene significance (GS) and module membership (MM). Notably, the top 10 genes are listed in Table 2. Given that the mechanism of the number 1 gene (BUB1B) in HCC has already been elucidated (Qiu et al., 2020), we chose the number 2 key gene (SGO2) for subsequent analyses. Notably, in vitro experiments were conducted to further investigate the biological functions of SGO2 in HCC. Specifically, the expression level of SGO2 was knocked down in Huh-7 cells by CRISPR/Cas 9 system as illustrated in Figure 8A. The success of SGO2 knockdown was confirmed using qRT-PCR analysis (Figure 8B). The CCK-8 assay results showed that decreased expression of SGO2 had a significant inhibitory effect on the proliferation of Huh-7 cells (Figure 8C).Clone formation experiments showed, Huh-7 cells with knockdown of SGO2 showed worse ability to form colonies.The difference was statistically significant (p < 0.05), as shown in Supplementary Figure S5. The results showed that SGO2 knockdown effectively inhibited the clone formation ability of HCC cells.

To make the prognostic signature more convenient for clinical application, we constructed a nomogram and tested its capacity to predict OS on the TCGA-LIHC cohort at 1-, 3-, and 5-years based on risk scores, and other signifcant independent risk factors (Figure 9A). The 1-, 3-, and 5-years OS calibration curves of the TCGA-LIHC data revealed that the nomogram had good predictive discrimination and accuracy (Figure 9B). Moreover,the nomogram had a higher consistency index compared to other clinical markers and the risk score (Figure 9C). In addition, Comparison of the net benefits of various models, such as none, risk score, all, nomogram, and clinical indicators revealed that the nomogram had a higher net income and a wider threshold probability (Figure 9D). The analysis results revealed that the nomogram had a better prognostic ROC value compared to the 3-CSC-related-genes signature and other four clinical indicators, and it could predict OS for 1-, 3-, and 5-years (Figures 9E–G).Finally, the predictive value of nomogram was further validated in the ICGC dataset (Supplementary Figure S6). Altogether, results of ROC, DCA, calibration curve, and C index analyses indicated that the nomogram had better clinical benefits than risk scores based on the four 3-genes signatures alone.