Comparison of genome-scale DNA methylation profiles in hepatocellular carcinoma by viral status

ABSTRACT

Hepatocellular carcinoma (HCC) incidence has steadily increased in the US over the past 30 years. Our understanding of epigenetic regulation in HCC is still limited, especially the impact of hepatitis B virus (HBV) or hepatitis C virus (HCV) infection on aberrant DNA methylation. We performed genome-wide DNA methylation profiling in 33 fresh frozen tumor samples, including 10 HBV-HCC, 13 HCV-HCC, and 10 non-infected (NIV-HCC) using the Illumina HumanMethylation450 BeadChip. Gene expression profiling was also performed using the Illumina whole-genome DASL HT Assay. Biological influences and gene networks of the differentially-methylated (DM) CpG loci were predicted using the Ingenuity Pathway Analysis. Genome-wide methylation analysis identified 7, 26, and 98 DM loci between HBV-HCC vs. HCV-HCC, HBV-HCC vs. NIV-HCC, and HCV-HCC vs. NIV-HCC, respectively, at P < 5 × 10⁻⁵ for each. Overall, the DM loci were highly enriched for enhancers (48%), promoters (37%), or CpG islands and surrounding regions (37%). Most DM loci were hypermethylated in HCV-HCC compared to HBV-HCC or NIV-HCC. The DM loci were associated with a variety of biological functions including Cell Morphology (HBV-HCC vs. NIV-HCC), Cell Death/ Survival (HBV-HCC vs. NIV-HCC), or Cellular Growth and Proliferation (HCV-HCC vs. NIV-HCC). A subset of the DM loci were correlated (either positively or negatively) with their gene expression or associated with alcohol consumption, BMI, cirrhosis, diabetes, and cigarette smoking. Our findings of differential methylation by viral infection lend insights into the potential effects of viral infection on the epigenetic regulation and further the development and progression of HCC.

KEYWORDS:

• Epigenetics
• Genome-wide DNA methylation
• hepatitis B virus
• hepatitis C virus
• hepatocellular carcinoma

Abbreviations

DM	=	differentially-methylated
HBV	=	hepatitis B virus
HCC	=	Hepatocellular carcinoma
HCV	=	hepatitis C virus
IARC	=	International Agency for Research on Cancer
IP	=	Ingenuity Pathway Analysis
NIV	=	non-infected
shelves	=	2 kb regions upstream and downstream of the CpG island shores
shores	=	2 kb regions upstream and downstream of the CpG islands
TSG	=	tumor suppressor genes
TSS	=	transcription start site
TSS200	=	200 bp of a transcription start site
TSS1500	=	1500 bp of a transcription start site
UBC	=	Ubiquitin C

Introduction

Liver cancer is the second leading cause of cancer deaths in the world.¹ In the US, 39,230 new primary liver cancer cases are expected to be diagnosed in 2016, with approximately 3-fourths being hepatocellular carcinoma (HCC).² The development of HCC is a complex and multifactorial process involving accumulation of genetic and epigenetic alterations that occur along with the progression of normal liver to chronic hepatitis, cirrhosis and, ultimately, HCC.^3,4 The major risk factor for HCC is infection with either hepatitis B virus (HBV) or hepatitis C virus (HCV), the prevalence of which varies from country to country.⁵ Both HBV and HCV are categorized as human carcinogens by the International Agency for Research on Cancer (IARC).⁶

Aberrant DNA methylation has been consistently observed in HCC. With the candidate gene approach, studies have found hypermethylation in various tumor suppressor genes (TSG) and hypomethylation in oncogenes in HCC.^7-9 Recent genome-wide DNA methylation profiling studies have noted substantial DNA methylation changes in HCC,^3,4,10,11 but most of these studies were not designed to specifically identify the changes associated with HBV or HCV infection.

While our prior study showed substantial genome-wide differences in DNA methylation between HCC and adjacent normal tissues,¹⁰ we did not examine methylation changes by viral status and other risk factors of HCC. The aim of this study is to evaluate if there are differentially-methylated (DM) loci that are associated with viral infection in HCC and how these differences in epigenetic regulation affect gene expression and relate to biological functions that may influence on HCC development or progression.

Materials and methods

Study enrollment and biospecimen and data collection

This study was approved by the Institutional Review Board of the University of Hawaii at Manoa and all study participants signed a written informed consent. Tumor samples from 33 HCC patients, including 27 from our previous study, ¹⁰ were analyzed for genome-wide DNA methylation. The 33 samples include 10 non-infected (NIV-HCC), 13 HCV-HCC, and 13 NIV-HCC. HCC was diagnosed histologically by percutaneous biopsy, liver biopsy at the time of surgery, or explanted liver at transplant. Patient information was self-reported and includes age, sex, and ethnicity as identified by the patients, a history of diabetes, a history of cancer besides HCC, previous or current tobacco use, and a family history of cancer. Risk factors for HCC including viral hepatitis, alcohol abuse (defined as greater than 2 drinks/daily for at least 10 years). The presence of other chronic liver diseases was obtained from the medical records. HBsAg, anti-HBc, and anti-HCV serologic status was derived from clinical laboratory testing conducted as part of patients' diagnoses and/or prior medical history. All operative procedures were performed at Hawaii Medical Center-East (2009–2011) by a hepatobiliary/liver transplant surgeon (LW). Final tumor diagnoses and grading of adjacent normal tissue were confirmed by a hepatobiliary pathologist. Adipose tissue was removed through blunt-dissection before the tissue specimen were treated RNA later® (Ambion, Inc., Austin, TX) and stored at −80°C until use.

DNA preparation and genome-wide methylation analysis

DNA was extracted from the tissues using the AllPrep DNA/RNA Mini kit (Qiagen Inc., Valencia, CA). About 500 ng of DNA based on NanoDrop-2000 reading (Thermo Scientific, Wilmington, DE) was bisulfite treated using the EZ DNA Methylation kit (Zymo Research, Irvine, CA) according to the manufacturer's recommendations. The Methylation450 BeadChips (HM450) (Illumina, San Diego, CA) were used to analyze the genome-wide DNA methylation profile across 485,512 CpG loci.⁹ GenomeStudio (v.2011.1) Methylation module (v.1.9.0) software was used to generate a β-value after normalization and background subtraction with the internal controls. The methylation level at each CpG locus represented by a β-value was calculated according to the normalized probe fluorescence intensity ratio between methylated and unmethylated signals. β-values vary between 0 (fully unmethylated) and 1 (fully methylated). Any probe's fluorescence intensity that was not significantly different from the background (P > 0.05) was removed from further analysis. In addition, probes that were designed for CpG sites on either the X or Y chromosome were excluded from analysis. The remaining CpG sites analyzed in the study were 472,009.

RNA preparation and genome-wide expression analysis

Total RNA was extracted from tissues using the AllPrep DNA/RNA Mini kit (Qiagen Inc., Valencia, CA). The quality of the RNA was checked by running 1 μl of each RNA sample on the Bioanalyzer (Agilent, Santa Clara, CA) using the RNA 6000 Nano chip. RNA (200 ng) was labeled, hybridized onto expression microarrays, and scanned. Briefly, the Illumina Whole Genome DASL (WG-DASL) method uses biotinylated random nonamer and oligo(dT) primers to convert input RNA to cDNA. Next, the biotinylated cDNA is annealed to DASL Assay Pool (DAP) probe groups for extension and ligation followed by PCR amplification with a biotinylated and a fluorophore-labeled universal primer. Last, the single stranded PCR products are eluted and hybridized onto an Illumina HT-12 BeadChip. Due to tissue availability, only a subset of samples (10 HBV-HCC, 9 HCV-HCC, and 9 NIV-HCC) was analyzed by the WG-DASL Assay.

Pyrosequencing

A subset of DM loci (n=6) found by the Methylation450 BeadChip were selected for technical validation in 7 HBV-HCC, 7 HCV-HCC, and 7 NIV-HCC samples. About 25 ng of bisulfite treated DNA was first amplified by PCR in the regions of interest using the PyroMark PCR kit (Qiagen, Valencia, CA) following the suggested protocol. All assays were custom designed using PyroMark Assay Design 2.0. Human Methylated and Non-methylated DNA Set from Zymo Research was used as controls (data not shown). All assay information is described in Supplementary Table 1. A PyroMark Q24 instrument (Qiagen) was used and the sequencing results were analyzed using the PyroMark Q24 software (Qiagen).

Statistical analysis

Distributions of clinical characteristics between HCC participants with different viral infection status were compared using the chi-square and Mann-Whitney rank tests for categorical and continuous variables, respectively.

Methylation β-values were converted to M-values for statistical analysis.¹² Partek Genomics Suite (Partek Inc., St. Louis, MO) was used for the analysis. A multivariate ANOVA was performed to identify DM loci by viral infection after adjusting for age and scan date (for batch effect removal) since these variables also contributed to the data variance (data not shown). Raw P-values lower than 5×10⁻⁵ were considered statistically significant since no CpGs were significantly different at the Benjamini and Hochberg's false discovery rate (FDR) q < 0.05. Raw P-value < 5×10⁻⁵ was corresponding to FDR q < 0.99, q < 0.70, and q < 0.26 in HBV-HCC vs. HCV-HCC, HBV-HCC vs. NIV-HCC, and HCV-HCC vs. NIV-HCC, respectively.

Unsupervised clustering analysis and Principal Component Analysis (PCA) were performed using the Partek Genomics Suite. PCA produced linear combinations of the variables (DM loci) to generate the axes known as principal components, or PCs. In this study, PC1 and PC2 were shown in the results. The Euclidian distance between the groups of samples was calculated with the average linkage.

To identify DM CpG islands, methylation levels (β-value) were averaged based on hg18 CpG islands downloaded from UCSC Genome Browser. This yielded 26,327 CpG islands. The averaged β-values were converted to M-values. A multivariate ANOVA was performed to identify CpG islands by viral infection after adjusting for age and scan date.

DM loci were classified by CpG island or its neighbors [2 kb regions upstream and downstream of the CpG islands (shores) or 2 kb regions upstream and downstream of the CpG island shores (shelves)] and functional promoters [within 1500 bp of a transcription start site (TSS) (TSS1500), within 200 bp of a TSS (TSS200), 5′ untranslated regions (5′UTR), first exon (1stExon)], and other regions (body, 3′UTR, or intergenic). To identify DM regions associated with gene function, loci were first stratified by functional locations described above, plus enhancers and non-enhancers. Next, methylation levels in all loci located each region were averaged. Mann-Whitney rank tests with P < 0.05 was considered statistically significant.

GenomeStudio was used to process the raw data generated by the Illumina WG -DASL Assay, and the processed data were analyzed in Partek Genomics Suite. Probes with detection P-value less than 0.05 were considered valid, and 14,167 genes were reliably measured on the expression assay. For the correlation analysis, Target IDs (Probe IDfrom HM450) of the DM loci were generated to match the UCSC reference genes. Since a CpG site could be localized in more than one gene (co-localization of the genes), a Target ID could be linked to multiple genes. For gene expression, intensities from the Illumina WG -DASL Assay were averaged if there were more than one probe for a gene. Then, the correlation analysis between methylation and gene expression was performed. If there were more than one probe for a gene, the intensities were averaged (gene level analysis). Spearman correlation with P < 0.05 was considered statistically significant when analyzing the relevance of the association of methylation to expression.

DM loci associated with viral infection were further analyzed in all 33 HCC samples for their associations with HCC risk factors, including alcohol consumption, BMI, cirrhosis, diabetes, gender, hyperlipidemia, and smoking. Mann-Whitney rank test was used for the analysis and P < 0.01 were considered statistically significant.

Ingenuity pathway analysis

Ingenuity Pathway Analysis (IPA) (Qiagen, Redwood City, CA) was used to predict gene networks and biological functions. Genes stored in the IPA Ingenuity Knowledge Base (IPKB) were cataloged based on reported or suggested biological and molecular functions as well as disease involvement.¹³ The identified genes were mapped to networks available in the IPA and then ranked by score. The score (score = −log₁₀[P-value]) computed by IPA is according to the probability of the set of identified genes and a list of biological functions stored in the IPKB. A score greater than 3 represents a corresponding P < 0.0001, indicating that a gene network was not generated by chance alone.¹⁴

The cancer genome atlas data (TCGA)

We used level 1 data from 380 liver hepatocellular carcinoma (LHCC) tumors in TCGA, including the idat files from the Methylation450 BeadChip. The data were normalized with internal controls and background subtraction as we did for our study. Among 380 HCC tumors, 115 had information on viral infection including 44 HBV/HCV-HCC, 57 HBV-HCC, and 14 HCV-HCC. The 71 samples (HBV-HCC or HCV-HCC) were used for biological validation of DM loci.

Results

Clinical and pathological characteristics of patients

Characteristics of study participants are shown in Table 1. The mean ages at HCC diagnosis for HBV-HCC, HCV-HCC, and NIV-HCC were 62.3 (SD ± 2.2), 64.08 (SD ± 8.5), and 72.0 (SD ± 11.6) years, respectively. The presence of cirrhosis and diabetes differed between HBV-HCC and NIV-HCC at P = 0.02 and P = 0.026, respectively. Age, alcohol consumption, BMI, smoking status, hyperlipidemia, AJCC stage, and race were not significantly different (P > 0.05).

Table 1. Characteristics of study participants.

				P-value^a
Characteristics	HBV-HCC	HCV-HCC	NIV-HCC	HBV- vs. HCV-HCC	HBV- vs. NIV-HCC	HCV- vs. NIV-HCC
Age, years, mean±SD	62.3 ± 12.2	64.1 ± 8.5	72.0 ± 11.6	0.852	0.105	0.107
BMI (kg/m2), mean±SD	24.6 ± 5.1	24.8 ± 5.5	26.5 ± 6.3	0.951	0.631	0.369
Gender, n (%)				0.400	0.400	1.000
Male	7 (70)	11 (85)	7 (30)
Female	3 (30)	2 (15)	3 70)
Cirrhosis, n (%)				0.163	0.250	0.020
Yes	1 (12.5)	5 (42)	0 (0)
No	7 (87.5)	7 (58)	10 (100)
Missing	2	1	0
Alcohol consumption^b, n (%)				0.231	0.305	0.054
Yes	1 (10)	4 (31)	0 (0)
No	9 (90)	9 (69)	10 (100)
Smoking status, n (%)				0.382	0.371	0.072
Ever smoker	6 (60)	10 (77)	4 (40)
Never smoker	4 (40)	3 (23)	6 (60)
Hyperlipidemia, n (%)				0.277	0.708	0.162
Yes	3 (23)	1 (8)	3 (30)
No	10 (77)	12 (92)	7 (70)
AJCC stage, n (%)				0.282	0.214	0.292
1	5 (50)	8 (62)	9 (90)
2	1 (10)	2 (15)	0 (0)
3a	2 (20)	1 (8)	1 (10)
3b	2 (20)	0 (0)	0 (0)
4	0 (0)	2 (15)	0 (0)
Diabetes, n (%)				0.183	0.371	0.026
Yes	4 (40)	2 (15)	6 (60)
No	6 (60)	11 (85)	4 (40)
Race				0.196	0.196	0.270
Asian	9 (90)	9 (69)	7 (70)
Black	0 (0)	1 (8)	0 (0)
Pacificist	1 (10)	0 (0)	2 (20)
White	0 (0)	3 (23)	1(10)

^aMann–Whitney rank test for continuous variables, Chi-squre for categorical variables.

^bif they drank more than 2 drinks/day for 10 years total.

DM loci in HBV-HCC vs. HCV-HCC

Seven DM loci localized in 6 genes and one intergenic locus (cg03994465) were identified between HBV-HCC and HCV-HCC (Fig. 1, Supplementary Table 3). Both PCA and heatmap showed clear separations by these DM loci (Fig. 1A and Fig. 1B), respectively. The loci were approximately 2.2-fold more methylated in HCV-HCC (average β-value: 0.236) relative to HBV-HCC (average β-value: 0.513) (P = 5.57×10⁻⁵) (Fig. 1C).

Figure 1. Differentially methylated CpGs between HBV-HCC vs. HCV-HCC. (A) Principal component analysis of the methylation data was plotted using the first 2 principal components (PC1 = 67.2% and PC2 = 9.04%). Each dot represents a sample (red for HBV-HCC and blue for HCV-HCC). (B) Unsupervised hierarchical clustering of β-values for differentially methylated loci. Red and blue blocks on the top of the maps represent HBV-HCC (n = 10) and HCV-HCC (n = 13), respectively. (C) Mean methylation level for all 7 differentially-methylated loci among HBV-HCC (left) and HCV-HCC (right). (D) The top IPA network involving differentially methylated genes. Red and blue genes indicate differentially methylated genes and the connected cancer related genes, respectively. (E) The top associated cellular functions with numbers of differentially-methylated genes.

Among 7 DM loci, 5 of those localized in genes (BMP8A, C1orf194, FLNC, GFl1, and KIAA1324) were associated with promoter regions, covering TSS1500, TSS200, 5′UTR, and 1stExon (Supplementary Table 3). Six of the DM loci (85.7%) were located in CpG islands and 4 (57%) in Enhancers (Supplementary Table 3). The DM locus in the BMP8A promoter was in a CpG island, and it had the most significant difference in methylation, Δβ = 0.36 (P = 8.42×10⁻⁶) (Supplementary Table 3).

Gene network analysis showed that “Organ Morphology, Skeletal and Muscular System Development and Function, and Connective Tissue Disorder” are the top functions or disease (Score = 8) (Fig. 1D). This network was connected with Ubiqutin C (UBC), and 88% of genes in this network were related to cancer development, including HCC (15 HCC related genes; blue colored in Fig. 1D). The top associated cellular function was cell morphology (associated with BMP8A, FLNC, GIF1, KIAA1324, and SMYD3) (Fig. 1E).

DM loci in HBV-HCC vs. NIV-HCC

We identified 26 DM loci between HBV-HCC and NIV-HCC (Supplementary Table 4). The two types of HCC showed a clear difference in the methylation pattern of those DM loci (Figs. 2A and 2B). All except one locus located in the promoter of THBS1 were significantly hypermethylated in HBV-HCC (range of β-value: 0.66–0.80) relative to NIV-HCC (range of β-value: 0.34–0.66) at P = 1.08×10⁻⁵ (Fig. 2C). Of the 26 DM loci, 14 (54%) were located in Enhancer regions (Supplementary Table 4). Among 25 hypermethylated loci in HBV-HCV compared to NIV-HCC, only one locus (TFAP2A) was located in a CpG Island (Supplementary Table 4); while 8 loci were in the surrounding regions of a CpG island (Shores or Shelves). Of the DM loci, 17 (65%) were associated with 17 unique genes that were involved in the top IPA network of “Cellular Development, Cellular Growth and Proliferation, and Embryonic Development” (Score = 19) (Fig. 2D). This network was connected by extracellular regulated kinase (ERK). In this network, 68% of the genes played roles in cancer development including HCC (14 HCC related genes, blue colored in Fig. 2D). The DM loci were enriched for the genes having roles in Cell Death and Survival, including ATF5, CPEB4, CSF2RB, HOXA3, MMP13, NUP62, TFAP2A, and THBS1 (Fig. 2E).

Figure 2. Differentially methylated CpGs between HBV-HCC vs. NIV-HCC. (A) Principal component analysis of the methylation data was plotted using the first 2 principal components (PC1 = 68.4% and PC2 = 6.55%). Each dot represents a sample (red for HBV-HCC and green for NIV-HCC). (B) Unsupervised hierarchical clustering of β-values for differentially methylated loci. Red and green blocks on the top of the maps represent HBV-HCC (n = 10) and NIV-HCC (n = 10), respectively. (C) Mean methylation level for all 7 differentially-methylated loci among HBV-HCC (left) and NIV-HCC (right). (D) The top IPA network involving differentially methylated genes. Red and blue genes indicate differentially methylated genes and the connected cancer related genes, respectively. (E) The top associated cellular functions with numbers of differentially-methylated genes.

DM loci in HCV-HCC vs. NIV-HCC

Ninety-eight DM loci were identified between HCV-HCC and NIV-HCC (Supplementary Table 5). This number (98 DM) was much higher than that identified between HBV-HCC and NIV-HCC (26 DM) (Fig. 2). The 98 DM loci clearly separated HCV-HCC from NIV-HCC (Figs. 3A and 3B). All these loci were hypermethylated in HCV-HCC (range of β-value: 0.66–0.78) relative to NIV-HCC (range of β-value: 0.32–0.68) at P = 1.21×10⁻⁴ (Figs. 3B and 3C).

Figure 3. Differentially methylated CpGs between HCV-HCC vs. NIV-HCC. (A) Principal component analysis of the methylation data was plotted using the first 2 principal components (PC1 = 64% and PC2 = 6.65%). Each dot represents a sample (blue for HCV-HCC and green for NIV-HCC). (B) Unsupervised hierarchical clustering of β-values for differentially methylated loci. Blue and green blocks on the top of the maps represent HCV-HCC (n = 13) and NIV-HCC (n = 10), respectively. (C) Mean methylation level for all 7 differentially-methylated loci among HCV-HCC (left) and NIV-HCC (right). (D) The top IPA network involving differentially methylated genes. Red and blue genes indicate differentially methylated genes and connected cancer related genes, respectively. (E) The top associated cellular functions with numbers of differentially-methylated genes.

Of the 98 DM loci, 45 (46%) were located in Enhancers, and 34 (35%) were in promoter regions (TSS1500, TSS200, 5′UTR, and 1stExon) (Supplementary Table 5). Sixty-eight of the DM loci (70%) were localized in 60 genes and the rest (30%) were in the intergenic regions (Supplementary Table 5). The DM genes were enriched for transcription regulators (20%, 13/60), including DNAJB6, HEXIM1, HIF3A, HOXA3, HOXA6, HR, NCOR1, PLAGL1, SHOX2, TBX15, TBX4, and TEAD2, and enzymes (17%, 10/60), including CES3, DHRS3, FBXL13, GAD1, SMYD3, ST5, TRAP1, TRIM71, TTLL1, and VARS2. Gene network analysis indicated that the DM genes were involved in “Embryonic Development, Organismal Development, and Cardiovascular Disease” (Score=43) (Fig. 3D), and were associated with Cellular Growth and Proliferation (AXIN2, CDHR2, DISC1, DNAJB6, GAD1, HEXIM1, HOXA3, HR, IGFBP4, ITGB4, MST1R, NCOR1, PLAGL1, PRKCDBP, RARA, ROR1, SHOX2, SMYD3, ST5, STEAP2, TEAD2, TFAP2A, TRAP1, TRIM71, and TSPO) (Fig. 3E).

Correlation between DM loci and gene expression

Spearman correlation analysis was performed to investigate correlations of the DM loci with expression of their corresponding genes. The analysis showed that 15 DM loci were significantly correlated with gene expression, 60% of which were located in the enhancer regions (Table 2). Two highly correlated CpG island DM loci between HBV-HCC and HCV-HCC were in inverse correlation (high methylation and low expression) and were located in the gene body of MFSD7 (r = −0.41, P = 0.031) and promoter of FLNC (r = −0.43, P = 0.022) (Table 2). Also, 2 DM loci that were hypermethylated in HBV-HCC compared to NIV-HCC were significantly correlated with gene expression. One was in the S_Shore of the ATF5 promoter and correlated with low expression of the gene (r = −0.64, P = 0.00025). Another located in the 5′UTR (not related to CpG island) of KIAA1199 was also positively correlated with its expression (r = 0.40, P = 0.036) (Table 2). In addition, we found 11 DM loci between HCV-HCC and NIV-HCC that were significantly correlated with the expression of 10 paired genes either negatively (CES3, DHRS3, IGFBP4, LRRC17, and SGIP1) or positively (AXIN2, HOXA6, RAB11FIP4, ROR1, and STEAP2) at P < 0.05 (Table 2).

Table 2. Correlations of DM loci and gene expression.

Comparison	Gene name (DASL)	TargetID (HM450)	Relation to UCSC CpG Island	Functional Location	Enhancer	Spearman, r	Spearman, P-value
HBV vs HCV	MFSD7	cg05827190	Island	Body	Enhancer	-0.41	0.0315
HBV vs HCV	FLNC	cg02661879	Island	5'UTR;1stExon	Non-enhancer	-0.43	0.0221
HBV vs NIV	ATF5	cg04315771	S_Shore	TSS1500;5'UTR	Non-enhancer	-0.64	0.0003
HBV vs NIV	KIAA1199	cg02585724	Opensea	5'UTR	Enhancer	0.40	0.0357
HCV vs NIV	AXIN2	cg09231741	Island	Body	Non-enhancer	0.45	0.0156
HCV vs NIV	CES3	cg26538442	Opensea	TSS200	Enhancer	-0.48	0.0099
HCV vs NIV	DHRS3	cg26029221	N_Shelf	Body	Enhancer	-0.51	0.0054
HCV vs NIV	HOXA6	cg23129930	Island	1stExon	Non-enhancer	0.45	0.0161
HCV vs NIV	IGFBP4	cg02627216	N_Shore	TSS1500	Non-enhancer	-0.44	0.0199
HCV vs NIV	LRRC17	cg20066782	Opensea	Body	Enhancer	-0.45	0.0169
HCV vs NIV	RAB11FIP4	cg07146073	Opensea	Body	Enhancer	0.51	0.0055
HCV vs NIV	ROR1	cg03937361	Opensea	Body	Enhancer	0.67	0.0001
HCV vs NIV	ROR1	cg13942627	Opensea	Body	Enhancer	0.44	0.0193
HCV vs NIV	SGIP1	cg09665332	Opensea	Body	Enhancer	-0.52	0.0046
HCV vs NIV	STEAP2	cg06112923	Opensea	Body;3'UTR	Non-enhancer	0.51	0.0057

DM loci associated with clinical characteristics of HCC

The DM loci from above analyses were further analyzed for their associations with clinical characteristics, including alcohol use, BMI, cirrhosis, diabetes, gender, hyperlipidemia, and smoking in all 33 patients. The most significant results are shown in Fig. 4, and a detailed list of correlations is provided in Supplementary Table 6. The DM loci in MMP13, FMNL3, and MST1R were associated with drinking and smoking habits, and history of cirrhosis, with higher methylation in ever drinkers compared to never drinkers (average β-value 0.95 vs. 0.75), in cirrhosis vs. non-cirrhosis (average β-value 0.81 vs. 0.67), and ever smokers vs. never smokers (average β-value 0.82 vs. 0.72), respectively. A significantly lower methylation in an interegenic locus (enhancer linked region, cg06817454) was found in diabetes relative to non-diabetes (average β-value 0.56 vs. 0.70). Another intergenic locus (cg05412459) located in the S_Shelf was significantly correlated with BMI (r = −0.55, P = 0.0009). None of the DM loci were correlated with gender or hyperlipidemia. A heatmap of methylation levels of all identified DMs by different viral infection as well as clinical characteristics was shown in all 33 samples (Fig. 5). Lower methylation levels of DM loci among NIV-HCC (more consistently blue) were found compared to the ones in HBV-HCC and HCV-HCC (more consistently red). Although overall the unsupervised clustering analysis showed well separated methylation profiles between HBV-HCC and HCV HCC, these groups still had quite similar methylation patterns regardless of clinical characteristics (Fig. 5).

Figure 4. Association between differentially methylated loci by virus infection and other HCC risk factors. Mann-Whitney rank test was used for alcohol consumption (ever vs. never), cirrhosis (yes vs. no), diabetes (yes vs. no), and smoking (yes vs. no) and Spearman correlation was used for BMI (continuous variable).

Figure 5. Unsupervised hierarchical clustering of β-values for all differentially methylated loci identified by different viral infection. Unsupervised hierarchical clustering of β-values for DM loci (rows) in 33 HCC samples (columns). Red, blue, green blocks on the top of the maps represent HBV-HCC (n = 10), HCV-HCC (n = 13) and NIV-HCC (n = 10), respectively. The clinical characteristics of the samples including alcohol use, BMI, cirrhosis, diabetes, and smoking were represented as other blocks below sample group.

Pyrosequencing analysis

To validate the HM450 results, we selected 6 DM loci for pyrosequencing analysis based on their P-values and availability for assay design. The loci include cg15931894 (CRIM1) and cg02585724 (KIAA1199) identified in HBV-HCC vs. NIV-HCC, and cg00727675 (HOXA3) and cg02149965 (VARS2) identified in HCV-HCC vs. NIV-HCC (Supplementary Table 1). Due to the limitation in assay design for top significant DM loci between HBV-HCC vs. HCV-HCC, we had to select 2 alternate DM loci at raw P-values lower than 1×10⁻⁴ for technical validation, cg00201414 (THSD1P) and cg13170600 (Supplementary Table 1). High concordance of methylation levels for all DM loci analyzed was found between the microarray and pyrosequencing analyses (Spearman correlation r = 0.94, P = 6.5 × 10⁻⁶⁰) (Supplementary Fig. 1).

Biological validation of DM loci using TCGA samples

Characteristics of study patients between our study and TCGA are shown in Supplementary Table 2. Most of the patients in TCGA were HBV-HCC (80%) compared to our samples (43%) at P = 0.0007. Race was significantly different for HBV-HCC and HCV-HCC, P = 0.0002 and P = 0.0001, respectively. A majority of the participants in our samples were Asian (78%) while 2 thirds of the patients in TCGA were White (66%). Stages were significantly different for HBV-HCC (P < 0.0001). All HBV-HCC in our samples was the early stage (1 or 2) while 30% of HBV-HCC in TCGAs was the late stage (3 or 4). Age and gender were not different (P > 0.05).

The seven DM loci identified between HBV-HCC and HCV-HCC in UHCC were also analyzed in the TCGA samples (Supplementary Fig. 2). A locus in BMP8A (cg02575697) was significantly methylated in HCV-HCC compared to HBV-HCC (P = 0.004). A borderline significance was observed for another locus co-localized in KIAA1324 and C1orf194 (cg09820792) (P = 0.06). Although other DM loci did not show statistical difference, we found the expected direction of methylation differences in all loci except 2 loci in an intergenic (cg03994465) and a body of SMYD3 (cg08106377).

Discussion

DNA methylation has been shown as significantly different between HCC and adjacent normal tissue.^3,4,10,11 HBV and HCV infection are the 2 main risk factors responsible for HCC development in humans,⁵ but our understanding of viral infection-associated epigenetic changes in HCC is still limited.^15,16 In this study, we investigated genome-wide methylation profiles in HCC in association with viral infection, and found many differentially methylated loci between virally infected and non-infected tumors as well as by different types of infection (HBV vs. HCV). We also discovered that most of the DM loci were in the enhancer regions, CpG islands or gene promoters suggesting a role in the regulation of gene expression (Supplementary Table 7). Indeed, a subset of the DM loci was significantly correlated with the expression of their corresponding genes. We also identified some of the DM loci associated with alcohol use, BMI, cirrhosis, diabetes, and smoking. Biological gene network analysis showed that the majority of the DM genes were related to cancer development including that of HCC and were associated with Cell Morphology, Cell Death/Survival, and Cellular Growth and Proliferation.

HBV and HCV are 2 different types of viruses (DNA vs. RNA), indicating that their molecular mechanisms underlying the liver tumorigenesis may be different.¹⁷ HBV may work either through viral DNA insertion into the human genome to activate oncogenes and/or induce genetic instability, or through immune responses to the viral infection which lead to fibrosis and cirrhosis.^18,19 HCV, however, is a single-stranded RNA virus the life cycle of which is thought to be predominantly cytoplasmic,²⁰ and the tumorigenic effects may be mediated through the alteration of cell signaling and metabolism as well as via increased immune responses.²¹

In the current study, we observed more hypermethylated DM loci in HCV-HCC and HBV-HCC as compared to NIV-HCC, particularly in HCV-HCC. This finding is supported by several previous studies that focused on tumor suppressor genes. Deng et al. showed hypermethylation of 15 differentially methylated genes in HCV-HCC relative to HBV-HCC.²² Similarly, Yang et al. found three tumor suppressor genes (APC, p15, SOCS1) to be more frequently methylated in HCV-HCC compared to HBV-HCC or NIV-HCC.²³ Moreover, Nishida et al. reported large abundance and higher levels of methylation in HCV-HCC compared to normal liver tissues and NIV-HCC, respectively.²⁴ Thus, evidence supports more prominent methylation changes in HCV-HCC, and our study identified several novel loci with the changes in comparison to HBC-HCC or NIV-HCC.

Altered methylation in regulatory regions such as enhancers or regions surrounding the transcription start sites is known to be important in regulating gene expression.²⁵ Understanding how altered methylation relates to gene regulation can be complex because it has been shown that promoter methylation regulates gene expression while gene body methylation affects alternative splicing, noncoding RNA expression, and disease-causing mutations.^26,27 Our studies showed approximately 2.5-fold higher frequency of DM loci in enhancers (46–57%) compared to total loci analyzed in enhancers (21%), P < 0.05 (Supplementary Table 7), indicating potential importance of these DM loci in regulation of distant genes. Indeed, most of the enhancer-related loci were demonstrated to be correlated with the expression of their associated genes in our study. While a majority of the DM loci in HBV-HCC vs. HCV-HCC were located in a CpG island, DM loci identified in HBV-HCC vs. NIV-HCC or in HCV-HCC vs. NIV-HCC were more likely to be in the surrounding regions of a CpG island (<4 kb from CpG island) or not related to CpG island (>4 kb from CpG island), respectively. These findings were consistent with a separate regional analysis of CpG islands based on UCSC Genome Browser to identify DM CpG island by different viral infection (Supplementary Table 8), where 17 DM CpG islands were identified between HBV-HCC vs. HCV-HCC, while 10 or none were found between HCV-HCC vs. NIV-HCC or HBV-HCC vs. NIV-HCC, respectively. This observation indicates that additional studies are needed to clarify whether different regulatory mechanisms exist for these DM loci with regard to their epigenetic regulation and if any, how this difference in regulation contributes to HCC development or progression in viral HCC.

In this study, gene network analyses revealed potential biological differences between HBV-HCC, HCV-HCC, and NIV-HCC. The top network using the DM genes identified in HBV-HCC vs. HCV-HCC was “Cell Morphology, Cell death/Survival, and Cellular Growth and Proliferation,” which was connected with Ubiquitin C (UBC). UBC was first described as a stress responsive gene ²⁸ that regulates various functions of a cell (e.g., protein degradation, DNA repair, cell cycle regulation, kinase modification, endocytosis, and regulation of other cell signaling pathways) (http://www.ncbi.nlm.nih.gov/gene/7316). Through an in silico network analysis based on protein-protein interaction, a recent study suggests potentially important roles of the UBC network in the pathogenesis of HCC.²⁹ Another network, “Cellular Development, Cellular Growth and Proliferation, Embryonic Development” was indicated by the DM loci identified in HBV-HCC vs. NIV-HCC. This network was connected by ERK which is a well-known oncogene frequently activated in many cancers including HCC.^30-34 Importantly, it has been known that hepatitis Bx, a transactivation protein to stimulate HBV gene expression, promotes angiogenesis by stimulating the ERK expression during HBV infection.¹⁷ The last top network, “Embryonic Development, Organismal Development, Cardiovascular Disease,” predicted through the DM loci in HCV-HCC vs. NIV-HCC, was connected by 5 major genes including AKT, ERK1/2, NFkB, VEGF, and Histone H3. First four of these genes have been well studied for their roles in invasion and metastasis and are suggested to be therapeutic targets for HCC.^35,36 Histone H3 is involved with the structure of nucleosomes, and is an important protein in epigenetic regulation of gene expression.³⁷ A recent study by Quan et al. suggests that HCV core-induced epigenetic silencing of SFRP1 may lead to the activation of the Wnt signaling pathway and the downregulation of SFRP1 was due to the loss of binding capacity of acetyl-Histone H3 in the SFRP1 promoter,³⁸ indicating an important role of HCV infection in gene regulation through Histone3.

Given that our findings of significant DM loci were based on locus-by-locus analysis, we questioned if there was significant differences in functional regions by viral infection. A separate analysis was performed to test for DM regions based on functional locations including enhancers. In this analysis, methylation levels of CpG sites were averaged by each region (Supplementary Fig. 3). We found significantly higher methylation in HCV-HCC compared to NIV-HCC only for TSS200 (P < 0.041) and for enhancer (P < 0.038). Although we did not see significant differences for those regions between HBV-HCC and HCV-HCC or between HBV-HCC and NIV-HCC, the methylation levels of TSS200 in HBV-HCC was 3.9% lower than HCV-HCC, but 2.4% higher than NIV-HCC. Similarly, the methylation levels of HCV-HCC in Enhancer were 7% and 4.4% higher than those of HBV-HCC and NIV-HCC, respectively. Taken together, the results of regional analysis were consistent with those of locus-locus analysis, showing more hypermethylation in HCV-HCC and HBV-HCC compared to NIV-HCC.

This study has several features that distinguish it from previous studies. The present study employed a genome-wide approach to examine DNA methylation differences by viral status in HCC. In contrast, many previous studies focused on selected tumor suppressor genes to evaluate the differences between HCV-HCC and HBV-HCC. We identified additional potentially novel DM loci by viral status using the genome-wide approach. In addition, we explored regional analyses by CpG islands, functional location as well as enhancers. We also evaluated the biological implications of the DM loci by assessing their correlation with gene expression and analyzing their regulatory networks. As strength of the study, we were able to validate a subset of DM loci identified between HBV-HCC and HCV-HCC using an independent dataset (TCGA) even though the differences were less significant in TCGA.

It is important to note some limitations of this study. First, the statistical power of our analyses is limited by a relatively small sample size. Therefore, there may be additional DM loci by viral infection that we were unable to identify. Another limitation is that the patients were primarily Asian, and so our finding may not be generalizable to other populations. However, we had fairly equal numbers of HBV-HCC, HCV-HCC, and NIV-HCC, which may provide valuable insight into different epigenetic alterations by viral infection within a population with high incidence of both viral hepatitis and HCC, while many previous studies focused only on HBV-HCC or HCV-HCC depending on the predominant viral incidence in that specific population.

In summary, we report differentially methylated CpG loci by viral infection in HCC using genome-wide methylation profiling. In spite of the limited sample size, we found clear differences in methylation coupled with altered gene expression among HBV-HCC, HCV-HCC, and NIV-HCC, suggesting viral infection to be an important contributor to the epigenetic changes seen in HCC.

Disclosure of potential conflicts of interest

No potential conflicts of interest were disclosed.

Acknowledgments

The authors thank the late Mr. Mike Loomis for building and maintaining the clinical database and Mr. Kyle Miyazaki for collecting clinical data. We thank the Pathology Shared Resource for the collection and processing of tumor tissue samples and the Genomics Shared Resource for performing Illumina HumanMethylation450 Beadchip assay at University of Hawaii Cancer Center. We thank Drs. Janos Molnar and Karolina Peplowska for performing the pyrosequencing.

Funding

This work was supported by the following grant: NIH 3P30CA071789-12S6.