Introduction

Recent projections suggest hepatocellular carcinoma (HCC) will be the fifth leading cancer-related cause of death in the US in 2026, accounting for ~31,000 deaths.1 Despite advances in cancer prevention and treatment, HCC remains associated with poor prognosis, and 5-year survival remains low (22%).1,2 The growing population burden of metabolic dysfunction, including obesity and type 2 diabetes mellitus, contributes to HCC occurrence, with epidemiological studies consistently showing obesity as an independent risk factor, even after accounting for viral hepatitis and alcohol use.3–7 The increase of harmful environmental exposures like viruses, mycotoxins, hepatocarcinogens, tobacco smoke, and workplace carcinogens also contribute to HCC incidence.4

In the United States, high mortality is due in part to HCC is typically being undiagnosed until advanced stages, when treatment options are few, with limited efficacy. Notably, HCC arises almost exclusively in the context of chronic liver disease (CLD), most commonly following a prolonged, clinically silent course of hepatic injury that culminates in cirrhosis, the end stage of CLD. This extended and often asymptomatic phase represents a substantial and underutilized opportunity for early intervention. While early-stage HCC is potentially curable, with 5-year survival exceeding 70%,5,8 and surveillance is associated with earlier detection and improved outcomes,6,9,10 existing screening strategies lack sufficient uptake, underscoring the urgent need for novel, reliable approaches to detect HCC earlier along the disease continuum.

Current HCC surveillance guidelines rely heavily on semiannual or annual ultrasound, often combined with serum alpha-fetoprotein (AFP) levels.8 However, ultrasound and AFP serum levels have suboptimal sensitivity and specificity for early-stage disease, limiting their utility as stand-alone risk stratifying markers.11,12 Many risk stratification tools with a variety of approaches have been proposed to make up for this.

For example, etiology-based risk stratification has its shortcomings, as not all HCC patients present with a recognized risk factor, like cirrhosis, and reliance on clinical factors can lead to both overdiagnosis and missed cases.13 Moreover, many models were developed in cohorts dominated by chronic viral hepatitis infection (ie., hepatitis B virus (HBV), hepatitis C virus (HCV)).14–16 With effective antiviral therapies reducing the burden of HBV and HCV, and the increasing prevalence of other CLD etiologies, especially metabolic dysfunction-associated steatotic liver disease (MASLD), these models may be less relevant or accurate; the epidemiology of HCC is changing.

Multi-omics approaches have begun to address this gap. For instance, proteomic approaches, including a measure of circulating biomarker aberration (Hepatoscore-14) and phosphoproteomic subclassifications, have identified molecular alteration subgroups linked to prognosis.17–22 Additionally, polygenic risk scores (PRS) such as PRS-HPC and NAFLD-PRS may be useful as predictors of progressive hepatic damage, including fibrosis, cirrhosis and HCC.23(p90),24–27 Composite models that integrate tumor-associated proteins with demographic factors offer improvements, but remain imperfect.16–18,28–32 For example, the GALAD score (age + sex + AFP ng/mL + AFP-L3% + DCP ng/mL) or newer Oncoguard® Liver test (sex + AFP ng/mL + 3 DNA methylation markers in HOXA1, TSPYL5, and B3GALT6), now in post-market surveillance, improve diagnostic accuracy compared to ultrasound alone. However, given their inclusion of tumor-derived signals (AFP in both), these approaches provide a snapshot of the molecular landscape of HCC, but their utility for early risk stratification remains uncertain, as abnormal AFP levels are not present consistently leading up to tumorigenesis. Additionally, few models incorporate environmental factors, and performance may vary across etiologies and racial/ethnic groups.29

Mounting evidence of the last two decades supports epigenetic modifications as substantial contributors to tumor initiation and progression, making these modifications particularly promising biomarkers for early detection.30–35 DNA methylation changes at CpG dinucleotides, which contribute to regulation of chromatin structure and gene expression, are among the earliest alterations in carcinogenesis and are readily measured in blood.36,37 Profiling liver disease in tissue- or leukocyte-derived DNA using Illumina Infinium Beadchip arrays has revealed differential methylation in CpG sites located near or within genes involved in oxidative stress and CLD, including MASLD, metabolic dysfunction associated steatohepatitis (MASH), and cirrhosis.38–46 Methylation-based classifiers for HCC detection, such as HepaQ and HelioLiver, have demonstrated higher sensitivity and specificity than AFP or GALAD.47,48 Yet, these assays have inherent limitations: both were designed with and for Asian populations and HelioLiver was designed with a control sample of liver disease patients, thus performance can vary based on underlying disease etiology, and limited validation in diverse populations restricts generalizability. Additionally, Illumina Infinium Bead Chip arrays, commonly used for these classifiers, cover <5% of CpG sites, leaving the majority of the genome unexplored.

However, the broad challenge for all of these tools is that although DNA methylation changes regularly arise in cancer, they vary widely across tissues and individuals, including race- and ethnicity-specific differences, and are often consequences of deregulation in tumors, rather than promoters of tumorigenesis. Moreover, deregulation of DNA methylation often occurs due to comorbid conditions such as various neurological, autoimmune, metabolic, infectious and cardiovascular diseases, confounding the identification of specific HCC associations. To effectively utilize methylation as a clinically reliable risk signal, markers must be both consistent across individuals and mechanistically linked to disease predisposition, regardless of time before disease onset.

Imprinting control regions (ICRs) are uniquely suited to meet these criteria for effective prognostic markers. Defined by parent-of-origin–dependent DNA methylation, ICR methylation is established in the gametes and pre-implantation zygote, resisting postfertilization epigenetic reprogramming. Therefore, these methylation profiles are perpetuated through germ layer specification and tissue differentiation and normally maintained throughout the life course. ICRs are also sensitive to environmental exposures during their window of establishment in early development, such as parental smoking, alcohol use, endocrine-disrupting chemicals, and non-chemical stressors (eg. physical, social). Such disruption of ICR methylation can then propagate through the germline, transferring altered gene expression to the next generation, and influence disease susceptibility over the life course,4,49,50 making them robust biomarkers that can reveal mitotically inherited epigenetic disruptions existing prior to tumorigenesis. Because imprinted genes regulate growth, metabolism, and cellular homeostasis, their dysregulation due to aberrant ICR methylation plausibly contributes to cancer initiation.30–32,51,52 By evaluating methylation at ICRs in blood, we are targeting loci that are both biologically constrained and clinically accessible; loci for which deviation from the expected pattern may serve as a strong indicator of mitotically inherited susceptibility, applicable for risk stratification.

Previously, aberrant ICR methylation or loss of imprinting (LOI) has been implicated in HCC and other cancers. Well-characterized loci including the ICRs of growth inhibitory long non-coding RNAs (lncRNA) H19, KCNQ1OT1, and MEG3 and imprinted growth promoters IGF2 and DLK1 have key roles in tumorigenesis and metabolic disfunction, with LOI of these loci reported in several cancers, including HCC, demonstrating key roles in tumorigenesis and metabolic dysfunction.33,34,53–56 More recently, methods such as quantitative chromogenic imprinted gene in-situ hybridization (QCIGISH) have demonstrated the feasibility of measuring imprinting dysregulation in cancer.35,52,57,58 Despite their importance, until the recent characterization of ~1,488 putative ICRs, only a small fraction of this set has been confirmed and thoroughly studied in human diseases.59

Herein, we used whole genome bisulfite sequencing (WGBS) to comprehensively profile the 5-methylcytosine landscape in DNA derived from peripheral leukocytes to identify previously uncharacterized HCC-associated DMRs. We augmented these with published methylation array data, and through cross referencing with previously reported putative ICRs, we identified 97 HCC-associated putative ICRs exhibiting differential methylation between HCC cases and otherwise healthy individuals.

Materials and Methods Participants and Samples

The characteristics of 39 primary HCC cases and 87 otherwise healthy controls are summarized in Supplementary Table 1. CpG methylation data used in these analyses were derived from two case-control comparisons: 1) WGBS of leukocyte-derived DNA from n = 10 primary HCC cases (mean age 66 ± 6 years) consented at the Lineberger Comprehensive Cancer Center at UNC Chapel Hill who underwent HCC resection surgery, and n = 51 otherwise healthy controls (mean age 63 ± 10 years) consented in primary care or accompanying neurology patients at Duke Health, 2) DNA methylation data using a custom Illumina methylation array (10,000 CpG sites at ICRs) from 29 incident HCC cases and 36 otherwise healthy controls identified using a combination of the rapid case ascertainment system of the North Carolina Cancer Registry, Duke Gastroenterology and the Duke Cancer Institute. The age of the latter participants ranged from 63 to 68 years. Data collection and analyses were approved by Institutional Review Boards at UNC Chapel Hill, Duke University School of Medicine and North Carolina State University.

WGBS Identification of DMRs Associated with Primary HCC

WGBS was performed in peripheral blood of HCC patients (n = 10) and otherwise healthy individuals (n = 51). Libraries for the Illumina NovaSeq 6000 Sequencing System were prepared from bisulfite converted DNA derived from blood using previously described methods.59 To generate long, high-quality reads from both ends of DNA fragments, sequencing was performed using S2 flow cells with 150bp paired-end (PE) reads. The reads were aligned to the hg38 reference genome using BSMAP60 as implemented in the MOABS framework.61 The average number of reads mapped was 132.7 million (75.6 M–258.1 M) and the average alignment percentage was 74.3% (65.1%–90.0%). The average percent of paired reads excluded as duplicates was 4.1% (0.20%–5.6%). Quality control metrics at the BAM level were obtained with samtools flagstat (v1.15.1). Both the control and HCC samples showed mapping rate ~100% and a high proportion of properly paired reads (97.0% and 95.2%, respectively), with very few singleton reads (≤4%) (Supplementary Table 2). The resulting data were processed, and MOABS software was used to call methylation sites using the “MCALL” function.61 Differentially methylated regions (DMRs) were defined using a 300 bp window, requiring at least 4 CpGs per window to all have a minimum credible methylation difference of 10%, in the same direction, between groups. CpG sites were included only if they had a minimum coverage depth of three reads (--minDepthForComp = 3). Approximately 24 million CpG sites were analyzed. Sequencing depth and average methylation levels were consistent across coverage bins, suggesting that both case and control datasets are of similar and sufficient quality for downstream methylation analysis (Supplementary Figures 1 and 2). The WGBS data has been deposited in NCBI’s Gene Expression Omnibus and are accessible through GEO Series accession number GSE302608 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE302608).

We additionally transitioned from the unadjusted and pooled MOABS DMR analysis to a multifactor generalized linear model using the DSS (Dispersion Shrinkage for Sequencing) R package to perform a stratified analysis with age as a covariate. DNA Methylation counts at each locus were obtained using MOABS mcall (the MOABS function was preferred for methylation calling as the initial read alignment performed with BSMAP was not suitable for Bismark pipeline). Spatial smoothing was applied to recover the statistical power at WGBS locus level with bsseq package BSmooth. The generalized linear model was fit using the formula ~ HCC_Status + Age. We then tested for differential methylation against the HCC status coefficients by performing the Wald test on each locus using DMLtest.multifactor to investigate the disease-specific association. The local inclusion parameter for DMR calling was p < 0.05, and false discovery rate was controlled by spatial constraints of minCG = 4, dis.merge = 100bp, and pct.sig = 0.05.

To confirm that methylation changes were not driven by genomic alterations, we performed copy number analysis on the WGBS data using CNVkit (v0.9.10).62 A pooled reference was constructed from control samples. Analysis was conducted using a moving average window of 100 kb to mitigate the inherent variance in bisulfite sequencing coverage.

To post-hoc estimate the power to identify significant DMRs, 500 simulations were run using base R software. Within each, 500 regions were simulated, with 25% randomly selected to be true DMRs. Within each region, a varying number of CpG sites was simulated using a beta-binomial approach to stay methodologically in line with MOABS software, and account for both variation due to coverage and between participants. Replicating our WGBS DMR calling requirements, DMRs were defined as requiring at least 4 CpGs per simulated region window each with a minimum methylation difference of 10%. Comparing the proportion of DMRs declared to the preset DMRs, this simulation gave a power estimation of ~72%.

Gene Annotation

Genomic coordinates were annotated with gene information using the RefSeq gene annotation track for the human genome (GRCh38/hg38). Refseq data were downloaded from the UCSC Genome Browser database (https://hgdownload.soe.ucsc.edu/goldenPath/hg38/database/).

Literature Review of Published HCC-Associated Loci

To comprehensively identify publications reporting altered methylation in association with liver damage, we used the National Library of Medicine’s PubMed database, with search terms including HCC and pre-HCC conditions – liver cirrhosis, non-alcoholic fatty liver disease (NAFLD), MASLD, fibrosis, MASH, steatosis, hepatic fat levels, HBV, HCV, and alcoholic liver disease – and “methylation”, “differential methylation” or “DMR[s]”, and “biomarkers”. With narrowing hits based on methods measuring differential methylation levels between etiologies using Illumina Infinium BeadChip arrays (27k, 450k, EPIC, EPIC v2.0), WGBS, or pyrosequencing, this search yielded n = 66 articles; n = 64 were primary research articles, including meta-analyses, and n = 2 were reviews. From the 64 primary research articles, 19 articles were excluded for missing chromosomal coordinate information. We identified regions reported as significantly differentially methylated between the cases of liver disease (HCC or other degenerative liver disease conditions that can lead to HCC, as listed above in the search terms) and non-cases. All chromosomal positions were converted to the hg38 build. A total of 45 publications were included and are summarized in Supplementary Table 3.38–41,43–46,55,63–97

Determining HCC-Associated Putative ICRs

Among the 1,519 DMRs associated with HCC in mixed leukocytes, we sought to identify overlaps with putative ICRs to facilitate early detection for interventions or triage. We used two sources of published putative ICRs. Briefly, as described in Jima et al (2022), 1,488 putative ICRs were identified in humans by performing WGBS of DNA derived from tissues representing the three germ layers. Additionally, Akbari et al (2022) determined candidate ICRs using publicly available Oxford Nanopore sequencing datasets. Using the NanoMethPhase package, that study identified 143 putative ICRs.98,99

To determine previously reported DMRs associated with HCC risk that are potential biomarkers prior to HCC onset, we determined direct overlaps between the 1,519 HCC-associated DMRs and putative ICRs from Jima et al and Akbari et al, resulting in 81 HCC-associated putative ICRs.

To ensure a comprehensive list of HCC-associated putative ICRs, we additionally intersected previously published loci with the compendium of putative ICRs reported by Jima et al and Akbari et al Analysis of the 45 published studies revealed 3,275 unique regions associated with liver pathology, 673 of which are single dinucleotide CpG sites. Of the 3,275 DMRs/CpG sites previously reported, 59 were identified using WGBS blood-derived DMRs and the remainder from methylation arrays (Supplementary Table 3). The overlap of these regions with the putative ICRs identified 16 additional putative ICRs that contain loci previously implicated in HCC development via published literature, totaling to n = 97 HCC-associated putative ICRs.

Replication via Human Imprintome Methylation Array Platform

Methylation levels in mixed leukocyte-derived DNA isolated from the whole blood of 29 HCC cases and 36 otherwise healthy controls were measured on the Human Imprintome Methylation array. The array, an Illumina Infinium HTS iSelect Custom BeadChip, is a targeted DNA methylation platform designed using the repertoire of 1,488 putative ICRs across the human genome characterized by Jima et al as target regions.59,100 This array contains 22,819 probes (704 control probes and 22,115 CpG probes, of which 10,438 CpG probes target unique CpG sites). Among the CpG probes, 9,757 probes passed design criteria for 1,088 of the 1,488 candidate ICRs; see Carreras-Gallo et al for additional details.100 Additional cgBackground probes enabled background-normalization. The IDAT files containing the raw intensity data for each probe on the microarray chip have been deposited in NCBI’s Gene Expression Omnibus and are accessible through GEO Series accession number GSE303108 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE303108).

The raw IDAT files produced from the array were processed using the openSesame function in the Sesame R package.101 The Sesame preprocessing code “0CDB” applies a non-linear dye bias correction and background subtraction using the noob method (Norm-Exp out of band signal), which deconvolves out-of-band fluorescent signal from true signal. The beta values are calculated by the ratio of the methylated signal to the total signal (methylated + unmethylated) for each probe, representing the methylation fraction as a value between 0 and 1. These processing steps produced a beta measurement for each probe, along with a corresponding Sesame signal detection p-value quality score. Signal detection p-values were assigned using the pOOBAH method (P-value with out-of-band (OOB) array hybridization).101

Probes that did not map to a CpG site in the Imprintome were discarded, and any remaining probes were also discarded if their signal detection p-value exceeded 0.2 for more than 10% of samples measured. For cases where multiple probes targeted the same CpG site, averaged beta measurements were used.

A series of Student’s t-tests were performed to test for significant differential methylation of beta values between cases and controls (null hypothesis of µ = 0, two-tailed, α = 0.05, unpaired, unequal variance assumed). Resulting p-values were adjusted based on false discovery rate for the number of ICRs tested. CpG sites were reported as significant based on an adjusted p-value < 0.05. Additionally, CpG beta values were aggregated if within the same ICR, producing beta measurements for each ICR site of each participant for ICR-level case-control comparison. See Supplementary Table 4 for all specialized R packages used.61,101–109

To estimate the power of this study post-hoc, we first estimated the true standard deviation of the data by calculating the median pooled standard deviation of methylation beta values for the subset of CpGs with coverage over the 70 ICRs. Standard deviations were calculated for each group before pooling (σpooled = 0.081). We set the minimum effect size to methylation beta = 0.06 (Cohen’s d = 0.73), marginally below the minimum of 0.1 to define a DMR from WGBS analysis. Statistical power for multiple comparisons within the 70 ICRs was estimated using the Sample Size Calculations for Two-Sample Microarray Experiments from the ssize.fdr package, with parameters (FDR = 0.05, power = 0.8, m = 70, pi0 = 0.5).110 The proportion of nulls (pi0) was set to 50% because this subset of ICRs was already found differentially methylated in the WGBS/previously reported DMR analysis. This test assumes a constant sample size, so we conservatively specified the smaller sample size between groups (n = 29) to calculate the minimum estimated power across all comparisons (βall > 0.73, function:ssize.twoSamp). This analysis confirmed that we had sufficient sample size for statistical power of >73%.

Results WGBS Identification of DMRs Associated with HCC in Peripheral Blood Leukocytes

Using WGBS on mixed leukocytes from 10 HCC cases and 51 controls, 1,519 DMRs were significantly associated with HCC (FDR q-value range 0–0.000486) (Supplementary Figure 3). Of the 1,519 blood-derived DMRs, the majority (73%; 1103/1519) displayed hypomethylation, and 27% (416/1519) exhibited hypermethylation associated with HCC (Figure 1). The range of the absolute values of the methylation difference (regardless of direction) is 0.100–0.454.

Figure 1 Distribution of hyper- and hypomethylation among the blood-derived DMRs. In peripheral blood between HCC cases and controls, the majority of DMRs (1103/1519) exhibited hypomethylation associated with HCC.

To assess functional significance, the identified DMR positions were overlapped with publicly available ENCODE chromatin immunoprecipitation sequencing (ChIP-seq) peak datasets for H3K4me3 and H3K27ac – regions enriched for active promoter and enhancer marks – in a B-lymphoblastoid cell line (GM12878) and liver tissue (ENCODE Experiment IDs ENCFF416SIL and ENCFF805YRQ). Overlap with these peak regions indicates localization of DMRs to putative regulatory elements. Additionally, CpG island enrichment in the DMRs was assessed using the UCSC Genome Browser CpG Islands track (hg38) (https://genome.ucsc.edu/cgi-bin/hgTrackUi?hgsid=3564790247_NWaSAOctGsVQ4vD5dafUzrKX1tk4&db=hg38andc=chr7andg=cpgIslandSuper) (Supplementary Table 5). In the liver tissue line, 17.2% (261/1519) of the HCC-associated DMRs were enriched for H3K4Me3 marks, and 16.0% (240/1519) were enriched in the B-lymphoblastoid cell line. Approximately the same proportion of DMRs were enriched for H3K27ac in the liver tissue and B-lymphoblastoid cell line (14.4% (218/1519) and 12.8% (195/1519), respectively). Of these, 72 DMRs were enriched in both liver and B-lymphoblastoid, suggestive of conservation of essential epigenetic regulation and that blood methylation captures HCC-relevant enhancers. Additionally, 26% (401/1519) of the DMRs overlapped with at least one CpG island.

The age-stratified generalized linear model analysis of differential methylation identified 955 DMRs associated with HCC status, accounting for approximately two-thirds of the 1,519 pooled DMRs. The substantial overlap with the 1,519 DMRs identified in the primary analysis suggests that the modest (~3 year) age difference between groups does not meaningfully confound observed associations. Due to potential deflation of significance in stratified/covariate-adjusted analyses, given the small sample size, further analyses were performed on the entire set of 1,519 HCC-associated DMRs to maximize sensitivity.

CNVkit was used to examine copy number variation in each of the ten HCC cases for which WGBS was performed (Supplementary Figure 4A–J). Analysis of log2 copy ratios across all autosomes demonstrated a stable diploid baseline centered at 0.0, indicating the observed depth at the sites in the HCC samples matches the baseline read depth calculated from the control samples. No large-scale aneuploidy or consistently amplified/deleted segments were detected.

The nearest genes of the 1,519 blood-derived DMRs were used in Qiagen Ingenuity Pathway Analysis (IPA) to determine if there was over-representation of HCC-related pathways in the set of DMR-associated genes. Using bedtools suite and RefSeq hg38 gene annotations, 1,321 unique nearest genes mapped to the DMRs, with several DMRs mapping to the same gene. After removing four genes not annotated in IPA (“LOC730100”, “LOC124903068”, “LOC102724219”, “LOC124905221”), 1,317 genes were processed for pathway enrichment. The top five “Diseases and Disorders” pathways are cancer (n = 838 genes, FDR adjusted p-value range 3.82E-03 – 1.05E-42), organismal injury and abnormalities (n = 878 genes, 3.82E-03 – 1.05E-42), endocrine system disorders (n = 748 genes, 3.82E-03 – 1.06E-39), hereditary disorder (n = 317 genes, 3.82E-03 – 4.1E-37), and dermatological diseases and conditions (n = 638 genes, 3.82E-03 – 2.94E-33). The top five hepatotoxicity pathways are liver hyperplasia/hyperproliferation (n = 412 genes, FDR adjusted p-value range 1.00E00 – 5.51E-11), liver inflammation/hepatitis (n = 39 genes, 1.00E00–4.3E-03), hepatocellular carcinoma (n = 147 genes, 1.00E00 – 7.2E-03), liver steatosis (n = 40 genes, 1.00E00 – 1.81E-02), and liver regeneration (n = 8 genes, 2.80E-01 – 3.46E-01). The top 5 significant canonical pathways are Protein Kinase A Signaling (n = 31/396 genes in the pathway, FDR adjusted p-value 1.49E-05), Hepatic Fibrosis Signaling Pathway (n = 30/416, 2.18E-04), 3-phosphoinositide Biosynthesis (n = 19/215, 2.18E-04), mRNA 3 Prime End Processing Signaling Pathway (n = 12/102, 2.18E-04), and Superpathway of Inositol Phosphate Compounds (n = 20/241, 2.18E-04). All canonical pathways are available in Supplementary Table 6.

Defining a Comprehensive List of HCC-Associated Putative ICRs

For comparison of our findings with results from published data, we created a compendium of previously reported DMRs or differentially methylated CpG sites found in multiple source tissues, including circulating cell-free DNA (ccfDNA) from whole blood or plasma and liver tissue biopsies. The overlap of these regions with the putative ICRs identified 16 additional putative ICRs that overlap differential methylation previously implicated in HCC development via published literature. Of the 16 ICR/DMR regions, 15 of the ICRs were first identified in liver tissue, and one was first identified in peripheral blood associated with alanine aminotransferase (ALT) levels.45,55,68–70

With 81 HCC-associated blood-derived DMRs from our WGBS and 16 DMRs previously published as associated with HCC or pre-HCC risk, we compiled a set of 97 HCC-associated DMRs that overlap with putative ICRs, of which 92 are unique Imprintome regions (Supplementary Table 7). This includes 21 regions co-located with previously identified ICRs (see Supplementary Table 8 for the list of previously identified ICRs). All but five ICRs from Akbari et al that overlapped with the DMR compendium were also reported in the 1,488 Imprintome ICRs (Figure 2).

Figure 2 Intersection of HCC-associated DMRs with Putative ICR targets. Our compendium of HCC-associated regions, containing the 1,519 blood-derived DMRs from WGBS and from our literature analysis, contained 4,794 DMRs. The 143 Nanopore-detected ICRs were reported by Akbari et al,98 and the 1,488 putative Imprintome ICRs were reported by Jima et al.59 For some ICRs, there were multiple overlapping DMRs, thus we note the number of unique ICRs included.

Table 1 HCC/Liver Disease-Associated Differentially Methylated Regions Detectable in Mixed Leukocytes, cfDNA and Liver Tissue

See Table 1 for further details of the 97 regions and Figure 3 for a schematic workflow for identifying putative ICR targets. These 97 regions represent a subset of putative ICRs that have exhibited differential methylation between HCC cases and controls in an accessible surrogate tissue and are therefore potential biomarkers for early detection prior to HCC onset due to the establishment of imprinting regulation in utero.

Figure 3 Schematic workflow for identifying putative ICR targets. We combined the 1,519 blood DMRs with the DMRs reported in previous literature from mixed blood and tissue sources. Using a Perl script, we intersected those with the 1,488 putative ICRs previously described111 as well as 143 putative ICRs reported by Akbari et al98 to find 97 direct overlaps, of which 81 DMRs came from our WGBS of HCC case/control blood methylation levels.

Determining Molecular Pathways of Liver Disease-Associated Targets

We used Qiagen IPA to analyze the genes overlapping or closest to the 97 putative ICR regions. The two most significant metabolic gene pathways are associated with HCC risk: PTEN signaling (FDR adjusted p = 9.12E-03) and IGF-1 signaling (p = 1.28E-02). Additional significant (p<0.05) HCC risk-associated pathways include growth hormone signaling, signaling by MET, toll-like receptor signaling, p38 MAPK signaling, PDGF signaling, signaling by IGF1R, acute phase response signaling, and the folate signaling pathway (Supplementary Table 9). Dysregulation of these pathways is associated with cellular proliferation and several cancers. The 10 most significant canonical pathways are shown in Table 2.

Table 2 Top Pathways Associated with HCC-Associated Genes by Ingenuity Pathway Analysis

Table 3 Independent Validation of Differentially Methylated ICRs Using the Imprintome Array

Table 4 Significantly Differentially Methylated ICRs at the CpG Level in the Human Imprintome Array

Replication of the Putative ICRs in an External HCC Case-Control Analysis

To externally replicate the identified ICR targets, an independent population of n = 29 cases and n = 36 otherwise healthy controls was used to measure DNA methylation on the Human Imprintome Methylation Array, which contains probes for CpGs within 1,088 of the 1,488 putative ICRs.100 Because this array was designed to the ICRs reported by Jima et al, five of the 97 regions that were reported by Akbari et al are not included.98 In addition, 12 ICRs lacked associated probes after quality control filtering, limiting this replication to 70 ICR sites overlapping HCC-associated DMRs (containing a total of 1,049 CpG sites). We observed significant differential methylation associated with HCC at both the CpG site-level and ICR-level. At the ICR level, 14 ICRs were significant after adjustment for false discovery rate (Table 3). These regions are in proximity to the transcripts LINC02774, RPN1, MCHR2, GRB10, PEG10, MEST, KCNQ1, MEG3, ZNF597, KCNJ18, ZNF331, BLCAP, ELK1, and CMC4. Nine of the 14 genes are known imprinted genes.

When considering CpGs individually, 46 ICRs had at least one significantly differentially methylated CpG site between the cases and controls (Table 4). Together, these results provide supporting evidence for the association of these ICRs with HCC risk.

CpG methylation summary plots were created for each ICR to visualize the difference between case and control methylation levels. Figure 4 displays the mean methylation of CpG sites in cases and controls for ICR 454 (GRB10, previously identified ICR), ICR 1417 (ELK1, putative ICR), ICR 1192 (BLCAP|NNAT, previously identified ICR), and ICR 1438 (NAP1L6P, putative ICR), as examples. Each ICR displays methylation levels consistent with hypomethylation in HCC cases. Full plots and CpG methylation summaries can be found at: https://mollyrbio.shinyapps.io/HCCIMPlots/.

Figure 4 Differential CpG methylation for four ICRs. Visualization of differential methylation between HCC cases and controls for two previously identified ICRs (A and B) and two putative ICRs (C and D). Mean methylation beta values for each CpG site across cases and controls were estimated as detailed in methods and plotted for each ICR (https://mollyrbio.shinyapps.io/HCCIMPlots/). X-axis: Mean methylation beta value, a continuous measure (0–1) representing the proportion of methylated CpGs at the site. Y-axis: Illumina array CG identifier number. * Adjusted p-value <0.05.

Discussion

We conducted WGBS in mixed leukocytes to characterize HCC-associated alterations in the 5mC landscape. This analysis was complemented by a comprehensive literature review of previously reported liver disease-associated DMRs to provide broader biological context. We focused specifically on loci overlapping putative ICRs. Given their established methylation stability across tissues and cell types and persistence throughout the life course, diseased-associated aberrant methylation at these regulatory regions has potential to serve as early-risk biomarkers.59,98 Findings were subsequently evaluated in an additional case-control study using leukocytes, accessible cells that could serve as surrogates for inaccessible liver cells, given the conservation of ICR methylation. Ultimately, these efforts were designed to identify epigenetic alterations detectable in leukocytes that may inform early risk assessment for HCC, by identification of susceptible individuals otherwise appearing healthy.

The WGBS of leukocytes derived from a cohort of HCC cases and otherwise healthy controls identified 1,519 blood-derived DMRs, 81 of which coincide with ICRs, both demonstrated and putative.

Functional enrichment analysis of these regions with histone marks suggests essential epigenetic regulatory function and that methylation analysis of blood-derived DNA captures HCC-relevant enhancers. Enrichment for H3K4me3, which is associated with active promoters, indicates potential transcriptional effects from methylation changes on active genes, with these effects conserved across tissues. Similarly, DMR overlap with H3K27ac suggests involvement of active enhancers with cell type-specific regulatory roles. Enrichment at CpG islands – regions typically exhibiting tightly regulated and stable methylation in blood – further supports the interpretation that these DMRs reflect biologically meaningful regulatory variation rather than shifts in cell composition.112

IPA analysis of genes nearest to the 1,519 HCC DMRs revealed enrichment across three major functional themes: (i) fibrosis and cirrhosis signaling, (ii) metabolic and obesity-related pathways, and (iii) canonical oncogenic signaling cascades implicated in HCC progression. Notably, several of these pathways include known or putative imprinted genes, reinforcing the biological relevance of our ICR-focused approach.

Among oncogenic signaling networks, Protein kinase A (PKA) signaling emerged as significant. Multiple PTPN family members were represented in this pathway. PTPN members regulate signal transduction and metabolic processes, have been implicated in obesity-related phenotypes, and are generally considered tumor suppressors in HCC, where reduced expression is associated with poor prognosis.113–118

Growth factor signaling genes were also prominent, including IGF2 and IGF2R. IGF2, a well-established imprinted gene, is frequently overexpressed in HCC and functions as a key mitogenic driver of tumor growth, while IGF2R, which modulates IGF2 availability, is often downregulated, further contributing to cancer progression.119–121 The representation of these genes within enriched pathways supports a model in which altered methylation near imprinted growth regulators may contribute to dysregulated proliferative signaling. Also, notably, emerging evidence suggests that epigenetic dysregulation in HCC intersects with metabolic reprogramming to promote tumor progression, including mechanisms such as lactylation- and acetylation-mediated activation of glycolytic and angiogenic pathways.122,123 In this context, hypermethylation at the IGF2R|AIRN ICR may have implications beyond growth factor signaling, as IGF2R has been linked to lysosomal function and iron homeostasis, raising the possibility that its dysregulation contributes to ferroptosis resistance. Future studies integrating iron-related biomarkers (eg. ferritin, hepcidin) and ferroptosis regulators such as GPX4 could help determine whether altered imprinting at this locus promotes metabolic vulnerabilities that support therapeutic evasion.

Pathways related to metabolic dysfunction and fibrosis were additionally represented, including hepatic fibrosis signaling. LEP (leptin), a hormone that regulates body weight, appeared among the significant genes, consistent with crucial role of leptin in obesity, cirrhosis, and MASLD – conditions that elevate HCC risk.124–126 Several imprinted genes are known regulators of leptin production or signaling, further linking imprinting biology to metabolic and inflammatory pathways implicated in hepatocarcinogenesis.127–131 High leptin levels are also linked to increased HCC risk due to its ability to promote tumor cell proliferation, invasion, and angiogenesis through various signaling pathways including JAK2/STAT3 (both of which were enriched pathways) and PIK3/AKT (a gene among the 1,519 DMRs).132,133

Collectively, these findings indicate that genes proximal to the HCC-associated DMRs converge on interconnected signaling networks governing growth regulation, metabolic homeostasis, and fibrotic progression – core biological processes underlying HCC development. Incorporating previously reported liver-disease associated DMRs, we identified 16 additional regions overlapping putative ICRs, bringing the total to 97 HCC-associated ICR-linked loci. Of these, 21 overlap established imprinted domains, including regions within or proximal to IGF2R|AIRN, MEST|MESTIT1, H19, KCNQ1|KCNQ1OT1, MEG3|DLK1 (including the intragenic IG-DMR and the MEG3-DMR), GRB10, PEG10, CHRNE, PRMT2, DHFR, DLGAP2, INPP5F, OVOS2, and ZNF597. These loci broadly cluster into growth factor signaling, metabolic regulation, and canonical tumor suppressor/oncogene networks – processes central to hepatocarcinogenesis.

Several of the identified imprinted regions converge on the IGF signaling axis. As discussed above, dysregulation of IGF2 and its receptors is a recurrent feature of HCC, with aberrant activation of the axis reported in a subset of early tumors.121,134,135 The presence of IGF1R and IGF2R within these putative ICR regions reinforces the relevance of altered imprinting and growth factor regulation in HCC development. Other loci, including MEST, H19, MEG3, DLK1, PEG10, and GRB10, have established roles in tumor growth, metabolic regulation, or both, further supporting the biological plausibility of imprinting-linked epigenetic disruption in liver disease and cancer progression.136–142 Indeed, single ICRs can regulate multiple tumor-associated genes; for example, the IG-DMR of MEG3 regulates not only MEG3 and DLK1, but also is responsible for establishing the imprinting status of >80 genes in the 14q32.2 imprinted domain.143PRMT2, INPP5F, and OVOS2 are also novel oncogenes.142,144–146

Additional genes proximal to putative ICRs, such as KCNK15, have not been established as imprinted, but have been associated with HCC prognosis in prior transcriptomic analyses.147 We identified a hypomethylated DMR with KCNK15, and hypothesize that the methylation at this locus may regulate gene expression and alterations can serve as a marker of disease progression.

IPA of genes nearest to these 97 DMRs revealed significant enrichment of pathways previously implicated in HCC, notably PTEN signaling and IGF-1 signaling. These findings mirror the broader IPA results from the 1,519 DMRs, again highlighting convergence on growth factor signaling and PI3K/AKT pathway regulation. PTEN signaling, a critical negative regulator of PI3K/AKT activity and cellular metabolism, has also been shown to influence the expression of multiple imprinted genes.148,149 A 2021 study in mice showed that when Pten is absent, the expression of several imprinted genes, including Igf2, Plagl1 (Zac1), Cdkn1c, Dlk1, Mkrn1, Magel2, and Dlx5, changes in mouse embryonic stem cells, embryoid bodies, and cardiomyocytes.150 The second top pathway, IGF-1 signaling, mediates cell proliferation, survival, migration, and blocks apoptosis in HCC, as detailed above.151,152

This may provide a potential mechanistic bridge between imprinting dysregulation and oncogenic metabolic reprogramming.

Despite the small sample size (N = 10 cases), we determined the reproducibility of our findings in mixed leukocyte-derived DNA of an additional HCC case-control study using a novel methylation array targeting a majority of ICRs comprising the Human Imprintome.59,100 Of the 97 HCC-associated putative ICRs, 46 (47%) had significantly associated CpGs, and further 14 of these ICRs, several of which had not previously been linked to HCC, were significant across their spans. While we do not report the direction of methylation change between cases and controls for these ICRs, both hyper- and hypomethylation could coincide with imprinting dysregulation. Whether this dysregulation happens via hypo- or hypermethylation is dependent on the method of imprinting regulation at the locus (ie. parent-of-origin specific methylation at the gene promoter, long-range chromatin mediated interactions such as enhancer blocking, insulator function, or lncRNA-mediated silencing). Future studies are needed to determine which mechanisms are in operation at each candidate ICR.

Additional imprinted genes implicated in HCC and liver disease were not represented among the 97 regions, but warrant consideration in future biomarker refinement. For example, PLAGL1 (ICR-404), involved in cell-proliferation control, was not included in our DMR compendium. However, this gene has been shown to have abnormal transcription in cell lines derived from HCC patients with HCC, and its murine ortholog Zac1 has been demonstrated to directly coordinate imprinted gene networks involved in MASLD pathophysiology.153–155 While genes such as PLAGL1 may play mechanistic roles in hepatocarcinogenesis, our biomarker strategy prioritizes loci demonstrating stable germline imprinting detectable in blood, as these are more suitable for non-invasive risk stratification. These excluded loci may show variable methylation, partial imprinting, or tissue specificity, hence not meeting our stringent criteria.

Other studies have investigated blood-based DNA methylation marks in HCC. Hernandez-Meza et al67 used Illumina HumanMethylation450K arrays (485,000 CpGs, <2% of known CpG islands) to profile a cohort of normal liver tissues, cirrhotic tissues, dysplastic nodules, and HCC tissues. They found an increased proportion of hypermethylated samples in the progression from cirrhotic tissue (<1%) to dysplastic nodules (≥25%) to HCC (>50%). The study also observed an inverse correlation between DNA methylation and gene expression for TSPYL5, KCNA3, LDBH, and SPINT2 (all p < 0.001).67 Interestingly, KCNA3 encodes a potassium voltage-gated channel, similar to KCNK15, which was identified in our study. Lubecka et al found hypermethylation at LSP1 in HCC relative to cirrhosis controls, a gene we have identified in our list of putative ICRs.91 While each of these panels show promise for differential methylation at CpG sites as biomarkers, they are not selective for ICRs, so aberrant methylation identified may be consequential to tumorigenesis and progression, and thus useful in early risk stratification. Thus, while multiple studies have identified epigenetic marks with high sensitivity and specificity to detect HCC in blood-derived DNA, such marks are difficult to interpret.15,16 Additionally, many previous studies were limited by small sample sizes lacking population diversity, inferred methylation levels using novel techniques, and/or lack of target validation.

Our findings identify epigenetic biomarkers that are present across the life course and independent of cirrhosis status, which can be incorporated into existing risk models, such as GALAD. The biomarker panel used as part of GALAD has recently shown 82% specificity and 62% sensitivity when retrospectively evaluated for cirrhosis patients within 12 months before clinical diagnosis of HCC.156 Merging epigenetic risk markers with protein markers of progressive liver disease can be evaluated for their incremental predictive performance.

Like other HCC studies discussed, we are limited by HCC-patient and control sample sizes (WGBS: n = 10 cases, 51 controls, replication case control: n = 29 cases, 36 controls). Because HCC remains a rare—though highly lethal—cancer despite rising incidence, assembling larger well-characterized cohorts remains inherently challenging. However, post-hoc power analyses show sufficient power to detect meaningful differences in both the WGBS discovery of DMRs and the Human Imprintome Methylation Array to replicate these DMRs (72% and 73% at our given sample sizes, respectively). Also, although several of the 97 putative ICRs did not replicate in the external validation, despite being previously implicated in HCC risk, likely due to population heterogeneity, these findings pave the way for experimental validation of the 97 regions and identifying potential mechanisms that could be targeted for dietary or pharmacological interventions. Additionally, our studies were performed in case-control settings with methylation measurements made in samples obtained at diagnosis. Although ICR disruption occurs early in development, methylation changes could still arise concurrently with HCC onset. From contemporary cohort studies being assembled such as All of Us or Cohorts for Environmental Exposure and Cancer Risk (CEECR), it will soon be possible to determine if methylation differences identified here precede cancer onset or arise secondarily from tumorigenesis.

Nonetheless, these results provide comprehensive evidence for association of a set of ICRs with HCC risk. These regions are thus strong candidates for further validation in diverse cohorts, to be developed as early-risk stratifying HCC biomarkers. Given the stability of ICRs, temporally and spatially, such biomarkers would be applicable to DNA from accessible cell types, such as mixed leukocytes, and at any time during the life course, providing efficient, non-invasive, early detection.

Conclusions

Here, we address the critical need to identify early-detection epigenetic markers of HCC risk that may allow providers to facilitate interventions while the liver has the capacity to regenerate and, thus, reduce overall HCC incidence. For public health screening, it is critically important that DNA methylation marks be consistent across individuals, and within individuals, both across age and cell types, so that accessible DNA sources, such as blood and saliva, are informative as proxies for affected tissues.

We have achieved a pivotal step towards this goal of creating a clinically applicable screening panel, with our key finding of 97 putative ICRs with demonstrated differential methylation between patients with HCC and otherwise healthy individuals in DNA from a clinically accessible source, blood. When replicated in larger prospective studies, these marks could be developed into a biomarker panel for early detection in any setting, including primary care.

Abbreviations

HCC, Hepatocellular carcinoma; ICRs, imprinting control regions; D