Transposon mutagenesis identifies genes driving hepatocellular carcinoma in a chronic hepatitis B mouse model

HBV-associated HCC mouse model for mutagenesis screening

Liver-SB/HBsAg transgenic mice develop chronic liver inflammation with associated reactive hyperplasia (Figure 1a, Supplementary Figure 1), hepatocytomegaly, and ground glass hepatocytes (Figure 1b) similar to the human chronic hepatitis B and to the previously characterized HBsAg mice ¹⁸. Liver-SB/HBsAg modifications lead to the appearance of preneoplastic foci, visible from 19.7 weeks of age, followed by hepatocellular adenoma (HCA) and trabecular HCC (Supplementary Figure 2). SB induced a cooperative tumorigenic effect with HBsAg as Liver-SB/HBsAg mice displayed reduced survival (Figure 1c). We also noted a tendency to larger numbers of tumors (Supplementary Figure 2) and to more advanced stage disease compared to HBsAg mice alone (Figure 1d).

To identify the genes mutated by SB that cooperate with HBsAg-associated inflammation in tumor induction, we PCR amplified and sequenced the transposon insertion sites from 250 tumors harvested from 34 mice¹⁹, which yielded 328,687 sequence reads and identified an average of 1315 unique transposon insertion sites per tumor. By comparing the location of the transposon insertions in all tumors we found a few tumors that were genetically related. These tumors were removed leaving 228 genetically unrelated tumors, which were subsequently used for downstream analysis (Supplementary Figures 3,4).

Large-scale screen identifies Common Insertion Sites (CISs)

We then screened for CISs, which are regions in the cancer genomes that contain a higher density of transposon insertions than predicted by chance and therefore are likely to contain a cancer gene (CISs on chromosome 1 could not be conclusively identified because of local transposon hopping and CISs on this chromosome were therefore excluded). A gene-centric CIS-calling method (gCIS) that looks for a higher density of transposon insertions within the coding regions of all RefSeq genes than predicted by chance²⁰ identified 2525 gCIS genes (Supplementary Table 1a), while a non-gene-centric Gaussian Kernel Convolution (GKC) method²¹ that looks for a higher density of insertions within fixed kernel widths of 15Kb to 240Kb²² identified 2041 CISs containing 2103 genes (Supplementary Table 1b). There was a 83% overlap between the CISs identified by these two methods (p<2×10⁻¹⁶, Chi-square with Yates correction) and together, they identified 2871 CIS loci containing 2881 genes (Supplementary Table 1c, Supplementary Figure 5a). In vitro transposition cell culture assays have suggested that SB is a random insertional mutagen and the only requirement for insertion is a TA dinucleotide¹². Although the GKC method scans the entire cancer genome for CISs, most CISs were located within or in close proximity to genes, providing additional evidence that SB targets the coding regions of genes that confer a selective advantage for tumor growth. gCISs not identified by GKC often contained very large (>300Kb) or small genes (<10Kb) and were probably missed in part because of the use of fixed kernel widths (Supplementary Figure 5b,c).

To estimate the genetic coverage of this screen we randomly selected subgroups of 10 to 220 tumors from the total of 228 tumors and then used GKC to identify the CISs for each subgroup. The number of tumors in each subgroup was then plotted against the genomic base pairs covered by all the CISs identified for each subgroup. The curve plateaued before reaching 100 tumors, indicating that the screen was approaching saturation with as few as a 100 tumors (Figure 2a,b). Adding more tumors merely identified additional CISs of lower frequency (Supplementary Figure 6,7). Consistent with this notion of saturation, this screen identified most of the CIS genes found previously in much smaller-scale HCC transposon screens performed in p53 mutant¹⁷ and Sav1-deficient mice (see below) (Supplementary Figure 8). To our knowledge this is the first saturating transposon screen for cancer genes reported for mice.

Nearly 68% of the CIS genes were mutated in <20% of the tumors and on average each tumor contained 387 mutated CIS genes (Figure 2c,d). This is indicative of extensive intratumor heterogeneity, likely resulting from branching tumor evolution, similar to that reported for several human cancers²³, combined with the large mutational load induced by SB. In this HCC model, every liver cell contains ~350 copies of the transposon¹³, which when induced to transpose, can transpose over and over again. In addition, because the transposon insertion sites are PCR amplified before being deep-sequenced, insertional mutations present in a small number of cells in the tumor can be identified. It is this ability to sample the extensive intratumor heterogeneity at such a deep level that we believe made it possible to approach saturation with only a 100 tumors.

General features of HCC CIS genes

The transposon¹³ contains a promoter used to deregulate oncogenes and stop cassettes in both transcriptional orientations to inactivate TS genes. Transposon insertions at oncogenes therefore tend to be located at the 5′ end of the gene in the same transcriptional orientation, while insertions in TS genes are usually located throughout the gene in either orientation. Insertion site patterns are thus often predictive of gene function. For the genes identified here and in other transposon screens performed in solid tumors^{13-15,17,19,22,24}, the majority are predicted to function as TS genes. Most genes also have insertions only in one allele, suggesting they may function as haploinsufficient TS genes, although it remains possible that the other allele is silenced by an epigenetic or other mutational event. This is consistent with recent studies performed in human cells and tumors that identified hundreds if not thousands of haploinsufficient genes, which are able to drive cellular proliferation. These genes are postulated to operate in cooperative networks to maximize proliferative fitness²⁵, and support other studies indicating that even partial inactivation of TS genes can contribute to tumorigenesis⁸.

Branching tumor evolution complicates efforts to implement personalized medicine and suggests that targeted therapies might be directed to genes mutated at the trunk of the evolutionary tree. Transposons provide powerful tools for identifying trunk genes since insertions in trunk genes are more likely to be associated with a higher number of sequencing reads. Despite the splinkerette PCR method preventing a perfect indication of read quantitation, previous studies show that CIS genes with high read depth are more likely to be highly penetrant, as seen with APC²². To identify potential trunk genes for HCC we again screened for CISs but this time using only the 25,197 SB insertions that were represented by at least four sequencing reads each. Among the 524 CISs identified (Supplementary Table 1d) were 21 outlier CISs that showed the highest sequencing read counts and frequency of occurrence (Figure 2e, Supplementary Figure 9a). Six genes including Arid1a^10,11, Gsk3b²⁶, Iqgap2²⁷, Magi1²⁸, Pten²⁹ and Sav1^30,31 have known tumor suppressive roles in HCC, while Snd1 is described as an oncogene promoting HCC angiogenesis^32,33. Two genes including Zbtb20^34,35 and Ankrd17^36,37 are involved in hepatocyte differentiation and maturation, while Zbtb20 is a transcription factor that represses the expression of α-fetoprotein, a widely used biomarker used for HCC surveillance. The Setd2 TS gene³⁸ and the putative TS genes Adk³⁹, Dpyd^40,41, Mll5⁴² and Nfia^43,44 were also identified as potential driver genes for HCC, although there is no published evidence they are involved in HCC. Six driver genes also regulate hepatic metabolism with Pten⁴⁵, Gsk3b⁴⁶, Adk^47,48, Zbtb20⁴⁹ and Ghr^50,51 regulating lipid and glucose metabolism, and Dpyd controling pyrimidine catabolism in hepatocytes. Seven driver genes including Mll5, Setd2, Wac, Arid1a, Nfia, Snd1 and Zbtb20 are transcription factors, cofactors or chromatin remodelers, highlighting the importance of the modulation of transcriptional programs in HCC. One driver CIS we named Rtl1 Chr12 locus (Figure 2e) is the merger of several CISs at the same imprinted locus⁵² (Rian, Rtl1, 6430526N21Rik genes and a microRNA cluster) and was targeted in 38.4% of the tumors. The majority (92%) of the transposon insertions at this CIS are located on the negative strand (Supplementary Figure 9b), indicating a potential oncogenic loss of imprinting and transcriptional activation of the Rtl1 gene, or one or more of the microRNAs located downstream. Interestingly, these microRNAs are among the most strongly up-regulated microRNAs in murine and human HCC⁵³.

Comparisons of CIS genes with human data and other indications

We next asked whether the 2881 HCC CIS genes are specific for HCC by comparing them to the genes identified in transposon screens for colorectal cancer (CRC)²² and pancreatic adenocarcinoma (PDAC)^19,24. These comparisons showed that >60% of the HCC CIS genes are also mutated in CRC and PDAC (p<2.10⁻¹⁶) (Supplementary Figure 10). This is not surprising since many cancer genes are known to function in more than one cancer cell type. HCC CIS genes are also significantly mutated in human HCC as shown by the very significant overlap with the genes somatically mutated in human HCC identified in two whole-genome HCC sequencing studies^10,11, and in a genome-wide analysis of HBV integrations⁵⁴ (Table 1). In addition, we found a significant enrichment among the genes listed in the Catalogue of Somatic Mutations in Cancer (COSMIC) database, which reports all of the somatic mutations identified in human HCC, and the genes listed in the Cancer Gene Census database, which aims to catalogue all of the driver genes in human cancers (Table 1). HCC CIS genes thus play important roles in human HCC and provide additional validation for the genes identified in human HCC.

The expression of HCC CIS genes is also specifically deregulated in human HCC. Among 2189 HCC CIS genes contained in a large microarray dataset obtained from 223 HBV-positive HCC patients⁴ (Supplementary Figure 11a), 45% were significantly deregulated in human HCC (absolute fold change >1.5 and adjusted p-value <0.00001) (Supplementary Figure 11b; Supplementary Table 2), compared to 29% misregulated genes on the array that are not CIS genes (p <2.10⁻¹⁶). Differentially expressed HCC CIS genes also displayed significantly lower p-values and more of them were up-regulated than down-regulated, although this difference was not large (Supplementary Figure 11c-d). Comparison of the HCC microarray data with other CIS gene lists from PDAC^19,24 and CRC²² screens showed the highest enrichment for HCC CIS genes (Supplementary Figure 11e), indicating better specificity when the data were from the same tumor type. Overall, 42.5% of the HCC CIS genes are mutated or mis-expressed in human HCC (Figure 3, Supplementary Table 3). The deregulation of 30 HCC CIS genes was subsequently confirmed by real time PCR using a microfluidic dynamic array containing 18 HBV-positive human tumors and 9 Liver-SB/HBsAg tumors (Supplementary Figure 11f; Supplementary Table 4). The direction of the expression change (up- or down-regulated) was positively correlated between human and mouse tumors (Supplementary Figure 11g).

Transposon drive HCC via conserved pathways

HCC CIS genes are also enriched in the major cancer signaling pathways known to be important for human HCC such as the Ras/Erk, p53, Akt, Wnt, and Tgfβ/Bmp pathways, together with the recently identified hepato-tumor suppressor Hippo pathway³⁰ (Figure 4a, Supplementary Table 5, Supplementary Figure 12). Another known HCC signaling pathway, IL6/Stat3, was less enriched, possibly because the tumor microenvironement locally secretes IL6, thus activating this pathway in hepatocytes. A similar analysis for PDAC^19,24 and CRC²² CIS genes showed that many of these signaling pathways are also enriched for mutations in these genes. Interestingly, the Hippo pathway was more highly mutated in HCC, whereas the Wnt pathway, which is critical for CRC²², was more highly mutated in CRC (Supplementary Figure 12a). Transposons also prominently targeted the Hippo pathway in tumors generated in a Sav1-deficient background (Methods and Supplementary Table 1e). Sav1 activates MST1, which results in the downstream activation of Hippo signalling⁵⁵. Together, these results suggest that transposon-induced mutations in Hippo pathway genes cooperate with Sav1-deficiency to further deregulate Hippo signaling, and are consistent with the results from a transposon screen performed in Apc mutant mice, which identified a large number of CRC CIS genes that targeted the Wnt pathway²². Liver-SB/Sav1-deficient mice also developed multiple large liver tumors nearly two months earlier than Liver-SB/HBsAg mice, providing additional evidence for an oncogenic role for deregulated Hippo signaling in HCC.

HCC CIS genes target metabolic events

When we classified the HCC CIS genes according to their gene ontology we noticed a striking over-representation of genes linked to cellular metabolic processes (p=1.4 × 10⁻⁶⁰, 1199 genes) (Figure 4b, Supplementary Table 6). A less striking enrichment for metabolic processes was detected for CRC and PDAC CIS genes (Supplementary Figure 13), while most other ontologies such as transcription or cell cycle displayed similar enrichments among all three tumor types. When we focused the gene classification on metabolic processes, we found 47 metabolic categories that were targeted by the transposon during HCC tumor development, against only 9 and 10 metabolic categories for CRC and PDAC CIS genes, respectively (Supplementary Figure 14). When we looked for enrichment for disease genes annotated in Ingenuity Pathway Analysis, we also found a highly significant enrichment for HCC CIS genes in metabolic diseases (265 genes, p = 4 × 10⁻⁶); genetic association to cancer was also highly evident (431 genes, p = 6 × 10⁻¹⁰). Overall, genes associated with protein, carbohydrate, lipid and nucleic acid metabolism were highly mutated by SB. Several other processes were also targeted including oxidoreduction and metabolism of ATP, organophosphate, hormones, isoprenoid and vitamins (Supplementary Table 6).

To determine whether similar genetic changes occur in human HCC we measured metabolic gene expression levels in 18 HBV-positive human HCC and 9 Liver-SB/HBsAg tumors using RT-qPCR. The deregulation of key metabolic genes was positively correlated in these two species, and similar gene expression signatures were found for most genes (Supplementary Figure 15). The classification of HCC CIS genes has thus highlighted the importance and specificity of targeting metabolic genes in HCC. The fact that normal hepatocytes display important general metabolic functions that most likely need to be reorganized during malignant transformation may explain the reprogramming by SB of many aspects of hepatocyte metabolism.

Metabolic profiling and mapping of mouse HCC tumors

To better understand how genetic changes in metabolic genes affect hepatic tumor metabolism we performed non-targeted metabolic profiling on eight food-restricted Liver-SB/HBsAg animals. Chemical derivatisation and gas chromatography coupled to time-of-flight mass spectrometry (GC/TOFMS) were utilized to identify the differences in metabotypes between malignant (tumor masses > 3mm) and adjacent normal liver tissue. The quality of acquired metabolomics data was verified by principal component analysis (PCA) and partial least squares discriminant analysis (PLS-DA) (Supplementary Figure 16). Even with eight matched samples, statistical analyses using the t-test and Welch-test were able to identify 56 annotated metabolites significantly changed in tumors compared to normal adjacent tissues (Supplementary Figure 17). Consistent with the transposon data we noticed significant alterations corresponding to specific metabolic pathways. For instance, carbohydrate metabolism was affected through increased levels of glucose and fructose and reduced ribitol, ribonic acid and allonic acid. Levels of several amino acids including isoleucine, valine, aspargine, tyrosine, methionine, serine, leucine, phenylalanine, threonine and alanine were also significantly downregulated in tumors, consistent with data from human tumors⁵⁶. Moreover, components of nucleic acid including ribose, pyridine, uracil and uridine were found at higher levels in tumors. We then integrated the genetic (Supplementary Table 1c), expression⁴, and metabolic results (Supplementary Table 7) in order to identify potential correlations. We noticed a marked consistency between genetic alterations, mRNA levels and actual amounts of metabolites, revealing a probable control of metabolic changes by genetic modifications at least to glycolysis, the tricarboxylic acid (TCA) cycle, glutaminolysis, the pentose phosphate (PP) pathway (Figure 5) and redox regulation (Supplementary Figure 18). Collectively, these metabolic alterations control bioenergetics, biosynthesis and redox regulation, which are known to optimize tumor cell metabolism for proliferation and survival.

We also tested our interpretation of the integrated data using biochemical assays and in vivo live metabolic imaging ⁵⁷. We found that intrahepatic lactate dehydrogenase (Ldh), alanine transaminase (Alt) and glutaminase (Gls) enzymatic activities were consistently more strongly activated in large hepatic tumors (Supplementary Figure 19), supporting the enhanced conversion of pyruvate to alanine and lactate, and the inclusion of glutamine possibly as an anaplerotic carbon source to supplement the increased metabolic fuel requirement in larger tumors⁵⁸. Metabolic imaging with hyperpolarized pyruvate⁵⁷ also showed a major conversion of pyruvate into alanine and lactate (Figure 6), while its conversion into bicarbonate and other TCA intermediates are unchanged, indicating a shift towards aerobic glycolysis and anabolic processes, relative to oxidative phosphorylation. It is to be noted however that pyruvate is only one of the fuel sources in HCC, and other processes generating ATP such as lipid catabolism have indeed been found to be significantly increased in HCC as well ^56,59. Overall, these experimental findings were consistent with our predictions based on the integration of multiple datasets.

Mice and histology

We utilized the following alleles to generate a mouse model of HBV-associated hepatocellular carcinoma: HBsAg transgenic¹⁸, Alb-Cre⁷⁵, T2/Onc2 (6113)¹³ or T2Onc3 (12740) and Rosa26-lsl-SB11¹⁴. The resulting cohorts were on a mixed B6.129 genetic background. Alb-Cre/+; T2Onc2/+; Rosa26-lsl-SB11/+; HBsAg/+ mice were used for the SB screen. Alb-Cre/+; T2Onc3/+; Rosa26-lsl-SB11/+; HBsAg/+ were used for metabolic imaging. All animals were monitored on a biweekly basis in accordance with IACUC guidelines. Full necropsies were performed. For DNA extraction, liver tumors were measured and snap-frozen. For histology, livers were 10% formalin fixed and paraffin-embedded. 5μm sections were processed for haematoxylin and eosin staining, or orcein staining. All slides were reviewed by our veterinary pathologist, Dr. Jerrold M. Ward. Histological classification of hepatic lesions and tumors were as previously reported^18,76.

Identification of transposon insertions sites

Identification of transposon insertion sites was performed using splinkerette PCR to produce barcoded PCR products that were pooled and sequenced on the 454 GS-Titanium (Roche) platform. Reads from sequenced tumors were mapped to the mouse genome assembly NCBI m37 and merged together to identify non-redundant Sleeping Beauty (SB) insertion sites. Cloning and mapping of the transposon insertion sites were performed as previously described^19,22. 328,687 non-redundant insertion sites (Supplementary Dataset 1) were used to identify CISs using a gene centric method²⁰ and a Gaussian Kernel Convolution (GKC) statistical framework method²¹.

Gaussian Kernel Convolution (GKC) method for common insertion site (CIS) identification

Any insertions sites on the transposon donor mouse chromosome 1 for tumors derived from T2/Onc2 were excluded from common insertion site (CIS) analysis. The likelihood of ‘local hopping’ of the transposon is increased where the transposon array is located⁷⁷. This phenomenon can significantly increase the background level of transposon insertion events, thereby complicating CIS analysis. 328,687 non-redundant insertion sites were used to identify CISs using a Gaussian kernel convolution (GKC) statistical framework²¹. The previous GKC analysis approach²² was enhanced by utilizing multiple kernel scales (widths of 15K, 30K, 50K, 75K, 120K and 240K nucleotides). CISs predicted across multiple scales and overlapping in their genomic locations were clustered together, such that the CIS with the smallest genomic “footprint” was reported as the representative CISs. For highly significant CISs with narrow spatial distributions of insertion sites, the 15K kernel is typically the scale on which CISs are identified. The P value for each CIS was adjusted by chromosome and a cut-off of p<0.05 was used.

Pruning genetically similar tumors from the same animals

We initially sequenced 250 tumors from 34 mice and called CISs using the GKC method on all tumors (with a non-stringent p-value significance threshold to capture many loci). Using the boundaries of the CIS loci, we obtained genomic positions that defined 2110 bins on the genome. We scored every tumor for every bin as having an insertion in the bin or not. We computed the Spearman correlation between all tumors and performed hierarchical clustering (complete linkage). When tumors from the same animal clustered together as pairs in the dendrogram, we retained only the tumor, which had the largest total number of insertions. The aim of the pruning is to ensure that there is at least a single tumor from another animal that is most similar to a tumor from a given animal. In theory, pairs of clustered tumors from the same animal can remain after the pruning step if, for example, three tumors originally clustered together and only one was removed. To resolve this problem the pruning step can be repeated iteratively. After the pruning step, we tested whether the median Spearman correlation between the 228 remaining same-animal-tumors and the median Spearman correlation between the non-same-animal tumors was statistically significant. If this was not the case, we terminated the pruning process. On our dataset a single pruning round was sufficient. We used the remaining 228 tumor samples further to call CISs with both GKC and gCIS methods.

Analysis of screen saturation

In order to determine whether the SB screen had reached saturation, and if so, what number of samples was sufficient to reach saturation, the Gaussian Kernel Convolution (GKC) CIS calls from all 228 samples were analyzed using the ACT software package⁷⁸. ACT considers genomic locations generated by multiple samples for specific biological phenomenon under study (e.g.ChIP-seq peaks) to determine the saturation of a screen. The program considers the various combinations in which samples can be added so that the increase in base pair coverage is a range of values based on all the samples. The results can be depicted as a series of boxplots showing the increase in base pair coverage, where the boxplot at each position n on the x-axis shows the coverage values of all combinations of n samples. Boxplots that approach a horizontal asymptote indicate that the coverage has reached saturation.

For the GKC CISs generated by all 228 samples, the insertion sites that contributed to CISs were extracted, resulting in a set of approximately 131,000 sites. The insertion sites were then selected per sample and pseudo-kernels of 7.5k nucleotides either side of each insertion were applied to mimic GKC kernels of 15k nucleotides. Overlapping kernels within each sample were merged into continuous genomic regions. These 228 modified insertion files were then analyzed using ACT. For each combination of samples the median values, and 25th and 75th percentiles were plotted using the ggplot2 visualization package for the R statistical analysis platform. As a control, the 228 samples were reanalyzed where the same number of insertion sites per sample was selected at random across the mouse genome (excluding chromosome 1, the donor chromosome). The pseudo-15k nucleotide kernels were applied. Figure 2a shows the saturation plot for all 228 samples, clearly indicating a ‘knee’ in the profile as the graph asymptotes with an increasing number of samples. Visual inspection of the graph indicates that below 50 samples, the screen does not appear to have reached saturation. From 100 samples upwards there appears to be sufficient coverage to report saturation. Conversely, Supplementary Figure 4a for the randomized samples is virtually a straight line, indicating no saturation. While the analysis does not produce a clear-cut asymptote this is to be expected due to the type of data under consideration. ACT was designed to analyze such data as ChIP-seq arrays for predicting transcription factor binding sites. In these scenarios ChIP-seq replicates should ideally report the same key binding sites/genomic locations. Hence across multiple samples the same locations should be reported. For SB screens however, while insertions in the same gene will be found from different samples, the locations of the insertion sites will not overlap perfectly, even with the addition of the 15k nucleotide pseudo kernels. Hence each sample will introduce novel regions, such that the overall coverage will continue to increase even if the screen has truly reached a ‘saturation’ point. Also, not all samples will contribute to all CISs. Different combinations of samples will thereby result in varying coverage’s, causing the coverage profile not to asymptote perfectly.

Cross-species oncogenomics analysis

Gene expression assays, conducted on 223 HBV-positive human patients diagnosed with HCC, was used to compare tumors and paired non-tumor tissues. The microarray data are publicly available and downloadable at the Gene Expression Omnibus database with accession number GSE14520⁴. Human hepatocellular carcinoma (HHC) expression profiles were downloaded from the GEO database (GSE14520). Only the 223-paired samples hybridized to the HG-U133A platform were selected. Then, CEL files were Robust Multi-array Average (RMA) normalized and differential gene expression between normal and tumor cells were computed using Affy and limma package, implemented in Bioconductor. As a result, adjusted P-value (FDR-based p-value adjustment according to Benjamini and Hochberg) and fold-change were obtained for each probset. Our mouse candidate cancer genes were converted to human genes using three different databases of homologous genes (MGI, HomoloGene and Ensembl)^79-81. The obtained human genes were then mapped to the list of Affymetrix platform probe sets. For genes detected with more than one probe set, the probe set with best-adjusted p-value was considered as gene representative. The threshold of significance used was the absolute value of fold change more than 1.5, and adjusted p value less than 0.00001.

Smaller scale liver mutagenesis screen in the Sav1 deficient background

This study was performed in a Sav1 deficient background with Alb-Cre/+, T2Onc2/+, Rosa26-lsl-SB11/+, Sav1^fl/fl animals. The Sav1 conditional allele was previously described³⁰. A total of 136 large tumors from 31 animals were utilized for analysis. Most animals were euthanized at 8 to 9 months of age, when they had multiple large tumors. The protocol for genomic DNA extraction, barcoded amplification, and sequencing on Illumina Hi-Seq2000 was previously described²⁰. 6,251,915 reads were mapped to the mouse genome (mm9) and filtered when they were represented by very few sequence reads²⁰. Finally, the gene-centric CIS identification method (gCIS) identified 73 common insertion sites (Supplementary Table 1e).

Patient Samples

Resected tumor samples were obtained from patients undergoing curative resection for HCC (n=172). Operations took place between 1991-2009 at the National Cancer Centre, Singapore. All samples were collected in accordance with the requirements of the local Singapore Ethical Committees and informed consent was obtained from all subjects. The patient demographics and clinical characteristics have already been described in a previous study⁸². Total RNA was isolated using TRIzol (Invitrogen) following the manufacturer’s instructions, and the RNA concentration was quantified by NanoDrop (Thermo Scientific). RNA samples were reverse transcribed to cDNA using the SuperScript III cDNA Synthesis Kit (Invitrogen) in 10 μl reactions containing 1 μg total RNA and oligo(dT) primers according to manufacturer’s instructions.

Gene expression analysis from patient samples

Specific primers targeting genes of interest were designed using Primer Express software version 3.0. Primer sequences used for specific target amplification (STA pre-amplification reaction) as well as qRTPCR are listed in Supplementary Table 4. For the STA reaction, each cDNA sample was pre-amplified with 200nM pooled STA primer mix and Tagman PreAmp Master Mix (Applied Biosystems) in a 5ul reaction, which was run for 14 cycles according to the manufacturer’s protocol. To remove unincorporated primers, each sample was treated with Exonuclease I (Fermentas) following incubation at 37°C for 30 minutes. For inactivation, the mix was in a second step, incubated at 80°C for 15 minutes. At the end of the Exonuclease I treatment, the reactions were diluted 1:5 in TE buffer (pH 8.0) prior to use for qRT-PCR. The Fluidigm BioMark™ real-time PCR system and 48.48 Microfluidic Dynamic Arrays were employed for high-throughput qRT-PCR analysis⁸³. As volume per inlet is 5 μl, the 6 μl volume per inlet with overage was prepared. For the samples, 2.7 μl of each STA and ExoI-treated sample were mixed with 20X DNA Binding Dye Sample Loading Reagent (Fluidigm) and 2X SsoFast EvaGreen SuperMix with Low ROX (Bio-Rad). For the gene expression assays, 0.3 μl of mix primer pairs (100 uM) was added with 2X Assay Loading Reagent (Fluidigm) following the addition of 1X TE buffer to 6 ul volume. Prior to loading the samples and assays into the inlets, the chip was primed in the NanoFlex 4-IFC Controller. The samples and assays were then loaded into the inlets of the dynamic array. Following loading and mixing of the samples and assays into the chip by the IFC Controller, PCR was run with the following reactions conditions: 50°C for 2 minutes, 95°C for 10 minutes, followed by 40 cycles of 95°C for 15 seconds and 60°C for 60 seconds. Global threshold and linear baseline correction were automatically calculated for the entire chip. The melting curve analysis and cycle threshold (Ct) were provided by Fluidigm Real-Time PCR Analysis software version 3. Fold change in expression of genes of interest between liver tumor and adjacent non-tumor samples were determined using the comparative Ct method following the formula: 2^{−ΔCt(tumor)}/2^{−ΔCt(non-tumor) 84}. GUSB and Gusb were used as the internal control gene in human and mouse samples, respectively. The -ΔCt data obtained from this calculation was then used to generate the heatmap as well as supervised hierarchical clustering of tumor and adjacent non-tumor samples by dChip software. In addition, statistical analysis was performed by a standard t-test to identify significant differentially expressed genes (p ≤ 0.05) between two studied groups for each species.

Pathway and gene ontology analyses

Gene-annotation enrichment analyses were completed with DAVID Bioinformatics v6.7⁸⁵. The list of all human genes was used as default background. For gene ontology analyses GOTERM and PANTHER were used. For pathway analyses, KEGG was selected. We obtained enrichment of genes that display a genetic association with cancer through DAVID using GENETIC_ASSOCIATION_DB_DISEASE and OMIM.

Metabolomics experiments

Ex vivo biochemical enzyme activity assays

For each enzyme activity assay, 100 mg of liver tissue was homogenized in 200 μL of ice-cold 100 mM tris HCl buffer, then centrifuged for 10 min at 13,000 × g to remove insoluble material. All assays were based on continuous spectrophotometric rate determination. Hepatic ALT, AST and GLS enzyme activities were measured with commercially available enzymatic assay kits (Biovision, catalog #752-100, #753-100; Sigma-Aldrich, EC 3.5.1.2).

Magnetic resonance spectroscopy (MRS) and metabolic imaging

This technology is partially described in a recent publication ⁵⁷.