Abstract The tumor microenvironment (TME) is a crucial mediator of tumor progression and treatment response. Here, we compare the immune microenvironments of HBV and non-HBV hepatocellular carcinoma (HCC) and investigate the reason for the persistence of HBV infection in the liver. We combine the Viral-Track method with scRNA sequencing and profile the transcriptomes of 70,056 cells from HBV and non-HBV-HCC patients. In addition to hepatocytes and macrophages, HBV transcripts were also detected in T and B cells using the Viral-Track method, confirming the lymphotropic nature of HBV in scRNA-sequencing data for the first time, to the best of our knowledge. HBV-HCC tumors have reduced levels of NK cells, macrophages, DCs, and increased malignant hepatocytes compared with those in non-HBV HCC. Notably, we report the enrichment of metallothioneins (MTs), particularly MT1G, in HBV-related HCC TAMs, which is associated with a worse prognosis. HBV-tumor-infiltrated CD8⁺ T cells exhibit a dysfunctional cytotoxic phenotype, characterized by upregulated MDK and CTLA4 expression and reduced IFN-γ production, unlike the non-HBV-HCC. Additionally, HBV-HCC exhibits immunosuppressive ligand-receptor interactions, whereas non-HBV-HCC exhibits antitumor ligand-receptor interactions. Our deeper understanding of the HBV-HCC ecosystem using Viral-Track integrated scRNA sequencing provides insights into immune evasion mechanisms and HBV lymphotropism associated with viral persistence. Subject terms: Tumour virus infections, Tumour heterogeneity __________________________________________________________________ Viral-Track integrated single-cell RNA-seq reveals HBV lymphotropism, TAM-specific MT1G upregulation and reduced IFN-γ production by CD8 + T cells, contributing to an immunosuppressive tumor microenvironment in HBV-associated hepatocellular carcinoma. Introduction Hepatitis B virus (HBV) infection is the most significant risk factor for the development of hepatocellular carcinoma (HCC)^[38]1. Despite sophisticated therapeutic interventions, the overall survival (OS) rate of HCC patients has not improved because viruses constantly evolve and adapt to overcome host defenses. Therefore, high-resolution analytical methods are needed to investigate viral dynamics, persistence, and virus-mediated tumor heterogeneity. The HCC tumor microenvironment (TME) is a complex ecosystem of tumor cells, activated hepatic stellate cells (HSCs), tumor-infiltrating lymphocytes (TILs), tumor-associated macrophages (TAMs), and stromal cells, promotes the immune escape, invasion, and metastasis of HCC^[39]2. CD8^+ cytotoxic T cells are indispensable for antitumor immunity^[40]3. In chronic HBV infection, severe dysfunction of HBV-specific CD8^+ T cells is frequently observed, characterized by a loss of cytokine secretion function (IL-2, TNF-α, and IFN-γ) due to prolonged antigen exposure^[41]4. Exhausted CD8^+ T cells frequently exhibited overexpression of various inhibitory receptors, including PD-1, LAG3, HAVCR2, and CTLA4, contributing to the failure of T-cell response^[42]5. Immune checkpoint inhibitors (ICIs) can reinvigorate T-cell function; however, some checkpoint inhibitors still exhibit acquired resistance during treatment^[43]6. It is necessary to investigate novel checkpoint inhibitors and other antineoplastic therapies that may benefit patients in the future. Several pioneering single-cell RNA sequencing (scRNA-seq) studies have investigated the immune or malignant cells of primary HCCs^[44]7, early-relapse HCCs^[45]8, and the differences between primary and metastatic HCCs^[46]9. However, these studies have not characterized a global landscape of TME combining HBV and non-HBV-HCC and their immunosuppressive characteristics of the cellular immune components. We utilized the computational tool Viral-Track, which has demonstrated the ability to comprehensively map scRNA-seq data onto a large database of known viral genomes, providing precise annotation of the cell types associated with viral infections^[47]10,[48]11. Few studies have suggested that lymphocytes are a site of HBV-host DNA integration, leading to persistent infection and potentially producing infectious virus particles, which supports the possibility of lymphotropism^[49]12,[50]13. In some viral infections, including HBV infection, lymphotropism may be essential, a primary mechanism of viral persistence, and a major contributor to disease pathogenesis^[51]14. However, evaluating HBV lymphotropism in liver tissues remains necessary using sophisticated molecular methods, such as single-cell RNA sequencing and the Viral-Track method. Here, we profiled the transcriptomes of 70,056 cells from HBV and non-HBV HCC patients. Using the Viral-Track method, HBV transcripts were also found in B and T cells in addition to the hepatocytes and macrophages, demonstrating the HBV-lymphotropic nature for the first time in scRNA-sequencing data, to the best of our knowledge. Notably, we found a subset of MT1G^+ M2 macrophages associated with a worse prognosis in HBV HCC patients. CD8^+ T cells in HBV tumors showed a cytotoxic-exhausted state, with lower IFN-γ expression. Our comprehensive understanding of the HCC ecosystem provides profound insights into immune evasion mechanisms, which may be beneficial in developing new immune therapy strategies for HBV and non-HBV tumors. Results Study design and Viral-Track method for identification and characterization of HBV infection in scRNA-Seq data To investigate the cellular diversity of the HBV and non-HBV-HCC tumor microenvironment, we collected surgical tumor specimens from five HBV and three non-HBV-HCC patients (discovery cohort) for scRNA-seq and the Viral-Track analysis (Fig. [52]1A; Supplementary Data [53]1). Moreover, we curated two publicly available scRNA sequencing datasets from the GEO database, [54]GSE234015 (validation cohort 1) and [55]GSE149614 (validation cohort 2), to substantiate our findings in the discovery cohort. In the discovery cohort, the quality control filtering was used to eliminate low-quality cells, indicated by a high correlation coefficient (r = 0.89) between the nCount_RNA and nFeature_RNA in the scatter plot, and 70,056 quality-assured singlet cells were recovered (Fig. [56]1B; Supplementary Fig. [57]1A, B). We observed a higher proportion of exonic reads, indicating our data’s good quality (Supplementary Fig. [58]1C). Moreover, each sample’s median of reads, genes, and UMI count distribution per cell favorably met the quality criterion (Supplementary Fig. [59]1D). Genes with substantial cell-cell variation (n  =  4000) were chosen for further downstream analysis (Fig. [60]1C). Principal component analysis (PCA) was performed on the cleaned dataset, and the first 12 PCs were included in the subsequent analysis (Supplementary Fig. [61]1E–G). Fig. 1. Study design and Viral-Track method for identification and characterization of HBV infection in scRNA-Seq data. [62]Fig. 1 [63]Open in a new tab A Schematic representation of the experimental strategy. B Results of doublet and singlet cells by scDblFinder. C The characteristic variance diagram was drawn for 4000 genes with significant cell-cell variations. D HBV transcripts identification and quality control metrics: (left panel) overall mapping quality of sequencing reads, including both host and viral transcripts; (middle panel). After removing human sequences, the percentage of mapped viral contigs and the entropy-based sequence complexity (right panel). E A viral genome coverage map of 5′ scRNAseq data from Viral-Track analysis shows the HBV gene locations from the NCBI database (left panel), and the HBV genome coverage map of 5′ scRNAseq data of validation cohort 1 (right panel). F Histogram shows the percentage (%) of infected and non-infected cells in each patient’s tumor tissue. G HBV reads in HCC tissues were detected by the Viral-Track method (left panel), and a pie chart of viral reads in each infected patient in the HBV tumour (right panel). H Uniform Manifold Approximation and Projection (UMAP) plot of the liver HBV data. Infected cells are red, and non-infected cells are green (left panel)-pie chart of infected and non-infected cells in HBV tumor (right panel). I The ratio of HBV-infected cells relative to the total cell populations defines the HBV transcript enrichment factor in each cell type. J Histogram of each analysed patient’s infected immune and non-immune cells in tumor tissue. See also Supplementary Figs. 1 and 2. In order to efficiently detect HBV-infected and non-infected cells from raw scRNA-seq data, we employed Viral-Track, an R-based computational pipeline. Briefly, Viral-Track relies on the STAR aligner^[64]15 to map the reads of scRNA-seq data to both the host reference genome and an extensive list of high-quality viral genomes^[65]10,[66]16. The majority of detected reads were mapped nearly 100% to the reference sequence 003977 (3182 nt) of the Hepatitis B virus based on two quality metrics: the percentage of mapped viral genome segments and the sequence complexity (Fig. [67]1D). Coverage plot analysis revealed a strong peak at the 5’ end of the C gene, encoding for the main core and HBx proteins, suggesting that the HBV virus actively produces virions, a similar pattern of peaks were also observed in the coverage plot of validation cohort 1 (Fig. [68]1E). The detection of HBV reads in tissue by Viral-Track further confirmed the HBV-HCC (n = 5, KC02A, KC06A, KC07A, KC08A, and KC09A) or non-HBV-HCC (n = 3, KC03A, KC05A, and KC10A) patients (Fig. [69]1F). The patient KC07A exhibited highest viral transcripts in HCC tissues, detected by Viral-Track method (Fig. [70]1G). We detected infected and non-infected cells using automated thresholding over the viral signal, and 23.57% of the cells in the HBV-related HCC scRNA-seq samples were found to be infected (Fig. [71]1H; Supplementary Fig. [72]2A–E). After calculating the frequency of HBV-infected cells in the total cell populations, the highest number of infected cells were detected in malignant hepatocytes, hepatic stellate cells (HSCs), and B cells, followed by myeloid and T/NK cells (Fig. [73]1I). A few studies using PCR-based techniques confirmed the productive replication of HBV in the immune cell compartments^[74]12,[75]14. Our findings indicated that the viral infection had a major effect on the transcriptional profile of malignant cells compared to T/NK and myeloid cells (Supplementary Fig. [76]2F). Next, the proportion of each infected immune and non-immune cell type was calculated in all HBV HCC patients (Fig. [77]1J). The landscape of the tumor ecosystem in HBV and non-HBV-HCC tumor tissues After the Viral-Track method annotation, we constructed a global tumor niche atlas of HCC tumors and performed cell classification and marker gene identification using Seurat^[78]17. We displayed 36 clusters by employing the uniform manifold approximation and projection (UMAP) algorithm^[79]18 (Fig. [80]2A; Supplementary Fig. [81]3A). All these clusters were shared among patients and between HBV and non-HBV-HCC samples (Fig. [82]2B, C; Supplementary Fig. [83]3B). The infection map was overlaid onto the annotated cell map. It showed that most HBV transcripts were observed within the malignant hepatocyte clusters, suggesting the efficiency of the Viral-Track method (Fig. [84]2D). All cells (both immune and non-immune fractions) expressed high levels of housekeeping genes, such as ACTB, B2M, GAPDH, and LDHA (Supplementary Fig. [85]3C), whereas immune cells expressed PTPRC (CD45)^[86]19 (Fig. [87]2E), indicating the accuracy of our data. Cluster 22 had the most up-regulated 1289 genes, while Cluster 1 had the most down-regulated 2,894 genes (Supplementary Fig. [88]3D). We then annotated all these clusters with canonical marker genes of major cell types. The identified non-immune cells mainly consist of eighteen clusters of malignant hepatocytes (ALB, TTR, APOB, and GPC3), one cluster of endothelial cells (CDH5, CLDN5, EPCAM, and VWF), and one cluster of HSCs (ACTA2, COL1A2, and PDGFRB). Immune cells consisted of seven clusters of T&NK cells (CD3D, CD3G, TRAC, TRBC1, and KLRD1), seven clusters of myeloid cells (LYZ, CD68, and C1QA/B), and two clusters of B cells (CD79A, MZB1, IGHG1, and JCHAIN) (Fig. [89]2E, F; Supplementary Fig. [90]3E). Additionally, we displayed each cluster’s top marker gene (0–35) (Supplementary Fig. [91]4A, B). Fig. 2. Landscape of the tumor ecosystem in HBV positive and HBV negative HCC tumor tissues. [92]Fig. 2 [93]Open in a new tab A UMAP plot, showing the annotation for cluster number and colour codes for cell type in the HCC ecosystem. B The UMAP plot, showing patient IDs by colour (left panel); the bar graph represents patient cell counts (right panel). C UMAP plot, showing tumor type by colour, HBV(blue), or non-HBV HCC(green). D UMAP plot coloured by viral status, HBV-infected cells (green), and non-infected cells (red). E UMAP plot showing the expression levels of marker genes, defined for all cell types. F Scaled heatmap of marker gene expression in the indicated cell types. The top bars label cell-type clusters, and the number in brackets matches the cluster number in (A). G Histogram indicating the immune and non-immune cell fraction of each patient’s tumor tissue. H Pie chart of immune and non-immune cell subtypes in HBV and non-HBV-HCC tissues. See also Supplementary Figs. [94]3, [95]4, and [96]5. The copy number variants (CNVs) revealed amplifications of chromosomes number 1 and 6 in non-HBV-HCC and deletion of chromosomes 5 and 6 in HBV-HCC. These findings align with previous studies^[97]20–[98]22 (Supplementary Fig. [99]5A). Next, the proportion of immune and non-immune cells in the HBV and non-HBV-HCC patients was calculated (Fig. [100]2G). Notably, HBV-associated HCC samples revealed a substantially higher number of malignant hepatocytes and a lower number of myeloid and T&NK cells (Fig. [101]2H). In summary, HBV and non-HBV-HCC samples shared the same major immune and non-immune cell subtypes but variable proportions of malignant, T/NK, and myeloid-derived cells (Fig. [102]2G, H; Supplementary Fig. [103]5B, C). Prevalence of metallothioneins enriched TAMs in HBV-mediated HCC microenvironment We reclustered myeloid cells and identified 16,466 cells divided into eight clusters: five for macrophages (Mø), one for monocytes, and one for DCs (Fig. [104]3A–C). By employing the Viral-Track method, HBV transcripts were detected in macrophages, which are known to engulf dead hepatocytes and potentially HBV virions^[105]11 (Fig. [106]3D). All myeloid cells exclusively exhibited LYZ, C1QB, and FTL marker genes. CD163 and CD68 expression identified macrophages, while MKI67 and HMGB2 expression identified proliferating macrophages (Fig. [107]3E; Supplementary Fig. [108]6A). Mø_c0 (40.23%), Mø_c1 (26.72%), and Mø_c2 (16.26%), which were non-HBV-HCC clusters, demonstrated increased expression of CCL3L1, CCL3, CCL4, and CXCL3, signatures of inflammatory TAMs^[109]23. Notably, Mø_c3 (12.18%) represented TAMs from HBV-HCC expressed higher levels of the metallothionein genes MT1G, MT1A, and MT2A, along with MMP-9 that were distinct from non-HBV subpopulations (Fig. [110]3E; Supplementary Fig. [111]6A, C). MMP-9 has been shown to interact with MTs in tumor progression and angiogenesis^[112]24. Mø_c7 (4.61%) were proliferating macrophages predominantly in non-HBV-related HCC samples (Supplementary Fig. [113]6B, C). Fig. 3. Prevalence of metallothioneins enriched TAMs in HBV-mediated HCC microenvironment. [114]Fig. 3 [115]Open in a new tab A UMAP plot of subclustered myeloid cells labelled in different colours. Cluster number and cell type annotations are provided in the Figure. B The UMAP plot, showing patient IDs by colour (left panel); the bar graph represents the cell counts of each individual (right panel). C UMAP plot, showing tumor type by colour, HBV(blue), or non-HBV-HCC(green) myeloid. D UMAP coloured by viral status, HBV infected (green), and non-infected cells (red). E Expression of marker genes correlated with myeloid subtypes (macrophages, monocytes, and DCs). F Scatterplots showing M1 and M2 polarization scores (top panel) and pro- and anti-inflammatory scores (bottom panel), colour-coded according to patient ID. G Boxplot showing the fraction of myeloid cell subtypes in HBV (purple) and non-HBV (green) (left panel), violin plot indicating the fraction of myeloid cell subtypes in HBV (purple) and non-HBV (green) samples from the Validation cohort 2 (right panel). Error bars represent mean ± standard deviation (SD). H Representative IHC staining, indicating CD68^+ cells, in paired HBV-HCC and non-HBV samples. Scale bar, 20 μm. I The volcano plot shows differentially expressed genes between HBV (red dots) and non-HBV myeloid (blue dots). The names of the most significant genes are indicated in the plot. J Violin plots indicating the expression of metallothionein genes (MT1G and MT1H) from macrophage clusters between HBV(blue) and non-HBV tumors (green) (top panel) violin plots showing the MT1G and MT1H genes expression in HBV (blue) and non-HBV (green) samples from the Validation cohort 2 (bottom panel). K Violin plots indicate the average metallothionein gene expression score from HBV (blue) and non-HBV (green) macrophage clusters (left panel) or infected and non-infected macrophages (middle panel); violin plots depicting the average metallothionein gene expression score in HBV (blue) and non-HBV (green) samples from the Validation cohort 2(right panel). L Representative IHC staining, indicating CD68^+MT1G^+ macrophages, in paired HBV and non-HBV HCC samples. Scale bar, 20 μm. The p-values in (J) and (K) were calculated using Student’s t-test; in (G), using the unpaired Wilcoxon test. All statistical comparisons were performed on biologically independent samples (n = 5 HBV-HCC and n = 3 non-HBV-HCC) from both the discovery and validation cohort 2. See also Supplementary Figs. [116]6 and [117]7. Established marker genes, such as FCGR3A (M1) and CD163 (M2), were expressed in both M1 and M2 macrophages, making it difficult to distinguish between them (Fig. [118]3E), as previously reported^[119]8. However, by calculating M1 and M2 polarization and pro- and anti-inflammatory scores, we observed a dominant M2-like phenotype concomitant with prior studies^[120]25 (Fig. [121]3F; Supplementary Fig. [122]6D; Supplementary Data [123]2). None of these subtypes was linked with liver-resident Kupffer cells (MARCO^+) or myeloid-derived suppressor cells (MDSCs, CD33^+)^[124]8 (Supplementary Fig. [125]6A). We found elevated ISG15 and ISG20 expression in HBV-HCC linked to HBV-related diseases, while elevated levels of TGF-β, IDO1, and ARG1 demonstrated M2 phenotype in non-HBV-HCC (Supplementary Fig. [126]6E). Notably, we observed macrophages, proliferating macrophages, and DCs were abundant in HBV-negative HCC tumours compared to those with HBV infection; similarly, the scRNA-seq data for the five additional HBV HCC and three non-HBV-HCC samples in the validation cohort 2 also demonstrated increased myeloid fraction in non-HBV-HCC, as previously depicted in a study on HBV-positive and HBV-negative HCC^[127]26 (Fig. [128]3G; Supplementary Fig. [129]6F; Supplementary Data [130]12). IHC analysis verified the abundance of macrophages (CD68^+) in non-HBV tumors (Fig. [131]3H). We identified CD1c^+/CD14^− type 2 conventional dendritic cells (cDC2) in cluster 5, mainly comprised of non-HBV HCC samples, previously reported to be responsible for the induction of IFNγ production in CD4^+ T cells^[132]27 (Fig. [133]3E; Supplementary Fig. [134]6A, B). Differentially expressed gene (DEG) analysis indicated that the top upregulated genes in HBV vs non-HBV-HCC myeloid cells included MT1G and MT1H (Fig. [135]3I; Supplementary Data [136]3). Similar results were obtained in infected vs non-infected myeloids (Supplementary Fig. [137]7A; Supplementary Data [138]4). KEGG pathway indicated that up-regulated DEGs in HBV-HCC were predominantly associated with metallothionein-mediated mineral absorption (Supplementary Fig. [139]7B). Further, we explicitly reported the enrichment of metallothioneins (MTs) expression in TAMs found in HBV-related HCC, supported by the validation cohort 2 (Fig. [140]3J, K; Supplementary Fig. [141]7C). Moreover, the Cancer Genome Atlas Liver Hepatocellular Carcinoma (TCGA-LIHC) dataset showed high expression of MT1G in HBV-HCC compared to non-HBV-HCC and correlated with adverse clinical outcomes in Hepatitis virus-mediated HCC patients (Supplementary Fig. [142]7D). A dual-colour IHC assay verified the enrichment of CD68^+MT1G^+ cells in HBV-related HCC samples (Fig. [143]3L). Investigating T/NK cells functional diversity and HBV lymphotropism in HBV-mediated HCC We identified a total of 10,919 T/NK cells that were divided into 14 clusters shared across patients between HBV and non-HBV-HCC tumor samples (Fig. [144]4A; Supplementary Fig. [145]8A). All eleven T cell subtypes were identified by the expression of CD3D, CD3E, CD3G, TRAC, and TRBC2 (Fig. [146]4B, C), including two subtypes of CD4^+FOXP3^+ Tregs: Treg TNFRSF4 in cluster 5 and Treg TNFRSF18 in cluster 11, which were tumor-infiltrating Tregs in non-HBV-HCC. One subtype of CD62L^−CCR7^+IL7R^+ phenotype memory T cells (cluster 0, 19.97%) comprised mainly HBV-HCC, one subtype of IL17A^+, IL17F^+, RORA^+ T helper type 17 (Th17) phenotype (cluster 10, 6.39%) predominantly composed of HBV-HCC and previously associated with a worse prognosis in the HBV-HCC^[147]28, one subtype of MKI67^+STMN1^+ T cycling (cluster 12, 4.82%) were abundant in non-HBV-HCC. Additionally, there were four subtypes of CD8^+ T cells, with two subtypes, CD8^+ MDK (cluster 1, 16.57%) and CD8^+ CCL3L1 (cluster 2, 15.19%), being primarily associated with HBV-HCC. MDK expression has been previously reported in activated T cells and HIV-infected cells^[148]29. In non-HBV-HCC, the remaining two subtypes, CD8^+ CLU (Cluster 4; 12.71%) and CD8^+ CCL4L2 (cluster 8, 7.17%), were identified as having enriched expression of IFNG (Fig. [149]4A, C; Supplementary Fig. [150]8B, C). Furthermore, a significant expression of αβ TCRs and a substantial reduction in γδ TCR markers (TRGV, TRDV, TRDC, and TRGC) were observed in all clusters of T cells in HBV and non-HBV-HCC. This observation is consistent with previous findings, where an assessment of TRGV, TRGC, and TRDV variant expression patterns in γδ T cell clusters demonstrated a universal loss of TCR diversity in HCC livers. The loss of TCR diversity is inversely correlated with tumor progression^[151]30 (Fig. [152]4B; Supplementary Fig. [153]8D–F). Fig. 4. Investigating T/NK cells functional diversity and HBV lymphotropism in HBV-mediated HCC. [154]Fig. 4 [155]Open in a new tab A UMAP plot of HBV and non-HBV-tumor T/NK subsets. Cluster number and cell type annotations are provided in the Figure. B Expression of marker genes correlated with T/NK subtypes (αβ T cells, γδ T cells, T memory, Tregs, CD8^+ T, and NK cells). C T/NK subtype selected gene expression heatmap for naive, resident, inhibitory, cytokines, co-stimulatory, transcriptional factors (TF), and cell type. D Boxplots showing the ratio of KLRF1^+, KLRD1^+, and KLRB1^+ NK cells in HBV(purple) and non-HBV (green) tumor patients (left panel), violin plots indicating the expression of KLRD1 in NK cells between HBV(blue) and non-HBV HCC (green) samples from the Validation cohort 2 (right panel). E UMAP colored by viral status, representing HBV-infected (green) and non-infected (red) T cells (left panel). F UMAP colored by viral status, representing HBV-infected (green) and non-infected (red) B cells (left panel); UMAP colored by viral status, representing only HBV-infected (green) B cells (right panel). G A coverage map of the HBV genome. The locations of the HBV genes, as listed in the NCBI database, are depicted at the bottom. H Representative IHC staining, showing CD3^+HBc^+ T cells in paired HBV and non-HBV HCC samples. Scale bar, 20 μm. I Boxplot showing the ratio of NK cells in HBV and non-HBV tumor patients. J Boxplot showing the ratio of CD4^+, CD8^+ T, and Th17 cells in HBV and non-HBV-HCC patients. Error bars represent mean ± standard deviation (SD). The p-values in (D, I, and J) were calculated using the unpaired Wilcoxon test. All statistical comparisons were performed on biologically independent samples (n = 5 HBV-HCC and n = 3 non-HBV-HCC) from both the discovery and validation cohort 2. See also Supplementary Figs. [156]8, [157]9 and [158]10. Two NK subtypes (CD16 and CD160, predominantly non-HBV-HCC) were identified by high expression of KLRF1, KLRC1, and KLRD1 in clusters 6 and 9 (Fig. [159]4B, C; Supplementary Fig. [160]8B). We observed higher fractions of KLRF1^+, KLRD1^+, and KLRB1^+ NK cells in non-HBV-HCC; similarly, higher fractions of KLRD^+ were observed in validation cohort 2, exhibited a mature phenotype^[161]31, consistent with the patients’ favourable survival in non-viral HCC using the TCGA database (Fig. [162]4D; Supplementary Data [163]12; Supplementary Fig. [164]9). NKT cells (cluster 2, non-HBV-HCC) were characterized by the expression of both NK and T cell markers (Fig. [165]4A–C). Although some early PCR-based studies on chronic Hepatitis B virus infection revealed that the virus might reside in B and T cells other than hepatocytes, suggesting HBV lymphotropism is a primary mechanism of viral persistence and pathogenesis^[166]12–[167]14. However, a need remains to evaluate HBV in the liver tissues using advanced molecular techniques such as scRNA sequencing. In this study, for the first time, HBV transcripts have been comprehensively detected in T/NK and B cells within HBV-HCC tumor samples using the Viral-Track method in scRNA sequencing data (Fig. [168]4E–G; Supplementary Fig. [169]10A). To ensure the robustness of our findings and rule out technical artifacts, we performed a mutually exclusive marker analysis within our HBV-HCC dataset. Specifically, we examined the co-expression of CD3D (T cell marker) with hepatocyte-specific markers (TTR and GPC3), as well as IGHG1 (B cell marker). This analysis confirmed that the vast majority of CD3D+ T cells did not co-express hepatocyte or B cell markers, and importantly, HBV-infected T cells showed negligible to no co-expression with these non-T cell lineage markers. These results validate the effectiveness of our doublet exclusion approach and support the conclusion that the detected HBV signal in lymphoid cells is not an artifact of cell doublets (Supplementary Fig. [170]10B–H). The coverage plot indicated that the HBV was actively producing virion in T/NK cells, as evidenced by a prominent peak at the 5′ end of the C gene, which encodes the main core and HBx proteins, and verified the existence of productive replication of HBV in the immune cell compartment, suggesting lymphoid cells act as a reservoir of replicating HBV likely contribute to the persistence of hepatitis B virus (Fig. [171]4G). Dual-colour IHC verified the enrichment of CD3^+HBc^+ T cells in HBV-HCC samples, further supporting the HBV lymphotropism (Fig. [172]4H). In summary, HBV-HCC exhibits a distinct immune ecosystem characterized by HBV lymphotropism and a diminished proportion of NK cells, T cycling, fractionally higher CD8^+ T cells, and significantly increased proportions of Th17 in comparison to the non-HBV-HCC microenvironment, whereas there was no significant difference in CD4^+ Tregs proportion in both HBV and non-HBV-HCC (Fig. [173]4I, J; Supplementary Fig. [174]10I; Supplementary Data [175]12). Cell-state transition trajectory analysis of CD8^+ tumor-infiltrating T cell populations in HBV and non-HBV-HCC We used Monocle to determine the pseudo-time cell trajectories of CD8^+ T cells to understand the cellular state evolution in HBV and non-HBV-HCC TME^[176]8 (Fig. [177]5A, B; Supplementary Fig. [178]11A). We then examined the transcriptional alterations using a branching heatmap; after combining the findings from the clustering and pseudo-time analyses, we observed intriguingly that CD8^+ T cells could differentiate into three pseudotime-dependent states and eventually shifted toward subpopulations with exhaustion status (Fig. [179]5C, D; Supplementary Fig. [180]11B). The independent trajectories analysis of CD8^+ T cells in HBV and non-HBV-HCC revealed that early-state CD8^+ T cells were found in non-HBV-HCC, whereas HBV-HCC CD8^+ T cells predominated during the second state. In the third state, CD8^+ T cells in HBV-HCC were predominant at the initial ends, while non-HBV-HCC CD8^+ T cells were primarily located at the terminal ends (Fig. [181]5E, F). CD8^+ CLU cells were predominantly state one cells in non-HBV-HCC, characterized by upregulated expression of FASN, FABP1, HMGCS2, and low expression of cytotoxic markers GZMB, NKG7, and GNLY, whereas second state comprised CD8^+ MDK and partially CD8^+ CCL3L1 cells were found in HBV-HCC (Fig. [182]5C). These CD8^+ MDK cells were deficient in interferon-gamma (IFN-γ), indicating a dysfunctional phenotype, consistent with prior findings of IFN-γ downregulation in MDK-enriched malignancies^[183]32 (Supplementary Fig. [184]11C). State 3 exhibited both cytotoxic (GZMB, GZMA, GZMH, NKG7, and IFNG) and exhausted markers (PDCD1, CXCL13, LAG3, TIGIT, CTLA4, and HAVCR2), demonstrating dysfunctional CD8^+ T cells (Fig. [185]5G, H; Supplementary Fig. [186]11C). Fig. 5. Cell-state transition trajectory analysis of CD8^+ tumor-infiltrating T cell populations in HBV and non-HBV-HCC. [187]Fig. 5 [188]Open in a new tab A Pseudotime-ordered analysis of CD8^+ T cells from HBV and non-HBV tumor samples. B CD8^+ T cell subtypes are labelled by colors. C Heatmap showing the dynamic changes in gene expression along the pseudotime (lower panel). The distribution of CD8 subtypes during the transition (divided into 3 states), along with the pseudo-time. Subtypes are labeled by colors (upper panel). D A 2D pseudo-time plot shows the dynamics of cytotoxic (upper panel) or exhausted signals (lower panel) in CD8^+ T cells in HBV and non-HBV tumor samples (left panel). Violin plots showing the average score of cytotoxic and exhausted signals from the clusters of CD8^+ T cells in HBV and non-HBV tumor samples (right panel). E 2D pseudo-time-ordered CD8^+ T cell graph from non-HBV (left panel) and HBV (right panel) tumor tissues. State-specific cell density distribution is shown at the top of the Figure. F Histogram of CD8^+ T lymphocytes in three states in HBV and non-HBV tumor tissues. CD8 subtypes labeled by colors. G Two-dimensional plots of average scores for cytotoxic (upper panel) and exhausted signals(lower panel) in non-HBV (green) and HBV (blue) HCC samples, dynamically along the pseudotime. H Two-dimensional plots of exhausted marker genes (CXCL13, LAG3, CTLA4, TIGIT, PDCD1, and HAVCR2) during the T cell pseudotime transitions. I Two-dimensional plots showing the expression scores for genes related to oxidative phosphorylation, fatty acid metabolism, and PPAR signalling pathway in non-HBV (green) and HBV (blue) HCC samples. J Two-dimensional plots showing the dynamic expression of exhaustion transcription factors (EOMES, TOX, HOPX, IRF2, and IRF9) along the pseudo time. See also Supplementary Fig. [189]11. According to GO and KEGG pathway analysis of differentially expressed genes (DEGs) along the pseudo-time, state 1 DEGs were strongly associated with oxidative phosphorylation, fatty acid, and cholesterol metabolism, and the PPAR signaling pathway (Fig. [190]5I; Supplementary Fig. [191]11D; Supplementary Data [192]5). Upregulation of the PPAR signaling pathway enhances FAO to fuel the CD8^+ T cell’s energy demands and T cell effector functions in the TME^[193]33. State 2 CD8^+ T cells involved in cysteine, methionine, glutathione, and carbon metabolism pathways linked to T cell proliferation and survival^[194]34 (Supplementary Fig. [195]11E). Glutathione metabolism proteins GPX2, GPX3, and GSTA1 increased significantly in activated T cells to buffer reactive oxygen species (ROS)^[196]35 (Fig. [197]5C). CD8^+ T cells in state 3 were associated with the chemokine signaling, PD-1 checkpoint, HIF-1, RAS, TNF, and VEGF pathways (Supplementary Fig. [198]11F; Supplementary Data [199]5). State 3 was characterized by increased expression levels of chemokines (CCL3L1, CCL4L2, CCL4, and CCL5) (Supplementary Fig. [200]11G), classical cytotoxic (GZMB, GZMA, GZMH, NKG7, and IFNG), and T cell exhaustion markers (PDCD1, CXCL13, LAG3, TIGIT, CTLA4, and HAVCR2), suggesting exhausted cytotoxic CD8^+ T cell state (Fig. [201]5C, H; Supplementary Fig. [202]11C). In addition, state 3 was characterized by high levels of T cell exhaustion-related transcriptional factors TFs, including EOMES, TOX, HOPX, IRF2, IRF9, and the programmed cell death protein 1 (PD-1) pathway, suggesting the exhausted state of these cells (Fig. [203]5C, J; Supplementary Fig. [204]11H). Overall, we conclude that the transition trajectories, metabolic pathways, and expression of immune checkpoints and cytotoxic markers in HBV or non-HBV-HCC CD8^+ T cells were not functionally homogeneous, as distinct T cell subsets have different degrees of dysfunction or impairment of effector functions. CD8^+ T cells in HBV-HCC exhibited a dysfunctional phenotype characterized by MDK expression and diminished IFN-γ production, whereas non-HBV-HCC retained higher cytotoxic activity. HBV-related HCC tumor infiltrating MDK^highIFN-γ^lowCD8^+ T cells displayed a dysfunctional immunosuppressive phenotype We further compared and analyzed the expression of the cytotoxic and exhausted markers of CD8^+ T cells in HBV and non-HBV-HCC. LAG3 and CTLA4 expression facilitated exhaustion in HBV-HCC, whereas CXCL13 contributed to exhaustion in non-HBV-HCC CD8^+ T cells. Additionally, cytotoxic marker genes such as GZMA, PRF1, GNLY, and CST7 substantially increased, while IFN-γ dramatically decreased in HBV-HCC; intriguingly, similar findings were observed in the validation cohort 2 (Fig. [205]6A; Supplementary Fig. [206]12A, B). The loss of IFN-γ production occurs in highly exhausted CD8^+ T cells^[207]36. Overall, the cytotoxic and exhausted marker expression scores were more pronounced in CD8^+ T cells of HBV-HCC samples, demonstrating an exhausted cytotoxic CD8^+ T cells phenotype validated by validation cohort 2 (Fig. [208]6B). Then, DEG analysis of HBV versus non-HBV CD8^+ T cells indicated that MDK and GNLY were highly up-regulated genes, while IFNG was down-regulated, indicating an immune-suppressive TME with a dysfunctional cytotoxic phenotype^[209]32,[210]37 (Fig. [211]6C; Supplementary [212]Data 6). Previous research suggests that low IFN-γ levels produced by host-infiltrating cells at the tumor site increase tumor cell survival and metastatic potential^[213]38. To assess the effector capacities of CD8^+ T cells, we examined tumor-infiltrated CD8^+ IFN-γ^high (log[2] expression>0.585) and IFN-γ^low (log[2] expression<0.585) cells in HBV and non-HBV-HCC samples. Our study confirmed that IFN-γ^high cells significantly expressed cytotoxic markers such as GZMA, PRF1, and GZMB, supporting earlier findings^[214]32 (Fig. [215]6D). Dual-colour IHC analysis confirmed the reduced level of CD8^+IFNG^+ in HBV compared to non-HBV tumors (Fig. [216]6E), and staining specificity was validated using non-tumor liver tissue, which showed minimal IFN-γ signal under the same conditions (Supplementary Fig. [217]12C). Further, we observed that HBV-HCC had a higher expression of MDK^+CD8^+ T cells, strongly supported by validation cohort 2 (Fig. [218]6F). According to TCGA-LIHC analyses, an MDK^highCD8^+ T cell signature (top 50 genes) predicted worse clinical outcomes in virus-mediated HCC patients (Fig. [219]6G; Supplementary Data [220]7). High MDK expression was also associated with poor clinical outcomes in the TCGA-based data analysis of hepatitis virus-mediated HCC patients (Fig. [221]6H). HCC patients with elevated levels of Midkine have a poor prognosis and immunosuppressive tumor microenvironment^[222]39. Further, we evaluated the expression of classical exhausted (CTLA4, HAVCR2, PDCD1, and LAG3) and cytotoxic markers (GZMA, PRF1, and IFNG) in CD8^+ CTLA4^high (log[2] expression>0.585) and CTLA4^low (log[2] expression<0.585) cells. In HBV-HCC, CD8^+CTLA4^high cells had higher exhaustion scores than non-HBV-HCC. However, the cytotoxic score was not statistically significant; similar findings were observed in the validation cohort 2 (Fig. [223]6I, J). Similarly, both cytotoxic and exhaustion markers were notably higher in HBV-HCC LAG3^high (log[2] expression >0.585) CD8^+ T cells, as supported by validation cohort 2 (Supplementary Fig. [224]12D, E). Fig. 6. HBV-related HCC tumor infiltrating MDK^highIFN-γ^lowCD8^+ T cells displayed a dysfunctional immunosuppressive phenotype. [225]Fig. 6 [226]Open in a new tab A Violin plots indicating expression of exhausted genes in CD8^+ T cells (LAG3, CTLA4, and CXCL13) (top panel) or cytotoxic genes (GZMA, PRF1, and IFNG) (middle panel) in HBV (blue) and non-HBV (green) tumor samples, violin plots depicting the expression of LAG3, CTLA4, CXCL13 and IFNG in HBV (blue) and non-HBV (green) samples from the Validation cohort 2 (bottom panel). B Violin plot indicating the cytotoxic and exhausted average scores of CD8^+ T cells in HBV(blue) and non-HBV (green) tumor samples in the discovery (left panel) and validation cohort 2 (right panel). C The volcano plot shows differentially expressed genes between HBV (red dots) and non-HBV tumor-infiltrated CD8^+ T cells (blue dots). The names of the most significant genes are indicated in the plot. D Violin plot indicating the expression of cytotoxic marker genes in IFNG^high (log[2] expression>0.585) and IFNG^low (log[2] expression[<]0.585) CD8^+ T cells across HBV and non-HBV tumor samples. E Representative IHC staining, showing CD8^+IFNG^+ T cells in paired HBV and non-HBV-HCC samples. Scale bar, 20 μm. F Violin plots indicating the expression of MDK cells in HBV (blue) and non-HBV (green) tumor samples in the discovery cohort (left panel) and validation cohort 2 (right panel). G The Kaplan-Meier plot of disease-free survival in the TCGA LIHC cohort compares low versus high MDK^high (signature of top 50 genes) in hepatitis virus HCC patients. The plot indicates the number of patients from each subgroup, the Log-Rank t-test results, and the number of patients (N). H The Kaplan-Meier plot of disease-free survival in the TCGA database compares low versus high MDK in Hepatitis (left panel) and non-hepatitis virus (right panel) HCC patients. The plot shows patient numbers from each subgroup, the Log-Rank t-test results, and the number of patients (N). I Violin plot indicating the expression of exhausted (upper panel) and cytotoxic marker genes(lower panel) in CTLA4^high (log[2] expression>0.585) and CTLA4^low (log[2] expression<0.585) CD8^+ T cells across HBV and non-HBV tumor samples. J Violin plots showing the average score of cytotoxic and exhausted signals in CTLA4^high (log[2] expression>0.585) and CTLA4^low (log[2] expression<0.585) CD8^+ T cells across HBV and non-HBV tumor samples in both the discovery cohort (left panel) and the validation cohort 2(right panel). The p-values in (A, B, D, F, I, and J) were calculated using Student’s t-test. All statistical comparisons were performed on biologically independent samples (n = 5 HBV-HCC and n = 3 non-HBV-HCC) from both the discovery and validation cohort 2. See also Supplementary Fig. [227]12. In conclusion, CD8^+ T cells in HBV-HCC samples showed a cytotoxic dysfunctional state, characterized by upregulated MDK and CTLA-4 expression and reduced IFN-γ production, unlike non-HBV-HCC. Highly proliferative malignant hepatocytes in HBV-related HCC foster more immunosuppressive TME than non-HBV-HCC Our results showed that HBV-HCC tumor microenvironments had significantly lower immune cell subtype fractions in Fig. [228]2H. This prompted us to hypothesize that HBV and non-HBV-HCC malignant hepatocytes might differ in immunosuppression and proliferation characteristics. To address this assumption, 40,376 hepatocytes were arranged into 18 clusters shared across HBV and non-HBV-HCC patients (Fig. [229]7A, B). The Viral-Track method comprehensively detected a higher number of HBV transcripts across all hepatocytes, which serve as infection reservoirs (Fig. [230]7C). Each subtype of malignant hepatocytes exhibited unique gene expression profiles (Supplementary Fig. [231]13A). Cluster 1, enriched with immunosuppressive MDK expression^[232]37,[233]39, and Cluster 12, abundant with cycling malignant hepatocytes (MKI67^+, STMN1^+, and HIST1H4C^+), were identified in HBV-HCC (Fig. [234]7A, D). Interestingly, proliferative malignant hepatocytes belonged to patient KC07A, which has the highest viral load, suggesting that viral load enormously affects malignant hepatocyte proliferation^[235]40 (Fig. [236]1D; Supplementary Fig. [237]13B, C; Supplementary Data [238]12). We identified malignant hepatocytes by the marker genes such as ALB, TTR, GPC3, and APOA2; however, AFP’s expression pattern indicates that it may not be a reliable HCC marker, consistent with prior studies^[239]41 (Fig. [240]7E). Notably, malignant hepatocytes were more prevalent in HBV-HCC (Fig. [241]7F; Supplementary Data [242]12). Next, we performed a DEG analysis in HBV versus non-HBV malignant hepatocytes. We observed an enrichment of genes involved in immune response pathways (e.g., IFN_gamma signaling) in non-HBV-HCC. In contrast, in HBV-HCC, the upregulated genes primarily belonged to cell-cycle (e.g., MYC_targets_V1) and hypoxia-related pathways, suggesting an immunosuppressive microenvironment (Fig. [243]7G; Supplementary Fig. [244]13D; Supplementary Data [245]8). Non-HBV-HCC presented higher immune surveillance scores, whereas HBV-HCC had higher metastatic potential, proliferative, and immune escape scores, as validated by validation cohort 2 (Fig. [246]7H; Supplementary Fig. [247]13E–G; Supplementary Data [248]9). High expression levels of genes related to immune surveillance, HLA/B/C, suggest that MHC-I antigen presentation was dominant in non-HBV-HCC tumors (Supplementary Fig. [249]13H). Further, we found enhanced expression of immune escape genes EBAG9, ADAM10, PTGES3, LGALS4, NECTIN-2, MDK, CD47, and TGFB1 in HBV-HCC, confirmed by validation cohort 2^[250]37,[251]42 (Fig. [252]7I; Supplementary Fig. [253]13I). Fig. 7. Highly proliferative malignant hepatocytes in HBV-related HCC foster more immunosuppressive TME than non-HBV-HCC. [254]Fig. 7 [255]Open in a new tab A UMAP plot of re-clustered malignant cells derived from HBV and non-HBV liver tumor. Cluster number annotations are provided in the Figure. B UMAP plot showing malignant cells, color-coded according to patient origin (left panel). The UMAP plot showing tumor type by color, HBV, or non-HBV HCC (right panel). C UMAP colored by viral status in malignant cells, HBV-infected (green), and non-infected cells (red). D Sankey diagram showing the percentages of malignant cell clusters among HCC patients and vice versa. E Expression of marker genes correlated with malignant hepatocytes (ALB, TTR, APOA2, GPC3, SERPINA1, and AFP). F Boxplot showing the fractions of malignant cells in HBV (blue) and non-HBV tumor patients (green). Error bars represent mean ± standard deviation (SD). G Based on the HALLMARK gene set of upregulated genes, the bar chart showing the enrichment of specific pathways in HBV (upper panel) and non-HBV (lower panel) malignant cells. H Violin plot indicating the average score of metastatic, immune escape, immune surveillance, and proliferation signals in malignant cells across HBV (blue) and non-HBV (green) tumor samples in discovery (left and middle panel) and validation cohort 2 (right panel). I Violin plot showing selected immune escape genes (MDK, CD47, and TGFB1) in tumor cells from HBV (blue) and non-HBV (green) HCC patients in discovery (upper panel) and validation cohort 2 (lower panel). The p-values in (H and I) were calculated using Student’s t-test, and in (F) by the unpaired Wilcoxon test. All statistical comparisons were performed on biologically independent samples (n = 5 HBV-HCC and n = 3 non-HBV-HCC) from both the discovery and validation cohort 2. See also Supplementary Fig. [256]13. In brief, our findings revealed that a greater proportion of malignant hepatocytes in the HBV-HCC possessed enhanced capacities for metastasis, cellular proliferation, immune evasion, and decreased immune surveillance. Deciphering ligands-receptors interactions in HBV and non-HBV tumor microenvironment To systematically dissect the ligand-receptor (L-R) interactions in the tumor ecosystem, we conducted a cell-cell communication network based on CellPhoneDB (Fig. [257]8A). Our results demonstrated that the HBV possessed a significantly greater number of ligand-receptor (L-R) pairs, 11,251 (p < 0.05), than the non-HBV-HCC, 6383 (p < 0.05) (Supplementary Data [258]10, [259]11). In HBV and non-HBV-HCC, myeloid, malignant, and T/NK cells had more potential ligand-receptor pairings than HSCs, endothelial, and B cells (Fig. [260]8B; Supplementary Fig. [261]14A). Fig. 8. Deciphering ligand-receptor interactions in HBV and non-HBV-HCC microenvironment. [262]Fig. 8 [263]Open in a new tab A A directed network demonstrating the differential ligand-receptor (L-R) interaction intensities among the six major cell types (colour-coded by cell type) in HBV (left panel) and non-HBV (right panel) tumor tissues. Edge width shows L-R interaction intensities. B Bar chart showing the number of significant ligand-receptor pairs among myeloid (top panel), malignant cells (bottom panel), and all other cells in HBV and non-HBV tumor samples. C Dotplot showing HBV and non-HBV-HCC prevalent L-R interactions among malignant, T/NK, and myeloid cell types, based on significant ligand and receptor pairs. These scores are normalized expression levels, dot size indicates the significance of the interaction, and dot color indicates interaction intensity, calculated by CellPhoneDB. D A Kaplan-Meier plot of disease-free survival in the TCGA LIHC cohort, comparing low versus high MDK-LRP1 ligand-receptor pair’s mean expression in hepatitis virus (right panel) and non-hepatitis virus (left panel) HCC patients. The number of patients from each subgroup and the Log-Rank t-test results are indicated in the plot. E Dotplot showing the interaction between T/NK clusters and myeloid cells in HBV and non-HBV-HCC, based on selected proinflammatory cytokines ligand and receptor pairs. These scores are normalized expression levels; dot size indicates the significance of the interaction, and dot color indicates interaction intensity, calculated by CellPhoneDB. F Kaplan-Meier plot of disease-free survival in the TCGA database, comparing low versus high IFNG-IFNGR1 (left panel) and TNF-TNFRSF1A (right panel) ligand-receptor pair’s mean expression in non-Hepatitis virus HCC patients. The number of patients from each subgroup and the Log-Rank t-test results are indicated in the plot. See also Supplementary Fig. [264]14. Next, we looked for the most prevalent L-R interactions in myeloid, malignant, and T/NK cells. HBV-HCC immunosuppressive interactions included SPP1 with PTGER4/CD44^[265]43, MDK with LRP1/SORL1^[266]44, and MIF-TNFRSE14^[267]45. In non-HBV-HCC cells, enrichment of IFNG-IFNGR1 and LGALS9 interactions with CD44, SORL1, DAG1, CD47, and MET were reported to induce tumor cell apoptosis^[268]46 (Fig. [269]8C). Using the TCGA-LIHC, we found that high MDK-LRP1 expression was associated with poor clinical outcomes in hepatitis virus-mediated HCC, while LGALS9-SORL1 expression was associated with better survival in non-viral HCC (Fig. [270]8D; Supplementary Fig. [271]14B). The immune checkpoints CTLA4-CD86 and PDCD1-FAM3C interacted strongly in HBV and non-HBV-HCC, respectively (Supplementary Fig. [272]14C). Most CD8^+ T cell investigations showed higher CTLA-4 expression in CHB patients. CTLA-4 binds to CD80/CD86 twenty-fold more than CD28, thus inhibiting CD28’s positive co-stimulatory signal and suppressing T cells^[273]47. The TCGA database showed that hepatitis virus-mediated HCC with elevated CTLA4-CD86 expression had worse clinical outcomes (Supplementary Fig. [274]14D). Next, we found that non-HBV T/NK cells had the strongest proinflammatory cytokine interactions with myeloid clusters via the IFNG-IFNGR1, TNF-TNFRSF1B, TNF-VSIR, TNF-TNFRSF1A, TNF-NOTCH1, and TNF-SEMA4C axes, suggesting that IFNG and TNF activated myeloid cells^[275]48 (Fig. [276]8E). A higher expression of proinflammatory cytokine L-R interactions, such as TNF-TNFRSF1B and IFNG-IFNGR1, were linked to favourable survival in non-viral HCC (Fig. [277]8F). In summary, our findings showed that HBV-HCC has more immunosuppressive ligand-receptor interactions. In contrast, non-HBV-HCC had stronger antitumor interactions. Discussion In recent years, the overall survival (OS) rate of hepatocellular carcinoma (HCC) patients has not improved despite advancements in treatments. This is largely due to treatment decisions being made without distinguishing between HBV-induced and non-HBV HCC, even though 50–80% of HCC cases worldwide are caused by HBV infection^[278]49,[279]50. Therefore, it is crucial to differentiate the immune microenvironments of these HCC subtypes to develop more targeted therapeutic strategies. Our study utilized single-cell transcriptome profiling, the Viral-Track method, CD8 cell lineage trajectory analysis, and cell-cell interaction analysis to analyze 70,056 cells from HBV-associated and non-HBV-associated HCC tumors. We employed Viral-Track, an R-based computational pipeline that utilizes the STAR aligner^[280]15, identifies HBV-infected and non-infected cells by mapping reads to both the host genome and a comprehensive viral genome database^[281]10,[282]16, providing valuable insights into the distinct tumor microenvironments of HBV and non-HBV-related HCC. We observed a lower number of myeloid cells, including macrophages and dendritic cells (DCs), in HBV-HCC compared to non-HBV-HCC, using canonical markers such as LYZ, CD68, CD163, C1QA, and C1QB. These findings are consistent with previous studies showing reduced LYZ expression in tumor tissues from the S-I/S-II subtypes of HCC^[283]51. Moreover, prior studies reporting decreased frequencies and functional impairments of plasmacytoid dendritic cells (pDCs) in chronic HBV patients further support our findings of an altered dendritic cell population in HBV-HCC^[284]52,[285]53. Regarding natural killer (NK) cells, our study demonstrated a reduction in NK cell numbers in HBV-HCC, which is consistent with previous reports showing that liver NK cells decrease in number and exhibit functional impairments as HCC progresses. This reduction, coupled with abnormal NK receptor expression, contributes to the progression of chronic HBV infection and HCC, leading to poor prognosis in HCC patients^[286]54,[287]55. We found two predominant NK subtypes (CD16^+ and CD160^+) in non-HBV-HCC. Previously high infiltration of CD160^+ NK cells was found to be associated with improved disease-free survival and overall survival in HCC patients^[288]56. These findings suggest that intratumoral CD160 expression is positively correlated with improved outcomes in HCC patients. Our analysis identified two subtypes of CD4^+FOXP3^+ Tregs (Treg TNFRSF4 and Treg TNFRSF1)-primarily found in tumor-infiltrating Tregs in non-HBV-HCC. Although some studies report a significant increase in tumor-infiltrating FoxP3^+ Tregs in HBV-HCC compared to non-HBV-HCC^[289]57,[290]58. Other studies also suggest that intratumoral Tregs decrease in advanced HBV-HCC compared to early-stage patients, with circulating Treg frequencies increasing as HCC progresses^[291]59. The analysis of copy number variants (CNVs) identified notable chromosomal alterations, including amplifications of chromosomes 1 and 6 in non-HBV-HCC and deletions of chromosomes 5 and 6 in HBV-HCC, consistent with prior studies^[292]20–[293]22. The reproducibility of these CNV alterations across multiple HCC cohorts reinforces their relevance as genomic signatures. Early PCR-based studies on chronic hepatitis B virus infection revealed that the virus may reside in B and T cells in addition to hepatocytes, suggesting lymphoid cells may act as HBV reservoirs^[294]12–[295]14. Since HBV is present in immune cells ∼10 to 100-fold lower than in serum or hepatocytes, its detection is challenging using standard methods^[296]12. Consequently, assessing HBV in liver tissues requires sophisticated molecular methods, such as scRNA sequencing and the Viral-Track method. Previously, the viral track was applied to bronchoalveolar lavage (BAL) samples from moderate and severe COVID-19 patients, revealing the infection landscape of SARS-CoV-2 and its interaction with the host tissue. Similarly, in another study, viral-track combined with scRNA sequencing effectively detected a reduction in HBV transcription in the liver of CHB patients at baseline and 12 and 24 weeks after starting tenofovir alafenamide treatment^[297]11. However, in this study, the Viral-Track method was applied for the first time in combination with scRNA-sequencing data to detect HBV-infected immune cells. Detection of HBV transcripts in immune cell compartments (myeloid, B, and T/NK cells) in liver tumor tissues emphasizes the importance of HBV lymphotropism in viral persistence and genomic diversification, which enables HBV replication to resume later in immunocompromised patients (such as liver transplant recipients) or subjects discontinuing effective antiviral therapy. Although Viral-Track is a novel and potent tool for analyzing host-viral interactions, the potential scarcity of viral readings and infected cells in the sample can limit its applicability. Therefore, alternative approaches are necessary to explore the mechanisms underlying HBV reads and transcripts in immune cells. Multiple HBV transcripts in macrophages and hepatocyte-specific APOB and APOA1 transcripts suggest that macrophages phagocytosed dying infected hepatocytes^[298]11. In HBV-HCC, macrophages displayed noticeable plasticity in phenotypic and functional properties. Notably, we identified a distinct subset of M2 macrophages with high expression of Metallothionein (MT) proteins, particularly MT1G and MT1H, in HBV-HCC TAMs compared to the non-HBV-HCC tumor microenvironment (TME), as supported by validation cohort 2 ([299]GSE149614). MTs are stress-responsive proteins with immune-modulating functions^[300]60. Previously, MT-positive TAMs were reported in scRNA sequencing of human fetal hepatocellular carcinoma^[301]61. MT-positive cell infiltration correlated with the degree of inflammation^[302]62. We also identified MMP-9 in HBV-HCC TAMs, which interact with MTs in tumor development and angiogenesis^[303]24. These findings reinforce the association between MTs and HBV-related HCC, suggesting that MT-enriched TAMs, particularly MT1G^+ TAMs, may serve as a potential immunotherapeutic target. However, further research is needed to elucidate the mechanisms by which MT1G contributes to HBV-mediated tumor progression. The transition trajectories, metabolic pathways, and expression of immune checkpoints and cytotoxic markers in HBV or non-HBV-HCC CD8^+ T cells were not functionally homogeneous due to impairment of effector functions in distinct T cell subsets. HBV-tumor-infiltrated CD8^+ T cells exhibited a cytotoxic dysfunctional T cell phenotype, characterized by upregulated MDK and CTLA4 expression and reduced IFN-γ production, unlike the non-HBV-HCC. While MDK is not typically associated with T cells, its expression has been reported in activated T cells and HIV-infected cells^[304]29. Additionally, studies in prostate cancer have shown that MDK abundance in the tumor microenvironment is linked to reduced cytolytic capacity of CD8^+ T cells, as evidenced by decreased expression of key effector molecules, such as GZMB, PRF1, and IFN-γ^[305]32. Further research is needed to clarify the role of MDK in T cell function and its impact on immune responses in HBV-HCC. IFN-γ produced by T cells plays a key role in eradicating HBV during acute infection by activating macrophages^[306]63. In HBV-HCC, T cells had weaker ligand-receptor interactions with macrophages via IFNG-IFNGR1 due to low IFN-γ levels, which may explain the reduced myeloid cell count in HBV-HCC, potentially contributing to HBV clearance failure. The effectiveness of immune checkpoint blockade (ICB) therapies varies by cancer type and patient characteristics. A better understanding of immune checkpoint expression in the HBV and non-HBV-HCC microenvironments can aid in the development of more effective ICB therapies. Our findings suggest that CTLA4 inhibitors and IFN-γ therapy may be beneficial for treating HBV-HCC. While anti-CTLA4 therapies, including Ipilimumab and nivolumab + Ipilimumab, are already being used in HCC treatment with promising results^[307]64,[308]65. The impact of HBV status on checkpoint blockade response requires further study. Clinical trials, such as CheckMate 459 and Phase 3 HIMALAYA, included HBV-mediated HCC patients with a viral load <100 IU/mL, as higher HBV loads are associated with poorer outcomes in patients receiving immune checkpoint inhibitors^[309]66,[310]67. These findings emphasize that HBV status significantly influences the response to immunotherapy, highlighting the potential benefit of antiviral therapy prior to immunotherapy for HBV^+ HCC patients. Our study found a significantly higher number of malignant hepatocytes in HBV-HCC compared to non-HBV-HCC, consistent with the mechanism of HBV cccDNA integration into hepatocytes, driving their transformation into malignant hepatocytes and progression to HCC^[311]68–[312]71. Moreover, we found that MHC-I antigen presentation was more prominent in non-HBV-HCC tumors, previously reported that HLA class I antigens are strongly expressed in HCC cells, particularly in primary tumors compared to metastatic ones, which is in line with our results that HBV-HCC had higher metastatic potential^[313]72–[314]75. Additionally, these malignant hepatocytes in HBV-HCC displayed increased capacities for metastasis, proliferation, and immune evasion^[315]76–[316]78. In contrast, non-HBV malignant hepatocytes displayed enhanced immune surveillance. In line with the heterogeneous nature of malignant hepatocytes, our data revealed that HBV-HCC tumors exhibited more immunosuppressive ligand–receptor interactions, whereas non-HBV-HCC tumors showed comparatively enriched antitumor interactions. These findings underscore the mechanistic impact of tumor heterogeneity on immune evasion and its implications for immunotherapy targeting. In summary, detecting HBV transcripts using the Viral-Track method offered an opportunity to truly measure host-pathogen interactions and HBV lymphotropism within the human liver. Our data can be a valuable resource, facilitating a deeper understanding of the immune and tumor phenotypes associated with HBV and non-HBV tumors and assisting in developing more effective therapeutic targets and biomarkers for immunotherapies in HCC patients. Methods Tissue samples from patients with HCC Eight surgically resected hepatocellular carcinoma (HCC) tissue samples were obtained and categorized into five HBV-associated HCC cases (KC02A, KC06A, KC07A, KC08A, and KC09A) and three non-HBV HCC cases (KC03A, KC05A, and KC10A). Classification was based on the presence or absence of HBV viral reads in tumor tissues, as detected by the Viral-Track method^[317]10,[318]11. HBV-associated HCCs were selected to have a uniform etiological background of HCC because >80% of our local HCC cases were caused by chronic HBV infection. Each HCC patient was newly diagnosed, pathologically confirmed, and proved not to have other types of cancer. The HCC tumour tissues were obtained after surgical resection at Hubei Cancer Hospital, Wuhan, Hubei, China. This study was approved by the Research and Ethical Committee of Hubei Cancer Hospital (study approval No: LLHBCH2020LW-026). Informed consent was obtained from each patient or their guardian. All ethical regulations relevant to human research participants were followed. The patients’ HBV reads in tissue as detected by the Viral-Track method, along with clinical information, were summarised in Supplementary Data [319]1. Tissue dissociation and single-cell sorting Following surgical resection, tissue dissociation and preparation of single-cell suspensions were performed as previously described^[320]8. The fresh tissues were immediately transferred to pre-cooled MACS Tissue Storage Solution (Cat#130-100-008, Miltenyi Biotech, Germany) and shipped at 4 °C. For each sample, ~1 g tissue was used for dissociation, and the remaining tissue was fixed in formalin for 48 h and embedded in paraffin for subsequent immunohistochemistry analysis. The tissue was dissociated in 5 mL RPMI 1640 medium (Cat#12633020, Thermo Fisher, USA) supplemented with 0.5 mg/ml collagenase II (Cat# C6885-100MG, Sigma-Aldrich, USA) and 0.02 mg/ml DNase I (Cat# 11284932001, Sigma-Aldrich, USA) at 37 °C for 25 min on a shaker. Subsequently, 15 mL RPMI 1640 medium was added to terminate the dissociation process, and the cell pellet was collected by centrifugation after passing through a 70 μm cell strainer (Cat#130-098-462, Miltenyi Biotech, Germany). The cell suspension was treated with 1.5 mL of 1X Red Blood Cell Lysis Solution (Cat#130-094-183, Miltenyi Biotech, Germany) and incubated on ice for 10 min to lyse the red blood cells. The cell pellet was washed with 10 mL of PBS containing 0.04% BSA at 4 °C and centrifuged at 300 g for 5 min. The supernatant was then discarded. The cell pellet was resuspended in 10 mL of PBS containing 0.04% BSA, passed through a 30 μm cell strainer (Cat#130-098-458, Miltenyi Biotech, Germany), and centrifuged at 400 g for 5 min at 4 °C. The supernatant was then discarded. Single-cell RNA reverse transcription, amplification, and sequencing The viable singlet cell suspension was prepared at a concentration range of 700–1200 cells/µl, as recommended by the user protocol of 10× Genomics. It was then subjected to single-cell capture using the Chromium platform ([321]https://www.10xgenomics.com/instruments/chromium-controller), 10× Genomics, CA, USA). The chromium platform is a droplet-based system in which HCC single cells, gel beads with barcoded oligos, and reagents were mixed and captured as droplets in oil emulsion. Gene expression application was performed. Pooled and barcoded single cells were subjected to subsequent manipulation procedures, including complete cell lysis, mRNA capture, reverse transcription, cDNA amplification, and Illumina sequencing library preparation (Cat# 1000121, Illumina, San Diego, USA). For each patient, all the samples (HBV and non-HBV tissues) were processed in parallel in the same thermal cycler. Sequencing libraries of pooled single cells were obtained and sequenced using the Illumina Novaseq 6000 platform (USA) to provide adequate coverage. Public scRNA-seq datasets Beyond the newly sequenced scRNA-seq data, we downloaded viral sequence data of five HBV-HCC samples([322]GSM7441022, [323]GSM7441024, [324]GSM7441027, [325]GSM7441029, and [326]GSM7441031) from the NCBI Gene Expression Omnibus (GEO) database under the accession number [327]GSE234015^[328]11 as Validation cohort 1. Furthermore, we additionally curated a previously published Lu’s study^[329]9 processed gene expression data of single-cell RNA sequencing 5 HBV-positive HCC (HCC01T, HCC02T, HCC06T, HCC09T, HCC10T) and three non-HBV HCC patients (HCC05T, HCC07T, HCC08T) available at the GEO database under the accession code [330]GSE149614, used as a validation cohort 2 to corroborate our findings and results robustness. Immunohistochemical staining in human HCC samples Formalin-fixed, paraffin-embedded (FFPE) tumor sections from eight hepatocellular carcinoma (HCC) patients were used for immunohistochemistry (IHC). The cohort included five HBV-associated HCC patients (KC02A, KC06A, KC07A, KC08A, and KC09A) and three non-HBV-HCC patients (KC03A, KC05A, and KC10A), as detailed in Supplementary Data [331]1. Tissue sections (5μm thick) were baked at 70 °C for 1 h, deparaffinized in xylene, and rehydrated through a graded ethanol series (100%, 95%, and 75%). Antigen retrieval was performed using citrate buffer (pH 6.0) in a 95 °C water bath for 20 min. Endogenous peroxidase activity was blocked using 3% hydrogen peroxide for 15 min, followed by blocking with 10% normal goat serum at room temperature for 30 min. Sections were incubated overnight at 4 °C with the following primary antibodies: CD68 (1:1000, Cat# ab955, Abcam, USA), MT1G (1:200, Cat# CSB-PA17384A0Rb, CUSABIO Technology LLC, USA), CD8 (1:1000, Cat# GB12068, Diagnocine, USA), IFN-γ (1:200, Cat# abs119966, Absin Bioscience Inc, China), CD3 (1:1000, Cat# ab237721, Abcam, USA), and anti-HBC (courtesy of Prof. Yuchen Xia, Wuhan University). Double staining was performed for CD68 and MT1G, CD8 and IFN-γ, and CD3 and HBC in paired HBV and non-HBV HCC samples. Following primary antibody incubation, sections were washed and incubated with appropriate HRP- or fluorescence-conjugated secondary antibodies. Nuclei were counterstained using DAPI. All slides were imaged under 200× magnification using a fluorescence microscope. Single-cell RNA-seq data processing and integration We used Cell Ranger (v3.1.0) from 10× Genomics for initial data demultiplexing, read alignment, UMI counting, and annotation on the raw read data. Cell Ranger filtered the low-quality reads with default parameters. Then, the clean reads were aligned to the human reference genome (GRCh38) to quantify gene expression in single cells. The Seurat R package (v3.1.5)^[332]17 was used to filter, preprocess, scale, and normalize the gene-barcode matrices of transcriptomic data (one for each HCC case). To address potential ambient RNA contamination, we applied the DecontX algorithm^[333]79 after initial quality control and doublet removal using DoubletFinder (v 2.0.3)^[334]80. For batch correction, we employed the SCVI and SCANVI models from the scVI-tools package ([335]https://github.com/scverse/scvi-tools), which were implemented in Python (v3.10.11). SCVI was initialized with 90 latent dimensions and trained for 100 epochs to generate batch-corrected embeddings, which were then used in SCANVI to refine cell-type annotations while preserving batch correction^[336]81. The integrated dataset was visualized using UMAP and stored as an AnnData object for further analysis. Cell annotation and HBV identification in single cells A total of 70,056 cells were selected for further analysis. Cell clustering and markers identification were performed on the single cells using the reads count matrix as input by the Seurat R package (v3.1.5). The Louvain method was used for clustering, and data were visualized using uniform manifold approximation and projection implemented by the Seurat package. We performed principal component analysis (PCA) on the normalized expression matrix using 4000 highly variable genes identified by the “FindVariableFeatures” function. Following the results of PCA, the appropriate principal components (PCs) were selected for clustering with the specific resolution parameters (see below). Finally, to detect cluster-specific expressed genes, the clusters were compared pairwise using the Seurat “FindAllMarkers” function to test genes with >0.25-fold difference (log-scale) on average between the two groups of cells and detectable expression in more than 50% of cells in either of the two populations. For the clustering of all cells, the top 12 PCs were selected with a resolution parameter of 0.5. For the clustering of T/NK lymphocytes, the top 10 PCs were chosen with a resolution parameter of 0.5. For the clustering of myeloid cells, the top 10 PCs were selected with a resolution parameter of 0.5. For the clustering of malignant hepatocyte cells, the top 10 PCs were selected with a 0.5 resolution parameter. The viral read filtering and alignment in scRNA-seq data were performed using the Viral-Track method ([337]https://github.com/PierreBSC/Viral-Track) using the viruSITE database (release 2020.3), containing 12,163 genome sequences from 9297 viruses^[338]16. Viral segments with more than 50 mapped reads were detected using the samtools (v1.16.1). The viral contigs from the Viral-Track analysis were combined for each patient and used as inputs for HBV tracking in different types of single cells. Infected and non-infected cells were visualized using uniform manifold approximation and projection implemented by the Seurat package. Single-cell chromosomal copy number variations (CNV) The inferCNV R package (v1.10.1) ([339]https://github.com/broadinstitute/infercnv) was used to infer the genetic profiles of each cell based on the average gene expression level of large gene sets (100 genes) within a chromosome across all cells to identify malignant cells with clonal large-scale chromosomal copy number. Cells with a deletion at a specific region will show an average lower expression level from that region than cells without the deletion. The posterior probabilities for each cell belonging to one of the components were then used to construct a heatmap that visualizes the copy number status of each cell. Identification of DEGs Differentially expressed gene analysis for HBV and non-HBV myeloid, infected and non-infected myeloid, CD8^+ T Cells, and malignant cells was carried out using the “FindMarkers” function, implemented in the Seurat package, with log-scaled fold change = log [0.25] and P value < 0.05 (Wilcoxon Rank Sum test). Differentially expressed genes passing the criteria were shown in Supplementary Data [340]3, [341]4, [342]6, and [343]8. Cell developmental trajectory of CD8^+ T cells The cell lineage trajectory of CD8+ T was inferred using Monocle2 (v2.30.0)^[344]82. We first used the “newCellDataSet” function in Monocle2 to convert the count’s matrix into normalized mRNA counts and created an object with the parameter “expressionFamily = negbinomial.size”. We used the “differentialGeneTest” function to derive DEG from each cluster, and genes with a q-value < 1e-3 were used to order the cells in pseudotime analysis. After the cell trajectories were constructed, differentially expressed genes along the pseudotime were detected using the “FindAllMarks” of the Seurat R package “differentialGeneTest” function. Definition of cell scores and signature We used the average expression (measured by log[2] (TPM + 1) of 5 metallothionein genes (MT1G, MT1H, MT1E, MT1X, and MT2A) to define the metallothionein score in macrophages. To assign M1/M2 polarization and pro-/anti-inflammatory ability estimates to macrophage cells, we applied Gene Set Variation Analysis (GSVA) (v1.48.3)^[345]83 using standard settings, as implemented in the GSVA package. The gene sets associated with the above functions were described previously (Supplementary Data [346]2). The features of CD8^+ T cells were defined using the mean log[2] (TPM + 1) normalized expression of 10 cytotoxicity-associated genes (GNLY, GZMB, PRF1, GZMK, IFNG, GZMA, NKG7, GZMH, GZMM, CST7) and 8 exhausted markers (CTLA4, TIGIT, LAG3, PDCD1, HAVCR2, TOX, EOMES, CXCL13) to define the cytotoxic and exhausted scores. To assign the metabolism and signalling pathways along the pseudotime, the cluster-specific expressed genes of CD8^+ T cell clusters from HBV/non-HBV HCC tumors were used to investigate the genes enriched in different metabolic and signalling pathways by METASCAPE^[347]84. To define metabolism and signaling pathways phenotypes in CD8^+ T cells in HBV and non-HBV, the signatures of metabolism were determined as the mean expression of gene signatures involved in oxidative phosphorylation, cholesterol metabolism, cysteine and methionine, glutathione, Fatty acid, and carbon metabolism, while the signatures of signaling pathways were determined as the mean expression of gene signatures involved in PD-1 checkpoint, chemokine, HIF-1, PPAR, TNF, VEGF and Ras signalling pathways (Supplementary Data [348]5). The features of malignant cells were defined with the mean log[2] (TPM + 1) normalized expression of metastatic, hepatic, proliferation, immune surveillance, and immune escape genes (Supplementary Data [349]9). Pathway analysis To investigate biological state or function differences, the differentially expressed genes (DEGs) between HBV/non-HBV malignant cells were used to investigate Hallmark gene set enrichment, and the DEGs of HBV HCC macrophages were used to analyze KEGG pathway enrichment analysis using the online METASCAPE tool^[350]84. Only the topmost enriched terms were reported, and the p-value was <0.05. Survival analysis We selected MDK^highCD8^+ T cells (log[2] expression>0.585) and MDK^lowCD8^+ T cells (log[2] expression<0.585), then identified the DEGs between MDK^high and MDK^lowCD8^+ T cells with the “FindMarkers” function in the Seurat R package. The MDK^highCD8^+ T cell feature genes included MDK, CFAP47, RBP2, CRABP2, EFCAB12, TRPC1, FANK1, EPO, RFX6, KCNS3, PTGDR2, SLC30A8, CCDC155, ZIC4, POTEI, EFR3B, ZNF497, SLC2A14, IL5RA, ARHGAP28, NPTX1, SAMD11, DENND2A, SH3RF3, FGF12, RASL11B, GLB1L3, TRPM6, SGCZ, ITGB1BP2, TBX4, CPAMD8, NME5, SYNGR4, LRRC66, PDZD3, FOXB1, CKMT1B, SELP, NFE2, GALR3, FGF20, SLC34A3, CCNB3, APELA, DPEP3, SLC25A2, BAALC, KIF18B, CACNA1I. Kaplan-Meier survival curves of MT1G, KLRF1, KLRD1, KLRB1, MDK, MDK^highCD8^+ T cell signatures (top 50 signatures as described above), and the mean expression of ligands-receptors such as MDK-LRP1, LGALS9-SORL1, IFNG-IFNGR1, TNF-TNFRSF1A and CTLA4-CD86 were plotted in the Cancer Genome Atlas (TCGA) Liver Cancer (LIHC) data using the online tool Kaplan-Meier plotter^[351]85 ([352]https://kmplot.com/analysis). Cell-cell interaction analysis To analyze cell-cell interactions between different cell types, we used CellPhoneDB (V 3.0.0)^[353]86 to identify significant ligand-receptor pairs within HBV and non-HBV HCC samples (Supplementary Data [354]10-[355]11). Cell-type-specific receptor-ligand interactions were found for HBV- and non-HBV-HCC based on cell-type-specific receptor and ligand expression. If a pair of ligands and receptors had an interaction intensity greater than 1 and an adjusted P value less than 0.01 between two cell clusters, we defined this L-R pair as a potential molecular axis mediating cellular interactions. Statistics and reproducibility The statistical approaches and experimental design used in this study are detailed in the Methods and figure legends. All statistical comparisons were conducted using either unpaired two-tailed Wilcoxon rank-sum tests or unpaired two-tailed Student’s t-tests, depending on the data distribution. Analyses were performed using GraphPad Prism 9 and R (version 4.2.2), with statistical significance defined as p < 0.05. Exact p-values are reported in the figure panels or legends where applicable. Single-cell RNA-seq analyses were based on biologically independent tumor samples, comprising five HBV-HCC and three non-HBV-HCC cases in both the discovery and validation cohorts. Replicates were defined as independent clinical samples from different patients. No technical replicates were used for scRNA-seq experiments. Immunohistochemistry (IHC) was performed on independent tumor sections, and representative images are provided. All sample sizes, statistical tests, and reproducibility details are indicated in the corresponding figure legends. Reporting summary Further information on research design is available in the [356]Nature Portfolio Reporting Summary linked to this article. Supplementary information [357]Supplementary Information^ (5.8MB, pdf) [358]Supplementary Data 1.^ (11.3KB, xlsx) [359]Supplementary Data 2.^ (13KB, xlsx) [360]Supplementary Data 3.^ (28.2KB, xlsx) [361]Supplementary Data 4.^ (39.4KB, xlsx) [362]Supplementary Data 5.^ (13.4KB, xlsx) [363]Supplementary Data 6.^ (43.8KB, xlsx) [364]Supplementary Data 7.^ (21.4KB, xlsx) [365]Supplementary Data 8.^ (91.7KB, xlsx) [366]Supplementary Data 9.^ (12.8KB, xlsx) [367]Supplementary Data 10.^ (568.1KB, xlsx) [368]Supplementary Data 11.^ (481KB, xlsx) [369]Supplementary Data 12.^ (14.5KB, xlsx) [370]42003_2025_8443_MOESM14_ESM.pdf^ (37.3KB, pdf) Description of Additional Supplementary Materials [371]Reporting Summary^ (2.4MB, pdf) [372]Transparent Peer Review file^ (1.4MB, pdf) Acknowledgements