Genome-wide methylation profiling of the different stages of hepatitis B virus-related hepatocellular carcinoma development in plasma cell-free DNA reveals potential biomarkers for early detection and high-risk monitoring of hepatocellular carcinoma

Genome-wide methylation profiling with plasma cell-free DNA from different stages of hepatitis B virus-related hepatocellular carcinoma development

Data mining resulted in 37,610,900; 37,072,952; 35,016,215; 33,286,609 and 33,002,633 raw reads for the samples from the HC, CHB, LC, HCC, and NSCLC patients, respectively. When these were mapped to the hg19 genome reference, 19,083,142; 18,285,110; 16,389,041; 14,422,110 and 14,894,363 reads were obtained, respectively, with a mapping rate of approximately 50% (Figure 2A, Additional file1: Table S4). Notably, very high background levels of methylation were observed for the HC plasma cfDNA.

Association of plasma cell-free DNA methylation disturbances with the progressive stages of hepatocellular carcinoma development

We used hierarchical clustering to trace the progression of methylation disturbances during the carcinogenic process for the HC, CHB, LC and HCC samples. The analysis revealed the closeness of aberrant DNA methylation patterns among different groups of samples (Figure 2B), suggesting the possible sequential direction of the development of epigenetic alterations for the samples. Moreover, we estimated the relationships between the aberrant methylation patterns in the various stages of HCC development by the pairwise overlap of the peak bases from each stage, which revealed the similarities between HC and CHB as well as between LC and HCC (Additional file2: Table S5). In addition, we investigated the DNA methylation densities surrounding the transcription start site (TSS), which is a crucial region that affects the gene expression of regulatory mechanisms, and we observed an accumulation of aberrant DNA methylation events as the disease stage progressed from HC to CHB to LC and, finally, to HCC (or NSCLC) (Figure 2C). These results suggest a stage-related stepwise process in which the gene methylation regulatory elements were sequentially affected during the carcinogenesis process. Collectively, the above observations indicate that the aberrant methylation of plasma cfDNA occurs in a dynamic, progressive and stage-responsive manner during HCC development. This creates the basis for our further analyses of HCC progression-related DNA methylation abnormalities.

Identification of differentially methylated loci among the various stages of hepatocellular carcinoma development

To reliably obtain the DNA methylation peaks, we used two methods, MACS and BALM. The peaks that were generated by these two methods were generally consistent (90 to 97% coincidence). We performed our analysis using these peaks, which totaled 226,898; 257,344; 282,203; 242,987 and 290,096 for the HC, CHB, LC, HCC and NSCLC groups, respectively (Additional file3: Table S6).

Because our long-term aim was to establish methylation markers that can be used for detecting and monitoring HCC development and progression, we were particularly interested in hypermethylated DNA alterations. To identify the stage-related DNA methylation disturbances, we investigated the alterations at each stage, including CHB, LC and HCC as well as NSCLC, which were compared to the normal background levels of the HC. The different methylation peaks that were subsequently generated are referred to as the differential methylated regions (DMRs). The most abundant DMRs were observed in the “total DMR” category, reflecting all of the hypermethylation peaks that were elevated compared with the HC (Figure 2D; Additional file3: Table S6). In addition, the “gene-related DMR” and “CGI-associated DMR” categories could be isolated, and the category that included both the CGI and gene-related DMRs became the focus of our further analyses (Figure 2D, Additional file4: Table S7).

The analysis of aberrantly methylated DMRs at different loci in association with the HCC developmental process also revealed that widespread genomic structural elements were affected, which included a variety of structural elements as defined by the UCSC Genome Bioinformatics (UCSC) (Figure 2E and Additional file5: Table S8).

Early-, middle- and late-stage differentially methylated regions and their related differentially methylated genes

To further elucidate the cfDNA methylation alterations that occur during the CHB, LC and HCC stages, we compared the DMRs that were associated with each stage on a chromosomal distribution basis against those from the HC (Figure 3A). The analysis showed that CHB is characterized by sparse and mild hypermethylation alterations, suggesting that certain HBV-induced epigenetic disturbances begin as early as this stage. However, these hypermethylation alterations obviously increased during the transition from CHB to LC, indicating the presence of more extensive epigenetic disturbances. The HCC stage was characterized by substantial hypermethylation, and a similar degree of epigenetic change was also observed for NSCLC, suggesting that such alternations are commonly associated with various malignancies in the late stages of carcinogenesis.

Provided that the introduction of DMRs was a continuous, dynamic process during the different stages of HCC development, we observed the following three categories of DMR involvement: 1) those that were common to CHB, LC and HCC, which may play important roles in the process of HCC development from CHB to HCC; because they occurred as early as the CHB stage, we referred to them as early DMRs; 2) those that were common to both LC and HCC, which may be essential for the transition from LC to HCC; we called these middle DMRs; and 3) those that were unique to HCC, which are most likely indispensible for maintaining the malignant status; we called these late DMRs. To identify the exact genes that are affected by the aforementioned early-, middle- and late-stage DMRs (the differentially methylated genes (DMGs)), we used “CGI + gene” involvement as the selection criterion and generated 240 early-, 272 middle- and 286 late-stage DMGs (Figure 3B, Additional file6: Table S9, Additional file7: Table S10 and Additional file8: Table S11).

Recent studies have reported that tissue and cancer-specific DMRs occur more frequently within CGI shores, which are regions of relatively low CpG densities that flank traditional CGIs (up to 2 kb in distance) than within CGIs themselves, indicating the importance of considering CGI shores in DNA methylome analyses[19–21]. Therefore, we generated 225 early-, 232 middle- and 259 late-stage DMGs according to CGI involvement (Additional file9: Table S12, Additional file10: Table S13 and Additional file11: Table S14), while a significantly expanded panel of 1,349 early-; 1,290 middle- and 967 late-stage DMGs were generated when CGI-shore involvement was considered (Figure 3C, Additional file12: Table S15, Additional file13: Table S16 and Additional file14: Table S17). This represents an obvious trend where CGI-associated DMGs gradually increase along with HCC progression, while CGI-shore-associated DMGs decrease, possibly suggesting that the epigenetic effects that occur during the late stages of HCC development are associated with CGI alone rather than CGI-shores.

Differential methylation profiling between hepatocellular carcinoma and non-small cell lung cancer

To locate plasma DNA methylation loci that are specific to HCC and that may serve as potential diagnostic markers, we compared the DMGs between HCC and NSCLC and found that 414 were shared by both malignancies, while 462 were specific to HCC and 178 were specific to NSCLC (Figure 3D and Additional file15: Table S18, Additional file16: Table S19 and Additional file17: Table S20). These results suggest that although carcinogenesis and tumor development may originate from different tissues, they share common DNA methylation alterations, as well as possess unique disturbances. Therefore, these differences can be used to determine the different tissue origins of carcinomas.

Functional involvement of plasma cell-free DNA differentially methylated genes associated with hepatocellular carcinoma developmental stages

We performed gene ontology and pathway analyses of the early-, middle- and late-stage HCC DMGs using the DAVID method. All of these DMGs were selected for analysis because their promoter regions were aberrantly methylated.

Methylation-specific PCR (MSP) screening, qMSP validation and gene expression of the library revealed differentially methylated genes in independent tissue DNA

Two techniques were used to assess the tissue DNA for validating the accuracy of the plasma DNAm profile library. First, 40 DMR targets that possessed the top P values were screened by MSP in 10 normal liver and 33 HCC (Additional file18: Figure S2) tissue DNA samples, generating 33 targets that were informative of HCC. Subsequently, the informative targets were investigated by qMSP in the same tissues sets that had supplemented with 29 LC tissues. As expected, a substantial portion of the targets exhibited hypermethylation in the LC, HCC or both the LC and HCC. The typical results in selected genes are shown in Figures 4A and B, with the UCSC scheme of the gene information and their quantitative methylation statues in HC, LC and HCC. These results suggest that the altered methylation patterns that were observed in the plasma are consistent with the respective liver tissues.

In addition, the expression levels of selected genes, ZNF300, SLC22A20 and SHISA7, were examined using real-time reverse transcription (RT-PCR) in the liver tissue of HC, LC and HCC. There appeared a coincidence that hypermethylated ZNF300 showed downregulation with progression of HC to LC and then to HCC (P <0.05). However, SLC22A20 and SHISA7 did not have hypermethylation coupled with gene silence pattern (Figure 4C).

Multiplex-BSP-seq validation of differentially methylated genes in independent plasma cell-free DNA

To further confirm the cfDNA methylation levels, we used Multiplex-BSP-seq on the targets that were selected by qMSP in the independent plasma cfDNA samples (Additional file19: Table S2). Of the 33 target genes that were studied, 22 were successfully amplified and deep-sequenced. The representative results are shown in Figure 5. For each gene, the methylation levels of seven to ten individual CpG sites were investigated. GAPDH and KCNV1 exhibited relatively fixed methylation patterns for the seven CpGs that were inspected, regardless of the HCC developmental stage. Therefore, both can serve as quality control measurements. In this study, the methylation levels in some of the CpGs of 3 particular genes were closely related to the developmental stage of HCC; therefore, they were predictive of the disease stage. These genes and their informative CpG sites included ZNF300 at CpG6, SLC22A20 at CpG3 and 5 and SHISA7 at CpG1, 5 and 6 (see Figure 6A, B), which could distinguish between HC + CHB and LC + HCC with a lower cut-off value, while distinguishing between HC + CHB + LC and HCC with a higher cut-off value of CpG methylation (Figure 6A). When these CpG methylation measurements were subjected to multiple univariate logistic regression analyses, they showed strong or very strong associations with stage discrimination (OR: 5 to approximately 14) in distinguishing between HC + CHB and LC + HCC or between HC + CHB + LC and HCC (Figure 6B). Taken collectively, Multiplex-BSP-seq using an independent set of plasma samples validated certain CpG methylation changes in a subset of genes that were revealed by MethylCap-seq during HCC development. This also suggests the important role of aberrant methylation in these CpG sites during HCC carcinogenesis.

Patient samples

In this study, group of patients were recruited by following criteria:

CHB (chronic HBV-related hepatitis) patient recruitment: more than 6 month history of diagnosed or suspected acute hepatitis; laboratory test evidence of HBV infection; mild to moderate increase of alanine aminotransferase (ALT) over years; no serious manifestations outside of hepatic system.

HC (healthy control) recruitment: individuals who attended routine healthy examinations; no histories of hepatitis; negative hepatitis B-related tests; no abnormal findings for hepatic biochemical tests.

HCC (HBV-related hepatocellular carcinoma) patient recruitment: pathological or clinical evidence indicating HCC diagnosis; chronic HBV related hepatitis or long-term HBV infection; normal or mildly abnormal liver biochemical examinations; no serious dysfunction in kidney, heart, lung, brain and other organs.

LC (HBV-related liver cirrhosis) patient recruitment: pathology or imaging diagnoses of cirrhosis patients; over a year of expected survival (liver function at level A or B by Child, MELD score ≤10); histories of chronic HBV infection or long-term hepatitis B infection; normal or mildly abnormal liver biochemical examinations; normal or slightly elevated serum AFP levels (≤20μg/L).

NSCLC (non-small cell lung cancer) patient recruitment: pathological or clinical diagnoses of lung cancer; no histories of hepatitis; no positive findings of HBV hepatitis infections by related blood tests; no abnormal findings in liver biochemical tests; no obvious dysfunctions in renal, heart, liver, brain or other organs.

Slash represents missing data.

DNA and RNA extraction

A peripheral venous blood sample was drawn into an EDTA tube for each subject, and the plasma sample was separated within 2 hours. The plasma was isolated by centrifugation of the whole blood at 2,500 × g for 10 min. To eliminate potential cellular debris, an additional centrifugation step was performed at 15,000 × g for 10 min[34]. The cfDNA was extracted from 800-μl aliquots of plasma using the Qiagen Blood DNA Kit (QIAGEN, Hilden, Germany). The tissue DNA was isolated from 100 mg of corresponding tissues using a conventional proteinase K/organic extraction method as previously described[35]. The extraction of total RNA from liver tissue and reverse transcription were performed as described in reference[36].

Processing of MethylCap-seq

To obtain adequate starting material for the MethylCap-seq library, 30 ng of plasma cfDNA from each participant were pooled for each group. Approximately 500 ng of DNA for each of the five libraries (HC, CHB, LC, HCC and NSCLC) was sonicated until all of the DNA fragments were within the desired size range (100 to 200 bp). The end-repair, adenosine base addition and adaptor ligation steps and methylated DNA enrichment were performed as previously described[37]. Finally, 1 μg of each product was placed in the Genome Analyzer II (Illumina, Inc., San Diego, CA) to generate 36 bp-long unpaired reads. The resulting data were processed to generate genomic methylation profiles and then uploaded to a public database (Gene Expression Omnibus: GSE54961 (2014)).

Processing of MethylCap-seq data and GO analysis

We used the BWA alignment tool[38] with the default settings to map the 36-bp unpaired reads to the hg19 human genome reference assembly[39]. After removing the PCR duplicates with Picard (http://broadinstitute.github.io/picard/), SAMtools and Picard were used to convert, sort and index the aligned data. The DMRs in the four observed groups (CHB, LC, HCC and NSCLC) for comparison with the HC were identified using two methods, MACS[40] and BALM[41], to increase the power of detection for the MethylCap-seq analysis. Methylation peaks (hypermethylated regions) were identified using MACS for the four observed groups and HC, respectively, as previously described[37]. For BALM, a dual-threshold strategy was applied to decrease the false-positive detection of the DMRs. A high-confidence threshold of 0.975 for the four observed hypermethylated regions and a low-confidence threshold of 0.95 for the HC hypermethylated regions were set. Then, the specific methylation peaks of the four observed groups were defined as the hypermethylated regions. Similarly, HC-specific methylation peaks were defined as hypomethylated regions using the reverse setting. Whole-genome methylation (methylation of each CpG) was inferred with BALM, which was processed to perform Pearson correlation analyses among all of the samples using the R software. The refSeq genes (UCSC genes) and corresponding CpG islands (CGIs) were downloaded from the table browser of the UCSC database[39]. The BED files were manipulated with Bedtools[42] in addition to assorted Perl scripts.

Screening with methylation-specific PCR (MSP) and validation with qMSP

The tissue DNA and cfDNA were converted with sodium bisulfite using an EpiTect Kit (QIAGEN, Hilden, Germany). All of the DNA (converted or unconverted) samples were immediately stored at -80°C for further analyses. Methylation-specific PCR (MSP) primer pairs were designed using online software (http://www.urogene.org/cgi-bin/methprimer/methprimer.cgi; Additional file21: Table S1). In total, 40 MSPs were performed, and the results were visualized on gels. The MSP products were cloned by T-Vectors and verified by sequencing. The in vitro methylated DNA from the HepG2 cells was obtained using the CpG (M. SssI) methyltransferase (NEB, MA,USA) treatment and used as a positive control. Water was used as a non-template control. Up to 33 gene targets were measured using methylation-specific real-time PCR (qMSP) as previously described[43]. In brief, the PCR products of the methylated and unmethylated sequences of each gene target were cloned into pGEM-T Easy Vectors (Tiangen, Beijing, China). The plasmids were diluted, and a standard curve was constructed to measure both the methylated and unmethylated targets. Real-time PCR was conducted with the Rotor-Gene Q (QIAGEN, Valencia, CA). Target gene methylation was calculated using the following formula: 100 × methylated reaction/(unmethylated reaction + methylated reaction).

Multiplex-BSP-seq

The BSP regions were situated within the regions that were assessed by the qMSP analysis. BSP amplicons that were 80 to 120 bp in size were obtained because of the fragmentation of the cfDNA template. The BSP primer pairs were designed using online software (http://www.urogene.org/cgi-bin/methprimer/methprimer.cgi; Additional file19: Table S2). Degenerate primers were utilized in the promoter regions in case there where inevitably contained CpG sites. To overcome the difficulties that are associated with the scarcity of cfDNA, Multiplex-BSP-seq was performed. In total, 35 BSP primers were divided into three panels according to their annealing temperatures. Three independent multiplex-PCR amplifications were conducted for each cfDNA sample, and the three PCR products were combined in equal amounts as determined by gel electrophoresis and Qubit for ligation with a barcode adaptor (oligonucleotide sequences © 2007-2012 Illumina, Inc. All rights reserved). A final pooling, which contained equal quantities of each participant, was performed with all of the ligation products, and 12 cycles of PCR were performed using library amplification primers to produce the templates for deep sequencing. Illumina HiSeq 2500 was used to generate 100-bp paired-end reads (Additional file22: Figure S1). The mapping of the raw reads was performed by Bismark (version 0.7.7) and Bowtie2 (version 2.0.5), and the resulting SAM data were sorted by SAMtools. The R methylKit package was utilized to assess the methylation levels at each CpG locus with default settings unless otherwise specified. The methylation level (%) that was measured at each individual CpG site was expressed as the percentage of methylated CpG versus unmethylated CpG.

Real-time reverse transcription PCR

The expression levels of selected genes were measured by real-time quantitative RT-PCR analysis. Three or more reactions were performed for each sample using a RealMasterMix (SYBR Green) kit (TianGen, Beijing, China) with gene-specific primers (Additional file23: Table S3) on a Rotor-Gene Q (QIAGEN, Valencia, CA). The RNA quantity was normalized to the ACTB content, and gene expression was quantified according to the 2^-ΔCt method.

Statistical analyses

To evaluate the assay’s potential for discriminating between different disease stages, the data that were obtained for the HC, CHB, LC and HCC samples were used to calculated the optimal cut-off value using the most significant area under the receiver operating characteristic (ROC) curve. Univariate models were used to examine the association between the gene CpG methylation and HCC stage. Multivariate models were then developed that adjusted for the most important covariates, including age and gender. The measurement data were analyzed using a one-way ANOVA. The sample means were compared with an unpaired t-test, assuming unequal variances, and all of the tests were two-tailed. P values are shown for statistically significant findings. All of the statistical analyses were performed using the SPSS statistical package (version 13.0, Chicago, IL, USA).