In this study, 12 HCC patients with their 47 blood samples were included before and after complete RFA therapy in You’An Hospital from October 2017 to December 2019. The diagnosis of HCC was preferably histologically confirmed. In cases when tumor biopsy results were unavailable, diagnosis was established by contrast-enhanced CT or MRI. Treatment in each patient was performed by an interventional radiologist with >5 years of experience, guided by the National Comprehensive Cancer Network (NCCN) and Chinese HCC treatment guidelines. Complete ablation is defined as the presence of an ablative margin of at least 5 mm around the entire tumor, no enhancement in the arterial phase, and no defect in the portal phase on enhanced CT scan. Patients underwent follow-up every 1 month in the first 3 months post-RFA and every 3 months thereafter. The longest follow-up time was 14 months. The definition of early recurrence was the relapse of HCC within 1 year after complete ablation. Patient characteristics are shown in Table 1 and Supplementary Table 1.

The analyses included in this study were performed in accordance with the Declaration of Helsinki and were approved by the Ethics Committee of the Beijing You’An Hospital (Jing You Ke Lun Zi [2022]055# in Chinese). Written informed consent from the patients for the research use of data was obtained before the investigation.

The fresh blood samples were collected in Ethylenediaminetetraacetic acid (EDTA) tubes, and then genomic DNA was extracted and bisulfite-converted (ZYMO Research, Irvine, CA, USA) for the DNA methylation assay. The DNA methylation test panel named HYGEIA was constructed by BioChain (Beijing) Science & Technology Inc. (Beijing, China). It contained 5,117,910 CpG sites on a single-stranded DNA chain, and the detailed information is shown in Supplementary Table 2. The single-chain library was constructed by BioChain and was used for the methylation sequencing of target interval capture.

First, the sequencing quality of the raw data was analyzed using the Fastp software (15) to remove low-quality bases/reads and obtain clean reads. The Bismark alignment software (16) was used to align clean reads with the reference genome (GRCh38). After deduplication, “Bismark_sthylation_detractor” was used to extract CpG sites and count the number of reads. The methylation level of a single CG site = number of reads of the C base/the total number of reads; the methylation level of a gene = the number of reads covering C bases within the gene/the total number of reads within the gene.

Total RNA was extracted using TRIzol (Invitrogen, Life Technologies, Carlsbad, USA) and assessed using Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA) and Qubit Fluorometer (Invitrogen). Total RNA samples that met the following requirements were used in subsequent experiments: RNA integrity number (RIN) > 7.0 and a 28S:18S ratio > 1.8. RNA-sequencing (RNA-Seq) libraries were generated and sequenced by CapitalBio Technology (Beijing, China). The final libraries were quantified using the KAPA Library Quantification kit (KAPA Biosystems, Cape Town, South Africa) and Agilent 2100 Bioanalyzer. After quantitative reverse transcription–polymerase chain reaction (RT-qPCR) validation, libraries were subjected to paired-end sequencing with paired-end 150-bp read length on an Illumina NovaSeq sequencer (Illumina, Inc., San Diego, CA).

The human genome version of hg38 was used as a reference. The sequencing quality was assessed using FastQC (v0.11.5) (17), and then low-quality data were filtered using NGSQC (v2.3.3) (18). The clean reads were then aligned to the reference genome using HISAT2 (v2.1.0) (19) with default parameters. The processed reads from each sample were aligned using HISAT2 against the reference genome. The gene expression analyses were performed using StringTie (v1.3.3b) (20).

Four paired samples (pre-treatment and post-treatment) of non-recurrent patients and three paired samples of recurrent patients were applied to differential expression/methylation analysis.

Differentially expressed genes were identified using the DESeq2 package in RNA-Seq count data. Differentially methylated CpG sites were identified using the Limma package in methylation rate data.

The calculation of the correlation coefficient of the methylation of each CpG site with its corresponding gene expression was performed using Spearman’s correlation in a nested for-loop. The gene regions were categorized into nine groups: promoter (2–3 kb), promoter (1–2 kb), promoter (≤1 kb), 5′UTR, exon, intron, distal intergenic region, 3′UTR, and downstream.

DIABLO mixOmics (version 6.30.0) (21) with sparse partial least-squares discriminant analysis (sPLS-DA) was applied for the integrative analysis of transcriptome with methylome. Gene expression levels were quantified as Transcripts Per Million (TPM) and log2-transformed [log₂(TPM + 1)]. To account for potential inter-sample technical variation, quantile normalization was further performed across all samples prior to integrative analysis. Similarly, quantile normalization was applied across all samples to mitigate batch effects. A total of nine pre-treatment and nine post-treatment samples with transcriptome and methylome datasets were taken for integration. DIABLO was used to identify highly correlated gene signatures across two omics datasets. To identify key immune processes that were differentially activated or suppressed between experimental groups, gene set enrichment analysis (GSEA) was conducted using the gseKEGG functions implemented in the “ClusterProfiler” package. The normalized enrichment score (NES) > 1 and adjusted p-value < 0.05 in GSEA results were defined as significant.

The ClusterGVis package was applied to cluster and visualize the time-series gene expression data of the transcriptome. TPM values were log₂-transformed by dividing the number of reads mapped to a transcript by the transcript length. TPM normalization of RNA-Seq data was used for the ClusterGVis package. The elbow method was used to define five as a suitable cluster value based on the total within sum of squares value. The mfuzz method could be chosen to cluster data and detect genes with a consensus trend. Gene symbol was converted to gene ID using the database “org.hs.eg.db”, and enrichment analysis was performed. The “enrichCluster” function was used to enrich immune-related pathways in the “BP” type with p-value < 0.05.

To evaluate the association between T-cell differentiation signatures and the risk of recurrence after radiofrequency ablation, the following analyses were performed. First, four T-cell differentiation-related pathway gene sets were obtained from the C5 collection of the MSigDB database (http://www.gsea-msigdb.org/gsea/msigdb): regulation of T-cell differentiation in the thymus, regulation of gamma-delta T-cell differentiation, positive regulation of CD8-positive alpha-beta T-cell differentiation, and positive regulation of CD4-positive alpha-beta T-cell differentiation. The corresponding signature scores for each patient’s pre-treatment sample were calculated using the gene set variation analysis (GSVA) algorithm. Then, stratified analysis was performed to compare RFS outcomes between the high-score and low-score groups. For survival analyses, the Kaplan–Meier curves were generated, and log-rank tests were used to derive p-values. Hazard ratios (HRs) and 95% confidence intervals (CIs) were estimated using Firth’s penalized likelihood (to address complete separation in subgroups with no events). To validate the independent prognostic value of T-cell differentiation signatures, univariate and multivariate Cox proportional hazards regression analyses were conducted. The multivariate models incorporated clinical confounders, including age, gender, Barcelona Clinic Liver Cancer (BCLC) stage, and the Child–Pugh score. For external validation, The Cancer Genome Atlas Liver Hepatocellular Carcinoma (TCGA-LIHC) cohort (n = 371) was analyzed using the same T-cell differentiation gene sets, with multivariate adjustment for clinical factors (sex, age, tumor stage, and grade).

After methylation HYGEIA panel and RNA-Seq assays, a total of 18 samples were used in this section (Table 1). Over 2.2 million CpG sites were measured, with the majority located in promoter regions, particularly within 1 kb of transcription start sites, where methylation levels were the lowest (Figures 1A, B). All detected CpG sites were mapped to over 18,000 genes, and most genes contained more than 30 measurable CpG sites. Tens of thousands of CpG sites whose methylation significantly correlated with expression (adj. p < 0.05, |r| > 0.5) were identified, with approximately 75% residing in promoters. Most promoter-proximal CpGs showed negative correlation with gene expression, and this inverse relationship became more pronounced closer to the transcription start site and with stronger correlation coefficients (Figure 1C). At the gene level, average CpG methylation intensity demonstrated a significant moderate negative correlation with expression across all samples (Spearman’s R = −0.36, p-value < 0.001) (Figure 2, Supplementary Table 3). In summary, these results indicated that the correlation between DNA methylation and mRNA expression was complex and not simply a negative correlation.

The sPLS-DA was applied to identify a subset of variables that could explain the variability between nine pre-treatment and nine post-treatment samples. It was demonstrated that the methylome data had a higher variance than the transcriptome data for the first two components (Figure 3A). A global overview of the correlation at the component level is displayed in Figure 3B. It revealed the correlation between the methylation data and transcription data, resulting in a correlation of 0.95 and 0.89 for the first and second components, respectively. All the correlations between pre-treatment and post-treatment samples assessed in multi-omics were shown in a circos plot (Figure 3C) with a cut-off level of 0.8. The result revealed a number of strong negative correlations. Furthermore, unsupervised analysis shown by heatmap clustering (Figure 3D) demonstrated that the highly correlated genes (>0.7) between multi-omics mostly distinguished pre-treatment and post-treatment samples. These gene signatures underwent GSEA to evaluate the median rank in all PBMC samples. The results showed that RFA treatment enhanced immune activity, leading to increased NES of antigen processing and presentation, Th1 and Th2 cell differentiation, and T- and B-cell receptor signaling pathway and reduced NES of PD-L1 expression and PD-1 checkpoint pathway in cancer (Figure 3E). In conclusion, the mixOmics method indicated that the gene signatures with strong correlation between multi-omics play important roles in exploring immune profiles of ablation treatment.

We analyzed RNA-sequencing count data from 14 paired samples (four non-recurrent and three recurrent patients) using the DESeq2 package. With an False Discovery Rate (FDR) threshold of < 0.05, we identified a set of differentially expressed (DE) genes in both patient groups: non-recurrent patients showed more upregulated than downregulated DE genes, while recurrent patients exhibited the opposite trend (Figures 4A, B). For DNA methylation analysis, we used the Limma package to assess methylation rates in the same paired samples. Using the same threshold, we detected numerous differentially methylated (DM) CpG sites in both cohorts, with upregulated DM CpG sites significantly outnumbering downregulated sites in both non-recurrent and recurrent patients (Figures 4C, D). To ensure the robustness of the results, we performed a sensitivity analysis (Supplementary Table 4). Although the number of identified differential features varied across different thresholds, the non-recurrent group had more upregulated than downregulated genes, while the recurrent group was dominated by downregulated genes, remaining consistent across all tested levels.

To analyze how RFA impacts HCC patients with different prognoses, we compared DE and DM genes between the non-recurrent and recurrent cohorts. In non-recurrent patients, RFA treatment suppressed the expression of a subset of genes (characterized by downregulated expression and upregulated methylation) while enhancing the expression of another gene subset (upregulated expression and downregulated methylation) (Figure 5A; Supplementary Table 5). A similar pattern of gene expression regulation was observed in recurrent patients, with RFA leading to both the downregulation and upregulation of distinct gene sets (Figure 5B; Supplementary Table 6).

The Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis of these DE and DM genes revealed the key functional differences between the two groups: antigen processing and presentation, as well as Th1/Th2 cell differentiation signaling pathways, were activated in non-recurrent patients post-RFA and partially in recurrent patients. Additionally, compared to recurrent patients, non-recurrent patients showed inhibition of the PD-L1 expression, PD-1 checkpoint, and TNF signaling pathways (Figure 5C).

We identified the DE genes of five periods, including pre-treatment and post-treatment at 1, 6, 9, and 12 months in HCC patients who underwent complete ablation therapy. We performed differential gene expression analysis between four post-treatment periods and the pre-treatment period. There were 269 upregulated and 319 downregulated genes in the intersection of differential genes in each time period of post-treatment compared with pre-treatment in the non-recurrent patient group. The selenium binding protein 1 (SELENBP1) was the most significantly downregulated in post-treatment samples of non-recurrent patients (Figure 6A; Supplementary Figure 1). Meanwhile, we identified 29 upregulated and 140 downregulated genes in the intersection of differential genes in each time period of post-treatment compared with pre-treatment in the recurrent patient group (Figure 6B; Supplementary Figure 1).

To comprehensively investigate the global temporal patterns of the immune-related impact of RFA, we prepared all genes for cluster analysis in different prognoses using the ClusterGVis package. We clustered those genes into five groups (C1–C5) according to their expression patterns in non-recurrent and recurrent patients. We defined the highly expressed genes in five periods using a line map (Supplementary Figure 2). Using the specific marker genes, we also performed gene ontology (GO) pathway enrichment analysis for each period. Genes specifically expressed in pre-treatment samples in non-recurrent patients (C3) (Figure 7A; Supplementary Table 7) enriched more T-cell immune-related pathways, such as T-cell differentiation and T-helper cell differentiation, than the recurrent group (C1) (Figure 7B; Supplementary Table 8). Ablation treatment induced and enhanced many types of immune responses 6 months in non-recurrent patients (C1 and C5) and recurrent patients (C2, C3, and C4). Nine months after treatment, the adaptive immune system exhibited sustained activation in non-recurrent patients (C2); however, no such immune pathways were enriched in recurrent patients at this time point. At post-treatment 12 months, more types of immune pathways were still activated in non-recurrent patients (C4) compared to recurrent patients (C5). In summary, the RFA therapy in HCC could improve various anti-tumor immune responses persisting for more than 9 months, which would be critical for recurrence prevention. In the meantime, the greater ability of T-cell differentiation in HCC patients was also one of the key factors determining whether recurrence occurred.

The univariable Cox regression analysis (Figure 4) showed that the T-cell differentiation signature score was significantly associated with RFS [p < 0.05 for gene ontology biological process (GOBP) regulation of T-cell differentiation in the thymus and GOBP positive regulation of CD4 alpha-beta T-cell differentiation]. The Kaplan–Meier analyses of pre-ablation treatment samples (n = 12) demonstrated that high expression of three distinct T-cell differentiation signatures correlated with improved RFS (Figures 8A–D). To further validate whether the T-cell differentiation signature was independently associated with recurrence, we performed multivariable Cox proportional hazards regression analyses (Supplementary Figure 3). After adjusting for age, gender, BCLC stage, and the Child–Pugh score, the four T-cell differentiation signatures did not reach statistical significance in the multivariable model. We acknowledged that statistical significance attenuated after multivariable adjustment, likely reflecting the small cohort size with limited events that reduced statistical power, and complete separation requiring Firth’s penalized likelihood correction, which yielded more conservative estimates. However, effect directions remained consistent, and the CD4 alpha-beta differentiation signature showed a particularly strong effect size (HR = 0.10, p = 0.060) approaching significance.

To externally validate the role of host T-cell fitness in predicting recurrence risk, we analyzed the four T-cell differentiation signatures in the TCGA-LIHC cohort (n = 371) via stratified survival analysis (Figure 9) and multivariate Cox regression (Figure 10). After multivariable adjustment for clinical confounders (sex, age, tumor stage, and grade), the gamma-delta T-cell, CD8 alpha-beta T-cell, and CD4 alpha-beta T-cell differentiation signatures maintained independent prognostic significance, confirming the robustness of these findings. Notably, while the primary cohort used PBMCs and the TCGA-LIHC cohort used tissue, requiring different cutoff values, the trends remained consistent.