Abstract A lack of relevant genetic models and cell lines hampers our understanding of hepatoblastoma pathogenesis and the development of new therapies for this neoplasm. Here, we report an improved MYC-driven hepatoblastoma-like murine model that recapitulates the pathological features of embryonal type of hepatoblastoma, with transcriptomics resembling the high-risk gene signatures of the human disease. Single-cell RNA-sequencing and spatial transcriptomics identify distinct subpopulations of hepatoblastoma cells. After deriving cell lines from the mouse model, we map cancer dependency genes using CRISPR-Cas9 screening and identify druggable targets shared with human hepatoblastoma (e.g., CDK7, CDK9, PRMT1, PRMT5). Our screen also reveals oncogenes and tumor suppressor genes in hepatoblastoma that engage multiple, druggable cancer signaling pathways. Chemotherapy is critical for human hepatoblastoma treatment. A genetic mapping of doxorubicin response by CRISPR-Cas9 screening identifies modifiers whose loss-of-function synergizes with (e.g., PRKDC) or antagonizes (e.g., apoptosis genes) the effect of chemotherapy. The combination of PRKDC inhibition and doxorubicin-based chemotherapy greatly enhances therapeutic efficacy. These studies provide a set of resources including disease models suitable for identifying and validating potential therapeutic targets in human high-risk hepatoblastoma. Subject terms: Paediatric cancer, Cancer models, Cancer genetics, Cancer therapy, Liver cancer __________________________________________________________________ The availability of relevant animal models that can recapitulate high-risk hepatoblastoma will help to better understand its pathogenesis. Here the authors report and characterize a hepatocyte-specific, MYC-driven hepatoblastoma mouse model and show it recapitulates the human hepatoblastoma pathophysiology. Introduction Hepatoblastoma and hepatocellular carcinoma (HCC) are the most common primary liver malignancies in children and adolescents/young adults. While primary liver cancers account for only 1–2% of all pediatric tumors^[80]1, the largest incidence increase has been observed for hepatoblastoma in children under 5 years in nearly all regions of the world^[81]2. The rate is increasing at more than 4.3% annually in the US^[82]3. Hepatoblastoma is an embryonal neoplasm that likely arises from hepatic cell precursors^[83]4,[84]5. Genetically, hepatoblastoma has the fewest somatic mutations among all human cancers^[85]6, suggesting that hepatic precursor cells during the early stage of liver development may be particularly susceptible to fundamental events resulting in oncogenic transformation. In line with previous findings as reviewed^[86]7,[87]8, genomic sequencing studies have confirmed that mutations in the Wnt-β-catenin signaling pathway are the most common genetic event in hepatoblastoma^[88]9–[89]17. The gene for the antioxidant transcription factor, NFE2L2, is also altered in a subpopulation of high-risk hepatoblastomas^[90]9,[91]10,[92]15,[93]17, suggesting that liver cells undergo oxidative stress during cellular transformation or disease progression. The Hippo signaling pathway plays a critical role in liver organogenesis and cancer^[94]18–[95]20. The dysregulated downstream effector molecule of Hippo signaling, YAP1, is involved in hepatoblastoma tumorigenesis^[96]21–[97]25. Combination of the activated form of YAP1 (YAP^S127A) with either hepatoblastoma relevant NFE2L2 mutant or CTNNB1 mutant promotes liver tumorigenesis although either alone is unable to transform normal liver cells in these mouse models^[98]26. Hepatic developmental pathways may determine the differentiation capacity of mutated liver progenitor/stem cells, and differentiation status may determine the aggressiveness of hepatoblastoma^[99]9. Hepatoblastomas with high expression levels of stem/progenitor cell markers (EpCAM, LIN28B, SALL4, HMGA2, AFP) are usually associated with poor prognosis^[100]9. Such liver stem/progenitor cells have the ability to accumulate mutations following chemotherapy, leading to the development of post-treatment residual disease resulting in relapse and metastasis^[101]17. The MYC oncogenes are involved in many cancers including hepatoblastoma^[102]11,[103]27–[104]31. Gain of chromosome 2 and 8 (with MYCN and MYC oncogenes, respectively) is common (25–50%) in hepatoblastoma^[105]7,[106]11,[107]32,[108]33. While β-catenin mutation alone (CTNNB1) is usually insufficient to transform liver progenitor cells into hepatoblastoma^[109]26, MYC cooperates with β-catenin and YAP to sustain tumorigenesis^[110]29 and is an essential requirement for tumor maintenance in a β-catenin-based hepatoblastoma mouse model^[111]28. β-catenin drives MYC expression^[112]30 and MYC silencing prevents tumor growth in human hepatoblastoma cancer cell line-based xenograft models^[113]11. These data indicate that MYC plays an essential role in hepatoblastoma growth. Due to a lack of targetable somatic mutations and a paucity of genetic animal disease models and cell lines^[114]34,[115]35, identification of therapeutic targets in hepatoblastoma remains challenging. Conventional chemotherapy is critical for most hepatoblastoma treatment. However, the genetic response of hepatoblastoma cells to chemotherapy is not well defined, which impedes development of more effective therapies because of an incomplete understanding of the mechanism of therapeutic response and resistance. Here we generate a hepatocyte-specific MYC-driven multifocal hepatoblastoma-like tumor model that resembles high-risk human hepatoblastoma. The transcriptomics of this transgenic hepatoblastoma-like model are characterized by bulk RNA-seq, single cell RNA-seq, and pathology-based spatial transcriptomics, all of which confirm its similarity with human hepatoblastoma. Cell lines generated from this model are readily passaged in vitro. Cancer dependency genes are mapped by a genome-wide CRISPR-Cas9 screening approach. We also perform genetic mapping of cellular responses to doxorubicin, a commonly used chemotherapeutic for hepatoblastoma treatment, with a genome-wide CRISPR-Cas9 screen and identify genes that synergize with and antagonize the effect of chemotherapy. Based upon this screen, a combination therapy is developed which shows better efficacy than doxorubicin treatment alone. Our studies characterize hepatoblastoma disease models (mouse and cell lines) that recapitulate pediatric hepatoblastoma and identify potential therapeutic targets of hepatoblastoma that are conserved across mouse and human species. Results Hepatocyte-specific MYC overexpression drives rapid hepatic oncogenesis Previous genetic hepatoblastoma mouse models have provided invaluable information toward our understanding of the role of oncogenic drivers in this cancer. However, most of these models have only addressed well differentiated hepatoblastoma, which has a relatively good clinical outcome, or they do not align with the onset of liver development in children. To overcome these limitations, we generated a model in a C57BL/6J genetic background by crossing hepatocyte-specific transgenic Alb-Cre mice (Cre recombinase under the control of the mouse albumin enhancer/promoter hybrid)^[116]36 with CAG-STOP^flox/flox-Myc mice (CAG promoter-driven human c-MYC, whose expression is prevented by a LoxP site flanked STOP cassette)^[117]37 (Fig. [118]1a). Hepatocyte-specific, Cre-mediated excision of the floxed STOP cassette allows expression of the CAG promoter-driven human Myc, leading to a typical phenotype with hepatomegaly and paraneoplastic alopecia in double transgenic Alb-Cre;CAG-Myc mice (ABC-Myc, Figs. [119]1b and [120]S1a). Strikingly, activation of one allele of the Myc oncogene led to rapid onset of liver tumors in neonatal mice; all these mice died within 1–10 weeks after birth (Fig. [121]1c). which is consistent with the known role for MYC in sustaining hepatoblastoma growth. Embryonic lethality was not induced by Myc activation as all possible genotypes were recovered at the expected Mendelian ratio (Table [122]S1). Western blot and immunohistochemistry validated MYC overexpression in non-neoplastic hepatocytes and liver tumor tissues at fetal (E17.5) and different postnatal stages (Figs. [123]1d, [124]S1b), suggesting that MYC is activated in the fetal stage. In parallel, we also developed an ABC-Myc;TdTomato model that had a similar tumor penetrance and lethality but that also expressed TdTomato as a lineage reporter (Fig. [125]1c). Together the data demonstrate that the introduction of human MYC alone is sufficient to quickly drive tumorigenesis in susceptible murine hepatic stem/progenitor cells in the fetal livers of transgenic mice. Fig. 1. ABC-Myc drives hepatoblastoma-like tumor development. [126]Fig. 1 [127]Open in a new tab a Breeding strategy to generate Alb-Cre;CAG-Myc (ABC-Myc) compound mice. b Hepatomegaly with tumor nodules in ABC-Myc liver in comparison with age matched normal mouse liver. c Inferior overall survival of ABC-Myc (green color, n = 33), and ABC-Myc/TdTomato (pink, n = 11) mice, respectively. Log-rank (Mantel-Cox) test used for statistical analysis in Kaplan-Meier survival. d Western blot showing overexpression of C-MYC in ABC-Myc livers at postnatal day 5 (P5), 10 (P10),15 (P15) and 36 (P36) in comparison with the normal controls. The blots are representative of three independent experiments. e Hematoxylin and Eosin (H&E) shows histology of ABC-Myc tumors at postnatal day 7 (P7), 25 (P25), 67 (P67). Sample number for each image n = 1. Scale bar = 25 μm. f Hematoxylin and Eosin (H&E) shows mixed histology of human and ABC-Myc tumors. Sample number for each image n = 1. Scale bar = 25 μm. g Hematoxylin and Eosin (H&E) shows pleomorphism of human and ABC-Myc tumors. Sample number for each image n = 1. Scale bar = 25 μm. h Immunostaining of alpha fetoprotein (AFP), glutamine synthetase (GLUL), Spalt Like Transcription Factor 4 (SALL4), Glypican 3 (GPC3), arginase (ARG1), β-catenin (BCAT), cytokeratin 19 (KRT19) and integrase interactor 1 (INI1). Sample number for each image n = 1. Scale bar = 25 μm. Source data are provided as a Source Data file. Pathological analyses define the ABC-Myc-driven liver neoplasm as a hepatoblastoma-like malignancy Hepatoblastoma is histologically heterogenous, with two main histologic patterns (epithelial, and epithelial mixed with mesenchymal components). Epithelial patterns are further delineated into fetal, embryonal, mixed fetal and embryonal, cholangioblastic, small cell undifferentiated, macrotrabecular, mixed and others^[128]38,[129]39. Tumors arising in the ABC-Myc model effaced most of the sampled liver tissues and had a highly resembling human hepatoblastoma histology with embryonal or combinations of both fetal and embryonal morphologies comprising the bulky tumors, as well as scattered foci of extramedullary hematopoiesis (Fig. [130]1e, Table [131]1). These multifocal tumors involving all liver lobes correspond to human PRE-Treatment EXTent of tumor (PRETEXT) stage IV disease^[132]40. Clinically, 35% of patients with hepatoblastoma present as multifocal tumors at diagnosis, and 43% of these are PRETEXT stage IV^[133]41, a poor prognostic factor that usually requires high-intensity, dose-dense cisplatin and doxorubicin-based chemotherapy, and often liver transplantation^[134]42. Samples at time points E14.5, E17.5, and postnatal day 7 (P7) were evaluated to assess the presence of pre-neoplastic lesions. Neoplastic transformation was first observed in E17.5 livers in low numbers of scattered developing hepatocytes with abnormal nuclear morphologies (Fig. [135]S1c). Nuclear changes consisted of karyomegaly, marginalization of chromatin, and a single, centralized, and prominent nucleolus that is consistent with other cancers where constitutive MYC activation is present. These dysplastic cells were interpreted as pre-neoplastic lesions based on the biological time course of the ABC-Myc mouse model described in this paper. Table 1. Pathological characterization of ABC-Myc tumors Accession Number Fetal Embryonal Crowded Macrotrabecular Cholangioblastic Mesenchymal INI1 SALL4 GPC3 AFP GLUL ARG Bcat Grouping RS19-2057 + + - - +/15–20% - +/100%; retained +/90% +/100% +/100% +/100% + +/100%; membraneous C2 RS22-578 - + - - +/<5% - +/100%; retained + + + +/patchy & strong + +/100%; membraneous C2 RS22-575 + + - + +/<5% - +/100%; retained +/15-20% + + +/87% + +/100%; membraneous C2 AP21-532 + + + + +/<5% - +/100%; retained + + + + + +/100%; membraneous C2 RS22-576 - + - - - - +/100%; retained + + + + + +/100%; membraneous C2 RS22-577 + + - - +/<5% - +/100%; retained + + + + + +/100%; membraneous C2 [136]Open in a new tab Neoplastic nodules were grossly visible in all liver sections of ABC-Myc mice starting at P7. Multifocal to coalescing neoplastic foci with an embryonal morphology could be observed in the livers of P7 ABC-Myc mice (Fig. [137]1e), consistent with the hypothesis that hepatoblastoma-like neoplasia may arise from epithelial-lineage committed hepatic stem progenitor cells with the introduction of human oncogenic Myc signaling resulting in impaired differentiation. Further evaluation of the hepatoblastoma-like tumors from time points P25 to P67 showed a coexistence of distinct subpopulations of neoplastic cells with embryonal, fetal, and rarer cholangioblastic-like morphologies (Fig. [138]1e). The co-existence of these morphologies in advanced hepatoblastoma-like tumors is most consistent with human pediatric hepatoblastoma with a mixed epithelial phenotype. Small cell undifferentiated, rhabdoid, teratoid, and mesenchymal morphologies were not observed. There were no definitive well-differentiated fetal morphologies identified in the sections except within P67 tumors (Fig. [139]1e). All tumors had combinations of primitive morphologies comparable to the previously described C2 morphologic phenotype as described by Cairo et al.^[140]11. Hepatoblastomas characterized as C2 are documented to have aggressive biological behavior and an unfavorable prognosis, which is observed in this model. The Myc-driven murine hepatoblastoma-like tumors demonstrate phenotypic plasticity of hepatocyte lineage committed stem/progenitor cells. While the co-existence of embryonic and fetal histological features of ABC-Myc tumors resemble the human hepatoblastoma (Fig. [141]1f), it is important to differentiate hepatoblastoma from hepatocellular carcinoma in pediatric patients, because of differing treatment and prognosis^[142]38. While the poorly differentiated histology is consistent with the pediatric C2 phenotype, some tumor areas also contain histologic features of the subclassification of pediatric hepatoblastomas with hepatocellular carcinoma features that were previously called transitional liver cell tumors (TLCT)^[143]9 (Fig. [144]1g), indicating that some ABC-Myc tumor cells have features of HCN-NOS (Hepatocellular Malignant Neoplasm, Not Otherwise Specified) that frequently presents phenotypic plasticity. We further determined the pathological features of this hepatoblastoma-like malignancy using immunohistochemical markers of human pediatric hepatoblastoma, and observed overexpression of hepatic stem/progenitor cell markers documented in C1 and C2 human pediatric hepatoblastomas^[145]11 (Fig. [146]1h). Murine hepatoblastoma-like neoplasms had diffuse immunopositivity for alpha fetoprotein (AFP) and glypican 3 (GPC3), two stem cell markers used to distinguish neoplastic hepatocellular cells^[147]39,[148]43,[149]44, as well as immunoreactivity for glutamine synthetase (GLUL or named as GS), a β-catenin target and a marker of β-catenin activated hepatocytes^[150]38,[151]39, SALL4, another embryonal type of hepatoblastoma marker^[152]45,[153]46, and Arginase-1 (ARG-1), a marker used to distinguish primary hepatocellular tumors from metastatic tumors^[154]47. Immunoreactivity was visually observed in greater than 75% of the bulky hepatoblastoma-like neoplasms and staining intensity for all markers was visually graded as moderate to strong in staining intensity for all markers (Table [155]1). Rare subpopulations of poorly differentiated neoplastic cells, visually quantified at less than 1% of the neoplasm, had immunoreactivity for cytokeratin 19 (KRT19), a marker for biliary cancer or small-cell undifferentiated type hepatoblastoma^[156]38, as well as non-neoplastic, entrapped bile ducts. INI1 (SMARCB1) was retained in all neoplasms, further demonstrating the hepatocellular origin of these cells. The strong cytoplasmic staining of β-catenin may suggest an activation of the Wnt/β-catenin signaling pathway in these tumors (Fig. [157]1h). Ki67 staining showed that 3–6% of cells were positive (Fig. [158]S1d, e). In summary, the ABC-Myc hepatoblastoma-like model overall recapitulates the embryonal or mixed fetal and embryonal histologic features of human hepatoblastoma, with some bearing HCN-NOS features, and has anatomic and molecular characteristics of human disease highly associated with the high-risk C2 subtype^[159]11. C1 and C2 hepatoblastoma subclasses were initially defined by gene expression profiling and can be delineated by epithelial cell type, proliferation differences, and expression of stem cell markers that may be assayed by IHC. C1 and C2 features have been further correlated with the phase of hepatic differentiation in which susceptible lineage-committed subpopulations may undergo tumorigenesis. C2 tumors have molecular features of non-neoplastic murine liver at E11.5 and E12.5, while C1 tumors have features of hepatic differentiation in late and postnatal stages^[160]11. Some retention of both C1 and C2 characteristics by our murine hepatoblastoma-like neoplasms may result from differences in transgene copy number expression in the embryonal liver; this difference may affect the timing of malignant transformation in susceptible hepatocyte specified stem-progenitor populations starting at E9.5, when albumin expression can first be detected^[161]48. Neoplastic transformation was first observable in E17.5 livers in small subpopulations of atypical appearing cells by histology, suggesting that tumorigenesis is occurring along a continuum of time in this model that is based on the increasing expression of albumin into adulthood. Therefore, hepatoblastoma-like neoplasms with hybrid features of C2 and C1 may be expected in the model. Thus, the ABC-Myc hepatoblastoma-like model recapitulates the morphologic features of human hepatoblastoma, histologically most similar to the high-risk C2 class of hepatoblastoma, with some bearing HCN-NOS features, but also retains some immunohistochemical and molecular characteristics of low-risk C1 neoplasms. These subtypes may occur sequentially or randomly. Nevertheless, we only observed well differentiated hepatoblastoma in P67 while embryonal and cholangioblastic subtypes occur at an early time (P7, P25) (Fig. [162]1e), suggesting that there could be a sequential event during MYC-mediated cellular transformation that coopts with liver developmental program. C2 aggressive type may be derived from stem/progenitor cells in hepatoblast, and thus appeared at an early developmental stage, while the C1 type may be derived from a more differentiated cells at late developmental stage. Serum chemistry panel analysis reveals liver dysfunction of ABC-Myc mice similar to that of human hepatoblastoma To assess the liver function of ABC-Myc mice, we performed serum chemical analysis (Fig. [163]2a, b). Not surprisingly, ABC-Myc mice showed abnormal elevation of AFP (Fig. [164]2a), alkaline phosphatase (ALP), alanine transaminase (ALT), and total bilirubin (Fig. [165]2b), the three commonly used biomarkers of liver function, indicating that the livers in ABC-Myc mice are damaged. One clinical study showed that 80% of hepatoblastoma patients had abnormal levels of ALP and 12.5% had increased ALT^[166]49. As liver is the major organ that produces glucose, liver cancer can cause hypoglycemia. Indeed, the serum glucose levels in ABC-Myc mice were remarkedly reduced (Fig. [167]2b). The serum levels of creatinine and blood urea nitrogen (BUN) in ABC-Myc mice were also declined in comparison with the age-matched normal mice although not statistically different. While creatinine and BUN are the commonly used chemical markers to assess kidney function, liver cancer can lead to less production of creatinine, a break-down product of creatine in liver through transamination of amino acids. Low levels of BUN may indicate liver disease in the clinic due to less production of urea. However, the albumin and globulin levels seemed to be in the normal range, and no abnormal levels of common electrolytes (Sodium, Potassium, Calcium and Phosphorous) were observed (Fig. [168]2b). We further performed complete blood count (CBC) measurements to assess if ABC-Myc mice had developed additional complications. While white blood cell counts showed no difference between normal mice and ABC-Myc mice, the absolute number of circulating eosinophils tended to be increased although the difference was not statistically significant (Fig. [169]S2). However, the ABC-Myc mice developed microcytic anemia, as indicated by a reduction in the proportion of red blood cells (hematocrit, HCT%), amount of hemoglobin (HB), mean corpuscular volume (MCV) and mean corpuscular hemoglobin (MCH), increase in size variation (percentage of red cell distribution width, RDW%), but normal range of total number of red blood cells and mean corpuscular hemoglobin concentration (MCHC) (Fig. [170]2c). Thrombocytosis also occurred in ABC-Myc mice, as indicated by an increase in total platelet counts, plateletcrit (PCT), and mean platelet volume (MPV) (Fig. [171]2d). One study reported that among hepatoblastoma patients, 75% had thrombocytosis and 37.5% had microcytic anemia, whereas only 23.1% of pediatric patients with hepatocellular carcinoma had thrombocytosis and none had microcytic anemia^[172]50. These chemistry and CBC parameters are consistent with the presence of hepatoblastoma-like disease in ABC-Myc mice. Fig. 2. Clinical chemistry analysis of serum from ABC-Myc mice. [173]Fig. 2 [174]Open in a new tab a Quantification of serum AFP levels in normal (n = 4 biologically independent animals) and ABC-Myc (n = 6 biologically independent animals) mice by ELISA. Data are presented as mean ± SD. Unpaired two-sided t-test,**p = 0.0046. b Chemistry panel markers in determination of liver function, kidney function and electrolytes in serum from normal (n = 3 biologically independent animals) and ABC-Myc (n = 4 biologically independent animals) mice. Data are presented as mean ± SD. Unpaired two-sided t-test,****p < 0.0001, **p = 0.0056, *p = 0.0171, ns not significant. ALP Alkaline phosphatase, ALT Alanine transaminase, BUN Blood urea nitrogen. c Complete blood count to determine the changes in red blood cells in blood from normal (n = 3 biologically independent animals) and ABC-Myc (n = 4 biologically independent animals) mice. Data are presented as mean ± SD. Unpaired two-sided t-test,****p < 0.0001, ***p < 0.001, ns not significant. HCT hematocrit, RBC red blood cell, HB hemoglobin, MCV mean corpuscular volume, MCH mean corpuscular hemoglobin, MCHC mean corpuscular hemoglobin concentration, RDW red cell distribution width, RSD red cell standard deviation. d Complete blood count to determine the changes in platelets in blood from normal (n = 3 biologically independent animals) and ABC-Myc (n = 4 biologically independent animals) mice. Data are presented as mean ± SD. Unpaired two-sided t-test,****p < 0.0001, ***p < 0.001, **p = 0.0017. PLT platelet, MPV mean platelet volume, PDW platelet distribution width, PCT Plateletcrit. Source data are provided as a Source Data file. Signaling pathways in ABC-Myc tumor cells resemble those in human hepatoblastoma with a poor outcome To understand the molecular mechanisms of ABC-Myc hepatoblastoma-like tumors, we identified the differentially expressed genes in tumors versus age-matched normal murine livers using bulk RNA-seq (Fig. [175]3a and Supplementary Data [176]1), followed by signaling pathway analysis. Interestingly, the Igf2 oncogene ranked first (log2 fold change = 11.6) among the upregulated genes in tumors (Supplementary Data [177]1). In humans, IGF2 is located in the 11p15.5 imprinted locus, which is the second most frequently altered locus in hepatoblastomas and hepatocellular carcinomas, mostly through copy-neutral loss of heterozygosity^[178]17. IGF2 induction by 11p15.5 alterations is likely the first genetic event in hepatoblastoma^[179]17. The most significantly downregulated genes in tumors were cytochrome P450 (CYP) family genes related to metabolic functions of mature hepatocytes (Fig. [180]3a). Gene set enrichment analysis (GSEA)^[181]51 showed that the genes upregulated and downregulated in ABC-Myc tumors were significantly associated with the corresponding human hepatoblastoma gene signatures reported by Cario^[182]11 (Fig. [183]3b). Since Cario gene sets consisted of hepatoblastoma tissue RNA samples including those resected after preoperative chemotherapy, we compared ABC-Myc gene expression with the gene datasets generated from biopsy or surgery prior to any chemotherapy (Ikeda dataset, [184]GSE131329)^[185]52, which included 14 noncancerous liver tissues and 53 tumor tissues. We used the top 200 differentially expressed genes from Ikeda genset for GSEA analysis, and again, we obtained very similar results (Fig. [186]S3a), which further strengthened our conclusion. In agreement with the immunostaining findings, GSEA demonstrated that the β-catenin pathway was also significantly upregulated in ABC-Myc tumors as indicated by its association with gene signatures derived from β-catenin transgenic liver tumors^[187]29 and β-catenin knockdown in HepG2 cells^[188]53 (Fig. [189]3c). To further determine if ABC-Myc induces transcriptomes similar to those in human hepatoblastoma, we performed a comparative analysis using the VENN diagram showing the number of deregulated genes (and their %) in the comparison between tumor vs. non-tumor liver samples from hepatoblastoma patients and from the ABC-Myc model at FDR < 0.05 (Fig. [190]3d). Briefly, we integrated the RNA-seq from ABC-Myc tumors and control livers with the RNA-seq from Carrillo-Reixach’s study that included tumor and non-tumor samples from 32 patients with hepatoblastoma^[191]15. As a result, we obtained a matrix of 11,393 ortholog genes. Then, we performed a supervised analysis by comparing tumor vs. non-tumor samples using human and mouse samples. The results showed that 50.1% and 42.5% of the up- and down-regulated genes in the ABC-Myc tumors vs. control liver samples were also deregulated in human hepatoblastoma in comparison with non-tumor samples, respectively (Fig. [192]3d). The statistical analysis clearly showed a significant overlapping in upregulated (p = 1.6 × 10^−96) and downregulated (p = 2.1 × 10^−153) genes and clearly supports the high similarity of our ABC-Myc tumor model with human hepatoblastoma on the transcriptomic level. To further confirm the high similarity of molecular features between human and mouse tumor samples, we used the integrative human and mouse ortholog genes to perform a Principal Component Analysis of RNA-seq. The results showed that tumor and non-tumor samples were clearly grouped into two categories independent of the species from which samples were obtained (Fig. [193]3e). Specifically, mouse tumor samples were grouped with human tumor samples and control mouse liver samples were grouped with adjacent non-tumor samples from patients with hepatoblastoma. Additionally, we cross-referenced our RNA-seq results for the top 500 genes upregulated and downregulated in ABC-Myc tumors with human hepatoblastoma RNA-seq analysis reported by Hooks et al.^[194]54 (Fig. [195]S3b). The results again revealed that the top differentially expressed genes in ABC-Myc tumors were similarly altered in human hepatoblastomas (Fig. [196]S3b), further supporting that the murine hepatoblastoma-like model resembles human disease at the transcriptomic level. Further comparison of the gene pathways between ABC-Myc tumors and human hepatoblastomas revealed that both shared altered metabolic pathways and those regulating the cell cycle, DNA replication and repair, and RNA splicing (Table [197]S2). Altogether, our data clearly indicate the high similarity of our ABC-Myc model and human hepatoblastoma and support its use as an experimental model for this extremely rare disease. Fig. 3. Signaling pathways in ABC-Myc tumor cells resemble those in human hepatoblastoma with poor outcome. [198]Fig. 3 [199]Open in a new tab a Volcano plot showing differentially expressed genes in ABC-Myc tumors (n = 3) vs normal mouse livers (n = 3). X-axis represents the expression changes in log2 (fold). Y axis represents the significance of expression change for each gene in -log10 (p value). b GSEA showing genes highly downregulated and upregulated in ABC-Myc hepatoblastoma are significantly associated with the signatures downregulated (left panel) and upregulated (right panel) in human hepatoblastoma reported by Cario et al.^[200]11. P Value calculated by one-sided Fisher’s exact test. The FDR is calculated by comparing the distribution of normalized enrichment scores from many different genesets. c GSEA showing genes highly upregulated in ABC-Myc hepatoblastoma are significantly associated with the β-catenin signatures derived from mouse livers overexpressing β-catenin in dataset ([201]GSE79084)^[202]29 (left panel), and β-catenin knockdown in HepG2 cells from dataset ([203]GSE94858)^[204]53 (right panel). d Proportional VENN diagrams of the up-regulated (top) and down-regulated genes (bottom) in the human HB (n = 34) vs. adjacent non-tumor liver (NL, n = 32) samples (left) and mice Myc-ABC tumor (n = 3) vs. control liver (CL, n = 3) samples. The numbers in the Venn diagrams represent the number of significant genes at FDR < 0.05. The comparisons were performed considering the total of 11,393 ortholog genes. RNA-seq database from patients with HB was obtained from Carrillo-Reixach et al. ([205]GSE133039)^[206]15. The P values (upregulated p = 1.6 × 10^−96, downregulated p = 2.1 × 10^−153) of the overlaps are calculated by the hypergeometric distribution. e Principal Component Analysis using the integrated dataset consisting in 11393 genes present in mouse and human tumor (HB, [207]GSE133039)^[208]15, non-tumor liver (NL) and control liver (CL) samples. f GSEA showing stem cell gene signatures highly upregulated in ABC-Myc hepatoblastoma. P Value calculated by one-sided Fisher’s exact test. The FDR is calculated by comparing the distribution of normalized enrichment scores from many different genesets. g Pearson correlation heatmap using the dendrogram of bootstrapping hierarchical clustering from tumoral samples including the 11,393 ortholog genes present in mouse and human tumor samples. Human hepatoblastomas were annotated with molecular features obtained from Carrillo-Reixach et al. ([209]GSE133039)^[210]15. h Heatmap of the 11 ortholog genes of gene the 16-gene signature in C1 and C2 human hepatoblastomas ([211]GSE133039)^[212]15 and mouse ABC-Myc tumors. Source data are provided as a Source Data file. The outcomes of hepatoblastoma can be distinguished by two molecular signatures, C1 and C2, which represent better and worse outcomes, respectively^[213]11. We cross-referenced C1 and C2 signatures to our RNA-seq data and found that ABC-Myc tumors expressed higher levels of C2 and lower levels of C1 signatures (Fig. [214]S3c). To further validate that ABC-Myc tumors resemble the C2 class, we applied seven different prediction algorithms^[215]11, and all of which showed that ABC-Myc hepatoblastomas were classified as C2 (Table [216]S3). Hirsch et al. reported that hepatoblastoma can be further classified into 4 molecular subtypes, ‘Hepatic differentiation’, ‘Liver progenitor’, ‘Mesenchymal’ and ‘Proliferation’^[217]17. ABC-Myc tumors exhibited low expression of ‘Hepatic differentiation’ signature, but high expression of the ‘Liver progenitor’ and ‘Proliferation’ signatures (Fig. [218]S3d). We also found three out of six ‘Mesenchymal’ markers were expressed in ABC-Myc tumors (Fig. [219]S3d), albeit to a lesser degree (Supplementary Data [220]1). While the ‘Hepatic differentiation’ group overlaps with C1, the ‘Liver progenitor’ signature is associated with a subclass of hepatoblastoma that has the worst outcome^[221]15,[222]17. Consistent with the primarily embryonal histological features, GSEA results showed that ABC-Myc tumors had significant upregulation of cancer stem cell signatures including “liver cancer with upregulated EpCAM” and “undifferentiated cancer” (Fig. [223]3f). Compared with age-matched normal livers, the hepatoblastoma embryonic gene markers (i.e., Lin28b, Sall4, EpCAM, Hmga, Afp) in ABC-Myc tumors, which are usually associated with a poor outcome^[224]9, were increased over 4–250 fold (Fig. [225]3a, Supplementary Data [226]1). Notably, Lin28b is an oncogene that can drive hepatoblastoma in a transgenic mouse model^[227]55, and is highly expressed in high-risk hepatoblastoma^[228]31. Clustering analysis of the correlation between the gene expression of human hepatoblastoma in Carrillo-Reixach’s study and the mouse samples showed four main groups of tumor samples (Fig. [229]3g). In line with the results in Fig. [230]3e, the gene expression profile of ABC-Myc tumor samples was highly correlated with that of the human primary hepatoblastoma and specifically, with tumors in the proliferative “C2- Pure subclass” (p = 0.019) and with a strong “14q32-gene signature” overexpression (p = 0.027). These molecular features have been already reported to be associated with clinical features of poor prognosis^[231]11,[232]15. In line with these findings, Dlk1 was highly expressed in ABC-Myc tumors (log2 fold change = 8.85) (Supplementary Data [233]1). DLK1 is a well-known hepatoblast marker and is highly expressed in hepatoblastoma^[234]56. Carrillo-Reixach et al. recently identified the DLK1-DIO3 locus genes on 14q32 as a new hallmark of human hepatoblastoma that is associated with Wnt/β-catenin signaling, and high expression of 14q32 gene signature being associated with a poor outcome^[235]15, supporting that hepatoblasts could be the cells of origin of ABC-Myc hepatoblastoma. Finally, integrative analysis of the expression profile of the 11 ortholog genes of the 16-gene signature further confirmed that the mouse ABC-Myc tumors had a similar profile to the human C2 hepatoblastomas in the Carrillo-Reixach cohort^[236]11 (Fig. [237]3h). Taken together, these data demonstrate that ABC-Myc tumors resemble human hepatoblastoma with molecular signatures of aggressive disease. scRNA-seq analysis of ABC-Myc tumors reveals the heterogeneity of hepatoblastoma-like cells Single cell RNA sequencing (scRNA-seq) studies have shown that the mammalian liver is composed of multiple cell lineages in addition to hepatocytes and cholangiocytes^[238]57–[239]59. The heterogeneity of liver cells is further complicated by the anatomical structure of liver zonation^[240]60,[241]61, which shows a distinct expression pattern of metabolic genes distributed from the central vein to the portal vein along the lobule axis^[242]62. While scRNA-seq analysis has provided insight into adult hepatocellular carcinoma and its tumor microenvironment^[243]63,[244]64, tumor heterogeneity in hepatoblastoma at the single cell level has just been recently appreciated^[245]65,[246]66. To investigate if the transcriptomic ecosystem of ABC-Myc-driven tumors recapitulates human hepatoblastomas, we performed scRNA-seq to define the distinct cellular populations of cells dissected from 4 tumors and 3 healthy livers. Cell Ranger Single-Cell Software Suite (version 6, 10X Genomics) was used to quality control and quantify the single-cell expression data to generate filtered gene-barcode matrices for 90,715 cells with an average of 3037 mRNA molecules (UMIs, median = 1963, range: 302–32,765). First, we characterized the transcriptional differences between tumor samples and the control group by applying the NBID algorithm^[247]67 that we developed for differential analysis of scRNA-seq. The top 10 genes upregulated in tumor samples, ordered by the fold change (adjusted P < 6.428e-323, log2FC > 5), were Camp, Ngp, Igf2, Ltf, Prtn3, Afp, Ermap, Rhd, Elane, and Mpo (Fig. [248]S4a). The high levels of Igf2 and Afp further verified the hepatoblastoma-like tumors arising from ABC-Myc mice. Highly elevated expression of granule genes such as Camp, Ngp, and Ltf are correlated with neutrophil development. In addition, genes Prtn3, Elane, and Mpo are functional activation markers of neutrophils involved in inflammation, infection, and tumor invasion. Activated tumor-associated neutrophils release enzymes, mainly proteinase 3 (encoded by Prtn3), neutrophil elastase (encoded by Elane), and myeloperoxidase (coded by Mpo), destroying surrounding tissues, which may lead to tumor invasion. Numerous studies have shown that neutrophils in the tumor microenvironment can promote rapid tumor development and growth^[249]68. We also noted that erythroid genes Ermap and Rhd are among the top 10 genes, and erythroid-like signature is present in human pediatric hepatoblastomas^[250]65. Pathway analysis with Hallmark genes showed that DNA repair, cell cycle, MYC targets, E2F targets, and heme metabolism were enriched in tumor samples (adjusted P < 4.5e-08, Fig. [251]S4b). After correction for batch effect, 16 clusters of cells in normal and tumor tissues were generated by unsupervised clustering of the global single-cell transcriptomic datasets with Latent Cellular State Analysis (Fig. [252]4a, b). The frequencies of low-quality cells (cells with low UMI counts ( ≤ 500) or more than 20% UMI of mitochondrial genes) in clusters 5 were greater than 50%, which was therefore removed, and the rest of the 15 clusters were kept for further analysis. Among the 15 clusters that remained, six (clusters 2, 3, 7, 9, 12, and 16) were dominated by cells from tumor tissues, with more than 97.5% of cells in each cluster from a tumor sample (Fig. [253]4a, b). However, the cells in four clusters (clusters 4, 13, 14 and 15) were predominantly from the control group, with more than 93.6% of cells in each cluster from a normal liver sample. The remaining five clusters (clusters 1, 6, 8, 10, and 11) were largely shared by both tumor and control groups (with 23.5–60.6% of cells in each cluster from the tumor group). To define the biological functions of each cluster, we used Seurat (version 4.3.0) to compare the average single-cell expression profiles from each cluster with previously annotated reference datasets^[254]69. After determining the similarity of each cluster to previously defined cell types^[255]70, we generated top 10 markers for each cell cluster (Supplementary Data [256]2), and then performed gene set enrichment analysis to associate the marker genes of each group to known functional pathways in KEGG (Supplementary Data [257]2). The clusters shared between the tumor and control groups (clusters 1, 6, 8, 10, and 11) were all enriched with immunity related genes. Genes involved in DNA damage, heme biosynthetic process, and oxidative phosphorylation were enriched in tumor-specific clusters (clusters 2, 3, 7, 9, 12, and 16) (Supplementary Data [258]2). Clusters specific to the control group (clusters 13 and 15) were enriched for genes involved in the endoplasmic reticulum, essential for hepatocytes’ protein synthesis function. Fig. 4. scRNA-seq analysis of ABC-Myc tumors reveals the heterogeneity of hepatoblastoma-like cells. [259]Fig. 4 [260]Open in a new tab a A t-SNE plot showing the source of cells partitioned by normal livers (n = 3) and ABC-Myc tumor (n = 4) samples. b A t-SNE plot showing the inferred cell clusters using Latent Cellular State Analysis for the normal livers (n = 3) and ABC-Myc tumor (n = 4) samples. Cluster 5 mainly consisted of low-quality cells with low (≤500) UMI counts or more than 20% UMI counts from mitochondrial genes. Therefore cluster 5 was not pursued further in our analysis. The cells in each cluster were colored and labeled with numbered annotations assigned by the Bioconductor package SingleR using normal mouse cell-type marker genes and orthologs of human hepatoblastoma tumor signature genes (Song et al.,)^[261]65 as references. c Confusion matrices of cell clusters,