A MicroRNA Cluster in the DLK1-DIO3 Imprinted Region on Chromosome 14q32.2 Is Dysregulated in Metastatic Hepatoblastomas

Hepatoblastoma (HB) is the most common malignant liver neoplasm in children. Despite progress in HB therapy, outcomes for patients with metastatic disease remain poor. Dysregulation of miRNA expression is one of the potential epigenetic mechanisms associated with pathogenesis of HB. However, miRNA profiles related to the different stages of HB tissues and cells, in particular of lung metastatic tumor cells, are unknown. In the present study, using array-based miRNA expression and DNA methylation analysis on formalin-fixed paraffin-embedded tissues, we aimed to identify miRNA changes that can discriminate between lung metastatic tumors, primary tumors (fetal and embryonal subtypes), and nontumorous surrounding livers. Our analysis demonstrated that a large cluster of microRNAs and snoRNAs located within the 14q32.2 DLK1-DIO3 region showed a strikingly upregulated expression pattern in HB tumors, especially metastatic tumors, compared to normal liver tissues. This revealed dysregulation of miRNAs similar to that seen in a malignant stem-like subtype of hepatocellular carcinoma associated with poor prognosis. These findings in HB mirror similar findings made in multiple other cancer types. With further analysis this may in future allow stratification of different stages and types of HB tumors based on their miRNA profiles, which could lead to new approaches to diagnosis and treatment in progressive HB patients.


INTRODUCTION
Hepatoblastoma (HB) is the most common malignant neoplasm of the liver in children. Despite progress in the therapy of HB, the outcome of patients with metastatic disease remains poor (1,2). Approximately one fifth of HB patients have pulmonary metastasis at diagnosis, and recurrence of HB most frequently occurs in the lung. Among those patients who had residual lung disease after induction chemotherapy, but were able to undergo complete resection of the liver tumor, the overall survival rates in those whose lung tumors were completely resected, and of those whose lung tumors were incompletely resected, were 63.6% and 41.8%, respectively (3). To improve the prognosis, identification of novel molecular-genetic markers predictive of treatment outcome, and innovations that would allow personalized treatment are needed.
It is well known that mutations in beta-catenin are a hallmark of HB (4,5). In contrast, results from recent whole-exome sequencing analyses showed other mutations are rarely seen in HB (3.9 -4.6 mutations per tumor) (6,7). The paucity of mutations has been reported in many pediatric tumors and may be correlated with their early age of onset (8). This paucity of mutations might also suggest that epigenetic aberrations are an important mechanism involved in the pathogenesis of HB.
Small non-coding RNAs, including microRNA (miRNA) and small nucleolar RNA (snoRNA), have been the focus of many studies in the last few decades and their fundamental role in cancer is currently well established. MiRNAs are a class of 19-25nt single-stranded RNA molecules that serve as major regulators of gene expression through their ability to bind to and post-transcriptionally inhibit the expression of specific target messenger RNAs. SnoRNAs are a class of 60-300 nt non-coding RNAs implicated in the chemical modification of ribosomal RNA, and their role in cellular regulation as well as a role in cancer development and progression has been highlighted (9).
We previously reported that aberrant DNA methylation of some tumor suppressor genes was related to poor outcome in HB patients and that disruption of imprinting status was implicated in its pathogenesis (10)(11)(12). Previous work on miRNA alterations in HB has shown that clustering of some miRNA expression profiles were related to several clinical aspects, including prognosis (7,(13)(14)(15). Therefore, dysregulation of epigenetic mechanisms, particularly miRNA expression and DNA methylation, could be a useful parameter for diagnosis as well as classification of HB. Primary liver cancers such as HB or hepatocellular carcinoma (HCC) often comprise heterogeneous types of cells with different histological features. HB tumors are classified as wholly epithelial, or of mixed epithelial and mesenchymal types. In the wholly epithelial type, there are two major subtypes, called the fetal subtype and the mixed fetal and embryonal subtype. Cairo et al. identified a 16-gene signature discriminating tumors with a fairly well-differentiated histology and a favorable prognosis against advanced and poorly differentiated tumors with a dismal outcome (16). Their signature defines HB prognostic subtypes that reflect liver developmental stage, which is in line with the ability of the miR signature to stratify patients with HB (13). A better understanding of the difference in molecular profiles related to the heterogeneity can aid in exploring molecular mechanisms of carcinogenesis and therapeutic options. Moreover, comparative molecular analyses between primary and metastatic tumors are useful to identify important molecules and pathways related to tumor progression. Therefore in the present study, we evaluated miRNA expression profiles related to different tissues and cell types in HB by using formalin-fixed, paraffin-embedded (FFPE) specimens. We previously confirmed that high quality data can be generated in HB FFPE sample analysis through miRNA array assays (17). The present study has identified remarkable dysregulation of microRNAs and snoRNAs located within the 14q32.2 DLK1-DIO3 imprinted region, especially in metastatic HB tissues. The pattern of the miRNA dysregulation observed was typical of that seen in a malignant stem-like subtype of HCC associated with poor prognosis.

Patients
FFPE specimens were obtained from 16 patients referred to Hokkaido University Hospital and Kanagawa Children Medical Center for surgical treatment in between 2000 and 2016. Ten patients were male and six were female. Fourteen patients with a median age of 3.9 (0-10.5) years underwent primary tumor resection. Nine patients had lung metastasis (7, synchronous; 2, metachronous), and the nine metastatic tumor samples were resected at a median age of 7.6 (2.8-10.9) years ( Table 1). The ethics committee at our institution approved the study protocol. For 6 of the patients, informed signed consent was provided by the parents. For the remaining 10 patients, the institutional ethics committee approved a web-based opt-out consent process, whereby each patient's parents had the option of withdrawing the patient's tissues and clinical information from the study.

RNA Extraction and miRNA Array
The FFPE specimens included tumor cells from 9 metastatic tumors, 21 primary tumors (11 fetal subtypes and 10 embryonal subtypes), and normal cells from 14 nontumorous surrounding liver samples. After dissecting the specimens macroscopically under a light microscope with a sterile scalpel, ensuring tissues were collected well away from tissue boundaries and stored systematically, which enabled us to avoid contamination by normal tissues or mesenchymal components, we extracted 44 RNA samples using a miRNeasy FFPE Kit (Qiagen, Hilden, Germany) according to the manufacturer's instructions. We conducted a quality check (QC) of the RNA samples measuring the concentration of total RNA and OD 260/280 ratio using DS-11 NanoPad (DeNovix, Wilmington, USA) or NanoDrop ND-2000 (Thermo Fisher Scientific, Waltham, USA), and RNA integrity number (RIN) using Agilent 2100 Bioanalyzer using Agilent RNA 6000 Nano Kit (Agilent Technologies, Santa Clara, USA). All the data for QC is shown in Supplementary Table 1. The OD 260/280 ratios for all the

DNA Extraction and Methylation Array
We extracted DNA samples from the tissue FFPE samples mentioned above. To extract the DNA, we used a QIAamp DNA FFPE Tissue Kit (Qiagen) according to the manufacturer's instructions. We conducted a quality check of the DNA samples using RT-PCR, following the Infinium HD FFPE QC Assay Protocol, and we confirmed that 12 samples from six nontumorous surrounding liver (8N, 9N, 10N, 11N, 12N, and 14N in Table 1) and six metastatic tumor (4M, 8M, 10M, 12M, 14M, and 15M in Table 1) were appropriate for the methylation assay. We next conducted genome-wide methylation analyses using an Infinium HumanMethylation450 BeadChip (Illumina) and the 12 DNA samples, following the Illumina Infinium HD Methylation protocol. This array includes 485,577 cytosine positions in the human genome (482,421 CpG sites (99.4%), 3,091 non-CpG sites and 65 random SNPs). We linked the UCSC Genome Browser annotation (version hg19 of the human reference genome) to each of the CpG sites on the array. Based on the UCSC chromosome annotation, we then screened for probes that were located in the DLK1-DIO3 imprinted region.

Data Analysis and Statistical Analysis of miRNA Data
A workflow describing the study design, which involved the miRNA and snoRNA microarray data analysis is shown in Supplementary Figure 1. All raw data were normalized by the use of a tool that uses Robust Multichip Analysis (RAM) plus detection above background (DABG) algorithms and analyzed using Transcriptome Analysis Console ver. 4.0 (Affymetrix, Santa Clara, CA). Statistical analysis and data visualization of heatmap with hierarchical clustering were performed by Qlucore Omics Explorer (QOE) (Qlucore AB, Lund, Sweden). The differential expression of miRNAs in the different tissue cell types was analyzed using multi-group comparisons with ANOVA (false discovery rate (FDR) < 0.05). Variance filtering was used to reduce the noise and the filtering threshold was set of 0.52 applied by QOE. To compare differentially expressed miRNAs between normal surrounding liver tissues and primary or metastatic tumors, a two-group comparison (with > 2 fold mean expression change and one-way ANOVA, FDR < 0.05) was used. Data visualization of box plots was performed using GraphPad Prism software (version 6, La Jolla, CA) and P-values were calculated by Tukey's test. The average Euclidean distance was adopted for hierarchical clustering. To perform pathway enrichment analysis and targets prediction, we used Diana mirPath v.3 web-based computational tool (http://microrna.gr/ mirpath). We used the "pathways union" option of the miRPath software. P-values were obtained by the Fisher's exact test as enrichment analysis method and the FDR was estimated using the Benjamini and Hochberg method. All pathways showing Pvalues < 0.05 were considered significantly enriched between the groups under comparison (22). To represent the HB risk stratification as described by Cairo et al. (2010), we classified the HB tumors into C1/C2 subclasses that are related to Cm1/Cm2 subclasses determined by miR expression profiles (13). In Table 1, C2 tumors were defined as tumors which showed mir-371 and/or 373 upregulation, compared to the matched nontumorous surrounding livers.

Identification of Differentially Expressed MiRNAs in Four Different Tissue Types in HB
Following the strategy outlined in the workflow diagram (Supplementary Figure 1), we compared miRNA expression levels in four different tissue and cell types in HB: nontumorous liver tissue surrounding HB (n=14), primary HB tumor tissues of two subtypes [fetal (n=11) and embryonal (n=10)], and metastatic tumor tissue (n=9), and detected 73 differentially expressed miRNAs in multi-group comparisons using ANOVA (FDR < 0.05). The genomic locations, fold changes, and statistical values of the miRNAs are shown in Table 2. The results of hierarchical clustering of the expression of these 73 miRNAs in the 44 tissue samples are presented in the heatmap diagram ( Figure 1), showing that, with a few exceptions, almost all the nontumorous surrounding liver samples expressed most of the identified miRNAs at relatively low levels, and that these 73 miRNAs became upregulated in primary and metastatic tumors. The exceptions included one nontumorous surrounding liver sample that expressed most of the miRNAs at relatively high levels, and one embryonal tumor plus two metastatic tumors that expressed a subgroup of the miRNAs at relatively low levels, while five nontumorous surrounding liver samples expressed a smaller but different subgroup of the miRNAs at relatively high levels. Hierarchical clustering showed that the metastatic tumors, and the fetal and embryonal subtypes of primary tumors did not segregate from each other. Overall, these results suggest that by hierarchical clustering the identified differentially expressed miRNAs were able to discriminate between normal surrounding liver tissue and HB tumor tissue (whether it be primary or metastatic), but that clustering of the miRNA expression was unable to discriminate between metastatic tumors and the primary tumors of the two subtypes ( Figure 1). In addition to hierarchical clustering, we compared the miRNA expression profiles of nontumorous surrounding liver vs. primary and metastatic tumors in a two-group comparison (with > 2-fold mean expression change and one-way ANOVA, FDR < 0.05), and identified 60, 51, and 81 differentially expressed miRNAs in the fetal subtype, embryonal subtype and metastatic tumor tissue compared to normal tissue, respectively (Supplementary Figure 2). However, except for several miRNAs that were expressed at relatively low levels, we ultimately were unable to identify single miRNAs using this approach that were significantly differentially expressed in metastatic tumors versus either fetal or embryonal subtype tumor.

Pathway Enrichment Analysis and Target Prediction
Next, we aimed to investigate whether specific pathways involving up-and down-regulated miRNAs were altered in primary and metastatic HB tumors. For this analysis, we calculated the enrichment of multiple miRNA target genes and compared each set of miRNA targets to all known Kyoto Encyclopedia of Genes and Genomes pathways to identify altered pathways. We first investigated whether co-expression of the top 21 upregulated miRNAs (highlighted as bold in Table 2, and selected on the basis of fold-difference in expression) in metastatic tumors compared to those in the normal surrounding liver tissue could affect some signaling pathways. The miRNA versus pathway heatmap showed that signaling pathways regulating Hippo (Pvalue, 0.024), also known as the Salvador/Warts/Hippo pathway, which controls organ size in animals (5), were significantly enriched in the case of the upregulated miRNAs ( Figure 2). Moreover, when we investigated all the upregulated genes in Table 2 in the pathway enrichment analysis, 18 KEGG pathways were significantly enriched (Supplementary Table 2), which included Hippo (P-value, 8.28e-10), TGF-beta (P-value, 0.0073), and p53 (P-value, 6.02e-05) signaling pathways. On the other hand, when we examined the targets of the four downregulated miRNAs, miR-378d, miR-378a-5p, miR-146b-3p and miR-203a, in metastatic tumors versus normal surrounding tissue for enriched functional pathways, this yielded biosynthesis of unsaturated fatty acids (hsa01040) at the top of the list (P-value, 0.033).
DLK1-DIO3 Genomic Imprinted miRNA Cluster at 14q32.2 In a previous study (17), we had noted that chromosome 14 contained the highest percentage (~15%) of the total complement of microRNAs expressed in HBs upon mapping of miRNAs to each chromosome, and that these miRNAs were concentrated in a relatively small region of chromosome 14, the DLK1-DIO3 imprinted region, which contains a miRNA cluster at 14q32.2. We therefore examined the chromosomal locations of the differentially expressed miRNAs in our tumor samples to identify whether this might reveal miRNA clusters enriched in a specific genomic location. Surprisingly, we found 55% (40 out of the 73) of all the identified miRNAs were located within the DLK1-DIO3 region on chromosome 14q32.2 ( Table 2). Moreover, 19 of the top 21 miRNAs investigated in Figure 2 are located within the DLK1-DIO3 region on chromosome 14q32.2. The miRNAs that cluster in the 14q32.2 region are separated by a cluster of putative C/D box small nucleolar RNAs, therefore, we evaluated the expression of 93 miRNAs and 54 small nucleolar RNAs (snoRNAs) located in this region in the tumor samples and compared that with the normal surrounding liver samples to identify patterns of differential expression. It showed that 54 of these miRNAs (58.0%) and 31 small nucleolar RNAs (57.4%) were upregulated with > 2-fold mean expression change in at least one type of tumor tissue. As shown in Figure 3, the increased miRNA and snoRNA expression levels, which were calculated by weighted average using Transcriptome Analysis Console ver. 4.0, were overall the highest in the metastatic     tumors. This was most apparent in the snoRNAs, which were highly upregulated with greater than 50-fold mean expression change in at least one type of tumor tissue. Comparing mean fold differences in expression, differential expression of miRNAs and snoRNAs occurred at the DLK1-DIO3 locus in HB metastatic and primary tissues versus normal liver tissues (Supplementary Figure 3).
To examine whether metastasis-associated changes in miRNA expression occurred, we used matched trios of metastatic, primary, and normal liver tissues derived from the same patients. Relative ratios of miRNA expression in the matched tissues were calculated in primary tumors (n=21) and metastatic tumors (n=6) divided by expression in the matched normal liver tissue. In this analysis, we excluded three metastatic  tumor samples (21M, 23M-1, and 23M-2) because of lack of the matched normal samples. Comparison of the ratios showed that 23 miRNAs and 8 snoRNAs were significantly upregulated in the metastatic compared to primary tumor tissues (Supplementary Table 3, and Supplementary Figure 4).
To determine whether the miRNAs that were upregulated in this region had any specific function, we investigated all of the 40 identified upregulated miRNAs using the pathway enrichment analysis tool. This showed that only two KEGG pathways, Prion diseases (hsa05020) (P-value, <1e-325) and ECM-receptor interaction (hsa04512) (P-value, <1e-325), were significantly enriched, both of which were similarly yielded at the top of the list in the previous analyses ( Figure 2 and Supplementary Table 1).
We then went on to further analyse whether differential methylation between metastatic tumor tissues and normal liver tissues occurred at two regulatory DMRs (Differentially Methylated Regions), which are responsible for controlling the imprinted expression of genes, including MEG3, in the DLK1-DIO3 imprinted region. We observed that the methylated CpG dinucleotides in these two DMRs all displayed an average loss of DNA methylation when analysed in a series of six metastatic tumour tissues, as compared to DNA methylation levels in six normal liver tissues (Figure 4). Three CpG probes located in the MEG3-DMR showed significantly lower methylation levels in metastatic tumor tissues compared to normal liver tissues using Mann-Whitney U analysis. As shown in Figure 5 and Supplementary Table 4, low methylation levels at the three probes were significantly correlated with the upregulation of the miRNAs and snoRNAs shown in Figure 3. While these results are not sufficient to demonstrate a causal association between hypomethylation and increased miRNA expression from the DLK1-DIO3 imprinted region, they are nevertheless consistent with elevation of miRNA and snoRNA expression levels in the DLK1-DIO3 imprinted region potentially resulting from an overall loss of DNA methylation in the locus. Also consistent with our results, Luk et al. (2011) had shown that miRNAs in the DLK1-DIO3  imprinted region were overexpressed in a stem-like subtype of hepatocellular carcinoma (23), and very recently Carrillo-Reixach et al. (2020) have reported that genes in the DLK1-DIO3 imprinted region are hypomethylated and overexpressed in HB, enabling stratification of patient outcome (24).

Expression of miRNAs Participating in Key Signaling Pathways Involved in the Pathogenesis of HB Tumors
Several studies have explored the role of miRNAs in some key signaling pathways altered in the pathogenesis of HB tumors (14) and references therein. The known miRNAs participating in a key signaling pathway, the Hippo pathway, include miR-186, miR-125a, and miR-206 (14). We compared the expression levels of two oncogenic miRNAs (miR-125a and miR-206) and one tumor-suppressive miRNA (mir-186), by comparing the probe intensities in each tumor tissue type to that of normal surrounding liver tissue in each of the 14 cases where this comparison could be done ( Table 1). As shown in the coloredcell patterns in Table 1, we found that miR-125a-5p was upregulated in 10 out of the 14 cases (71%) in at least one tumor tissue type compared to normal surrounding liver. In contrast, changes in the expression of miRNAs participating in the Wnt (miR-4510, let-7i-3p, miR-624-5p, and miR-885-5p), and Myc (miR-371 cluster and miR-100/let-7a-2/miR-125b-1 cluster) signaling pathways, which have also been highlighted in HB (5), were unremarkable, except miR-885-5p, which was downregulated in 10 and upregulated in 4 out of the 14 cases in at least one tumor tissue type compared to normal surrounding liver, respectively. In the present analysis, we found no apparent associations between any alterations in those miRNA expressions and risk stratification factors, such as C1/C2 classification, or age at diagnosis and the presence of metastasis ( Table 1). On the other hand, when we compared the expression levels of 14q32 transcripts in the DLK1-DIO3 imprinted region with these same prognostic factors in the 21 primary tumor samples, we identified seven miRNAs and two snoRNAs, which were significantly upregulated in the six primary tumors that were classified into the C2 classification (Supplementary Figure 5). Furthermore, in the seven primary tumors which had metachronous or synchronous metastasis, three miRNAs and one snoRNA were significantly upregulated (Supplementary Figure 6). There were no differentially upregulated 14q32 small RNAs when comparing age at diagnosis of more than versus less than 8 years of age.

DISCUSSION
HB tumors demonstrate a diverse range of histological phenotypes and clinical behaviors, which may arise as a result of the proliferation of transformed liver stem cells or of early hepatic precursors with varying degrees of differentiation (7). Furthermore, the presence of metastasis is the most important prognostic factor for HB patients (1)(2)(3). The use of FFPE specimens has enabled us to separately collect different types and stages of HB tumor tissues. Our analysis has demonstrated that HB shows distinct miRNA profiles related to both normal and HB tumor type and stage.
Our data provide the first description, to our knowledge, of miRNA profiles investigated in primary tumour and lung metastatic HB tumors. However, the main limitation of our data was that, due to HB being a very rare tumor, only a relatively small number of the clinical samples were examined. Nonetheless, miRNAs and snoRNAs located in the 14q32.2 cluster showed significantly upregulated expression in metastatic tumors. Interestingly, our results resemble recent findings within a stemlike subtype of HCC, which is known to be associated with overexpression of miRNAs located in this cluster (23), and they are also very similar to a recent report showing that genes in the DLK1-DIO3 imprinted region are hypomethylated and overexpressed in HB, enabling stratification of patient outcomes (24). These studies suggest that certain molecular mechanisms, such as Wnt signaling pathways or TGF-beta, might be related to the tumorigenesis of this unique stem-like subtype of HCC. Indeed, in their very recent paper Carrillo-Reixach et al. (24) investigated a large cohort of HB patients, and demonstrated that overexpression and hypomethylation of genes in the 14q32 DLK1-DIO3 imprinted region, which includes the largest known cluster of miRNAs and snoRNAs in the human genome, is associated with HB patient outcome, and activation of the Wnt/beta-catenin pathway. Activation of the Wnt/b-catenin is a hallmark in HB (4,25), and the Hippo pathway has also been reported in HB (26,27). Hippo signaling transcription factors YAP and TAZ are known to regulate the b-catenin destruction complex and to orchestrate Wnt responses (28,29). These findings suggest that dysregulation of 14q32 clusters might play an important role in HB pathogenesis via interaction with Wnt/b-catenin/Hippo signaling pathways. The strong upregulation of imprinted genes in the 14q32 locus, DLK1, and MEG3, which are abundantly expressed in fetal liver, has also been shown in HB (16). Products of genes targeted by miRNAs from this locus are expected to regulate receptor tyrosine kinase-activated pathways in the pathogenesis of relapsing-remitting multiple sclerosis (30).
Several previous studies investigating copy-number variations in HB showed that copy-number alterations in the chromosomal region are rarely seen (7,31), suggesting that it is unlikely that gene amplification was the cause of upregulation of the genes and miRNAs at this locus. This is also supported by our findings that reduced average methylation was identified at two differentially methylated regions (DMRs), which have previously been shown to control transcription from imprinted genes in the locus, including MEG3. In contrast to our findings, Oshima G, et al., showed that expression of 14q32-encoded miRNAs was a favorable prognostic factor in patients with metastatic colorectal cancer (32). Further investigations of aberrant methylation at imprinting control regions are warranted to determine whether disruption of the imprinting status of this locus is involved in the pathogenesis of HB. In addition, HCC patients with overexpressed 14q32.2 miRNAs exhibited shorter survival times (23), therefore, miRNAs from this locus might also become prognostic biomarkers and candidates for targeted therapies designed to treat HB patients with progressive disease. These miRNAs should therefore be further explored using more clinical specimens for their diagnostic and prognostic use in HB.
Among miRNAs that participate in key signaling pathways previously identified as being involved in the pathogenesis of HB tumors, miR-125a-5p might frequently function as an oncogene in the development of both primary and metastatic HB tumors. Suppression of miR-125a-5p upregulates the Tafazzin (TAZ)-EGFR signaling pathway, and has been shown to promote retinoblastoma proliferation. It has also been suggested to have a potential role in regulating the Hippo pathway (19). In summary, new findings involving miRNAs could serve to provide alternative therapeutic strategies for the treatment of HB, especially in patients with metastatic disease.

CONCLUSION
Stratifying different types and stages of primary and metastatic HB tumors based on their miRNA profiles could lead to new approaches to diagnosis and treatment in HB patients with progressive disease. The present work has identified a surprisingly high proportion of upregulated miRNAs located within the 14q32.2 locus in HB tumors, and particularly in metastatic tumors. Discovering more about the upregulation of both miRNAs and snoRNAs on the 14q32.2 locus could improve our understanding towards innovating better therapeutic strategies to improve clinical outcomes of progressive HB, as well as potentially elucidating their role in stem-like subtypes of HCC.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: the NCBI Gene Expression Omnibus (GSE153089).

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by Hokkaido University Hopital Ethics Committee. Approval code is 010-0202. Written informed consent to participate in this study was provided by the participants' legal guardian/next of kin.