To reveal the proportional difference in the two subsets of γδ T cells from peripheral blood mononuclear cells (PBMCs) between the healthy population and HCC patients, the raw immune-phenotyping data of these two populations provided by the Shuangzhi Purui Medical Laboratory Co., Ltd. (Wuhan, China) were obtained and analyzed. The healthy population consisted of 170 disease-free adults (age range, 26–72; gender includes both female and male). The HCC population consisted of 42 patients (age range, 30–71; diagnosis, HCC I to IV including both primary and recurrence; gender includes both female and male). The antibodies used included anti-human CD3, anti-human TCR Vδ1, and anti-human TCR Vδ2 antibodies (BioLegend, San Diego, CA, USA). All examinations were conducted using BD FACSCanto cytometry (BD Biosciences, San Jose, CA, USA).

Tumoral mRNA-sequencing data (level 3) and corresponding clinical information of 371 HCC patients were obtained from The Cancer Genome Atlas (TCGA) dataset (https://portal.gdc.cancer.gov/) by following the related guidelines and policies. The mRNA-seq data of 50 paired peri-tumor tissue samples, together with 226 healthy liver tissue samples from the GTEx database (GTEx V8, https://gtexportal.org/home/datasets), were used as the control group, which described donors’ clinical information including gender, race, age, clinical stages, therapy history, and biospecimen collection. Protein–protein interaction was analyzed using the online database (https://string-db.org/). The R Project for Statistical Computing version 4.0.3 was used for graph plotting and analysis unless specified otherwise. In our analysis, HCC patients were divided into three groups (I, II, and III–IV) according to the clinical stage (pTNM) due to the limited amount of stage IV patients. To analyze the multi-gene correlation, the mRNA-seq dataset from NIH National Cancer Institute (https://tcga-data.nci.nih.gov/tcga/) was also included. The R software packages ggstatsplot and pheatmap were respectively used to produce two-gene and multi-gene correlation maps. Spearman’s analysis was used to describe the correlation.

As for ferroptosis analysis, the related marker genes were derived from previously published data (21). The marker genes for apoptosis were adopted from http://biocc.hrbmu.edu.cn/CancerSEA/home.jsp. The marker genes for pyroptosis were from published work (22). Then R packages ggplot2 and pheatmap were used for data visualization. Based on marker genes of apoptosis, ferroptosis, and pyroptosis, we wanted to know whether samples could be clustered into subgroups. The R packages ConsensusClusterPlus and pheatmap were used to generate clustering heatmaps. Moreover, the Kaplan–Meier survival analysis with log-rank test was used to plot survival curves among groups; p-values and hazard ratio (HR) with 95% CI were generated by log-rank tests and univariate Cox proportional hazards regression. For correlation between single gene and survival, R packages ggrisk, survival, survminer, glmnet, and timeROC were applied for analysis.

To evaluate immune cell infiltration in the context of apoptosis, ferroptosis, and pyroptosis, the R package immunedeconv was applied, which integrates six classes of algorithms, including TIMER, xCell, MCP-counter, CIBERSORT, EPIC, and quanTIseq. Meanwhile, the gene expression of inhibitory checkpoint molecules PDCD1(PD1), CD274(PDL1), PDCD1LG2(PDL2), CTLA4, LAG3, TIGIT, HAVCR2(TIM3), and SIGLEC15 in HCC tissues were intensively inspected. To analyze the correlation between “-optosis” related marker genes and 15 checkpoint molecules, both stimulatory checkpoint molecules (CD27, CD28, CD40, DNAM1, ICOS, and TNFRSF9) and inhibitory checkpoint molecules (CD274, CTLA4, HAVCR2, LAG3, PDCD1, PDCD1LG2, TIGIT, SIGLEC15, and BTLA) were included. To assay ICB among HCC subgroup samples, the potential ICB response was predicted with the Tumor Immune Dysfunction and Exclusion (TIDE) algorithm (23). p-Value <0.05 was considered statistically significant.

To comprehensively inspect the immune microenvironment of HCC TME, we firstly compared gene expression of immune inhibitory checkpoint molecules between tumor tissues and normal (peri-tumor) tissues, including CD274, CTLA4, HAVCR2, LAG3, PDCD1, PDCD1LG2, TIGIT, and SIGLEC15. It showed that all inhibitory checkpoint molecules except for CD274 are significantly upregulated in the HCC tumor group (Figures 1A, B). Then, 371 HCC samples were divided into 3 groups based on clinical stage (pTNM), and the checkpoint molecules were further compared among the 3 subgroups. It turned out that all checkpoint molecules are not statistically different among subgroups (Figure 1C). Further analysis revealed that HCC patients with different clinical stages respond to ICB therapy similarly (Figure 1D). We further checked the difference of immune cell infiltration between the HCC and normal groups by CIBERSORT score comparisons. We found that in the HCC group, scores of regulatory T cells (Tregs), resting macrophages M0, resting myeloid dendritic cells (DCs), and activated mast cells are strikingly increased. On the contrary, CIBERSORT scores of infiltrated γδ T cells, monocytes, macrophages M2, and resting mast cells are all depressed (Figure 1E). Furthermore, we inspected tissue infiltration of these immune cells among subgroups according to clinical stages (I, II, III–IV, and normal) and observed a similar pattern (Figure 1F) to Figure 1E.

Due to the excessive proliferation of HCC tumor cells and corresponding antitumor immunity and because cancer cell apoptosis routinely occurs in the TME, we first checked the expression differences of 63 apoptosis marker genes between HCC and normal tissues. We can see that 48 apoptosis-related marker genes are significantly upregulated in the HCC group (9 downregulated genes highlighted by blue arrows and 6 genes with no statistical difference) (Figure 2A). We then analyzed the correlation between each of the 63 marker genes and the overall survival (OS) rate and found that only 22 apoptotic genes are significantly correlated with the OS (p < 0.05). Specifically, 21 genes have negative correlations, while only one gene (PPP2R1B) has a positive correlation (Supplementary Figure 2A). Afterward, the multi-gene correlation among these 22 genes was calculated, as shown in Figure 2B. Moreover, we intensively analyzed those genes that have a better correlation with OS (p ≤ 0.01), including BCL10, E2F2, BAK1, BCAP31, CASP2, SNW1, and STK4. It shows that higher expressions of these 7 genes lead to poorer prognoses (Figure 2C). Further analysis revealed that expressions of these 7 genes among different clinical stages (I, II, and III–IV) are consistently and greatly upregulated as compared with the normal group, and the expression level dramatically increases from stage I to stage III–IV (Figure 2D). We then found that six of these seven genes (except BCAP31) are positively correlated with the expression level of checkpoint molecules and immune cell infiltration (Figures 2E, F; Supplementary Figure 2B).

Based on the expressions of those 22 genes that are correlated with OS, 371 HCC samples were further divided into 3 subgroups using the R package ConsensusClusterPlus (Supplementary Figure 2C), and the clinical information of these 3 groups is shown in the table. It shows that the expressions of most checkpoint molecules (except SIGLEC15) are significantly elevated in “G2” and “G3” compared with “G1.” CIBERSORT scores indicated that infiltrations of naïve B cells, CD8⁺ T cells, naïve CD4⁺ T cells, resting NK cells, monocytes, macrophage M1, and activated mast cells are statistically different among the three groups (p < 0.01). Further TIDE score analysis revealed that “G2” and “G3” have better ICB responses. However, OS probability shows no statistical difference among the three subgroups. Furthermore, 7 genes with a p ≤ 0.01 correlation with OS were used to construct a prognostic signature model using the least absolute shrinkage and selection operator (LASSO) Cox regression model. The obtained Riskscore equals (0.072) * BAK1 + (0.1257) * BCAP31 + (0.0895) * BCL10 + (0.2596) * E2F2 + (0.1424) * SNW1 (Supplementary Figure 2D). Based on this prognostic model, the correlation between the Riskscore and immune cell infiltration was further visualized (Figure 2G).

Given that ferroptosis can be conditionally triggered by chemotherapy or can spontaneously occur under oxidative stress in the HCC microenvironment, we then checked the overall expressions of ferroptosis marker genes between HCC and normal samples. It shows 19 upregulated and 3 downregulated genes in HCC (Figure 3A), and the gene correlation in HCC and normal tissues is shown in Figure 3B and Supplementary Figure 3A. We then calculated how expressions of these genes would affect patient prognosis, and we found that higher expressions of genes SLC1A5, CARS1, SLC7A11, RPL8, and TFRC lead to poorer OS (p < 0.01) (Figure 3C and Supplementary Figure 3B). Moreover, the expressions of these 5 genes are gradually enhanced from tumor stage (pTNM) I to III/IV (Figure 3D). Subsequently, the correlation between these 5 genes and checkpoint molecules or immune cells was analyzed, as shown in Figures 3E, F. Moreover, HCC samples were also sub-clustered into 4 groups based on the expression of 22 ferroptosis marker genes, showing the difference in immune cell infiltration as well as expressions of inhibitory checkpoint molecules among subgroups (Supplementary Figure 3C). TIDE score revealed that these 4 subgroups respond to ICB differently, and the OS probability among groups is statistically different (p < 0.001). Furthermore, a prognostic signature model was constructed, and the obtained Riskscore equals (0.1218) * CARS1 + (0.0359) * RPL8 + (0.0095) * TFRC + (0.1585) * SLC7A11 + (0.1652) * SLC1A5 (Supplementary Figure 3D), which was used to calculate the correlation between the Riskscore and immune cell infiltration (Figure 3G).

In addition to ferroptosis, pyroptosis is another type of inflammatory cell death. It is thus interesting to unveil the relationship between pyroptosis and the immune response in the HCC microenvironment. We found that, among 33 marker genes, 21 genes are upregulated, 6 genes are downregulated, and 6 genes have no statistical difference (Figure 4A). Further analysis revealed that only 6 genes (CASP3, GSDME, NLRC4, NLRP6, NOD1, and PLCG1) are statistically correlated with the OS rate of HCC patients (p < 0.05) (Figure 4B and Supplementary Figure 4A), and only GSDME, NLRP6, and NOD1 are significantly correlated with the OS (p < 0.01). Moreover, genes GSDME and NOD1 expressed higher in pTNM stage III–IV compared with stages I and II (Figure 4C). Further evaluation revealed that these 6 genes except NLRP6 are all positively correlated with checkpoint molecule expression and immune cell infiltration (Figures 4D, E). Additionally, all HCC samples can be grouped into two subclusters based on pyroptosis marker genes (Supplementary Figure 4B), and then expressions of checkpoint molecules, immune cell infiltration (TIMER and CIBERSORT scores), ICB, and the OS probability between two subgroups were all deciphered (Supplementary Figures 4C, D). The prognostic signature model was also generated by using the genes with p < 0.01, and the Riskscore equals (0.2413) * GSDME + (−0.2142) * NLRP6 + (0.0892) * NOD1, and the correlation between the Riskscore and immune cell infiltration was obtained (Figure 4F and Supplementary Figure 4E).

In the context of HCC TME, apoptosis, ferroptosis, and pyroptosis occur in an interconnecting manner, and their cross-talk may exist at both the gene and protein levels. We thus evaluated those marker genes that have a p < 0.01 correlation with the OS rate of HCC patients and found that these three types of “-optosis” genes do show cross-talk at the genomic level (Figure 5A and Supplementary Figure 5A). For example, NOD1 (pyroptosis gene) associates with both apoptosis genes (BCL10, CASP2, STK4, E2F2, and BAK1) and ferroptosis genes (TFRC, SLC1A5, and CARS1), and apoptosis gene BCL10 connects with pyroptosis genes (NOD1 and GSDME) and ferroptosis genes (SLC1A5, TFRC, CARS1, and SLC7A11). The protein–protein correlation analysis (STRING database, https://string-db.org/) further supports the analysis at the mRNA level, clearly showing the existing interactions among apoptosis, ferroptosis, and pyroptosis at the protein levels (Figure 5B and Supplementary Figure 5B). Next, we examined the protein–protein interactions between 15 checkpoint molecules and the “-optosis” genes that significantly correlate (p < 0.01 and <0.05, respectively) with the OS rate (Supplementary Figures 5C, D); the result indicated that only 5 checkpoint molecules (CD27, CD28, CD40, CTLA4, and PDL1) correlate with “-optosis” proteins. Afterward, 15 “-optosis” genes (significantly correlating with the OS rate, p < 0.01) were used to construct the prognostic signature model, and the obtained Riskscore (Supplementary Figure 5E) was correlated with immune cell infiltration (Figure 5C), indicating positive correlations with 6 main types of immune cells and specifically relative higher correlations with neutrophils, macrophages, and myeloid DCs (p_Spearman > 0.4).

Since cross-talk of apoptosis, ferroptosis, and pyroptosis was demonstrated at both the gene and protein levels, it is interesting to decode whether they conspiringly participate in inducing the immune-suppressiveness of HCC TME. We thus analyzed how each marker gene (p < 0.01 correlation with the OS) correlated with immune cell infiltration (Figure 5D). The result unequivocally indicates that tumor-infiltrating Tregs, macrophage M0, resting myeloid DCs, and activated mast cells are all significantly enriched in the context of expressions of the listed marker genes. Meanwhile, γδ T cell, monocytes, macrophage M2, resting mast cells, and neutrophils are all significantly suppressed. Importantly, both enriched and reduced immune cells significantly (p < 0.001) correlate with the “-optosis” marker genes regardless of their high or low expression in HCC tissue compared with normal tissue. Additionally, OS analysis generated by the prognostic signature model shows that higher expressions of marker genes (p < 0.01 correlation with the OS) of apoptosis, ferroptosis, and pyroptosis lead to poorer prognosis both individually (Figures 6A–C) and collaboratively (Figure 6D) (Supplementary Figures 2D, 3D, 4E, and 5E). These results again suggest that apoptosis, ferroptosis, and pyroptosis altogether participate in creating the immunosuppressive HCC TME.

Since the liver is one of the richest sources of tissue-resident γδ T cells, and our previous works unequivocally proved that adoptive transfer therapy of allogeneic Vδ2⁺ γδ T cells could safely and efficiently help control liver cancer progression or even achieve complete remission (4, 5), we thus focused on exploring how γδ T cell is repressed in the HCC TME (Figure 5D). We tried to inspect the underlying difference of circulating γδ T cells between the healthy population and HCC patients. We found that even though the global proportion of γδ T cells in CD3⁺ T cells is not statistically different between healthy and HCC populations, the Vδ1⁺ subset dramatically elevated, while the Vδ2⁺ subset strikingly reduced in the HCC population (p < 0.0001) (Figure 6E). This leads to a striking augmentation of the Vδ1⁺/Vδ2⁺ ratio in the HCC population. Further analysis revealed the significant reduction of NKG2D expression in the Vδ1⁺ subset, but not expressions of PD1, NKP30, and NKP46. For the Vδ2⁺ subset, expression of PD1 is greatly increased in the HCC group, implying depressed antitumor immunity of Vδ2⁺ γδ T cells. For other makers, NKG2D, NKP30, and NKP46 show no statistical difference between the healthy and HCC groups (Figures 6F, G).