Research Article Hepatic and Biliary Cancer Epigenetic footprint enables molecular risk stratification of hepatoblastoma with clinical implications Graphical abstract Authors 113 HB patients Juan Carrillo-Reixach, Laura Torrens, Marina Simon-Coma, ., Maria Rosa Sarrias, Josep M. Llovet, Carolina Armengol GENOMIC STUDY TRANSCRIPTOMIC STUDY METHYLATION STUDY RNA-sequencing/ HTA® -array/ CytoScan HD® -array ddPCR RT-qPCR 850K (EPIC)-array/ QUAlu Correspondence carmengol@igtp.cat Dysregulation Overexpression 2 epigenomic HB subtypes (C. Armengol). of global RNA of 14q32 & BLCAP editing DLK1-DIO3 genes (Epi-CA & Epi-CB) 16 + VIM-gene signature (C1/C2/C2B) Lay summary Hepatoblastoma is a rare CLINICAL PARAMETERS: prognostic marker identification childhood liver cancer that has been understudied. CHKA new therapeutic Molecular risk stratification We have used cutting- target Poor prognostic factors: MRS1 MRS2 MRS3 edge technologies to Survival -4q, -18, 17q11.2 AI (NF1) expand our molecular Strong 14q32 Epi-CB knowledge of this cancer. Time Our biological findings can be used to improve clinical management and pave the Highlights way for the development  Hepatoblastoma (HB) involves global dysregulation of RNA editing, including of novel therapies for this in the tumor suppressor BLCAP. cancer.  Overexpression of a 300 kb region within the 14q32 DLK1/DIO3 locus is a new hallmark of HB.  We identified 2 epigenomic HB subtypes -Epi-CA and Epi-CB- with distinct degrees of DNA hypomethylation and CpG island hypermethylation.  The molecular risk stratification of HB, based on the 14q32-signature and epigenomic subtypes, is associated with patient outcomes.  The enzyme CHKA could be a novel therapeutic target for patients with HB. https://doi.org/10.1016/j.jhep.2020.03.025 © 2020 European Association for the Study of the Liver. Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/). J. Hepatol. 2020, 73, 328–341 Research Article Hepatic and Biliary Cancer Epigenetic footprint enables molecular risk stratification of hepatoblastoma with clinical implications Juan Carrillo-Reixach1, Laura Torrens2,3, Marina Simon-Coma1,4, Laura Royo1, Montserrat Domingo-Sàbat1, Jordi Abril-Fornaguera2,3, Nicholas Akers3,5, Margarita Sala4,6,7, Sonia Ragull1, Magdalena Arnal8, Núria Villalmanzo9, Stefano Cairo10,11, Alberto Villanueva12, Roland Kappler13, Marta Garrido14, Laura Guerra15, Constantino Sábado16, Gabriela Guillén17, Mar Mallo18, David Piñeyro19, María Vázquez-Vitali1, Olga Kuchuk2, María Elena Mateos20, Gema Ramírez21, Manuel López Santamaría22, Yasmina Mozo23, Aroa Soriano24, Michael Grotzer25, Sophie Branchereau26, Nagore García de Andoin27, Blanca López-Ibor28, Ricardo López-Almaraz29, José Antonio Salinas30, Bárbara Torres31, Francisco Hernández22, José Javier Uriz27, Monique Fabre32, Julià Blanco33,34, Claudia Paris35, Viera Bajciová36, Geneviève Laureys37, Helena Masnou6, Ariadna Clos6, Cristina Belendez38, Catherine Guettier26, Lauro Sumoy19, Ramón Planas4,6, Mireia Jordà9,39, Lara Nonell8, Piotr Czauderna40, Bruce Morland41, Daniela Sia2, Bojan Losic5,42,43,44, Marie Annick Buendia45, Maria Rosa Sarrias4,46, Josep M. Llovet2,3,47, Carolina Armengol1,4,* 1 Childhood Liver Oncology Group, Germans Trias i Pujol Research Institute (IGTP), Program for Predictive and Personalized Medicine of Cancer (PMPPC), Badalona, Spain; 2Mount Sinai Liver Cancer Program, Divisions of Liver Diseases, Tisch Cancer Institute, Icahn School of Medicine at Mount Sinai, New York, USA; 3Translational research in Hepatic Oncology, Liver Unit, Institut d’Investigacions Biomèdiques August Pi i Sunyer (IDIBAPS), Hospital Clínic, Universitat de Barcelona, Barcelona, Spain; 4CIBER, Hepatic and Digestive Diseases, Barcelona, Spain; 5Department of Genetics and Genomic Sciences, The Icahn Institute for Genomics and Multiscale Biology, Icahn School of Medicine at Mount Sinai, New York, USA; 6Gastroenterology Department, Hospital Universitari Germans Trias i Pujol Hospital, Badalona, Spain; 7Gastroenterology Department, Hospital Universitari Josep Trueta, Girona, Spain; 8MARGenomics, IMIM (Hospital del Mar Medical Research Institute), Barcelona, Spain; 9PMPPC, IGTP, Badalona, Spain; 10XenTech, Evry, France; 11Instituto di Ricerca Pediatrica, Corso Stati Uniti 4, Padova, Italy; 12 Chemoresistance and Predictive Factors Group, Program Against Cancer Therapeutic Resistance, Catalan Institute of Oncology (ICO), Bellvitge Biomedical Research Institute, L’Hospitalet del Llobregat, Barcelona, Spain; 13Department of Pediatric Surgery, Dr. von Hauner Children’s Hospital, Ludwig-Maximilians-University Munich, Lindwurmstr. 2a, 80337 Munich, Germany; 14Hospital Vall d’Hebron, Pathology Department, Barcelona, Spain; 15University Hospital La Paz, Pathology Department, Madrid, Spain; 16Hospital Vall d’Hebron, Pediatric Oncology Department, Barcelona, Spain; 17Hospital Vall d’Hebron, Pediatric Surgery Department, Barcelona, Spain; 18MDS Research Group, Josep Carreras Leukaemia Research Institute, ICO-Hospital Germans Trias i Pujol, Universitat Autònoma de Barcelona, Badalona, Spain; 19High Content Genomics and Bioinformatics Unit, PMPPC, IGTP, Badalona, Spain; 20Pediatric Oncology Unit, Department of Pediatrics, University Hospital Reina Sofía, Córdoba, Spain; 21University Hospital Universitario Virgen del Rocío, Pediatric Oncology Department, Sevilla, Spain; 22 University Hospital La Paz, Pediatric Surgery Department, Madrid, Spain; 23University Hospital La Paz, Pediatric Oncology Department, Madrid, Spain; 24Biomedical Research in Cancer Stem Cells Group, Pathology Department, Institut de Recerca Hospital Vall d’Hebron (VHIR), Barcelona, Spain; 25Department of Pediatric Oncology, University Children’s Hospital Zurich, University of Zurich, Zurich, Switzerland; 26 Bicêtre Hospital, Le Kremlin-Bicêtre, France; 27Pediatric Oncology, Hospital Universitario Donostia, Doctor Begiristain Kalea, 117, 20080, Donostia, Spain; 28Department of Pediatric Hematology and Oncology, HM Montepríncipe Hospital, Boadilla del Monte, Madrid, Spain; 29 Pediatric Oncology and Hematology, Hospital Universitario Cruces, Bilbao, Spain; 30Division of Hematology-Oncology, Department of Pediatrics, Hospital Universitari Son Espases, Palma de Mallorca, Spain; 31Medical Oncology Department, Pediatric Oncology Department, University Hospital La Fe, Valencia, Spain; 32Department of Pathology, Hôpital Universitaire Necker-Enfants Malades, Assistance Publique- Hôpitaux de Paris, and Université Paris Descartes, Paris, France; 33IrsiCaixa AIDS Research Institute, IGTP, Badalona, Catalonia, Spain; 34 University of Vic - Central University of Catalonia, Vic, Spain; 35Stem Cell Transplant Unit, Hospital Luis Calvo Mackenna, Santiago, Chile; 36 Department of Pediatric Oncology, Childrens University Hospital Brno, Brno, Czech; 37Department of Pediatric Hematology, Oncology and Hematopoietic Stem Cell Transplantation, Ghent University Hospital, Ghent, Belgium; 38Oncohematology Service, Hospital Gregorio Marañón, Madrid, Spain; 39Consortium for the Study of Thyroid Cancer, CECaT, Barcelona, Spain; 40Department of Surgery and Urology for Children and Adolescents, Medical University of Gdansk, Gdansk, Poland; 41Department of Oncology, Birmingham Women’s and Children’s Hospital, Birmingham, UK; 42Graduate School of Biomedical Sciences, One Gustave L. Levy Place, Box 1022, New York, NY; 43Tisch Cancer Institute, Cancer Immunology Program, Icahn School of Medicine at Mount Sinai, New York, NY; 44Icahn Institute for Data Science and Genomic Keywords: Hepatoblastoma (HB); RNA editing; BLCAP; 14q32; DLK1-DIO3 locus; de Can Ruti. Camí de les Escoles, s/n 08916 Badalona, Spain. Tel.: + 34 93 554 30 72, Molecular risk stratification; Prognostic biomarker; CHKA. fax: +34 93 497 86 54. Received 5 September 2019; received in revised form 11 March 2020; accepted 13 March E-mail address: carmengol@igtp.cat (C. Armengol). 2020; available online 30 March 2020 https://doi.org/10.1016/j.jhep.2020.03.025 *Corresponding author. Address: Germans Trias i Pujol Research Institute (IGTP), Edifici Muntanya, Campus Can Ruti, Childhood Liver Oncology Group (c-LOG), Ctra. Journal of Hepatology 2020 vol. 73 j 328–341 Technology, Icahn School of Medicine at Mount Sinai, New York, NY; 45INSERM, UMR 1193, Paul-Brousse Hospital, Hepatobiliary Centre, Villejuif, France; 46Innate Immunity Group, IGTP, Badalona, Spain; 47Institució Catalana de Recerca i Estudis Avançats (ICREA), Barcelona, Spain Background & Aims: Hepatoblastoma (HB) is a rare disease. stratification in children with HB.6 However, the rarity of the dis- Nevertheless, it is the predominant pediatric liver cancer, ease has impaired the incorporation of molecular data into this with limited therapeutic options for patients with aggressive clinical classification. In this context, there is a need to increase our tumors. Herein, we aimed to uncover the mechanisms of HB understanding of the biology of this rare tumor and its prognostic pathobiology and to identify new biomarkers and therapeutic determinants to be able to move towards biology-driven precision targets in a move towards precision medicine for patients with medicine, which includes biomarkers for therapeutic tailoring. advanced HB. The origin of HB is largely unknown. Most tumors are spo- Methods: We performed a comprehensive genomic, tran- radic, and their extreme rarity has limited our understanding of scriptomic and epigenomic characterization of 159 clinically their underlying molecular mechanisms. Regarding the genetic annotated samples from 113 patients with HB, using high- alterations identified to date, the most significant are activating throughput technologies. mutations of the catenin beta 1 (CTNNB1) gene, which encodes Results: We discovered a widespread epigenetic footprint of HB b-catenin, in more than 70% of HBs.7 b-catenin is a key regulator that includes hyperediting of the tumor suppressor BLCAP of cell fate and proliferation during liver development and concomitant with a genome-wide dysregulation of RNA editing regeneration. CTNNB1 mutations in cancer impair the proteoso- and the overexpression of mainly non-coding genes of the onco- mal degradation of this protein and lead to the constitutive genic 14q32 DLK1-DIO3 locus. By unsupervised analysis, we iden- activation of the Wnt pathway.8 tified 2 epigenomic clusters (Epi-CA, Epi-CB) with distinct degrees High-throughput technologies now enable us to identify the of DNA hypomethylation and CpG island hypermethylation that are molecular subtypes of diverse cancers and their associated associated with the C1/C2/C2B transcriptomic subtypes. Based on oncogenic aberrations. Based on transcriptomic studies, we these findings, we defined the first molecular risk stratification of identified 2 HB subclasses—C1 and C2—that resemble late and HB (MRS-HB), which encompasses 3 main prognostic categories early stages of liver development, and a discriminating 16-gene and improves the current clinical risk stratification approach. The signature.9 The recent studies led by French10 and American11 MRS-3 category (28%), defined by strong 14q32 locus expression teams described a third HB subclass not detected by the 16- and Epi-CB methylation features, was characterized by CTNNB1 and gene signature. This subclass, called C2B in the paper by Hooks NFE2L2 mutations, a progenitor-like phenotype and clinical et al.,10 is characterized by increased expression of epithelial- aggressiveness. Finally, we identified choline kinase alpha as a mesenchymal transition markers such as vimentin (VIM). A promising therapeutic target for intermediate and high-risk HBs, as recent pan-cancer analysis showed that HB is the tumor with the its inhibition in HB cell lines and patient-derived xenografts lowest rate of somatic mutations (1–7 mutations per tumor strongly abrogated tumor growth. genome).12 However, exome sequencing studies of HB have Conclusions: These findings provide a detailed insight into the revealed nuclear factor erythroid 2-related factor 2 (NFE2L2), a molecular features of HB and could be used to improve current regulator of critical antioxidant and stress-responsive genes, as clinical stratification approaches and to develop treatments for the second most mutated gene in ~10% of cases.13 In comparison patients with HB. with hepatocellular carcinoma (HCC), the main liver cancer in Lay summary: Hepatoblastoma is a rare childhood liver cancer that adults, HB has more than 10-fold fewer mutations. This obser- has been understudied. We have used cutting-edge technologies to vation thus suggests that childhood liver tumorigenesis is driven expand our molecular knowledge of this cancer. Our biological by mechanisms other than DNA mutations, such as epigenetic findings can be used to improve clinical management and pave the modifications.14 To date, genome-wide epigenetic studies on HB way for the development of novel therapies for this cancer. are scarce and included a limited number of cases.15–17 © 2020 European Association for the Study of the Liver. Published by Through a high-throughput genomic, transcriptomic and epi- Elsevier B.V. This is an open access article under the CC BY-NC-ND genomic study of unprecedented size, we have discovered and license (http://creativecommons.org/licenses/by-nc-nd/4.0/). validated a profound epigenetic footprint in HB, spanning RNA editing dysregulation to specific DNA methylation profiles linked Introduction to strong overexpression of 14q32 genes, Wnt signaling, and a Hepatoblastoma (HB) is the predominant pediatric liver tumor, progenitor-like phenotype. Based on our findings, which includes mainly affecting infants under 3 years of age.1 Although its inci- an updated 16-gene signature, we present the first molecular risk dence has increased markedly over the last 30 years, HB is a rare stratification of HB, which seeks to improve on the current clinical disease (1.8 cases per million children per year).2,3 Clinical studies CHIC risk staging system, and we identify choline kinase alpha combining chemotherapy and efficient surgical approaches have (CHKA) as a potential therapeutic target for patients with HB. led to dramatic improvements in outcomes for patients with HB, with a 3-year event-free survival (EFS) above 80%.4 However, there are limited treatment options for clinically advanced tumors, with Patients and methods a 3-year EFS of only 34%.5 Furthermore, patient survivors can suffer Patients and samples severe and lifelong side effects derived from chemotherapy and The study included 113 patients with HB (discovery set: 67 immunosuppression. A recent unified analysis from the Children's samples, 33 patients; validation set: 92 samples, 80 patients). In Hepatic tumors International Collaboration (CHIC) led to the total, we analyzed 112 primary tumors, 3 recurrences and 44 development of a new international clinical staging system for risk paired non-tumor samples (Table S1). The main clinical Journal of Hepatology 2020 vol. 73 j 328–341 329 Research Article Hepatic and Biliary Cancer Table 1. Main clinical and pathological features of the 113 patients with detected in the complete set (total incidence <1%), thereby ruling hepatoblastoma included in the study. out their relevance for HB tumorigenesis. Discovery set Validation set We analyzed the genomic profiling of the same tumors with the (n = 33) (n = 80) high-resolution CytoScan HD array. The recurrent altered chro- Age, months (median, [range]) 16 [1–180] 18 [0.2–204] mosomal regions in HB are shown in Fig. 1, confirming previous Gender (M/F) 20/13 52/28 single-nucleotide polymorphism array- or karyotype-based re- Serum AFP, ng/ml (range) 341–2,186,461 300–12,299,925 ports.9,11,19 The most frequent chromosomal alterations included Clinical classification: CHIC-HS (VL-L/I/H) 16/4/13 48/7/25 broad and focal copy number gains in 1q, 2q, 5p, 6, 7, 8, 12, 13q34, Pre-operative chemotherapy 31/2 (94%) 73/7 (91%) 15q, 17q and 20, and losses in 1p, 4q and 18p11.32 (Fig. S2). The (Y/N, %) most recurrent, already reported,20 allelic imbalance involved 3.9 Tumor characteristics: Mb of the 11p15 locus (13/32, 41% cases). Other additional recur- PRETEXT stage (I/II/III/IV/n.a.) 2/12/13/6/0 4/31/28/16/1 rent allelic imbalances were found in 1p, 2q, 2p, 3p, 7q,11p and 17q Vascular Invasion (Y/N, %) 13/20 (39%) 16/64 (20%) Multifocality (Y/N, %) 12/21 (36%) 19/61 (24%) in 13–22% of the primary tumors, of which the last two, to our Metastasis at diagnosis 9/24 (27%) 18/62 (22.5%) knowledge, have not been found in previous studies. Tumors from (Y/N, %) recurrences or with a HCN-NOS (hepatocellular malignant neo- HB histology: plasms not otherwise specified21) histology showed an increased Epithelial/Mixed/n.a. 17/16/0 54/25/1 number of chromosomal aberrations (Fig. S2). MEC: Fetal/Non-Fetala/n.a. 19/13/1 65/14/1 HCN-NOS 2 – Follow-up, months (mean, 41.76 [1–100] 41.45 [0,2–100] Genome-wide dysregulation of RNA editing and BLCAP [range]) hyperediting in HB Outcome: cancer-related 11/22 (33%) 17/63 (21%) Because BLCAP RNA is a highly conserved edited transcript,22 we deaths or tumor recurrence (Y/N, %) examined whether the nt 5 A>G substitution (which confers a a Non-fetal includes crowded fetal, macrotrabecular and embryonal histological Y2C change) observed in the RNA-seq data is due to an editing subtypes. CHIC-HS, Children's Hepatic tumors International Collaboration-Hepato- blastoma Stratification (VL-L, very low; or low; I, intermediate and H, high risk)6; event. RT-PCR-Sanger sequencing revealed that the nt 5 substi- HCN-NOS, hepatocellular neoplasm not otherwise specified; MEC, main epithelial tution was present only in the RNA and not in the DNA of the component; n.a., non-available; PRETEXT, PRETreatment EXTent of disease. tumors (Fig. 2A). This observation strongly points to the dysre- gulation of RNA editing in the BLCAP transcript. The nt 5 editing characteristics of these patients as well as the pathological and of BLCAP was further confirmed by droplet digital PCR. To this molecular (CTNNB1 status and C1/C2 classification) features of end, we used probes to measure the fractional abundance of the tumors are summarized in Table 1. wild-type and edited nt 5 of BLCAP and found that the latter was 1.85-fold higher in tumor than in non-tumor samples (p <0.0001, Fig. 2B). Molecular profiling These findings on BLCAP prompted us to study whether RNA RNA-sequencing (RNA-seq), human transcriptome array (HTA), editing is globally disrupted in HB. Genome-wide analysis of RNA CytoScan HD and methylation 850K-arrays were performed on changes using RNA-seq data revealed that tumor samples had a the discovery set. The main findings were confirmed in the lower overall editing index than non-tumor samples in both Alu validation set and 5 human fetal livers. Sample-assay overlap is and non-Alu regions (p <0.0001, Fig. 2C). detailed in Table S1. The omics data generated in this study have Adenosine deaminases acting on RNA (ADARs) are respon- been deposited in NCBI's Gene Expression Omnibus18 and are sible for converting A to I in nuclear-encoded RNAs, leading to available through GEO Series accession number GSE132219. A>G substitutions.23 We therefore studied whether these en- Whole-genome sequencing data are available under accession zymes were aberrantly expressed in HB and found a significant number PRJNA548663 at the Sequence Read Archive of the NCBI. overexpression of both ADAR1 and ADAR2 genes in tumor Additional detailed protocols are provided in supplementary compared to non-tumor samples (p < −0.0005; Fig. 2D). In sum- information and CTAT table. mary, we discovered an unprecedented dysregulation of global editing and regulatory enzymes and identified BLCAP as the first Results hyperedited gene in HB. Genomic profiling reveals a recurrent altered sequence of BLCAP in HB Overexpression of 14q32 genes is a new hallmark of HB RNA-seq data were examined for nucleotide alterations that lead By performing an unsupervised hierarchical clustering of the to amino acid changes and fusion transcripts. The study of transcriptomic data, we identified 3 groups of tumors according missense changes found the same point CTNNB1 mutations as to the co-clustering study of 12 dendrograms (Fig. S3A). Co- identified by RT-Sanger-sequencing (Table S1), and revealed cluster 1 (CC1) and co-cluster 2 (CC2) were enriched in C1 and changes (A/G) in nucleotide (nt) positions 5, 14 and/or 44 of the C2 tumors (p = 0.001) whereas the third co-cluster (CC3) apoptosis-inducing factor BLCAP (bladder cancer-associated composed of C2 tumors, significantly overlapped with the protein) transcript in 9 cases of the discovery set (28%) (Fig. 1). recently identified C2B subclass10 (Fig. S3B) which had high NFE2L2 mutations were found in 3 cases (9%). Analysis of fusion expression of VIM (Fig. S3C). The tumor gene expression profile transcripts identified 15 events with perfect alignment in 12 showed upregulation of the Wnt/b-catenin pathway, and distinct tumors (Table S2). Four of these transcripts were imprinted and stem cell-related genes (fold change [FC] >2, false selected and validated in tumor samples and their corresponding discovery rate [FDR] <0.001, Table S3A,B). Among the most patient-derived xenografts (PDX) by RT-PCR-Sanger sequencing strongly dysregulated genes in tumors, as identified by HTA (Fig. S1). No additional tumors with these fusion transcripts were (FC >30, FDR <10−7), we also found a previously undescribed, 330 Journal of Hepatology 2020 vol. 73 j 328–341 Frequency (%) A R19* T30* R27 R32 T15 T18 T10 T14 T19 T12 T29 T21 T16 T28 T22 T20 T13 T33 T31 T25 T34 T17 T23 T11 T32 T26 0 20 40 60 80 T9 T4 T6 T8 T5 T2 T7 T1 CTNNB1 mutations (69%) BLCAP nt 5 A>G (28%)* 1 1 1 1 1 1 1 1 1 1 1 Fusion transcripts (12%) 1 1 1 1 NFE2L2 mutation (9%) 1 1 1 1 G 20p11.1-q13.33 (53%) G 1q25.1-44 (50%) G 2q12.1-37.3 (47%) 1 G 8p23.1-q24.3 (47%) AI 11p15.4-15.5 (41%) 1 1 1 1 1 1 1 1 1 1 1 1 1 1 G 12p11.1-q14.3 (38%) G 17q11.1-q25.3 (38%) AI 1p32.3-43 (28%) 1 1 1 1 1 1 1 1 1 1 AI 11p12-11.2 (28%) 1 1 1 1 1 1 1 1 1 1 G 6p11.2-q14.1 (28%) L 4q32.3-q35.2 (25%) 1 1 1 1 1 1 1 1 1 L 1p12-p36.33 (22%) 1 1 1 1 1 1 1 1 G 15q12-q22.2 (19%) G 7p11.2-q36.3 (19%) G 5p12-35.3 (16%) G 13q34 (16%) L 18p11.32 (16%) 1 1 1 1 1 1 AI 3p21.2-p21.31 (16%) 1 1 1 1 1 1 1 AI 2p11.2-q11.2 (13%) 1 1 1 1 1 1 AI 7q11.21 (13%) 1 1 1 1 1 1 AI 17q11.2 (13%) 1 1 1 1 1 AI 2q36.3-q37.1 (13%) 1 1 1 1 1 B 20 N. alterations/sample Presence of gene alterations No Yes 10 Copy number alterations Gain (G) Loss (L) Allelic imbalance (AI) 0 Fig. 1. The genomic landscape of hepatoblastoma. (A) Main gene mutations, fusion transcripts, and chromosomal alterations in the 34 tumor (T) samples from the 32 patients (discovery set) ranked by their frequency (right plot) and considering only 1 sample per patient. CTNNB1 and NFE2L2 mutations and BLCAP nucleotide 5 change were identified by RNA-seq. Only copy number alterations and allelic imbalances identified by the CytoScan HD array in at least 4 primary tumors are shown. HBs with an HCN-NOS histology are marked with an asterisk. (B) Total number of main aberrations in the above tumor samples. X-axis indicates the individual tumors; Y-axis indicates the number of total gene mutations (black), copy number alterations (gains, blue; losses, red) and allelic imbalances (pink) per sample. HB, hepatoblastoma; HCN-NOS, hepatocellular neoplasm not otherwise specified; R, recurrence; RNA-Seq, RNA sequencing. highly upregulated 330 kb region within the DLK1 (delta-like (p <0.0001, Fig. 3C). We assessed the mRNA expression of 14q32 non-canonical notch ligand)-DIO3 (iodothyronine deiodinase 3) in the validation set and confirmed the overexpression of 14q32 locus, spanning from 14q32.2 to q32.31, called 14q32 henceforth genes. Its correlation with the Wnt/b-catenin pathway activation (Fig. 3A, see Fig. S4A,B for more details). This is an imprinted was also confirmed in the validation set by measuring LGR5, a region with a key role in human development and cancer.24,25 It well-known wnt/b-catenin target gene27 (Fig. S5). Moreover, the contains more than 100 transcripts, including DLK1 (a well- study of fetal liver samples also indicated an elevated expression known hepatoblast marker highly expressed in HB9,26), MEG3 of 14q32 genes, thereby reinforcing the idea that HB re- and MEG8 (2 maternally expressed non-coding genes), small capitulates pathological and molecular features of developing nucleolar RNAs of the C/D box family (namely SNORD113 and livers.9,28 SNORD114), and the largest microRNA cluster in the human To gain insight into the possible mechanisms conferring genome (Fig. 3A). Tumor overexpression of 14q32 genes was strong 14q32 gene expression in HB, we examined the 14q32 further validated by gene set enrichment analysis (GSEA) (Fig. 3B, region at the genomic and epigenomic level. Since 14q32 gene Table S4). overexpression has previously been linked to adeno-associated Hierarchical clustering based on the gene expression profile virus integration in this locus and hepatocarcinogenesis,29,30 of all the genes localized at 14q32 showed 2 main groups of we used whole genome sequencing to search for viral inte- tumors, thus revealing variability among HBs (Fig. 3C). This gration in a tumor with a strong 14q32-gene signature. No heterogeneity was observed at the level of tumor/non-tumor viral integration site was detected (Table S5, Fig. S6). Neither expression of 14q32 genes and in the number of overexpressed did the CytoScan HD array reveal focal chromosomal rear- genes (Fig. S4C). Among these genes, 4 (DLK1, MEG3, SNORD113- rangements. In contrast, the comparison of tumor and non- 3, SNORD114-22) were selected to classify tumors on the basis of tumor methylation profiles using the 850K-array identified the degree of 14q32 gene expression (strong/moderate), and 32 significant differently methylated CpGs localized at the they are referred to as the 14q32-gene signature hereafter 14q32 locus with predominant tumor DNA hypomethylation (Fig. S4D). Strikingly, the resulting 14q32 classification was also (FDR <0.05), specifically, the methylation levels of a CpG associated with mutations in the Wnt/b-catenin pathway (cg02412314) localized within the intragenic MEG3 region Journal of Hepatology 2020 vol. 73 j 328–341 331 Research Article Hepatic and Biliary Cancer A DNA RNA B C C A T G T A T T G C A T G T A T T G ddPCR RNA-seq 30 0.4 p <0.0001 p <0.0001 NT BLCAP nt 5 editing (FA) 0.3 Global editing index 20 wt Tyr2Cys (c.5A>G) 0.2 C A T G T A T T G C A T G T A/G T T G 10 0.1 T 0 0.0 NT T NT T wt Tyr2Cys (c.5A>G) D 8,000 HTA 300 400 RNA-seq 60 p = 0.0003 p <0.0001 p = 0.0005 p <0.0001 gene expression gene expression gene expression gene expression 6,000 300 200 40 ADAR2 ADAR2 ADAR1 ADAR1 4,000 200 100 20 2,000 100 0 0 0 0 NT T NT T NT T NT T Fig. 2. Dysregulation of RNA editing in hepatoblastoma. (A) Chromatogram of DNA and RNA Sanger sequences of the BLCAP gene in tumor (T) and non-tumor (NT) samples of a representative case with RNA editing of nucleotide 5 (highlighted in yellow). The black arrow indicates the ATG start codon. (B) FA of nucleotide 5 edited vs. non-edited BLCAP assessed by ddPCR in the 31 paired T and NT samples (discovery set) for which RNA was available (paired t test). (C) Global editing index in the 32 cases of the discovery set determined by RNA-seq (paired t test). (D) Gene expression of ADAR1 and ADAR2 genes in T and NT samples (Mann- Whitney U test). HTA plot includes data of 18 NT and 32 T samples; RNA-seq plot includes data of 32 NT and 32 T samples. Gene expression is given in normalized arbitrary units (HTA array) or counts (RNA-seq). BLCAP, bladder cancer associated protein; ddPCR, droplet digital PCR; FA, fractional abundance; HTA, human transcriptome array; RNA-seq, RNA sequencing. showed a strong correlation with the mean gene expression of classification9 (p = 0.6882) but were strongly associated with our the 14q32 region (r = −0.61, p <0.0001, Fig. 3D). Interestingly, CC1/CC2/CC3 transcriptomic co-clusters (p <−0.0005) and with the tumors with a strong 14q32-gene signature displayed lower Hooks signature10 (p <0.005). Moreover, the Epi-CB cluster, which methylation of the 32 CpGs than tumors with a moderate was enriched with tumors of the C2 subtype, exhibited consti- signature (p = 0.0082, Fig. 3E). Moreover, fetal liver samples tutive activation of Wnt/b-catenin signaling (p = 0.0391) and a showed the lowest levels of 14q32 methylation compared with strong 14q32-gene signature (p = 0.0010) (Fig. 4B). The study of tumor samples (Fig. 3E). the methylation profiles between the 2 tumor clusters revealed that Epi-CB tumors had a sharp global hypomethylation Identification of two distinct epigenetic profiles in HB compared to Epi-CA tumors in all epigenomic structures, except Next, we used the 850K-array data to extend our methylation for CpG islands, which were hypermethylated (Fig. 4C). We next study to the complete genome. Principal component analysis of studied the impact of this specific CpG island hypermethylation the methylation data showed that tumor samples were clearly on the transcriptome of Epi-CB tumors and identified KLF6, ITGB3, distinct from non-tumor samples (Fig. 4A). In general, tumors NFIC, TRANK1 and TSPYL5 as possible tumor suppressor genes, as were characterized by genome-wide DNA hypomethylation the hypermethylation of their CpG islands (b >0.2 and FDR (p <0.0001). The supervised analysis comparing tumor and non- <0.0001) was associated with a switch-off of their expression tumor samples identified 30,165 differently methylated CpGs (RNA-seq/HTA: FC <−2 and FDR <0.001, Fig. 4D). Next, we sought regulating 7,234 genes (|Db| >0.20, FDR <0.0001), including to investigate whether the dysregulation of methylation observed hypermethylation of the RASSF1 promoter.31 Using HTA and in HB was associated with changes in the expression of the RNA-seq to associate these data with gene expression, we found enzymes regulating DNA methylation. The expression of tet 21 hypermethylated (Db >0.20, FDR <0.0001) genes in tumors methylcytosine (TET) family genes — specifically TET1 and TET3 — with concomitant reduced gene expression (FC <−2, FDR <0.05) involved in DNA demethylation, was significantly higher in tu- (Table S6). Among them, we recognized genes endowed with mors compared to non-tumor tissue and correlated with the tumor suppressor functions such as AKR7A3,32 EDNRB,33 ESRP2,34 degree of hypomethylation (Fig. S7A). Therefore, Epi-CB tumors PEMT35 and PER3.36 with strong global hypomethylation had significantly higher The unsupervised analysis showed separate clusters of tumor levels of TET1 and TET3 than Epi-CA tumors (p <0.0025). Similarly, and non-tumor samples (p <0.0001). The epigenetic clustering the expression levels of DNA methyltransferases (DNMT1, also disclosed two distinct tumor clusters, which we called DNMT3A and DNMT3B) were higher in tumors than non-tumor Epigenetic-Cluster A and Epigenetic-Cluster B (Epi-CA and Epi- samples (p <0.0001) and mainly in Epi-CB tumors characterized CB) (Fig. 4B) which were not associated with the 16-gene C1/C2 by CpG island methylation. Moreover, expression of DNMTs was 332 Journal of Hepatology 2020 vol. 73 j 328–341 also correlated with high level of CpG island methylation CHKA as a new therapeutic target for intermediate and high- (Fig. S7B). risk HBs As we reported that C1/C2 subtypes resembled late and early To identify therapeutic targets for aggressive HBs, we performed a liver developmental stages,9 we examined whether the 2 distinct supervised analysis comparing the 3 main molecular risk cate- tumor methylation profiles also mimicked the different gories. Among the 392 differentially expressed genes (FDR methylation profiles of adjacent non-tumor and fetal liver sam- <0.0001, Table S9), we observed overexpression of 14q32 tran- ples based of different ages. In line with our previous findings,9 scripts (DLK1, MEG3, SNORD113-4 and SNORD114-13) and liver we observed that the global methylation value of Epi-CB tu- progenitor markers (GPC3, KRT19, AFP, EPCAM) in tumors mors, enriched with C2 tumors, was similar to that of early belonging to the high-risk category (MRS-3) compared to tumors embryonal/fetal phases of liver development at 8.4 ± 7.7 weeks of in the low-risk (MRS-1), and to a lesser extent, to those in the gestation, while the Epi-CA tumors, enriched with C1 and C2B intermediate-risk (MRS-2) categories. The most widely overex- tumors, had a global methylation value similar to that of late fetal pressed coding gene in high-risk and intermediate-risk tumors or postnatal liver phases at ~5 weeks after birth (Fig. 4E). was CHKA (Table S9, Fig. 6A), the main regulator of the biosyn- thesis of phosphatidylcholine via the CDP-choline pathway, which plays a key role in regulating cell growth and carcinogenesis.37 The Molecular risk stratification of HB differential expression of CHKA between MRS categories (Fig. 6B), To address the relevance of our molecular findings in the clinical as well as proliferation (Ki67), 14q32 (DLK1) and liver progenitor setting, we studied their association with clinical parameters in (EpCAM, GPC3, KRT19, AFP) markers was also seen by immuno- the whole set of 113 patients (77% CTNNB1 mutations, 4% NFE2L2 histochemistry (IHC) in 20 tumors (Fig. S9). mutations, 25% BLCAP nt 5 hyperediting, 63% strong 14q32-gene To address whether CHKA could be used as a therapeutic signature, 33% Epi-CB, and 44% C2; Table S7). Moreover, we target for HB, we tested the anti-tumoral ability of two CHKA measured VIM expression in order to determine its impact on our inhibitors (MN58b and TCD-717)38 at growing concentrations previous 16-gene signature9 and defined an updated 16+VIM- (2–8 lM) in two HB cell lines, HepG2 and Huh6. Both CHKA gene signature that classified the C2 tumors as either C2B10 (11%) inhibitors exerted a dose-dependent reduction of cell viability in or C2-Pure (34%) on the basis of high or low levels of VIM, the two cell lines (p = 0.0005 and p <0.0001, respectively) respectively. The association of these molecular features with (Fig. 6C). Similarly, MN58b and TCD-717 completely inhibited clinical data revealed that the losses of chromosome 4q or 18 and colony formation in HepG2 and Huh6 cells (p <0.0001, Fig. 6C). the 17q11.2 allelic imbalance (where the tumor suppressor We also investigated the anti-tumoral effects of silencing CHKA neurofibromin 1, NF1, is localized) were associated with poor gene expression via small interfering RNA in the Huh6 cell line, prognostic parameters (see more details Table S8). Moreover, which has the lower expression of this enzyme. After depleting patients with tumors with a strong 14q32-gene signature or CHKA by ~4-fold, cell viability of Huh6 cells was reduced by ~15% classified as Epi-CB or C2-Pure had a poorer outcome than those (p <0.0001; Fig. 6D). with tumors with a moderate 14q32-gene signature or classified Next, we assessed the in vivo anti-tumor effects of CHKA in- as Epi-CA or C1/C2B (Fig. S7). On the contrary, CTNNB1, VIM and hibition using MN58b in a PDX established from a high-risk HB BLCAP editing were not associated with any parameter of poor whose CHKA mRNA and protein levels are representative of in- outcome. termediate and high-risk tumors (Fig. S10). Interestingly, CHKA Next, based on the presence of the novel biomarkers of poor inhibition fully abrogated tumor growth throughout treatment prognosis, we defined the first molecular risk stratification of HB compared with the control arm (vehicle) (p = 0.028, Fig. 6E). The (MRS-HB) (Fig. 5A). The low-risk category (MRS-1) included tu- IHC study revealed that MN58b-treated tumors showed a signif- mors without any biomarker of poor prognosis (i.e., a moderate icantly lower proliferation rate, as determined by CCND1 (cyclin 14q32-gene signature and Epi-CA) and was enriched for wild- D1) and Ki67, as well as more commonly reverting from the tumor type CTNNB1 tumors (p = 0.012). Tumors in the intermediate- progenitor-like phenotype than vehicle-treated PDXs (Fig. S11). In risk category (MRS-2) were defined by having only one addition, MN58b-treated tumors showed a significant increase of biomarker (i.e., a strong 14q32-gene signature or Epi-CB), necrotic areas and a trend of having higher levels of the active form whereas those in the high-risk category (MRS-3) had two poor of caspase-3 than tumors in the control arm (Fig. S12). prognostic biomarkers (i.e., a strong 14q32-gene signature and Epi-CB) and were enriched for NFE2L2 mutations (p = 0.005). Discussion Kaplan-Meier survival curves showed that the 3-year EFS was Through comprehensive molecular profiling, we herein identi- 91%, 82%, and 52% for patients with MRS-1, MRS-2 and MRS-3 fied an unprecedented widespread epigenomic footprint of HB tumors, respectively (p <0.0001, Fig. 5B). To identify the most that includes RNA editing dysregulation, overexpression of aggressive tumors, we integrated the 16+VIM-gene signature to mainly non-coding genes in the oncogenic 14q32 DLK1-DIO3 the MRS and subdivided the high-risk category (MRS-3) into locus, and two distinct epigenomic tumor profiles that associate MRS-3a (C1 and C2B) and MRS-3b (C2-Pure); we defined the with the transcriptomic C1/C2/C2B subtypes previously defined latter as a very high-risk category, which was associated with a by Cairo-Armengol et al.9 and Hooks et al.10 The integration of 3-year EFS of only 37% (p <0.0001, Fig. 5B). Importantly, the these epigenetic hallmarks together with an updated 16-gene multivariate analysis indicated that this novel risk stratification signature allowed us to develop the first molecular risk stratifi- based on molecular parameters is an independent prognostic cation of HB, which improves on current clinical patient risk factor of the current clinical CHIC hepatoblastoma stratification6 classification, and to identify CHKA as a potential therapeutic (Fig. 5C,D). Accordingly, the combination of clinical and molec- target for patients with HB. ular staging systems (Fig. 5E) resulted in improved performance Our study revealed a genome-wide dysregulation of RNA at discriminating low- and high-risk patients (p <0.0001, Fig. 5F). editing in HB for the first time. RNA editing is a widespread Journal of Hepatology 2020 vol. 73 j 328–341 333 Research Article Hepatic and Biliary Cancer A Chr 14 2 p13 p1 14p11.2 q11.2 14q12 q21.1 q24.3 31.1 31.3 q32.2 chr14: 101,193,164-101,539,274 Ref. seq curated Gene symbol DLK1 MEG3 RTL1 MEG8 SNORD113’s SNORD114’s miR’s MEG9 miR’s Mir370 14q32.2-q32.31 10 5 gene expression Log2FC T/NT 0 SNORD114-26 DLK1 SNORD114-27 MEG3 -5 MIR2392 SNORD114-28 SNORD114-29 MIR770 SNORD114-30 MIR493 SNORD114-31 -10 MIR337 MIR379 MIR665 MIR441 RTL1 B C MIR431 MIR433 MIR299 MIR380 0.7 MIR1197 MIR127 MIR323A Enrichment score 0.6 MIR432 MIR758 NT26 NT15 NT16 NT19 NT76 NT30 NT18 NT78 NT17 NT13 NT14 NT12 MIR136 NT5 NT1 NT4 NT9 NT2 NT7 R27 T16 T15 T45 T10 T17 T10 T18 T30 T26 T29 T31 T33 T23 T21 T22 T28 T24 T12 T19 T32 T13 T25 T14 T20 T11 MIR329-1 T6 T5 T7 T1 T5 T3 0.5 MEG8 MIR329-2 0.4 MIR370 (ES) MIR494 0.3 FWER = 0.008 SNORD113-1 MIR1193 FDR = 0.017 SNORD113-2 0.2 MIR546 NES = 1.70 SNORD113-3 MIR495 0.1 SNORD113-4 MIR376C 0.0 SNORD113-5 MIR376A2 SNORD113-6 MIR654 SNORD113-7 MIR376B SNORD113-8 MIR376A1 SNORD113-9 MIR300 Rank list metric (Signal2Noise) SNORD114-1 MIR1185-1 1.0 T NT SNORD114-2 MIR1185-2 0.0 SNORD114-3 MIR381 SNORD114-4 -1.0 MIR487B SNORD114-5 MIR539 -2.0 SNORD114-6 MIR889 SNORD114-7 MIR544A SNORD114-8 MIR655 0 5,000 10,000 15,000 20,000 25,000 SNORD114-9 MIR487A Rank in ordered dataset SNORD114-10 MIR382 SNORD114-11 -3.0 -2.6 -2.3 -1.9 -1.6 -1.2 -0.9 -0.5 -0.2 0.2 0.5 0.9 1.2 1.6 1.9 2.3 2.6 3.0 MIR134 SNORD114-12 Enrichment profile Hits Ranking metric scores MIR668 Sample type: Non-tumor Tumor SNORD114-13 MIR485 SNORD114-14 14q32-gene signature: Moderate Strong MIR323B SNORD114-15 Wnt/β-catenin pathway mutations: No Yes MIR154 SNORD114-16 MIR496 SNORD114-18 MIR377 SNORD114-19 MIR541 SNORD114-20 MIR409 SNORD114-21 MIR412 SNORD114-22 MIR369 SNORD114-23 D E SNORD114-24 SNORD114-25 MIR410 MIR656 Non-tumor MEG9 Sample type Tumor with moderate 14q32-gene signature <0.0001 Chr 14 Tumor with strong 14q32-gene signature position 14q32 mean expression <0.0001 101,192,913 Fetal Liver 0 0 0 1 0 1 0 1 0 1 0 1 0 1 0 1 1 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 0 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 1 1 0 1 0 1 0 0 1 1 0 1 1 1 0 1 1 1 0 0 0 1 0 1 0 1 0 1 0 1 0 1 0 1 cg10528576 cg09459542 0.85 <0.0001 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 1 1 1 1 1 1 1 1 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 cg16562127 cg18603483 cg02412314 0.0279 cg12197579 20 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 0 1 0 0 0 0 1 0 1 0 1 1 1 1 1 1 1 cg13551098 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 cg02313552 Mean 14q32 β-value expression (log2) 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 0 1 1 cg02412314 0.0038 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 0 1 1 cg04828958 0.80 14q32 mean 15 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 0 0 1 0 0 0 1 0 1 0 0 0 0 1 0 0 1 1 0 1 0 0 1 1 0 0 1 1 1 1 cg13484881 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 0 1 1 1 1 1 1 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 cg08051604 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 0 1 1 1 1 1 1 1 1 cg04997124 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 0 1 1 1 1 1 1 1 1 cg23032838 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 cg05218696 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 cg27370351 cg04010923 10 1 0 1 0 1 0 1 0 1 1 1 0 1 1 1 0 1 1 1 1 1 1 1 0 1 0 1 0 1 1 1 0 1 1 1 0 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 0 1 0 1 0 1 0 1 1 1 0 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 cg00910338 cg08821973 cg09848749 0.75 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 0 1 1 0 1 1 0 1 0 0 1 1 0 1 1 1 1 1 0 1 1 1 1 1 0 1 1 0 1 0 0 1 1 0 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 0 1 1 0 0 0 1 1 1 0 0 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 0 1 0 1 1 0 1 cg05899471 cg18910298 5 r = -0.63 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 0 1 0 1 1 1 1 1 1 1 1 1 0 0 1 1 1 1 cg20425627 cg21178618 p <0.0001 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 cg19603100 cg00175487 0 0.70 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 0 1 0 1 1 1 0 0 1 1 1 1 1 1 1 cg25562958 0.2 0.4 0.6 0.8 1.0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 0 1 1 0 0 0 1 1 0 1 1 1 1 1 1 1 1 cg09243104 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 0 1 0 0 0 1 1 1 0 0 0 1 0 1 1 1 1 cg21146428 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 cg10137253 cg03019112 β-value 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 cg23029159 101,528,354 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 0 1 1 0 1 1 0 1 1 1 1 1 1 1 cg11134801 β values + - 0.65 Fig. 3. Overexpression of genes of the 14q32 DLK1-DIO3 locus is a hallmark of hepatoblastoma. (A) Top, scheme of chromosome 14 in which the overex- pressed 14q32 region in HB is marked with a red square. Middle, detail of the overexpressed 14q32 genes of the DLK1-DIO3 locus. Bottom, the X-axis indicates the gene chromosomal localization from 14q32.13 to 14qter. Y-axis is the mean tumor (T)/non-tumor (NT) FC of the genes of the 14q32 region present in the HTA array in normalized arbitrary units. The box and whiskers represent the 25th to 75th percentiles and ± min to max values, respectively, of log2 gene expression of 32 tumors. (B) GSEA enrichment plot using the 96 genes of the HTA array out of the total of 101 genes presenting the 14q32 region (chr14:101193164-101539274) based on the Genome Browser Database hg19. (C) Unsupervised clustering and heatmap of the 133 HTA array probesets (right) at the overexpressed 14q32 chromosomal region of 32 T and 18 NT samples. (D) Left, heatmap of the methylation degree of the 32 significant CpGs (FDR <0.05) localized in the overexpressed 14q32 region. The color range in the heatmap indicates the methylation degree of each CpG. Right, correlation between mean gene expression of 14q32 genes (133 probesets) and the methylation (b-value) of the most significant 14q32 CpG, which is localized in a MEG3 intronic region (Spearman test). (E) Mean 334 Journal of Hepatology 2020 vol. 73 j 328–341 epigenetic mechanism that confers specific nucleotide changes is linked to high expression of liver progenitor cell markers and on RNA transcripts without altering the sequence of genomic Wnt/b-catenin targets and that it influences the survival of pa- DNA; thereby contributing to transcriptomic diversity in normal tients with HB. Collectively, our study pinpoints the over- but also cancer cells.39 The functional impact of RNA editing on expression of 14q32 genes of the DLK1-DIO3 locus as a novel cell biology ranges from protein recoding to alterations in oncogenic hallmark of HB. alternative splicing, miRNA specificity and RNA stability.40 Dysregulation of DNA methylation might be considered a third Herein, in HBs, we identified global hypoediting, specific epigenetic hallmark of HB. DNA methylation plays an important hyperediting of BLCAP, a highly conserved gene with potential role in cell differentiation and cancer, influencing the regulation tumor suppressor functions,22,41,42 and also an imbalance in the of gene expression networks.52,53 The genome-wide DNA hypo- expression of ADAR enzymes; the main regulators of RNA edit- methylation that we found in HB is consistent with the findings of ing.23 Our findings agree with the observed dysregulation of RNA previous reports.15–17 As a novelty, we discovered 2 distinct editing across different cancer types with over- and under- epigenetic profiles in HB based on the degree of DNA hypo- editing patterns relative to non-tumor samples.39,40,43 Similar methylation and CpG island hypermethylation. We also associ- to our data, global RNA editing dysregulation in cancer has also ated these epigenomic traits with the previously defined C1/C2/ been associated with the hyperediting of key genes.40 In that C2B molecular subclasses.9,10 Our results demonstrate how the regard, increased hyperediting of BLCAP has already been re- interplay between the epigenome and the transcriptome de- ported in HCC44 and to a minor degree in brain, cervical, oral termines distinct tumor molecular entities. Moreover, we inves- cavity and lung tumors.45,46 Previous investigations revealed that tigated the impact of CpG island hypermethylation, an additional BLCAP editing could affect the functions of several binding pro- level of epigenetic dysregulation in Epi-CB tumors belonging teins such as RB47 or STAT346; thereby, influencing proliferation mainly to the C2 subtype, and identified novel tumor suppressor and apoptotic signaling pathways. In HCC, experiments per- candidates whose expression was strongly repressed in aggres- formed in SMMC_7721 and Focus liver cancer cell lines revealed sive HB, which could be explored in further functional studies. that nt 5 editing confers a growth advantage, modulating the In an attempt to translate our findings into the clinical setting, activation of AKT/mTOR signaling.44 Overall, our study reveals an we propose the first molecular risk stratification system, called unexplored field related to RNA editing in HB. Future studies will MRS, for HB. This system is based on the presence of the 2 novel need to clarify the functional effects of BLCAP editing in these prognostic biomarkers identified in the current study (i.e. 14q32- tumors. gene signature and Epi-CA/B). Of note, the prognostic impact of The second epigenetic alteration we observed pertains to the the MRS is improved by incorporating the updated 16-gene pronounced overexpression of mainly non-coding genes local- signature described here, which includes VIM expression, to ized in a small 14q32 region of the DLK1-DIO3 locus in almost all distinguish the recently reported C2B subclass.10 The benefit of the HBs examined and that is highly correlated with their degree the MRS compared to our previously published 16-gene signa- of methylation. In addition to DLK1, a well-known hepatoblast ture is that it is able to better differentiate patients according to marker overexpressed in HB,9,26 this locus is characterized by a their prognosis. This could be explained by the fact that MRS cluster of imprinted genes whose altered dosage is associated integrates both epigenetic and transcriptomic classifiers, with developmental defects and liver oncogenesis.24,29,48 Inter- providing a better representation of the molecular complexity of estingly, the expression of the DLK1-DIO3 locus has been pro- HB. Moreover, like the clinical CHIC hepatoblastoma stratifica- posed as a marker of induced pluripotent stem cells, thereby tion,6 the integration of multiple molecular prognostic factors supporting its role in early development.49 In agreement with may have an additive effect in terms of risk prediction. By these previous studies, our data suggest a fine-tuned regulation combining the clinical CHIC and the molecular MRS systems, we of the 14q32 region during liver development, since its tran- have been able to further improve patient risk prediction. In this scripts were highly expressed in fetal livers but strongly regard, our findings highlight the importance of incorporating repressed in postnatal ones. In that regard, we found that HBs molecular factors into the clinical setting, thereby facilitating present an aberrant expression of 14q32 genes, an observation future precision medicine. The main limitations of the current that supports the notion that these genes are involved in hepatic study lie in the use of a retrospective cohort of patients treated tumorigenesis. The oncogenic role of 14q32 genes in the liver with different chemotherapeutic protocols and the study of post- was initially identified through research into the mechanisms chemotherapy specimens. Thus, the implementation of our involved in the spontaneous development of HCC in mice treated findings into the clinical setting requires a further validation in with adeno-associated viruses (AAVs).29 Moreover, the gene- diagnostic biopsies from homogeneously treated patients and targeting frequency of this locus by AAVs was shown to be suf- probably the definition of a new algorithm integrating clinical ficient to initiate multiple foci of HCC in mice characterized by and molecular parameters; the prospective cohort of patients Dlk1-Glt2 overexpression linked to CpG hypomethylation.50 In enrolled in the ongoing Paediatric Hepatic International Tumour agreement with these experimental data, 2 independent studies Trial (PHITT, NCT03017326) provides a unique opportunity to reported the overexpression of 14q32 genes in a subset of 6–19% conduct such validation. of patients with HCC and poor prognosis.48,51 In line with these Finally, we identified CHKA as a novel potential therapeutic studies on HCC, we found that the overexpression of 14q32 genes target for HB patients. CHKA, a key gene for membrane = methylation levels (b-value) of all 568 CpGs localized in the 14q32 region in the 3 different sample types (NT, n = 19; T with moderate and strong 14q32-gene signature, n = 12 and n = 15, respectively; fetal liver samples, n = 5). p values were calculated using the ANOVA test with the Tukey post-test. FC, fold change; FDR, false discovery rate; GSEA, gene set enrichment analysis; HB, hepatoblastoma; HTA, human transcriptome array. Journal of Hepatology 2020 vol. 73 j 328–341 335 Research Article Hepatic and Biliary Cancer A B 400 600 300 400 PC3 10% 200 600 200 0 400 -200 200 % 0 100 12 -400 -200 -400 2 -600 PC -600 -800 0 00 -6 0 00 00 0 00 0 0 0 0 20 40 60 ,0 -8 -2 -4 NT13 NT15 NT30 NT10 NT12 NT6 NT14 NT1 NT28 NT26 NT16 NT23 NT2 NT4 NT17 NT5 NT11 NT33 NT7 T7 T8 T28 T9 T18 T17 T4 T6 T15 T10 T2 T16 T14 T5 T20 T21 T13 T1 T12 T31 T22 T25 T11 T33 T30 T23 T32 T26 -1 PC1 20% 0.65 Epi-CA HB Epi-CB HB p value Sample type <0.0001# 0.60 16-gene classification 0.6882 Mean β-values HTA 0.0003# Transcriptomic co-clusters RNA-seq 0.0005# 0.55 Wnt/β-catenin pathway mutations 0.0391 14q32-gene signature 0.0010 Sample type: NT T 0.50 16-gene classification: C1 C2 p <0.0001 Transcriptomic co-clusters: CC1 CC2 CC3 Wnt/β-catenin pathway mutations: No Yes 0.45 14q32-gene signature: Moderate Strong NT T C Open sea N shores N shelves CpG islands S shelves S shores **** **** *** 0.80 0.55 0.80 0.26 0.85 **** 0.55 *** *** * ** 0.80 *** *** Mean β-values **** 0.50 0.75 **** 0.24 0.50 0.70 0.75 **** 0.70 0.22 0.45 0.70 0.45 0.60 0.65 0.20 0.65 0.40 0.40 0.60 0.18 0.60 0.50 0.35 0.55 0.16 0.55 0.35 NT Epi-CA HB subclass Epi-CB HB subclass D KLF6 ITGB3 NFIC TRANK1 TSPYL5 30,000 1,500 4,000 2,000 500 Gene expression 3,000 400 1,500 20,000 1,000 300 2,000 1,000 200 10,000 500 1,000 500 100 0 0 0 0 0 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 0 0.2 0.4 0.6 0.8 1 cg06048750 β-value cg12683944 β-value cg03536983 β-value cg07696485 β-value cg04917181 β-value NT Epi-CA HB subclass Epi-CB HB subclass E 0.65 NT Epi-CA HB subclass FL Epi-CB HB subclass 0.60 Mean β-values 0.55 0.50 r2 = 0.5749 0.45 -9 0 10 20 30 40 50 60 70 190 Birth Time (months) Fig. 4. Identification of two epigenetic hepatoblastoma subtypes resembling the transcriptomic co-clusters. (A) Top, principal component analysis of tumor (T, n = 28) and non-tumor (NT, n = 19) samples using normalized methylation data of the 685,375 CpGs of the 850K-array obtained after filtering those probes containing single nucleotide polymorphisms. Bottom, mean global methylation levels (b-value) of the same samples (unpaired t test). The box and whiskers represent the 25th to 75th percentiles and ± min to max values, respectively. (B) Representative unsupervised clustering (Euclidean Ward method) of same CpGs used for the principal component analysis. p values indicate the associations between the listed molecular features and the 2 tumor epigenetic clusters, Epi-CA and Epi-CB. p values were calculated using Fisher's exact and #Chi-square tests. (C) Methylation levels of the epigenetic substructures in NT (n = 19), Epi-CA (n = 13) and Epi-CB (n = 13) tumors. The box and whiskers represent the 25th to 75th percentiles and ± min to max values, respectively. p values were calculated using the ANOVA test with the Tukey post-test (*p <0.05; **p < −0.001; ***p <0.0001). (D) Representation of the 5 most strongly repressed genes by CpG island 336 Journal of Hepatology 2020 vol. 73 j 328–341 Very high-risk A MRS-3 MRS-1 MRS-2 Low-risk Intermediate-risk MRS-3a MRS-3b High-risk Very high-risk T104 T101 T108 T107 T103 T100 T106 T102 T105 T109 T110 T36 T40 T34 T51 T86 T85 T15 T18 T52 T59 T16 T17 T10 T75 T70 T77 T46 T53 T44 T43 T41 T35 T76 T93 T79 T57 T80 T90 T55 T68 T74 T47 T62 T54 T49 T24 T19 T78 T96 T71 T66 T67 T39 T14 T13 T28 T58 T20 T61 T60 T83 T42 T63 T48 T84 T37 T38 T89 T87 T26 T81 T72 T95 T98 T88 T33 T25 T12 T45 T29 T99 T82 T50 T92 T27 T22 T69 T94 T97 T91 T32 T23 T30 T31 T21 T11 Classifier T7 T2 T4 T6 T9 T8 T1 T5 T3 Strong 14q32 Epi-CB C2-Pure p value Wnt/β-catenin activation 0.0124 NFE2L2 mutation 0.0046 BLCAP nt5 hyperediting 0.4908 Histology: mixed 0.2281 Main ep. comp: non-fetal 0.0748 CHIC-HS 0.0242 Age >8 yr 0 0.0192 AFP >1,000 ng/ml 1 0.6217 PRETEXT IV 0.5139 Metastasis 0.1813 Vascular invasion 0.0033 Multifocality 0.0388 Poor outcome <0.0001 NA No Yes Clinical risk stratification (CHIC-HS): Very low/low risk Intermediate risk High risk B Molecular risk stratification (MRS-HB) 1.0 MRS-1 (n = 35) 1.0 91% 91% MRS-1 (n = 35) 82% MRS-2 (n = 40) 82% MRS-2 (n = 40) 0.8 0.8 C 67% 0.6 0.6 95% CI MRS-3a (n = 15) EFS MRS-3 (n = 31) Lower Higher p-value EFS 52% HR 0.4 0.4 CHIC-HS (L/VL-I-H) 3.17 1.83 5.48 <0.0001 37% MRS-3b (n = 16) MRS-HB (1-2-3) 2.21 1.26 3.88 0.006 0.2 0.2 Log rank = 1.67x10-4 Log rank = 2.10x10-5 0.0 0.0 0 20 40 60 80 100 0 20 40 60 80 100 Time (months) Time (months) D CHIC-HS E CHIC-HS F CHIC-HS + MRS-HB 1.0 95% VLR/LR (n = 60) CHIC-HS MRS-HB 1.0 + MRS-HB 100% LR (n = 22) 0.8 73% IR (n = 11) 0.8 MRS-1 LR 81% IR (n = 68) VLR/LR 0.6 0.6 MRS-2/3 EFS EFS 46% HR (n = 35) IR IR 0.4 0.4 MRS-1/2 HR HR (n = 16) 0.2 MRS-3 HR 0.2 25% Log rank =3x10-6 Log rank =1.03x10-7 0.0 0.0 0 20 40 60 80 100 0 20 40 60 80 100 Time (months) Time (months) Fig. 5. Molecular Risk Stratification of Hepatoblastoma (MRS-HB). (A) Top, representation of the MRS-HB categories for the 106 patients from which all biomarkers could be assessed. Patient stratification was based on transcriptomic and epigenetic classifiers (right): 14q32- and 16+VIM-gene signatures and epigenomic Epi-CA/Epi-CB classification. Bottom, main genomic, pathologic and clinical features and their association with the MRS (right, Chi-Square test). (B) EFS Kaplan-Meier plots of the same patients stratified according to the MRS-HB into 3 (left) or 4 (right) categories. Vertical line indicates 3-year EFS probability (in %). (C) Multivariate Cox regression analysis comparing MRS with clinical CHIC-HS Stratification.6 (D) EFS Kaplan-Meier plots of the same patients with HB classified following the CHIC-HS stratification (VL/L, Very Low and Low; L, Low; I, Intermediate; H, High risk). (E) Scheme used to combine CHIC-HS and MRS-HB classifications. (F) EFS Kaplan-Meier plots of the 106 patients with HB classified by combining clinical and molecular risk stratification systems. CHIC-HS, Children's Hepatic tumors International Collaboration-Hepatoblastoma Stratification; Epi-CA/B, Epigenetic-Cluster; EFS, event-free survival; HB, hepato- blastoma; HR, hazard ratio; MRS, Molecular Risk Stratification. = hypermethylation in Epi-CB tumors (RNA-seq/HTA criteria: FC <−2 and FDR <0.001; 850K-array criteria: b >0.2 and FDR <0.0001). The X-axis indicates methylation levels of the most hypermethylated CpGs islands for each gene and the Y-axis the linear gene expression levels (HTA). The grey shadow indicated the low CpG island methylation levels associated to high gene expression. (E) Global methylation levels (Y-axis) of 19 NT (white dots) and 5 fetal liver samples (FL; green dots) at different gestational and postnatal ages (X-axis). The light and dark purple shadows indicate the 25th and 75th percentiles of the methylation levels of Epi-CA and Epi-CB tumors, respectively, and their extrapolation over time. Note the plateau of methylation levels at ~ 12–24 months of age. Red, exponential curve. Epi- CA/B, Epigenetic-Cluster A/B; FC, fold change; FDR, false discovery rate; HTA, human transcriptome array; ITGB3, integrin subunit beta 3; KLF6, Kruppel like factor 6; NFIC, nuclear factor I C; TRANK1, tetratricopeptide repeat and ankyrin repeat containing 1; TSPYL5, TSPY like 5. Journal of Hepatology 2020 vol. 73 j 328–341 337 Research Article Hepatic and Biliary Cancer 0.0012 A HTA RNA-Seq B <0.0001 0.0020 <0.0001 MRS-1 MRS-2 MRS-3 4,000 100 <0.0001 80 0.0256 % stained area (CHKA) 0.0001 CHKA expression CHKA expression 80 3,000 60 60 2,000 40 40 1,000 20 20 0 0 0 -1 -2 -3 -1 -2 -3 -1 -2 -3 RS RS RS RS RS RS RS RS RS M M M M M M M M M C HepG2 Huh6 150 150 150 150 % viability % viability % viability 100 100 100 100 % viability * ** *** ** *** *** **** **** **** 50 50 50 **** 50 **** **** **** **** 0 0 0 0 0 2 4 6 8 0 2 4 6 8 0 2 4 6 8 0 2 4 6 8 [MN58b] (μM) [TCT-717] (μM) [MN58b] (μM) [TCT-717] (μM) 150 150 150 150 % colony formation rate % colony formation rate % colony formation rate % colony formation rate 100 100 100 100 50 50 50 50 **** **** **** **** **** **** **** 0 0 **** **** 0 **** 0 **** **** 0 2 4 6 0 2 4 6 0 2 4 6 0 2 4 6 [MN58b] (μM) [TCT-717] (μM) [MN58b] (μM) [TCT-717] (μM) MN58b MN58b TCD-717 TCD-717 0 μM 2 μM 4 μM 6 μM 0 μM 2 μM 4 μM 6 μM D Sham siRNA E 4,000 p = 0.028 CHKA siRNA 3,500 Tumor volume (mm3) 1.5 300 3,000 *** 2,500 CHKA mRNA % viability 1.0 200 2,000 Control * 1,500 MN58b 0.5 100 1,000 *** 500 0 0 0 24 h 48 h 72 h 2 5 9 12 16 19 23 26 Time after treatment (days) Fig. 6. The effect of inhibition of CHKA in hepatoblastoma. (A) CHKA gene expression assessed by HTA and RNA-seq in the different risk molecular risk categories (MRS-1, n = 11; MRS-2, n = 8; MRS-3, n = 13)╪. (B) Left, Images of CHKA immunostaining for representative tumors according to MRS. Right, quantification of stained areas of CHKA immunohistochemistry in formalin-fixed paraffin-embedded tumor samples (MRS-1, n = 9; MRS-2, n = 6; MRS-3, n = 5)╪. (C) Cell viability assay (MTT)╪ (top) and colony-formation assay╪ (middle) of HepG2 (left) and Huh6 (right) cells treated by at different concentrations of CHKA inhibitors (MN58b and TCD-717) for 48 h. 0 lM of TCD-717 was prepared with DMSO at the same concentration of 8 lM of TCD-717. Representative colony-formation assay images of the different conditions (bottom). Data given are from a minimum of 3 independent experiments. p values were calculated vs. control. (D) CHKA gene expressionY assessed by RT-qPCR of control (Sham siRNA, white) and CHKA knock-down (CHKA-siRNA, black) in Huh6 cells. Cell viability assay (MTT) assay╪ of the same cells at different time points of culture. Data given are from four independent experiments. p values were calculated vs. control. (E) Tumor growth curve of a PDX established from a high-risk tumor (MRS-3, case T50) following intraperitoneal injection of vehicle (phosphate buffered saline, control; n = 5) or 3 mg/kg/day of CHKA inhibitor (MN58b; n = 6). p values was calculated using the two-tailed t test with Welch's correction 338 Journal of Hepatology 2020 vol. 73 j 328–341 biosynthesis, is overexpressed in different neoplasms such as (HUNTER, C9380/A26813), National Cancer Institute, Tisch Can- breast, lung, prostate, and HCC.54–57 The complete abrogation of cer Institute (P30-CA196521), Samuel Waxman Cancer Research tumor growth that we observed in vitro and in vivo using two HB Foundation, Spanish National Health Institute (SAF2016-76390) cell lines and a PDX model is achieved by an inhibition of pro- and AGAUR (SGR-1358). CA and MRS were supported by Ramón y liferation and induction of cell death. Of note, CHKA has been Cajal (RYC-2010-07249) and Miguel Servet (CPII14/00021) pro- proposed as a therapeutic target in different tumor types57,58 and grams of the Ministry of Science and Innovation of Spain and TCD-717 has already been evaluated in a phase I clinical trial ISCIII, respectively. CA, MRS and MS received funding from (NCT01215864) in advanced solid tumors. Our findings thereby CIBERehd (CB06/04/0033) and AGAUR (2017-SGR-490). IGTP is a support further attention to CHKA inhibition as a potential member of the CERCA network of institutes. The funders had no therapy for patients with HB. role in study design, data collection and analysis, decision to The similarities observed between HB (i.e. BLCAP hyper- publish, or preparation of the manuscript. editing, 14q32 locus overexpression, CHKA overexpression) and HCC also point to common underlying mechanisms between the Conflict of interest main hepatic tumors in childhood and adulthood, shedding light Prof. Josep M. Llovet is receiving research support from Bayer on overlapping molecular mechanisms and common therapeutic HealthCare Pharmaceuticals, Eisai Inc, Bristol-Myers Squibb and approaches. Ipsen, and consulting fees from Bayer HealthCare Pharmaceuti- In summary, our data provide novel epigenetic insights into cals, Bristol-Myers Squibb, Eisai Inc, Celsion Corporation, Eli Lilly, HB and establish the rationale to advance towards precision Exelixis, Merck, Ipsen, Glycotest, Navigant, Leerink Swann LLC, medicine by identifying new biology-driven therapies and Midatech Ltd, Fortress Biotech, Sprink Pharmaceuticals and incorporating molecular data into patient stratification. Nucleix and CANFITE. CA has a research contract with CHIOME Biosciences Inc. The other authors report no conflicts of interest Abbreviations in this work. AAV, adeno-associated viruses; BLCAP, bladder cancer associated Please refer to the accompanying ICMJE disclosure forms for protein; CHIC, Children's Hepatic tumors International Collabo- further details. ration; CHIC-HS, CHIC-Hepatoblastoma Stratification; CHKA, choline kinase alpha; ddPCR, droplet digital PCR; EFS, event-free Authors' contributions survival; Epi-CA/B, Epigenetic-Cluster A/B; FA, fractional abun- JCR designed, performed research, analyzed all data, and wrote dance; FC, fold change; FDR, false discovery rate; GSEA, gene set the results; LT designed, performed research and reviewed the enrichment analysis; HB, hepatoblastoma; HCC, hepatocellular manuscript; MSC designed and performed research; JAF, LR, MV, carcinoma; HCN-NOS, hepatocellular neoplasm not otherwise NV and SR performed research; LN, MA and MM performed specified; HR, hazard ratio; HTA, human transcriptome array; bioinformatics analyses of transcriptomic, genomic and epige- ITGB3, integrin subunit beta 3; KLF6, Kruppel like factor 6; MEC, netic data; BL and NA developed and performed RNA variant and main epithelial component; MRS, Molecular Risk Stratification; editing analysis and BL, DS and NA performed expression ana- NFIC, nuclear factor I C; NF1, neurofibromin 1; PDX, patient- lyses of RNA-seq data; OK, performed research; AV performed derived xenograft; PRETEXT, PRETreatment EXTent of disease; PDX experiment; SC, RK, MAB, BM and PC provided clinical RNA-seq, RNA sequencing; RT-qPCR, reverse transcription samples, the associated relevant clinical and pathological infor- quantitative PCR; siRNA, small interfering RNA; TRANK1, tetra- mation and critically review the manuscript; MDS; performed tricopeptide repeat and ankyrin repeat containing 1; TSPYL5, research and critically review the manuscript; AS, CS, LG, MGa, TSPY like 5. MGo, MEM, SB, GR, MLS, GG, YM, NGA, JJU, BLI, BT, MF, RLA, JAS, CP, VB, GL, CB, CG, FH, RP, HM, AC and MS provided clinical samples as well as the associated relevant clinical and patho- Financial support logical information; JB, DP, LS and MJ designed and performed This article was possible thanks to the inputs from the Instituto research; MRS, designed, analyzed data, and reviewed the de Salud Carlos III, ISCIII (PI09/00751, PI10/02082, PI13/02340). manuscript; JML designed research and supervised the study. CA The project has received funding from the European Union's designed research, analyzed data, supervised the whole study Horizon 2020 research and innovation programme under grant and wrote the manuscript. agreement No 668596 (ChiLTERN) and grant agreement No 826121 (iPC). JCR is supported by the Catalan Agency for Man- agement of University and Research Grants (AGAUR, 2019 Acknowledgements FI_B01024). LT is supported by an Accelerator Award (CRUCK, We thank patients and their families for allowing the collection AECC, AIRC) (HUNTER, C9380/A26813). DS is supported by the of biological samples. We also would like to acknowledge the Gilead Research Scholar in Liver Disease. JML is supported by the support of the Spanish Group for the Study Of Childhood Liver European Union's Horizon 2020 research and innovation pro- Cancer of the Spanish Society of Paediatric Haematology and gramme (HEPCAR, 667273-2), Institució Catalana de Recerca i Oncology (SEHOP) and the clinical staff of Spanish hospitals that Estudis Avançats (ICREA), U.S. Department of Defense provided us with samples and clinical data (Maciej Murawski, (CA150272P3), an Accelerator Award (CRUCK, AECC, AIRC) Marta Villa, Sandra Pisa, Montserrat Torrent, Pablo Trujillo, Ester = at 26 days of treatment. Other statistical analysis performed were ANOVA test with the Tukey post-test╪ and unpaired t testY (*p <0.05; **p < −0.001; ***p <0.0001). CHKA, choline kinase alpha; HTA, Human Transcriptome Array; MRS, Molecular Risk Stratification; PDX, patient-derived xenograft; RNA-seq, RNA sequencing; RT- qPCR, reverse transcription quantitative PCR; siRNA, small interfering RNA. Journal of Hepatology 2020 vol. 73 j 328–341 339 Research Article Hepatic and Biliary Cancer Anton, Jose Andres Molino, Hermenegildo González, Ma Esther [16] Cui X, Liu B, Zheng S, Dong K, Dong R. Genome-wide analysis of Llinares, José Luis Fuster and Isabel Ojanguren). We thank the DNA methylation in hepatoblastoma tissues. Oncol Lett 2016;12:1529– 1534. IGTP Biobank and Genomic Units and the HUVR-IBiS Biobank [17] Honda S, Minato M, Suzuki H, Fujiyoshi M, Miyagi H, Haruta M, et al. (Andalusian Public Health System Biobank and ISCIII-Red de Clinical prognostic value of DNA methylation in hepatoblastoma: four Biobancos PT17/0015/0041) for the assessment and technical novel tumor suppressor candidates. Cancer Sci 2016;107:812–819. support provided. We also thank Dr. Eduard Cabré for his advice [18] Edgar R, Domrachev M, Lash AE. Gene Expression Omnibus: NCBI gene expression and hybridization array data repository. Nucleic Acids Res on statistical analysis. We acknowledge the personnel at the 2002;30:207–210. Genome Analysis Platform of CIC bioGUNE and CIBERehd for [19] Tomlinson GE, Douglass EC, Pollock BH, Finegold MJ, Schneider NR. Cy- their help on experimental design and for processing the togenetic evaluation of a large series of hepatoblastomas: numerical ab- methylation Beadchips. We acknowledge Koji Nakamura normalities with recurring aberrations involving 1q12-q21. Genes (CHIOME Biosciences Inc) for kindly providing DLK1 antibody. Chromosomes Cancer 2005;44:177–184. [20] Albrecht S, von Schweinitz D, Waha A, Kraus JA, von Detailing A, Pietsch T. We thank Ugne Balaseviciute for her assistance in experimental Loss of maternal Alleles on chromosome arm 11p in hepatoblastoma. procedures. Cancer Res 1994;54:5041–5044. [21] López-Terrada D, Alaggio R, de Dávila MT, Czauderna P, Hiyama E, Katzenstein H, et al. Towards an international pediatric liver tumor Supplementary data consensus classification: proceedings of the Los Angeles COG liver tumors symposium. Mod Pathol 2014;27:472–491. Supplementary data to this article can be found online at https:// [22] Levanon EY, Hallegger M, Kinar Y, Shemesh R, Djinovic-Carugo K, doi.org/10.1016/j.jhep.2020.03.025. Rechavi G, et al. Evolutionarily conserved human targets of adenosine to inosine RNA editing. Nucleic Acids Res 2005;33:1162–1168. References [23] Nishikura K. Functions and regulation of RNA editing by ADAR Author names in bold designate shared co-first authorship deaminases. Annu Rev Biochem 2010;79:321–349. [24] da Rocha ST, Edwards CA, Ito M, Ogata T, Ferguson-Smith AC. Genomic imprinting at the mammalian Dlk1-Dio3 domain. Trends Genet [1] Schnater JM, Köhler SE, Lamers WH, Von Schweinitz D, Aronson DC. 2008;24:306–316. Where do we stand with hepatoblastoma? A review. Cancer [25] Benetatos L, Hatzimichael E, Londin E, Vartholomatos G. The microRNAs 2003;98:668–678. within the DLK1-DIO3 genomic region: involvement in disease patho- [2] Linabery AM, Ross JA. Trends in childhood cancer incidence in the U.S. genesis. Cell Mol Life Sci 2013;70:795–814. (1992-2004). Cancer 2008;112:416–432. [26] López-Terrada D, Gunaratne PH, Adesina AM, Pulliam J, Hoang DM, [3] Kehm RD, Spector LG, Osypuk TL, Poynter JN, Vock DM. Do pregnancy Nguyen Y, et al. Histologic subtypes of hepatoblastoma are characterized characteristics contribute to rising childhood cancer incidence rates in the by differential canonical Wnt and Notch pathway activation in DLK+ United States? Pediatr Blood Cancer 2018;65:1–9. precursors. Hum Pathol 2009;40:783–794. [4] Aronson DC, Czauderna P, Maibach R, Perilongo G, Morland B. The [27] Van de Wetering M, Sancho E, Verweij C, De Lau W, Oving I, Hurlstone A, treatment of hepatoblastoma: its evolution and the current status as per et al. The b-catenin/TCF-4 complex imposes a crypt progenitor phenotype the SIOPEL trials. J Indian Assoc Pediatr Surg 2014;19:201–207. on colorectal cancer cells. Cell 2002;111:241–250. [5] Semeraro M, Branchereau S, Maibach R, Zsiros J, Casanova M, Brock P, [28] Zimmermann A. The emerging family of hepatoblastoma tumours: from et al. Relapses in hepatoblastoma patients: clinical characteristics and ontogenesis to oncogenesis. Eur J Cancer 2005;41:1503–1514. outcome-Experience of the International Childhood Liver Tumour Strat- [29] Donsante A, Miller DG, Li Y, Vogler C, Brunt EM, Russell DW, et al. AAV egy Group (SIOPEL). Eur J Cancer 2013;49:915–922. vector integration sites in Mouse hepatocellular carcinoma. Science [6] Meyers RL, Maibach R, Hiyama E, Häberle B, Krailo M, Rangaswami A, 2007;317:477. et al. Risk-stratified staging in paediatric hepatoblastoma: a unified [30] Wang P-R, Xu M, Toffanin S, Li Y, Llovet JM, Russell DW. Induction of analysis from the Children's Hepatic tumors International Collaboration. hepatocellular carcinoma by in vivo gene targeting. Proc Natl Acad Sci Lancet Oncol 2017;18:122–131. U S A 2012;109:11264–11269. [7] Armengol C, Cairo S, Fabre M, Buendia MA. Wnt signaling and hep- [31] Harada K, Toyooka S, Maitra A, Maruyama R, Toyooka KO, Timmons CF, atocarcinogenesis: the hepatoblastoma model. Int J Biochem Cell Biol et al. Aberrant promoter methylation and silencing of the RASSF1A gene 2011;43:265–270. in pediatric tumors and cell lines. Oncogene 2002;21:4345–4349. [8] Perugorria MJ, Olaizola P, Labiano I, Esparza-baquer A, Marzioni M, [32] Chow RKK, Sin ST-K, Liu M, Li Y, Chan THM, Song Y, et al. AKR7A3 sup- Marin JJG, et al. Wnt – b -catenin signalling in liver development, health presses tumorigenicity and chemoresistance in hepatocellular carcinoma and disease. Nat Rev Gastroenterol Hepatol 2019;16:121–136. through attenuation of ERK, c-Jun and NF-jB signaling pathways. Onco- [9] Cairo S, Armengol C, De Reyniès A, Wei Y, Thomas E, Renard CA, et al. target 2017;8:83469–83479. Hepatic stem-like phenotype and interplay of wnt/b-catenin and [33] Chen C, Wang L, Liao Q, Huang Y, Ye H, Chen F, et al. Hypermethylation of Myc signaling in aggressive childhood liver cancer. Cancer Cell EDNRB promoter contributes to the risk of colorectal cancer. Diagn Pathol 2008;14:471–484. 2013;8:909. [10] Hooks KB, Audoux J, Fazli H, Lesjean S, Ernault T, Senant N-D, et al. New [34] Deloria AJ, Höflmayer D, Kienzl P, Łopatecka J, Sampl S, Klimpfinger M, insights into diagnosis and therapeutic options for proliferative hepato- et al. Epithelial splicing regulatory protein 1 and 2 paralogues correlate blastoma. Hepatology 2018;68:89–102. with splice signatures and favorable outcome in human colorectal cancer. [11] Sumazin P, Chen Y, Treviño LR, Sarabia SF, Hampton OA, Patel K, et al. Oncotarget 2016;7:73800–73816. Genomic analysis of hepatoblastoma identifies distinct molecular and [35] Tessitore L, Sesca E, Vance DE. Inactivation of phosphatidylethanolamine prognostic subgroups. Hepatology 2017;65:104–121. N-methyltransferase-2 in aflatoxin-induced liver cancer and partial [12] Gröbner SN, Worst BC, Weischenfeldt J, Buchhalter I, Kleinheinz K, reversion of the neoplastic phenotype by PEMT transfection of hepatoma Rudneva VA, et al. The landscape of genomic alterations across childhood cells. Int J Cancer 2000;86:362–367. cancers. Nature 2018;555:321–327. [36] Zhang F, Sun H, Zhang S, Yang X, Zhang G, Su T. Overexpression of PER3 [13] Eichenmüller M, Trippel F, Kreuder M, Beck A, Schwarzmayr T, Häberle B, inhibits self-renewal capability and chemoresistance of colorectal cancer et al. The genomic landscape of hepatoblastoma and their progenies with stem-like cells via inhibition of notch and b-catenin signaling. Oncol Res HCC-like features. J Hepatol 2014;61:1312–1320. Featur Preclin Clin Cancer Ther 2017;25:709–719. [14] Buendia MA. Unravelling the genetics of hepatoblastoma: few mutations, [37] Aoyama C, Liao H, Ishidate K. Structure and function of choline kinase what else? J Hepatol 2014;61:1202–1204. isoforms in mammalian cells. Prog Lipid Res 2004;43:266–281. [15] Maschietto M, Rodrigues TC, Kashiwabara AY, de Araujo ÉSS, Aguiar TFM, [38] de la Cueva A, Ramírez de Molina A, Álvarez-Ayerza N, Ramos MA, da Costa CML, et al. DNA methylation landscape of hepatoblastomas Cebrián A, del Pulgar TG, et al. Combined 5-FU and ChoKa inhibitors reveals arrest at early stages of liver differentiation and cancer-related as a new alternative therapy of colorectal cancer: evidence in alterations. Oncotarget 2017;8:97871–97889. 340 Journal of Hepatology 2020 vol. 73 j 328–341 human tumor-derived cell lines and mouse xenografts. PLoS One [49] Benetatos L, Vartholomatos G, Stem ÁMÁ. DLK1-DIO3 imprinted cluster 2013;8:1–13. in induced pluripotency: landscape in the mist. Cell Mol Life Sci [39] Paz-Yaacov N, Bazak L, Buchumenski I, Porath HT, Danan-Gotthold M, 2014;71:4421–4430. Knisbacher BA, et al. Elevated RNA editing activity is a major contributor [50] Wang P, Xu M, Toffanin S, Li Y, Llovet JM, Russell DW. Induction of he- to transcriptomic diversity in tumors. Cell Rep 2015;13:267–276. patocellular carcinoma by in vivo gene targeting. Proc Natl Acad Sci U S A [40] Kung C-P, Maggi LB, Weber JD. The role of RNA editing in cancer 2012;109:11264–11269. development and Metabolic Disorders. Front Endocrinol (Lausanne) [51] Toffanin S, Hoshida Y, Lachenmayer A, Villanueva A, Cabellos L, 2018;9:1–21. Minguez B, et al. MicroRNA-based classification of hepatocellular carci- [41] Galeano F, Leroy A, Rossetti C, Gromova I, Gautier P, Keegan LP, et al. noma and oncogenic role of miR-517a. Gastroenterology 2011;140:1618– Human BLCAP transcript: new editing events in normal and cancerous 1628. tissues. Int J Cancer 2010;127:127–137. [52] Pujadas E, Feinberg AP. Regulated noise in the epigenetic landscape of [42] Gromova I, Romov P, Celis JE. bc10 : a novel human bladder cancer- development and disease. Cell 2012;148:1123–1131. associated protein with a conserved genomic structure downregulated [53] Kulis M, Esteller M. DNA methylation and cancer. Adv Genet 2010;70:27–56. in invasive cancer. Int J Cancer 2002;546:539–546. [54] Ramírez de Molina A, Gutiérrez R, Ramos MA, Silva JM, Silva J, Bonilla F, [43] Han L, Diao L, Yu S, Xu X, Li JBJJ, Zhang R, et al. The genomic landscape et al. Increased choline kinase activity in human breast carcinomas: and clinical relevance of A-to-I RNA editing in human cancers. Cancer Cell clinical evidence for a potential novel antitumor strategy. Oncogene 2015;28:515–528. 2002;21:4317–4322. [44] Hu X, Wan S, Ou Y, Zhou B, Zhu J, Yi X, et al. RNA over-editing of BLCAP [55] Lin X-M, Hu L, Gu J, Wang R-Y, Li L, Tang J, et al. Choline kinase a mediates contributes to hepatocarcinogenesis identified by whole-genome and interactions between the epidermal growth factor receptor and mecha- transcriptome sequencing. Cancer Lett 2015;357:510–519. nistic target of rapamycin complex 2 in hepatocellular carcinoma cells to [45] Paz N, Levanon EY, Amariglio N, Heimberger AB, Ram Z, Constantini S, promote drug resistance and xenograft tumor progression. Gastroenter- et al. Altered adenosine-to-inosine RNA editing in human cancer. Genome ology 2017;152:1187–1202. Res 2007;17:1586–1595. [56] de Molina AR, Sarmentero-Estrada J, Belda-Iniesta C, Tarón M, de [46] Chen W, He W, Cai H, Hu B, Zheng C, Ke X, et al. A-to-I RNA editing of Molina VR, Cejas P, et al. Expression of choline kinase alpha to predict BLCAP lost the inhibition to STAT3 activation in cervical cancer. Onco- outcome in patients with early-stage non-small-cell lung cancer: a target 2017;8:39417–39429. retrospective study. Lancet Oncol 2007;8:889–897. [47] Zhao M, Zhang L, Qiu X, Zeng F, Chen W, An Y, et al. BLCAP arrests G1/S [57] Asim M, Massie CE, Orafidiya F, Pértega-Gomes N, Warren AY, Esmaeili M, checkpoint and induces apoptosis through downregulation of pRb1 in et al. Choline kinase alpha as an androgen receptor chaperone and HeLa cells. Oncol Rep 2016;35:3050–3058. prostate cancer therapeutic target. J Natl Cancer Inst 2015;108:djv371. [48] Luk JM, Burchard J, Zhang C, Liu AM, Wong KF, Shek FH, et al. DLK1-DIO3 [58] Hernández-Alcoceba R, Fernández F, Lacal JC. In vivo antitumor activity of genomic imprinted microRNA cluster at 14q32.2 defines a stemlike sub- choline kinase inhibitors: a novel target for anticancer drug discovery. type of hepatocellular carcinoma associated with poor survival. J Biol Cancer Res 1999;59:3112–3118. Chem 2011;286:30706–30713. Journal of Hepatology 2020 vol. 73 j 328–341 341