A transgenic zebrafish liver tumor model with inducible Myc expression reveals conserved Myc signatures with mammalian liver tumors

Liver cancer, mainly hepatocellular carcinoma (HCC), is one of the leading causes of cancer-related death worldwide (Nordenstedt et al., 2010). It generally has a poor prognosis as it is often diagnosed at an advanced stage when treatment is not effective. Tremendous efforts have been made to decipher the molecular mechanisms of HCC, and gene expression profiling of human HCCs has been used to identify subgroups of patients according to etiological factors, early pre-neoplastic lesions, stages of the disease, rate of recurrence and survival (Hoshida et al., 2009; Lee and Thorgeirsson, 2004; Roessler et al., 2010; Wurmbach et al., 2007). In human HCC, MYC is commonly amplified and associated with unfavorable prognosis (Abou-Elella et al., 1996). A higher expression level of MYC is associated with more advanced status of HCC (Gan et al., 1993). However, at which stage and to what extent Myc contributes to human HCC is still unclear. Recently, Myc has been suggested to be at the center of human liver tumor malignant conversion, on the basis of genome-wide gene expression profiling (Kaposi-Novak et al., 2009), but this conclusion still lacks experimental confirmation and thus an in vivo animal model is highly desired.

In the past decade, the zebrafish has become an increasingly popular experimental model for human diseases (Lieschke and Currie, 2007; Liu and Leach, 2011; Payne and Look, 2009). We have demonstrated that the zebrafish and human liver tumors have molecular conservation at various levels, indicating the potential of zebrafish for modeling human liver cancer (Lam et al., 2006). Although several Myc transgenic mouse models for liver cancers have been reported (Lee et al., 2004), we envision that comparative studies of Myc tumors in an evolutionary distant model such as zebrafish should disclose highly conserved features and biomarkers of Myc tumors. Thus, in the present study, we have generated an inducible Myc transgenic zebrafish line and demonstrated that overexpression of Myc in the liver results in obvious liver tumors. Furthermore, by transcriptomic analyses, we found that the zebrafish Myc tumor model not only captured the human liver cancer signature but also showed high similarity with mouse Myc-induced liver tumor models. We also identified a short list of Myc targets that could be used to predict Myc-causing malignancies. Our analyses suggested that this transgenic zebrafish line could serve as a model for human liver tumors induced by over-activation of MYC.

In order to investigate the role of Myc in liver tumorigenesis, a transgenic line, named TO(Myc), was generated by introducing mouse Myc cDNA in a Tet-On vector under a liver-specific zebrafish fabp10 promoter (Her et al., 2003). By whole mount in situ hybridization, we confirmed that transgenic Myc expression was specifically induced in the liver after doxycycline (Dox) treatment (Fig. 1A,B). To investigate the effect of Myc overexpression on liver growth, TO(Myc) fish were crossed with the LiPan line, which has red fluorescent protein (RFP) expression in the liver to facilitate observation of liver development (Korzh et al., 2008). As shown in Fig. 1D, overgrowth of the liver was apparent after 4 days of Dox (40 μg/ml) treatment starting from 3 dpf (days postfertilization), in comparison with that of untreated controls (Fig. 1C). By using different concentrations of Dox, we observed a clear dose-dependent increase of the liver size in TO(Myc) embryos, as indicated by two-dimensional (2D) measurement (Fig. 1E). By contrast, no obvious change was noted in the non-TO(Myc) siblings (M-D+) even under high dosage (40 μg/ml) of Dox treatment. Consistent with the observations, the dose-dependent induction of Myc transgene expression was also confirmed by reverse transcription polymerase chain reaction (RT-PCR), in which the transgenic embryos were treated with increasing concentrations of Dox (10, 30 and 60 μg/ml) for 2 days from 2 dpf (Fig. 1F). Although there was a dosage-dependent increase of proliferation in the TO(Myc) livers after Dox induction, there was no detectable apoptosis in the TO(Myc) livers at all concentrations of Dox from 10-60 μg/ml (data not shown).

To examine the tumorigenic potency of the TO(Myc) line, TO(Myc) transgenic zebrafish and their non-transgenic siblings were treated with Dox (60 μg/ml) from 21 dpf. All Dox-treated transgenic fish (M+D+) showed enlarged abdomens (Fig. 2D) and obvious liver overgrowth (Fig. 2E) after only 2-3 weeks of Dox induction (100%, n=50) compared to Dox-treated non-transgenic siblings (M-D+) (Fig. 2A,B). Histological analysis of the M+D+ fish showed hyperplasia (HP; 6/10), mixed hyperplasia/ hepatocellular adenoma (HP/HCA; 1/10) and adenoma (3/10) (Fig. 2F; Table 1). In general, the entire liver was uniformly transformed into tumor phenotype without the existence of normal liver tissue. For observation of the effect of prolonged activation of Myc, we treated the fish continuously for up to 16 weeks. The treated fish (M+D+) did not show significant increase in mortality (Fig. 2M), but they did show more prominent liver enlargement (Fig. 2J,K) and significantly smaller body size (Fig. 2J,N) than control fish. Hematoxylin and eosin (H&E) examination of samples from 8 and 16 weeks post-treatment (wpt) revealed a progression of liver tumors from predominantly hyperplasia to adenoma (Table 1). We also noticed one case of early HCC, as shown in Fig. 2L, confirming the liver tumor progression with prolonged Myc activation.

In order to investigate the molecular mechanism of Myc-induced liver tumors, livers from all four groups (transgenic and non-transgenic fish, treated or untreated with Dox) were collected after 4 months of Dox-treatment at 60 μg/ml. Two biological replicates (each with four or five pooled livers) were included for each group. About 40 million tags were sequenced for each sample mapped to the zebrafish RefSeq database for quantitative analysis of transcript expression (supplementary material Table S1). To verify the dynamic range of gene expression, reverse-transcription quantitative PCR (RT-qPCR) was performed for several selected genes with a range of 249,696 to 0.12 TPM (transcripts per million) in the control sample 1 (M-D-_1). As shown in Fig. 3A, there was a good correlation of the RNA sequencing data and RT-qPCR results.

To compare the molecular similarity of the tumor and control samples, whole transcriptomic data from all eight samples were subjected to hierarchical clustering using MeV software (Saeed et al., 2006; Saeed et al., 2003). As shown in Fig. 3B, the two tumor samples (M+D+) were clustered as an independent branch and were well separated from all control samples (M-D-, M-D+ and M+D-), suggesting that there were substantial changes in the tumor transcriptome.

In order to identify the genes important for liver cell transformation, a t-test was performed between the six control and two tumor samples. Differentially expressed genes were selected by accepting a P-value of <0.05 and fold change >1.5 for further analyses. In order to focus on physiologically relevant genes, we also set an expression level cutoff at an average of 10 TPM in either control or tumor samples, which is equivalent to about three transcripts per cell, based on an estimation of 300,000 transcripts per cell (Mortazavi et al., 2008). Finally, 148 up- and 353 downregulated transcript entries were selected using these criteria (supplementary material Tables S2, S3). These up- and downregulated zebrafish genes were subsequently mapped to 125 and 207 human gene symbols. As illustrated in the volcano plot (Fig. 3C), there were fewer upregulated genes than downregulated genes, but the upregulated genes were more deregulated as judged by their distances from the base point.

To understand the composition of these differentially expressed genes, Gene Ontology Slim classification analysis was performed (Fig. 4A-C). Interestingly, all genes in the macromolecular complex and ribosome groups were upregulated (Fig. 4B). By contrast, more genes involved in nucleotide binding, ion binding, protein binding, transferase activity, hydrolase activity, transporter activity, membrane and nucleus were downregulated. It is also interesting to note that there were more uncharacterized genes in the downregulated group than the upregulated group, probably implying that less attention was devoted to the downregulated genes in the Myc pathway and cancers in previous studies.

To gain more detailed functional insights into the differentially expressed gene list, gene ontology enrichment analysis was conducted. As shown in Table 2, ribosome and translation are the most prominent groups in the upregulation division. Consistent with this, all of the 59 ribosomal protein genes identified showed upregulation in the tumors (Fig. 4D). Of these, 34 genes (rpl4/5a/7a/8/9/10a/11/12/13/13a/18a/19/21/22/26/27a/34/36/36a and rps3/5/7/8/11/12/13/14/15/16/20/21/23/25/29) and several major translation factors (eif2s1/3, eif3c/d/ea/ha, eef1g and ef1a) were present in the significantly upregulated group in the tumors (supplementary material Table S2).

We also compared the expression of endogenous myc genes. Out of the six members of the myc gene family in zebrafish (myca, mycb, mych, mycl1a, mycl1b and mycn), only mycb and mych were relatively highly expressed in the liver, with averages of 54.9 and 13.3 TPM, respectively, in the control samples (M-D-, M-D+ and M+D-) although these decreased to 13.3 and 3.5 TPM, respectively, in the liver tumor samples (M+D+). Thus, it seems that high level of transgenic expression of the mouse Myc negatively affects the expression of endogenous myc genes. By contrast, all the rest of the myc genes had extremely low expression level (0.0 to 2.3 TPM) in all the liver samples (control and tumor). Thus, it is impossible to evaluate their change during tumorigenesis.

We next examined the similarity of Myc zebrafish model to various mouse HCC models, including two chemically induced [ciprofibrate (CIP) and diethylnitrosamine (DENA)], four transgenic (overexpression of Myc, E2f1, Myc/E2f or Myc/Tgfa in the liver) and one knockout (Acox1^−/−) models (Lee et al., 2004). Hierarchical clustering analysis of gene expression patterns was used to assess the relative similarities between different mouse HCC models. Three distinctive HCC clusters were identified. Our Myc overexpression-induced liver cancer clearly distinguished the Myc-E2f1-Myc/E2f1 cluster from the other models, and it showed the highest relatedness to Myc/E2f1 transgenic mice, suggesting conservation of Myc signatures between zebrafish and mice liver tumors induced by overexpression of Myc (Fig. 5A; supplementary material Table S4). Thus, our zebrafish model successfully resembles mouse liver tumor under similar a inducer (here Myc) and confirmed conservation between zebrafish and mouse liver cancers.

Hepatocarcinogenesis is a multistage process following a dysplasia-adenoma-carcinoma sequence with distinct morphological and molecular signatures. Therefore, we compared the transcriptome of Myc-induced zebrafish liver tumors with those from different stages of human liver conditions: low grade dysplastic nodules (LGDN), high grade dysplastic nodules (HGDN), and very early, early, advanced and very advanced hepatocellular carcinoma (veHCC, eHCC, aHCC and vaHCC) (Table 3). In the first human HCC dataset (GSE12443) (Kaposi-Novak et al., 2009), the upregulated genes of Myc-induced zebrafish tumors clearly distinguished HGDN and eHCCs from others, and they showed significant similarity to eHCC. Another human HCC dataset (GSE6764) included more advanced HCC stages, in which the pathological HCC stages were defined by tumor size, differentiation status and metastasis level (Wurmbach et al., 2007). Similarly, the upregulated genes of the zebrafish tumors successfully captured human HCC but not the less severe groups including HCA (Fig. 5B; Table 3). Within the HCC group, it showed the highest similarity with very advanced human HCC. Analyses of the leading edge genes revealed that most of them are involved in mRNA transcription (SND1, BTF3, MED30, EXOSC3), rRNA processing (EXOSC3, NAS2), protein translation and folding (RPSs, RPLs, EIF3, EEF1, CCT2, CCT6A, CCT8), cell cycle regulation (CAPRIN1, GNB2L1, CGREF1, DTYMK, PCNP) and mitochondrial function (SND1, CHCHD3, PDSS, SLC25A39) (supplementary material Table S5).

Myc is a global transcription factor, therefore it is interesting to examine the Myc target genes involved in hepatocarcinogenesis. Thus, the upregulated genes from the Myc-induced zebrafish liver tumors were compared with Myc target genes retrieved from the Molecular Signatures Database (MSigDB) (Subramanian et al., 2005), where 382 Myc targeted genes could be mapped by the zebrafish transcripts in the RefSeq database. We found that 16 of them presented in the 125 upregulated genes from the Myc-induced zebrafish tumors, including genes encoding ribosomal proteins (RPL13A, RPL19, RPL22, RPS2, RPS11, RPS28), mitochondria (FPGS, SLC7A3, ALDOA), translation (PABPC4, ANKHD1), cell cycle inhibition (DUSP1, CGREF1), glycogen biosynthesis (PPP1R3C), transcription factor (TEF) and a member of the arrestin family (ARRDC3). Notably, DUSP1 and ARRDC3 have been suggested as tumor suppressors (Calvisi et al., 2008; Draheim et al., 2010). We also performed gene set enrichment analysis (GSEA) using these 16 genes as Myc signature to classify different mouse HCC models and various stages of human HCCs (Fig. 5B; Table 3). The 16-gene signature showed a better prediction of Myc transgenic mice (Fig. 5A). Interestingly, these 16 genes also clearly distinguish HCC with less severe human liver tumors, with better correlation to more advanced HCC (Fig. 5B; Table 2). Therefore, the 16-gene signature could potentially be used as a molecular prognosis signature for human HCC patients.

Similarly, we also compared the 125 upregulated genes from our data with another set of Myc target genes (Zeller et al., 2003) (http://www.myc-cancer-gene.org/index.asp), whereby 681 genes were identified with zebrafish homologs and 24 of them were overlapped (indicated by asterisks in supplementary material Table S2), including five genes (RPL13A, RPL19, RPL22, RPS11 and FPGS) in the earlier-identified 16-gene list. We found that the 24 genes also showed a better prediction of Myc transgenic mice (supplementary material Table S4) and high correlation with advanced and very advance HCC (supplementary material Table S5) as compared with the entire upregulated gene list.

By generating a transgenic zebrafish model and utilizing RNA-sequencing technology, we have examined the transcriptome of zebrafish liver tumors driven by the oncogene Myc. Our analyses have shown that ribosome biogenesis constitutes the key features of Myc-mediated liver tumors. Cross-species analyses showed that the zebrafish Myc model correlated well with Myc transgenic mouse models for liver cancer. We also illustrated that the Myc-induced zebrafish liver tumors possessed molecular signatures common to human HCC. Interestingly, a small set of 16 Myc target genes was identified that could be further tested as HCC prognosis markers in clinical samples.

Consistent with previous reports that MYC is the key regulator of ribosome gene expression and protein translation (van Riggelen et al., 2010), we observed that essentially all ribosomal protein genes as well as several translation factor genes were upregulated. The ribosome is the machinery for protein synthesis; therefore, upregulation of ribosome gene expression could reflect the requirement for increased cell proliferation. This is in contrast to our two other liver tumor models induced by overexpression of kras^v12 or xmrk in transgenic zebrafish (Li et al., 2012; Nguyen et al., 2011; Nguyen et al., 2012), which also showed prominent cell proliferation and liver overgrowth but did not show such obvious upregulation of ribosome genes (data not shown). Thus, the upregulation of ribosome genes in the Myc transgenic zebrafish model should be a specific event elicited by the Myc pathway and thus suggests extra-ribosomal functions for these ribosomal proteins. Extra-ribosomal functions have been found in both prokaryotes and eukaryotes and include regulation of their own transcription and translation, DNA repair endonuclease activity and modulation of gene-specific transcription and translation, such as sequestering MDM, protecting p53 and sequestering their activator Myc (Lindström, 2009; Ruggero and Pandolfi, 2003). Consistent with these reports, we observed obvious upregulation of rpl11, which can bind to Myc and inhibit its function. In addition, rpl26 and rps7, whose homologs have been shown to inhibit MDM2 and to protect p53, were also upregulated. Moreover, among the significantly upregulated ribosome genes, rpl12, rpl13, rpl36a and rps8 have also been reported to be upregulated in human HCC (Kondoh et al., 2001).

Although liver tumorigenesis progressed with time in the TO(Myc) line, we observed mostly hepatocellular adenoma after 4 months of treatment, despite the fact that the same mouse Myc gene expressed in lymphocytes causes leukemia in transgenic zebrafish (Langenau et al., 2005; Langenau et al., 2003). By contrast, transgenic expression of kras^v12 or xmrk alone is sufficient to induce HCC in transgenic zebrafish (Li et al., 2012; Nguyen et al., 2011; Nguyen et al., 2012), suggesting that these oncogenes function distinctively in different tissues and possess different potencies to drive liver tumorigenesis. Similar to our model, Sawey et al. reported that an embryonic hepatoblast cell line lacking P53 and overexpressing Myc was not tumorigenic in vivo in the mouse model. By introducing other mutations or transgenes, these cells provide a sensitized background where a single additional oncogenic factor can trigger tumorigenesis (Sawey et al., 2011). In contrast to the histopathological diagnosis, our comparative transcriptomic analyses showed that the Myc-induced zebrafish liver tumors captured the molecular signatures of advanced HCCs, suggesting that upregulation of Myc signatures are often correlated with malignancy status. In accordance with this speculation, a global gene expression analysis showed that the activation of the MYC transcription signature was strongly associated with malignant conversion of liver cancer (Felsher and Bishop, 1999; Jain et al., 2002; Kaposi-Novak et al., 2009; Shachaf et al., 2004). Thus, our TO(Myc) line could serve as a good driving line for hepatocarcinogenesis. Upregulated Myc can affect multiple hallmark capabilities, such as proliferative signaling, energy metabolism, angiogenesis, invasion and survival (Hanahan and Weinberg, 2011). By crossing the Myc transgenic line with other transgenic or mutant lines, the progress of hepatocarcinogenesis could be accelerated, as we observed in a p53^−/− background (our unpublished data). It could also provide a sensitized background in which a single additional deregulated oncogene can trigger tumorigenesis, thus enabling new oncogenes to be uncovered and studied in the economical and potentially high-throughput zebrafish model. Recently, myc transgenic medaka with liver hyperplasia have also been reported (Menescal et al., 2012) and consistent observations were produced from the two complementing transgenic fish models. It is also interesting to note that Myc-induced mouse HCCs could be rapidly regressed after inactivation of Myc expression (Felsher, 2010; Jain et al., 2002; Shachaf et al., 2004). We also observed a rapid liver tumor regression from our TO(Myc) transgenic zebrafish after removal of Dox (data not shown), similar to earlier reported liver tumor regression from our TO(xmrk) transgenic zebrafish (Li et al., 2012).

Zebrafish were maintained in compliance with Institutional Animal Care and Use Committee (IACUC) guidelines. Transgenic founders (F0) were generated by co-injection into one-cell stage embryos of three plasmid constructs: pFABP10-rtTA2s-M2 [a 2-kb zebrafish fabp10 promoter (Her et al., 2003) in prtTA2s-M2 (ClonTech)], pTRE-Myc [full-length mouse Myc cDNA in pTRE2 (ClonTech)] and pKRT4-GFP (Gong et al., 2002) (GFP marker for transgenic screening). Injected embryos were raised to adulthood and F0 founders were crossed with wild-type fish for transgenic screening. PCR was conducted on GFP-positive embryos to confirm the co-segregation of these three plasmids. A positive F0 fish was used to establish a stable Tg(fabp10:TA; TRE:Myc; krt4:GFP) transgenic line, known as TO(Myc). We always observed a co-segregation of the three plasmids and standard Mendelian inheritance ratios for a single transgenic locus in subsequent generations.

Dox treatment was conducted in six-well plates for embryos and in 6-l tanks for juveniles and adults, with different concentrations of Dox (Sigma-Aldrich). Water was changed every other day. Fish density was ∼60 juveniles or ∼25 adults per tank. Fish samples were fixed in either Bouin's fixative or formalin solution (Sigma-Aldrich). Sections (5 μm) were processed using the Reichert-Jung 2030 BIOCUT Microtome and stained with Mayer's H&E.

Liver tumors (M+D+) and three non-tumor control livers (M-D-, M-D+ and M-D+) were prepared for RNA extraction. Four to five livers were pooled for each RNA extraction. Total RNA was extracted using TRIzol Reagent (Invitrogen) and treated with DNase I (Invitrogen) to remove genomic DNA contamination. Total RNA samples were then submitted to Mission Biotech, Taiwan, for SOLiD-based SAGE sequencing. Briefly, mRNA was purified using Dynabeads Oligo(dT) EcoP (Invitrogen) and subjected to cDNA synthesis. Resultant cDNA was digested by NlaIII and then sequencing adapters were added to the cDNA fragments.

The tags were mapped to the zebrafish Reference Sequence database with a criterion of allowing maximum two nucleotide mismatches. Uniquely mapped tag counts for each transcript were normalized to TPM to facilitate comparison between different samples. To eliminate transcripts with only marginal expression, the expression level of the transcripts had to be above 10 TPM in at least one group of samples (control or treated) in order to be included in subsequent bioinformatics analyses. Because sequence tags for vitellogenins were highly variable in different samples, these vitellogenin tags were excluded in all subsequent transcriptomic analyses. To facilitate functional implications of zebrafish transcriptome, all zebrafish genes were mapped to annotated human and mouse genes in order to use existing online software developed for human genes. Thus, Unigene annotation of zebrafish transcript entries were retrieved from the UniGene database and human and mouse homology mapping of zebrafish Unigene clusters were retrieved from the HomoloGene database. For Unigene clusters mapped by more than one transcript entries, the highest TPM value was used to represent the expression level of the Unigene cluster (van Ruissen et al., 2005). Functional characterization of Unigene clusters was based on Gene Ontology classification and can be obtained from Stanford's SOURCE database (Diehn et al., 2003).

Gene ontology enrichment analysis was performed using DAVID (Database for Annotation, Visualization and Integrated Discovery) (Huang et al., 2009) with the total zebrafish genome information as the background and P-values representing a modified Fisher's exact t-test. Gene Ontology Fat categories were used for this analysis. The cutoff for P-values was 0.01. We used GSEA to establish the relatedness between zebrafish and mammalian liver cancers (Subramanian et al., 2005). Mammalian liver cancer transcriptome data were retrieved from the Gene Expression Omnibus (GEO) database. These data include different stages of human HCCs (GSE12443, GSE6764), human HCAs (GSE11819, GSE7473) and several mouse HCC models (GSE1897) (Kaposi-Novak et al., 2009; Lee et al., 2004; Rebouissou et al., 2009; Rebouissou et al., 2007; Wurmbach et al., 2007). The statistical significance of the normalized enrichment score (NES) was estimated by using an empirical phenotype-based permutation test procedure. A false discovery rate (FDR) value was provided by introducing adjustment of multiple hypothesis testing.

The RNA sequencing data reported in this study have been submitted to the Gene Expression Omnibus with access number GSE40745.

FUNDING: This work was supported by grants from the Biomedical Research Council of Singapore and National Medical Research Council of Singapore [R154000547511].