Epigenome-wide DNA methylation profiling of portal vein tumor thrombosis (PVTT) tissues in hepatocellular carcinoma patients.

Aberrant methylation is a hallmark of hepatocellular carcinoma and plays an important role in tumor initiation and progression. However, the epigenome-wide methylation patterns of portal vein tumor thrombosis (PVTTs) have not been fully explored. Here, we performed epigenome-wide DNA methylation of adjacent normal tissues (ANTs), paired tumor tissues and paired PVTTs using an Infinium HumanMethylation450 array (n = 11) and conducted the Sequenom EpiTYPER assays to confirm the aberrantly methylated genes. MTS and apoptosis assay were used to assess the synergistic effect of two drugs on the HCC cell lines. We found the mean global methylation levels of HCC tissues and PVTTs were significantly lower than ANTs (P < 0.01). A total of 864 differentially methylated CpG sites annotated in 532 genes were identified between HCC tissues and paired PVTTs (|mean methylation difference|>10%, P < 0.005). The pathway analysis based on hypermethylated genes in PVTT tissues was interestingly enriched in regulation of actin cytoskeleton pathway (P = 4.48E-5). We found 23 genes whose methylation levels were gradually alternated in HCC tissues and PVTTs. Aberrant methylation status of TNFRSF10A, ZC3H3 and SLC9A3R2 were confirmed in a validation cohort (n = 48). The functional experiments demonstrated the combination of decitabine (DAC) and tumor necrosis factor-related apoptosis-inducing ligand (rh-TRAIL) could synergistically suppress the proliferation and induce apoptosis in SK-Hep-1 and Huh7 cell lines. Together, our findings indicated that DNA methylation plays an important role in the PVTT formation through regulating the metastasis-related pathways. The combination of DAC and rh-TRAIL might be a promising treatment strategy for HCC.

Introduction

Hepatocellular carcinoma (HCC) is one of the most aggressive human malignancies and one of the leading causes of cancer death [1]. Portal vein tumor thrombosis (PVTT) occurs in about 35–50% of HCC patients and represents an extremely poor prognosis, due to the high frequency of intrahepatic metastasis and the lethal major bleeding caused by the increased portal pressure [2]. So far, the mechanism of PVTT formation remains obscure. To fully understand the underlying mechanism of PVTT formation at cellular and molecular levels has become one of the major goals in HCC research.

DNA methylation is one of most studied epigenetic mechanisms and mostly occurs in the context of CG dinucleotides [3]. Aberrant DNA methylation status has emerged as an important hallmark in various cancers [4]. Abnormal methylation of tumor suppressor genes and oncogenes has been proven to play an important role in many malignancies, such as liver cancer [5], [6], [7], [8], [9], gastric cancer [10] and colorectal cancer [11], [12]. DNA methylation is a dynamic epigenetic marker accompanied with occurrence and development of cancer, which has widely effects on the prevention, diagnosis, prediction and therapy of the disease [13]. Although many previous researches have explored the methylation patterns of HCC [5], [6], [7], [8], [9], [14], the study exploring the DNA methylomes of PVTT is limited. Fully exploring the DNA methylomes of PVTT allowed us to develop a better understanding of special biological process involved in the PVTT formation.

In the current study, we aimed to identify the genes with alternated methylation status from normal tissue to tumor to PVTT, so we performed an epigenetic genome-wide analysis to evaluate the DNA methylation profile of PVTT in 11 HCC patients using the Infinium Human Methylation450 (HM450) BeadChip. Then, we replicated the candidates methylated genes in an independent validation cohort. Because the methylation of TNFRSF10A was identified to be gradually decreased accompanied with PVTT formation, we used MTS and apoptosis assay to investigate whether the rh-TRAIL displayed a combination effect with the DNA methyltransferase inhibitor (DAC) on the suppressing the proliferation of HCC cell lines.

Results

Clinical characteristics

Consistent with the previous RNA-sequencing study [15], 11 HCC patients with PVTT were included in the first stage of genome-wide epigenetic study and the clinical characteristics were presented in Table 1. All the included patients were male, and two patients drank over 50 g per day of alcohol. The median age was 54 years old. Ten patients were HBV-infected. Two patients were found with more than one tumor. All the PVTTs were in right portal veins.

Table 1

Clinicopathological characteristics of HCC patients with PVTT.

Character	Patients (n = 11)	Validation cohort(n = 48)	P value (χ2)
Gender (male/female)	11/0	45/3	1
Age (median(range), years)	54 (39–71)	46(26–77)	0.1190
Heavy alcohol consumption (Yes/No)	2/9	13/35	0.7119
Diabetes (Yes/No)	1/10	1/47	0.3407
Hypertension (Yes/No)	1/10	2/46	0.4680
AFP (Positive/Negative)	10/1	42/6	1
HBV infected (Yes/No)	10/1	45/3	0.5725
Tumor number (1/2)	9/2	32/16	0.4756
Max tumor diameters (median(range), cm)	8.4 (3–17.2)	8.7(2–29.4)	0.9666
Tumor necrosis (Yes/No)	9/2	33/15	0.4938
Liver cirrhosis (Yes/No)	11/0	42/6	0.5820

Heavy alcohol consumption: patients who drinks over 50 g per day of alcohol. χ2.

To validate the finding in epigenome-wide DNA methylation study, we collected a validation cohort. The validation cohort consisted of 48 HCC patients. All the HCC patients companied with PVTT. The clinical characteristics were also listed in the Table 1.

Epigenome-wide DNA methylation profiling of HCC tissues and PVTT tissues

We analyzed the epigenome-wide DNA methylation status in adjacent non-cancerous tissues (ANT), HCC tissues and PVTT tissues using the Infinium Human Methylation450 (HM450) BeadChip. Both in HCC tissues and PVTT tissues, the mean global methylation level was significantly lower than that in the ANTs (HCC vs. ANT P = 0.0039, PVTT vs. ANT P = 0.0062; Fig. 1A). We found the mean methylation levels of the body regions, 3′UTR regions and IGRs (intergenic region) in HCC tissues and PVTT tissues were significantly lower than that in ANTs (Fig. 1B). However, the mean methylation levels of TSS200 regions, 5′UTR regions and 1st Exon regions were similar among the three different tissues. Regarding the broader CpG context, the mean methylation levels of the HCC tissues and PVTT tissues in N-shelf region, S-shelf region and open sea region were significantly lower than that in the ANTs (Fig. 1C).

Fig. 1

The methylation status of ANT, HCC tissues and PVTT tissues (n = 11) (A) The mean global methylation levels of ANT, HCC tissues and PVTT tissues. The mean methylation levels of ANT, HCC tissues and PVTT tissues regarding (B) gene regions and (C) the broader CpG context. Manhattan plots show epigenetic profiles of all DMPs (D, ANT vs. HCC; E, ANT vs. PVTT; F, HCC vs. PVTT) (G) Heatmap of an unsupervised cluster analysis based on the differentially methylated CpG sites among those three tissues (PVTT vs. ANT |Δ methy diff. | > 10%, P < 1E−6; HCC vs. ANT |Δ methy diff. | > 10%, P < 5E−4; PVTT vs. HCC |Δ methy diff. | > 10%, P < 0.005).

Compared to ANTs, 1771 differentially methylated positions (DMPs) were detected in HCC tissues (|delta methylation difference (Δ methy diff.)| > 10%, P < 5E−4, Fig. 1D), and 41,614 DMPs were found in PVTT tissues (|Δ methy diff.| > 10%, P < 5E−4, Fig. 1E). The rate of hypomethylated CpG sites (62.9%) was obviously higher than the rate of hypermethylated CpG sites (37.1%) in HCC. The hypomethylated CpG sites (89.2%) was also more than the hypermethylated CpG sites (10.2%) in PVTT tissues. In addition, we found 864 DMPs between PVTT tissues and HCC tissues (|Δ methy diff.| > 10%, P < 0.005, Fig. 1F). However, the rate of hypomethylated sites was lower than that of hypermethylated sites (hypomethylated vs. hypermethylated, 14.9% vs. 85.1%). The heatmap, using all the DMPs among the different groups, showed clustering of the ANTs and a separate cluster of HCC tissues and PVTT tissues (Fig. 1G).

Besides, the hypomethylated CpG sites in HCC tissues compared to ANTs were mainly located at chromosomes 1, 2, 6, 7, 8 and 19 and the hypermethylated CpG sites were mainly located at chromosomes 1, 2, 6, 7, 11 and 19 (proportion >6%, Supplementary Fig. 1A and Supplementary Fig. 1B). When we compared PVTT tissues to ANTs, we found that the hypomethylated CpG sites mainly located at chromosomes 1, 2, 5, 7, 8 and 11 and the hypermethylated CpG sites mainly located at 1, 2, 6, 7, 11 and 19 (proportion >6%, Supplementary Fig. 1A and Supplementary Fig. 1B). In addition, the hypomethylated CpG sites in PVTT tissues were mainly located at chromosomes 2, 4, 5, 6, 11, 14, 16 and 22 and the hypermethylated CpG sites were mainly located at chromosomes 1, 2, 6, 11 and 17 when compared to the HCC (proportion >6%, Supplementary Fig. 1A and Supplementary Fig. 1B).

As for the distribution of the DMPs based on the relation to CpG island, we found, both in the HCC tissues and PVTT tissues, the hypomethylated CpG sites were mainly located at open sea regions whereas the hypermethylated CpG sites were mainly located at the CpG island (Supplementary Fig. 1C). We further analyzed the distribution of the DMPs based on the relation to the nearest gene regions. The aberrant methylated CpG sites which located at IGRs and body regions were the major proportion among all the aberrant CpG sites (Supplementary Fig. 1D).

Properties of the differentially methylated genes

To identify the potential functions, biological processes and networks, we performed the GO analysis, KEGG pathway and protein–protein interaction (PPI) analysis of the aberrant methylated CpG sites in PVTT tissues when compared to HCC tissues. For hypermethylated CpG sites, a total of 457 genes were finally identified. We found these genes mainly involved in embryonic morphogenesis (P = 1.22E−6), intracellular signaling cascade (P = 2.63E−5) and cell projection organization (P = 6.71E−5) (Fig. 2A). For the hypomethylated CpG sites, 75 genes were identified and used for further analysis. The GO analysis revealed that these genes were involved in embryonic development ending in birth or egg hatching (P = 0.0047), chordate embryonic development (P = 0.0048) and in utero embryonic development (P = 0.0206) (Fig. 2B). It is interesting to find that the formation of PVTT might be associated with the genes which were involved in embryonic developmental process.

Fig. 2

GO analysis and KEGG pathway analysis based on the differentially methylated CpG sites between PVTTs and HCC tissues. GO analysis based on (A) hypermethylated CpG sites (the 9 most enriched are shown, P < 0.05) and (B) hypomethylated CpG sites (P < 0.05). (C) The KEGG pathway analysis based on the hypomethylated CpG sites (the 9 most enriched pathways are shown, P < 0.05). (D) A protein–protein interaction (PPI) network of genes that are differentially methylated between PVTTs and HCC tissues.

The KEGG pathway analysis based on the genes annotated in the differentially methylated CpG sites were conducted using DAVID (https://david.ncifcrf.gov/summary.jsp). For the genes annotated in the hypomethylated CpG sites in PVTT tissues, the pathway was marginally significantly enriched in cell adhesion molecules (P = 0.0705). The genes annotated in the hypermethylated CpG sites were mostly involved in the regulation of actin cytoskeleton (FGFR2, ARHGEF4, PDGFA, GNA12, PIK3CD, ITGB2, GNG12, VAV2, CDC42, ITGA9, TIAM2, PDGFRB, PAK1, PIK3R3, MYLK, F2R) (P = 4.48E−5), T cell receptor signaling pathway (P = 0.0021), focal adhesion (P = 0.0036) and pathways in cancer (P = 0.0040) (Fig. 2C).

PPI analysis indicated the differentially methylated genes showed significantly direct connectivity (P = 1.52E−11, Fig. 2D). We observed a total of 374 direct connections with high confidence (combined score >0.7) among 221 genes by using STRING online tool (https://string-db.org). An extensive network included 168 genes and 344 edges was identified, in which CDC42 and FYN were the center. CDC42 regulates cell cycle signaling pathways that control diverse cellular functions including cell morphology, migration, endocytosis and cell cycle progression [16]. FYN plays important role in the regulation of cell growth and survival, cell adhesion, integrin-mediated signaling, cytoskeletal remodeling, cell motility, immune response and axon guidance [17].

Genes with gradually altered methylation from ANT to PVTT tissues

Firstly, we compared their differentially methylated sites (DMSs) log2 fold change (log2FC) signature and differential gene expression (DGE) log2 fold change (log2FC) signature to analyze the consistence between tumor and PVTT. We found the regression slopes between tumor log2-FC effect sizes compared with PVTT of both methylation (0.65, Fig. 3A) and mRNA expression (0.77, Fig. 3B), which indicated gradient of both transcriptomic severity and methylome severity with PVTT > tumor (slope < 1, P < 0.0001).

Fig. 3

The consistency between tumor tissues and PVTTs (A) The methylation consistence and (B) the expression consistence. (C) The flow diagram of our screening process. (D) The KEGG pathway analysis based on the genes whose methylation levels gradually increased or decreased.

We then proceeded to identify the DMPs which were altered in HCC tissues when compared to ANTs and the methylation levels were further changed in PVTT tissues. In order to obtain more potential gradient altered DMPs, we used |Δ methy diff.| > 5% and P < 0.01 as our screening standards. Finally, we identified a total of 340 gradually increased DMPs and 160 gradually decreased DMPs (Fig. 3C).

The 340 gradually increased DMPs were annotated in 194 genes. The KEGG pathway analysis revealed that these genes were enriched in the valine, leucine and isoleucine biosynthesis pathway and renal cell carcinoma pathway (P < 0.05, Fig. 3C and D). In addition, the 160 gradually decreased DMPs were detected and annotated in 90 genes (Fig. 3C). The KEGG pathway analysis found that the genes were involved in adherens junction, focal adhesion and extracellular matrix receptor interaction (ECM-receptor interaction) (P < 0.05, Fig. 3D and Supplementary Table 1). Consisted with the methylation analysis, the transcriptome analysis based on the dysregulated genes between HCC and PVTT also suggested that the ECM-receptor interaction and focal adhesion were significantly associated with the venous invasion [15].

To identify potentially novel signatures of DNA methylation in the process of the PVTT formation, we further focused on those genes with at least two CpG sites located at promoter regions (TSS1500, TSS200, 5′UTR and 1st Exon) or at least three CpG sites located at body regions. According to these criteria, 23 genes including 64 CpG sites were finally identified (Table 2). 17 genes (MICALL2, ZC3H3, NDRG2, TNFAIP2, SLC9A3R2, GNA11, COL16A1, TBX15, SLC25A2, HISTIH2AL, EVX1, DLX5, PAK1, NRN1L, TBKBP1, EGLN2, SHROOM4) showed gradually hypermethylated and 6 genes (ADORA2A, KIAA1143, NAPEPLD, PARVG, TNFRSF10A, TSTD1) were gradually hypomethylated (Fig. 3C and Table 3).

Table 2

The genes whose methylation levels gradually increased from ANT to PVTT tissues.

Probe ID	CHR	Position	Gene name	Gene region	CpG island	Methylation levels			Mean difference		P value
						PVTT	HCC	ANT	HCC-ANT	PVTT-HCC	HCC vs. ANT	PVTT vs. HCC
cg23305408	1	32,169,701	COL16A1	1stExon	Island	0.635	0.453	0.271	0.182	0.182	1.96E−03	2.51E−03
cg05003322	1	32,169,706	COL16A1	1stExon	Island	0.586	0.436	0.256	0.180	0.151	1.35E−03	6.16E−03
cg24842086	1	119,531,857	TBX15	5′UTR	N_Shore	0.544	0.395	0.167	0.227	0.149	1.50E−03	1.34E−03
cg05940231	1	119,532,189	TBX15	TSS200	Island	0.596	0.428	0.228	0.200	0.168	5.68E−04	2.14E−03
cg25340966	1	119,532,195	TBX15	TSS200	Island	0.594	0.420	0.134	0.287	0.174	6.41E−04	3.32E−03
cg25212131	5	140,683,670	SLC25A2	TSS200	Island	0.867	0.789	0.691	0.098	0.077	8.17E−03	2.78E−03
cg07039560	5	140,683,681	SLC25A2	TSS200	Island	0.796	0.705	0.548	0.157	0.091	1.65E−03	2.08E−03
cg01171339	6	27,831,894	HIST1H2AL	TSS1500	N_Shore	0.830	0.739	0.659	0.080	0.091	8.69E−03	6.53E−03
cg16051742	6	27,832,098	HIST1H2AL	TSS1500	N_Shore	0.702	0.590	0.473	0.118	0.111	2.39E−03	7.82E−03
cg17525357	7	1,491,380	MICALL2	Body	N_Shore	0.868	0.779	0.679	0.100	0.089	8.14E−04	2.69E−03
cg11787789	7	1,491,464	MICALL2	Body	N_Shore	0.814	0.760	0.694	0.066	0.054	1.19E−03	7.72E−03
cg14520214	7	1,491,831	MICALL2	Body	Island	0.843	0.775	0.671	0.104	0.069	2.23E−05	1.33E−03
cg16857641	7	1,492,452	MICALL2	Body	S_Shore	0.951	0.872	0.793	0.079	0.079	1.14E−03	1.20E−03
cg11251690	7	1,493,126	MICALL2	Body	S_Shore	0.943	0.881	0.767	0.115	0.062	4.20E−05	1.54E−03
cg02492818	7	1,493,153	MICALL2	Body	S_Shore	0.866	0.803	0.739	0.065	0.062	6.15E−04	1.17E−03
cg09071155	7	27,280,914	EVX1	TSS1500	N_Shore	0.434	0.287	0.160	0.127	0.147	1.59E−03	6.51E−03
cg17610800	7	27,281,561	EVX1	TSS1500	N_Shore	0.654	0.503	0.283	0.220	0.151	2.98E−05	7.47E−03
cg15798385	7	27,281,569	EVX1	TSS1500	N_Shore	0.763	0.601	0.342	0.259	0.162	8.12E−05	8.77E−03
cg09426103	7	96,654,554	DLX5	TSS1500	S_Shore	0.388	0.270	0.122	0.148	0.117	7.26E−03	6.38E−03
cg10243939	7	96,654,788	DLX5	TSS1500	S_Shore	0.562	0.400	0.188	0.212	0.162	7.85E−03	9.92E−03
cg17083494	7	96,654,832	DLX5	TSS1500	S_Shore	0.630	0.470	0.284	0.186	0.160	4.49E−03	2.38E−03
cg13191508	8	144,600,686	ZC3H3	Body	N_Shore	0.942	0.884	0.798	0.087	0.058	1.76E−03	2.27E−03
cg19598713	8	144,601,734	ZC3H3	Body	Island	0.850	0.798	0.687	0.111	0.052	7.97E−03	2.93E−03
cg25900150	8	144,601,851	ZC3H3	Body	Island	0.962	0.903	0.760	0.143	0.059	3.59E−03	2.93E−03
cg17202086	11	77,122,839	PAK1	5′UTR;5′UTR	Island	0.754	0.595	0.400	0.194	0.160	5.25E−04	5.19E−03
cg26996201	11	77,122,864	PAK1	5′UTR;5′UTR	Island	0.695	0.581	0.489	0.092	0.114	9.45E−03	6.36E−03
cg12269002	11	77,122,929	PAK1	5′UTR;5′UTR	Island	0.646	0.451	0.240	0.211	0.195	2.54E−03	9.91E−03
cg18309286	11	77,123,047	PAK1	5′UTR;5′UTR	Island	0.762	0.546	0.352	0.194	0.215	7.22E−03	5.91E−03
cg10078898	14	21,491,045	NDRG2	Body	N_Shore	0.753	0.695	0.525	0.170	0.057	4.55E−03	1.95E−03
cg23288563	14	21,491,469	NDRG2	Body	N_Shore	0.650	0.528	0.379	0.148	0.122	3.79E−03	9.97E−03
cg18081258	14	21,494,161	NDRG2	TSS1500	Island	0.708	0.557	0.365	0.191	0.152	2.13E−03	4.58E−03
cg09814127	14	103,593,235	TNFAIP2	Body	Island	0.789	0.707	0.561	0.147	0.082	3.81E−03	9.99E−03
cg18587137	14	103,593,503	TNFAIP2	Body	Island	0.857	0.734	0.603	0.131	0.123	9.20E−03	2.18E−03
cg18620571	14	103,593,505	TNFAIP2	Body	Island	0.878	0.772	0.639	0.134	0.105	9.46E−03	7.81E−03
cg02668822	16	2,086,611	SLC9A3R2	Body	Island	0.791	0.704	0.557	0.147	0.087	4.65E−03	7.99E−03
cg01142668	16	2,086,846	SLC9A3R2	Body	Island	0.969	0.901	0.675	0.226	0.068	2.31E−04	9.77E−04
cg09321109	16	2,087,555	SLC9A3R2	Body	N_Shore	0.832	0.760	0.566	0.195	0.072	8.65E−04	4.88E−03
cg09756419	16	67,918,680	NRN1L	TSS200	Island	0.787	0.685	0.517	0.168	0.102	4.62E−03	6.29E−03
cg02308079	16	67,918,686	NRN1L	TSS200	Island	0.757	0.656	0.497	0.159	0.101	3.53E−03	2.61E−03
cg01068601	16	67,918,699	NRN1L	TSS200	Island	0.894	0.770	0.592	0.179	0.123	6.47E−03	1.92E−03
cg03294833	16	67,918,813	NRN1L	1stExon	Island	0.813	0.699	0.508	0.191	0.114	2.30E−03	2.96E−03
cg04896381	17	45,772,467	TBKBP1	TSS200	Island	0.247	0.113	0.006	0.107	0.134	9.77E−03	3.91E−03
cg01822573	17	45,772,473	TBKBP1	TSS200	Island	0.256	0.106	0.005	0.101	0.149	5.92E−03	4.88E−03
cg25613170	17	45,772,523	TBKBP1	TSS200	Island	0.255	0.107	0.039	0.068	0.148	4.88E−03	9.77E−04
cg00168191	19	3,097,565	GNA11	Body	Island	0.853	0.796	0.736	0.060	0.058	4.76E−03	6.91E−03
cg24366168	19	3,097,712	GNA11	Body	Island	0.910	0.832	0.746	0.086	0.078	6.24E−03	3.11E−03
cg02261780	19	3,097,728	GNA11	Body	Island	0.863	0.754	0.660	0.094	0.110	6.83E−03	2.42E−03
cg08080060	19	41,307,090	EGLN2; EGLN2	1stExon; Body	Island	0.873	0.808	0.740	0.068	0.065	7.00E−03	3.32E−03
cg08078058	19	41,307,140	EGLN2; EGLN2	1stExon; Body	Island	0.924	0.861	0.804	0.057	0.062	4.85E−03	4.30E−04
cg12642693	X	50,556,917	SHROOM4	1stExon	Island	0.593	0.412	0.136	0.277	0.180	4.37E−03	7.94E−03
cg18110168	X	50,557,055	SHROOM4	TSS200	Island	0.696	0.466	0.116	0.350	0.230	5.52E−03	8.65E−03
cg05152874	X	50,557,122	SHROOM4	TSS200	Island	0.653	0.470	0.134	0.336	0.182	3.00E−03	5.57E−03

CHR chromosome, UTR untranslated region.

Table 3

The genes whose methylation levels gradually decreased from ANT to PVTT tissues.

Probe ID	CHR	Position	Gene name	Gene region	CpG island	Methylation levels			Mean difference		P value
						PVTT	HCC	ANT	HCC-ANT	PVTT-HCC	HCC vs. ANT	PVTT vs. HCC
cg24411043	1	161,008,297	TSTD1	Body;	N_Shore	0.156	0.231	0.292	−0.061	−0.074	0.0082	0.0085
cg27180868	1	161,008,750	TSTD1	5′UTR	Island	0.110	0.186	0.280	−0.094	−0.077	0.0076	0.0053
cg24161057	1	161,008,811	TSTD1	TSS200	Island	0.084	0.156	0.252	−0.096	−0.072	0.0046	0.0062
cg12426196	1	161,008,826	TSTD1	TSS200	Island	0.080	0.162	0.250	−0.087	−0.082	0.0029	0.0028
cg12141030	3	44,803,447	KIAA1143; KIF15	TSS1500; Body	Island	0.084	0.141	0.333	−0.192	−0.057	3.03E−5	0.0010
cg07303143	3	44,803,452	KIAA1143; KIF15	TSS1500; Body	Island	0.058	0.131	0.365	−0.234	−0.073	0.0001	0.0072
cg19273773	7	102,790,112	NAPEPLD	TSS1500	Island	0.031	0.088	0.215	−0.128	−0.056	1.1E−5	0.0028
cg08904363	7	102,790,119	NAPEPLD	TSS1500	Island	0.057	0.113	0.212	−0.099	−0.057	0.0002	0.0056
cg23882019	8	23,082,951	TNFRSF10A	TSS1500	Island	0.073	0.142	0.260	−0.118	−0.068	0.0027	0.0032
cg22843797	8	23,082,961	TNFRSF10A	TSS1500	Island	0.097	0.184	0.323	−0.138	−0.087	0.0013	0.0014
cg26530341	8	23,083,353	TNFRSF10A	TSS1500	S_Shore	0.428	0.594	0.758	−0.163	−0.166	0.0004	0.0073
cg26001125	22	24,823,050	ADORA2A	TSS1500	S_Shelf	0.352	0.461	0.561	−0.100	−0.109	0.0024	0.0077
cg08025954	22	24,823,455	ADORA2A	TSS200	S_Shelf	0.267	0.361	0.491	−0.130	−0.093	0.0001	0.0012
cg12793123	22	24,823,554	ADORA2A; ADORA2A	1stExon; 5′UTR	S_Shelf	0.335	0.416	0.484	−0.068	−0.081	0.0056	0.0012
cg05014660	22	44,576,654	PARVG	TSS200		0.076	0.151	0.223	−0.072	−0.074	0.0025	0.0075
cg01311181	22	44,577,181	PARVG	5′UTR		0.105	0.173	0.379	−0.206	−0.067	0.0018	0.0065

CHR chromosome, UTR untranslated region.

The key genes with gradually altered methylation

Previous studies demonstrated that methylation of promoter, gene body and enhancer regions always led to alternation of gene expression. We further performed association analysis between mRNA expression and DNA methylation in all the 33 samples from the 11 HCC patients [15].

The methylation levels of seventeen genes were increased from the ANT to PVTT. Among the seventeen genes, four genes (PAK1, ZC3H3, NDRG2 and SLC9A3R2) were detected to be associated with changes in mRNA expression. The positive correlations were found in PAK1 and ZC3H3, and the negative correlations were found in SLC9A3R2 and NDRG2 (Fig. 4, Supplementary Table 2). We can also observe significant overall survival difference between high ZC3H3 expression group and low ZC3H3 group in the TCGA data set (HR = 1.5, P = 0.018, Supplementary Fig. 2). For NDRG2 gene, the high expression group had a better overall survival outcome, even though the difference was not statistically significant (HR = 0.72, P = 0.064, Supplementary Fig. 3).

Fig. 4

The methylation status of (A) PAK1, (B) NDRG2, (C) SLC9A3R2 and (D) ZC3H3 in promoter regions or in body regions (blue: ANT; purple: HCC tissues; red: PVTT tissues; blue rectangle: CpG island region). The expression levels of (E) PAK1, (F) NDRG2, (G) SLC9A3R2 and (H) ZC3H3 from our inhouse data set. The correlation analysis between methylation levels and expression levels in (I) PAK, (J) NDRG2, (K) SLC9A3R2 and (L) ZC3H3.

Among the six gradually hypermethylated genes, we found the expression levels of TNFRSF10A were significantly increased in HCC tissues and PVTT tissues when compared to ANTs (Fig. 5A). We could also find the expression of TNFRSF10A in PVTT tissues were slightly higher than that in HCC tissues, but the difference was not significant (Fig. 5B). The Pearson correlation analysis demonstrated a significantly negative correlation between the expression levels and methylation levels of TNFRSF10A (Fig. 5C, Supplementary Table 2). In other words, the hypomethylation of TNFRSF10A caused the increased expression in HCC. This finding implied that TNFRSF10A might contribute to the venous invasion process in HCC and was regulated by the methylation alteration. We could also observe the increased expression of TSTD1 in HCC tissues and PVTT tissues, even though no significantly statistical difference was detected (Supplementary Fig. 4).

Fig. 5

(A) The methylation status of TNFRSF10A in promoter region in different samples (blue: ANT; purple: HCC tissues; red: PVTT tissues; blue rectangle: CpG island region) (B) The expression levels of TNFRSF10A from our RNA-seq data set. (C) The Pearson correlation analysis shows the significantly negative correlation between methylation levels and expression levels (r = −0.4107, P = 0.0176).

Replication of the five aberrantly methylated genes in PVTT tissues in an independent cohort by using Sequenom EpiTYPER assays

To validate the above-mentioned five differently methylation genes (SLC9A3R2, ZC3H3, NDRG2, PAK1 and TNFRSF10A), we conducted the Sequenom EpiTYPER assays in an independent set which contained 48 HCC patients companied with PVTT. The location of the five genes and candidate amplicons were listed in the Supplementary Table 3. Three genes were located at the promoter regions and two genes were located at the body regions (Fig. 6A and Supplementary Fig. 6A). The candidate amplicons contained or closed to the differently methylated CpG sites in the epigenome-wide DNA methylation analysis, and the methylation levels CpG sites in the amplicons are significantly correlated with each other (Supplementary Fig. 5). The results from the validation cohort were consisted with the results in the epigenome-wide analysis.

Fig. 6

(A) Schematic diagram of the amplicons’ location in the five genes. MassARRAY-based validation of differential methylation of the candidates (TNFRSF10A, SLC9A3R2 and ZC3H3). The methylation levels of the CpG sites in (B) TNFRSF10A promoter region, (C) SLC9A3R2 gene body region, (D) ZC3H3 gene body region. Yellow rectangle: Amplicon; Red dot: Failed genotyping; Blue dot: successful genotyping. * P < 0.05.

For TNFRSF10A gene, the methylation levels of all the successfully genotyped CpG sites (except for CpG 3.4) in HCC tissues and PVTT tissues were significantly lower than that in ANTs (Fig. 6B, Supplementary Table 6). More importantly, we detected the methylation levels of CpG 5 were further decreased in the PVTT tissues when compared to HCC tissues (P < 0.05, Fig. 6B). Compared to the HCC tissues, we also found the methylation levels of CpG 7 tend to be further decreased in the PVTT tissues and the difference was marginally significant (P = 0.0756).

In the four hypermethylated genes, the methylation levels of almost all the CpG sites were significantly increased in both HCC and PVTT tissues when compared to ANTs (Supplementary Table 6). We further analyzed the methylation differences between the HCC and PVTT tissues. Several CpG sites in the amplicons of ZC3H3 and SLC9A3R2 presented significantly higher methylation levels in the PVTT tissues (P < 0.05, Fig. 6C and D). In addition, the mean methylation levels of ZC3H3 gene presented the gradually increased from the ANT to PVTT tissues (Fig. 6D). However, we failed to detect the methylation difference in NDRG2 and PAK1 between the PVTT and HCC tissues (Supplementary Fig. 6B and C). All the patients in our discovery cohort had liver cirrhosis, but 6 patients in our validation cohort had no liver cirrhosis. Considering liver cirrhosis might have some impacts on methylation patterns, we further validated our findings after excluding those six patients and found similar results in TNFRSF10A and SLC9A3R2 genes (Supplementary Fig. 7). Interestingly, Interestingly, we found more CpG sites showed significant differences between the HCC and PVTT after excluding those 6 patients in ZC3H3 gene (Supplementary Fig. 7).

DAC and rh-TRAIL synergistically suppress the proliferation and induce apoptosis in HCC cell lines

TNFRSF10A, also named TNF-Related Apoptosis-Inducing Ligand (TRAIL) Receptor 1, can be activated by rh-TRAIL and further transduces cell death signal and induces cell apoptosis. It has been developed as the target for monotherapies or combination therapies with existing chemotherapeutic drugs [18]. Based on our above finding, the expression level of TNFRSF10A presented significantly negative correlation with its methylation level. It implied that the expression levels of TNFRSF10A could be up-regulated by treatment with epigenetic drugs, such as DAC. Therefore, we speculated there might be a combination therapy effect of TRAIL and DAC on HCC.

Firstly, we treated two HCC cell lines (Huh7 and SK-Hep-1) with DAC to verify whether the expression of TNFRSF10A was regulated by DNA methylation. We found the expression levels of TNFRSF10A were increased in both Huh7 and SK-Hep-1 after DAC treatment (Fig. 7A and Fig. 7B). Then, the different combination groups of two drugs were performed (Fig. 7C). The rh-TRAIL and DAC combination treatment had significantly enhanced efficacy to suppress both Huh7 cells and SK-Hep-1 cells (Fig. 7D and Fig. 7E). The CompuSyn software was used to assess the synergistic effect, and all the CI values of different combination groups were below 1 (Fig. 7D and Fig. 7E). Because rh-TRAIL could introduce apoptosis in cancer cells, we further carried out apoptosis assay to assess the combination effects on two HCC cell lines. In SK-Hep-1 cell line, both DAC and rh-TRAIL could induce obvious apoptosis (24.0% and 17.2%, respectively). The apoptosis rate induced by the combination of DAC and rh-TRAIL drugs (39%) were significantly higher than that treated with single drug (Fig. 7F and Fig. 7G). More interestingly, DAC displayed little effect on cell apoptosis in Huh7 cell line when compared to control group (control group vs. DAC group, 4.5% vs. 4.9%, Fig. 7H and Fig. 7I). However, we found the combination group (38.4%) significantly induced more apoptosis than the rh-TRAIL alone group (23.4%) (Fig. 7I). In addition, we stained the viable cells with crystal violet after drugs treatment for three days and found the viable cell count was least in the combination group (Fig. 7J).

Fig. 7

The expression of TNFRSF10A was increased after DAC treatment in (A) SK-Hep-1 cell line and (B) Huh7 cell line. (C) Diagram for the different combination groups. DAC combined with rh-TRAIL synergistically suppressed HCC cell viability. (D, E) SK-Hep-1 cells and Huh7 cells were pre-cultured in 24-well plates (20000 cell/well) and then treated with indicated doses of DAC and rh-TRAIL for 72 h. MTS was conducted to assess the impact of DAC, rh-TRAIL and the combination on two cell lines. The CI values, calculated by the CompuSyn software, <1 indicated synergism. (F, H, G, I) SK-Hep-1 cells and Huh7 cells were pre-cultured in 6-well plates to 20–25% confluence and treated with drugs for 48 h. PI/Annexin V apoptosis assay of SK-Hep-1 (F, G; DAC 10 μM, rh-TRAIL 25 ng/ml) and Huh7 (H, I; DAC 15 μM, rh-TRAIL 20 ng/ml) demonstrated the combination of DAC and rh-TRAIL could significantly induce more apoptosis than the two drugs alone (P < 0.05). (J) SK-Hep-1 cells and Huh7 cells were pre-cultured in 6-well plates to 20–25% confluence and treated with drugs for 72 h. Then, the attached cells were washed with PBS and stained with crystal violet staining. Three different independent experiments were performed.

Discussion

In our study, we analyzed the genome-wide methylation profiles of PVTT to identify the aberrant methylation in PVTT tissues and the potential differences between the PVTT tissues and its primary HCC tissues. We observed decreased global DNA methylation in PVTT tissues and HCC tissue compared to normal liver tissue. The biological function analysis reveals the aberrantly methylated genes in PVTT tissues were significantly associated with embryonic development and morphogenesis. The pathway analysis based on aberrantly methylated genes in PVTT tissues was interestingly enriched in cell adhesion molecules pathway and regulation of actin cytoskeleton pathway. To identify more reliable genes which regulated by DNA methylation and involved in PVTT formation, we analyzed the genes whose methylation levels were gradually increased or decreased from ANTs, to HCC tissues and to PVTT tissues. Strikingly, both the methylation analysis and transcriptome analysis based on the hypomethylated genes indicated that ECM-receptor interaction pathway and focal adhesion pathway were significantly involved in the process of venous invasion. The similar phenomenon was also detected in colorectal cancers (CRC). The genes in ECM remodeling pathway were silenced in CRC and were regulated by comprehensive DNA methylation [19]. Finally, we also found several genes in which significantly methylated alterations were detected associated with expression change. The current study highlights the impact of DNA methylation dysregulation on the HCC initiation and venous invasion and provided several genes which are worthy to explore in the future research.

Recently, an inspiring study which compared the methylation landscapes of single circulating tumor cells (CTCs) and CTC clusters in breast cancer was published [20]. Compared to single CTCs, CTC clusters features a higher ability to seed metastasis [21]. Gkountela et al. found the specially hypomethylated status of binding sites for stemness-associated transcription factors in CTC clusters, paralleling embryonic stem cell biology [20]. We could regard the PVTT as a special pattern of CTC cluster and also found that the aberrant methylated genes in PVTT tissues were also significantly associated embryonic development and morphogenesis in GO analysis.

The top enriched pathway based on the eleven genes annotated in the hypermethylated CpG sites in PVTT tissues was the regulation of actin cytoskeleton pathway. The results of the pathway analysis emphasized the potential central role of aberrant methylation status in the regulation of actin cytoskeleton pathway during the process of PVTT formation. Actin, an important component of the cytoskeleton, is closely related to epithelial-mesenchymal transition (EMT) [22]. Reassembly of the actin cytoskeleton is essential for invasive cell migration [23]. Previous studies have demonstrated the significance of actin cytoskeleton dysregulation to the tumorigenesis and metastasis in colorectal cancer [24], gastric cancer [25] and renal cell cancer [26]. Recently, Peng et al. [27] reported that the actin cytoskeleton remodeling drove EMT transition in HCC cells and was critical in cell density-dependent induction of EMT. The epigenetic alterations of the genes involved in this pathway provide an alternative mechanism for dysregulation of actin cytoskeleton in HCC. In addition, parts of identified genes in this pathway, like FGFR2 and PIK3CD, have been proven to be methylated-differentially expressed genes in liver disease (e.g. nonalcoholic fatty liver disease) and liver cancer [28], [29].

We identified six genes which the methylation levels were lower in the HCC tissues and even lower in the PVTT tissues. Companied with the hypomethylation status, the expression levels of TNFRSF10A were increased in HCC tissues and PVTT tissues. In addition, we confirmed this result in our independent validation cohort. This finding suggested TNFRSF10A might contribute to the process of HCC initiation and portal vein invasion. TNFRSF10A, a member of the TNF-receptor superfamily, plays an important role in the regulation of cell apoptosis process. It binds with TNF-related apoptosis-inducing ligand (TRAIL) and introduces apoptosis in cells. The dysregulation of this gene is significantly associated with tumorigenesis [30]. Even though TNFRSF10A is one of key factors for inducing apoptosis in cancer cells, it is overexpressed in various tumors [30]. In colon cancer patients, van Geelen et al. [31] found that the expression of TNFRSF10A is higher in colon cancer tissues when compared to normal tissues and the high TNFRSF10A expression was associated with worse disease-free and overall survival. Lim et al. [32] reported that the TNFRSF10A could trigger the protective autophagy in TRAIL-resistant HepG2 cells and lead to the insensitivity to the TRAIL. In addition, previous study indicated DR5 rather than DR4 (TNFRSF10A) contribute to the lysosomal permeabilization and further induce apoptosis in liver cancer cell lines [33]. Therefore, the role of TNFRSF10A in HCC, especially in metastasis, is worthy to be clarified in the future research. Moreover, because TNFRSF10A is regulated by the DNA methylation, it is reasonable to speculate the DNA methylation transferase inhibitor can enhance the efficacy of human recombinant TRAIL in HCC treatment. In the study of small cell lung carcinomas, Vitaliy et al. reported combination of DAC, a DNMT inhibitor, with a histone deacetylase inhibitor could sensitize small cell lung carcinomas (SCLC) cells to TRAIL-introduced apoptosis [34]. In vitro experiment, we could observe DAC and rh-TRAIL synergistically suppressed the proliferation and induced apoptosis in HCC cell lines. Our study provides potentially promising drug treatment strategy for HCC.

As for the genes whose methylation status were gradually increased, four genes (PAK1, ZC3H3, NDRG2 and SLC9A3R2) were identified with expression alteration. In the validation cohort, the gradually altered methylation changes were detected in ZC3H3 and SLC9A3R2. For SLC9A3R2 gene, previous studies have proven that this gene is involved in the colon cancer [35] and breast cancer [36]. In addition, the aberrant methylation status of SLC9A3R2 in cell free-DNA could be a potential biomarker for lung cancer [37]. However, the study of both SLC9A3R2 and ZC3H3 in liver cancer is lacked. More interestingly, the two genes were regulated by the aberrant methylation of gene body regions. The detailed mechanism is worthy of further exploration.

In addition, detecting the methylation status of circulating tumor DNA (ctDNA) has been applied in clinical practice to diagnose early HCC [38]. In the current study, we explored the specific methylation patterns of PVTT and several aberrant methylated genes in PVTT and validated parts of our results in the validation cohort. Those genes might serve as promising ctDNA candidates to diagnose HCC with the existing of PVTT and help to draw a more individual treatment strategy.

The main limitation of this study is the limited patient number in the epigenome-wide methylation analysis. We therefore conducted a screen without using the common genome-wide significant threshold (10−8). Because of this reason, parts of false positive and false negative DMPs would be identified. Larger scale studies are encouraged in the future to identify the methylation patterns of PVTT tissues. In addition, the potentially important role of TNFRSF10A gene in PVTT formation has been validated in the expanded samples cohort. In the in vitro experiment, we proved that the DAC and rh-TRAIL could synergistically suppress the proliferation and induce apoptosis in HCC cells, and proposed a new strategy for HCC treatment. However, we merely provided the evidence in vitro. The in vivo experiment and the more detailed mechanism are expected in the future research.

In summary, our findings confirmed that DNA methylation alterations were significantly accompanied with the PVTT formation and provided us novel insights into this process. The pathway analysis based on differentially methylated genes in PVTT tissues was interestingly enriched in regulation of actin cytoskeleton pathway, focal adhesion and ECM-receptor interaction. We also reported several crucial genes (eg. TNFRSF10A, SLC9A3R2, ZC3H3) which was regulated by the DNA methylation and might participate in the PVTT formation. More interestingly, we found the combination of DNA methylation inhibitor (DAC) and rh-TRAIL might be a potential treatment strategy for HCC.

Materials and method

Patients and samples

Eleven HCC patients with PVTT from the Shanghai Eastern Hepatobiliary Surgery Hospital were included in epigenome-wide DNA methylation profiling analysis. For the independent validation cohort, we collected 144 samples from the forty-eight HCC patients. HCC, ANT tissues and PVTT tissues were collected and immediately frozen in liquid nitrogen and stored at −80 °C until DNA and RNA were extracted. All patients were diagnosed as HCC according to the TNM staging system of the American Joint Committee on Cancer (AJCC 8th edition). The patient who received pre-operative treatment for HCC was excluded from this study, because the pre-operative treatment might have some impact on the methylation patterns. The clinical and pathological information of all the patients was collected and presented in Table 1.

This study protocol was approved by the Clinical Research Ethics Boards of the Shanghai Eastern Hepatobiliary Surgery Hospital, and informed consents were obtained from all included patients.

DNA extraction, bisulfite conversion and methylation analysis

According to the manufacturer’s instructions, genomic DNA was extracted from HCC, ANT and PVTT tissues using the QIAmp DNA Mini Kit (QIAGEN). The concentrations of extracted DNA were measured by a NanoDrop 2000 (Thermo, Wilmington, USA). Genomic DNA (500 ng–1 μg) was converted by sodium bisulfite according to the manufacturer’s protocol of the EpiTect Fast DNA Bisulfite Kit (QIAGEN, Hilden, Germany). Then, 500 ng of each bisulfite-treated DNA sample was used to analyze the epigenome-wide methylation by The Infinium Human Methylation450 (HM450) BeadChip (Illumina, San Diego, CA, USA). Normalization and filtration of raw methylation matrices were performed by using R software (version 3.4.3; https://www.r-project.org/) with the Chip Analysis Methylation Pipeline (ChAMP) package [39]. As a result, the remaining 431,362 probes were utilized for data analysis. Paired T-test or Wilcoxon matched-pairs signed-ranks test was used to determine whether there was a difference in methylation.

HCC cell lines and cell culture

Two HCC cell lines, Huh7 and SK-hep1, were selected for decitabine treatment. The cell lines were purchased from Type Culture Collection of the Chinese Academy of Sciences (Shanghai, China) (Huh7 Cat# TCHu182, SK-hep1 Cat# TCHu109). Cell lines were cultured in DMEM (Gibco, Cat# C11995500BT) with 10% fetal bovine serum (FBS) (Cellmax, Cat# SA102.02) and maintained at 37 °C in 5% CO2.

Decitabine treatment

Decitabine (DAC), one of the DNMT inhibitors, was purchased from Targetmol (Cat # T1508). HCC cells were pre-cultured to 10–25% confluence and then cultured with medium containing DAC for three days at different concentrations (0 μM, 5 μM and 10 μM). The DAC was diluted in dimethylsulfoxide (DMSO) (Sigma, Cat# D2650). We refreshed the medium for three days. The experiments were repeated for 3 times in both HCC cell lines.

RNA extraction and real-time-quantification PCR (qPCR)

Total RNA was extracted after DAC treatment for three days using RNA-Quick Purification Kit (ESscience Cat# RN001). Then, 1 μg RNA was reversely transcribed to cDNA using Hifair™ III 1st Strand cDNA Synthesis Kit (Yeasen, Cat# 11123ES10). qPCR was performed with Hifair™ qPCR SYBR Green Master Mix (Yeasen, Cat# 11201ES08). We selected GAPDH as a normalizing control. The qPCR primers for TNFRSF10A, forward primer, 5′-TGTTGCATCGGCTCAGGTTGTG-3′ and reverse primer, 5′-GTGGACAGCGAGTCTGCGTTG-3′; for GAPDH, forward primer, 5′-GTGAAGCAGGCGTCGGA-3′ and reverse primer, 5′-AGCCCCAGCGTCAAAGG-3′. The 2−ΔΔCT was used as the calculation method to analyze the expression levels of TNFRSF10A.

MTS assay and synergism evaluation

HCC cells were seeded at 20,000 cells/well in 24-well plates and cultured overnight. Then, the cells were treated with rh-TRAIL (0, 2.5, 5, 10, 25, 50 ng/ml; CAT# TP750004; OriGene Technologies, Rockvile, MD) or DAC (0, 5, 10, 20 μM) alone and in combination for 72 h. MTS method was conducted to assess the effects of drugs on cell viability using CellTiter 96 AQueous One Solution Cell Proliferation Assay (G3580; Promega, Madison, WI) according to the manufacturer’s protocol. After incubating for 3 h, we further transferred the media with MTS to a 96-well plate and measured the absorbance at 490 nm using spectrophotometry (Multiskan GO, Thermo Scientific). The combination index (CI) was calculated by Compusyn software (ComboSyn, Inc.), and a CI value <1 represents synergism [40].

PI/annexin V apoptosis assay

HCC cells were pre-cultured to 20–25% confluence in 6-well plates. For SK-Hep-1 cell line, the cells were treated with rh-TRAIL (25 ng/ml) or DAC (10 μM) alone and in combination for 48 h. For Huh7 cell line, the cells were treated with rh-TRAIL (20 ng/ml) or DAC (15 μM) alone and in combination for 48 h. The blank control group was treated with DMSO. Following the drugs treatment for 48 h, both floating and trypsinized attached cells were collected and used for apoptosis detection. The apoptotic cells were detected using Annexin V-FITC/PI apoptosis kit (Multi Sciences, CAT# AP101-100-kit) according to the manufacturer’s instructions and analyzed by flow cytometry (BD Biosciences).

Crystal violet staining

HCC cells were pre-cultured to 20–25% confluence in 6-well plates and treated with rh-TRAIL or DAC alone and in combination as described above for 72 h. We refreshed the medium which contained drugs for three days. After treatment, we washed the attached cells with PBS and stained with crystal violet solution (Beyotime Biotechnology, cat#C0121) to visualize the viable cells.

RNA-seq data analysis

The corresponding HCC tissues, including HCC tissues, ANT and PVTT tissues, had been performed RNA-seq, and the data were published in our previous study [15]. We re-ran our transcriptome data, and the Pearson correlation analysis was used to assess the correlation between methylation status and expression levels.

The heatmap and the unsupervised hierarchical clustering analysis were analyzed using MORPHEUS (https://software.broadinstitute.org/morpheus/), an online heatmap making and analysis website.

Gene bioinformatics and primers design

The information of the five genes analyzed in our study were obtained from the UCSC genome database (http://genome.ucsc.edu/). Two online websites, the Gene Expression Profiling Interactive Analysis website (GEPIA, http://gepia.cancer-pku.cn) and cBioportal (http://www.cbioportal.org/), which analyze bioinformation based on the public data sets, were used to get more gene data. The location information of the five genes were summarized in the Supplementary Table 3.

We used the DBCAT software (http://dbcat.cgm.ntu.edu.tw/) to scan the corresponding sequence and identify the candidate CpG islands in the five genes. The location information of the five candidate amplicons were listed in the Supplementary Table 3. The primers used for quantitative methylation analysis of the five genes were listed in the Supplementary Table 4.

Mass array quantitative methylation analysis

The Sequenom EpiTYPER assay was used to detect the methylation levels of the converted genomic DNA. Three main steps were included in this process: PCR amplification, SAP cleanup and T cleavage. After SAP cleanup, we usually conducted the gel electrophoresis to verify the efficiency of PCR amplification. T cleavage was further performed after confirming high efficiency of amplification. The products were then transferred to a SpectroCHIP® array and analyzed on the MassARRAY® Analyzer 4 platform.

Gene ontology (GO) and KEGG pathway analysis

Functional and pathway enrichment of the aberrant methylated gene were performed by GO analysis and KEGG pathway analysis using DAVID Bioinformatics Resources 6.7 (http://david.abcc.ncifcrf.gov/) [41].

Statistical analysis

We used the SPSS 20.0 version to perform the statistical analysis. Student’s t-test and χ2 were used to compare the variables between the discovery cohort and validation cohort. Paired two-tailed Student’s t-test was conducted to compare the methylation differences among the ANT, HCC tissues and PVTT tissues. Pearson correlation analysis was used to assess the methylation correlations between each CpG sites.

Ethics approval and consent to participate

This study protocol was approved by the Clinical Research Ethics Boards of the Shanghai Eastern Hepatobiliary Surgery Hospital (EHBHKY2019-K-007), and informed consents were obtained from all included patients.

Consent for publication

Not applicable.

Availability of data and materials

The datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request.

Competing interests

The authors declare no conflict of interest.

Funding

We are greatly indebted to the subjects enrolled in our study. This work was supported by the (2016YFC0906400), (81400656, 81872297, and 81874059), Zhejiang province analysis and test technology project (2018C37062), the (2016XZZX002-05), Medical leading project of Shanghai Science and Technology Commission Foundation (17411960300).

CRediT authorship contribution statement

Xiaoxiao Fan: Methodology, Formal analysis, Writing - original draft, Writing - review & editing, Visualization. Yirun Li: Formal analysis, Writing - original draft. Xin Yi: Methodology. Guoqiao Chen: Software. Shengxi Jin: Methodology. Yili Dai: Formal analysis. Bin Cui: Resources. Binghua Dai: Conceptualization, Supervision, Funding acquisition. Hui Lin: Conceptualization, Supervision, Funding acquisition. Daizhan Zhou: Conceptualization, Methodology, Writing - review & editing, Visualization, Supervision, Funding acquisition.