Chromatin accessibility module identified by single-cell sequencing underlies the diagnosis and prognosis of hepatocellular carcinoma

Abstract

BACKGROUND

Hepatocellular carcinoma (HCC) is notorious for its aggressive progression and dismal prognosis, with chromatin accessibility dynamics emerging as pivotal yet poorly understood drivers.

AIM

To dissect how multilayered chromatin regulation sustains oncogenic transcription and tumor-stroma crosstalk in HCC, we combined multiomics single cell analysis.

METHODS

We integrated single-cell RNA sequencing and paired single-cell assay for transposase-accessible chromatin with sequencing data of HCC samples, complemented by bulk RNA sequencing validation across The Cancer Genome Atlas, Liver Cancer Institute, and GSE25907 cohorts. Cell type-specific chromatin architectures were resolved via ArchR, with regulatory hubs identified through peak-to-gene linkages and coaccessibility networks. Functional validation employed A485-mediated histone 3 lysine 27 acetylation suppression and small interfering RNA targeting DGAT1.

RESULTS

Malignant hepatocytes exhibited expanded chromatin accessibility profiles, characterized by increased numbers of accessible peaks and larger physical regions despite reduced peak intensity. Enhancer-like peaks enriched in malignant regulation, forming long-range hubs. Eighteen enhancer-like peak-related genes showed tumor-specific overexpression and diagnostic accuracy, correlating with poor prognosis. Intercellular coaccessibility analysis revealed tumor-stroma symbiosis via shared chromatin states. Pharmacological histone 3 lysine 27 acetylation inhibition paradoxically downregulated DGAT1, the hub gene most strongly regulated by chromatin accessibility. DGAT1 knockdown suppressed cell proliferation.

CONCLUSION

Multilayered chromatin reprogramming sustains HCC progression through tumor-stroma crosstalk and DGAT1-related oncogenic transcription, defining targetable epigenetic vulnerabilities.

Key Words: Hepatocellular carcinoma; Single-cell multiomics; Chromatin accessibility; Enhancer-like peak; Chromatin coaccessibility

Core Tip: This study unveils multilayered chromatin accessibility dynamics in hepatocellular carcinoma. Single-cell multiomics analysis discovered malignant hepatocytes adopting increased numbers of accessible peaks and larger physical regions despite reduced peak intensity, while tumor-stroma intercellular chromatin coaccessibility maintains tumor ecological symbiosis. Functional validation identified DGAT1 as an epigenetic vulnerability: Pharmacological histone 3 lysine 27 acetylation inhibition and relating decreased chromatin accessibility suppressed DGAT1 expression, and DGAT1 knockdown further blocked tumor growth. Our findings redefine hepatocellular carcinoma progression as a chromatin-adaptive disorder governed by tumor-intrinsic epigenetic prioritization and microenvironmental crosstalk, proposing simultaneous targeting of chromatin hub genes and compensatory pathways as precision strategies.

INTRODUCTION

Hepatocellular carcinoma (HCC), as the common cancer worldwide, caused more than 900000 new cases to be diagnosed and more than 800000 deaths annually, with low 5-year survival rate due to limited therapeutic options[1]. Emerging evidence highlights chromatin accessibility dynamics as critical epigenetic regulators of HCC progression[2-4]. Chromatin accessibility, which determines the physical access of DNA to the transcriptional machinery, affects determines the binding capacity of transcription factors (TFs) and RNA polymerase recruitment, thereby shaping cellular identity and oncogenic plasticity[5].

Early chromatin studies prioritized promoter-proximal regions, where accessibility directly facilitates transcription initiation. However, bulk epigenomic studies, including The Encyclopedia of DNA Elements and Roadmap Epigenomics, revealed enhanced distal enhancer activity over promoters during hepatocarcinogenesis[6]. Traditional enhancer biology posits physical looping via mediator complexes to amplify transcriptional output, yet the regulatory hierarchy between promoter-centric and enhancer-driven control in HCC remains unresolved.

Furthermore, the heterogeneous tumor microenvironment (TME) functions as a dynamic ecosystem. Key cellular constituents, such as malignant hepatocytes, immune cells, fibroblasts, and endothelial cells, each contribute distinct mechanisms to tumor progression[7]. Bulk analyses identified global chromatin decompaction as a hallmark of HCC[3], but mask cell-subtype specificity. Single-cell assay for transposase-accessible chromatin with sequencing (scATAC-seq) now delineates epigenetic divergence between malignant and non-malignant cells at cellular resolution. Strikingly, recent work demonstrated tumor exploitation of arginine starvation to induce global chromatin compaction and impair T cell activation, further emphasis on approaches to obscure cell type-specific regulatory programs[8]. Few studies emphasize on the intercellular chromatin accessibility networks and the relating mediate tumor-stroma symbiosis in HCC. Whether similar intercellular chromatin coordination exists in HCC and how it regulates bidirectional transcription remain unexplored.

Recent advances in single-cell epigenomic technologies have revealed hierarchical chromatin regulation in other malignancies, operating through local promoter control, distal enhancer wiring, and microenvironmental coordination[3,7-13]. While chromatin accessibility alterations are increasingly recognized as central to cancer biology[11], their role in HCC remains underexplored. Specifically, the functional contribution of enhancer-driven chromatin remodeling is poorly characterized in HCC. The integrative application of single-cell RNA sequencing (scRNA-seq) and scATAC-seq enables comprehensive investigation of intricate epigenetic mechanisms involved in tumor development at cellular resolution, offering not only enhanced capability to identify key regulatory pathways driving oncogenesis but also surpassing conventional limitations of traditional cell type classification methods.

Here, we present an integrated single-cell multi-layered data of HCC that combines bulk ATAC-sequencing (ATAC-seq), bulk RNA sequencing (RNA-seq), and single cell ATAC-seq and single cell RNA-seq. This approach overcomes the limitations of bulk analyses and provides the comprehensive atlas of chromatin accessibility dynamics in the HCC TME. Our study addresses three pivotal findings. First, malignant cells prioritize chromatin accessibility quantity over quality to maintain oncogenic transcription. Second, we uncovered clinically relevant long-range chromatin interactions compensate for promoter-centric regulation in HCC, defining diagnostic and prognostic signatures. And validated DGAT1, gene regulated by chromatin accessibility hub, showed sensitive to chromatin accessibility alterations at SNU449. Knocking down DGAT1 with small interfering RNA (siRNA) showed decreased cell proliferation. Third, intercellular chromatin coaccessibility mediates tumor-stroma symbiosis. These findings establish a paradigm shift from ‘quality to quantity’ in cancer epigenetics, with direct implications for the development of chromatin-focused therapies for HCC.

MATERIALS AND METHODS

The Cancer Genome Atlas and Liver Cancer Institute related data acquisition and processing

We obtained chromatin accessibility region data from The Cancer Genome Atlas (TCGA) and analyzed under the standard protocal RNA-seq and related clinical data from the TCGA-liver hepatocellular carcinoma (LIHC) cohort were downloaded from the Genomics Data Commons portal. The TCGA-LIHC dataset comprises 361 primary tumor cases. We focused on those specifically classified as nonrecurrent primary HCC. Within this cohort, matched peri-tumoral histologically normal hepatic tissue specimens with complete transcriptomic profiles were obtained from 59 participants.

To validate TCGA-derived findings, we incorporated a multi-institutional HCC transcriptomic cohort (GSE14520[14], GSE25907[15] and GSE10186[16]). GSE14520 comprised hepatitis B virus (HBV)-associated cases from Fudan University’s Liver Cancer Institute (LCI) and Zhongshan Hospital. This retrospective cohort provided curated gene expression profiles (n = 221), matched survival annotation, and histologically normal peri-tumoral specimens (n = 209). GSE25907 included 268 HCC tumors, 243 adjacent non-tumor tissues, 40 cirrhotic tissues, and 6 healthy liver samples. GSE10186 established the consensus molecular classification into three subtypes: S1 (WNT pathway hyperactivation driven by transforming growth factor β, not β-catenin mutation), S2 (proliferative MYC proto-oncogene/protein kinase B activation), and S3 (hepatocyte differentiation). Clinical-pathological parameters with survival outcomes were extracted from supplementary repositories, while pre-processed microarray datasets enabled immediate bioinformatic interrogation without additional normalization.

scRNA-seq data collection and processing

Our single-cell analytic framework integrated hepatic transcriptional profiles from eight HCC cases across two established Gene Expression Omnibus (GEO) cohorts: GSE125449 and GSE151530)[7,17]. Raw scRNA-seq data retrieval from GEO repositories were followed by meticulous processing: (1) Leveraged pre-curated droplet matrices from originating studies; and (2) Implemented Seurat-based computational workflows (v4.0) encompassing stringent quality filters (cellular doublet exclusion, mitochondrial content thresholding) and nonlinear normalization protocols, culminating in a refined transcriptomic atlas of 15660 high-confidence cells.

To analyze the single-cell level data across the cohorts, we first removed cohort and experiment-specific effects by performing data integration by Seurat. Counts were first normalized using log normalize. These 8 samples across 2 cohorts underwent independent normalization and finding variable features using 4000 genes for the number of variable features. After anchor identification, all samples were integrated with the Seurat’s scVIIntegration in using 40 dimensions. To minimize batch effects, we further used Harmony method. Finally, we performed clustering on these corrected counts for downstream cell-type assignment. We conducted principal component analysis and constructed the shared nearest neighbor graph using 20 principal components. This graph served as the input for Louvain clustering, executed via the FindClusters function with a resolution of 0.4.

Cell typing standards: Cluster annotation employed a dual-method validation framework: Logistic regression-based marker gene detection (Seurat’s FindAllMarkers: |avg_log₂FC| > 0.3, Bonferroni-corrected P < 0.05) Lineage confirmation through pathway enrichment analysis via cluster-specific differential expression signatures. Hierarchical clustering reconciled 13 computationally derived subgroups into six consensus lineages informed by lineage-specific transcriptional signatures and ontology analyses. Raw sequence demultiplexing utilized sample-specific cellular barcodes for experimental batch tracking. Copy number variation (CNV) analysis was further applied to testify the robustness of malignant cells clustering.

Data collection and processing for scATAC-seq of HCC

Our single-cell chromatin accessibility analysis in HCC leveraged paired scRNA-seq and scATAC-seq data from eight treatment-naïve HCC specimens (GEO: GSE227265 for scATAC-seq; GSE125449 and GSE151530 for scRNA-seq)[7,11,17]. This matched multiomics design ensured coordinated analysis of chromatin accessibility and transcriptional states within the same tumors, providing sufficient statistical power (n = 8) to resolve cell type-specific regulatory networks.

Quality control: For each patient tumor sample, a list of unique ATAC-seq fragments with associated barcodes was imported into the R package ArchR for quality control and doublet removal. To ensure that each cell was both adequately sequenced and had a high signal-to-background ratio, we filtered cells with less than 1000 unique fragments and enrichment at TSS below 4 in ArchR’s createArrowFiles function. Cells with doublet enrichment scores (ESs) > 1 were removed using the ArchR’s filterDoublets function.

Reduction and clustering: ScATAC-seq clustering was performed by adapting the approach of previous publications. Firstly, we used the R package ArchR to construct a matrix for each cell by feature in each patient. Then we added a column with the cell barcode id to the overlaps object and creating a sparse matrix. To reduce the dimensionality of these genomic tile features, we applied the iterative latent semantic indexing (LSI) procedure implemented in ArchR. In brief, this procedure first performs term frequency-inverse document frequency normalization to upweight more informative features, followed by an initial LSI reduction on the most accessible tiles.

Batch correction: Batch effect correction was implemented through Harmony integration. The algorithm refined Iterative LSI embeddings through mutual nearest neighbor alignment, generating harmonized low-dimensional representations that retained biological variability while suppressing technical noise for downstream analysis. FindClusters is then used to identify low-resolution clusters, within which feature counts are summed across all cells in each cluster to identify the top 14000 most variable features at resolution 0.4. This process is repeated 4 times by taking the top 14000 most variable tiles from the first iteration and using them as the most accessible tiles for the second iteration. The iterative LSI procedure produced 50 LSI dimensions, which were then collapsed into a two-dimensional uniform manifold approximation and projection (UMAP) embedding using ArchR’s addUMAP.

Peak calling: For each cluster, peak calling was performed on Tn5-corrected insertions using the MACS2 callpeak command with the following parameters: “-shift-75-extsize 150-nomodel-call-summits-nolambda-keep-dup all -q 0.05”. The top 20000 extended peak summits for each cluster were quantile normalized. Then these cluster peak sets were merged, and the iteration procedure was applied again. This process is repeated 4 times by taking the top 7000 most variable peaks from the first iteration and using them as the most accessible tiles for the second iteration at resolution 0.8. The iterative LSI procedure produced 50 LSI dimensions, which were then collapsed into a two-dimensional UMAP embedding using ArchR’s addUMAP. These reduced dimensions were input into a Seurat object, and final clusters were identified using Seurat’s FindClusters with a resolution of 0.8 and resulting 12 initial scATAC-seq clusters. Gene signature enrichment and visualization was calculated using Seurat’s AddModuleScore.

Integration of scRNA-seq and scATAC-seq identifyd cell type

The multimodal integration of scRNA-seq and scATAC-seq data employed ArchR’s addGeneScoreMatrix algorithm to compute chromatin-derived gene activity metrics, leveraging fragment-to-gene distance weighting (exponential decay model prioritizing proximal chromatin fragments) with strict genomic boundary definitions to prevent cross-gene interference. To transfer cell type labels from scRNA-seq, Seurat’s canonical correlation analysis implementation was used in the ArchR (FindTransferAnchors, TransferData), which assigned a cell type subcluster identity to each scATAC-seq cell from the matching scRNA-seq data, along with an associated label prediction score. This label transfer was limited to cells from the same patient dataset. Only those with a label prediction score > 0.5 were included in downstream analyses. Additionally, only inferred cell type subclusters with > 10 cells were considered for downstream analysis to ensure sufficient cells for peak calling in each cluster. To biologically validate cluster, infer CNV analysis and functional enrichment analysis were further applied.

Post-annotation pseudobulk profile reconstruction involved MACS2-based chromatin peak re-identification (parameters: -shift-75-extsize 150-nomodel-call-summits-nolambda-keep-dup all -q 0.05) followed by ArchR’s iterative merging protocol to unify peak matrices across cellular cohorts. Differential accessibility landscapes were resolved via ArchR’s getMarkerFeatures implementing Wilcoxon rank-sum testing with TSS enrichment and fragment count covariate adjustments, defining regulatory elements under stringent thresholds (| log₂[fold change] |≥ 4, Benjamini-Hochberg adjusted P < 0.001).

scRNA-seq CNV analysis

The R package inferCNV was used to infer putative copy number events for each cluster. Since all tissues involved in this study were tumor derived from tumors, we chose immune cells to mimic the normal background. The computational analysis was performed using parameter settings: Denoising = True, CNV-cutoff = 0.1, scale-normalization = off, and hidden Markov model disabled. For visual comparison of CNV patterns between sample groups, comparative boxplot visualization was constructed using ggplot2 to display the distribution of predicted CNV events across different biological clusters.

Gene ontology/pathway enrichment of cluster signal analysis

To characterize hepatic cellular heterogeneity at single-cell resolution, functional enrichment profiling of cluster-specific marker genes was implemented according to established methodology[18]. The cluster Profiler package was employed for overrepresentation analysis with Kyoto Encyclopedia of Genes and Genomes pathway annotations as the reference database. Statistical confidence was determined through Benjamini-Hochberg multiple testing correction at false discovery rate < 0.05.

Peak annotation

Peaks were annotated using ChIPseeker. Granges objects were annotated to transcription start sites (TSSs ± 1500 bp) and functional regions using the "https://TxDb.Hsapiens.UCSC.hg38.knownGene" database. Gene symbols were linked via the "https://org.Hs.eg.db" database.

Decoding cellular composition and cell-type-specific gene correlation networks in HCC bulk RNA seq

To investigate cellular composition and assess marker gene correlations in HCC bulk RNA-seq datasets, we performed cell-type deconvolution using the co-expression-based decomposition algorithm implemented by Bisque. For the TCGA-LIHC and LCI cohorts, normalized RNA-seq (TCGA) and microarray (LCI) expression matrices were processed with uniform analytical parameters and consistent marker gene sets derived from our cluster-specific signature definitions. This cross-platform standardization enabled comparative assessment of cell-type proportion estimates between the two cohorts. Validation of decomposition reliability was achieved through systematic evaluation of marker gene co-expression patterns across all cellular subtypes.

Comparative analysis of HCC vs adjacent non-tumor tissues was conducted through paired Wilcoxon signed-rank testing to identify tumor-associated alterations in cellular prevalence. Cell-type enrichment scoring preceded statistical testing to ensure biological relevance of interrogated markers. Multiple hypothesis correction accounted for simultaneous evaluation of six distinct cell populations, employing false discovery rate (FDR) adjustment at 0.05. Our combined analysis incorporated 268 matched tumor-normal pairs across both cohorts (TCGA: n = 59; LCI: n = 209), providing sufficient statistical power to detect subtle changes in cellular distribution.

This integrative approach combined molecular deconvolution with robust statistical frameworks to characterize TME alterations. Importantly, the implementation of platform-specific normalization protocols while maintaining constant marker gene sets ensured methodological consistency across distinct expression profiling technologies. The validation through co-expression pattern analysis further strengthened confidence in the biological plausibility of the estimated cellular proportions.

TF motif enrichment analysis in the scATAC-seq data

We assessed TF activities using two computational methods in ArchR’s ChromVAR deviation scores and TF foot printing. For motif enrichment, Hypergeometric Optimization of Motif EnRichment motif annotations were added via ArchR’saddMotifAnnotations. Combing 8771 cluster differential peaks (Log fold change > 4 and FDR < 0.05), ArchR’speakAnnoEnrichment were further used to calculate celltype related motifs. To compute GC bias-corrected motif deviation scores, ArchR’s addDeviationsMatrix and getVarDeviations were used to compute GC bias-corrected motif deviation scores. Top 5 TFs were involved, with FDR < 0.001.

For TF footprinting analysis, we first retrieved the positions of Hypergeometric Optimization of Motif EnRichment motifs using the ArchR’s getPositions. Next, we computed the footprints for selected motifs with the ArchR’s getFootprints function. To account for Tn5 bias and for visualization purposes, ArchR’s plotFootprints was employed under the parameter “normMethod = Divide” for normalization. To visualize the sequence logo for TFs, firstly, position weight matrix were obtained through TFBSTools and ArchR’sgetPeakAnnotation. Then, TFBSTools’s bg transferred position weight matrix to ProbMatrix, which could be identified by ggseqlogo in R and drawed TF sequence logo.

Prediction of gene - peak correlations and peak coaccessibility

First, ArchR’s addGeneScoreMatrix was utilized to calculate gene activity scores, which reflect the transcriptional activity of genes based on chromatin accessibility within gene and regulatory regions. Next, ArchR’saddPeak2GeneLinks was used to add peak-to-gene (P2G) links with the parameters “corCutOff = 0.75, and addEmpiricalPval = TRUE” to calculate the correlation between chromatin accessibility at specific peaks and the gene activity scores of genes. Furthermore, ArchR’s getPeak2GeneLinks was used to obtain P2G links with the parameters “corCutOff = 0.25, FDRCutOff = 0.001, varCutOffATAC = 0.5, and varCutOffRNA = 0.5”. We obtained links between 49720 peaks and 8007 genes.

To calculate peak coaccessibility, we utilized ArchR’s addCoAccessibility and getCoAccessibility to add coaccessibility scores to peaks with default parameters “dimsToUse = 1:30, corCutOff = 0.8, and varCutOffRNA = 0.5”. Results showed all significant associations survived rigorous FDR filtering thresholds (< 0.001), precluding the need for additional statistical adjustments. Coaccessibility analysis utilized ArchR’s default framework for identifying spatially coordinated chromatin regions, with a stringent correlation threshold (r > 0.8) informed by validated practices in single-cell epigenomics. Sensitivity testing at lower thresholds (r > 0.5-0.7) increased network complexity but introduced unvalidated linkages, aligning with recommendations from ArchR’s developers[19].

Inference and analysis of cell-cell communication

Intercellular signaling networks were systematically investigated by applying CellChat to single-cell transcriptomic datasets. This framework integrates mass action kinetic principles with differential ligand-receptor expression profiling and permutation-based hypothesis testing to reconstruct context-dependent cell signaling events. The algorithm prioritizes biologically relevant interactions by identifying overexpressed ligands and receptors through pairwise comparisons across cellular subpopulations. Signaling flux between distinct cell clusters was mathematically modeled by assigning probabilistic weights to ligand-receptor pairs, reflecting interaction likelihoods derived from expression correlation and stoichiometric balance. Through this integrative computational paradigm, CellChat delineates potential oncogenic signaling axes mediating crosstalk between malignant and stromal compartments, thereby elucidating tumor microenvironmental communication architectures. All analytical workflows incorporated default quality control parameters and bootstrap resampling procedures to ensure statistical rigor.

Gene ESs calculated by gene set variation analysis analysis

The ESs for the gene set were calculated using the gene set variation analysis package (v1.48.3) in R, with the following rationalized parameters: MaxDiff = TRUE, and kernel-based estimation. The output ESs represent the relative activity of the chromatin accessibility hub genes in each sample. Higher ES values reflect stronger coordinated expression of the hub gene set, suggesting chromatin accessibility-driven transcriptional synergy. Statistical associations between the ESs and clinical outcomes (e.g., survival, diagnosis) were analyzed using standard regression frameworks.

Diagnostic performance evaluation using ROC analysis

The receiver operating characteristic (ROC) curve’s integrated value area under the curve (AUC) was computed to assess the differentiation capacity of central gene expression profiles in matched tumor-normal tissue specimens. To ensure scientific rigor, the 95% confidence interval (CI) of the AUC was estimated using 5000 bootstrap re-samplings, a nonparametric method robust to sample size variability. Statistical significance against random prediction (AUC = 0.5) was computed using the DeLong test, which evaluates differences in AUC values via asymptotic covariance matrices. For subgroup analyses within a single cohort [e.g., low alpha fetoprotein (AFP] and high-AFP subgroups in the LCI cohort), P values were adjusted for multiple comparisons using the Benjamini-Hochberg procedure. Cross-cohort validations (e.g., TCGA and GSE25907) were considered independent replication tests and thus not subject to multiple testing correction. Sensitivity and specificity at the optimal cutoff were determined using the Youden index (J = sensitivity + specificity - 1) to balance clinical applicability. All analyses were implemented in R with the pROC package to ensure reproducibility and transparency.

Prognostic significance assessment via survival analysis

We assessed the prognostic relevance of hub gene expression using time-to-event survival analysis. Kaplan-Meier curves were generated to visualize differences in overall survival between patients stratified into high- and low-expression groups. The optimal expression cutoff for stratification was objectively determined using the survminer package with the maxstat algorithm, which maximizes the log-rank statistic without prior assumptions. Survival disparities between groups were quantified by the log-rank test, a nonparametric method robust to censored data. Additionally, univariate Cox proportional hazards regression models were employed to compute hazard ratios (HRs) with 95%CIs, reflecting the mortality risk associated with gene overexpression. Proportional hazards assumptions were verified using Schoenfeld residual-based tests (global P > 0.05) to ensure model validity.

Cell culture, drug treatment and siRNA transfection

Cell culture: SNU449 cells were obtained from the American Type Culture Collection. SNU449 was culturedin RPMI-1640 medium supplemented with 10% fetal bovine serum (FBS) and 1% penicillin/streptomycin and incubate at 37 °C with 5% CO₂.

Drug treatment: SNU449 cells was treated with 1 μM A485 (p300/CREB-binding protein (CBP) inhibitor) in DMSO for 48 hours[20]. Histone 3 lysine 27 acetylation (H3K27Ac) validation was performed through western blotting with H3K27Ac antibody (PTMBIO # PTM-116RM, 1:1000) and β-actin (Cell Signaling Technology #4967S, 1:5000) for normalization. Gene expression of midkine (MDK) (1:1000; Proteintech #11009-1-AP) and nucleolin (NCL) (1:1000; Proteintech #83380-1-RR) was further identified in biological triplicates with statistical significance assessed by student’s t-test (P < 0.05).

siRNA transfection: siRNA was designed using the Dharmacon™ siRNA Design Tool. SNU449 cells were transfected using jetPRIME Transfection Reagent (Polyplus). Briefly, cells were seeded in 6-well plates and cultured overnight to reach 60–70% confluency. siRNA complexes were prepared by mixing 25 nM DGAT1-targeting siRNA (sense: 5′-UGG ACUACUCACGCAUCAUTT-3′) or scramble control siRNA (5′-UUCUCCGAACGUGUCACGUTT-3′) with 4 μL jetPRIME reagent in 200 μL jetPRIME Buffer, followed by 15-minute incubation at room temperature. The complexes were added dropwise to cells and incubated for 24 hours without medium replacement. Knockdown efficiency was validated by Western blotting using anti-DGAT1 (1:1000; Proteintech #82945-1-RR) and β-actin (1:5000; loading control). All experiments were performed in biological triplicates with statistical significance assessed by student’s t-test (P < 0.05).

Cell proliferation assays

EdU assay: Cell proliferation was evaluated using the 5-ethynyl-2’-deoxyuridine (EdU) Cell Proliferation Kit (Elabscience, #E-CK-A377). Briefly, SNU449 cells transfected with DRG1-targeting siRNA or scramble control were cultured for 24 hours, followed by incubation with 10 μM EdU for 2 hours at 37 °C. Cells were fixed with 4% paraformaldehyde, permeabilized with 0.5% Triton X-100, and stained with Alexa Fluor 594-conjugated EdU detection reagent as per the manufacturer’s instructions. Nuclei were counterstained with 4’, 6-diamidino-2-phenylindole, dihydrochloride (DAPI) (1:5000, Thermo Fisher Scientific, #D1306). Images were captured using a fluorescence microscope (Zeiss Axio Observer) and analyzed with ImageJ software. The proportion of EdU-positive cells was calculated from ≥ 3 independent biological replicates and compared using Student’s t-test (P < 0.05).

Cell counting kit-8 assay: Cell viability was measured using the cell counting kit-8 (Biosharp, # BS350A). At 24, 48, 72 ,96, 108 hours post-transfection, 10 μL of cell counting kit-8 reagent was added to each well (seeded at 3 × 10³ cells/well in 96-well plates). Following 2 hours of incubation at 37 °C, absorbance at 450 nm was recorded using a microplate reader (Tecan). Background absorbance (medium only wells) was subtracted, and viability was normalized to the scramble control. Results are shown as mean ± SD of triplicate technical replicates across three independent experiments. Significance was calculated using a Student’s t-test (P < 0.05).

RESULTS

Chromatin accessibility facilitates tumorigeneses-related pathways

We obtained HCC tissues from 17 patients through TCGA database (Supplementary Table 1). ATAC-seq signal accumulated within ± 2 kb of the TSSs, which was consistent with canonical chromatin accessibility patterns (Figure 1A). Annotation of HCC-specific chromatin accessibility peaks (defined by Corces[3]) revealed their predominant localization in promoter (36.03%), intronic (35.32%), and distal regulatory regions (21.21%) (Figure 1B). Compared with the other genomic elements, the promoter-localized peaks presented significantly higher chromatin accessibility scores (P < 0.001, Figure 1C), indicating preferential regulatory activity.

Figure 1 Bulk the assay for transposase-accessible chromatin with sequencing dataset identifies chromatin accessibility dynamics and tumor-associated pathways in hepatocellular carcinoma. A: Transcription start site enrichment profiles (± 5 kb from transcription start site) validating sequencing quality across hepatocellular carcinoma (HCC) and paired normal samples; B: Genomic distribution (proportions of promoter, intronic, and distal regions) of chromatin accessibility peaks unique to HCC; C: Chromatin accessibility scores across genomic compartments. Promoter regions show the highest accessibility (Wilcoxon test, P < 0.001); D: Top enriched REACTOM pathways for genes with recurrent promoter-accessible peaks (adjusted P < 0.05); E: Weak but significant correlation between promoter accessibility and corresponding gene expression (Pearson r = 0.24, P < 0.001); F: Immune infiltration-associated chromatin features across HCC-specific peaks. Criteria for immune-associated peaks: Enhanced accessibility in immune cells; correlation with cytolytic activity. TSS: Transcription start site; UTR: Untranslated region; PE: Phosphatidyl ethanolamine; DC: Dendritic cells; NK: Natural killer cells; HSC: Hematopoietic stem cells.

Therefore, we calculated the frequency of genes in which the promoter region exhibited an HCC-specific peak. Through negative binomial modeling, we identified 65 genes with recurrent chromatin accessibility alterations. Pathway analysis of these genes revealed multiple liver- and tumor- related pathways, including nuclear receptor subfamily 1 group H member 2 (NR1H2) and nuclear receptor subfamily 1 group H member 3 (NR1H3) mediated signaling, Bile acid and bile salt metabolism and so on (q < 0.001, Figure 1D). This multipathway convergence suggested coordinated chromatin remodeling during hepatocarcinogenesis.

We further evaluated the promoter accessibility-expression correlations. Surprisingly, promoter accessibility was only moderately correlated with transcriptional outputs (Figure 1E, Pearson r = 0.24, P < 0.001). Given that immune-infiltrating cells contribute > 50% of the variance in chromatin accessibility according to Corces et al’s criteria[3] (enhanced accessibility in immune cells; Pearson correlation with cytolytic activity), we subclassed these peaks by immune-association status. Notably, immune-related peaks accounted for 52.7% of the total peaks (Figure 1F), underscoring the confounding effects of TME heterogeneity in bulk tissue analyses. Overcoming this technical constraint, we utilized single-cell assay for transposase-accessible chromatin sequencing to delineate lineage-specific chromatin organization-transcriptional interplay mechanisms.

Single-cell multiomics integration resolves TME heterogeneity

We performed integrative analysis of matched scRNA-seq and scATAC-seq datasets from eight treatment-naïve HCC specimens (Figure 2A)[7,11,17]. The cohort included patients aged 42-81 years (median age 63.5 years) with diverse clinical characteristics: Histological grades (well-differentiated to undifferentiated), etiologies (HBV, hepatitis C virus, non-alcoholic steatohepatitis, and autoimmune hepatitis), and tumor stages (Supplementary Table 2).

Figure 2 Integrated single-cell multi-omics characterization of hepatocellular carcinoma tumor from the same patients. A: Uniform manifold approximation and projection plot of 15660 single-cell RNA-sequencing (scRNA-seq) cells (left) and UMAP plot of 10201 single-cell assay for transposase-accessible chromatin with sequencing (scATAC-seq) cells (right) color-coded by patient origin; B: UMAP plot of 15660 scRNA-seq cells (left) and UMAP plot of 10201 scATAC-seq cells (right) color-coded by annotated cell types; C: Dot plot showing canonical marker gene expression (rows) across cell types (columns) in scRNA-seq. Dot size: Fraction of expressing cells; color: Mean normalized expression; D: Dot plot showing chromatin accessibility of marker gene promoters (rows) across cell types (columns) scATAC-seq. Dot size: Fraction of expressing cells; color: Mean normalized expression; E: Transcription start site enrichment profiles stratified by major cell types. Malignant hepatocytes exhibit reduced transcription start site signal compared to non-malignant cells; F: Fragment size distributions across cell types. Malignant cells display longer nucleosome-free fragments. UMAP: Uniform manifold approximation and projection; scRNA-seq: Single-cell RNA-sequencing; ATAC-seq: Assay for transposase-accessible chromatin using sequencing; scATAC-seq: Single-cell assay for transposase-accessible chromatin with sequencing; TSS: Transcription start site.

Following rigorous quality control and normalization procedures (including mitochondrial gene filtering, doublet removal, and library size normalization; Methods section), we performed unsupervised clustering on the integrated scRNA-seq dataset. Cell type annotation was achieved through canonical marker gene expression profiling, revealing six major cell types (Figure 2B).

The cell type labels were transferred from the scRNA-seq dataset to the scATAC-seq dataset via the ArchR (v1.0.2), with a predicted score over 0.5 (Figure 2B, Supplementary Figure 1A). Notable interpatient variability in cellular subpopulation proportions was observed (Supplementary Figure 1B and C), reflecting the intrinsic biological diversity of HCC. To validate clustering robustness, we generated marker gene expression dot plots (Figure 2C and D). Malignant cell identification was further validated via CNV analysis. Both the scRNA-seq and scATAC-seq data consistently revealed elevated CNV scores in the malignant clusters than in the nontumor populations (P < 0.001, Supplementary Figure 1D-F), confirming the robustness of classification.

Distinct chromatin accessibility patterns revealing transcriptional plasticity in tumor cells

Systematic comparison of chromatin accessibility features across annotated cell types revealed malignancy-specific epigenetic remodeling. After controlling for intra-patient variability using Harmony-integrated LSI dimensions, TSS enrichment analysis revealed significantly attenuated signals in malignant hepatocytes compared with nonneoplastic populations from the same patients (P < 0.001, Figure 2E, Supplementary Figure 2A). This pattern was particularly pronounced in undifferentiated tumors compared with their differentiated counterparts (a 17% reduction, P < 0.001; Supplementary Figure 2B). Since such a reduction could be contributed to the overall active chromatin regions in tumor cells[3], we further verified the peak average scores per cell, which revealed that both endothelial cells and malignant cells presented relatively high peak scores (Supplementary Figure 2C), while TSS enrichment was relatively high in endothelial cells (1.5-fold, P < 0.001, Supplementary Figure 2A).

Compared with nonmalignant populations, malignant cells presented both increased fragment size (164.9 vs 155.8, P < 0.001, Supplementary Figure 2D) and quantity (1.5-fold, P < 0.001; Supplementary Figure 2E) compared to non-malignant populations. However, analysis of nucleosome-free regions revealed a paradoxical pattern (Figure 2F). We further validated this finding in undifferentiated tumors, which presented 1.4-fold more fragments than well-differentiated samples did (P < 0.001, Supplementary Figure 2F).

These findings suggested that in malignant cells, the quantity of accessible regions was chosen instead of the quality of accessible regions. The increased number of nucleosome-free fragments, particularly in undifferentiated tumors, indicates that developmental plasticity enables microenvironmental adaptation[4,21].

Cell type-specific chromatin architecture reveals conserved malignant signatures

Comparative chromatin accessibility analysis identified 8780 cell type-specific accessible regions (FDR < 0.001, log₂FC > 4; Figure 3A). Peaks in malignant cells (24.11%) showed preferential promoter accessibility (Supplementary Figure 3A), while the accessibility score did not rank highest (Supplementary Figure 3B).

Figure 3 Cell type-specific chromatin accessibility underlies functional diversity in hepatocellular carcinoma. A: Heatmap of 8780 cell type-specific accessible regions (columns) across cell types (rows). Z-score indicates accessibility deviation from mean; B: Pathway analysis of the single-cell RNA-sequencing marker for each cell type using the Kyoto Encyclopedia of Genes and Genomes pathway database. The score indicated normalized enrichment scores and q value denote Benjamini-Hochberg-adjusted P values; C: Pathway analysis of the single-cell assay for transposase-accessible chromatin with sequencing (scATAC-seq) marker peak enriched genes for each cell type using the Kyoto Encyclopedia of Genes and Genomes pathway database. The score indicated normalized enrichment scores and q value denote Benjamini-Hochberg-adjusted P values; D: Expression of scATAC-seq marker peak enriched genes in Liver Cancer Institute cohort grouped by cell peak. The upper panel shows the T-statistic from a paired t-test between tumor and adjacent non-tumor samples; E: Expression of scATAC-seq marker peak enriched genes in The Cancer Genome Atlas cohort grouped by cell peak. The upper panel shows the T-statistic from a paired t-test between tumor and adjacent non-tumor samples. FDR: False discovery rate; FC: Fold change; scRNA-seq: Single-cell RNA sequencing; ATAC-seq: Assay for transposase-accessible chromatin using sequencing; scATAC-seq: Single-cell assay for transposase-accessible chromatin with sequencing; KEGG: Kyoto Encyclopedia of Genes and Genomes; LCI: Liver Cancer Institute; TCGA: The Cancer Genome Atlas.

Marker genes identified by scRNA-seq and promoter-localized peaks (n = 5622) were prioritized to identify proximal transcriptional regulators (Figure 3B and C). Gene ontology analysis of the top 50 cell type-enriched genes via two single-cell omics analyses revealed differences across cell types. The signals in fibroblasts were related to both the organization and function of the cytoskeleton. The characteristics of the endothelial cell signals were accordingly enriched in focal adhesion regions. Genes in B cells and T cells from the two omics analyses were consistently enriched in cell type-related receptor signaling pathways. With respect to myeloid cells, in addition to their role in regulating the immune response, other functions differed from those of the scATAC-seq and scRNA-seq data. Modality-specific pathway divergence in malignant cells emphasized the transcriptional regulatory network complexity and cellular heterogeneity.

Deconvolution of the bulk RNA-seq data (LCI and TCGA cohorts; via Bisque) revealed conserved expansion of malignant hepatocytes in tumor vs adjacent tissues (P < 0.001, Figure 3D) and the TCGA cohort (P < 0.001, Figure 3E). We additionally conducted a parallel analysis by selecting all tumor-upregulated genes with a log fold change greater than 1 from the bulk RNA-seq data and evaluated their average expression in the single-cell datasets. Consistently, we observed that malignant cells presented the highest bulk scores (Supplementary Figure 4A and B). In summary, the signatures enriched in the scATAC-seq datasets presented malignant chromatin-transcriptional coordination in the bulk sequencing data.

Cell type-specific transcription regulators determine cell function

The analysis of cell type-specific chromatin accessibility peaks revealed lineage-defining TF landscapes (Supplementary Figure 5). Surprisingly, HNF1A, a master regulator of hepatocyte differentiation, was enriched in malignant cells, suggesting partial retention of hepatic functions for survival advantage[22]. In addition, the SRY-box transcription factors (SOX) family (SOX9, SOX11, SOX12, and SOX13) was expressed in endothelial cells, implicating mesenchymal genes and steering endothelial cells toward a mesenchymal fate[23]. POU factors (POU2F2, POU5F1B and POU2F3), which regulate antibody-secreting plasma cell differentiation[24], were identified in B cells. Spi-1 proto-oncogene and spi-2 proto-oncogene were enriched in myeloid cells, indicating their role in immune modulation, associated with M2-like polarization and immune tolerance[25]. These findings delineate how cell-specific TFs orchestrate cellular functions within the HCC TME, from immune modulation to malignant adaptation.

Malignant hepatocytes exhibit extensive network rewiring

Extending beyond promoter analysis, we integrated the scATAC-seq and scRNA-seq data via the Peak2 GeneLinkage module to map distal regulatory elements. In total, 33,769 P2G links were identified, with 75.6% colocalizing with cell type-specific accessible regions (Figure 4A and B). Most of these peaks were in promoters (44%), introns (27%) and distal elements (21%, Supplementary Figure 6A). Malignant cells presented the highest P2G connectivity (Figure 4C), with no significant differences in genomic distribution (Supplementary Figure 6B).

Figure 4 Malignant cells harbored multiple enhancer-like peaks driving oncogenic transcription. A: Row scaled heatmaps of statistically significant peak-to-gene (P2G) links. Each row represents expression of a gene (left) correlated to chromatin accessibility of a peak (right); B: Overlap between P2G links and cell type-specific accessible regions; C: Stacked bar plot showing P2G uniquely related to specific marker across cell type. Malignant cells exhibit the highest regulatory complexity; D: Stacked bar plot showing number of genes per peak linked to across cell type. Malignant peaks target more genes (P < 0.001); E: Correlation coefficients (Pearson r) of P2G links across cell types. Malignant cells display stronger correlations (P < 0.001); F: Gene ontology enrichment analysis of the P2G enriched genes for each cell type. The score indicated normalized enrichment scores and q value denote benjamini-hochberg-adjusted P values; G: Browser track showing the accessibility profile at the strongest peak locus across all cell types (chr8 144286909-144366910). Peak was further highlighted by light red shallows; H: Uniform manifold approximation and projection plots of 15660 single-cell RNA-sequencing cells color-coded expression of the strongest peak hub. scRNA-seq: Single-cell RNA-sequencing; ATAC-seq: Assay for transposase-accessible chromatin using sequencing; scATAC-seq: Single-cell assay for transposase-accessible chromatin with sequencing; P2G: Peak-to-gene; UMAP: Uniform manifold approximation and projection;

Malignant hepatocytes displayed increased regulatory network complexity. In addition to the quantity of P2G links, increasing connections were also identified in malignant cells, with an average of 1.97 Linked genes per peak while for non-malignant cells, the number was 1.32 (P < 0.001, Figure 4D). Besides, the expansion of regulatory hubs showed that 54 peaks, predominantly located in promoters (P < 0.001, Supplementary Figure 6C), were linked to more than 8 genes, which served as the maximum number in non-tumor cells. Furthermore, stronger correlation coefficients were observed (P < 0.001, Figure 4E).

Therefore, to further understand the biological significance of these observed results, we conducted functional annotation of cell type-specific enhancer networks. Cell type-specific peak-linked genes were identified, and we used possion distribution modeling to screen them by cell type group (P < 0.05). The top 50 genes were used for enrichment analysis (Figure 4F). Genes in immune cells had similar intrinsic characteristics, such as antigen processing and activation. Fibroblast genes showed effects on hepatocyte differentiation, and genes enriched in endothelial cells affected angiogenesis. Malignant cells were enriched in the electron transport chain, suggesting increased metabolic activity and energy demands. These findings highlight the interplay between cell type-specific chromatin accessibility hubs and their functional roles in shaping the HCC tumor ecosystem.

A super enhancer-like chromatin accessibility peak modulated genes predicts HCC diagnosis and prognosis

A malignancy-specific super enhancer-like peak regulating 18 genes was identified (Figure 4G, Supplementary Table 3). These genes were enriched in malignant cells according to the scRNA-seq data (Figure 4H) and were consistently upregulated in HCC vs non-tumor adjacent tissues across cohorts. Due to the limits of the micro assay, although only 9 genes located in the LCI cohort, they all upregulated in HCC samples (Figure 5A). Additionally, 18 genes were all upregulated in HCC (Figure 5B). To further assess the diagnostic potential of chromatin accessibility-regulated genes, we conducted ROC analysis. In the LCI cohort, this evaluation demonstrated strong discriminative capacity with an AUC of 0.853 (95%CI: 0.813-0.889; Figure 5C). Subsequent survival analysis revealed elevated expression levels of these core genes correlated with unfavorable survival outcomes, showing statistical significance in the LCI cohort (log rank P = 0.02; HR = 1.64, 95%CI: 0.993-2.7; Figure 5D).

Figure 5 Enhancer-like chromatin peaks exhibited diagnostic and prognostic potential in hepatocellular carcinoma. A: Overexpression of enhancer-hub genes in tumors vs adjacent normal tissues in Liver Cancer Institute cohort (n = 209). The upper panels showed the T-statistic from a paired t-test between tumor and non-tumor adjacent non-tumor tissues (NATs); B: Overexpression of enhancer-hub genes in tumors vs adjacent normal tissues in The Cancer Genome Atlascohort (n = 59). The upper panels showed the T-statistic from a paired t-test between tumor and NATs; C: Receiver operating characteristic (ROC) curve analysis of hub genes for diagnosing hepatocellular carcinoma in the LCI cohort (n = 209, area under the curve (AUC) = 0.853, 95% confidence interval (CI): 0.813-0.889); D: Kaplan-Meier survival analysis of hub genes in the LCI cohort (n = 209). High hub gene expression correlates with poor overall survival (P = 0.0233, hazard ratio = 1.64, 95%CI: 0.993-2.70); E: ROC curve analysis of hub genes for diagnosing hepatocellular carcinoma in the LCI cohort in low alpha fetoprotein subgroup (AFP < 300 ng/mL, n = 59, AUC = 0.857, 95%CI: 0.813-0.89); F: ROC curve analysis of hub genes for diagnosing hepatocellular carcinoma in the LCI cohort in high AFP subgroup (AFP ≥ 300 ng/mL, n = 50, AUC = 0.857, 95%CI: 0.813-0.89); G: ROC curve analysis of hub genes for diagnosing hepatocellular carcinoma in the TCGA cohort (n = 59, AUC = 0.923, 95%CI: 0.865-0.968); H: Kaplan-Meier survival analysis of hub genes in the TCGA cohort (n = 59). High hub gene expression correlates with poor overall survival (P < 0.001, hazard ratio = 1.92, 95%CI: 1.02-3.63); I: Gene set variation analysis score of 18 gene in GSE25097 cohort. Expression in hepatocellular carcinoma (n = 268) was higher than in NATs (n = 243), cirrhotic tissues (n = 40), and healthy liver samples (n = 6), P value was adjusted with the Benjamini-Hochberg method; J: ROC curve analysis of hub genes for diagnosing hepatocellular carcinoma in GSE25097 (n = 557, AUC = 0.809, 95%CI: 0.882-0.844). LCI: Liver Cancer Institute; TCGA: The Cancer Genome Atlas; ROC: Receiver operating characteristic; KM: Kaplan-Meier; AFP: Alpha fetoprotein; GSVA: Gene set variation analysis.

To evaluate the differential performance of the 18-gene biomarker signature, the LCI cohort was stratified into AFP-low (AFP < 300 ng/mL) and AFP-high (AFP ≥ 300 ng/mL) subgroups. Remarkably, the biomarker demonstrated robust discriminative ability in both subgroups: AFP-low patients achieved an AUC of 0.86 (95%CI: 0.808-0.915, Figure 5E), while AFP-high patients showed similar diagnostic efficacy (AUC = 0.86, 95%CI: 0.807-0.913, Figure 5F). DeLong’s test comparing the two AUC curves revealed no statistically significant difference (P = 0.87), underscoring the consistent and complementary role of the 18-gene panel alongside established clinical biomarkers such as AFP. These results suggest that the signature provides additive diagnostic value independent of AFP levels.

Consistent with these findings, the TCGA cohort demonstrated even stronger diagnostic and prognostic associations. The AUC for hub genes in the TCGA cohort was 0.923 (95%CI: 0.865-0.968, Figure 5G), and survival analysis confirmed a significant correlation between higher expression of gene set and reduced overall survival (P < 0.001, HR = 1.92, 95%CI: 1.02-3.63; Figure 5H).

To address its specificity for HCC vs other liver pathologies, the 18-gene enhancer signature was validated in the GSE25907cohort (n = 557), including 268 HCC tumors, 243 adjacent non-tumor tissues, 40 cirrhotic tissues, and 6 healthy livers. Consistent with prior cohorts (TCGA and LCI), the biomarker exhibited significantly higher expression in HCC tumors compared to non-tumor, cirrhotic, and healthy tissues (P < 0.001, Figure 5I). ROC analysis demonstrated robust discriminative performance (AUC = 0.81; 95%CI: 0.844-0.882, Figure 5J) across all tissue subtypes, supporting its diagnostic utility as a standalone classifier for HCC even in the presence of heterogeneous liver disease states.

To explore the expression interactions among the hub genes, we constructed co-expression networks in both tumor and non-tumor adjacent tissues with in the LCI cohort and TCGA cohort with Pearson estimates (Supplementary Figure 6C and D). Notably, the hub genes presented significantly stronger co-expression patterns in tumor tissues than in normal tissues, suggesting their potential role in driving tumor-specific regulatory networks (Supplementary Figure 6E and F). These results collectively highlight the diagnostic and prognostic value of chromatin accessibility modulated genes in HCC, as well as their potential involvement in tumor-specific regulatory mechanisms.

To address the biological significance of the identified 18-gene enhancer-like hub, we conducted functional validation through pharmacological and genetic approaches. Pharmacological inhibition of chromatin accessibility using the CBP/P300 inhibitor A485 significantly reduced H3K27Ac levels (P < 0.001, Figure 6A and B), a marker of active enhancers, and concurrently downregulated DGAT1 expression (P < 0.01, Figure 6A and B), directly linking enhancer activity to DGAT1 transcription. siRNA-mediated knockdown of DGAT1 in two HCC cell lines (P < 0.01, Figure 6C and D) resulted in marked reductions in proliferation rates (P < 0.01, Figure 6E and F) and cell viability (P < 0.05, Figure 6G). The consistent outcomes across independent methods and cell lines substantiate the proposed mechanism, underscoring the hub’s functional relevance to tumor growth. This evidence has been integrated into revised figures and discussions, with limitations acknowledged but balanced by the reproducibility and statistical rigor of the current findings.

Figure 6 DGRT1 regulated by chromatin accessibility promoted cell proliferation. A: Western blot analysis of DGAT1 and histone 3 lysine 27 acetylation expression in SNU449 cells treated with 1 μmol/L A485 (CREB-binding protein/P300 inhibitor) for 48 hours. Glyceraldehyde-3-phosphate dehydrogenase served as a loading control; B: Quantification of DGAT1 and H3K27Ac protein levels normalized to GAPDH (n = 3); C: Western blot of DGAT1 and proliferation marker in SNU449 cells transfected with DGAT1-targeting small interfering RNA. Small interfering RNA with non-targeting gene was used as control; D: Quantification of DGAT1 and proliferation marker levels normalized to GAPDH (n = 3); E: Representative images of 5-ethynyl-2’-deoxyuridine assay showing proliferating cells (EdU-positive, red) in DGAT1-knockdown SNU449 cells. Nuclei were stained with 4’, 6-diamidino-2-phenylindole, dihydrochloride (blue); F: Quantitative analysis of EdU-positive cells (n = 3), P value < 0.05; G: Cell counting kit-8 assay demonstrating reduced cell viability in DGAT1-knockdown SNU449 cells over 96 hours (n = 3), P value < 0.05. EdU: 5-ethynyl-2’-deoxyuridine; DAPI: 4’, 6-diamidino-2-phenylindole, dihydrochloride; DGAT1: Diacylglycerol acyltransferase-1; siDGAT1: Small interfering RNA targeting DGAT1; H3K27Ac: Histone 3 lysine 27 acetylation; GAPDH: Glyceraldehyde-3-phosphate dehydrogenase; PCNA: Proliferating cell nuclear antigen.

Chromatin coaccessibility and intercellular signaling shape the tumor-promoting microenvironment

To elucidate the functional interplay between cell type-specific regulatory elements and intercellular communication, we performed chromatin coaccessibility and cell-cell interaction analyses. Coaccessibility is defined as the correlation of chromatin accessibility between two peaks across cells. Peaks that are coaccessible exhibit similar patterns of accessibility in cells, suggesting that they might be functionally related or physically interacting. Using stringent correlation thresholds (r > 0.8, resolution = 1), we identified 26,448 significant peak coaccessibility relationships, including 2550 intercellular interactions (Supplementary Table 4). Malignant cells and endothelial cells harbored the most coaccessible peaks (Figure 7A, Supplementary Table 4), suggesting extensive epigenetic coordination.

Figure 7 Epigenetic coordination between malignant and stromal cells. A: Co-accessible chromatin peaks between tumor and non-tumor cells (log10 transformed counts); B: Ligand-receptor interactions initiated by malignant cells in tumor microenvironments, quantified via CellChat; C: A bar plot showing the log10 of co-accessible peaks number among tumor cells and non-tumor cells across differentiated tumors and nondifferentiated tumors; D: A bar plot showing the number of ligand-receptor interactions when setting malignant cells as source among tumor cells and non-tumor cells across differentiated tumors and nondifferentiated tumors; E: Gene ontology enrichment for genes linked to malignant co-accessible regions. The score indicated normalized enrichment scores and q value denote Benjamini-Hochberg-adjusted P values; F: A circle network plot showing communication networks of the three common signaling pathways, Collagen, midkine (MDK) and peptidylprolyl isomerase A across differentiated tumors and nondifferentiated tumors. Line width indicates interaction strength; G: Expression of MDK pathway ligands (MDK, orange) and receptors (green) across differentiation states. Undifferentiated tumors exhibit upregulated MDK; H: MDK and nucleolin expression was lowest in well-differentiated samples grouped by Hoshida subtypes; I: Western blot of MDK, NCL and H3K27Ac expression in SNU449 cells transfected with 1 μmol/L A485 (CREB-binding protein/P300 inhibitor) for 48 hours. glyceraldehyde-3-phosphate dehydrogenaseserved as a loading control; J: Quantification of MDK, NCL and H3K27Ac levels normalized to glyceraldehyde-3-phosphate dehydrogenase (n = 3). scRNA-seq: Single-cell RNA-sequencing; MK: Midkine related pathway; CypA: Peptidylprolyl isomerase A; MDK: Midkine; SDC: Syndecan; LRP: LDL Receptor Related Protein; NCL: Nucleolin; H3K27Ac: Histone 3 lysine 27 acetylation; CBP: CREB-binding protein; GAPDH: Glyceraldehyde-3-phosphate dehydrogenase.

To determining the functionally-related signaling of these cellular coaccessible regions, we applied Kyoto Encyclopedia of Genes and Genomes analysis and found that the peaks enriched in tumor cells were related to P450 metabolism-related genes, endothelial cells in extracellular traps, monocytes in cholinergic synapse-related genes, and fibroblasts in the cytoskeleton in muscle cells (Supplementary Figure 7A).

Subsequently, we applied CellChat, which uses mass action models to infer cellular communication (Figure 7B). By setting signals in malignant cells as source targets, upon cellular interactions, fibroblasts (n = 38) and endothelial cells (n = 34) were found to be enriched. To better depict the signal changes that occur during tumor progression, upon the historical grading, we divided the samples into a differentiated group and an undifferentiated group. Equivalent coaccessible regions but more cellular communications occured in the undifferentiated group (Figure 7C and D).

Upon co-accessible analysis, we extracted genes that were subject to malignant cells and unique to the remaining cell types. Gene ontology analysis revealed that coaccessible regions occurred among malignant cells and endothelial cells, affecting transcription activity (Figure 7E). For the CellChat results, a total of 123 ligands and receptors were identified upon comparison, of which the MDK pathway was the most enriched in undifferentiated tumors with both incoming and outgoing strengths. In addition, Collagen was the bottom incoming signaling pathway while peptidylprolyl isomerase A was the top pathway (Supplementary Figure 7B, Figure 7F). The pathway enrichment of the ligands and receptors identified MDK and NCL, which were both further increased in the undifferentiated group. Apart from that, the remaining cellular signaling pathways differed between the cell types (Figure 7G). SDC2 was enriched in malignant cells, LRP1 in fibroblasts and SDC1 in malignant cells, but all of these genes were elevated in the undifferentiated group. These findings established a differentiation hierarchy in tumor-stroma communication, where undifferentiated tumors develop enhanced epigenetic coordination and signaling complexity to sustain their aggressive phenotype.

To address etiology-driven heterogeneity, we reclassified tumors using the Hoshida molecular subtypes (S1/S2/S3), revealing that our chromatin-defined hub genes MDK and NCL exhibited significantly lower expression in well-differentiated S3 tumors than in aggressive S1/S2 subtypes (P < 0.01, Wilcoxon test, Figure 7H), aligning with their differentiation-dependent regulatory roles. Furthermore, pharmacological inhibition of H3K27Ac deposition (via A485) reduced chromatin accessibility and suppressed both MDL and NCL expression (P < 0.01, Figure 7I and J). These findings underscore the etiologically agnostic yet differentiation-dependent functional roles of MDK/NCL in maintaining the tumor phenotype, independent of primary disease drivers such as HBV or non-alcoholic steatohepatitis.

DISCUSSION

Earlier single-cell transcriptomic study mapped HCC evolution under therapeutic pressure, revealing transcriptional plasticity as a survival strategy[17]. The dynamic regulation of chromatin accessibility is a cornerstone of cancer progression[3], yet its hierarchical organization and functional consequences in HCC remain incompletely understood. Here, we dissect the chromatin accessibility mechanisms underlying such plasticity, integrating multi-omics profiling to resolve enhancer-driven oncogenesis. This study extends our discovery of tumor-stroma metabolic symbiosis[7] by identifying shared chromatin coaccessibility as its epigenetic substrate. Our study’s sample size (n = 8 paired tumors) reflects current limitations in matched single-cell multiomics data availability for HCC. While larger cohorts are ideal for generalization, the integration of scRNA-seq and scATAC-seq at cellular resolution provides high-dimensional insights into tumor heterogeneity. Importantly, all datasets were derived from the same patient cohort to minimize batch effects, enabling robust tumor-specific regulatory inferences.

Our multi-scale analysis identified chromatin accessibility remodeling as a critical molecular mechanism driving HCC pathogenesis[5,26]. Bulk ATAC-seq revealed 65 genes with enriched accessible peaks, many linked to oncogenic pathways (Figure 1D). At single-cell resolution, malignant cells exhibited chromatin landscapes distinct to stromal and immune populations (Figure 3B and C), highlighting compartmentalized regulatory programs. For example, T-specific peaks were linked to immune modulation pathways (e.g., Th cell differentiation and the T cell receptor signaling pathway), potentially explaining the immunosuppressive phenotype observed in advanced HCC[8]. By integrating bulk and single-cell multiomics sequencing, we dissected chromatin accessibility landscapes in HCC, with single-cell resolution unmasking regulatory complexity obscured by TME heterogeneity (e.g., immune contamination in bulk peaks). Quantitative expansion of chromatin accessible regions, enhancer-like regulatory reprogramming, and TME remodeling. These findings collectively redefine our understanding of epigenetic regulation in liver cancer and suggest novel therapeutic strategies.

Notably, our study establishes enhancer-like chromatin hubs as pivotal contributors to HCC pathogenesis. As for peak identification, our stringent statistical thresholds (log₂FC ≥ 4, FDR < 0.001) prioritized high-confidence cell type-specific chromatin accessibility changes but may have overlooked regions with moderate yet biologically meaningful differential regulation. Future studies integrating dynamic statistical models (e.g., model-based analysis of single-cell transcriptomics) or chromatin velocity frameworks could refine this balance between sensitivity and specificity. Malignant-specific enhancer connectivity (e.g., 33769 peak-gene links; Figure 4A) may mirror cohesin-mediated looping in cancers like breast or ovarian tumors, yet driven by HCC-specific TF networks (e.g., HNF4A, Supplementary Figure 5)[27,28]. Recent studies in breast and gynecologic cancers have established enhancer-driven oncogene activation and intratumoral chromatin heterogeneity as recurring themes in solid tumors. For example, the regulatory elements that lead to the upregulation of clinically relevant oncogenes in breast cancer[28], while chromatin state diversity aligns with ovarian tumor heterogeneity[27]. However, HCC’s enrichment for liver lineage TFs (e.g., HNF4A) and metabolic regulatory elements suggests tissue-contextualized dependencies distinct from hormonally driven cancers[22]. Rigorous cross-tumor comparisons are needed to test whether these features represent HCC-specific mechanisms or variations of broader oncogenic principles.

Our findings align with pan-cancer chromatin accessibility studies[3], which identified immune microenvironment contributions to bulk epigenomic variability. However, by resolving cell type-specific enhancer-like hubs via scATAC-seq, we uncovered HCC-specific dependencies distinct from global chromatin decompaction trends. In malignant cells, 4.3% of enhancer-like peaks regulated ≥ 5 genes (vs 0.47% in non-tumor cells; Figure 4F), forming superhubs that amplify oncogenic transcription. A strongest enhancer-like peak was linked to 18 genes (Figure 5), and it served as biomarker demonstrated robust diagnostic accuracy across three cohorts. Furthermore, the 18-gene enhancer biomarker exhibited AFP-independent diagnostic accuracy (Figure 5F), highlighting its potential to complement AFP in clinical practice, especially for AFP-lower HCC cases. Besides, the 18-gene enhancer signature exhibited significantly higher activity in HCC tumor tissues compared to cirrhotic tissues, adjacent non-tumor tissues, and healthy liver samples (P < 0.001, Kruskal-Wallis test). This confirms its specificity for HCC over non-malignant liver diseases (e.g., cirrhosis) and underscores its diagnostic utility in distinguishing HCC from background liver pathology. However, lack of cholangiocarcinoma cohort validation weakened its specificity to HCC vs other liver diseases. Furthermore, given the scarcity of normal liver tissue samples in current public databases (our study has exhausted available resources from major repositories), the statistical power of this analysis remains limited. We proposed establishing multicenter collaborative cohorts to systematically investigate the impact of varying liver disease stages on HCC molecular features, which will be prioritized in our follow-up studies.

DGAT1 exhibiting the strongest correlation (r = 0.88, P < 0.001). Functional validations targeting this hub, via A485-mediated chromatin suppression and siRNA-based DGAT1 knockdown, demonstrated its critical role in driving proliferation. Although this study delineates chromatin accessibility-mediated regulation of DGAT1 in SNU449 cells, the universality of this mechanism across tissues requires further verification. While our study validated the chromatin accessibility-regulated hub gene DGAT1 in vitro, further investigation is warranted to confirm its oncogenic role and therapeutic potential in vivo. Orthotopic xenograft models, patient-derived organoids, or genetically engineered mouse models would provide critical pathophysiological context. Besides, the remaining 18 genes remain mechanistically unexplored, their upregulation in HCC and survival associations (Figure 5D and E) suggest coordinated oncogenic potential. Enhancer-like peaks were further enriched electron transport chain pathways (Figure 4F), resolving the Warburg effect paradox by epigenetically synchronizing glycolysis and mitochondrial metabolism[29]. These patterns mirror enhancer-driven oncogenesis in breast/ovarian cancers[27,28], yet HCC uniquely retains hepatocyte lineage TFs (HNF4A,FOXA1) and metabolic elements, exposing tissue-contextualized vulnerabilities. Though our study focuses on HCC, future work will explore whether chromatin-driven adaptation is conserved across cancer types using pan-cancer scATAC-seq datasets.

Chromatin coaccessibility and CellChat analyses also revealed differentiation-dependent escalation of tumor-stroma signaling. We further validated this result in Hoshida subtypes (Figure 7H). In undifferentiated tumors, the MDK-NCL axis emerged as a central communication hub - a finding with dual translational implications. Mirroring canonical drivers like Wnt/β-catenin and telomerase reverse transcriptase, MDK-NCL signaling was enriched in aggressive subtypes (Figure 7F-H), stabilizing tumor-stroma ecosystems through epigenetic coordination. And decreased chromatin accessibility of A485 reduced MDK and NCL expression level. While Wnt/β-catenin regulates cell-autonomous proliferation, MDK-NCL may orchestrate microenvironmental reprogramming, suggesting combinatorial targeting strategies[29].

While our study provides valuable insights, several limitations warrant attention. First, while our analysis of cell type-specific chromatin accessibility peaks identified lineage-defining TF motifs, we acknowledge key limitations in equating motif enrichment to functional TF activity[30]. ScATAC-seq captures chromatin accessibility but not post-translational modifications or cofactor interactions that regulate TF activation. motif presence does not guarantee TF binding or transcriptional regulation Future studies integrating proteomic, phospho-specific, or single-cell multiomic datasets will refine these predictions by resolving TF activation states and their dynamic regulation.

Besides, while clustered regularly interspaced short palindromic repeat/clustered regularly interspaced short palindromic repeat-associated nuclease 9 (CRISPR/Cas9) knockout of enhancer-like peaks or overexpression of hub genes would further solidify our conclusions, technical constraints (e.g., complexity of targeting multiple enhancers and time limitations) precluded these experiments. Importantly, the concordant results from A485-mediated chromatin suppression and siRNA-based DGAT1 knockdown functionally connect enhancer activity to DGAT1 expression and proliferation, aligning with our hub network analysis.

Furthermore, while our intra-sample design controls for patient-specific technical variability, we acknowledge that systematic biases inherent to scATAC-seq (e.g., chromatin fragmentation efficiency, TSS capture bias) could influence absolute signal quantification. Experimental validations such as ATAC-qPCR or technical replicates would strengthen these findings but were not feasible within the current study framework. Future work incorporating orthogonal assays will be critical to confirm these observations. Furthermore, the cohort size (n = 8) constrains robust detection of rare cell states or patient subtype-specific regulatory programs. While our coaccessibility networks were constructed using batch-harmonized chromatin embeddings (e.g., Harmony-integrated LSI dimensions), we acknowledge that residual technical variability could still subtly influence network topology. For instance, patient-specific biases in sequencing depth or chromatin fragmentation might propagate into coaccessibility score calculations, despite explicit corrections during dimensionality reduction. Besides, Spatial resolution of chromatin accessibility and transcriptional states, essential for dissecting tumor-niche crosstalk, is absent, highlighting a need for spatially resolved multi-omics in future work. Additionally, we cannot assess how therapeutic interventions reshape chromatin landscapes, as our cohort comprises untreated tumors. Addressing these gaps will clarify how epigenetic heterogeneity influences clinical outcomes and therapeutic vulnerabilities.

In summary, our study delineates the chromatin accessibility of HCC at single-cell resolution, revealing key regulatory hubs and intercellular chromatin cross-talk. We identified the enhancer-like peak, which was validated as a potent diagnostic and prognostic biomarker. This framework redefines the epigenetic vulnerability landscape of HCC, guiding precision strategies to disrupt chromatin-mediated oncogenesis.

CONCLUSION

This study elucidates how chromatin accessibility rewiring drives HCC progression through multiomics integration. Bulk analyses linked promoter accessibility to tumorigenic pathways (e.g., bile acid metabolism), with immune microenvironment heterogeneity explaining limited accessibility-expression correlations. Single-cell resolution revealed malignancy-specific chromatin remodeling: Malignant hepatocytes exhibited attenuated TSS signals but expanded nucleosome-free regions, reflecting transcriptional plasticity. Lineage-specific TFs (e.g., HNF1A in tumors, SOX family in endothelial cells) and complex regulatory hubs (1.97 genes/peak in malignant cells) underscored cell-type specific chromatin accessibility reprogramming. A super enhancer-like peak regulating 18 genes emerged as a robust diagnostic/prognostic biomarker (AUC = 0.81-0.923 across cohorts; survival HR = 1.64-1.92), functionally validated via CBP/P300 inhibition (DGAT1 suppression, P < 0.01) and siRNA knockdown (proliferation impairment, P < 0.05). Undifferentiated tumors showed enhanced epigenetic coordination (e.g., MDK/NCL signaling) and tumor-stroma crosstalk. These findings establish chromatin plasticity as pivotal to HCC pathogenesis and clinical progression.

ACKNOWLEDGEMENTS

We extend our gratitude to all the participants involved in this study.