HB-related data were obtained and downloaded from the Gene Expression Omnibus (GEO¹) database portal using the keyword “hepatoblastoma.” The inclusion criteria for expression profile data were as follows: (a) the organism was Homo sapiens, (b) samples used for gene expression analysis included both HB tissue and normal liver tissue, (c) data for all samples were complete, and (d) HB and normal liver tissue samples could be clearly separated by principal component analysis (PCA). Only datasets that met all of the above criteria were included. Two datasets, GSE75271 (Sumazin et al., 2017) and GSE131329 (Kanawa et al., 2019), were therefore included for further analysis. GSE75271, consisting of 50 HB samples and five normal liver samples, was analyzed using GPL570 platform (Affymetrix Human Genome U133 Plus 2.0 Array), while GSE131329, consisting of 53 HB samples and 14 normal liver samples, was analyzed via GPL6244 platform (Affymetrix Human Gene 1.0 ST Array).

Raw data files (^∗.CEL) from GSE75271 and GSE131329 were downloaded and processed. Data from GPL570 and GPL6244 platforms were imported using R packages affy (Gautier et al., 2004) and oligo (Carvalho and Irizarry, 2010), respectively. The gene expression profile probe names were transformed to gene symbols and Entrez IDs using the hgu133plus2.db R Bioconductor package and the hugene10sttranscriptcluster.db R Bioconductor package, respectively. If one gene symbol corresponded to different probes, combined average levels were considered for gene expression values (Barrett et al., 2013). All raw data were processed using data filtering, a base 2 log transformation, and quantile normalization. The imput.knn function in the impute R Bioconductor package was used for data filtering. After data merging, the ComBat function in the sva R package (Leek et al., 2012) was used for batch effect correction, and the results were verified by PCA. DEGs between HB and normal liver tissue samples were detected using the limma R package (Ritchie et al., 2015), with the following cut-off criteria for significance: adjusted P < 0.05 and |log2FC| > 1.

In order to explore the functional annotation of candidate genes, GO terms and KEGG pathway analyses were performed using the clusterProfiler R package (Yu et al., 2012). GO terms included biological process (BP), cellular component (CC), and molecular function (MF). Adjusted P values below 0.05 were deemed significantly enriched.

We constructed a PPI network using the Search Tool for the Retrieval of Interacting Genes (STRING) online database (Szklarczyk et al., 2015), with an interaction score >0.4 set as the cut-off value. Subsequently, the Cytoscape software was utilized to visualize the PPI network (Shannon et al., 2003). The Molecular Complex Detection (MCODE) (Bandettini et al., 2012) plugin in Cytoscape was applied in order to extract densely connected modules from PPI network, with degree cut-off = 2, node score cut-off = 0.2, K-score = 2, and max depth = 100. The Cytoscape plugin CytoHubba (Chin et al., 2014) was utilized for the identification of key genes from the PPI analysis. We extracted the top 20 genes from both approaches, and the intersecting genes of all four approaches of CytoHubba ranking were deemed as hub genes. The four approaches of CytoHubba ranking used here were maximal clique centrality (MCC), edge percolated component (EPC), maximum neighborhood component (MNC), and node connect degree.

We constructed an unsigned weighted gene co-expression network using the WGCNA R package (Langfelder and Horvath, 2008). After data merging, batch effect correction, and exclusion of outlier samples, the complete gene expression matrix contained 8,204 genes across 116 samples. An expression matrix of 2,051 genes with the top 25% highest variance was used for WGCNA. We conducted hierarchical clustering of samples to remove outliers with a cut-off value of 80 to produce two stable clusters. Then, the soft threshold power β was determined in order to ensure a scale-free network. The resulting Pearson correlation matrix was converted to adjacency matrix via the power function, followed by transformation into a topological overlap matrix (TOM). The TOM was used to calculate corresponding dissimilarity. We carried out hierarchical clustering in order to cluster similar genes into the same module. The dynamic cutting algorithm was then used to detect the gene modules. Subsequently, we clustered the eigengenes according to the relationship and merged them into modules with the association >0.75. Module–trait association between each module and the phenotype was evaluated based on Pearson correlation. For each gene, module membership (MM) was characterized according to the association between module eigengene (ME) and its expression level. The association between gene expression and clinical phenotype represented gene significance (GS). After identifying a module of interest, GS and MM for each gene were computed in the given module. Finally, we performed GO and KEGG pathway analyses to illustrate potential biological functions of the identified module.

Key genes were closely correlated genes in one module with a MM > 0.8 and a GS > 0.2. For subsequent analysis, intersecting genes identified from both PPI and the most significant modules were assessed. Based on GSE75271 and GSE131329 datasets, the expression values of key genes between HB and normal liver tissue samples were then compared.

FBS (cat. no. 10099141), DMEM (cat. no. 11995065), PBS (cat. no. 10010023), and 0.25% Trypsin-EDTA (cat. no. 25200072) were purchased from Gibco (Grand Island, NY). Antibodies against CDK1 (cat. no. ab133327), CCNA2 (cat. no. ab181591), and β-actin (cat. no. ab8226) were obtained from Abcam (Cambridge, MA, United States). Antibodies against CDC20 (cat. no. 4823) and GAPDH (cat. no. 5174S) were procured from CST (Beverly, MA, United States).

Human HB cell lines (HepG2 and HuH-6) were purchased from Shanghai Institutes for Biological Sciences, Chinese Academy of Sciences. Cells were cultured in DMEM supplemented with 10% FBS at 37°C/5% CO₂.

Total proteins were extracted from HepG2 or HuH-6 cells using the RIPA buffer supplemented with a protease inhibitor cocktail. Lysates were separated using sodium dodecyl sulfate polyacrylamide gel electrophoresis and transferred onto a polyvinylidene fluoride membrane. The membrane was then blocked for 1 h using western blocking buffer and subsequently incubated using a primary antibody, followed by incubation for 2 h with IgG HRP-conjugated secondary antibody (Jackson ImmunoResearch, PA, United States). Proteins were detected using ChemiDoc-It system (Tanon, Shanghai, China). Band intensities were assessed using ImageJ. GAPDH or β-actin served as the loading control.

Small interfering RNAs (siRNAs) targeting CDK1, CCNA2, or CDC20, as well as non-targeting control siRNAs, were obtained from RiboBio (Guangzhou, China). siRNAs were transfected into HepG2 or HuH-6 cell lines according to the manufacturer’s guidelines using Lipofectamine 2000. Transfection efficiency was confirmed via western blot (WB) 2 days after siRNA transfection.

HB cells were cultured in six-well plates containing media supplemented with 10% FBS to a density of 3 × 10³ cells/well. The culture media were replaced by media containing 5% FBS the following day, and cells were cultured for 2 weeks. This step was followed by paraformaldehyde (PFA) fixing and staining with crystal violet. Subsequently, photos were taken. Cells were subsequently fixed using PFA, stained with crystal violet, and microscopic images were acquired.

Cell viability was assessed using Cell Counting Kit-8 (CCK-8, Dojindo, Japan). Briefly, cells were seeded into 96-well plates (1 × 10³ cells per well) and cultured for 4 h until adherence. The CCK-8 agent was added in each well at the indicated time-point, and the optical density at 450 nm was assessed after 1 h using a plate reader.

The transwell invasion assay was carried out using six-well plates containing transwell inserts (8-μm pore size; BD Biosciences) according to the manufacturer’s guidelines. Matrigel purchased from BD Biosciences was added to serum-free media, transferred to the top chamber, and incubated for 5 h. Subsequently, cells were cultured in the top chamber supplemented with serum-free medium. The lower chamber was supplemented with 10% FBS, and cells were removed from the top chamber after 36-h incubation. For quantification of the cells in the lower chamber, membranes were PFA-fixed, stained with crystal violet, and invading cells were quantified using microscopy image analysis.

Data analysis was carried out using R (version 3.6.3) and GraphPad Prism (version 8.0.1). Gene expression levels between HB and normal liver tissue samples were compared using the Student’s t-test. To evaluate the predictive value of each hub gene for the distinction between HB and normal liver tissue, we applied the receiver operating characteristic (ROC) curve. An area under curve (AUC) > 0.90 and P < 0.05 indicated statistical significance.

A detailed outline of our study is summarized in Figure 1. For our analysis, we combined two publicly available microarray gene expression datasets of HB and normal liver tissue samples. We carried out PCA to visualize data before and after batch effect correction, during which four outlier samples (GSM1948577, GSM1948562, GSM1948566, and GSM3770543) were removed (Figures 2A–C), resulting in a total of 99 HB samples and 19 normal liver samples after data preprocessing and quality control. We applied a filtering step (P value < 0.05 and |log2FC| > 1) for the identification of DEGs, which resulted in a total of 856 DEGs. Among these DEGs, 350 were up-regulated, while 506 were down-regulated, with a volcano plot presented in Figure 2D. The heatmap of the top 100 genes is shown in Figure 2E. We conducted GO and KEGG pathway analyses to elucidate the biological functions and potential signaling pathways these genes may be involved in. GO analysis suggested that DEGs predominantly consisted of genes involved in small molecule catabolic processes, the collagen-containing extracellular matrix, and coenzyme binding (Figures 3A–C). KEGG analysis identified enrichment for PI3K-AKT signaling, cell cycle, and FoxO signaling (Figure 3D). Critical pathways, including the cell cycle, FoxO pathway, NF-kappa B signaling pathway, amoebiasis, and carbon metabolism, are presented along with their related genes in Figure 3E. For the term cell cycle, all enriched DEGs were up-regulated with the exception of GADD45B and GADD45D. It should be noted that we observed two genes concurrently enriched in three critical pathways: TGFB2 was enriched in the cell cycle, FoxO signaling, and amoebiasis pathways, while GADD45B was enriched in the cell cycle, FoxO signaling, and NF-kappa B signaling pathways.

A PPI network containing 791 nodes and 9,054 edges was conducted using the Cytoscape software based on results of the STRING online database (Figure 4A). All four methods within the CytoHubba plugin were adopted, and the top 20 genes of each method were listed (Table 1). The intersecting genes that were concurrently listed in the four methods were regarded as hub genes (AURKA, AURKB, CDK1, CCNA2, CDC20, and PLK1) for PPI analysis. Nineteen clusters were obtained after module analysis using the MCODE plugin of Cytoscape, and we selected the top three modules as hub modules based on MCODE scores (Figures 4B–D). Notably, all six hub genes were found in module 1, which played an essential role in the constructed PPI network. Specifically, module 1 contained 59 nodes and 1,600 edges and had the highest MCODE score (55.172) of all modules. Another notable observation from module analysis was that all genes from module 1 exhibited up-regulation. Subsequently, we conducted GO and KEGG analyses of genes in module 1 using the R clusterProfiler package. For BP within the GO analysis, we found that genes in module 1 played a critical role in nuclear division, organelle fission, cell cycle transition, mitotic cell cycle transition, as well as chromosome segregation (Figure 5A). For CC within the GO analysis, we found that up-regulated genes were significantly enriched in the chromosomal region, condensed chromosome, and spindle (Figure 5B). The MF of GO analysis showed that genes were associated with ATPase activity, catalytic activity, action on DNA, protein serine/threonine kinase activity, and single-stranded DNA binding (Figure 5C). KEGG analysis showed that genes from module 1 were enriched for the cell cycle, DNA replication, as well as oocyte meiosis pathways (Figure 5D).

During sample clustering, two samples were regarded as outliers and thus excluded (GSM1948574 and GSM3770517; Supplementary Figure 1A). Besides, we identified β = 10 and R² = 0.88 as the optimal soft threshold parameters to guarantee a scale-free network (Supplementary Figures 1B,C). We set clustering height cut-off to 0.25 in order to merge similar modules, which resulted in seven modules (Figure 6A). Specifically, blue, black, brown, pink, green, and magenta modules were identified as significant modules (Figure 6B). The blue module containing 259 genes appeared to be the most relevant module involved in HB. The top 100 genes of the blue module, ranked by gene significance for cancer, are listed in Supplementary Table 1. Subsequently, the module eigengenes and associations between eigengenes and sample types were computed. The module eigengene dendrogram was plotted, and the seven modules were divided into two clusters. Similar results were obtained from eigengene network heatmap (Figure 6C). Interestingly, the blue module was not only located close to cancer but also had a markedly positive association with cancer, meaning that genes in the blue module may be essential for tumor progression. Moreover, module–trait relationship analysis confirmed the highly positive correlation between the blue module and cancer (r = 0.64, P = 1e-14) (Figure 6D). When focusing on the blue module (Figure 6E), we substantiated a significantly positive association between MM and GS (r = 0.64, P = 6e-38). Consequently, the blue module was chosen for functional enrichment analysis, during which we aimed to elucidate potential biological processes involved in HB.

In order to explore potential genes and pathways associated with HB growth, we conducted GO and KEGG analyses on the blue module identified by WGCNA. KEGG analysis indicated that genes in the blue module were markedly enriched for the cell cycle, oocyte meiosis, and DNA replication pathways (Figure 6F). Key pathways and their associated genes are shown in the heatmap (Figure 6G). Additionally, GO analysis demonstrated that genes in the blue module were primarily associated with organelle fission, nuclear division, chromosomal region, tubulin binding, and ATPase activity (Figures 7A–C). For the description of functionally enriched GO clusters, we utilized cnetplots to highlight the relationships between genes and critical pathways (Figures 7D–F).

Several key genes identified using PPI analysis were also included in the WGCNA blue module, with the intersecting genes being AURKA, AURKB, CDK1, CCNA2, and CDC20. The scoring of each hub gene in PPI and WGCNA is summarized in Table 2. For further validation of these potential key genes, we compared their expression values between HB and normal liver samples in the GSE75271 and GSE131329 datasets. Expression levels of these five key genes were markedly elevated in HB samples compared with normal liver samples (Figure 8A). ROC curve was utilized to evaluate the predictive value of each hub gene for the distinction between HB and normal liver tissue. The AUC of expression levels for four of the genes exceeded 0.90 in the ROC analysis (Figure 8B). Specifically, the AUC was 0.918 (95% CI, 0.865–0.970) for AURKA, 0.964 (95% CI, 0.933–0.994) for CDK1, 0.952 (95% CI, 0.915–0.990) for CCNA2, and 0.928 (95% CI, 0.882–0.973) for CDC20. A literature search revealed AURKA as a previously reported oncogenic gene in HB, and elevated expression levels of AURKA have been associated with an advanced COG stage as well as metastasis (Zhang et al., 2018; Tan et al., 2020). However, the role of CDK1, CCNA2, or CDC20 in HB growth has not been reported to date, and thus, we selected these three genes for subsequent experimental validation.

In order to investigate the influence of CDK1, CCNA2, or CDC20 expression in HB cells, we knocked down CDK1, CCNA2, or CDC20 in HepG2 and HuH-6 cells via siRNAs. The knockdown efficiency of each hub gene was validated by WB analysis (Figure 9A). After transfection of CDK1-siRNA into HepG2 and HuH-6 cell lines, the effect of CDK1 knockdown on cell proliferation was explored using both a CCK-8 assay (Figure 9B) and a colony formation assay (Figures 9C,D). These assays indicated a significantly lower proliferative ability of the CDK1-siRNA group compared to the control siRNA group. Similar effects were observed for CCNA2 or CDC20 knockdown in HB cells (Figures 9B–D). Next, we evaluated the effect of CDK1, CCNA2, or CDC20 knockdown on the invasive ability of HB cells using a transwell invasion assay (Figures 9C,D), which revealed a significantly decreased rate of the relative invasive cells relative to controls for all three knockdown models. Lastly, wound-healing assays revealed that CDK1, CCNA2, or CDC20 siRNA knockdown groups exhibited a markedly lower relative migration distance than the control did (Figures 9E,F). Taken together, these results demonstrated that knockdown of CDK1, CCNA2, or CDC20 inhibited proliferative, migratory, and invasive capabilities in both HepG2 and HuH-6 cell lines.