Integrated single-cell and transcriptomic analysis of CD8⁺CD101⁺TIM3⁺ T cells in hepatocellular carcinoma: implications for tumor microenvironment and prognostic modeling

Introduction

Hepatocellular carcinoma (HCC) accounts for the third most common contributor to global cancer-induced mortality, with geographic disparities evident in its incidence—70.1% of cases originate in Asia, predominantly in China (60.5%) (1). Despite therapeutic advances including surgical resection and tyrosine kinase inhibitors, the 5-year survival rate stagnates below 20% due to frequent recurrence and chemoresistance in advanced stages (2,3). Thus, this unmet clinical need underscores the urgency for innovative strategies targeting HCC’s intrinsic molecular drivers.

Immunotherapy has emerged as a transformative paradigm in oncology, with immune checkpoint inhibitors demonstrating certain success in HCC (4-6). However, the therapeutic efficacy remains suboptimal, primarily because of the intricate nature of tumor microenvironment (TME) (7). Recent advancements in single-cell RNA sequencing (scRNA-seq) technology have created an unparalleled chance to dissect the traits of different immunocyte subgroups present within the TME. Prior studies suggest that scRNA-seq-derived signatures might predict immunotherapy responses in solid tumors (8,9). Among various immune cell types, CD8⁺ T lymphocytes are key players of antitumor immune surveillance (10-12). However, the tumor-induced immune suppression often facilitates the exhaustion of CD8⁺ T cells, which contributes to the evasion of immune responses and is related to an unfavorable prognosis in HCC (13-15). CCD8⁺CD101⁺TIM3⁺ T cells (CCT) constitute a distinct subpopulation of exhausted CD8⁺ T lymphocytes (16), which are linked to immune dysfunction in various cancers (17-20). The identification and functional understanding of CCT in HCC have opened up new avenues for potential immunotherapy. However, their mechanistic contributions to HCC progression and therapeutic resistance remain unexplored.

Our study seeks to explore the specific phenotypes of immune cell subsets within the HCC’s TME, with a particular emphasis on examining the functional properties of CD8⁺ T cells. Additionally, we aim to identify innovative therapeutic targets that could enhance the efficacy of current treatment modalities and assist in the progression of personalized medicine in HCC management. By addressing these objectives, our research aims to bridge the gap between current treatment limitations and the emerging potential of targeted immunotherapies, ultimately aspires to enhance the outcomes and propel advancements in clinical treatment strategy of HCC. We present this article in accordance with the TRIPOD reporting checklist (available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-1639/rc).

Methods

Acquisition of bulk transcriptome data

This study primarily used datasets from The Cancer Genome Atlas (TCGA) and Gene Expression Omnibus (GEO) databases. We retrieved the comprehensive genome expression data for liver hepatocellular carcinoma (LIHC) by utilizing “TCGAbiolinks (version 2.25.0)” (21). Through the R package GEOquery, expression profiling and clinical data from the GEO were obtained. We focused our analysis on the TCGA-LIHC dataset, which includes 374 LIHC samples and 50 normal tissue samples. The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.

Acquisition and processing of scRNA-seq data

A scRNA-seq datasets was obtained through the GEO database. GSE149614 included eight HCC tissues in this study. Raw data in GSE149614 were imported. The GSE149614 dataset was processed utilizing the Seurat package (version 4.2.0) (22). Initially, the inferior quality cells and genes were eliminated according to several criteria: removal of cells expressing fewer than 200 genes; exclusion of genes detected in less than one cell; cells displaying a gene expression range between 200 and 8,000, and cells with less than 20% of mitochondrial gene expression cells were preserved; retain the unique molecular identifier (UMI) reading is less than 10,000 cells. After standardizing the data, we identified highly variable genes. These genes guided our clustering approach, which used principal component analysis and significant principal components. To manage batch effects from different samples, we applied the harmony method. We then implemented the FindClusters function to cluster. Then the “RunUMAP” function was applied to facilitate unified manifold approximation and projection (UMAP). To characterize distinct cellular populations, differentially expressed genes (DEGs) across identified clusters were determined by employing the FindAllMarkers function. Eventually, cell clusters were identified by distinct biomarkers associated with cell types, followed by the calculation and evaluation of the proportions of these cell types.

Single-cell trajectories construction in pseudotime

Pseudo-time analysis was employed through Monocle 2 (23). Genes exhibiting high dispersion and elevated expression levels were selected for the establishment of the pseudotime trajectory (24). Default parameters were applied in accordance with the DDRTree algorithm. Furthermore, we employed the Branched Expression Analysis Modeling feature available in Monocle 2, which facilitates the identification of genes exhibiting significant expression changes dependent on branching events (23). The visualization of branch-dependent expression patterns was accomplished through the generation of a heatmap using Monocle 2.

Analysis of cell communication and expression of ligand-receptor interactions

The “CellChat” R software package generates CellChat objects from the UMI count matrix for both HCC and normal groups (https://www.github.com/sqjin/CellChat) (25). By referencing the “CellChatDB. Human” database, which specializes in ligand/receptor interactions, and using preset parameters, we conducted an examination of intercellular communication. To highlight differences among various cell types within each group, we utilized the “netVisual_diffInteraction” function. With computational visualization tools such as “netVisual_bubble” and “netVisual_aggregate”, we systematically characterized the spatial expression patterns of signaling molecules across diverse cellular subpopulations.

Functional enrichment

DEGs were subsequently carried into Gene Ontology (GO) functional analysis (26) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis (27). The “clusterProfiler (version 4.2.2)” (27) aided in executing GO and KEGG analysis for DEGs in HCC.

Model construction and validation

Initially, univariate Cox regression was used to filter CCT-associated DEGs correlated with overall survival (OS). Clinical data from all tumor samples were randomly split into a training set (n=237) and a validation set (n=101) in a 7:3 ratio. We then applied the LASSO Cox regression model via the “glmnet” package (28) to refine the candidate gene list and construct a model. The prognostic risk score was determined using a predefined mathematical formula:

Riskscore=∑i=1nCoef(genei)∗Expression(genei)[1]

An optimal cutoff value was employed to categorize patients in the training cohort into distinct prognostic groups: low-risk and high-risk. Survival analysis was conducted using the Kaplan-Meier method, with differences evaluated through the log-rank test. Additionally, the validation group served to assess the model’s efficacy using the receiver operating characteristic (ROC) curve.

Nomogram development and validation

Clinical data and pertinent variables were obtained from the TCGA. Both univariate and multivariate Cox regression analyses were conducted to assess prognostic significance. Furthermore, a comprehensive nomogram was devised using the “RMS” R package, facilitating individualized survival probability estimation. Its predictive performance was subsequently validated through an ROC curve

Gene set enrichment analysis (GSEA)

GSEA serves as a method to evaluate significant variations in predefined gene sets across distinct biological conditions (29). In this study, we conducted a differential expression analysis comparing short-term and long-term peritoneal dialysis cohorts using “limma” (version 3.50.0) (30). We calculated gene expression fold changes (FCs) between the two groups and performed GSEA using “clusterProfiler” (version 4.2.2), ranking genes based on their log₂ FC values. The gene set utilized was sourced from the Molecular Signatures Database (31,32).

Gene set variation analysis (GSVA)

To explore functional differences between the HCC and control groups, we implemented GSVA using “GSVA” (version 1.42.0) and visualized the resulting data with “pheatmap” (version 1.0.12).

Immune infiltration analysis

We sourced 28 distinct immunocyte types from the Tumor and Immune System Interactions Database (33) and computed the enrichment score for each immune cell in HCC patients using single-sample GSEA (ssGSEA) (34). Discrepancies in immune cell infiltration levels were visualized using “ggplot2” (version 3.3.6) (35).

Evaluation of the drug susceptibility

Drug response prediction for HCC samples was carried out using “oncoPredict” (version 0.2), which analyzed half maximal inhibitory concentration (IC50) values from the Genomics of Drug Sensitivity in Cancer platform (36,37). This approach facilitated a systematic assessment of potential therapeutic efficacy against HCC based on pharmacological profiling data.

Somatic mutation analysis

The genomic alteration profile was comprehensively characterized through the analysis of mutation and copy number variation data. The “maftools” package aided in the visualization of somatic mutations across different molecular subgroups (38). Genes exhibiting the highest mutation frequency, known as frequently mutated genes, were identified as key oncogenic drivers in tumorigenesis (39).

Prediction of the immunotherapy responses

To predict patients’ immunotherapy responses, the Tumor Immune Dysfunction and Exclusion (TIDE) analysis was conducted in this study (40).

Immune checkpoint-related genes expression

Furthermore, immune regulatory molecules expressed on the surface of immunocytes play a critical role in determining the intensity of immune activation. These checkpoints are essential for preventing excessive immune responses. Our study focused on the expression of genes related to the immune checkpoint response.

Statistical analysis

Statistical analysis was performed using R software. We used “survival” and “survminer” for Kaplan-Meier survival analysis. Correlation between continuous variable were evaluated through Spearman correlation test. Nonparametric Kruskal-Wallis and Wilcoxon tests were used for continuous variables. Results were considered significant when P value was below 0.05.

Results

Figure 1 shows flow chart of this research.

Figure 1 The flow chart for this research. GEO, Gene Expression Omnibus; GSEA, Gene Set Enrichment Analysis; GSVA, Gene Set Variation Analysis; LASSO, least absolute shrinkage and selection operator; PCA, principal component analysis; ROC, receiver operating characteristic; scRNA-seq, single-cell RNA sequencing; TCGA, The Cancer Genome Atlas; TIDE, Tumor Immune Dysfunction and Exclusion; TMB, tumor mutation burden; UMAP, unified manifold approximation and projection.

Single-cell sequencing analysis

Utilizing the single-cell sequencing dataset GSE149614, we conducted an analysis of the HCC cell population to investigate the origin of highly expressed genes. Following a preliminary evaluation and quality control, 25,486 cells were obtained. Figure 2A demonstrates the clustering of these cells into 24 distinct clusters. Cell-specific biomarkers, determined by examining gene expression profiles characteristic of each cluster, were used to annotate the various cell types. Additionally, further detailed annotation led to the classification of seven primary cell types: T/natural killer (NK) cells, fibroblasts, myeloid cells, B cells, and endothelial cells (Figure 2B). Cumulative histograms depict the proportions of these cell types among the groups (Figure 2C), while a heat map illustrates gene expression characteristics in each cell group (Figure 2D).

Figure 2 Identification of cell subsets. (A) UMAP shows the distribution of HCC cellular subsets. (B) UMAP figure shows the annotation of HCC cells. (C) The cumulative histogram shows the distribution of cell types in each sample. (D) Heat map shows the cell group of specific gene expression. HCC, hepatocellular carcinoma; NK, natural killer; UMAP, unified manifold approximation and projection.

T cell subsets were reclustered

In our study, we focused on assessing the impact of T cells on immunotherapy efficacy. Initially, we isolated T cells within HCC and then reclustered these cells (Figure 3A). As shown in Figure 3B, our analysis revealed that T cells could be categorized into four distinct types: CD8⁺ T cells, CCT, CD8⁺CD101⁻TIM3⁺ T cells, and CD4⁺ T cells. Each cell type exhibits unique gene expression profiles, visualized through a spot diagram (Figure 3C). Additionally, we conducted a GSVA enrichment analysis using the MsigDB database [online table 1 (available at https://cdn.amegroups.cn/static/public/tcr-2025-1639-1.csv)]. The most significantly different pathway among the T cell subsets was selected to create a heat map representing pathway activity (Figure 3D).

Figure 3 T cell subsets. (A) The distribution and (B) the annotation results of T cell subpopulations in UMAP plots. (C) Marker genes expression in distinct cellular populations. (D) The T cells subgroup of the most significant enrichment pathways in the heat map. Avg exp, average expression; KEGG, Kyoto Encyclopedia of Genes and Genomes; NES, normalized enrichment score; Pct exp, percentage of cells expressing; UMAP, unified manifold approximation and projection.

Pseudotime trajectory analysis

A quasichronological cell trajectory was established to identify the fundamental gene expression program determining HCC progression. This was achieved by utilizing identified CD8⁺CD101⁻TIM3⁺ T cells and CCT. CD8⁺CD101⁻TIM3⁺ T cells (cell types in pseudotime =0) were positioned at the beginning of the trajectory, while CCT was placed at the trajectory’s end. CD8⁺CD101⁻TIM3⁺ T cells served as CCT precursors, and the results revealed the differentiation process from CCT precursors to CCT (Figure 4).

Figure 4 Transcription trajectory analysis reveals the transcription pattern of CCT. (A) Pseudotime analysis visualized by a blue gradient (dark to light). (B) The states of the proposed temporal trajectory divided by Monocle 2. (C) Pseudotime trajectories showing cells distribution based on cell type. CCT, CD8⁺CD101⁺TIM3⁺ T cells.

Cellular communication patterns in HCC microenvironment

To further elucidate the combined effects between T cells, particularly those related to CD8⁺CD101⁻TIM3⁺T, we conducted intercellular communication analyses to pinpoint the key communications that shape CCT differentiation. Figure 5A,5B illustrate the interactions among cell types. We then delved deeper into the signal pathways and specific receptor pairs with CD8⁺CD101⁻TIM3⁺ T as communication receivers, noting significant changes between the two groups. Figure 5C,5D highlight that CD8⁺ T cells and CCT are the primary signal providers and receptors, engaging potential signal pathways including SPP1, MIF, MK, VISFATIN, GALECTIN, and PARs. Our subsequent investigation into specific signal pairs within these cell types revealed that the MIF-CD74/CXCR4 and SPP1-CD44 signaling pathways are the primary signals received by CD8⁺CD101⁻TIM3⁺ T cells. This suggests that the MIF and SPP1 pathways might play a pivotal role in the differentiation of CD8⁺CD101⁻TIM3⁺T cells (Figure 5E). These preliminary findings underscore the dynamic interplay among these cell populations, paving the way for further exploration of CCT’s comprehensive function in HCC. Figure 5F,5G demonstrate that the ligand MIF, receptors CD74 and CXCR4 are abundantly expressed in nearly all cells, while SPP1 is predominantly expressed in CCT and myeloid cells. However, CD44 is present in cells other than endothelial and hepatocyte cells (Figure 5F,5G).

Figure 5 Number and intensity of communication between normal and HCC cells. Network diagrams show the quantity (A) and strength (B) of interactions occurring among cells. (C) Heat maps visualize possible afferent signaling pathways between cells. (D) Heat maps visualize possible efferent signaling pathways between cells. (E) String diagram that visualizes pairs of possible incoming or outgoing signals. (F) Distribution of MIF signal ligand receptors. (G) Expression distribution of SPP1 signal ligand receptor. HCC, hepatocellular carcinoma.

Enrichment analysis of DEGs associated with CCT in HCC

By comparing CCT with other cell types, we identified 1,281 DEGs exhibiting statistically significant variations [|log₂ FC| >0.25, Figure 6A, online table 2 (available at https://cdn.amegroups.cn/static/public/tcr-2025-1639-2.xls)]. To unravel the biological roles associated with these genes, we conducted GO [online table 3 (available at https://cdn.amegroups.cn/static/public/tcr-2025-1639-3.csv)] and KEGG [online table 4 (available at https://cdn.amegroups.cn/static/public/tcr-2025-1639-4.csv)] pathway enrichment analyses. The results indicate that these genes are enriched in biological processes such as “aerobic respiration” and “electron transport chain”, cell components like “respirasome” and “mitochondrial respirasome”, and molecular functions including “electron transfer activity” (Figure 6B). The KEGG pathway enrichment includes “oxidative phosphorylation”, “non-alcoholic fatty liver disease”, and others (Figure 6C).

Figure 6 Enrichment analysis of DEGs associated with CCT in HCC. (A) Heat maps depict DEGs in HCC. (B) GO enrichment analysis of DEGs. (C) KEGG enrichment analysis of DEGs. BP, biological process; CC, cellular component; CCT, CD8⁺CD101⁺TIM3⁺ T cells; DEGs, differentially expressed genes; GO, Gene Ontology; HCC, hepatocellular carcinoma; KEGG, Kyoto Encyclopedia of Genes and Genomes; MF, molecular function.

Construction and verification of prognostic risk model

We further explored key genetic targets potentially influencing patient outcomes. Following univariate Cox analysis, 188 genes tied to the prognostic outcome of LIHC were identified [online table 5 (available at https://cdn.amegroups.cn/static/public/tcr-2025-1639-5.csv)]. We randomly assigned 7 out of 10 LIHC samples (n=338) to a training group (n=237) and the remainder to a validation group (n=101). As illustrated in Figure 7A,7B, 26 genes were detected as prognostic indicators for LIHC [online table 6 (available at https://cdn.amegroups.cn/static/public/tcr-2025-1639-6.csv)]. Based on an optimal cut-off value, samples were categorized into high- and low-risk groups. Kaplan-Meier analysis revealed a more favorable prognosis for the low-risk samples (Figure 7C,7D). In the training cohort, the model exhibited strong predictive power with AUC values of 0.825, 0.873, and 0.865 for 1-, 3-, and 5-year survival predictions, respectively (Figure 7E). In the validation set, the AUC values were 0.654, 0.636, and 0.644, respectively (Figure 7F).

Figure 7 Construction and verification of the prognostic model for LIHC. (A) Change pattern of the predictor variable regressed by LASSO. (B) The confirmation interval for LASSO regression under each lambda. Survival curves of training set (C) and validation set (D). ROC curves of the training set (E) and validation set (F). AUC, area under the curve; LASSO, least absolute shrinkage and selection operator; LIHC, liver hepatocellular carcinoma; ROC, receiver operating characteristic.

By stratifying various clinical characteristics, we also validated the stability of our prognostic model, revealing that age, gender, stage, and T stage significantly impact prognosis (Figure 8). Furthermore, we explored the influence of prognosis-related gene expression levels on outcomes, finding that all evaluated genes significantly affect prognosis (Figure 9).

Figure 8 Influence of clinical characteristics on prognosis. Survival difference between the two risk subgroups for age >65 years (A) and age ≤65 years (B). Survival differences between the two risk subgroups among males (C) and females (D). Survival differences between the two risk subgroups in stage I–II (E) and stage III–IV (F) . Survival differences between the two risk subgroups in T stage 1–2 (G) and T stage 3–4 (H).

Figure 9 The impact of prognosis-related gene expression levels on survival. ANXA11 (A), YIF1A (B), UBE2E3 (C), TPP1 (D), TPM1 (E), TM4SF1 (F), SNRPE (G), SMS (H), SMIM14 (I), SERINC2 (J), S100A16 (K), S100A13 (L), RTN3 (M), PSMA7 (N), PON1 (O), POLR2L (P), IGFBP3 (Q), HMGCS2 (R), FEZ2 (S), CYB5R3 (T), CXCR3 (U), C3 (V), PLPP1 (W), PHGDH (X), METTL9 (Y), and MAN1A1 (Z) significantly affect prognosis.

Development and verification of nomogram

The Cox regression analyses demonstrated that the risk score serves as a statistically significant independent predictor of prognosis (Figure 10A,10B). Based on these findings, we constructed a comprehensive nomogram that integrates the risk score with other clinical variables. This nomogram exhibited robust performance in predicting patient survival outcomes (Figure 10C). The predictive accuracy of this model was further validated through ROC analysis, with AUC values of 0.796, 0.818, and 0.796 for 1-, 3-, and 5-year survival, respectively, indicating high predictive accuracy (Figure 10D).

Figure 10 Development and verification of a LIHC diagnostic nomogram. (A) Forest plots illustrate the findings derived from univariate Cox regression analyses. (B) The clinical characteristics of multivariate Cox regression analysis. (C) Nomogram of prediction model. (D) The column chart displays the dependence of the ROC curve at 1, 3, and 5 years. AUC, area under the curve; CI, confidence interval; HR, hazard ratio; LIHC, liver hepatocellular carcinoma; M, metastasis; N, node; ROC, receiver operating characteristic; T, tumor.

GSEA and GSVA

To gain deeper insights into the mechanisms of DEGs, we conducted GSEA using data from MsigDB. The most significant pathways were identified [online table 7 (available at https://cdn.amegroups.cn/static/public/tcr-2025-1639-7.csv)]. GSEA results showed enrichment in pathways such as ribosome, spliceosome, and DNA replication, among others (Figure 11A-11F). Additionally, GSVA revealed upregulation of certain pathways in the low-risk group and others in the high-risk group [Figure 11G, online table 8 (available at https://cdn.amegroups.cn/static/public/tcr-2025-1639-8.csv)].

Figure 11 GSEA and GSVA. Ribosome (A), spliceosome (B), DNA replication (C), linoleic acid metabolism (D), glycine serine and threonine metabolism (E), and fatty acid metabolism (F) were significantly enriched. (G) GSVA of differentially enriched pathways between the two risk subgroups. FDR, false discovery score; GSEA, gene set enrichment analysis; GSVA, gene set variation analysis; KEGG, Kyoto Encyclopedia of Genes and Genomes; NES, normalized enrichment score.

Immune infiltration

An analysis of the infiltration levels in 28 distinct immunocyte subsets utilizing ssGSEA methods showed statistically significant variations in the infiltration densities of multiple immunocyte populations, including gamma-delta T cell, activated dendritic cell, activated CD8⁺ T cell, activated CD4⁺ T cell and other immunocytes in the two risk subgroups (P<0.05, Figure 12A). Moreover, most immunocyte populations exhibited significant positive interrelationships, although a limited subset demonstrated inverse associations; specifically, activated CD8⁺ T cell and CD56dim natural killer cells (Figure 12B). The correlation analysis revealed that the risk score showed: positive relationships with central memory CD8⁺ T cell (R=0.1942), activated dendritic cells (R=0.2792), regulatory T (Treg) cell (R=0.207), and macrophages (R=0.2302); inverse relationships with effector memory CD8⁺ T cell (R=0.3605), natural killer cell (R=−0.2605), activated CD8⁺ T cell (R=−0.2352), type 1 T helper cell (R=0.3502) and activated B cell (R=−0.2378) (Figure 13).

Figure 12 Immune infiltration level in the two groups. (A) Boxplot shows infiltration level of immunocytes. (B) The correlations between immunocytes.

Figure 13 The association between risk scores and immunocytes. The risk score showed relationships with activated dendritic cells (A), central memory CD8⁺ T cell (B), regulatory T cell (C), macrophages (D), type 1 T helper cell (E), effector memory CD8⁺ T cell (F), natural killer cell (G), activated CD8⁺ T cell (H), and activated B cell (I). TPM, transcripts per million.

Tumor mutation burden (TMB) and drug susceptibility analysis

Comprehensive mutational analysis in LIHC identified the 20 most frequently altered genes, with TP53 having the highest mutation prevalence, followed by that of CTNNB1 (Figure 14A,14B). TMB represents a pivotal biomarker for predicting response to immunotherapy. Consequently, an analysis of somatic mutations linked to HCC was performed, revealing no noteworthy differences in TMB levels in two cohorts (Figure 14C). Additionally, the risk score demonstrated predictive performance in chemotherapy sensitivity, indicating potential enhanced therapeutic sensitivity in the high-risk group to certain drugs [Figure 14D-14L, online table 9 (available at https://cdn.amegroups.cn/static/public/tcr-2025-1639-9.csv)].

Figure 14 The analysis of TMB and drug sensitivity. (A,B) The highly mutated genes are depicted in the two groups. (C) The TMB difference between the two groups. (D-L) The difference of the drug sensitivity between the two groups. TMB, tumor mutation burden.

Immune checkpoints gene expression and immunotherapy response

Significant differences in immune checkpoint gene expression were observed between risk levels, including TNFRSF4, CD274, LAG3, SIGLEC10, HAVCR2, and CTLA4 (Figure 15A). TIDE analysis revealed that lower TIDE scores in the low-risk group correlate with increased responsiveness to immunotherapy (Figure 15B,15C).

Figure 15 Immune checkpoints and TIDE analysis. (A) Boxplots showing the difference in immune checkpoint levels. (B) The difference in T cell dysfunction, exclusion and TIDE score. (C) Sankey plot shows the response to immunotherapy predicted by TIDE analysis. *, P<0.05; **, P<0.01; ***, P<0.001; ****, P<0.0001. TIDE, Tumor Immune Dysfunction and Exclusion.

Discussion

HCC is one of leading cause of cancer-related mortality worldwide. Despite remarkable progress on surgical techniques and adjuvant therapy, the prognosis of patients with HCC remains poor (1). Immunotherapy has demonstrated certain success in HCC. However, the intricate nature of TME reduces its anticancer efficacy (5). Hence, it is imperative to elucidate the mechanisms behind immunotherapy efficacy and develop more effective strategies to combat cancer. In our study, T cells were further subclassified into CD4⁺ T cells, CD8⁺ T cells, CD8⁺CD101⁻TIM3⁺ T cells (precursor cells) and CCT in HCC. This comprehensive cellular profiling highlights the diverse array of cell types present within HCC tissues. Our results are consistent with prior research highlighting the pivotal function of CD8⁺ T cells in antitumor immunity (41). These cells are crucial for establishing effective antitumor responses and are linked to improved outcomes in HCC patients (42). However, immunosuppressive cells, including tumor-associated macrophages, could trigger CD8⁺ T-cell exhaustion, thereby restricting the overall antitumor effect and contributing to tumor advancement (43-45). This detailed cellular mapping underscores the diversity of cell populations within HCC tissues, reflecting their potential impact on tumor progression and response to treatment.

The pseudotime analysis revealed a seamless progression of transcriptional states from CD8⁺CD101⁻TIM3⁺ T cells to CCT, suggesting a sequential differentiation pathway that might play a pivotal role in the development of HCC. This finding aligns with previous studies that have demonstrated distinct functions of specific CD8⁺ T cell subsets in tumor growth and immunotherapy response (18,19,46). Notably, the irreversible expression of TIM3 and CD101 marks a key feature of CD8⁺ T cell exhaustion in cancer (47). In the context of cancer, CCT cells signify terminal differentiation and might significantly contribute to suppressing antitumor T cell activities within the TME, leading to persistent immunosuppression (48,49). Furthermore, analysis of cell-to-cell communication indicated that CD8⁺CD101⁻TIM3⁺ T cells predominantly receive signals from the MIF-CD74/CXCR4 and SPP1-CD44 pathways, primarily originating from CCT. This reflects distinct functional differences between the two cell types. Previous research has established both MIF and SPP1 involved in regulating tumor growth and metastasis, including HCC (50-54). By synthesizing our findings with this prior knowledge, we postulate that the MIF and SPP1 signaling pathways could serve as regulators of CD8⁺CD101⁻TIM3⁺ T cell differentiation. Consequently, targeting TIM3, CD101 and these pathways in HCC might revitalize these exhausted T cells, potentially boosting antitumor immunity

Differential expression analysis revealed 1,281 CCT-specific genes that are significantly enriched in pathways pertaining to electron transport chains, respiratory chain complexes, mitochondrial damage, and oxidative stress. The electron transport chain and respiratory chain complexes generate energy via oxidative phosphorylation within mitochondria (55). Mitochondrial function profoundly affects the activity and persistence of T cells by modulating the oxidative stress response and metabolic activity, thereby influencing T cell depletion, which is central to immunotherapy (56-59). By regulating mitochondrial function, especially through metabolic intervention or oxidative stress control, the efficacy of immunotherapy could be potentially enhanced. These pathway enrichment findings suggest that CCT-specific genes might mediate mitochondrial function and impact the outcomes of immunotherapy.

We further investigated the functional disparities between annotated CCT and its precursors, revealing that CCT precursors were predominantly involved in linoleic acid metabolism, retinol metabolism, and steroid hormone biosynthesis pathways. Conversely, CCT exhibited notable enrichment in extracellular matrix (ECM) receptor interactions and tyrosine metabolism pathways. Notably, glucocorticoid signaling modulates the differentiation and dysfunction of CD8⁺ T cells within the TME, whereas linoleic acid enhances the metabolic capacity of CD8⁺ T cells, prevents T cell exhaustion, and promotes their differentiation into a memory-like phenotype with potent effector functions. The enrichment in linoleic acid metabolism pathways in precursor cells suggests that dietary fatty acid intervention may block the differentiation of CCT (60). In contrast, the ECM and tyrosine metabolism pathways enriched in CCT are associated with immune evasion, suggesting a link between CCT and unfavorable immunotherapy outcomes (61,62), which provides a theoretical basis for the combination of anti-fibrotic drugs (such as LOXL2 inhibitors) and immunotherapy.

The prognostic model we constructed, incorporating 26 genes, exhibited strong predictive capabilities across diverse clinical cohorts, highlighting the potential of these genes in guiding individualized immunotherapy regimens. Among these 26 signature-associated genes (ANXA11, YIF1A, UBE2E3, TPP1, TPM1, TM4SF1, SNRPE, SMS, SMIM14, SERINC2, S100A16, S100A13, RTN3, PSMA7, PON1, POLR2L, IGFBP3, HMGCS2, FEZ2, CYB5R3, CXCR3, C3, PLPP1, PHGDH, METTL9, MAN1A1), CXCR3 emerged as a key player, with elevated expression on Treg cells in various human cancers. This suggests that CXCR3⁺ Treg cells could represent a promising therapeutic target for immunotherapy (63-66). A study have shown that the absence of CXCR3 enhances CD8⁺ T cell activity, ultimately inhibiting cancer progression (67). Additionally, CXCR3⁺ B cells abundance is positively associated with early tumor recurrence in HCC (68). However, the precise role of CXCR3 in HCC remains unclear. It might increase protumorigenic macrophage polarization and CD8⁺ T cell priming and reactivation in tumor environment, promoting tumor growth.

In addition, the aberrant SNRPE gene expression has been linked to various malignancies and might serve as a predictive indicator for poor OS in HCC (69,70). While its specific functions in HCC are still unclear, a study by Li et al. have shown that SNRPE mRNA expression levels positively correlate with the degree of immunocyte infiltration, such as CD4⁺ T cells, B cells, and CD8⁺ T cells (71). SNRPE might promote tumor growth by influencing the infiltration, activation, and immunosuppressive effects of these immunocytes. In particular, within the immune evasion mechanism of HCC, SNRPE could regulate key pathways in the immune response, thereby affecting immunocyte function and the establishment of the immune microenvironment. Our study also indicated an inverse correlation between the expression levels of most signature-associated genes and the infiltration of immunocytes. Immune checkpoint analysis further revealed significant differences in the expression levels of CD274 and CTLA4 between two cohorts. Overall, among the 26 gene markers, genes such as CXCR3 and SNRPE influenced the TME by altering Treg activity and immune cell infiltration patterns. The elevated CD274/CTLA4 expression and higher TIDE scores observed in the high-risk group align with these regulatory effects. These results indicated that such genes may serve as key mediators of immunotherapy response and represent viable therapeutic targets in HCC.

Despite the findings in our study, there are still limitations and questions that warrant further exploration. Firstly, the research relied solely on bioinformatics approaches. To deepen our understanding of CCT’s biological role within the TME of HCC, future studies should incorporate data from in vitro functional assays. Additionally, although the cohort size was drawn from respected databases such as TCGA and GEO, it remains limited and was not stratified by key molecular subtypes, such as tumors with TP53 or CTNNB1 mutations. This may restrict the generalizability of our findings across distinct HCC subgroups, particularly given the known clinical relevance of these mutations in therapy response and prognosis. Lastly, the prognostic model has not undergone clinical validation, which is essential for translating these results into clinical practice.

Conclusions

In conclusion, this study comprehensively investigated CCT and their precursor cells in the context of HCC, emphasizing their altered functions and correlations with immunotherapy responses. A prognostic model, constructed using 26 genes, demonstrated strong prognostic power in multiple patient populations. Furthermore, the established nomogram offers a practical tool for predicting individualized survival. Our analysis indicated heightened immunotherapy sensitivity among patients categorized as low-risk. These findings provide valuable guidance for personalized immunotherapeutic strategies for HCC patients.

Acknowledgments

We thank all members of our research team for their wholehearted cooperation and their excellent work.

Footnote

Reporting Checklist: The authors have completed the TRIPOD reporting checklist. Available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-1639/rc

Peer Review File: Available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-1639/prf

Funding: This study was supported the National Natural Science Foundation of China (No. 82260627), the Health Commission of Guangxi Zhuang autonomous region (No. ZK20221817) and Bethune Charitable Foundation (No. BJ-YWZJ202211-09).

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-1639/coif). The authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Bray F, Laversanne M, Sung H, et al. Global cancer statistics 2022: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin 2024;74:229-263. [Crossref] [PubMed]
2. Zeng H, Chen W, Zheng R, et al. Changing cancer survival in China during 2003-15: a pooled analysis of 17 population-based cancer registries. Lancet Glob Health 2018;6:e555-67. [Crossref] [PubMed]
3. Vogel A, Meyer T, Sapisochin G, et al. Hepatocellular carcinoma. Lancet 2022;400:1345-62. [Crossref] [PubMed]
4. Chami P, Diab Y, Khalil DN, et al. Radiation and Immune Checkpoint Inhibitors: Combination Therapy for Treatment of Hepatocellular Carcinoma. Int J Mol Sci 2023;24:16773. [Crossref] [PubMed]
5. Langhans B, Nischalke HD, Krämer B, et al. Role of regulatory T cells and checkpoint inhibition in hepatocellular carcinoma. Cancer Immunol Immunother 2019;68:2055-66. [Crossref] [PubMed]
6. Zhu L, Liu J, Chen J, et al. The developing landscape of combinatorial therapies of immune checkpoint blockade with DNA damage repair inhibitors for the treatment of breast and ovarian cancers. J Hematol Oncol 2021;14:206. [Crossref] [PubMed]
7. Chen H, Ye F, Guo G. Revolutionizing immunology with single-cell RNA sequencing. Cell Mol Immunol 2019;16:242-9. [Crossref] [PubMed]
8. Li Y, Zhao X, Liu Q, et al. Bioinformatics reveal macrophages marker genes signature in breast cancer to predict prognosis. Ann Med 2021;53:1019-31. [Crossref] [PubMed]
9. Liang L, Yu J, Li J, et al. Integration of scRNA-Seq and Bulk RNA-Seq to Analyse the Heterogeneity of Ovarian Cancer Immune Cells and Establish a Molecular Risk Model. Front Oncol 2021;11:711020. [Crossref] [PubMed]
10. Mumm JB, Emmerich J, Zhang X, et al. IL-10 elicits IFNγ-dependent tumor immune surveillance. Cancer Cell 2011;20:781-96. [Crossref] [PubMed]
11. Zhang N, Bevan MJ. CD8(+) T cells: foot soldiers of the immune system. Immunity 2011;35:161-8. [Crossref] [PubMed]
12. Hu J, Zhu S, Xia X, et al. CD8+T cell-specific induction of NKG2D receptor by doxorubicin plus interleukin-12 and its contribution to CD8+T cell accumulation in tumors. Mol Cancer 2014;13:34. [Crossref] [PubMed]
13. Yu F, Feng Y, Wang Q, et al. N6-methyladenosine (m6A) Writer WTAP Potentiates Hepatocellular Carcinoma Immune Evasion and Aerobic Glycolysis. Cell Biochem Biophys 2024;82:2321-2331. [Crossref] [PubMed]
14. Gabrielson A, Wu Y, Wang H, et al. Intratumoral CD3 and CD8 T-cell Densities Associated with Relapse-Free Survival in HCC. Cancer Immunol Res 2016;4:419-30. [Crossref] [PubMed]
15. Zhou N, Guo C, Li X, et al. USP24 promotes hepatocellular carcinoma tumorigenesis through deubiquitinating and stabilizing TRAF2. Biochem Pharmacol 2024;229:116473. [Crossref] [PubMed]
16. Onyshchenko K, Luo R, Rao X, et al. Hypofractionated radiotherapy combined with lenalidomide improves systemic antitumor activity in mouse solid tumor models. Theranostics 2024;14:2573-88. [Crossref] [PubMed]
17. Yang G, Cai S, Hu M, et al. Functional status and spatial architecture of tumor-infiltrating CD8+ T cells are associated with lymph node metastases in non-small cell lung cancer. J Transl Med 2023;21:320. [Crossref] [PubMed]
18. Chen C, Zhao F, Peng J, et al. Soluble Tim-3 serves as a tumor prognostic marker and therapeutic target for CD8 T cell exhaustion and anti-PD-1 resistance. Cell Rep Med 2024;5:101686. [Crossref] [PubMed]
19. Hudson WH, Gensheimer J, Hashimoto M, et al. Proliferating Transitory T Cells with an Effector-like Transcriptional Signature Emerge from PD-1(+) Stem-like CD8(+) T Cells during Chronic Infection. Immunity 2019;51:1043-1058.e4. [Crossref] [PubMed]
20. Rosenberg SA, Yang JC, Sherry RM, et al. Durable complete responses in heavily pretreated patients with metastatic melanoma using T-cell transfer immunotherapy. Clin Cancer Res 2011;17:4550-7. [Crossref] [PubMed]
21. Colaprico A, Silva TC, Olsen C, et al. TCGAbiolinks: an R/Bioconductor package for integrative analysis of TCGA data. Nucleic Acids Res 2016;44:e71. [Crossref] [PubMed]
22. Butler A, Hoffman P, Smibert P, et al. Integrating single-cell transcriptomic data across different conditions, technologies, and species. Nat Biotechnol 2018;36:411-20. [Crossref] [PubMed]
23. Qiu X, Mao Q, Tang Y, et al. Reversed graph embedding resolves complex single-cell trajectories. Nat Methods 2017;14:979-82. [Crossref] [PubMed]
24. Karmaus PWF, Chen X, Lim SA, et al. Metabolic heterogeneity underlies reciprocal fates of T(H)17 cell stemness and plasticity. Nature 2019;565:101-5. [Crossref] [PubMed]
25. Fang Z, Tian Y, Sui C, et al. Single-Cell Transcriptomics of Proliferative Phase Endometrium: Systems Analysis of Cell-Cell Communication Network Using CellChat. Front Cell Dev Biol 2022;10:919731. [Crossref] [PubMed]
26. Gene Ontology Consortium. going forward. Nucleic Acids Res 2015;43:D1049-56. [Crossref] [PubMed]
27. Kanehisa M, Goto S. KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res 2000;28:27-30. [Crossref] [PubMed]
28. Friedman J, Hastie T, Tibshirani R. Regularization Paths for Generalized Linear Models via Coordinate Descent. J Stat Softw 2010;33:1-22. [Crossref] [PubMed]
29. Subramanian A, Tamayo P, Mootha VK, et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U S A 2005;102:15545-50. [Crossref] [PubMed]
30. Ritchie ME, Phipson B, Wu D, et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res 2015;43:e47. [Crossref] [PubMed]
31. Liberzon A, Birger C, Thorvaldsdóttir H, et al. The Molecular Signatures Database (MSigDB) hallmark gene set collection. Cell Syst 2015;1:417-25. [Crossref] [PubMed]
32. Liberzon A, Subramanian A, Pinchback R, et al. Molecular signatures database (MSigDB) 3.0. Bioinformatics 2011;27:1739-40. [Crossref] [PubMed]
33. Ru B, Wong CN, Tong Y, et al. TISIDB: an integrated repository portal for tumor-immune system interactions. Bioinformatics 2019;35:4200-2. [Crossref] [PubMed]
34. Wu S, Lv X, Li Y, et al. Integrated Machine Learning and Single-Sample Gene Set Enrichment Analysis Identifies a TGF-Beta Signaling Pathway Derived Score in Headneck Squamous Cell Carcinoma. J Oncol 2022;2022:3140263. [Crossref] [PubMed]
35. Ito K, Murphy D. Application of ggplot2 to Pharmacometric Graphics. CPT Pharmacometrics Syst Pharmacol 2013;2:e79. [Crossref] [PubMed]
36. Yang W, Soares J, Greninger P, et al. Genomics of Drug Sensitivity in Cancer (GDSC): a resource for therapeutic biomarker discovery in cancer cells. Nucleic Acids Res 2013;41:D955-61. [Crossref] [PubMed]
37. Maeser D, Gruener RF, Huang RS. oncoPredict: an R package for predicting in vivo or cancer patient drug response and biomarkers from cell line screening data. Brief Bioinform 2021;22:bbab260. [Crossref] [PubMed]
38. Mayakonda A, Lin DC, Assenov Y, et al. Maftools: efficient and comprehensive analysis of somatic variants in cancer. Genome Res 2018;28:1747-56. [Crossref] [PubMed]
39. Liu Z, Wang L, Guo C, et al. TTN/OBSCN 'Double-Hit' predicts favourable prognosis, 'immune-hot' subtype and potentially better immunotherapeutic efficacy in colorectal cancer. J Cell Mol Med 2021;25:3239-51. [Crossref] [PubMed]
40. Jiang P, Gu S, Pan D, et al. Signatures of T cell dysfunction and exclusion predict cancer immunotherapy response. Nat Med 2018;24:1550-8. [Crossref] [PubMed]
41. Reiser J, Banerjee A. Effector, Memory, and Dysfunctional CD8(+) T Cell Fates in the Antitumor Immune Response. J Immunol Res 2016;2016:8941260. [Crossref] [PubMed]
42. Lin Y, Ruze R, Zhang R, et al. Immunometabolic Targets in CD8(+) T Cells within the Tumor Microenvironment of Hepatocellular Carcinoma. Liver Cancer 2025;14:474-96. [Crossref] [PubMed]
43. Zhou J, Liu M, Sun H, et al. Hepatoma-intrinsic CCRK inhibition diminishes myeloid-derived suppressor cell immunosuppression and enhances immune-checkpoint blockade efficacy. Gut 2018;67:931-44. [Crossref] [PubMed]
44. Li JJ, Wang JH, Tian T, et al. The liver microenvironment orchestrates FGL1-mediated immune escape and progression of metastatic colorectal cancer. Nat Commun 2023;14:6690. [Crossref] [PubMed]
45. Lu Y, Sun Q, Guan Q, et al. The XOR-IDH3α axis controls macrophage polarization in hepatocellular carcinoma. J Hepatol 2023;79:1172-84. [Crossref] [PubMed]
46. St Paul M, Ohashi PS. The Roles of CD8(+) T Cell Subsets in Antitumor Immunity. Trends Cell Biol 2020;30:695-704. [Crossref] [PubMed]
47. Leung EL, Li RZ, Fan XX, et al. Longitudinal high-dimensional analysis identifies immune features associating with response to anti-PD-1 immunotherapy. Nat Commun 2023;14:5115. [Crossref] [PubMed]
48. Barsch M, Salié H, Schlaak AE, et al. T-cell exhaustion and residency dynamics inform clinical outcomes in hepatocellular carcinoma. J Hepatol 2022;77:397-409. [Crossref] [PubMed]
49. Liu Y, Yao R, Shi Y, et al. Identification of CD101 in Glioma: A Novel Prognostic Indicator Expressed on M2 Macrophages. Front Immunol 2022;13:845223. [Crossref] [PubMed]
50. Sumaiya K, Langford D, Natarajaseenivasan K, et al. Macrophage migration inhibitory factor (MIF): A multifaceted cytokine regulated by genetic and physiological strategies. Pharmacol Ther 2022;233:108024. [Crossref] [PubMed]
51. Wang Q, Zhao D, Xian M, et al. MIF as a biomarker and therapeutic target for overcoming resistance to proteasome inhibitors in human myeloma. Blood 2020;136:2557-73. [Crossref] [PubMed]
52. Zhu GQ, Tang Z, Huang R, et al. CD36+ cancer-associated fibroblasts provide immunosuppressive microenvironment for hepatocellular carcinoma via secretion of macrophage migration inhibitory factor. Cell Discov 2023;9:25. [Crossref] [PubMed]
53. Zhao H, Chen Q, Alam A, et al. The role of osteopontin in the progression of solid organ tumour. Cell Death Dis 2018;9:356. [Crossref] [PubMed]
54. Hirama M, Takahashi F, Takahashi K, et al. Osteopontin overproduced by tumor cells acts as a potent angiogenic factor contributing to tumor growth. Cancer Lett 2003;198:107-17. [Crossref] [PubMed]
55. Reuter S, Gupta SC, Chaturvedi MM, et al. Oxidative stress, inflammation, and cancer: how are they linked? Free Radic Biol Med 2010;49:1603-16. [Crossref] [PubMed]
56. Scharping NE, Rivadeneira DB, Menk AV, et al. Mitochondrial stress induced by continuous stimulation under hypoxia rapidly drives T cell exhaustion. Nat Immunol 2021;22:205-15. [Crossref] [PubMed]
57. Bantug GR, Hess C. The immunometabolic ecosystem in cancer. Nat Immunol 2023;24:2008-20. [Crossref] [PubMed]
58. Roth KG, Mambetsariev I, Kulkarni P, et al. The Mitochondrion as an Emerging Therapeutic Target in Cancer. Trends Mol Med 2020;26:119-34. [Crossref] [PubMed]
59. Weinberg SE, Chandel NS. Targeting mitochondria metabolism for cancer therapy. Nat Chem Biol 2015;11:9-15. [Crossref] [PubMed]
60. Nava Lauson CB, Tiberti S, Corsetto PA, et al. Linoleic acid potentiates CD8+ T cell metabolic fitness and antitumor immunity. Cell Metab. 2023;35:633-650.e9. [Crossref] [PubMed]
61. Yuan Z, Li Y, Zhang S, et al. Extracellular matrix remodeling in tumor progression and immune escape: from mechanisms to treatments. Mol Cancer 2023;22:48. [Crossref] [PubMed]
62. Du L, Su Z, Wang S, et al. EGFR-Induced and c-Src-Mediated CD47 Phosphorylation Inhibits TRIM21-Dependent Polyubiquitylation and Degradation of CD47 to Promote Tumor Immune Evasion. Adv Sci (Weinh) 2023;10:e2206380. [Crossref] [PubMed]
63. Li CX, Ling CC, Shao Y, et al. CXCL10/CXCR3 signaling mobilized-regulatory T cells promote liver tumor recurrence after transplantation. J Hepatol 2016;65:944-52. [Crossref] [PubMed]
64. Redjimi N, Raffin C, Raimbaud I, et al. CXCR3+ T regulatory cells selectively accumulate in human ovarian carcinomas to limit type I immunity. Cancer Res 2012;72:4351-60. [Crossref] [PubMed]
65. Kachler K, Holzinger C, Trufa DI, et al. The role of Foxp3 and Tbet co-expressing Treg cells in lung carcinoma. Oncoimmunology 2018;7:e1456612. [Crossref] [PubMed]
66. Santegoets SJ, Duurland CL, Jordanova ES, et al. Tbet-positive regulatory T cells accumulate in oropharyngeal cancers with ongoing tumor-specific type 1 T cell responses. J Immunother Cancer 2019;7:14. [Crossref] [PubMed]
67. Moreno Ayala MA, Campbell TF, Zhang C, et al. CXCR3 expression in regulatory T cells drives interactions with type I dendritic cells in tumors to restrict CD8(+) T cell antitumor immunity. Immunity 2023;56:1613-1630.e5. [Crossref] [PubMed]
68. Liu RX, Wei Y, Zeng QH, et al. Chemokine (C-X-C motif) receptor 3-positive B cells link interleukin-17 inflammation to protumorigenic macrophage polarization in human hepatocellular carcinoma. Hepatology 2015;62:1779-90. [Crossref] [PubMed]
69. Jia D, Wei L, Guo W, et al. Genome-wide copy number analyses identified novel cancer genes in hepatocellular carcinoma. Hepatology 2011;54:1227-36. [Crossref] [PubMed]
70. Quidville V, Alsafadi S, Goubar A, et al. Targeting the deregulated spliceosome core machinery in cancer cells triggers mTOR blockade and autophagy. Cancer Res 2013;73:2247-58. [Crossref] [PubMed]
71. Li FP, Liu GH, Zhang XQ, et al. Overexpressed SNRPB/D1/D3/E/F/G correlate with poor survival and immune infiltration in hepatocellular carcinoma. Am J Transl Res 2022;14:4207-28. [PubMed]