Exploring DNA Methylation Profiles Altered in Cryptogenic Hepatocellular Carcinomas by High-Throughput Targeted DNA Methylation Sequencing: A Preliminary Study for Cryptogenic Hepatocellular Carcinoma.

BACKGROUND: Hepatocellular carcinoma (HCC) includes cryptogenic hepatocellular carcinomas (CR-HCC) that lack a defined cause. Specific DNA methylation patterns and comparisons of the aberrant alterations in DNA methylation between CR-HCC and adjacent peritumor tissues (APTs) have not yet been reported.
METHODS: The SureSelectXT Methyl-Seq Target Enrichment System was used to sequence targeted DNA methylation in three paired CR-HCC tissues and APTs. Gene Ontology (GO) enrichment and KEGG pathway analysis were performed to investigate the DNA methylation mechanism of CR-HCC. The mRNA expression levels of HOXB-AS3, HOXB6, HOXB3, USP18, MAP3K6, TIRAP, TNNI2, SHC3, CTTN, and TFAP2A, selected from the identified signaling pathways, were evaluated by quantitative real-time PCR (qPCR).
RESULTS: A total of 1728 differentially methylated regions (DMRs) were identified in tumor tissues compared with non-tumor tissues, of which 868 DMRs were hypermethylated and 860 were hypomethylated. The DMRs were mapped within 2091 DMR-associated genes (DMGs). The mRNA expression of HOXB-AS3, HOXB3, and MAP3K6 was downregulated in CR-HCC tissues compared to the APTs. However, the mRNA expression of TIRAP, SHC3, and CTTN was upregulated in the CR-HCC tissues. Differences between the mRNA expression of HOXB6, USP18, TNNI2, and TFAP2A in the CR-HCC and APTS tissues were not statistically significant. GO analysis showed that the molecular functions of "binding", "protein binding", and "cytoskeletal protein binding" were the main categories for the hypermethylated DMGs. The hypomethylated DMGs were mostly enriched in the molecular functions "binding", "protein binding", "calcium ion binding", among others. KEGG pathway analysis showed that the hypermethylated DMGs were enriched in several pathways such as "estrogen signaling pathway", while hypomethylated DMGs were enriched in several pathways such as "proteoglycans in cancer", suggesting that epigenetic modifications play important roles in the cryptogenic hepatocarcinogenesis.
CONCLUSION: These results provide useful information for future work to characterize the functions of epigenetic mechanisms on CR-HCC.

Introduction

Hepatocellular carcinoma (HCC) is the second leading cause of cancer deaths worldwide.1 Improvement in HCC treatment and prognosis is, therefore, urgent, but is challenging. To our knowledge, various factors, including infection with hepatitis B (HBV) or C viruses (HCV) and alcoholic and non-alcoholic steatohepatitis (NASH),2–4 are primarily responsible for HCC. HCC patients with active HBV or HCV infection have worse prognoses for survival time than uninfected patients, namely non-B non-C HCC (nBnC-HCC).5,6 Although there have been many recent studies on the development of HCC the mechanisms promoting HCC initiation and development remain largely unknown. Thus, it is important to understand the carcinogenic mechanism of HCC and to explore new therapeutic targets for the disease.

We defined cryptic hepatocellular carcinoma (CR-HCC) as HCC patients without any background liver etiology after excluding all other testable liver disease etiologies, including occult HBV infection, NASH, primary biliary cirrhosis, hemochromatosis, Wilson’s disease, alcoholic liver damage, α1 antitrypsin deficiency, severe steatosis, and autoimmune hepatitis. Because there are no guidelines for early detection of CR-HCC, these tumors to be greater in size and at advanced stages when diagnosed compared with hepatitis virus-related HCCs.7 It has been reported that the survival time of advanced CR-HCC is even shorter than that of hepatitis B/C-HCC.8 Therefore, there is an urgent need for a surveillance program for the early detection of CR-HCC. Although studies have identified a few of the tumor characteristics and the postoperative prognosis of CR-HCC patients, the main mechanism underlying CR-HCC is still unclear because of the rarity of the disease.

Epigenetic mechanisms can play important roles as interfaces between the genome and environmental influences, including lifestyle and other risk factors.9,10 Some investigations have suggested that epigenetic multi-regulation causes different epigenetic changes related to liver disease and may account for the progression from liver disease to HCC. A number of studies have investigated novel HCC therapies from the perspective of epigenetic alterations, targeting epigenetic regulation.11,12 Accumulating evidence on genome-wide DNA methylation provides strong support for the significance of differences in methylation of CpG sites between HCC tissues and adjacent peritumor tissues (APTs).13,14 We have consistently observed aberrant DNA methylation in virus-associated HCC.15,16 Some studies have found abnormal DNA methylation in the carcinogenic progression of non-alcoholic fatty liver disease and NASH to HCC.17 One of the key molecular events in the development of NAFLD-derived HCC was found to be hypermethylation in the promoter region of the glycine N-methyltransferase (GNMT) gene resulting in reduced expression.18

In this study, we focused on cryptogenic HCC (CR-HCC), one of the HCCs of unknown etiology, investigating differences in DNA methylation between CR-HCC tissues and APTs with the aim of increasing our understanding of the role of epigenetic mechanisms in the carcinogenesis of CR-HCC.

Methods

Sample Collection and Preparation

Cryptogenic HCC patients were identified from patients undergoing hepatectomy in the First Affiliated Hospital of Anhui Medical University from January 2017 to May 2018. Patients with HBV or HCV infections or alcoholism-related HCC, or with chronic liver disease of known etiology, such as autoimmune hepatitis, severe steatosis, and α1 antitrypsin deficiency, were excluded A total of three CR-HCC patients were identified.

Specimens from these three patients were collected during the operation. The APTs were collected from a location at least 3 cm away from the tumor boundaries, and all specimens were stored at −80°C in liquid nitrogen. All the collections followed the same protocol. All of the cancerous tissues were diagnosed as primary hepatocellular carcinoma by two independent and experienced pathologists.

Library Preparation and Data Analysis

Agarose gels (1%) were used to monitor the genomic DNA degradation and contamination. The Nano Photometer® spectrophotometer (IMPLEN, CA, USA) was used to check DNA purity. The Qubit® DNA Assay Kit in the Qubit® 2.0 Fluorimeter (Life Technologies, CA, USA) was used to measure DNA concentrations. According to manufacturer’s protocols, the SureSelectXT Methyl-Seq Target Enrichment System (Agilent Technologies, Santa Clara, CA, USA) and EZ-DNA Methylation-Gold Kit (Zymo Research, CA, USA) were used to build whole-genome bisulfite DNA libraries. SureSelectXT Human Methyl-Seq is the first comprehensive target enrichment system to enable researchers to focus on the regions where methylation is known to impact gene regulation: CpG islands, CpG island shores, promoters, and differentially methylated regions (DMRs). After generating clusters, we sequenced the library preparations on the Illumina Hiseq 2000/2500 platform (San Diego, CA, USA) to obtain 125/150bp paired-end reads. The Illumina CASAVA pipeline was used to perform image analysis and base calling. Lastly, we collected the 125/150bp paired-end reads.

Differentially Methylated Sites (DMS) and Differentially Methylated Regions (DMR)

The BiSeq software (version 1.6.0; Hebestreit K, 2013) was used to identify differentially methylated sites (DMSs) and differentially methylated regions (DMRs). Only the target regions were analyzed. The genomic regions termed CpG clusters are covered by CpG sites with a high spatial density. The definition of a DMR is the adjacent rejected CpG region in a CpG cluster. The status of all CpG sites that are hyper- or hypo-methylated is determined by the DMRs. The CpG site is hypomethylated when the difference in methylation switches from positive to negative. Finally, we identified DMRs with more than two CpG sites where the difference of methylation level was greater than 0.1. For experiments without replicates, DMSs and DMRs were identified using the methylKit software (version 1.0.0; Akalin A, 2012). For DMR identification, methylKit tiles the genome with windows of 1000bp in length and a step size of 1000bp and summarizes the methylation information on the tiles. MethylKit uses Fisher’s exact test to calculate the P-values for the differential methylation which are then adjusted to Q-values using the SLIM method (Wang, 2011). The differentially methylated regions/bases are selected from the calculated Q-values and percent methylation difference cutoffs. The threshold value was a Q-value <0.05 and a percent methylation difference larger than 10%.

GO, KEGG Enrichment Analysis of DMR-Related Genes

GO enrichment analysis was performed with the GOseq function in R after correction for the gene length bias. DMR-related genes were considered significantly enriched when the corrected P-value was less than 0.05 in GO terms. In this analysis, KOBAS software (Mao X, 2005) was used to analyze the statistical enrichment of DMG with the assistance of the KEGG database, which helps to understand the high-level functions and utilities of biological systems from the molecular level, particularly from large-scale molecular datasets generated by genome sequencing and other high-throughput experimental techniques ( M, 2008).

Quantitative Real‐Time Polymerase Chain Reaction (qPCR)

Total RNA from the CR-HCC tissues and APTs was extracted using the TRIzol reagent according to the manual (Invitrogen, USA) and the cDNA was reverse‐transcribed using the HiScript IIQ RT SuperMix for qPCR kit (Nanjing Nuoweizan Biological Technology Co., Ltd., China) according to the manufacturer’s instructions. The real-time qPCR was performed using the AceQTM qPCR SYBR Green Master Mix Kit (Nanjing Nuoweizan Biological Technology Co., Ltd., China). GAPDH served as an internal control and the relative mRNA expressions were calculated using the equation 2−∆∆Ct. All reactions were performed at least three times. The PCR primer sequences are listed in Table 1.

Table 1

Nucleotide Sequences of Primers for Real-Time qPCR

Genes	Sequence(5′to 3′)	Genes	Sequence(5′to 3′)
GAPDH(F)	CTAAGGCTGTGGGGAAGGT	TNNI2(F)	AGATGTCAGAGTAAGTATTCAGATGTTTTG
GAPDH(R)	GTTGTTGACCTGACCTGCC	TNNI2(R)	CATCTTTTTTCTGTTCATGAACATGAGTTC
HOXB-AS3(F)	CCCTCCAAGTCCAGTAAGAAGT	TIRAP(F)	TGGTGCAAGTACCAGATGCT
HOXB-AS3(R)	AGATCCTAAGAGGTGCGAGTTTA	TIRAP(R)	GACTTGACGAAAGCCACCAT
HOXB6(F)	CACTCCGGTCTACCCGTGGATGCA	SHC3(F)	AAGATGCTGGAGGAACTGCAA
HOXB6(R)	CATATCTTGATCTGCCTCTCCGTCAG	SHC3(R)	CCAGGAAGTCTCCGTCTTTCTC
HOXB3(F)	TGGGAAAGGGCGAACGAGGATTAGGG	CTTN(F)	AAAGCTTCAGCAGGCCAC
HOXB3(R)	TGTATAGATATGAAATGTTCTAGA	CTTN(R)	TTTGGTCCTGTTTCAAGTTCC
USP18(F)	CCTGAGGCAAATCTGTCAGTC	TFAP2A(F)	TTTTCAGCCATGGACCGTCA
USP18(R)	CGAACACCTGAATCAAGGAGTTA	TFAP2A(R)	CGTTGACGTGGGAGTAAGGA
MAP3K6(F)	CCTGGAGGTGATTCTGCCCATT
MAP3K6(R)	GGTCTCGACTGCGATGGTCTTTAA

Abbreviations: F, forward primers; R, reverse primers.

Ethics Statement

All protocols in this study were approved by the Ethics Committee of the First Affiliated Hospital of Anhui Medical University. Written informed consent was gained from all patients in this study at the time of surgery.

Statistical Methods

Statistical analyses were performed with GraphPad Prism v. 8 (GraphPad Software 160 Inc., San Diego, CA, USA). The unpaired Student’s t-test was used to evaluate statistical significance. Results with P-values (P) < 0.05 were considered statistically significant.

Results

DNA Methylation Changes Between CR-HCC Tissues and APTs

To describe the methylome in CR-HCC and unravel the epigenetic basis of CR-HCC, we firstly examined if any individual sites exhibited differential DNA methylation in CR-HCC tissues compared with adjacent peritumor tissues. We showed substantial differences in DNA methylation between the CR-HCC tissues and the APTs and comprehensively identified a large subset of CpG sites/DMRs/genes that correlated with CR-HCC. A total of 38,554 DMSs were analyzed in CR-HCC tissues and APTs by using the human SureSelect XT arrays (Agilent Technologies, Santa Clara, CA, USA). Compared with the APTs, 21,853 DMSs in the CR-HCC tissues showed higher CpG methylation levels while 16,701 DMSs showed lower CpG methylation levels.

A DMR is a region in the genome that shows differences in methylation levels between samples. It is both more accurate and biologically relevant to use DMRs rather than CpG sites. The use of DMRs can also decrease potential artifacts due to random methylation in single CpG sites. The criterion used for DMR inclusion was a change in DNA methylation at a locus-specific level with an α≥10% methylation change with a p-value<0.05. Compared with the APTs, 1728 DMRs were identified in tumor tissues by using the Differential DNA Methylation Database (MethDiff) methods, of which 868 DMRs were hypermethylated and 860 DMRs were hypomethylated in the CR-HCC tissues. For characterization of the DMRs, we identified differences in the distribution of DMRs related to gene structures via the percentage of DMRs located within the promoter, exon, 5′ untranslated region (5′ UTR), and 3′untranslated region (3′ UTR). Similarly, we calculated the percentages of DMRs within a CpG island, CpG island shore, and CpG island shelf to determine differences in the distribution related to CpG islands. Our analyses showed prominent distribution of DMRs within promoter regions (1626), CpG islands (1821), CpG island shores (1792), CpG island shelves (982), and, more rarely, within 5′ UTRs (145) and 3′ UTRs (127)(Figure 1).

Figure 1

The distribution DMRs in the various genomic elements. The horizontal coordinate represents different functional areas, and the vertical coordinate represents the proportion of DMR in different functional areas.

DMR Analysis During CR-HCC Tissues and APTs

After identifying the DMRs between the CR-HCC and APT groups, we used the physical position in the genome and the annotation information of the species to annotate the genes associated with the DMRs. A total of 2091 genes were identified from the DMR datasets, calculated by the MethDiff method. These1129 genes were related to DMRs with high methylation levels in the CR-HCC group and included 242 genes with methylation located in the gene promoter regions. A total of 962 genes were related to DMRs with low methylation levels in the CR-HCC group, of which 202 genes were differentially methylated the gene promoter. The data show that it is significantly more common for DMRs to be hypermethylated rather than hypomethylated in the CR-HCC tissues. Tables 2 and 3 show the top 30 hyper- or hypomethylated DMRs, ranked by statistical significance. After reviewing the GenBank database and the relevant literature, we selected the hypermethylated genes HOXB-AS3, HOXB6, HOXB3, USP18, and MAP3K6 and the hypomethylated genes TIRAP, TNNI2, SHC3, CTTN, and TFAP2A for further analysis of their expression by real-time qPCR. The results showed that the expression of HOXB-AS3, HOXB3, and MAP3K6 was downregulated in the CR-HCC tissues compared to the APTs, while the mRNA expression levels of TIRAP, SHC3, and CTTN were upregulated. The mRNA expression levels of HOXB6, USP18, TNNI2, and TFAP2A were not significantly different (Figure 2).

Table 2

Top30 Significant HyperMethylated DMRs in CR-HCC Tumor Tissues Compared with Adjacent Non-Tumor Tissues

Chrom	No.CG	Meth.CR-HCC	Meth.APTs	Meth.Diff	Meth.p	Direction	Gene_Name
Chr13	6	0.63922	0.412935	0.226284	2.20E-113	Hyper
Chr13	6	0.63922	0.412935	0.226284	2.20E-113	Hyper	LINC00398
Chr17	9	0.341825	0.18653	0.155295	6.46E-90	Hyper	HOXB-AS3
Chr17	9	0.341825	0.18653	0.155295	6.46E-90	Hyper	HOXB6
Chr17	9	0.341825	0.18653	0.155295	6.46E-90	Hyper	HOXB3
Chr11	7	0.495095	0.351013	0.144788	1.13E-64	Hyper
Chr11	7	0.495095	0.351013	0.144788	1.13E-64	Hyper	ORAOV1
Chr11	7	0.495095	0.351013	0.144788	1.13E-64	Hyper
Chr19	6	0.210449	0.104174	0.106275	2.21E-47	Hyper
Chr19	6	0.210449	0.104174	0.106275	2.21E-47	Hyper
Chr19	6	0.210449	0.104174	0.106275	2.21E-47	Hyper
Chr19	6	0.210449	0.104174	0.106275	2.21E-47	Hyper	ICAM5
Chr19	6	0.210449	0.104174	0.106275	2.21E-47	Hyper
Chr1	8	0.711672	0.49796	0.213712	1.20E-41	Hyper
Chr1	8	0.711672	0.49796	0.213712	1.20E-41	Hyper
Chr1	8	0.711672	0.49796	0.213712	1.20E-41	Hyper
Chr1	8	0.711672	0.49796	0.213712	1.20E-41	Hyper	RP4-794H19.2
Chr1	8	0.711672	0.49796	0.213712	1.20E-41	Hyper	RP11-63G10.3
Chr6	8	0.93899	0.832476	0.104991	1.40E-41	Hyper
Chr6	8	0.93899	0.832476	0.104991	1.40E-41	Hyper	XXbac-BPG308K3.6
Chr6	8	0.93899	0.832476	0.104991	1.40E-41	Hyper	RPL13P
Chr22	6	0.837569	0.539437	0.296242	1.45E-41	Hyper
Chr22	6	0.837569	0.539437	0.296242	1.45E-41	Hyper
Chr22	6	0.837569	0.539437	0.296242	1.45E-41	Hyper	USP18
Chr22	6	0.837569	0.539437	0.296242	1.45E-41	Hyper
Chr1	2	0.572219	0.447216	0.125002	1.58E-40	Hyper
Chr1	2	0.572219	0.447216	0.125002	1.58E-40	Hyper
Chr1	2	0.572219	0.447216	0.125002	1.58E-40	Hyper	MAP3K6
Chr11	14	0.669233	0.534945	0.124809	2.47E-32	Hyper

Table 3

Top 30 Significant hypoMethylated DMRs in CR-HCC Tumor Tissues Compared with Adjacent Non-Tumor Tissues

Chrom	No.CG	Meth.CR-HCC	Meth.APTs	Meth.Diff	Meth.p	Direction	Gene_Name
Chr11	3	0.773586	0.895064	−0.12148	3.54E-31	Hypo
Chr11	3	0.773586	0.895064	−0.12148	3.54E-31	Hypo
Chr11	3	0.773586	0.895064	−0.12148	3.54E-31	Hypo	TIRAP
Chr8	26	0.486188	0.855095	−0.36891	4.99E-30	Hypo
Chr8	26	0.486188	0.855095	−0.36891	4.99E-30	Hypo
Chr11	15	0.488482	0.765907	−0.27743	2.59E-28	Hypo	TNNI2
Chr11	15	0.488482	0.765907	−0.27743	2.59E-28	Hypo	SYT8
Chr9	2	0.101594	0.213265	−0.11167	2.69E-26	Hypo	SHC3
Chr9	2	0.101594	0.213265	−0.11167	2.69E-26	Hypo
Chr6	5	0.253251	0.379558	−0.12631	5.70E-26	Hypo
Chr6	5	0.253251	0.379558	−0.12631	5.70E-26	Hypo
Chr6	5	0.253251	0.379558	−0.12631	5.70E-26	Hypo	TFAP2A
Chr6	5	0.253251	0.379558	−0.12631	5.70E-26	Hypo	TFAP2A-AS1
Chr11	10	0.022787	0.49504	−0.4729	2.34E-24	Hypo
Chr11	10	0.022787	0.49504	−0.4729	2.34E-24	Hypo	PHLDB1
Chr1	7	0.738378	0.866507	−0.12511	1.51E-22	Hypo
Chr1	7	0.738378	0.866507	−0.12511	1.51E-22	Hypo	BAI2
Chrx	16	0.309552	0.631625	−0.32207	2.37E-21	Hypo
Chrx	16	0.309552	0.631625	−0.32207	2.37E-21	Hypo	ERAS
Chrx	16	0.309552	0.631625	−0.32207	2.37E-21	Hypo
Chrx	16	0.309552	0.631625	−0.32207	2.37E-21	Hypo
Chrx	33	0.200289	0.435692	−0.25939	5.65E-21	Hypo	PCSK1N
Chrx	33	0.200289	0.435692	−0.25939	5.65E-21	Hypo
Chrx	33	0.200289	0.435692	−0.25939	5.65E-21	Hypo
Chrx	33	0.200289	0.435692	−0.25939	5.65E-21	Hypo
Chr17	15	0.33953	0.550398	−0.21731	7.80E-20	Hypo	ALOX12P2
Chr17	15	0.33953	0.550398	−0.21731	7.80E-20	Hypo	AC027763.2
Chr17	15	0.33953	0.550398	−0.21731	7.80E-20	Hypo
Chr11	11	0.721749	0.836882	−0.10649	1.44E-18	Hypo	CTTN
Chr1	4	0.808155	0.910256	−0.1021	2.90E-18	Hypo
Chr1	4	0.808155	0.910256	−0.1021	2.90E-18	Hypo
Chr1	4	0.808155	0.910256	−0.1021	2.90E-18	Hypo	SLC45A1

Abbreviations: Chrom, chromosome number; DMR start, the starting position of DMR; DMR end, end position of DMR; no.CG, number of CpG sites in this region; meth.CR-HCC, the average methylation level of CR-HCC in this region; meth.APTs, average methylation level of APTs in this region; meth-diff, difference of mean methylation level between CR-HCC and APTs in the DMR region; meth-p, p value of the significance difference between CR-HCC and APTs in the DMR region; Direction, the mean methylation levels of CR-HCC in the DMR region, hyper (high), hypo (low), relative to APTs.

Figure 2

The relative mRNA expressions of HOXB-AS3, HOXB6, HOXB3, USP18, MAP3K6, TNNI2, TIRAP, SHC3, CTTN, and TFAP2A were detected in CR-HCC tissues compared with APTs (control). The mRNA levels of these genes were analyzed using real-time qPCR and normalized to the GAPDH. (*P<0.05, **P <0 0.01, ***P<0.001).

Bioinformatic Analysis of Methylation Data

The hypergeometric distribution of these genes in relation to specific GO classifications was calculated for selected DMR genes. The number and functional classifications of the DMGs in relation to the GO categories, including the major enriched categories of molecular function, biological processes, and cellular components were displayed as a histogram. Further analysis of the enriched GO terms for both the hypermethylated and hypomethylated DMGs suggested that the molecular functions of “binding”, “protein binding”, and “cytoskeletal protein binding” were the main categories for the hypermethylated DMGs. The hypomethylated DMGs were mostly enriched in the molecular functions “binding”, “protein binding”, “calcium ion binding”, among others (Figure 3A and B). These molecular functions of the GO enrichment analysis in DMGs are shown in the Directed Acyclic Graph (DAG) graphical display (Figure 4A and B).

Figure 3

(A) The significant GO categories of hypermethylated genes in CR-HCC tissues compared with APTs; (P<0.05). (B) The significant GO categories of hypomethylated genes in CR-HCC tissues compared with APTs (P<0.05). “*” means the GO term is significantly enriched.

Figure 4

(A) Molecular function DAG of hypermethylated DMG. (B) Molecular function DAG of hypomethylated DMG. S each node represents A GO term, and the box represents the Top10 GO term with the degree of enrichment. The darker the color, the higher the degree of enrichment. The name of the TERM and the p-values of the enrichment analysis are displayed on each node.

KEGG pathway analysis suggested that the significantly differentially hypermethylated genes were enriched in several pathways including the “Estrogen signaling pathway” (P=5.91E-03)(including the genes ADCY4, ADCY2, ADCY8, ADCY7, EGFR, AK72, FOS, GPER1, KCNJ6, KCNJ9, NRAS, PIK3R2, PLCB2, and MAP2K2), the “Notch signaling pathway” (P=9.04E-03) (including the genes RBPJL, CTBP1, DTX1, HES5, JAG2, NOTCH4, NUMBL, NCOR2, and MAML1), and ‘Signaling pathways regulating pluripotency of stem cells’ (P=1.26E-02)(including the genes DLX5, APC2, AKT2, FGFR3, JAK2, JAK3, ZFHX3, NRAS, MAPK11, REST, MAP2K2, IK3R2, BMPR11A, WNT10A, FZD1, WNT5B, and HAND1), the “Rap1 signaling pathway” (P=1.41E-02) (including the genesADCY4, EPHA2, RASGRP2, EFNA5, RAPGEF3, ADCY8, RGS14, EGFR, ADCY2, ADCY7, AKT2, FGFR3, PFN3, NRNAS, RALGDS, FIK3R2, RALB, RALA, MAP2K2, MAPK11, PLCB2, and CDH1) while differentially hypomethylated genes were enriched in several pathways such as “Proteoglycans in cancer” (P=5.04E-03) (including the genes FLNC, AKT1, CTTN, ANK1, GPC1, GFR, IGF2, HOXD10, MAPK11, MAP2K2, PIK3CD, PLCE1, CCND1, TGFB1, FZD3, WNT5A, ACTG1, WNT7B, KDR, CAMK2B, WNT5B, WNT10A, and HPSE2) and “Calcium signaling pathway” (P=7.79E-03) (including the genes ADORA2A, ADORA2B, CRIN2A, CRIN2D, GNAS, GNAL, GRIN2C, GRIN1, GRM5, HTR5A, NTSR1, PPP3CC, PLCE1, CALML5, P2RX1, RYR1, CACNA1A, CAMK2B, EGFR, CACNA1H, and ORAI3). Table 4(A–B) shows all the significant KEGG pathways (p<0.05). Both hypermethylated and hypomethylated genes with hypermethylation and hypomethylation were mostly related to specific pathways including the “Rap1 signaling pathway” and ‘Signaling pathways regulating pluripotency of stem cells’.

Table 4

KEGG Pathway Analysis of Hypermethylated and Hypomethylated Genes

Pathyway Name	ID	Gene	EntrezGene	Statistic
A, KEGG pathway analysis of hypermethylated genes
Oxytocin signaling pathway	04921	22	ADCY4,ADCY2,CACNA2D4,MEF2C,PPP1CC,PIK3R2,TRPM2,PPP1R12C,MAP2K2,ADCY8,NRAS,MAPK7,PLCB2,EGFR,ADCY7,KCNJ9,KCNJ6,OXT,CACNG3,CACNG7,CACNG8,FOX	C=159O=22P=1.19E-03
Estrogen signaling pathway	04915	14	ADCY4,ADCY2,ADCY8,ADCY7,EGFR,AK72,FOS,GPER1,KCNJ6,KCNJ9,NRAS,PIK3R2,PLCB2,MAP2K2	C=100O=14P=5.91E-03
Notch signaling pathway	04330	9	RBPJL,CTBP1,DTX1,HES5,JAG2,NOTCH4,NUMBL,NCOR2,MAML1	C=48O=9P=9.04E-03
Signaling pathways regulating pluripotency of stem cells	04550	17	DLX5,APC2,AKT2,FGFR3,JAK2,JAK3,ZFHX3,NRAS,MAPK11,REST,MAP2K2,FIK3R2,BMPR11A,WNT10A,FZD1,WNT5B,HAND1	C=142O=17P=1.26E-02
Endometrial cancer	05213	8	CTNNA1,APC2,EGFR,AKT2,NRAS,FIK3R2,MAP2K2,CDH1	C=52O=8P=1.39E-02
Rap1 signaling pathway	04015	22	ADCY4,EPHA2,RASGRP2,EFNA5,RAPGEF3,ADCY8,RGS14,EGFR,ADCY2,ADCY7,AKT2,FGFR3,PFN3,NRNAS,RALGDS,FIK3R2,RALB,RALA,MAP2K2,MAPK11,PLCB2,CDH1	C=211O=22P=1.41E-02
Prolactin signaling pathway	04917	10	AKT2,FOS,JAK3,NRAS,FIK3R2,MAPK11,MAP2K2,TNFRSF11A,SOCS6,SOCS3	C=72O=10P=1.52E-02
Adrenergic signaling in cardiomyocytes	04261	17	ADCY7,ADCY4,CACNG3,RAPGEF3,ADCY2,ADCY8,ADRA1A,AKT2,ATP2B1,PPP1CC,MAPK11,CACNG8,PIK3R2,PLCB2,CACNG7,TPM1,CACNA2D4	C=149O=1 71.88E-02
Progesterone-mediated oocyte maturation	04914	11	ADCY4,ADCY2,ADCY8,ADCY7,AKT2,RPS6KA6,MAPK11,FIK3R2,RPS6KA2,MAPD1L1,CCNA1	C=86O=11P=2.06E-02
Inflammatory mediator regulation of TRP channels	04750	11	ADCY2,ADCY8,ADCY7,ALOX12,ADCY4,PRKCH,PPP1CC,MAPK11,FIK3R2,PLCB2,PRKCE	C=99O=11P=2.45E-02
Cholinergic synapse	04725	13	ADCY4,CHAT,ADCY2,ADCY8,ADCY7,AKT2,FOS,KCNJ6,JAK3,NRAS,FIK3R2,PLCB2,SLC18A3	C=113O=13P=2.98E-02
Adherens junction	04520	9	SORBS1,BAIAP2,EGFR,CTNNA1,LMO7,PTPRM,ACP1,ACTN4,CDH1	C=73O=9P=3.61E-02
SNARE interactions in vesicular transport	04130	4	GOSR2,STX2,VAMP2,STX5	C=34O=4P=3.83E-02
Transcriptional misregulation in cancer	05202	18	FLI1,HOXA11,HOXA9,HOXA10,MLLT1,H3F3C,MEF2C,PAX7,PAX5,SPINT1,RXRA,HIST2H3D,PAX8,ZBTB16,WT1,ZBTB17,HIST1H3D,CCNA1	C=179O=18P=4.59E-02
Long-term potentiation	04720	9	RAPGEF3,ADCY8,RPS6KA6,GRIN2D,NRAS,PPP1CC,PLCB2,MAP2K2,RPS6KA2	C=67O=9P=4.97E-02
C, the number of reference genes in the category; O is the number of genes in the gene set and also in the category, p-value from hypergeometric test
B, KEGG pathway analysis of hypomethylated genes
Endometrial cancer	05213	9	EGFR,CTNNA2,FOXO3,AKT1,CCND1,PIK3CD,MAP2K2,TCF7	C=52O=9P=8.51E-04
Nicotine addiction	05033	10	GRIN3B,GRIN3A,GRIN2A,GRIN2D,GRIN2C,GRIN1,GABRD,GABRA5,GABRB3,CACNA1A	C=40O=10P=1.35E-03
Glioma	05214	10	CDKN2A,EGFR,AKT1,CALML5,PDGFA,MAP2K2,SHC3,CCND1,PIK3CD,CAMK2B	C=65O=10P=4.33E-03
Acute myeloid leukemia	05221	8	CEBPA,AKT1,MAP2K2,CCND1,PIK3CD,TCF7,RUNX1,CCNA1	C=57O=8P=4.62E-03
Signaling pathways regulating pluripotency of stem cells	04550	17	DUSP9,AKT1,NEUROG1,ISL1,ID4,OTX1,MAPK11,MAP2K2,PIK3CD,LHX5,ZIC3,FZD3,WNT5A,WNT7B,AXIN2,WNT5B,WNT10A	C=142O=17P=4.96E-03
Proteoglycans in cancer	05205	23	FLNC,AKT1,CTTN,ANK1,GPC1,EGFR,IGF2,HOXD10,MAPK11,MAP2K2,PIK3CD,PLCE1,CCND1,TGFB1,FZD3,WNT5A,ACTG1,WNT7B,KDR,CAMK2B,WNT5B,WNT10A,HPSE2	C=204O=23P=5.04E-03
Amphetamine addiction	05031	12	GRIN3B,GRIN3A,GRIN2A,GRIN2D,ARC,GRIN2C,GNAS,GRIN1,CALML5,PPP3CC,CAMK2B,CREB5	C=68O=12P=5.93E-03
Calcium signaling pathway	04020	21	ADORA2A,ADORA2B,CRIN2A,CRIN2D,GNAS,GNAL,GRIN2C,GRIN1,GRM5,HTR5A,NTSR1,PPP3CC,PLCE1,CALML5,P2RX1,RYR1,CACNA1A,CAMK2B,EGFR,CACNA1H,ORAI3	C=180O=21P=7.79E-03
Melanoma	05218	7	CNKN2A,EGFR,AKT1,PDGFA,MAP2K2,CCND1,PIK3CD	C=71O=7P=7.96E-03
Chronic myeloid leukemia	05220	10	CTBP1,CNKN2A,CTBP2,SHC3,AKT1,MAP2K2,CCND1,PIK3CD,CAMK2B,TGFB1	C=73O=10P=9.59E-03
Rap1 signaling pathway	04015	20	ADORA2A,EGFR,EFNA2,ADORA2B,GRIN2A,GNAS,AKT1,GRIN1,KDR,INSR,NGF,MAPK11,MAP2K2,PIK3CD,PLCE1,CALML5,PRKCI,PDGFA,ACTG1,LPAR2	C=211O=20P=1.29E-02
Long-term potentiation	04720	10	GRIN2A,GRIN2D,GRIN2C,GRIN1,GRM5,CALML5,PPP3CC,MAP2K2,RPS6KA2,CAMK2B	C=67O=10P=1.29E-02
Neuroactive ligand-receptor interaction	04080	29	GRIN3B,DRD4,CHRNE,ADORA2A,HEH3,GPR156,ADRA2A,ADORA2B,GRIN3A,GRM5,GRIN1,GABRD,GABRB3,GRIN2C,P2RY8,GRM6,NPBWR1,GRIN2D,GRIN2A,GABRA5,GRIK4,HTR5A,OPRL1,NTSR1,P2RY11,PRSS3,P2RX1,S1PR4,LPAR2	C=275O=29P=1.84E-02
Non-small cell lung cancer	05223	7	EGFR,CNKN2A,FOXO3,AKT1,MAP2K2,CCND1,PIK3CD	C=56O=7P=2.59E-02
Wnt signaling pathway	04310	16	CSNK1A1,CTBP2,WIF1,CTBP1,NFATC1,LRP6,CCND1,PPP3CC,TCF7,FZD3,WNT5A,WNT7B,WNT5B,WNT10A,CAMK2B,AXIN2	C=140O=16P=2.97E-02
Pathways in cancer	05200	26	CEBPA,CNKN2A,CTNNA2,EGFR,CTBP2,DAPK3,CTBP1,EPAS1,AKT1,GSTP1,MAH2,PIK3CD,MAP2K2,CCND1,PDGFA,TGFB1,TCF7,WNT5A,TRAF3,FZD3,WNT7B,CAMK2B,WNT5B,AXIN2,CCNA1,WNT10A	C=327O=26P=3.44E-02
Hippo signaling pathway	04390	17	WTIP,GDF7,DLG1,CTNNA2,AMH,SCRIB,CCND1,TCF7,WNT5A,TGFB1,FZD3,WNT7B,ACTG1,AXIN2,WNT10A,WNT5B,PRKCI	C=154O=17P=3.47E-02
MAPK signaling pathway	04010	22	DUSP5,EGFR,DUSP9,DUSP8,AKT1,MAPK8IP3,FLNC,ARRB1,MAPT,NFATC3,MEF2C,NGF,MAPK11,PPP3CC,MAP2K2,ARRB2,PDGFA,RPS6KA2,TGFB1,CACNA1A,CACNB1,CACNA1H	C=257O=22P=3.53E-02
Cocaine addiction	05030	8	GRIN3B,GRIN3A,GRIN1,GNAS,GRIN2C,GRIN2D,GRIN2A,CREB5	C=50O=8P=3.54E-02
cAMP signaling pathway	04024	20	GRIN3B,HCN4,ADORA2A,GRIN3A,GRIN1,AKT1,GRIN2C,AMH,GNAS,GRIN2D,GRIN2A,NFATC3,MAP2K2,PIK3CD,PLCE1,CALML5,PDE4D,ACOX3,CAMK2B,CREB5	C=200O=20P=3.64E-02
Amyotrophic lateral sclerosis (ALS)	05014	7	GRIN2D,GRIN2C,GRIN1,GRIN2A,MAPK11,PRPH,PPP3CC	C=51O=7P=3.86E-02
Prolactin signaling pathway	04917	8	CISH,FOXO3,AKT1,MAPK11,MAP2K2,SHC3,CCND1,PIK3CD	C=72O=8P=4.25E-02
Colorectal cancer	05210	7	AKT1,MSH2,CCND1,PIK3CD,TGFB1,AXIN2,TCF7	C=62O=7P=4.27E-02
Oxytocin signaling pathway	04921	15	EEF2,EGFR,GNAS,MEF2C,NFATC3,PIK3CD,CCND1,MAP2K2,PPP3CC,CALML5,OXT,RYR1,ACTG1,CACNB1,CAMK2B	C=159O=15P=4.41E-02
Glutamatergic synapse	04724	14	GRIN3B,GRIN3A,GRM6,GRM5,GRIN1,SHANK2,GNAS,GRIN2C,GRIN2D,GRIN2A,GRIK4,PPP3CC,CACNA1A,SHANK3	C=116O=14P=4.72E-02
Axon guidance	04360	14	SEMA4D,DPYSL2,SEMA6B,EFNA1,EPHA7,LRRC4C,PPP3CC,PLXNA3,SEMA5B,PLXNB3,SLIT1,CXCL12,SEMA3B,NTNG2	C=127O=14P=4.73E-02

Abbreviations: C, the number of reference genes in the category; O is the number of genes in the gene set and also in the category, p-value from hypergeometric test.

Discussion

DNA methylation that can regulate gene expression is one of the most important epigenetic modifications in humans. Previous studies have discovered aberrant DNA methylation in HCC tissues rather than APTs.19 The main aims of this research were to study the associations between methylation and CR-HCC and to find new potential biomarkers to diagnose CR-HCC and differs from previous studies in the design and enrolment criteria in the selection of patients for the study.

The two most common forms of aberrant CpG methylation in cancer, namely promoter hypo- or hyper-methylation, may lead to the inappropriate activation of oncogenes or the silencing of tumor suppressor genes (TSGs) while global hypo-methylation causing chromosomal instability, have been widely studied in HCC.20,21 Several recent studies have suggested that aberrant CpG methylation of several tumor suppressor genes such as Ras-association domain family 1 isoform A (RASSF1A), Runt-related transcription factor 3 (RUNX3), E-cadherin, and P16, among others, and cause tumor suppressor gene inactivation in HBV- and HCV-related HCC, implying that epigenetic modification plays an important role in hepatocarcinogenesis.22–25 In the present study, we detected 38,554 differentially methylated CpG sites across the whole genome of CR-HCC and APTs, including 16,701 hypermethylated CpG sites (42.9%) and 21,853 hypomethylated CpG sites (56%). In addition, we selected significant DMRs using stringent criteria. We identified 1728 DMRs, revealing 868 significantly hypermethylated DMRs and 860 significantly hypomethylated DMRs in the CR-HCC group compared with the APT group using MethDiff methods. Our analyses also identified prominent DMR distributions within CpG islands (958), CpG island shores (1051), CpG island shelves (509), and, more rarely, within the 5′ UTR (127) and 3′ UTR (145). Our data also showed that the DMRs were predominantly located within promoter regions (871). Since the promoters in CpG islands generally regulate gene transcription through DNA methylation, DMRs may regulate CR-HCC pathogenesis mainly via the transcriptional inactivation of tumor suppressor genes. It is known that aberrant DNA methylation occurs in CR-HCC tissues, and our results showed that hypermethylation of CpG sites occurs more frequently than hypomethylation. Our results are consistent with those of previous HCC genome-wide methylation studies.26–29

A total of 1129 genes were found to be related to DMRs with high methylation levels in the CR-HCC group, among which methylation occurred in 242 genes in the gene promoter region. A further 962 genes were found to be related to DMRs with low methylation levels in the CR-HCC group, among which methylation occurred in 202 genes in the gene promoter region. Moreover, this study identified the top 30 significantly hyper-and hypo-methylated DMRs and DMR-related genes in CR-HCC tissues compared with APTs, including TIRAP, TFAP2A, CTTN, USP18, ICAM5, and MAP3K6. The real-time qPCR results showed that the expression of HOXB-AS3, HOXB3, and MAP3K6 was downregulated in the CR-HCC tissues compared to the APTs. However, the mRNA expression of TIRAP, SHC3, and CTTN was upregulated in the CR-HCC tissues. The results confirm that the methylation status within the gene structure affects the gene expression in CR-HCC tissues. An increasing number of studies have demonstrated that TLR signaling pathways are modulated by TIR domain-containing adaptors, such as MyD88, TIRAP, and TRIF. TIRAP is specifically involved in the MyD88-dependent pathway via TLR2 and TLR4.30,31 In NSCLC cells, knockdown of TIRAP expression by siRNA suppressed proliferation.32 Cortactin (CTTN) has been identified as a biomarker of cancer metastasis due to its overexpression in various cancers, such as colorectal cancer, head and neck squamous cell carcinoma, among others. Results have suggested that CTTN inhibits the degradation of EGF and promotes the proliferation of CRC cells.33 A previous study has demonstrated that CTTN is a novel HBx-interacting protein and that the HBx-CTTN interaction inhibits cell cycle arrest to enhance cellular proliferation and maintain cellular migration.34 Ubiquitin-specific peptidase 18, also known as UBP43, is an enzyme belonging to the specific ubiquitin protease family; altering USP18 levels affects the stability of the cyclin D1 protein in lung cancer35 with higher levels of USP18 expression influencing KRAS expression in K-Ras-driven murine lung cancers in vivo.36 Of all the identified genes, homeobox A9 (HOXA9) was the most frequently methylated in plasma samples from HCC patients and in HCC tissues, suggesting its potential as a useful HCC biomarker.37 The discovery of specifically hypermethylated genes leading to gene silencing or hypomethylated genes leading to increased transcription is helpful in identifying new factors that are vital for CR-HCC carcinogenesis and progression. Our study of DMRs provides an important association between the aberrant DNA methylation within these sites, which is a vital epigenetic mechanism, and CR-HCC carcinogenesis. These results suggest that the development of CR-HCC is regulated by different genes together with different signaling pathways and their crosstalk and that aberrant DNA methylation is an important event in the pathogenesis of CR-HCC.

The pathways associated with DNA methylation and CR-HCC carcinogenesis were screened by bioinformatic analysis. GO analysis showed that the hypermethylated DMGs were mostly associated with “binding”- and “adhesion” -related molecular functions, including cell adhesion, binding, cell-cell adhesion, plasma-membrane adhesion molecules. “Binding” is another enrichment category of the hypermethylated genes, and it contains “Ion binding” (GO:0043167 P=1.27E-09), “cell adhesion” (GO:0007155, P=1.75E-09), and “DNA binding” (GO:0003677, P= 1.43E-09). Other research has indicated that a suppressive oncogene, namely CDH1, encodes the epithelial cell adhesion molecule E-cadherin which is involved in cell-cell adhesion and cell polarity. The CHD1 gene is known to be an invasion and tumor suppressor gene and its decreased expression due to hypermethylation could promote tumor cell invasion and metastasis.38–40 Jiang L demonstrated that calcium-binding protein 39 (CAB39) promoted the carcinogenesis and development of HCC when CAB39 activated the ERK signaling pathway, suggesting that CAB39 may be useful as a biomarker for the prediction of prognosis as well as the potential of being a novel therapeutic target in clinical practice.41 The hypomethylated DMGs were enriched in the GO categories including “protein binding”, “extracellular ligand-gated ion channel activity”, “regulation of signal transduction”, “cytoskeletal protein binding”. These data show that genes enriched in a number of GO categories are related to the carcinogenesis of CR-HCC and may be useful as specific biomarkers for the diagnosis of CR-HCC.

Furthermore, the results of the KEGG analysis suggested that the hypermethylation-related genes were mainly enriched in the “Rap1 signaling pathway”, “signaling pathways regulating pluripotency of stem cells”, and “Estrogen signaling pathway” in CR-HCC tissues. The hypomethylation-related genes were mainly involved in “colorectal cancer”, “calcium signaling pathway”, “nicotine addiction”, among others. Rap1, a small GTPase, controls various processes, including cell polarity, cell-cell junction formation, and cell adhesion.42 As undifferentiated cells, pluripotent stem cells (PSCs), including embryonic stem cells (ES) and induced pluripotent stem cells (iPS), have the capacity for indefinite self-renewal and are able to generate all cell types of the three germinal layers.43 Human hepatic stem/progenitor cells (hHPCs) from adult liver that resides in the canals of Hering may have a great capacity to proliferate and the ability to differentiate into mature hepatocytes.44 Activation and dysregulated proliferation of hHPCs commonly contribute to the carcinogenesis of HCC in adverse environments. Furthermore, an increasing number of studies have demonstrated that HPCs show strong associations with HCC,45,46 which may provide new directions to study the relationship between methylation and CR-HCC. Other identified pathways such as “Cell cycle” and “Chemokine signaling pathway” are also considered to be important for the development and progression of malignant carcinoma, and the genes involved, together with their methylation patterns, may provide potential novel biomarkers for CR-HCC.

Conclusion

This preliminary work may pave the way for future studies to identify the vital functions of epigenetic mechanisms on CR-HCC carcinogenesis. The present study has several potential limitations, including the relatively small sample size. The sample size was restricted mainly due to the high cost of Methyl-seq, which precludes its application to a large-scale study. Nevertheless, our findings indicate the important role of epigenetic modification on the carcinogenesis of CR-HCC and pave the way for future analyses of a larger sample using a series of biochemical approaches.