Integrated Genomic Analysis Identifies Driver Genes and Cisplatin-Resistant Progenitor Phenotype in Pediatric Liver Cancer

Pediatric liver cancers (PLC) are rare tumors and therefore have not been molecularly characterized in large series. Hepatoblastomas represent about 67% to 80% of all PLCs worldwide, generally developing before five years of age on nonfibrotic liver (1, 2). Although some rare syndromes such as familial adenomatous polyposis, Beckwith–Wiedemann syndrome, or Simpson-Golabi Behmel syndrome predispose to hepatoblastoma development, the etiology of hepatoblastomas is poorly understood because most of these tumors are sporadic. Hepatoblastomas are characterized by their histologic heterogeneity, with three main histology patterns—fetal, embryonal, and mesenchymal—that often coexist within a single tumor. A handful of genomics studies (3–7) have established hepatoblastomas as genetically simple tumors, with the smallest mutation burden among 24 pediatric cancer types (8). Beyond activating CTNNB1 (β-catenin) alterations found in most hepatoblastomas (70%–90%), only few recurrent driver mutations have been described, including NFE2L2 (5%–10%) and TERT promoter mutations (2%–5%). Pediatric hepatocellular carcinomas (HCC; incidence = 0.24–0.65 per 1,000,000) resemble adult HCC in their histology and frequently develop on fibrotic/cirrhotic liver as a consequence of hepatitis B virus infection or rare congenital disorders (1, 9). Previous genomic studies identified recurrent alterations in WNT signaling and telomerase pathways, but no oncogenic driver in pediatric HCC with underlying liver disease (10, 11). Fibrolamellar carcinomas (FLC), a rare subtype of HCC, develop in adolescents and young adults on healthy liver and have been characterized by a recurrent DNAJB1–PRKACA driver gene fusion (12). Finally, benign lesions such as hepatocellular adenomas (HCA) and focal nodular hyperplasias are usually related to congenital malformation associated with vascular abnormalities or metabolic genetic diseases, or occur after chemotherapy (2, 13).

Hepatoblastomas are usually treated with cisplatin-based neoadjuvant chemotherapy and subsequent surgical removal of the tumor, leading to >80% five-year survival (14, 15). However, some hepatoblastomas develop resistance to chemotherapy during the initial neoadjuvant chemotherapy or after tumor recurrence, and the molecular determinants of cisplatin resistance are yet to be discovered. In contrast to hepatoblastomas, pediatric HCCs respond poorly to chemotherapy, and as in adults, they have a poor prognosis if not completely removed by surgery.

We aimed to establish the detailed driver landscape of 126 pediatric liver tumors, then we analyzed the plasticity of hepatoblastoma tumors in relation to cisplatin resistance, and explored new therapeutic strategies to overcome this resistance.

We analyzed a cohort of 126 pediatric liver tumors comprising 104 hepatoblastomas that developed in 65 patients (87 primary tumors including 18 prechemotherapy and 17 relapses/metastases), 10 HCCs (9 patients), 7 FLCs and 5 HCAs sequenced by whole-genome sequencing (WGS; n = 65), whole-exome sequencing (WES; n = 57), RNA sequencing (RNA-seq; n = 120), and reduced representation bisulfite sequencing (RRBS; n = 92; Supplementary Fig. S1; Supplementary Table S1). Among the 65 patients with hepatoblastoma, 14 were older than 5 years of age at diagnosis. HCC developed in fibrotic or cirrhotic liver related to various constitutional liver diseases (tyrosinemia, mitochondrial cytopathy, progressive familial intrahepatic cholestasis, associated with germline mutations of FAH, NDUFA11, NDUFB9, TJP2, and ABCB11; Supplementary Table S2) whereas the other tumors arose on normal liver. Germline truncating mutations of BRCA1 and BRCA2 were also identified in two hepatoblastomas and one FLC. In WGS analyses, primary hepatoblastomas displayed a small number of somatic mutations (median = 886, 0.3 mutations/Mb), whereas hepatoblastoma metastases and relapses showed a massive mutation load (median = 12,824, 4.3 mutations/Mb) with a high proportion of doublet-base substitutions (8%; Fig. 1A). Overall, HCC displayed a higher mutation rate than hepatoblastomas (median 5,318, P = 0.002).

The WNT/β-catenin pathway was the most frequently altered oncogenic pathway in PLCs (84.5%), with different activating mechanisms across diagnoses. In hepatoblastomas, CTNNB1 alterations activating β-catenin were identified in 92% of the tumors (Supplementary Tables S3 and S4), with missense mutations exclusively observed in young patients (<4 years), whereas exon 3 in-frame deletions were observed later in life (Fig. 1B; Supplementary Fig. S2). The remaining five patients with hepatoblastoma without CTNNB1 mutation showed germline truncating mutations of APC (n = 3) or AXIN1 (n = 1) with somatic inactivation of the second allele in the tumor, and only one hepatoblastoma remained without an identified alteration in the pathway. CTNNB1 was also altered in 40% of HCAs. In HCC, no CTNNB1 mutations were detected, but biallelic inactivation of AXIN1 and/or AMER1 occurred in six of nine patients (67%).

The 11p15.5 imprinted locus, containing the IGF2 oncogene, was the second most frequently altered locus in hepatoblastomas and HCC (84% and 89%, respectively), mostly through copy-neutral loss of heterozygosity (cn-LOH; 51%–56%; Fig. 1C). Adding to LOH, we found epimutations of the imprinting control regions IC1 [gain of methylation (GOM), 22% of hepatoblastomas and 33% of HCC] and IC2 [loss of methylation (LOM), 5% of hepatoblastomas], and recurrent somatic mutations of CDKN1C in four patients with hepatoblastoma. In six patients with hepatoblastoma (10%), imbalanced B allele frequency (BAF) profiles revealed that the cn-LOH found in tumor cells was present as a mosaic in the adjacent normal liver (Fig. 2A and B). These young patients (median = 8.4 months, P = 0.045) were not diagnosed with Beckwith–Wiedemann syndrome (BWS). Yet, the cn-LOH of 11p15.5 was detected in a significant fraction of normal liver cells (6%–58%) without the other driver alterations identified in the corresponding tumor, indicating premalignant clonal expansions of hepatocytes. Oncogenic transformation involved the acquisition of CTNNB1 missense activating mutation in the six cases (Fig. 2B). Accordingly, IGF2 was highly overexpressed in the nontumor liver samples of two patients harboring a mosaic 11p15.5 alteration, whereas WNT/β-catenin target genes were overexpressed only in tumor cells (Fig. 2C). In patient #3559 (30% mosaic cn-LOH in the nontumor liver), RNAscope in situ hybridization revealed a massive overexpression of IGF2 in both the tumor and adjacent nontumor tissue, whereas glutamine synthetase and β-catenin immunostainings demonstrated an oncogenic activation of the WNT/β-catenin pathway specifically in tumor cells (Fig. 2D). These data support the idea that premalignant clonal expansions of “normal” hepatocytes with 11p15.5 cn-LOH overexpress IGF2 and can lead to hepatoblastoma formation after oncogenic CTNNB1 mutation. Finally, one patient with BWS (#3180) displayed IC2 mosaic epimutation in both blood cells and the liver (27% of cells). Overall, the mechanism of 11p15.5 alteration in PLC was related to age, with cn-LOH occurring more frequently in young patients and gain of methylation of IC1 mostly observed in older patients (Fig. 1B; Supplementary Fig. S2).

A paucity of structural variants (SV) was observed in primary hepatoblastomas (median = 6/tumor) except for six tumors with chromoplexy, which mostly developed in older patients (Fig. 1A; Supplementary Table S5). In HCC, we identified an unusual number of focal deletions (median = 47/tumor) of small size (median = 27 kb). This remarkable deletor phenotype led to recurrent cancer driver alterations, particularly on chromosome X with complete inactivation of AMER1 (44%), GPC3 (44%), RPS6KA3 (33%), SMARCA1 (22%), and BCORL1 (11%), but also on chromosome 16 with homozygous deletion or combined deletion and truncating mutation of AXIN1 (44%) and CREBBP (11%; Fig. 1B and C; Supplementary Fig. S3). Of note, GPC3, BCORL1, CREBBP, RPS6KA3, and AXIN1 were also inactivated in hepatoblastomas through truncating or damaging mutations, including one germline GPC3 mutation in a patient with Simpson-Golabi syndrome. Recurrent copy-number alterations were identified in two genes controlling p53 degradation: focal deletions at 4q35 pinpointed inactivation of the tumor suppressor IRF2 in 30% of hepatoblastomas and 56% of HCC, whereas focal amplifications at 1q32.1 led to a high expression of MDM4 in four patients with hepatoblastoma (Fig. 1C; Supplementary Fig. S3; Supplementary Table S6). Interestingly, gains of the entirety of 1q were also observed in 50% of hepatoblastomas and HCC (Fig. 1C). At the chromosome arm level, HCC and hepatoblastomas had a roughly similar profile of gains but HCC harbored more losses, including loss of the 13q arm (encompassing RB1) in 33% of cases (Supplementary Fig. S4). Finally, focal amplifications of CCND1/FGF19 were found in three patients with hepatoblastoma, leading to the overexpression of both genes. Other recurrent driver mutations in hepatoblastomas involved TERT promoter in nine patients (all older than 40 months, P = 1.5 × 10⁻⁵; Supplementary Fig. S2), NFE2L2 (four patients) and ARID1A (two patients). No recurrent gene fusion was identified in the cohort except the PRKACA–DNAJB1 fusion pathognomonic of FLC. We also identified recurrent mutations of ERBB4 in 2/7 FLCs (29%) and of HNF1A in 3 HCAs (60%), including one germline mutation. Overall, hepatoblastomas and HCC shared common pathways altered by diverse genes and mechanisms, while HCA and FLC had specific driver alterations (Fig. 1B).

Unsupervised transcriptomic classification of 100 hepatoblastoma samples from 64 patients revealed 4 robust molecular groups characterized by diverse differentiation states, cell proliferation, and immune infiltration levels (Fig. 3A). The differentiated “hepatocytic” hepatoblastoma group comprised 44 samples, 73% of which had only “fetal” histologic component in their corresponding mirror block (vs. 10%,P = 1.5 × 10⁻¹⁰). They strongly expressed transcription factors (TF) involved in hepatic differentiation (HNF1A and HNF4A; Fig. 3A). “Hepatocytic” hepatoblastomas were divided into two clusters: the “hepatocytic hot” subgroup was defined by a strong signature of polymorphic immune infiltration, including a mixture of T, B, and natural killer (NK) cells and macrophages, with 20/21 of these samples collected after chemotherapy; the “hepatocytic cold” subgroup was characterized by lower levels of immune infiltrates, and it was enriched in prechemotherapy samples (39%,P = 0.010; Fig. 3A; Supplementary Fig. S5A–S5C). The “liver progenitor” group comprised 44 samples enriched in highly proliferative, immune-cold tumors with embryonal compartments, 11p15.5 alterations, and MDM4 amplifications. They expressed TFs involved in hepatic differentiation but also TFs involved in self-renewal and pluripotency maintenance (MYCN, MIXL1, and LIN28B). The last group, “mesenchymal” hepatoblastomas, included all samples with a mesenchymal histology that were derived from mixed hepatoblastomas with both epithelial and mesenchymal components, and was found mostly in younger patients (Supplementary Fig. S2). These mesenchymal samples displayed a distinct differentiation program, with no expression of liver differentiation genes but a strong expression of mesenchymal stem cell TFs (TWIST1, TBX5, and MSX2), and variable levels of immune infiltration. We validated these four transcriptomic subgroups and their associations with histology and immune infiltration in an independent RNA-seq cohort of 34 hepatoblastoma samples (ref. 16; Supplementary Fig. S6A–S6C).

Strikingly, 14/24 patients (58%) with multiple synchronous and/or metachronous samples displayed transcriptomic group switches (Fig. 3B; Supplementary Fig. S7). Spatial transcriptomic heterogeneity was identified in the primary tumors of seven patients and matched histologic heterogeneity between fetal and embryonal parts of the tumor (Fig. 3B). Longitudinal transcriptomic group changes were also identified between pre- and post-chemotherapy samples (4 out of 6 patients) and between primary tumors and relapses/metastases (5 out of 8 patients). Phylogenetic trees revealed private driver alterations in some cases, but no recurrent gene associated with specific transcriptional group changes. In contrast, DNA methylation profiles were closely associated with differentiation states, and transcriptomic group changes were associated with DNA methylation reprogramming (Fig. 3C; Supplementary Fig. S8A). In particular, mesenchymal samples displayed coordinated hypermethylation of HNF4 and PPAR binding sites, both in our cohort (Supplementary Fig. S8B and C) and in an independent data set (ref. 16; Supplementary Fig. S6B).

The molecular plasticity of hepatoblastomas recapitulated the heterogeneity observed at the histologic level (Fig. 3B). A systematic histologic review at the sample (mirror block) and whole tumor levels revealed that 80% of primary hepatoblastomas displayed spatial heterogeneity with a mixture of embryonal, fetal, or mesenchymal areas. The “liver progenitor” molecular group was associated with increased intrasample histologic heterogeneity (P = 1.0 × 10⁻⁶; Fig. 3D), with frequent coexistence of fetal and embryonal cells (Fig. 3A and B). Primary tumors with alterations in the 11p15.5 locus also displayed more histologic heterogeneity (89% vs. 46%,P = 0.0021; Fig. 3E). Thus, the phenotypic plasticity of hepatoblastomas may relate to the multipotency of progenitor cells with 11p15.5 alteration.

Consistent with this plasticity, the proportion of molecular groups evolved across disease stages (P = 0.0014; Fig. 3B). The “hepatocytic hot” group was enriched in post- versus prechemotherapy samples (30% vs. 6%), whereas the “liver progenitor” group was enriched in metastases and relapses (64% vs. 40%). WNT/β-catenin alterations were always trunk in phylogenetic trees and already present in 100% of prechemotherapy primary samples. In contrast, 11p15 alterations occurred late in 8 out of the 20 affected phylogenetic trees and their frequency increased in later disease stages, which was also the case for IRF2 deletions (Fig. 3B).

To validate this phenotypic plasticity at the single-cell level, we analyzed matched nontumor liver, primary hepatoblastoma, and lung metastasis from one patient (#3981) by single-nucleus RNA-seq (snRNA-seq). Virtual copy-number profiles were consistent with the copy-number changes observed in bulk WES data, with gains of 1q, 2q, 5q, and 15q identified in all tumor cells, and an extra gain of Xq in all metastatic cells (Fig. 4A). In contrast to this genetic homogeneity, we identified clusters of tumor cells with distinct transcriptomic profiles, corresponding to the “mesenchymal,” “liver progenitor,” and “hepatocytic” signatures (Fig. 4B). The three populations were present in different proportions in each sample, with a larger “mesenchymal” contingent in the primary tumor, consistent with the molecular classification of bulk samples (“mesenchymal” for the primary, “liver progenitor” for the metastasis). These data demonstrate the existence of the three differentiation states at the single-cell level, and the plasticity of hepatoblastoma cells across these differentiation states. Indeed, both the last common ancestor (LCA) of the primary tumor and the LCA of the metastasis (with additional Xq gain) were able to generate the three cell types repeatedly during tumor evolution in this patient (Fig. 4C).

Thus, dynamic switches between three differentiation states operate in hepatoblastoma cells, resulting in high intrapatient heterogeneity in space and time. This phenotypic plasticity is more frequent in tumors with 11p15.5 alterations and involves transcriptional and epigenetic reprogramming of TF modules. Molecular phenotypes match well with histologic cell types and display drastically different immune infiltrates.

Mutational signature analysis of 65 pediatric liver cancer genomes identified four single-base substitution (SBS), two doublet-base substitution (DBS), and five insertion/deletion (indel) signatures, most of which matched signatures previously described in pan-cancer studies (refs. 17, 18; Fig. 5A; Supplementary Fig. S9A and S9B; Supplementary Table S7). The majority of signatures, including the clock-like signatures SBS1 and SBS5, were found ubiquitously in the tumors, while the signature SBS18, commonly found in neuroblastoma (17) and related to oxidative DNA damage (19), was identified in a subset of hepatoblastomas. Finally, we identified a massive load of mutations due to the signatures SBS35, DBS5, and ID3 in 20/66 primary hepatoblastomas resected after neoadjuvant chemotherapy (vs. 0/17 hepatoblastomas sampled before cisplatin treatment, P = 0.009). Consistently, SBS35 and DBS5 are known to reflect the diverse mutation types caused by cisplatin adducts on DNA (18, 20). In primary hepatoblastomas after neoadjuvant chemotherapy, cisplatin mutations were subclonal and almost exclusively found in the “liver progenitor” tumor subgroup (18/27 vs. 2/38; P = 1.2 × 10⁻⁷; Fig. 5B and C). Furthermore, in four tumors with spatial heterogeneity, cisplatin signature was restricted to the “liver progenitor” component even though the non-“progenitor” samples were exposed to the same chemotherapy regimen (Supplementary Fig. S7).

Interestingly, all 16 hepatoblastoma relapses and metastases that developed after chemotherapy displayed a deluge of clonal cisplatin-induced mutations (Fig. 5C), leading to high mutation loads (median = 4.3 mutations/Mb), comparable to those of adult tumors and 15 times higher than those of primary liver tumors of the same age (Supplementary Fig. S9C–S9E). We analyzed in detail eight patients with chemoresistant hepatoblastomas at different stages of the disease including primary tumors, relapses, and metastases operated on between 4 months and 12 years after chemotherapy. All relapses were clonally related to their matched primary tumors, even the lung metastasis of patient #3538 (germline APC), detected 12 years after initial surgery (Fig. 5D). Phylogenetic trees revealed that a sudden burst of mutations is acquired during cisplatin treatment. Every relapse/metastasis was derived from a single common ancestor cell that acquired between 5,000 and 13,000 mutations during cisplatin treatment, with various modes of metastatic seeding. In patients #3370 and #3981, SBS35 mutations were private to each sample, indicating that metastatic seeding involved independent resistant cells. In patient #3694, the two relapses displayed a mixture of shared and private SBS35 mutations, indicating successive rounds of cisplatin mutagenesis followed by clonal expansion of resistant cells. Finally, patients #3529 and #3949 displayed late branching of relapse/metastasis samples derived from a same cisplatin-resistant clone. The seeding abilities of these clones are clearly illustrated by patient #3529 whose primary tumor was treated with liver transplantation: the same cisplatin-resistant clone gave rise first to a metastasis in the spleen, followed by two subsequent metastases on the grafted liver (Fig. 5D).

Overall, these data suggest that “liver progenitor” cells are able to bypass cisplatin–DNA adducts and proliferate under chemotherapy while accumulating SBS35 mutations. Accordingly, we found an enrichment of genes related to cisplatin resistance (21) or DNA repair among genes overexpressed in “liver progenitor” hepatoblastomas (63 genes, P < 2.2 × 10⁻¹⁶). Among those, 20 genes were also upregulated in post- versus prechemotherapy “liver progenitor” hepatoblastomas (Supplementary Fig. S10A–S10D), including genes involved in inhibition of apoptosis (BIRC5, coding for Survivin) and DNA repair through homologous recombination (BRCA1, RAD54L, and EXO1), Fanconi anemia (FANCA, FANCB, FANCD2, and FANCI) or base-excision repair pathways (LIG1, LIG3, POLE, and POLE2). Of note, most of these genes were also overexpressed in fetal liver samples (13th–30th weeks of amenorrhea) versus postnatal liver.

The subclonal presence of the signature SBS35 in primary tumors after neoadjuvant chemotherapy is thus a marker of cisplatin-resistant cell proliferation and was associated with poor progression-free and overall survival (log-rank P = 0.012 and P = 0.032, respectively; Supplementary Fig. S11). Heavily mutated resistant cells later give rise to relapses and metastases. A median of 70 coding sequence mutations per relapse/metastasis occurred due to the extra mutation load attributed to cisplatin. Although some of these mutations affected known cancer genes including NF1, BRAF^L485W, KMT2C, KMT2D, BCORL1, and NOTCH1, no recurrent driver gene associated with tumor progression was identified.

To investigate new therapeutic options, we characterized a panel of nine pediatric liver cancer cell lines (PL-CCL; eight hepatoblastomas and one HCC; Fig. 6A; Supplementary Table S8). Our panel of cell lines, which are mainly derived from older patients and of “liver progenitor” type with multiple alterations in cancer driver genes, reproduced the major driver events identified in PLC. In particular, three cell lines with a matched primary or relapsed hepatoblastoma displayed genomic alterations globally similar to the original tumors (Supplementary Table S9). All hepatoblastoma cell lines carried CTNNB1 alteration whereas Hep3B, derived from an HCC, carried TP53, AXIN1, and RPS6KA3 alterations with frequent homozygous deletions in agreement with the deletor phenotype that we observed in pediatric HCC. Furthermore, all cell lines exhibited 11p15.5 locus alterations through cn-LOH (5/9), GOM IC1 (3/9), or CDKN1C mutation (1/9).

We tested our panel of cell lines with conventional therapeutic regimens (cisplatin, carboplatin, and doxorubicin treatments; Supplementary Table S10). In the nine cell lines analyzed, resistance to cisplatin and carboplatin correlated with a high expression of “liver progenitor” markers (Fig. 6B; Supplementary Fig. S12A). Additionally, we treated five cell lines with long-term exposure to low doses of cisplatin (0.5 μmol/L, 4 weeks; Supplementary Fig. S12B). Only B6-2 was completely sensitive after two weeks. In the four remaining cell lines, WES of resistant clones revealed an accumulation of the SBS35 cisplatin mutational signature. These results demonstrate that the SBS35 signature is directly induced by cisplatin exposure in resistant cells (Fig. 6C). Interestingly, doxorubicin, which is commonly used in resistant or high-risk patients with hepatoblastoma, exhibited a significant antitumor effect in all cell lines.

Next, we tested treatments targeting specific genomic alterations and their consequences. Three out of five tested cell lines were sensitive to anti-IGF2 antibodies, and sensitivity to this drug correlated with the level of IGF2 expression (P = 0.03, R = −0.91, Pearson correlation; Fig. 6D). We also tested the monoclonal therapeutic antibody xentuzumab, and we obtained between 14% and 47% growth inhibition (Supplementary Fig. S12C), consistent with a previous report in mice xenografts (22). Surprisingly, sensitivity to xentuzumab was not correlated with IGF2 expression, possibly because it neutralizes both IGF1 and IGF2. We then tested our cell lines with the MEK1/2 inhibitor trametinib (Fig. 6E). Remarkably, five of nine cell lines were sensitive to trametinib and harbored alterations within the MAP kinase pathway through RPS6KA3 homozygous deletions, NRAS-activating mutations, or FGF19 amplifications as previously shown in adult HCC cell lines (23). Of note, one of our patients with a metastatic hepatoblastoma (#3538) displayed a hotspotMAPK1^E322K trunk mutation and acquired a BRAF^L485W mutation in the metastasis that has previously been associated with response to an ERK inhibitor (24).

To identify new targets in cisplatin-resistant “progenitor” cells, we performed a differential expression analysis between the “liver progenitor” and other hepatoblastoma subgroups. Of the genes overexpressed in “liver progenitor” samples, we selected three key targetable genes: PLK1, BIRC5 (Survivin), and CHEK1, which are involved in cell-cycle, apoptosis, and DNA damage response, respectively (Fig. 6F and G). In patients with hepatoblastoma, high expression of PLK1, BIRC5, and CHEK1 was associated with (i) poor response to chemotherapy as assessed with AFP reduction in the sera (P = 4.1 × 10⁻⁴, 2.0 × 10⁻⁴, and 4.7 × 10⁻⁴, respectively), (ii) embryonal histology (P = 1.2 × 10⁻⁵, 5.1 × 10⁻⁶, and 7.1 × 10⁻⁶), (iii) the cisplatin-resistance signature SBS35 (P = 8.8 × 10⁻¹⁰, 1.3 × 10⁻⁹, 7.6 × 10⁻⁹), and (iv) a high expression of proliferation genes (P ≤ 2.2 × 10⁻¹⁶; Supplementary Fig. S13A–S13D). We thus tested inhibitors of PLK1 (BI-2536), Survivin (YM155), and CHEK1 (AZD-7762) in our collection of nine cell lines. Interestingly, the three drugs were more efficient than both cisplatin and carboplatin in all the cell lines except B6-2, which has the least “liver progenitor”-like phenotype. All cell lines were sensitive to BI-2536, including those that are resistant to cisplatin (Fig. 6H). Notably, the cell lines most resistant to cisplatin were also the most sensitive to YM155 and AZD-7762 (Supplementary Fig. S14A–S14C). Finally, we tested the most effective drug, BI-2536, in vivo in 24 nude mice xenografted with HepG2 cells randomized in four groups (cisplatin, cisplatin + doxorubicin, BI-2536, and control). HepG2 xenografts were resistant to cisplatin treatment, both in isolation and in combination with doxorubicin. In contrast, tumors in mice receiving BI-2536 responded to treatment (P = 0.003) without significant toxicity, unlike with cisplatin (Fig. 6I and J).

In conclusion, the PLK1 inhibitor BI-2536 showed efficacy in both in vitro and in vivo models and appears to be a good candidate drug for the treatment of PLC. Trametinib and anti-IGF2 antibodies show promise for the treatment of resistant hepatoblastomas with CCND1/FGF19 amplification or high IGF2 expression, respectively.

Altogether, our findings suggest a major role of cell plasticity in the spatiotemporal evolution of hepatoblastomas and their resistance to cisplatin-based chemotherapy. In this study, genomic analysis of multiple parts of primary tumors and of longitudinal samples at various time points allowed us to profile the clonal evolution of hepatoblastoma along treatment. Cisplatin induced hundreds of mutations in progenitor tumor compartments, a phenomenon reproduced in cell lines, and targeting specific oncogenic processes with different inhibitors could be used to target resistant cells.

Primary liver cancers showed a low mutation burden, as seen in other pediatric cancers (8). In PLC, β-catenin activation, mainly caused by CTNNB1 mutations, acted as the major early and common mechanism of liver tumorigenesis. However, we identified an earliest recurrent premalignant alteration with mosaic cn-LOH of 11p15 in the liver of very young patients with hepatoblastoma (10%). This finding is in line with similar premalignant clonal expansions recently identified as precursors of Wilms tumors, associated with hypermethylation of IC1 at the 11p15 locus (25). Our data reinforce the role of 11p15 alterations in premalignant stages that occur during pregnancy or early in life. Systematic searches of 11p15 locus alterations in the liver and other organs of patients with pediatric cancer could reveal additional mosaic cases.

The present molecular transcriptomic classification revealed four robust hepatoblastoma subgroups defined by both their differentiation states and immune infiltration levels. We verified that our clustering was not driven by patients with multiple samples (Supplementary Fig. S15), and we validated our classification in an independent RNA-seq cohort of 34 hepatoblastomas (Supplementary Fig. S6). Our transcriptomic classification overlaps the previously published subtypes C1/C2 (26), C1/C2A/C2B (27), and MRS1/MRS2/MRS3a/MRS3b (16), with a good match between “hepatocytic” and C1, “liver progenitor” and C2A, and “mesenchymal” and C2B subgroups (Fig. 3A; Supplementary Fig. S16). Our “liver progenitor” subgroup also overlaps with the very high-risk MRS-3b subgroup identified by Carrillo-Reixach and colleagues (ref. 16; P = 2.0 × 10⁻⁶). The most original findings of our classification are (i) the definition of clear transcriptomic signature of “mesenchymal” histologic cell type validated by careful histologic review of mirror blocks, and (ii) the identification of two subgroups of “hepatocytic” hepatoblastomas showing varying levels of immune infiltration, which is potentially enhanced by cisplatin-based chemotherapy. Hepatoblastoma subgroups are not defined by specific driver events, as opposed to medulloblastoma subgroups (28), for example. In turn, each state of hepatoblastoma differentiation is defined by specific TFs and a characteristic DNA methylation landscape. Our methylation-based classification revealed two main hepatoblastoma clusters reminiscent of the previously described Epi-CA/CB groups (16): a normal-like cluster matching the Epi-CA group comprising most “hepatocytic” samples, and a differentially methylated cluster matching the Epi-CB group, which we further divided into three subgroups associated with “liver progenitor,” “mesenchymal” phenotype, and older age. Importantly, hepatoblastoma cells are able to switch between differentiation states along the course of the disease, leading to striking spatial and temporal heterogeneity. Our single-cell data support the notion that hepatoblastoma cells have the inherent ability to switch across three differentiation states. Yet, clinical and molecular features associated with transcriptomic subgroups suggest that the cell of origin, patient age, chemotherapy, and driver alterations can favor the predominance of a given subtype in each tumor (Supplementary Table S11; Supplementary Fig. S17). This molecular plasticity, driven by specific TF modules and epigenetic landscapes, is reminiscent of neuroblastoma cells that can transdifferentiate between committed adrenergic cells and undifferentiated mesenchymal cells (29, 30), the latter being chemoresistant and enriched in relapse tumors. Similarly, we showed that hepatoblastoma “progenitor” cells proliferate under neoadjuvant chemotherapy and are enriched in relapses and metastases.

Moreover, PLC development follows evolutionary trajectories that vary according to age, involving similar pathways but diverse mechanisms of activation. In contrast to hepatoblastomas, pediatric HCCs occur in older children, on cirrhotic liver due to various causes. Like hepatoblastomas, HCC displayed frequent 11p15.5 alterations, but also a characteristic chromosome deletor phenotype leading to WNT/β-catenin activation through AXIN1 or AMER1 deletions, as well as frequent GPC3, RPS6KA3, and SMARCA1 alterations. Because a deletor phenotype is seen in only 3% of adult HCC (31), this oncogenic mechanism could be associated with young age. Interestingly, many HCC driver genes are located on the X chromosome (AMER1, RPS6KA3, GPC3, and SMARCA1), which may partly explain the enrichment of pediatric HCC in male patients described previously (9) and retrieved in the present cohort (7/9 males, 78%). Pediatric HCCs are commonly resistant to cisplatin, but their mechanism of resistance likely differs from that involved in hepatoblastomas because they do not exhibit similar progenitor and stem cell features.

“Progenitor” cells are able to proliferate under treatment, bypass cisplatin-induced adduct formation, acquire a huge number of cisplatin-induced mutations (median 11,142) with a characteristic SBS35 signature, and give rise to highly mutated relapse tumors. Various DNA polymerases and DNA repair genes involved in double-strand break repair, homologous recombination, and Fanconi pathway are already overexpressed both in “liver progenitor” hepatoblastomas before chemotherapy and in the fetal liver, showing that these mechanisms are constitutively active in liver progenitor cells. Relatedly, a previous work had shown that Fanconi anemia inhibitors were able to block the growth of hepatoblastoma cells in vitro and in vivo (27). Pich and colleagues estimated that ∼15 of every 1,000 cisplatin-induced mutations affect the sequence of coding genes, of which ∼0.7 are expected to affect the sequence of known cancer genes (32). Thus, an average hepatoblastoma relapse or metastasis is at risk of acquiring 167 coding mutations, including 7.8 in cancer genes due to cisplatin treatment. Although most cancer gene mutations are unlikely to play a functional role in hepatoblastoma cells, this increased mutational burden provides an opportunity for resistant cells to acquire additional oncogenic capabilities. These data highlight the necessity to identify alternative treatments for cisplatin-resistant hepatoblastomas, in which more aggressive chemotherapy regimens may just promote the selection of aggressive progenitor cells with a massive extra mutational burden.

Our findings revealed new potential therapies to combat hepatoblastoma resistant to cisplatin. Targeting the different mechanisms involved in cisplatin resistance (21) such as apoptosis (YM155 targeting Survivin), DNA repair (AZD-7762 targeting CHEK1), or cell-cycle control (BI-2536 as a PLK1 inhibitor) is very promising because an efficient antitumor effect was observed in PLC cell lines with a progenitor molecular signature. In the same line, other drugs targeting the proteasome (bortezomib) could also be efficient in resistant hepatoblastomas (27). Immunotherapy could be another appealing type of treatment for hepatoblastomas because we showed that cisplatin treatment can induce an intratumor polymorphic immune response leading to “hepatocytic hot” hepatoblastomas, consistent with our previous observation in hepatoblastomas with germline APC mutations (33). Reinforcement of an efficient intratumor immune response could be improved by the use of oxaliplatin to induce an immune cell death (34) or by treatments with immunomodulators. However, “liver progenitor” hepatoblastomas are immune- cold despite showing multiple cisplatin-induced mutations. In these cases, penetrating the progenitor compartments with immune cells will be challenging. Ongoing and planned clinical trials aim to test immunomodulators or targeted therapies in high-risk hepatoblastoma; however, these trials are limited by the small number of patients.

In conclusion, PLC showed various mechanisms of tumorigenesis related to age and cell of origin. Hepatoblastomas demonstrated a striking spatial and longitudinal phenotypic plasticity related to the progenitor compartment associated with cisplatin resistance and the mutational signature of DNA adduct bypass. And, finally, the identification of drugs targeting “progenitor” cells opens new avenues to treat children at high risk of resistance.

A series of 126 liver tumor samples and their nontumor counterparts were collected from 86 patients surgically treated in various French hospitals. The study was approved by the local Ethics Committee (CCPRB Paris Saint-Louis). Written informed consent was obtained in accordance with French legislation. All samples were immediately frozen in liquid nitrogen and stored at −80°C. Tumors included in this study comprised 87 primary hepatoblastoma (18 prior chemotherapy and 69 after chemotherapy) and 17 hepatoblastoma relapses and metastases (from 65 patients with hepatoblastoma including 14 older than five years old), five HCAs (five patients), seven FLCs (seven patients), and ten HCCs (nine patients). HCCs were developed in various etiologic contexts including tyrosinemia (two patients), mitochondrial cytopathy (three patients), and progressive familial intrahepatic cholestasis (three patients). Samples were analyzed by a combination of WGS (n = 65), WES (n = 57), RNA-seq (n = 120) and RRBS (n = 92). A summary of the cohort is provided in Supplementary Fig. S1, and detailed clinical characteristics of each sample are provided in Supplementary Table S1. Also, seven human fetal liver samples between the 13th and 30th weeks of amenorrhea were analyzed by RNA-seq.

All tumors were reviewed by three expert pathologists specializing in pediatric liver tumors. For hepatoblastoma, fractions of histologic components (fetal, embryonal, mesenchymal, cholangioblastic, and small cell undifferentiated) were estimated according to the consensus classification (14) for the whole tumor as well as for mirror blocks corresponding to frozen samples when available. Thus, spatial heterogeneity was defined by the coexistence of at least two histologic components, either at the whole tumor level or at the intrasample level (mirror block).

We extracted DNA using the Maxwell DNA extraction kit (Promega) or the AllPrep DNA/RNA/miRNA Universal Kit (Qiagen). Sixty-five tumors and matched nontumor samples were sequenced at the Centre National de Génotypage (Evry, France) on an Illumina HiSeqX5 as paired-end 150 bp reads, with an average depth of 90× for tumors and 30× for nontumor liver samples. Sequence reads were aligned to the hg19 version of the human genome using BWA (35). We used Picard tools (http://broadinstitute.github.io/picard/) to remove PCR duplicates and GATK (36) for local indel realignment and base quality recalibration, as recommended in GATK best practices (37).

Whole-exome data from 57 tumors and matched nontumor samples were analyzed in this study, as well as nine pediatric liver tumor cell lines. Sequence capture, enrichment, and elution of genomic DNA were performed by IntegraGen (Evry). Agilent in-solution enrichment was used with the manufacturer's biotinylated oligonucleotide probe library SureSelect Clinical Research Exome V2 or Twist Bioscience Human Core Exome Enrichment System, according to the manufacturer's instructions. The eluted enriched DNA samples were sequenced on an Illumina HiSeq4000 as paired-end 75 bp reads (n = 35) or Illumina NovaSeq as paired-end 100 bp reads (n = 31), with an average depth of 100× for tumors and 65× for nontumor liver samples. We used BWA to align reads on the hg38 version of the human genome and sambamba to remove duplicate reads.

We used MuTect2 to call somatic mutations from WES and WGS data by comparing each tumor sample with its matched nontumor counterpart. The three cell lines derived from patient tumors within our cohort were compared with the matched nontumor samples, while the six other cell lines were compared with a panel of normal samples. We excluded mutations belonging to the ENCODE Data Analysis Consortium blacklisted regions (http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeMapability/wgEncodeDacMapabilityConsensusExcludable.bed.gz) and regions covered by < 6 reads in the tumor or normal sample. We then selected only single-nucleotide variants (SNV) with a MuTect2 flag among “PASS,” “clustered_events,” “t_lod_fstar,” “alt_allele_in_normal,” or “homologous_mapping_event” and small insertions and deletions (indels) with a MuTect2 flag among “PASS,” “clustered_events,” or “str_contraction.” To improve specificity in the calling of mutations with low variant allele frequency (VAF), we quantified the number of high-quality variant reads in the tumor (mapping quality ≥ 20, base quality ≥ 20) and the number of variant reads in the nontumor sample with no quality threshold using bamreadcount (https://github.com/genome/bam-readcount). Only variants matching the following criteria were finally retained: VAF ≥ 2% in the tumor with ≥ 3 variant reads, VAF ≤ 5% in the nontumor samples with ≤ 2 variant reads, and a VAF ratio ≥ 5 between the tumor and nontumor sample.

Germline variant calling was performed independently for WGS and WES data according to GATK (version 4.0) best practices recommendations (38). For each data set, we performed SNV and indel discovery using HaplotypeCaller and joint genotyping across all non-tumor samples simultaneously using GenomicsDBImport followed by GenotypeGVCFs. We used hard filtering for WES data with different parameters for SNVs (ExcessHet > 54.69, FS > 10.0, MQ < 50.0, MQRankSum < −5.0, QD < 5.0, QUAL < 50.0, ReadPosRankSum < −5.0, SOR > 2.0, DP < 8.0, GQ < 20.0) and indels (ExcessHet > 54.69, FS > 200.0, QD < 2.0, QUAL < 30.0, ReadPosRankSum < −20.0, DP < 8.0, GQ < 20.0), and we used variant quality score recalibration for WGS (truth sensitivity level: SNVs = 99.6%, indels = 95.0%). Resulting high-quality variants were annotated using the Variant Effect Predictor toolset (39). We selected variants with an allele frequency lower than 0.01 or lacking in the gnomAD database (40). All candidate pathogenic variants were manually verified using the Integrative Genomics Viewer (41).

We used MANTA (42) software to identify somatic structural rearrangements in WGS data. To keep only the most reliable events, we selected only rearrangements supported by ≥ 15 reads and with a variant allele fraction ≥ 10% in the tumor, and not more than 1 variant read in the nontumor counterpart. We used the cgpBattenberg (43) algorithm to reconstruct absolute copy-number profiles from WGS data and estimate tumor purity. We used the Genome Alteration Print method (44) to infer absolute copy-number profiles and tumor purity from WES data, and the circular binary segmentation algorithm implemented in the Bioconductor package DNAcopy (45) to identify focal homozygous deletions and high-level amplifications, as previously described (46).

We used two different approaches to identify driver genes in PLCs, using both WES and WGS data. First, the MutSigCV tool (47) was used to identify genes with significantly recurrent mutations while taking into account gene size and genomic covariates. High-level amplifications and homozygous deletions were included in this analysis as additional mutation categories. Second, we used Oncodrive (48) to identify genes with a significant enrichment of mutations with functional impact. Finally, we defined as putative drivers genes with damaging alterations in ≥2 patients and either (i) a P value < 0.05 with both MutSigCV and Oncodrive tools or (ii) identified as driver genes in previous pediatric pan-cancer studies (25, 49). After removing three genes not expressed in normal liver nor in tumors (95th FPKM quantile < 0.1), we obtained 22 candidate driver genes (Supplementary Table S2). Only genes with damaging alterations in ≥3 patients were represented in Fig. 1B. The frequency of driver alterations per tumor type was estimated by patient and not by sample, in order to avoid biases due to patients with multiple samples. To this end, a patient was considered as altered for a gene as long as one of its samples harbored an alteration in this gene.

For each tumor with a cn-LOH of 11p15.5 locus, we searched for the presence of the same cn-LOH in matched nontumor liver tissue. To that aim, we calculated the BAF of common single-nucleotide polymorphisms (SNP), obtained from WGS or WES data, in the tumor (BAF_T) and in the nontumor (BAF_NT) sample. We considered a cn-LOH to be present in the nontumor liver if there was a significant overlap between SNPs with the B allele retained in the tumor (BAF_T > 0.5) and those with a BAF_NT > 0.5 (binomial test). We then estimated the proportion of nontumor liver cells carrying the cn-LOH as 2*RAF_NT-1 with RAF_NT the median BAF of retained alleles (those with a BAF_T > 0.5) in the nontumor sample.

We used Palimpsest (50) to extract signatures of SBS, DBS, and indels in WGS data of 65 PLCs, and to compare them with known signatures from the COSMIC database (v3). De novo analysis revealed four SBS signatures corresponding to the COSMIC SBS1, 5, 18, and 35 (cosine similarity scores > 0.8), two DBS signatures corresponding to the COSMIC DBS5 and a new signature (DBSnew), and five indel signatures corresponding to the COSMIC ID1, 2, 3, 5, and 8. To ensure comparability with other data sets, signature exposures in each tumor were recalculated using the COSMIC versions of signatures, except for signature SBS35 (we kept our own version as it was based on a larger number of mutations than the COSMIC version) and DBSnew (absent in COSMIC). SBS signature exposures were also calculated in WES data. Finally, we estimated the cancer cell fraction (CCF), i.e., the proportion of tumor cells carrying each mutation (see the “Clonality” section below). Mutations with an upper bound of their CCF confidence interval < 0.95 were considered subclonal. We then quantified SBS signature exposures among clonal and subclonal mutations separately.

We used Palimpsest to estimate the CCF of each mutation from its variant allele fraction, taking into account tumor purity and absolute copy number, as previously described (31). For 23 patients with multiple tumor samples, we used a Bayesian Dirichlet process in multiple dimensions (51) to identify clusters of mutations with a similar CCF distribution across all samples, hence belonging to the same branch of the phylogenetic tree. We then manually reconstructed the phylogenetic tree of each patient by organizing branches to fit the observed clonal composition of each sample. Damaging mutations affecting candidate driver genes were annotated on the trees, and we quantified the contribution of each mutational process on each branch with Palimpsest.

We performed RNA-seq for 120 tumors, 4 nontumor liver samples, 7 fetal liver samples, and 9 pediatric liver tumor cell lines. RNA samples were enriched for polyadenylated RNA from 1 μg of total RNA, and the enriched samples were used to generate sequencing libraries with the Illumina TruSeq Stranded mRNA Kit or the NEBNext Ultra II Directional RNA Library Prep Kit and associated protocol as provided by the manufacturer. Libraries were sequenced by IntegraGen (Evry) on an Illumina HiSeq 4000 as paired-end 75 bp reads or Illumina NovaSeq as paired-end 100 bp reads. Full Fastq files were aligned to the reference human genome hg38 using TopHat2 (52). We removed reads mapping to multiple locations, and we used HTSeq (53) to obtain the number of reads associated with each gene in the Gencode database. We used DESeq2 (54) to import raw read counts into R statistical software and apply variance stabilizing transformation to the raw count matrix. FPKM scores (number of fragments per kilobase of exon model and millions of mapped reads) were calculated by normalizing the count matrix for the library size and the coding length of each gene.

Fusions detected by TopHat2 (ref. 52; –fusion-search –fusion-min-dist 2000 –fusion-anchor-length 13 –fusion-ignore-chromosomes chrM) were filtered using the TopHatFusion-post algorithm and validated using FusionInspector (ref. 55; https://github.com/FusionInspector). We kept only fusions validated by BLAST, with at least 10 split-reads or read pairs spanning the fusion event, and an FFPM (fusion fragments per million reads) ≥ 0.1. We removed fusions identified recurrently in a cohort of 36 normal liver samples, involving genes with inconsistent orientations, noncoding genes, and putative read-through transcripts.

Large in-frame CTNNB1 deletions were screened in WGS, WES, and RNA-seq data using dedicated approaches. We used MANTA software (42) to identify structural variations spanning CTNNB1 region in WGS and WES (with the –exome option) data, with the following parameters to increase sensitivity: minEdgeObservations = 2, minScoredVariantSize = 10. For RNA-seq, we analyzed the junctions.bed output files from TopHat2 (52) to identify abnormal junctions involving CTNNB1 exons 3 or 4. We also screened for paired reads with one mate in exon 2 and the other in exon 3/4, and an abnormally long insert size, greater than the mean insert size of the library + 1 standard deviation.

We restricted the transcriptomic classification to hepatoblastoma samples, excluding other diagnoses (HCC, FLC, and HCA) with a limited number of samples, in order to identify robust subgroups. Hierarchical clustering was performed on 100 hepatoblastomas and 4 nontumor liver samples, based on the vst-normalized expression of the 3,000 autosomal genes with sufficient expression (95th FPKM percentile ≥ 0.1) and the highest standard deviation. The hclust R function was used with Euclidean distance and Ward.D2 linkage method, after median-centering the data. To verify that hepatoblastoma subgroups were not driven by patients with multiple samples, we reproduced the clustering with only one sample by patient and we obtained similar results (Supplementary Fig. S15). For external validation of hepatoblastoma clusters, we downloaded an RNA-seq data set from Carrillo-Reixach and colleagues (16) comprising 34 hepatoblastomas (Gene Expression Omnibus accession number GSE132219) and performed a hierarchical clustering based on the expression of the 19 markers of hepatic differentiation, progenitor, and mesenchymal subgroups represented in Fig. 3, with Cosine distance and Ward.D2 linkage method (Supplementary Fig. S6). We used the Bioconductor limma (56) package to identify differentially expressed genes between hepatoblastoma subgroups. We applied a q-value threshold of ≤ 0.05 to define differentially expressed genes.

We used previously established molecular signatures from the MSigDB database to quantify the level of hepatic differentiation (“Hsiao liver-specific genes”), cell proliferation (“Hallmark REACTOME cell cycle”), and inflammation (“Hallmark inflammatory response”) in each sample. We used MCPcounter (57) to estimate the infiltration by diverse immune cell types and classify hepatoblastoma samples accordingly. We also classified each hepatoblastoma sample according to the previously described C1/C2 (27), C1/C2A/C2B (28), and RMS-1/2/3a/3b (16) classifications. For C1/C2 (respectively C1/C2A/C2B) classification, we generated a consensus clustering (1,000 resampling iterations of hierarchical clustering with Euclidean distance and Ward.D2 linkage) of our cohort in two (respectively three) groups based on the 16 marker genes defined by Cairo and colleagues (ref. 27; respectively the four markers VIM, HSD17B6, TOP2A, and ITGA6 defined by Hooks and colleagues; ref. 28), and we assigned each consensus cluster to the relevant subgroup based on the expression of their respective markers. For the RMS classification, we first performed a hierarchical clustering based on the expression of 68 genes of the 14q32 imprinted locus to define two groups (moderate vs. strong 14q32-gene signature). We verified that our two main DNA methylation-based clusters were consistent with Carrillo-Reixach's Epi-CA/CB groups (16). Finally, we combined the 14q32 signature, Epi-CA/CB methylation groups, and C2A transcriptomic subgroup to define the RMS classification as in Carrillo-Reixach and colleagues (16). A signature of the “Liver progenitor” subgroup of hepatoblastomas comprising the seven representative genes indicated in Fig. 3A was used to quantify this phenotype in cell lines.

Single nuclei were isolated from four samples for snRNA-seq as previously described (58), using EZ Lysis buffer workflow with slight modifications. Briefly, tissue samples were thawed in PBS and cut into pieces < 0.5 cm. Approximately 35 mg of tissue was poured in a glass Dounce tissue grinder (Sigma, cat. no. D8938) and homogenized 25 times with pestle A and 25 times with pestle B in 1.5 mL of ice-cold nuclei EZ lysis buffer. Samples were then incubated on ice for 5 minutes with an additional 3 mL of cold EZ lysis buffer. Nuclei were centrifuged at 500 × g for 5 minutes at 4°C, washed with 5 mL ice-cold EZ lysis buffer, and incubated on ice for 5 minutes. After centrifugation, the nucleus pellet was washed with 5 mL of Nuclei Wash buffer containing 1 × PBS, 0.1%, nonacetylated BSA (Thermo AM2618), and 200 units/mL RNase inhibitor (NEB M0307L). Isolated nuclei were resuspended in 2 mL of Nuclei Suspension Buffer containing 1 × PBS, 1% nonacetylated BSA (Thermo AM2618), and 200 units/mL RNase inhibitor (NEB M0307L), filtered through a 70-μm and then a 30-μm MACS SmartStrainer (Miltenyibiotec 130-098-462 and 130-098-458), and counted under microscope using C-chip disposable hemocytometer. A final concentration of 1,000 nuclei per μL was used for loading on a 10× channel.

snRNA-seq was performed by Integragen SA (Evry) on matched nontumor liver, primary hepatoblastoma, and lung metastasis samples from one patient (#3981), following the Chromium Next GEL Single Cell 3′ V3.1 protocol. In short, about 8,800 single nuclei were loaded into each channel of a Chromium single-cell 3′ chip. Single nuclei were partitioned into droplets with gel beads in the Chromium Controller. After emulsions were formed, barcoded reverse transcription of RNA took place, followed by cDNA amplification, fragmentation, and adapter and sample index ligation, according to the manufacturer's recommendations. Libraries from the 10X channels were pooled together and sequenced as paired-end 100 bp reads on an Illumina NovaSeq. We used 10X Genomics Cell Ranger 5.0 (59) to align snRNA-seq reads to the human genome (GrCh38/hg38) and generate UMI counts for each sample, including intronic reads. We obtained, respectively, 9,142, 5,859, and 2,875 nuclei, with a median of 2,042, 2,334, and 3,446 genes per nuclei from samples #3981N (nontumor liver), #3984T (primary hepatoblastoma), and #3988T (lung metastasis). We then filtered the feature-barcode matrix to retain only good-quality nuclei and reliable genes. We removed nuclei with <1,000 read counts, <500 genes detected, or >5% of UMI counts in mitochondrial genes, and we removed genes detected in <3 nuclei as well as ERCC and mitochondrial genes. After QC we kept, respectively, 9,140, 5,685, and 2,826 nuclei, with a total of 16,470, 18,384, and 18,294 genes for samples #3981N, #3984T, and #3988T. All secondary analyses were performed using Seurat v3 (60). We normalized each data set using the SCTransform function with default parameters, performed principal component analysis on the 3,000 most variable genes and ran Louvain graph-based clustering on the 30 principal components with a resolution of 0.5. We used a Uniform Manifold Approximation and Projection (UMAP) with default settings to visualize the results. We used inferCNV v1.6 (Tickle T and colleagues, available from https://github.com/broadinstitute/inferCNV) to reconstruct virtual copy-number profiles, using healthy hepatocytes from the nontumor liver sample as reference, and genes with an average read count >0.1 in reference nuclei. We performed a clustering of copy-number profiles in each sample separately. We identified tumor cells using both inferCNV clusters and Seurat classification, and we redid the whole analysis (SCTransform normalization, dimensionality reduction, clustering and UMAP visualization) restricted to tumor cells. To characterize tumor cell clusters, we selected marker genes of each differentiation state (“hepatocytic,” “liver progenitor,” and “mesenchymal”). To that end, we selected the top 70 genes overexpressed in the corresponding group in bulk RNA-seq data, restricted to genes identified in snRNA-seq data. We obtained ∼50 marker genes for each subgroup, and we computed for each nucleus the mean log-normalized expression (NormalizeData pipeline) over each set of markers. We then assigned a status to each cluster based on the expression of these marker genes.

RRBS was performed in two distinct projects (PJ17, n = 82 and PJ20, n = 23) by Integragen SA (Evry) as described by Gu and colleagues (61), with the Diagenode Premium RRBS kit. In brief, 100 ng of qualified genomic DNA were digested with MspI. After end-repair, A-tailing and ligation to methylated and indexed adapters, the size-selected library fragments were subjected to bisulfite conversion and PCR-amplified. Samples of the PJ17 project (respectively PJ20) were then sequenced on an Illumina HiSeq4000 (respectively NovaSeq) sequencer as paired-end 75 bp (respectively 100 bp) reads. Image analysis and base calling were performed using Illumina Real-Time Analysis with default parameters. Reads were aligned to the hg38 version of the human genome using BS_Seeker2 (62). Sorted bam files were converted into CGmap files using CGmaptools (63), and the methylation level of each CpG site was defined as the ratio between the number of effective cytosines after bisulfite conversion (= methylated cytosines) and the total number of cytosines and thymines after bisulfite conversion (= methylated + unmethylated cytosines). On average, approximately 7 million CpG sites were covered in each sample after discarding CpG sites located in ENCODE-blacklisted genomic regions (wgEncodeDacMapabilityConsensusExcludable tract from the UCSC genome browser). We next integrated methylation levels across 100-bp-long genomic regions (tiles) using the tileMethylCounts function from the Methylkit package (64). Tile coverage was defined as the coverage sum of all CpGs inside the tile, and methylation was defined as the ratio between the total number of methylated CpGs and tile coverage. On average, approximately 1.3 million tiles were covered in each sample. Due to the different read lengths in PJ17 and PJ20, some tiles displayed heterogeneous coverage leading to systematic biases between the two projects. We thus compared the nontumor liver samples from PJ20 (n = 4) and PJ17 (n = 9), and we removed 58,179 tiles (∼3%) with a methylation difference ≥ 0.05 or ≤ −0.05 between all PJ17 and all PJ20 nontumor samples.

To determine the methylation status of 11p15.5 locus, we computed for each sample the mean methylation over imprinted regions IC1 (chr11:1998745–2003509, hg38) and IC2 (chr11:2697587–2700983, hg38). We then used a K-means clustering to identify samples with GOM of IC1 and/or LOMof IC2.

We generated a DNA methylation–based classification of 84 hepatoblastomas and 13 nontumor liver samples. To do so, we selected the 5,000 tiles with coverage >50× in all samples and the highest standard deviation, and we used the hclust R function to perform a hierarchical clustering with Pearson distance and Ward.D linkage method. We used Independent Component analysis as previously described (65) to characterize independent sources of DNA methylation changes in hepatoblastomas, based on ∼212,000 tiles with coverage >10× in all samples and standard deviation >0.035. We identified five methylation components (MC) including components related to the “liver progenitor” (MC1) and “mesenchymal” (MC2) subgroups, and to patient age (MC3). We explored various (epi)genomic features associated with differentially methylated tiles and MCs, including gene features (promoter/gene body), CpG island features (island/shore/shelf/open sea), chromatin states (66), and DNA methylation domains (ref. 67; highly methylated domains vs. partially methylated domains) in normal liver, and replication timing in the HepG2 cell line (68). We used ELMER v2 package (69) to identify transcription factor binding sites enriched within tiles hypermethylated in the “mesenchymal” component. We used an in-house adaptation of the gene set enrichment analysis method to identify overrepresented gene sets (MSigDB v6 database) among genes paired with tiles hypermethylated in “mesenchymal” samples. We compared the vst-normalized expression of TFs involved in liver differentiation (HNF1A, HNF4A, PPARG, and NR2F6) with the methylation of tiles containing their respective binding motifs and hypermethylated in the “mesenchymal” component. These correlations were validated in Carrillo-Reixach's data set comprising 26 hepatoblastomas analyzed with Illumina 850k methylation arrays (16).

Nine pediatric liver tumor cell lines were collected from collaborations or obtained from commercial sources (Supplementary Table S8). Cells were grown in either DMEM or Advanced DMEM F-12 supplemented with 10% fetal bovine serum and using the usual conditions (100 U/mL penicillin/streptomycin, 1% glutamine at 37°C, 5% CO₂, identity confirmed using WES, Mycoplasma-free verified with the MycoAlert Mycoplasma PLUS detection kit; Lonza).

Doxorubicin (S1208, Selleck Chemicals), YM155 (S1130, Selleck Chemicals), BI-2536 (S1109, Selleck Chemicals), AZD7762 (S1532, Selleck Chemicals), trametinib (S2673, Selleck Chemicals) were dissolved in DMSO at 10 mmol/L final concentration. Cisplatin (S1166, Selleck Chemicals) was dissolved in H₂O, 0.9% NaCl, 0.3% Tween20 at a 0.5 g/L final concentration, carboplatin (S1215, Selleck Chemicals) was dissolved in H₂O, 0.1% Triton. Cells were seeded at 2,500–4,500 cells/well (Supplementary Table S10). After overnight incubation at 37°C and 5% CO₂, cells were treated with five different concentrations of drugs (0.001, 0.01, 0.1, 1, and 10 μmol/L) using HP D300 digital dispenser (Tecan). Growth inhibition was measured 72 hours after treatment with MTS diluted 1:6 in fresh culture medium. Cell viability was assessed by recording absorbance at 490 nm using a FLUOstar microplate reader. Dose–response curves were performed using GraphPad Prism 6 Software to determine two parameters reflecting drug sensitivity: GI₅₀ and the area under the curve (AUC). When the GI₅₀ was not reached, the values were set to the highest concentration tested (10 μmol/L). Each concentration was tested in duplicate. Polyclonal goat anti-IGF2 antibody (Ref AF-292-NA; Biotechne) and control anti-Igg goat antibody (Ref AB-108-C; Bio-techne) were tested using five concentrations (0.1, 0.5, 1, 5, and 10 μg/mL). Monoclonal human anti-IGF1/IGF2 therapeutic antibody (Xentuzumab, #TAB-475CQ, Creative Biolabs) was resuspended in PBS at 1 mg/mL final concentration and tested at 1 μmol/L.

Mice were housed in a specific pathogen-free facility, and experiments were conducted using protocols and conditions approved by the institutional animal ethical committee (Authorization no. 2015082610113065.01, Ethics Committee Paris-Nord C2EA 121). At day 0, 5 × 10⁶ HepG2 cells with 1:1 matrigel were inoculated subcutaneously on each flank of female BALB/cAnNRj-Foxn1 five-week-old mice. At day 20, when tumors started growing (volume >200 mm³), mice were randomized in four groups of six mice allocated to the following treatment arms: cisplatin, cisplatin/doxorubicin, BI-2536, and vehicle. Mice in the cisplatin group were injected intraperitoneally at days 20 and 27 with 5 mg/kg cisplatin. Mice from the cisplatin/doxorubicin group received one injection of cisplatin at day 20 (5 mg/kg) and one intraperitoneal doxorubicin at day 27 (2 mg/kg). Mice in the BI-2536 group were intravenously injected with 40 mg/kg BI-2536 at days 20 and 27, whereas mice in the vehicle control group received injections of H₂O, 0.9% NaCl. Mice were weighted two times a week, and tumor volume was measured at the same time using a caliper and the following formula: (length*width²)/2 and was expressed as a percentage of initial tumor volume at day 20.

To test cisplatin ability to induce SBS35 mutational signature in hepatoblastoma, six pediatric liver tumor cell lines were treated with 0.5 μmol/L cisplatin for four weeks. At day 0, 60,000 cells were seeded in a 6-well plate and reseeded each week at the same concentration. When mortality was too high, cells were not split and fresh medium supplemented with 0.5 μmol/L cisplatin was added. At day 28, cells were grown in fresh DMEM and limit dilutions were performed in order to expand clonal cell lines. After amplification, samples were extracted using the AllPrep DNA/RNA/miRNA universal extraction kit (Qiagen). Finally, mutational signature analysis was performed on cells derived from limit dilution and compared with bulk nontreated baseline mutational profile.

Multiplex ligation-dependent probe amplification (MLPA) was used to determine the status of locus 11p15.5 in 27 samples without RRBS data available. MLPA was carried out with 50–100 ng of DNA diluted in Tris-HCl, 0.1 mmol/L EDTA. DNA samples were screened using ME030-C3 BWS/RSS kit (MRC-Holland) containing 42 (methylation-specific) MLPA probes: 26 probes specific of BWS 11p15 region, 2 probes targeting NSD1 gene, 13 reference probes targeting other chromosomes, and 1 digestion control probe. Two reference probes targeting, respectively, 2p25 and 2q24 were excluded because of the presence of very frequent copy-number variation in hepatoblastoma tumors. Among the 26 probes targeting the 11p15 region, 10 were methylation-specific, allowing for the detection of methylation abnormalities. To determine thresholds for IC1 GOM and IC2 LOM, we used the K-means clustering method as well. Samples with a gain of methylation in at least two of three probes targeting IC1 were annotated as GOM IC1, whereas samples with a loss of methylation in at least two of four probes were considered LOM IC2. One probe covering IC1 (H19.11.001.976583) was excluded from analysis because its distribution among samples was not discriminating samples with and without GOM IC1.

In situ hybridization experiments were performed using RNAscope 2.5 HD Detection Reagent BROWN Kit (cat. #322300) and IGF2 probe (cat. #594361) according to the manufacturer's protocol. FFPE slides underwent target retrieval under standard conditions (15 minutes in target retrieval reagent >98°C) and a 20-minute protease digestion at 40°C in the HybEZ oven.

IHC analyses of anti–β-catenin (BD Transduction, clone 14, 610154, 1/300) and anti-GS (Bioscience, 1/500) were performed as previously described (33).

TERT promoter mutations were identified with WGS data when available and completed with Sanger sequencing as previously described (70) for other samples.

The sequencing data reported in this paper have been deposited to the EGA (European Genome–Phenome Archive) database (accession numbers EGAS EGAS00001005108 and EGAS00001003536).

T.Z. Hirsch reports grants from Cancéropôle Ile-de-France, Fondation d'Enterprise Bristol-Myers Squibb pour la Recherche en Immuno-Oncologie, and CARPEM during the conduct of the study. J. Pilet reports grants and personal fees from Fondation pour la recherche médicale during the conduct of the study. G. Morcrette reports grants from Assistance Publique Hôpitaux de Paris, CARPEM, and Ligue Nationale Contre le Cancer during the conduct of the study. A. Roehrig reports grants from Fondation pour la Recherche médicale during the conduct of the study. L. Molina reports grants from Chateaubriand Fellowship, Embassy of France in the US, Office of Science and Technology during the conduct of the study. S. Caruso reports grants from Labex OncoImmunology and CARPEM during the conduct of the study. E. Jacquemin reports personal fees from Laboratoire CTRS, France, and Vivet Therapeutics, France, outside the submitted work. J. Zucman-Rossi reports grants from France Génomique, Institut National du Cancer, SFCE, association Etoile de Martin, Fédération Enfants et Santé, association Hubert Gouin “Enfance et Cancer,” INSERM Plan Cancer, Cancéropôle Ile-de-France, Inserm ITMO Cancer 3R, Ligue Nationale contre le Cancer, Labex OncoImmunology, IREB, Fondation Bettencourt-Shueller, SIRIC CARPEM, Fondation pour la recherche médicale (FRM), Ligue Contre le Cancer Comité de Paris, Fondation Mérieux, Fondation d'Entreprise Bristol-Myers Squibb, and INCA during the conduct of the study; personal fees from CAIXA Fundation and EASL JHEP Report outside the submitted work. No disclosures were reported by the other authors.

T.Z. Hirsch: Conceptualization, investigation, visualization, writing–original draft. J. Pilet: Conceptualization, investigation, visualization, writing–review and editing. G. Morcrette: Conceptualization, visualization, resources, and investigation. A. Roehrig: Investigation and visualization. B.J.E. Monteiro: Investigation and visualization. L. Molina: Investigation. Q. Bayard: Investigation. E. Trepo: Data curation and investigation. L. Meunier: Investigation and visualization. S. Caruso: Investigation. V. Renault: Resources and data curation. J.-F. Deleuze: Resources. B. Fresneau: Resources. C. Chardot:Resources. E. Gonzales: Resources. E. Jacquemin: Resources. F. Guerin: Resources. M. Fabre: Resources. I. Aerts: Resources. S. Taque: Resources. V. Laithier: Resources. S. Branchereau: Resources. C. Guettier: Resources. L. Brugieres: Resources. S. Rebouissou: Conceptualization, supervision, and investigation. E. Letouze: Conceptualization, supervision, investigation, visualization, writing–original draft. J. Zucman-Rossi: Conceptualization, supervision, investigation, writing–original draft.

We thank all the clinicians, radiologists, pathologists, laboratory technicians, and scientists who participated in this work: Jayendra Shinde, Sandrine Imbeaud, Samantha Schaeffer, Barkha Gupta (Inserm UMR 1138); Yves Lebouc (Inserm UMR 938), Isabelle Janoueix (Inserm UMR 830); Anne Boland, Robert Olaso (CEA, CNRGH); Pr. Frédéric Gauthier, Pr. Hélène Martelli, Cécile Jeannot, Dr. Danièle Pariente, Dr. Stéphanie Franchi-Abella, Isabelle Friteau, Dr. Charlotte Mussini, Dr. Jeanine Quillard, Katia Posseme, Olivier Trassard, Martine Prsle, Véronique Bruna, Marie-Josée Redon, Véronique Mamesier, Aurélie Jove-Garcia, Damien Schmidt, Ludivine Meni, Violette Gendre, Delphine Colmant, Christophe Avignon, Julie Tisserand (CHU Bicêtre, APHP, Le Kremlin Bicêtre); Dr. Dominique Debray, Gisèle Le Gall (CHU Necker, APHP, Paris); Pr. Michel Peuchmaur, Pr. Dominique Berrebi, Pascal Blain, Muriel Mateo-Gajardo, Elodie Dauvergne, Maxette Pierin, Louise Le Guen (CHU Robert Debré, APHP, Paris); Pr. Aurore Coulomb L'hermine, Dr. Sabah Boudjemaa, Dr. Hélène Dariane (CHU Trousseau, APHP, Paris); Dr. Julien Calderaro (CHU Henri Mondor, APHP, Créteil); Pr. Jean Yves Scoazec, Brenda Mallone (Institut Gustave Roussy, Villejuif); Dr. Paul Freneaux, Dr. Daniel Orbach (Institut Curie, Paris); Dr. Carole Coze, Pr. Dominique Figarella Branger, Dr. Andre Maues de Paula (CHU de la Timone, Marseille); Pr. Valérie Rigau, Dr. Yuri Musizzano (CHU Gui de Chauliac, Montpellier); Pr. Jannick Selves (CHU Toulouse); Dr. Estelle Thebaud, Dr. Louise Galmiche, Pr. Jean-François Mosnier (CHU Nantes); Dr. Cécile Dumesnil, Dr. Arnaud François (CHU Charles Nicolle, Rouen); Dr. Anne Lutin, Dr. Antoine Gourmel, Dr. Gilles Morin (CHU Amiens); Dr. Claire Berger (CHU Saint-Etienne); Dr. Frédérique Dijoud, Dr. Cécile Picard (CHU Lyon); Dr Clémence Yguel (CHU Nancy); Charles Marcaillou (Integragen); Josh Waterfall (Institut Curie, INSERM U830), Andrei Zinovyev (Institut Curie, INSERM U900); HEPATOBIO database. We warmly thank Xin Chen, USCF (USA) and Karl-Dimiter Bissig, Baylor College, for their kind gift of the Hep293TT and B6-2 cell lines. This work was supported by France Génomique (GEPELIN project), Institut National du Cancer (INCa, PELICAN.resist project), SFCE, association Etoile de Martin, Fédération Enfants et Santé (FES), association Hubert Gouin “Enfance et Cancer,” INSERM with the HTE (plan Cancer), Cancéropôle Ile-de-France “ThéBioCaPe” project, Inserm ITMO Cancer 3R “PEILICANS” project. The group is supported by the Ligue Nationale contre le Cancer (Equipe Labellisée), Labex OncoImmunology (investissement d'avenir), grant IREB, Coup d'Elan de la Fondation Bettencourt-Shueller, the SIRIC CARPEM, FRM prix Rosen, Ligue Contre le Cancer Comité de Paris (prix René et Andre Duquesne), and Fondation Mérieux. T.Z. Hirsch was supported by a fellowship from Cancéropôle Ile-de-France, Fondation d'Entreprise Bristol-Myers Squibb pour la Recherche en Immuno-Oncologie, and CARPEM; G. Morcrette was supported by CARPEM, the Ligue Nationale contre le Cancer and Assistance Publique Hôpitaux de Paris. Q. Bayard was funded by a fellowship from the Ministry of Education and Research. This work was supported by the Fondation pour la Recherche Médicale, grant number ECO201906008977 to A. Roehrig and grant number ECO20170637540 to J. Pilet. L. Molina was funded by the Chateaubriand Fellowship from the Embassy of France in the USA. S. Caruso was funded by Labex OncoImmunology and CARPEM. HEPATOBIO was supported by SFCE, Etoile de Martin, Enfant Sans Cancer, INCA/DGOS (Grant INCA-TRANSLA N° 2013-209).

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.