Multiomics Integration Reveals the Landscape of Prometastasis Metabolism in Hepatocellular Carcinoma*

The systematic investigation of gene mutation and expression is important to discover novel biomarkers and therapeutic targets in cancers. Here, we integrated genomics, transcriptomics, proteomics, and metabolomics to analyze three hepatocellular carcinoma (HCC) cell lines with differential metastatic potentials. The results revealed the profile of the prometastasis metabolism potentially associated with HCC metastasis. The multiomic analysis identified 12 genes with variations at multiple levels from three metabolic pathways, including glycolysis, starch, and sucrose metabolism, and glutathione metabolism. Furthermore, uridine diphosphate (UDP)-glucose pyrophosphorylase 2 (UGP2), was observed to be persistently up-regulated with increased metastatic potential. UGP2 overexpression promoted cell migration and invasion and enhanced glycogenesis in vitro. The role of UGP2 in metastasis was further confirmed using a tumor xenograft mouse model. Taken together, the compendium of multiomic data provides valuable insights in understanding the roles of shifted cellular metabolism in HCC metastasis.

Hepatocellular carcinoma (HCC) 1 is the second leading cause of cancer-related death in the world (1). Liver transplan-tation and resection are believed to be the best approaches to treat HCC and to deliver long-term survival for HCC patients (2,3). However, these methods are not applicable for advanced HCC, and the overall prognosis for HCC remains poor mostly due to tumor metastasis and recurrence (4). Metastasis, the spread of tumor from its primary site to other parts of the body, defines the switch between benign tumor and malignant cancer. Metastasis is estimated to be responsible for ϳ90% of cancer-associated deaths (5). To explore the underlying mechanisms of HCC metastasis is crucial for developing novel and more efficient HCC treatments, especially for advanced HCCs.
The reprogramming of cellular metabolism has been recognized as a key hallmark of cancer (6). Since Otto Warburg first reported that some cancer cells use the glycolysis pathway for energy production even in the presence of oxygen, it has become more and more clear that cancer cells rewire the metabolic fluxes to cope with various microenvironmental situations in order to sustain proliferation and invasion (7). Increasing evidence suggests that metabolism is also a major driver for cancer metastasis (8,9). For example, a glycolytic enzyme phosphoglucoseisomerase has long been known as the autocrine motility factor that promotes tumor cell migration and invasion (10). The glycolytic end-product lactate is reported to be positively associated with metastasis in many types of cancer (11). Lipid metabolism has also been implicated in tumor metastasis (12). A recent study shows that blocking lipid synthesis can overcome tumor metastasis after antiangiogenic therapy (13). However, the detailed metabolism shift associated with metastasis is still largely unclear. As for HCC, most studies have focused on the investigation of glucose metabolism (14). For instance, the up-regulation of several enzymes in the glycolysis pathway, including pyruvate kinase M2 (PKM2), glucose transporters (GLUTs), lactate dehydrogenase, etc., have been reported to be associated with HCC progression and poor prognosis (15)(16)(17). However, less is known about the other metabolic pathways. More importantly, the link between metabolism and HCC metastasis is still missing.
The metabolic reprogramming in cancer is a highly complicated process that requires the coordination of diverse intertwined metabolic pathways. These pathways form a dynamic network that is regulated by multiple levels of gene expression. Therefore, a large-scale and comprehensive analysis of cancer cell metabolism is required to understand the mechanisms and functional consequences of metabolic alterations associated with metastasis. In this study, we integrated data of genomics, transcriptomics, proteomics, and metabolomics from three HCC cell lines, including a low-metastatic cell line, Huh7; a medium-metastatic cell line, MHCC97L; and a highly metastatic cell line, HCCLM3, to mine potential genes and pathways contributing to HCC metastasis. Based on the multiomic analysis and functional study, UDP-glucose pyrophosphorylase 2 (UGP2), an enzyme critical for glycogen synthesis, was found to play an essential role in promoting HCC cell migration and tumor metastasis. Overall, our study described a systematic view of the cellular metabolism associated with HCC metastasis, providing valuable information for developing novel prognostic tools and therapeutic strategies for HCC.
Whole-exome Sequencing and Variation Calling-Genomic DNA was extracted from Huh7, MHCC97L, and HCCLM3 cells using a TIANamp Genomic DNA Kit (Tiangen, Beijing, China) according to the manufacturer instructions. Deep-coverage exome sequencing for HCC cell lines of Huh7 (111ϫ), MHCC97L (131ϫ), and HCCLM3 (124ϫ) were performed at Shanghai Biotechnology Corporation (Shanghai, China) using the illumina 2500 platform (2 ϫ 125 bp). Adapters and low-quality sequences were cleaned by using Trimmomatics (18). Cleaned reads were mapped to the reference genome (GRCh38) with Burrows-Wheeler Aligner (BWA) (19). On average, 99.9% of the exon positions in the reference genome were covered by the studied samples. We then removed duplicated reads and sorted remaining reads with SAMtools (20).
VarScan 2 was used to call candidate single-nucleotide polymorphisms (SNPs) (21). The putative SNPs in cell lines were determined using the additional steps: (1) the SNPs showing read depth Ͻ 20 and phred score Ͻ 15 were removed. (2) We required that a SNP at a certain nucleotide position should be supported by two or more reads on forward strand and two or more reads on reverse strand, as well as be supported by at least 50% of the reads covering this position. Finally, the specific SNPs in MHCC97L and HCCLM3 were identified as these distinguishing from the Huh7 cell line. We removed SNPs from the Single Nucleotide Polymorphism database (dbSNP147).
A sliding window approach (CNV-seq) (22) was adopted by using the 100 kb windows sliding in 50 kb increments to estimate read number (coverage) of each window in the three cell lines. The coverage of each window was normalized by the mean coverage of windows for each cell line. Copy numbers of MHCC97L and HCCLM3 were estimated by comparing the normalized coverage of the paired window between an MHCC97L cell and Huh7 cell, as well as between an HCCLM3 cell and Huh7 cell.
Microarray Analysis-Total RNAs from Huh7, MHCC97L, and HCCLM3 cells were isolated using Trizol reagent (Life Technologies, Carlsbad, CA, USA). A hybridization-based microarray assay was performed at Shanghai Biotechnology Corporation using the Human lncRNA expression microarray (4 ϫ 180K, Agilent). We used three biological replicates for each cell line. Over 20,000 coding genes were covered by Agilent probes in each cell sample. The raw data of nine samples were normalized using the R package limma (quantile algorithm).
Proteolysis and Triplex Dimethylation Isotopic Labeling-Cells were lysed with 8 M urea, and the protein concentration was measured using BCA assay. The samples were reduced by incubating with 10 mM DTT at 37°C for 1 h. The reduced proteins were alkylated for 1 h in darkness with 40 mM iodoacetamide. The alkylation reaction was quenched by adding DTT to a final concentration of 50 mM. The urea in the solution was exchanged to 50 mM sodium bicarbonate buffer by centrifugation using 3 kDa ultrafiltration devices (Millipore). The samples were incubated with trypsin at 37°C overnight for the digestion to complete. For triplex dimethylation isotopic labeling, sodium cyanoborohydride was added to the protein digest for a final concentration of 50 mM, and the deuterated sodium borocyanohydride was used for the heavy labeled samples. Samples were incubated with 0.2 mM formaldehyde, deuterated-formaldehyde or C13-labeled deuterated-formaldehyde, respectively, at 37°C for 1 h. The reaction was quenched with 2 M NH 4 OH, and the samples were mixed and separated using high pH RP-HPLC.
Proteomics-Tryptic digested samples were injected onto an HPLC system (Waters, Milford, MA, USA) coupled with a high pH stable C18 column (Phenomenex Gemini C18, 150 ϫ 2.1 mm, 3 m) at a flow rate of 150 l/min. The peptides were eluted with a 40-min gradient 5-45% buffer B (Buffer A: 50 mM ammonium formate, pH 10; Buffer B: acetonitrile). Fractions were collected every 3 min for 60 min. The tryptic peptide samples eluted from the first-dimensional HPLC were desalted using C18 ziptip and loaded on a nanoUPLC system (Waters) equipped with a self-packed C18 column (C18, 150 ϫ 0.075 mm, 1.7 m). The peptides were eluted using a 5-40% B gradient (0.1% formic acid in acetonitrile) over 90 min into a nanoelectrospray ionization Q Exactive mass spectrometer (ThermoFisher Scientific). The mass spectrometer was operated in data-dependent mode in which an initial Fourier transform (FT) scan recorded the mass range of m/z 350 -1,500. The spray voltage was set between 1.8 and 2.0 kV, and the mass resolution used for MS scan was 70,000. The dynamic exclusion was set to 45 s. The top 20 most intense masses were selected for higher-energy collision dissociation fragmentations. MS/MS spectra were acquired starting at m/z 200 with a resolution of 17,500. The automatic gain control (AGC) target value and maximum injection time were set as 1e6 and 100 ms, respectively, for MS scans, as well as 5e4 and 110 ms for MS/MS scans.
Raw data were searched against the Uniprot (release December 2016) human protein database containing 129,499 sequence entries using the SEQUEST database search algorithm embedded in the Protein Discoverer 1.4 Software (ThermoFisher Scientific). The following parameters were applied during the database search: 10 ppm precursor and fragment mass error tolerance, static modifications of carbamidomethylation for all cysteine residues, dimethylation for the formaldehyde labeling (⌬28 Da), deuterated-formaldehyde labeling (⌬32 Da), or C13-labeled deuterated-formaldehyde labeling (⌬36 Da) on lysines and the N terminus, flexible modification of oxidation modifications for methionine residues, and one missed cleavage site of trypsin were allowed. To determine the confidence of identification, false discovery rate was calculated by searching a decoy database generated by reversing all the protein sequences, and false discovery rate Ͻ0.01 was used as filtering criterion for all identified peptides. In addition, proteins identified with two or more peptides were considered, and proteins identified with the same set of peptides were grouped and treated as one. Quantification analysis was conducted using the Protein Discoverer 1.4 software. Ratios between the three HCC cell lysates were calculated based on the extracted chromatography areas and normalized using the median.
Metabolomics-To analyze the metabolites in the three HCC cell lines, we employed a quantitative polar metabolomics profiling platform by using selected reaction monitoring with a 5500 QTRAP hybrid triple quadrupole mass spectrometer (AB/SCIEX, Framingham, MA) that covered all major metabolic pathways by using a protocol reported by Yuan et al. (23). Briefly, we extracted the metabolites from 10 7 cells with 80% (v/v) methanol (cooled to Ϫ80°C). The metabolites were separated with hydrophilic interaction liquid chromatography (3.5 m; 4.6 mm inner diameter ϫ 100 mm length; Waters) and detected with positive/negative ion switching to analyze 287 metabolites (315 Q1/Q3 transitions) from a single 15-min LC-MS acquisition with a 3-ms dwell time and a 1.55 s duty cycle time. Once the selected reaction monitoring data were acquired, peaks were integrated to generate chromatographic peak areas used for quantification across the sample set by using MultiQuant 2.0 (AB/SCIEX).
Raw Data Processing-For the raw data of transcriptome, proteome, and metabolome, the redundant data were merged, and then each sample was scaled and centered. We marked outliers with normalized values larger than three. There were three, six, and three replicates in each cell line for transcriptome, proteome, and metabolome. Here, we allowed one, three, and two abnormal or missed replicates for each cell line in the transcriptome, proteome, and metabolome, respectively. The mRNAs, proteins, and metabolites that did not fit the criteria in replicate counts were removed. All the outliers and missed replicates were supplied by multiple interpolation method using R package mice.
Differential Modules Analysis-To investigate the association between metastatic capabilities and the levels of mRNA, protein, and metabolite, weighted correlation network analysis (WGCNA) (24) was applied for finding clusters (modules) of highly correlated genes or metabolites. We quantitated the metastatic capability as 1, 2, 3 for Huh7, MHCC97L, and HCCLM3 (trait 1) or as 1, 2, 2 for Huh7, MHCC97L, and HCCLM3 (trait 2). Soft threshold powers were set as 10, 10, 20 for the interpolated data of transcriptome, proteome, and metabolome. Here, we kept modules that significantly correlated to traits (r Ͼ 0.8 or r Ͻ -0.8; p Ͻ 0.05), and then kept mRNAs, proteins or metabolites that significantly correlated to the traits in these modules (r Ͼ 0.8 or r Ͻ -0.8; p Ͻ 0.05).
Quantitative Reverse Transcription PCR (qRT-PCR)-Total RNA was transcripted to complementary DNA using a FastQuant RT kit (TianGen) according to the manufacturer's protocols. Quantitative mRNA expression analysis was performed on a 7500 Fast Real-Time PCR System (ABI, Foster City, CA, USA) using the SuperReal SYBR Green PreMix (TianGen) following the manufacturer's protocols. The mean Ct for each sample was normalized to 18s-rRNA as the reference gene (for primer sequences, see Supplemental Table  S1).
Western Blotting-Protein concentration of each cell extracts was assayed using BCA protein assay kit (ThermoFisher Scientific), and the same amount of proteins was separated by SDS-PAGE and transferred onto the polyvinylidene fluoride membranes using a wet electro-blotter. The membranes were incubated with primary antibodies at 4°C overnight, followed by incubation with secondary antibodies at room temperature for 1 h. Bound antibodies were detected by the enhanced chemiluminescence (ECL) immunoblotting detection reagent.
Plasmid and Transfection-DNA with the complete coding sequence of UGP2 was amplified by PCR using the Premix TaqDNA Polymerase (TaKaRa, Otsu, Japan). The flanking NheI and BamHI restriction sites were created, and the UGP2 DNA was cloned in the pCDH-GFP lentivector expression vector (System Biosciences, Palo Alto, CA, USA). The UGP2 construct was transfected into cells using Lipofectamine 2000 (Life Technologies, Paisley, Scotland). The overexpression efficiency of UGP2 was measured by qRT-PCR and Western blotting. The primer sequences used for cloning the full-length UGP2 are listed in Supplemental Table S2.
Wound-healing, Migration, and Invasion Assays-In a wound-healing assay, 5 ϫ 10 5 cells/well were seeded in six-well plate, allowed to grow for 24 h to 90 -100% confluence, and starved overnight. A scratch was created through the confluent monolayer using a sterile pipette tip. The floating cells were removed with serum-free medium. Then, the cells were cultured with medium containing reduced fetal bovine serum (FBS, 2%) for another 24 h. The remaining width of the scratch was recorded from five randomly selected fields.
Migration and invasion assays were performed using 24-well transwell chamber filters (Millicell Hanging Cell Culture Insert, polyethylene terephthalate 8.0 m, Millipore). For the invasion assays, the membrane was prepared with Matrigel (BD Biosciences, San Jose, CA, USA) following the manufacturer's protocols. After starvation overnight, 1 ϫ 10 5 cells in 200 l of serum-free medium were added to the upper chamber for incubation at 37°C. Next, 600 l DMEM with 10% FBS were added to the lower chamber. Then, nonmigrated or noninvaded cells on the upper membrane surface were removed with a cotton swab, and the migrated and invasive cells on the lower membrane surface were fixed, stained with 0.01% crystal violet solution for 10 min, imaged, and counted in five random 200ϫ microscopic fields.
Glycogen Content Measurement-The glycogen was quantified using a glycogen detection kit (Jiancheng, Nanjing, China) according to the manufacturer's protocol. First, 50 mg of tumor tissues or cells were washed with normal saline, mixed with alkali solution, and boiled for 20 min. Then, the solution was mixed with detection solution, vortexed, and boiled for 5 min before the absorbance was measured at 620 nm wave length. The contents of glycogen were determined by using the standard curve measured at the same time.
In Vivo Metastasis Assay-In vivo metastasis assays were performed using five-week-old male BALB/c-nude mice (Chinese Academy of Sciences, Beijing, China) (26). Briefly, 2 ϫ 10 6 cells were mixed with 20 l of serum-free DMEM and 20 l of Matrigel, then orthotopically inoculated in the left hepatic lobe by a microsyringe through an Experimental Design and Statistical Rationale-For proteomics, three biological replicates and two technical replicates were investigated for each type of cell, resulting in six data points for each cell line for quantification. Statistical and correlation analysis was performed in the same way as the transcriptomics and metabolomics data by WGCNA as described in details in the sections of "Raw Data Processing" and "Differential Modules Analysis." Data Availability-All sequencing data that support this study have been deposited in the NCBI (Bioproject: PRJNA399198). These raw data of cell lines are available under the following experiment accessions: Huh7 (SRX3108375), MHCC97L (SRX3108376), and HCCLM3 (SRX3108377). Microarray data have been deposited in the NCBI (Bioproject: PRJNA382487, GEO accessions: GSE97626). These raw data of cell lines are available under the following experiment accessions: Huh7 (GSM2573327, GSM2573328, GSM2573329), MHCC97L (GSM2573324, GSM2573325, GSM2573326), and HCCLM3 (GSM2573330, GSM2573331, GSM2573332). The proteomic data have been deposited to the ProteomeXchange Consortium via the proteomics identifications database PRIDE (27) partner repository with the dataset identifier PXD005647.

RESULTS
To understand the relationship between cellular metabolism and metastasis in HCC, we employed a multiomic strategy to compare the genome, transcriptome, proteome, and metabolome in three different HCC cell lines with increasing metastatic capabilities, including Huh7, MHCC97L, and HCCLM3 cells. Huh7 is a well differentiated hepatocyte-derived cellular carcinoma cell line that was originally taken from a primary liver tumor, and the MHCC97L and HCCLM3 cell lines were both derived from a metastatic tumor with differential metastatic capabilities (28,29). First, whole-exome sequencing was performed to analyze the SNPs and copy number alterations (CNAs); then, microarray analysis was performed to examine the levels of mRNAs; next, a high-resolution Q Exactive mass spectrometer was employed to quantitatively analyze the protein expressions in the three cell lines; finally, a targeted LC-MS/MS strategy was applied to examine the relative levels of metabolites. The data from genomic, transcriptomic, proteomic, and metabolomic analyses were integrated to describe the landscape of HCC metabolism in the context of metastasis (Supplemental Fig. S1).
Genomic Variations in HCC Cells-First, we profiled the mutations in the three cell lines by whole-exome sequencing analysis, and identified a total of 2,726, 1,562 and 1,374 variants (SNPs and indels) in the coding regions or exon/ intron boundaries for Huh7, MHCC97L, and HCCLM3, respectively. The mutations were enriched in two types of transitions, i.e. A:T-to-G:C and G:C-to-A:T (Fig. 1A). This mutation spectrum was consistent with a previous genomic study using clinical samples of liver cancer (30), as well as the analysis results based on the genomic data of 266 HCC samples from The Cancer Genome Atlas (https://cancergenome.nih.gov/) (Supplemental Fig. S2). The SNP spectra of MHCC97L and HCCLM3 were quite similar due to the same genetic backgrounds between these two cell lines. As compared with the two metastatic cell lines, Huh7 showed a higher percentage of G:C-to-T:A transversion (Fig. 1A). This difference was also confirmed by comparing the mutation spectra between nonmetastatic (M0) and metastatic (M1) liver cancer samples from The Cancer Genome Atlas records (Supplemental Fig. S2).
The majority of the mutations were in the noncoding regions of the genome (Fig. 1B). We here focused on nonsynonymous SNPs (nsSNPs). There were 12 nsSNPs commonly seen in all the three cell lines (Fig. 1B), including NPC1, PRAMEF13, DRB1, MUC3A, COX3, ND4, etc. Compared with the Huh7 cell line, 169 and 65 nsSNPs were specifically observed in MHCC97L and HCCLM3, respectively. These metastatic-cellspecific nsSNPs involved 185 genes. Twenty-one nsSNPs were detected in both of the two metastatic cell lines, and nine of them were predicted to be damaging mutations by PolyPhen-2 (score Ͼ 0.8) (31). These nine nsSNPs were from six genes, including GOLGABA, HDAC6, RIBC2, CDC27, etc. (Supplemental Fig. S3).
Using a sliding window method, the copy number patterns of the two metastatic cell lines were profiled by using Huh7 cells as the control. MHCC97L and HCCLM3 cells displayed very similar relative copy number patterns (r ϭ 0.976, p Ͻ 2.2 ϫ 10 Ϫ16 ; Fig. 1C and 1D). CNAs of the two metastatic cell lines were determined by using 40% variation ratio as the cutoff threshold. Up-regulated CNAs were observed in 19.8% and 18.8% of the MHCC97L and HCCLM3 regions, while down-regulated CNAs were observed in 13.6% and 7.1% of the MHCC97L and HCCLM3 regions, respectively. The abun- dant CNAs suggested that the copy number patterns of MHCC97L and HCCLM3 cells were distinct from Huh7 cells (Fig. 1D). To investigate the correlations between the CNA patterns and metastatic potential, we extracted significant differential CNAs (dCNAs) that continuously up-regulated or down-regulated in the order of Huh7, MHCC97L, and HCCLM3. As a result, 3,321 significant dCNAs were detected, including 2,309 up-regulated dCNAs (orange dots, Fig. 1C), and 1,012 down-regulated dCNAs (blue dots, Fig. 1C). These dCNAs were from 2,449 genes and most frequently occurred in Chr. 3 (21.3%) and Chr. 19 (24.8%).
Next, we investigated the potential biological functions of the 185 genes with metastatic-cell-specific nsSNPs and the 2,449 genes with dCNAs associated with increased metastatic capability. GO analysis showed that the top five enriched biological processes included cellular process, metabolic process, response to stimulus, localization, and biological regulation (Fig. 1E). KEGG pathway analysis indicated that these altered genes were involved with a large spectrum of pathways that fell into seven categories, e.g. metabolism, cellular community, signal transduction, and signaling interaction (Fig. 1F). Some of these signaling pathways were well-known for their critical roles in tumorigenesis and cancer progression, such as the MAPK, Wnt, and ErbB signaling pathways. Moreover, three of the enriched pathways were involved with cell metabolism, including purine metabolism, biosynthesis of amino acids, and ether lipid metabolism. The results suggested that the nucleotide and structure variations in metabolic genes may contribute to the change of metastatic potential. Gene Expression-To determine the gene expression changes among the three cell lines, a hybridization-based microarray analysis was performed to investigate the transcriptome. In addition, quantitative mass spectrometry analysis by triplex stable isotopic dimethylation labeling was employed to study the proteome. WGCNA analysis was applied to investigate the gene co-expression. Two significant modules were detected both in the data of transcriptome and proteome, describing the up-regulated or the down-regulated gene expressions along with the metastatic capability ( Fig.  2A). Overall, 9,635 differentially expressed mRNAs were detected by microarray (Fig. 2B, upper panel). For protein analysis, 338 proteins were observed to be up-regulated in the HCC cell lines with increased metastatic capability, and 399 proteins were down-regulated (Fig. 2B, lower panel).
Overall, 179 genes were observed to be up-regulated in metastatic cells at both mRNA and protein levels, and 132 genes were down-regulated (Fig. 2C). The overlapped genes were searched for KEGG pathways, and the top 30 relevant pathways with p Ͻ 0.05 were shown in Fig. 2D. Nineteen of the 30 pathways were involved with cellular metabolism, e.g. glycolysis, pentose phosphate pathway, biosynthesis of amino acids. The transcriptomic and proteomic analyses showed that the differentially expressed genes were significantly associated with cellular metabolism, pointing to a potential link between metabolism and metastasis. We also mapped 43 metabolic genes with genomic and/or expression changes from six major metabolism processes that were most intensively studied in cancer research, i.e. glycolysis, pentose phosphate metabolism, TCA cycle, glutamine metabolism, nucleotide metabolism, and lipid metabolism (Supplemental Fig. S4). Gene expression differences were observed for many of these genes; however, genetic variations were rarely seen, suggesting that the different expressions of these genes were likely caused by epigenetic or transcriptional regulations.
Multiomics Integration-The genomic, transcriptomic, and proteomic analyses have unveiled the genetic alterations and differential gene expressions in metastatic HCC cells. To understand the outcome of the cellular metabolism changes, a quantitative polar metabolomics profiling platform using selected reaction monitoring with a 5500 QTRAP hybrid triple quadrupole mass spectrometer was employed to analyze the quantitative changes of 287 metabolites. WGCNA was adopted to analyze the modules of the metabolites potentially related with metastasis. The levels of 21 metabolites were detected to be up-regulated along with the metastatic capability, and 32 metabolites were decreased (Fig. 3A). The metabolites with different concentrations between the cells were enriched in 19 metabolic pathways, such as TCA-cycle, glycolysis, purine and pyrimidine metabolism, and amino acid metabolism (Supplemental Table S3). Next, we integrated the data of transcriptome, proteome, and metabolome for KEGG pathway analysis to identify the metabolic genes with variations at multiple levels. By integrating the KEGG pathway analysis results of gene expression (Fig. 2D) and metabolomics (Supplemental Table S3), genes from several pathways were observed with variations at multiple levels (Fig. 3B). qRT-PCR was used to examine the mRNA levels of the identified genes in the three studied cell lines, and the results showed that ten genes were up-regulated in metastatic cells, and two genes were down-regulated, consistent with the microarray analysis results (Figs. 3B and 4). Furthermore, the protein expressions of several identified genes were validated by Western blotting, and the results were consistent with the mass spectrometric measurements (Fig. 5). Among the 12 identified genes, GSTO1 was the only one accompanied by CNA with the same trend as the expression change (Fig. 3B). In addition, ALDH3A1, PYGB, and PGD were mutated in single-nucleotide sites in the MHCC97L genome, and PYGB also showed a SNP mutation in the HCCLM3 genome (Fig.  3B). Furthermore, the identified genes were from three different metabolic pathways, including glycolysis, starch and sucrose metabolism, and glutathione metabolism (Fig. 3B). Taken together, the integrated multiomic analysis identified multiple genes and metabolic pathways that are potentially linked to HCC metastasis.
A Glycogen Regulator UGP2 Regulates HCC Cell Migration and Metastasis-To compare the potential impacts of the 12 identified genes on the metastatic capability, we calculated the coefficients of determination (R 2 ) between gene expression levels and the quantified metastatic capability (trait 1, see "Methods") for the 12 genes (Supplemental Fig. S5). UGP2, a key enzyme regulating glycogen biosynthesis, was observed to be positively correlated with the metastatic capability at both of the mRNA (r ϭ 0.954) and protein levels (r ϭ 0.965) with the highest coefficient of determination among these genes. At the same time, the up-regulation of UGP2 in metastatic HCC tumor samples was also confirmed by comparing the M1 and M0 liver cancer samples from the The Cancer Genome Atlas records (Fig. 6A). UGP2 transfers a glucose moiety from glucose-1-phosphate to MgUTP and forms UDPglucose and MgPPi (32). In liver and muscle tissue, UDPglucose is a direct precursor of glycogen. UGP2 also functions as a critical regulator in the starch and sucrose metabolism pathway. The multiomic analysis showed that UGP2 was up-regulated in the metastatic cells, and three related metabolites were also detected with different levels, suggesting that UGP2 might play an important role in promoting metastasis through the regulation of glycogen synthesis.
Despite the recent revival of interest in cancer metabolism, very little is known about the role of glycogen synthesis. To further explore the implication of this interesting observation, we first examined the mRNA levels of UGP2 in seven different HCC cell lines and a normal hepatic cell line LO2. The levels of UGP2 were higher in the HCC cells as compared with the LO2 cells, and its expression correlated with the migratory potentials of the investigated cells (Fig. 6B). Next, we overexpressed UGP2 in HCC cells and investigated its impact on the cell phenotype (Supplemental Fig. S6). In the PLC and HepG2 cells, the overexpression of UGP2 greatly increased cell migration in the wound-healing assay (Fig. 6C). In addition, Transwell-chamber and Matrigel invasion assays showed that UGP2 overexpression significantly enhanced cell migration (Fig. 6d) and invasion (Fig. 6e). Consistent with these observations, enhanced migration and invasion was also observed in the UGP2-overexpressing Huh7 cells (Supplemental Fig.  S7). The effect of UGP2 on cell proliferation was examined by using colony formation and cell counting kit-8 (CCK8) assays, and interestingly, cell proliferation was not affected by UGP2 in the PLC, HepG2, and Huh7 cells (Supplemental Fig. S8), suggesting a specific role of UGP2 in the regulation of cell migration and invasion. Furthermore, increased glycogenesis was observed with UGP2 overexpression (Fig. 6F). The results suggest that UGP2 may promote cell migration and invasion through elevating glycogen production.
To determine the in vivo impact of UGP2 in the regulation of HCC metastasis, we studied the effects of UGP2 in an orthotopic HCC mouse model. UGP2-overexpressing Huh7 cells showed significant increase in the intrahepatic metastasis as compared with the control cells (Fig. 6G). In addition, higher levels of glycogen were detected in the tumors from mice injected with UGP2-overexpressing cells as compared with the mice injected with the corresponding control cells (Fig.  6H). Taken together, these results indicate that higher level of glycogen could be beneficial for HCC tumor cells to survive the nutrient-deprived condition during metastasis, and UGP2 may play an important role in cell migration and metastasis through the regulation of glycogen synthesis in HCC. DISCUSSION The integrative multiomic analysis of three HCC cell lines with different metastatic potentials revealed the link between metabolism and metastasis. We mapped the gene expression patterns of canonical metabolic pathways and also investigated the downstream metabolites. Three metabolic pathways were detected to be altered at multiple levels, including glycolysis, starch and sucrose metabolism, and glutathione metabolism. Several enzymes from the glycolysis pathway were up-regulated in the metastatic cells, such as PKM2, glucose-6-phosphate isomerase (GPI), and PFKP. It has been well-known that cancer cells utilize enhanced glycolysis to produce building blocks for biosynthesis in order to sustain cell proliferation (33). However, the role of glycolysis in metastasis is controversial under different circumstances (34,35). Our results suggest that the up-regulation of glycolysis may be prometastasis in HCC cells, which is consistent with the fact that many of these glycolysis enzymes are overexpressed in advanced HCC tissues and usually correlate with poor prognosis (36). Furthermore, the overexpression of PKM2 has been reported to enhance metastasis both in vitro and in vivo (37). This study further confirmed the role of key glycolysis enzymes in the regulation of HCC metastasis.
It has been long known that carbohydrate metabolism is highly involved in tumorigenesis and cancer progression, and most studies have focused on glucose and the subsequent Warburg effect (6 -8). On the other hand, it has been reported that different carbohydrates may have distinct impacts on cancer (34,35,38). A recent study demonstrated that sucrose intake in mice could lead to increased tumor growth and metastasis of breast cancer as compared with a nonsugar starch diet (38). However, the underlying molecular mechanism is still largely unclear. No prior studies have investigated the involvement of the starch and sucrose metabolism in HCC. In this study, we observed that three genes in the starch and sucrose metabolism pathway were up-regulated in the highly metastatic cell lines, including UGP2, GPI, and PYGB. In addition, sucrose was decreased in the metastatic cell lines even though all the three cell lines were cultured under the same condition, possibly indicating an elevated sucrose metabolism level in metastatic cells.
UGP2 is the hub molecule connecting starch and sucrose metabolism with glycogen synthesis (32). UGP2 exchanges UDP, from uridine triphosphate (UTP), at the carbon-1 position of the glucose 1-phosphate, generating UDP-glucose and releasing pyrophosphate, PPi. UDP-glucose is the activated glucose required for glycogen synthesis. A recent study indicates that UGP2 is up-regulated in gallbladder cancer and can be used as a biomarker for cancer progression and poor prognosis (39). Hypoxia could induce the expression of UGP2 in a hypoxia inducible factor-dependent manner (40). Our study suggested that overexpression of UGP2 resulted in elevated glycogen level and enhanced cell migration in vitro and increased intraheptic metastasis of HCC in vivo. The most intuitive explanation is that UGP2 helps fill up the glycogen reservoir that can be directly catalyzed to provide sufficient energy for cell migration. Glycogen accumulation improves survival of cancer cells in response to impoverishment of nutrients in the microenvironment, especially during metastasis. Taken together, these data pointed out the important role of UGP2 and glycogen synthesis in the regulation of HCC metastasis.
Furthermore, the integrative-omics analysis also revealed the different status of glutathione metabolism among the investigated three cells. Glutathione is one of the major endogenous antioxidants that directly participate in the neutralization of free radicals and reactive oxygen compounds, as well as maintaining exogenous antioxidants such as vitamins C and E in their reduced forms (41). Here, five genes, i.e. GCLC, GCLM, PDG, GSTO1, and TXNDC12, from the glutathione metabolism pathway were observed with different expression levels in the metastatic cells. For example, glutamate-cysteine ligase catalytic subunit GCLC and GCLM, two isoforms of glutamyl cysteine synthetase, catalyzing the first and ratelimiting step of glutathoine synthesis, are both up-regulated in the metastatic cell lines. Overexpression of GCLC has been reported to increase glutathione production and counteract enhanced oxidative stress in human malignant mesothelioma (42). Meanwhile, the inhibition of GCLM impairs tumor initiation of breast cancer stem cells (43). In this study, we detected higher concentrations of glutathione disulfide (GSSG) and NADPϩ by metabolomics (Fig. 4B). In addition, the levels of intracellular ROS gradually increased in order of Huh7, MHCC97L, and HCCLM3 (Supplemental Fig. S9), indicating higher levels of oxidative stress in the metastatic cells. This observation is consistent with the literature that oxidative stress acts as a key driver of the malignant transformation observed in primary tumors that enhances their metastatic potential (44). Taken together, the results suggest that the metastatic HCC cells elevate the glutathione metabolism to protect themselves from the enhanced oxidative stress and to sustain a highly metastatic phenotype.
In conclusion, this integrative multiomic study provides evidences for the correlation between cellular metabolism and metastasis. The comprehensive analysis defines the dysregulations of metabolic genes at different levels that may impact the metastatic traits of HCC cells. Based on the multiomic analysis, a set of genes were revealed for their potential roles in HCC metastasis through the regulation of metabolism. Furthermore, this study provides an example of integrating omic data from different levels for systematic description of biological systems, paving the way toward omics-based clinical management and personalized medicine.

DATA AVAILABILITY
The proteomic data have been deposited to the Proteome-Xchange Consortium via the proteomics identifications database PRIDE (27) partner repository with the dataset identifier PXD005647.