HB-related data were obtained from the Gene Expression Omnibus (GEO¹) database portal via the keyword “hepatoblastoma.” The dataset was included when all four items of the following criteria were met: (1) there were both HB and normal liver samples; (2) the dataset had miRNA or mRNA transcriptome data; (3) data for all samples were completely presented; (4) HB and normal liver samples could be clearly distinguished using principal component analysis (PCA). After screening, we chose one mRNA dataset (accession number: GSE131329) and one miRNA dataset (accession number: GSE153089) for further analysis. GSE131329 (Hiyama et al., 2019), consisting of 53 HB samples and 14 normal liver samples, was analyzed via GPL6244 platform (Affymetrix Human Gene 1.0 ST Array). GSE153089 (Honda et al., 2020), comprising of 30 HB samples and 14 normal liver samples, was analyzed via GPL21572 platform (Affymetrix Multispecies miRNA-4 Array). General information of the two datasets used for the present study is shown in Supplementary Table 1.

The GSE153089 dataset included nine specimens from metastatic tumor, 21 specimens from primary tumor (11 fetal subtypes and 10 embryonal subtypes), and 14 specimens from surrounding normal liver (Honda et al., 2020). Due to the lack of metastatic tumor samples in the GSE131329 dataset (Hiyama et al., 2019), we excluded all specimens from metastatic tumor in the GSE153089 dataset before further analysis. The remaining specimens in the GSE153089 dataset were subsequently divided into three groups, namely, normal surrounding liver, fetal-type tumor, and embryonal-type tumor groups. Each patient in the GSE153089 dataset possessed no more than one specimen from the same group except for patient 7. There were two fetal-type tumor specimens for patient 7 (Sample ID: 25F-1 and 25F-2), one (Sample ID: 25F-2) of which was randomly excluded for further analysis. Detailed information of samples in the GSE153089 dataset used for the present study is listed in Supplementary Table 2.

Raw data files (^∗.CEL) of GSE153089 and GSE131329 were imported using the oligo (Carvalho and Irizarry, 2010) R package. The data were sequentially filtered, background corrected, log base 2 transformed, and normalized. Based on the platform annotation information, gene symbol was obtained via conversion of the probe. If one gene symbol corresponded to two or more probes, the mean expression level of these corresponding mRNAs or miRNAs was treated as the final expression value. Before and after clustering and removing outliers, we detected the distribution patterns of HB and normal liver samples via PCA. DE-mRNAs and DE-miRNAs were then detected using the limma R package (17). An adjusted P < 0.05 and | log2FC| > 1 indicated statistical significance. Benjamini–Hochberg (BH) method was used to adjust the P value. Regarding the GSE131329 dataset, DE-mRNAs were obtained based on the comparison between HB and normal liver samples. As for the GSE153089 dataset, the first cluster of DE-miRNAs were detected between fetal-type tumor and normal liver groups, while the second cluster of DE-miRNAs were detected between embryonal-type tumor and normal liver groups. Through the intersection of these two clusters of DE-miRNAs, the upregulated or downregulated hub miRNAs were obtained.

Based on the hub miRNAs, we predicted the upstream TFs via TransmiR v2.0 (Tong et al., 2019), an easy-accessible public tool integrating experimentally verified TF-miRNA regulatory relationships from the publications. The Cytoscape software was subsequently utilized to visualize TF-miRNA regulatory relationships (Shannon et al., 2003). In addition, miRNet database was used for the prediction of the downstream target genes of hub miRNAs (Fan et al., 2016).

To further explore functional annotation of the candidate genes, we performed Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses via the clusterProfiler R package (Yu et al., 2012). GO terms consisted of biological process (BP), cellular component (CC), and molecular function (MF). An adjusted P < 0.05 was considered significantly enriched, and BH method was used to adjust the P value.

To unveil the relationships between the candidate genes, we established the PPI network via the STRING database (Szklarczyk et al., 2015). PPI pairs were considered significant with a combined score ≥ 0.4. Cytoscape software was subsequently applied to visualize the network (Shannon et al., 2003). On the basis of the degree obtained through Cytoscape plugin CytoHubba (Chin et al., 2014), top 10 hub genes were detected in the PPI network.

We applied Oncopression database² to validate expression levels of top 10 up-regulated hub genes and top 10 down-regulated hub genes. Oncopression is a web-based integrated gene expression profile using single sample normalization method UPC (Lee and Choi, 2017).

Hepatoblastoma and matched normal liver tissue samples from eight children undergoing surgical excision for primary HB were obtained from our hospital between 2014 and 2019. None of the patients received adjuvant radiotherapy or chemotherapy prior to surgery. Tissues were stored at −80°C immediately after harvest until further use. The pathological diagnosis of the tissue adjacent to each frozen tissue specimen was confirmed by at least two independent pathologists.

Total RNAs were isolated from the tissues using TRIzol reagent (Life Technologies, Carlsbad, CA, United States). Total mRNA was subsequently reverse-transcribed to produce complementary DNA (cDNA) using TaKaRa reverse transcription kit (TaKaRa Bio, Shiga, Japan). The SYBR Green fluorescence system (Roche, IN, United States) was used, and mRNA qRT-PCR was performed using a quantitative mRNA kit (TaKaRa Bio, Shiga, Japan). Based on the 2^–ΔΔCt method, the relative mRNA levels were normalized to GAPDH mRNA levels. All primers were synthesized by Sangon (Shanghai, China). The sequence of primers is summarized in Supplementary Table 3.

We conducted data analysis and visualization using R software (version 3.6.3) and GraphPad Prism (version 8.0.1). The expression levels of mRNAs or miRNAs between groups in the datasets were compared via a moderated t-test. For differential expression analysis of mRNAs or miRNAs in the datasets, a P value < 0.05 and | log2FC| > 1 were considered statistically significant. The mRNA expression levels of hub genes in HB and matched normal liver tissues from our hospital were statistically analyzed by a paired Student t test, and P values below 0.05 were considered significant.

The expression levels of fetal-type tumor and normal liver samples in the GSE153089 dataset prior to and after normalization are shown (Supplementary Figures 1A,B). PCA results before and after removing outliers (GSM4633970, GSM4633988, and GSM4633998) are also presented (Supplementary Figures 1C,D). Based on the differential expression analysis, we detected 41 upregulated DE-miRNAs and 36 downregulated DE-miRNAs, which are presented via volcano plot in Figure 1A. In addition, these DE-miRNAs between fetal-type tumor and normal liver samples were regarded as the first cluster of DE-miRNAs.

The expression values of embryonal-type tumor and normal liver samples in the GSE153089 dataset prior to and after normalization are shown (Supplementary Figures 2A,B). PCA results before and after removing outliers (GSM4633970 and GSM4633988) are also presented (Supplementary Figures 2C,D). Through the differential expression analysis, we detected 37 upregulated DE-miRNAs and 33 downregulated DE-miRNAs, which are presented via volcano plot in Figure 1B. Additionally, these DE-miRNAs between embryonal-type tumor and normal liver samples were regarded as the second cluster of DE-miRNAs.

Through the intersection of the aforementioned two clusters of DE-miRNAs, a total of 33 upregulated DE-miRNAs and 12 downregulated hub miRNAs were obtained (Figures 2A,B). Detailed information of respective hub miRNAs is also listed (Supplementary Tables 4,5).

As for the upregulated hub miRNAs, the predicted TFs included HNF4A, GTF2I, CEBPB, CREB1, MAZ, NR3B3, SHP, KLF4, PKM, and EED (Figure 2C). Regarding the downregulated hub miRNAs, the predicted TFs included NR3B3, SHP, CDKN1A, KLF3, USP7, and HSF2 (Figure 2D). Detailed information of the TFs predicted for the upregulated or downregulated DE-miRNAs is also listed (Supplementary Tables 6,7). Apart from the predicted TFs, we also predicted 5,772 target genes of the upregulated hub miRNAs and 5,600 target genes of the downregulated hub miRNAs. Upregulated hub miRNA-target gene network and downregulated hub miRNA-target gene network are presented in Figures 2E,F, respectively.

The expression levels of all samples in the GSE131329 dataset before and after normalization are visualized in Supplementary Figures 3A,B, respectively. PCA results prior to and after excluding the outlier (GSM3770543) are also shown (Supplementary Figures 3C,D). We then obtained 580 upregulated DE-mRNAs and 790 downregulated DE-mRNAs, which are presented via volcano plot in Figure 3A. The detailed information of these respective DE-mRNAs is listed (Supplementary Tables 8,9). Subsequently, we intersected the target genes of upregulated hub miRNAs and downregulated DE-mRNAs, resulting in a total of 250 downregulated candidate genes (Figure 3B). In addition, the intersection of target genes of downregulated hub miRNAs and upregulated DE-mRNAs resulted in 202 upregulated candidate genes (Figure 3C). Detailed information of these downregulated and upregulated candidate genes is also listed (Supplementary Tables 10, 11).

Biological process analysis indicated that enriched GO terms for downregulated candidate genes included response to nutrient levels, response to metal ion, small molecule catabolic process, and steroid metabolic process (Figure 4A). CC analysis showed that the candidate genes were markedly enriched in collagen-containing extracellular matrix, mitochondrial matrix, vesicle lumen, cytoplasmic vesicle lumen, secretory granule lumen, and blood microparticle (Figure 4C). In the process of MF analysis, the candidate genes were markedly enriched in coenzyme binding, cytokine activity, heme binding, tetrapyrrole binding, and oxidoreductase activity (Figure 4E). The complex relationships between these candidate genes and their related GO terms were visualized using the cnetplot R package (Figures 4B,D,F). Moreover, KEGG analysis identified complement and coagulation cascades, TNF signaling pathway, mineral absorption, and valine, leucine and isoleucine degradation as markedly enriched pathways (Figure 4G). Besides, the enriched pathways and their associated candidate genes were unveiled, which are shown as heatmap in Figure 4H.

We then conducted GO terms analysis based on the upregulated candidate genes. BP analysis revealed that nuclear division, organelle fission, mitotic nuclear division, and chromosome segregation served as the top enriched GO terms (Figure 5A). CC analysis identified chromosomal region, spindle, microtubule, and condensed chromosome as significantly enriched GO terms (Figure 5C). In addition, MF analytic results revealed microtubule binding, single-stranded DNA binding, DNA-dependent ATPase activity, DNA helicase activity, and histone kinase activity as markedly enriched GO terms (Figure 5E). The complex relationships between the aforementioned enriched GO terms and their associated candidate genes are also shown (Figures 5B,D,F). KEGG analysis identified cell cycle, cellular senescence and PI3K-AKT signaling as significantly enriched pathways (Figure 5G), and the associations of these pathways and their related candidate genes were visualized using heatmap (Figure 5H).

The downregulated or upregulated candidate genes were loaded into the STRING database, resulting in the construction of respective PPI network. A total of 211 node pairs were obtained for the downregulated candidate genes (Figure 6A), while 157 node pairs were obtained for the upregulated candidate genes (Figure 6C). The node pairs were input into Cytoscape software to visualize genes in respective PPI network. The respective top 10 hub genes were detected via Cytoscape plugin CytoHubba (Figures 6B,D). Specifically, the top 10 upregulated hub genes were CDK1, CCNB1, KIF11, PLK1, NCAPG, TOP2A, AURKA, TP53, ASPM, and TPX2, while the top 10 downregulated hub genes included IL6, DECR1, EGFR, CXCL8, CAT, IGF1, IL1B, F2, PTGS2, and FOS.

Oncopression database was utilized to validate the expression values of respective top 10 hub genes in the PPI network. As shown in Figure 7, eight of the top 10 downregulated hub genes (IL6, EGFR, CXCL8, CAT, IGF1, IL1B, PTGS2, and FOS) had significantly lower expression levels in HB tissue samples compared to normal liver tissue samples, while nine of the top 10 upregulated hub genes (CDK1, CCNB1, KIF11, NCAPG, TOP2A, AURKA, TP53, ASPM, and TPX2) had markedly higher expression levels in HB tissues in comparison to normal liver tissues.

Through literature search, we identified AURKA (Zhang Y. et al., 2018; Tan et al., 2020) and CDK1 (Tian et al., 2021) as previously reported oncogenic genes in HB. In contrast, the role of the other 15 hub genes in HB has not been reported to date or remains controversial. Based on HB and matched normal liver tissue samples in eight children with HB, the mRNA expression levels of these 15 hub genes were validated by qRT-PCR. The expression levels of EGFR, CAT, IGF1, IL1B, PTGS2, and FOS were significantly lower for HB tissues when compared with normal liver tissues. On the other hand, the expression values of CCNB1, KIF11, NCAPG, TOP2A, TP53, ASPM, and TPX2 in HB tissues were significantly higher than those in the normal liver tissues (Figure 8). Lastly, according to the predicted miRNA–mRNA pairs and the final verification results, we constructed the potential miRNA–mRNA interaction network involved in HB (Figure 9).