Uncovering immune cell heterogeneity in hepatocellular carcinoma by combining single-cell RNA sequencing with T-cell receptor sequencing

Abstract

BACKGROUND

Understanding the status and function of tumor-infiltrating immune cells is essential for improving immunotherapeutic effects and predicting the clinical response in human patients with carcinoma. However, little is known about tumor-infiltrating immune cells, and the corresponding research results in hepatocellular carcinoma (HCC) are limited.

AIM

To investigate potential biomarker genes that are important for the development of HCC and to understand how immune cell subsets react throughout this process.

METHODS

Using single-cell RNA sequencing and T-cell receptor sequencing, the heterogeneity and potential functions of immune cell subpopulations from HCC tissue and normal tissue adjacent to carcinoma, as well as their possible interactions, were analyzed.

RESULTS

Eight T-cell clusters from patients were analyzed and identified using bioinformatics, including six typical major T-cell clusters and two newly identified T-cell clusters, among which Fc epsilon receptor 1G⁺ T cells were characterized by the upregulation of Fc epsilon receptor 1G, tyrosine kinase binding protein, and T cell receptor delta constant, whereas metallothionein 1E⁺ T cells proliferated significantly in tumors. Differentially expressed genes, such as regulator of cell cycle, cysteine and serine rich nuclear protein 1, SMAD7 and metallothionein 1E, were identified as significantly upregulated in tumors and have potential as biomarkers. In association with T-cell receptor analysis, we inferred the clonal expansion characteristics of each T-cell cluster in HCC patients.

CONCLUSION

We identified lymphocyte subpopulations and potential biomarker genes critical for HCC development and revealed the clonal amplification of infiltrating T cells. These data provide valuable resources for understanding the response of immune cell subsets in HCC.

Key Words: Single-cell RNA sequencing; Paired T-cell receptor sequencing; Hepatocellular carcinoma; Immune cell subpopulations; Tumor-infiltrating immune cells

Core Tip: In this study, the heterogeneity of immune cell subpopulations, their potential functions, and their possible interactions in hepatocellular carcinoma (HCC) tissues and cancer-adjacent normal tissues were analyzed using single-cell RNA sequencing and T-cell receptor sequencing. Genes that may serve as biomarkers were identified by characterizing the major T-cell populations and 2 newly identified T-cell populations. In combination with the T-cell receptor analysis, we also inferred the clonal expansion characteristics of each T-cell cluster in HCC patients. In conclusion, we identified lymphocyte subpopulations and potential biomarker genes critical for HCC development and revealed the clonal expansion of infiltrating T cells.

INTRODUCTION

Primary liver cancer (PLC) exhibits significant heterogeneity in its development and progression, poses a global health challenge and is the second leading cause of cancer-related deaths worldwide[1]. China has the highest number of PLC cases primarily due to its large population and high prevalence of chronic hepatitis B virus infections[2]. Hepatocellular carcinoma (HCC) is the most common type of PLC, accounting for approximately 75% of cases[3]. The absence of early symptoms and the lack of highly sensitive and specific diagnostic techniques often result in patients being diagnosed at an advanced stage, missing the optimal window for treatment. Relevant data show that the 5-year survival rate of these patients is less than 12.5%[4,5]. Treatment failure, resistance to multiple drugs, and high mortality rates in HCC are significantly linked to the diversity of malignant cells and the tumor microenvironment (TME)[6].

Further research has provided new treatment options. Among these new options, immunotherapy has emerged as a promising approach that can enhance the body’s immune response and promote specific tumor immunity, aiming to treat tumors while minimizing recurrence and metastasis. However, the effectiveness of immunotherapy against solid tumors such as HCC still needs to be improved because of the challenges posed by the TME and the intricate mechanisms of tumor immune evasion[7]. Therefore, there is an urgent need to understand the cellular atlas and molecular composition of HCC and to identify key genetic targets to help research and develop new therapeutic strategies.

Notably, HCC, as a type of solid tumor, has a more complex TME composed of cellular (tumor-infiltrating immune cells and stromal cells), chemical (chemokines), and physical (extracellular matrix) components. These components interact to facilitate HCC progression and diminish the effectiveness of existing treatments[8]. Previous research on HCC has relied primarily on pathological assessments or bulk transcriptome sequencing, which provides only average gene expression levels across a cell population, fails to accurately identify the specific cell types and often misses important, less common subpopulations. In contrast, single-cell sequencing has high resolution and preserves information about heterogeneity within tissues as well as the transcriptional information of some rare cell populations[9]. In addition, it can specifically identify a class of cells in the TME and their corresponding gene expression profiles[10]. Single-cell sequencing is currently being extensively utilized in research related to the TME of HCC, mechanisms of tumor development, clonal evolution, and immunotherapy. For example, Zheng et al[11] identified novel T-cell exhaustion signature genes, such as layilin, pleckstrin homology like domain family A member 1, and synaptosome associated protein 47, using single-cell RNA sequencing (scRNA-seq) data of T cells from HCC patients. Moreover, they reported that exhausted CD8⁺ T cells and regulatory T cells (Tregs) were preferentially enriched and possibly clonally amplified in tumor tissues. However, a deeper understanding of the immune microenvironment composition and the mechanisms behind tumor immunosuppression in HCC is still needed. Therefore, identifying new tumor immunotherapies and finding new targets and effective biomarkers are highly important for the diagnosis and treatment of HCC through an in-depth understanding of the TME.

In this study, we combined scRNA-seq data from HCC with T-cell receptor (TCR) immunome library data. We examined the differences between normal tissue and tumor samples and analyzed the number of tumor-infiltrating lymphocytes within the TME along with their gene expression profiles. We identified lymphocyte subpopulations critical to HCC development, identified several potential biomarker genes, and revealed tumor immune heterogeneity in HCC. These findings will aid in clinical diagnosis and improve the effectiveness of immunotherapy and combination therapy for HCC.

MATERIALS AND METHODS

Tissue samples

One male patient and one female patient with HCC diagnosed by pathology were included in this study. All patients were diagnosed with stage I disease. Patient P0301 was positive for hepatitis B virus according to the hepatitis B surface antigen test. All patients did not receive chemotherapy or radiotherapy before tumor resection. The clinical characteristics of these patients are shown in Supplementary Table 1. We obtained paired fresh HCC tissues and adjacent normal liver tissues from two patients after cell isolation and single-cell suspension preparation. The adjacent normal tissues were at least 2 cm from the adjacent matched tumor tissue. This project was approved by the Research Ethics Committee of Nanjing Hospital, which is affiliated with Nanjing University of Traditional Chinese Medicine. All patients who participated in the experiment signed a written informed consent form.

Dissociation of tissue samples and preparation of single-cell suspensions

Tissues excised during surgery were immersed in Dulbecco’s modified eagle medium/nutrient mixture F-12 medium and transferred to the laboratory via refrigerated transfer. The fresh tissues were washed twice with 1 × phosphate buffered saline and cut into approximately 1-3 mm³ blocks in a 60 mm cell culture dish using sterile surgical scissors. The tissue blocks were transferred to lyophilization tubes. Next, 1 mL of preconfigured lyophilization solution (90% foetal bovine serum +10% dimethyl sulfoxide) was added, and the mixture was mixed gently to ensure that the tissue blocks were in full contact with the lyophilization solution. The samples were placed in a gradient cooler at room temperature and -80 °C overnight, covered with sufficient dry ice, and sent to Wuhan BGI Genomics. The tissue blocks were transferred into gentleMACS C tubes (130-093-237, Miltenyi Biotec, Germany) containing the enzyme mixture according to the protocol of the tumor dissociation kit (130-095-929, Miltenyi Biotec, Germany). The tissue separation process was performed using the gentleMACS dissociator (130-093-235, Miltenyi Biotec, Germany) by selecting the built-in program ‘h_tumor_01’. During this time, the C tube was rotated continuously for 30 minutes at 37 °C. Afterward, the cell suspension was filtered through a 70-μm filter placed over a 50-mL tube. The cell suspension was centrifuged at 300 × g for 7 minutes, and the entire supernatant was aspirated. Finally, the cell pellet was resuspended in Dulbecco’s modified eagle medium/nutrient mixture F-12.

Single-cell capture, library construction and transcriptome sequencing

The cells subjected to treatment were stained with 0.4% trypan blue. Following the enumeration of cells and the assessment of cell viability to ensure compliance with quality control standards, the library was constructed. The single-cell library was constructed using the Chromium™ Controller and Chromium™ single cell 3’ reagent version 2 kit (10x Genomics, Pleasanton, CA, United States). Using a microfluidic chip, beads and cells with cell barcodes are wrapped in a droplet. The cells are dissected so that the mRNA in the cells is attached to the cell barcode on top of the bead, forming single-cell gel bead-in-emulsion, in which the reverse transcription reaction is carried out, followed by cDNA library construction. The library fragment size was determined using an Agilent 2100 bioanalyzer instrument and quantified with quantitative polymerase chain reaction. The libraries were sequenced on the machine after passing quality control.

ScRNA-seq and data processing

The gene expression matrix was generated with 10x Genomics Cell Ranger v3.1.0 software (http://10xgenomics.com). Seurat (version 3.0.1) was used for scRNA-seq data preprocessing[12]. Cells with total unique molecular identifiers (UMIs) less than 200 and more significant than 7000, gene numbers less than 750, and a percentage of mitochondrion-derived UMI counts > 10% in a single cell were considered low-quality cells and removed. Principal component analysis was conducted, and the top 15 principal components and the first 2000 variable genes were selected for subsequent clustering and visualization steps. Seurat’s “FindClusters” function was used to identify the main cell clusters. The Louvain algorithm embedded in Seurat software was used for clustering, and the uniform manifold approximation and projection method was used to visualize the clustering results. Differential expressed genes (DEGs) between each cell type were calculated with the function ‘FindMarkers’ in Seurat, and genes with fold changes greater than 0.5 and adjusted P values lower than 0.05 were selected as the significant DEGs. Gene enrichment analysis of the DEGs was conducted using the R package clusterProfiler. Cell communication analysis was performed with CellChat.

Trajectory analysis

Trajectory analysis was performed using Monocle2 (version 2.8.0). We then analyzed differential gene expression using the differential gene test function to identify significant genes (Benjamini-Hochberg-corrected P value < 0.01). The cellular ordering of these genes was performed in an unsupervised manner. Trajectory construction was performed after dimensionality reduction and cell sorting with default parameters.

TCR analysis

The raw fastq files were processed with the CellRanger (3.0.1) vdj pipeline with default settings. For TCR analysis, cells with a single TCR α chain (TRA) and a single TCR β chain (TRB) or a heavy chain and a light chain were initially employed to construct cell lineages. Cells with the same nucleic acid sequence of the α- and β-chains were merged into the meta-cell, and the number of cells in each sample was calculated. Then, cells with the same V and J gene families and the same complementary determining region 3 (CDR3) length in both the α and β chains were included in the same group. The hamming distance for the nucleic acid sequence of CDR3 between each cell pair was calculated for each group. The quality threshold clustering method was used, and cells with a hamming distance of less than 2 were grouped into the same lineage. The clonotype frequency for each lineage was calculated. All lineages were then categorized into five different groups according to the cell number in each lineage: Hyperexpanded (100 < n ≤ 500), large (20 < n ≤ 100), medium (5 < n ≤ 20), small (1 < n ≤ 5), and single (0 < n ≤ 1). To view the relationships of different lineages in each sample, a network was constructed from the distance matrix calculated above with Gephi version 0.9.7. Specifically, the distance matrix for all cells in the same lineage was calculated first. Cell pairs with a distance less than 2 were considered closely connected cell pairs, and an edge was added between them.

RESULTS

Single-cell sequencing of tumor and adjacent normal samples

To ascertain the cellular composition of liver tumor tissue, the 10x Genomics Chromium platform was employed to generate scRNA-seq libraries from four liver samples, comprising two samples from HCC tissues and two from adjacent normal tissues (Figure 1A). Approximately 30000 cells were captured using the 10x Genomics Chromium microfluidic technique and subjected to RNA sequencing. All sequenced cells were classified into 17 cell clusters with Seurat[13] and manually annotated into 13 different cell types according to the expression of some classical marker genes (Figure 1B and C). The marker genes for each cluster are shown in Figure 1D. Cells exhibiting low CD3 expression levels and elevated natural killer (NK) cell granule protein 7 expression levels were classified as NK cells. Conversely, cells characterized by high human leukocyte antigen, CD14, CD80, and CD86 expression were identified as dendritic cells (DCs). B and plasma cells were annotated according to the expression of CD19, CD79A, CD79B, immunoglobulin heavy constant G, immunoglobulin heavy constant A, and immunoglobulin heavy constant M. As indicated in Figure 1A, the most common cell types were CD8⁺ T cells and CD4⁺ T cells, which accounted for nearly 80% of the total sequenced cells. A slight difference in cell composition was observed when the tumor tissue was compared with the adjacent normal tissue (Figure 1E). The adjacent normal tissue had a greater proportion of NK cells, granzyme H⁺ (GZMH⁺) CD8⁺ T cells, and activated CD8⁺ T cells, whereas the tumor tissue had a greater percentage of Tregs, metallothionein 1E⁺ (MT1E⁺) CD4⁺ T cells, and naive CD4⁺ T cells. In addition, a minor difference in the cellular composition was also observed between the two patients. In P0302, a larger fraction of B cells, plasma cells, and Tregs was observed in adjacent normal tissue compared with tumor tissue. In P0301, γδ T cells were more common in adjacent normal tissue.

Figure 1 Single-cell transcriptome clustering of all cells from hepatocellular carcinoma patients. A: Schematic diagram of the single-cell RNA sequencing experimental design; B: Uniform manifold approximation and projections of the combined single-cell RNA sequencing profiles of immune cells from control and tumor samples are plotted and colored according to samples; C: Cells are plotted and colored according to cell clusters. Each point corresponds to a single cell colored according to the cell cluster; D: A bubble plot of typical marker genes in each cluster is displayed. Darker colors indicate higher average expression; larger bubbles indicate a greater percentage of genes; E: Histogram of the abundance distribution of each cell subpopulation in the hepatocellular carcinoma and control groups. HCC: Hepatocellular carcinoma; TCR: T-cell receptor; UMAP: Uniform manifold approximation and projections; GZMH: Granzyme H; MT1E: Metallothionein 1E; NK: Natural killer; Treg: Regulatory T cells; DC: Dendritic cell.

Clustering and analysis of tumor-infiltrating T-cell subpopulations in HCC

Tumor-infiltrating T lymphocytes play a crucial role in tumor immunotherapy and constitute a key component of the HCC TME[14]. To further explore the ecology of T lymphocytes in the tumor immune microenvironment, we focused on analyzing the distribution and potential functions of these subpopulations. All T cells were reclustered into different subgroups according to the expression of the corresponding marker genes. Finally, eight different T-cell clusters were identified (Figure 2A), including activated CD8⁺ T cells, CD8⁺ tissue-resident memory T (TRM) cells, NK-like CD8⁺ T cells, naive CD4⁺ T cells, Tregs, Fc epsilon receptor 1G⁺ (FCER1G⁺) T cells, MT1E⁺ T cells and γδ T cells. TRM cells with higher expression of marker genes, such as interleukin (IL) 7 receptor and CD69, were predominantly observed in adjacent normal tissue (Figure 2B). Activated CD8⁺ T cells were characterized by high GZMK, CD8A, and CD8B expression (Supplementary Figure 1A).

Figure 2 Comprehensive analysis of cancer and paracancerous T cells in hepatocellular carcinoma patients. A: Cluster distribution of T-cell subsets on uniform manifold approximation and projections after subdivision; B: Marker gene heatmap of the subpopulations. The horizontal line represents the marker gene name, with the cells listed vertically and the upper color bar representing different subpopulations; C: Histogram showing the difference in the proportions of T-cell subsets between cancer patients and paracancerous patients; D: Scatter plot of differential expressed genes showing genes highly expressed in tumor 1_Activated CD8⁺ T cells, 2_TRM, 4_Activated CD8⁺ T cells and 5_Activated CD8⁺ T cells at P0301. Each dot represents a gene; E: Gene ontology (GO) enrichment analysis results showing the functional enrichment of upregulated differential expressed genes s in tumor tissues. The vertical axis is the GO entry, and the horizontal axis is the number of genes enriched in the GO entry. The redder color indicates greater significance. UMAP: Uniform manifold approximation and projections; TRM: Tissue-resident memory T cells; NK: Natural killer; Treg: Regulatory T cells; FCER1G: Fc epsilon receptor 1G; MT1E: Metallothionein 1E; BP: Biological process; CC: Cellular composition; MF: Molecular function.

Overall, the relative proportions of different CD8⁺ T-cell subtypes between the tumor and normal tissues did not differ much between the two patients (Figure 2C). Notably, TRM cells are not found in the circulatory system but are stable and widespread in human and mouse tissues[15]. Our findings illustrate the expression and localization of TRM cells within HCC tissues, and this information will facilitate subsequent specific functional investigations. The numbers of activated CD8⁺ T cells and NK-like CD8⁺ T cells were greater in adjacent normal tissues. In contrast, tumor tissue was mostly converted to an exhausted state, resulting in a decreased proportion, which is consistent with previous findings[16]. The findings of our comprehensive investigation support the varied distribution of CD8⁺ T cells within the TME of HCC. We further classified CD4⁺ T cells into two CD4⁺ T-cell subpopulations. Naive CD4⁺ T cells, which are derived mainly from tumor tissue, exhibit low of cytokine and effector gene expression levels. Given that certain inhibitory molecules may hinder the antitumor immune response, we subsequently assessed the expression levels of some known immune checkpoint-related genes in each cell type. In addition, Tregs overexpress the classical marker genes cytotoxic T-lymphocyte associated protein 4 (CTLA4) and T cell immunoreceptor with Ig and ITIM domains, as well as the costimulatory marker inducible T cell costimulator, and their expression is significantly greater in tumor tissues (Supplementary Figure 1B)[11].

In addition to CD8⁺ and CD4⁺ T cells, some nonclassical T cells, including γδ T cells, FCER1G⁺ T cells, and MT1E⁺ T cells, were also identified in our sequencing data. Although γδ T cells, which are characterized by elevated T cell receptor delta variable 2 and killer cell lectin like receptor C1 (KLRC1) expression, constitute a minor fraction of the overall T-cell population, their significant role in antitumor immunity is well documented. The preferential expression of GZMK in some γδ T cells also corresponds to previous findings in peripheral blood, but its role in tumor tissue still requires further exploration[17]. FCER1G⁺ T cells and MT1E⁺ T cells, with lower expression of both CD4⁺ and CD8⁺ T cells, also accounted for less than 10% of the total population. FCER1G⁺ T cells express NKT cell-related genes such as FCER1G, tyrosine kinase binding protein, and T cell receptor delta constant, suggesting that this group of cells may play similar functional roles.

Transcriptome heterogeneity of T-cell subsets in different HCC patients

As T cells usually play a special role in antitumor immune reactions[18], the differences in gene expression between tumor-derived T cells and adjacent tissues were determined. DEGs in tumors and adjacent normal tissues were highly variable between different patients. Regulator of cell cycle (RGCC), cysteine and serine rich nuclear protein 1 (CSRNP1), SMAD7, lamin A and basic helix-loop-helix family member E40 expression was commonly upregulated in P0301 tumor tissues, whereas the expression of MT2A, MT1E, and MT1X was upregulated in P0302 tumor tissues (Figure 2D, Supplementary Figure 1D). As shown in Figure 2D, the upregulation of dual specificity phosphatase 4 expression in 1_activated CD8⁺ T cells inactivated the mitogen-activated protein kinase signaling pathway, specifically the ERK, p38, and JNK pathways, via dephosphorylation. This mechanism is important for both cellular physiology and pathology.

Furthermore, circumstantial evidence suggests that dual specificity phosphatase 4 is involved in the progression of various malignancies, including HCC, ovarian carcinoma, esophageal carcinoma with rib metastasis, and pancreatic tumors[19]. The SMAD7 gene, which is highly expressed in all four types of T cells, functions as an inhibitor of transforming growth factor-β (TGF-β) by binding to type I TGF-β receptors, thereby obstructing TGF-β-mediated Smad signaling. Additionally, SMAD7 interacts with various intracellular proteins to modulate the TGF-β-independent signaling pathway, which is important for maintaining cellular homeostasis and regulating immune responses[20,21]. To examine the signaling pathways through which these genes exert their effects, gene ontology analysis was conducted. The findings indicated that the functions related to DEGs upregulated in the tumor tissue of P0301 were mostly enriched in “T-cell-mediated immunity” and “regulation of cell killing”. In contrast, the DEGs identified in P0302 were related to the “chemokine-mediated signaling pathway”, “regulation of lymphocyte migration” and “cytokine activity” (Figure 2E).

Currently, immune checkpoint inhibitors have achieved remarkable results in treating melanoma[22] and lymphoma[23], demonstrating their promising future in tumor immunotherapy. Additionally, we compared the differences in the expression of some known checkpoint-related genes between tumor and normal tissues in each cell cluster. Among these genes, only CTLA4, a gene associated with exhaustion, was highly expressed in the tumor-derived Tregs, indicating a depleted immune environment (Supplementary Figure 1B). Overall, we not only identified DEGs that were highly expressed in tumor-infiltrating T cells from two different patients but also confirmed the high expression of the immune checkpoint gene CTLA4 in tumor tissues.

Analysis of pseudotime trajectories for different T-cell subsets

Pseudotime trajectory analysis was conducted using Monocle 2 to elucidate the developmental pathways of CD8⁺ T cells and CD4⁺ T cells[24]. In P0301, the trajectory of CD8⁺ T-cell differentiation progressed from 5_activated CD8⁺ T cells to 1_activated CD8⁺ T cells, ultimately culminating in 4_activated CD8⁺ T cells (Figure 3A). Among them, 1_activated CD8⁺ T cells appeared to be in an intermediate stage of differentiation and extended into a new branch. The cells on this branch had the same differentiation period as the 1_activated CD8⁺ T cells. However, at P0302, 1_activated CD8⁺ T cells were connected to 4_activated CD8⁺ T cells and 5_activated CD8⁺ T cells in two different directions (Figure 3A), suggesting that these cells may be in different developmental locations and be functionally differentiated. Similarly, we analyzed the differentiation trajectory of CD4⁺ T cells. Most of the 0_Naive CD4⁺ T cells are located on the backbone of the pseudotime trajectory. Then, the 6_Tregs extend in different directions at the branching points of the trajectory plot (Figure 3B).

Figure 3 Analysis of the development trajectories of the CD8⁺ T and CD4⁺ T-cell subsets. A: Schematic diagram of the pseudotime trajectory of CD8⁺ T-cell subpopulations at P0301 (left) and P0302 (right). Each dot corresponds to a single cell, and each color represents a T-cell cluster; B: The same trajectory analysis was applied to CD4⁺ T cells; C: Heatmap of differential gene expression using branch point expression analysis comparing different cell fates of CD8⁺ T cells at P0301. Each row represents a gene, and the column represents a total of 100 bins that divide the pseudotime value from start to finish; D: Branch point expression analysis of the same CD4⁺ T cells in P0301.

The classification of three distinct activated CD8⁺ T-cell subsets, designated 1_Activated CD8⁺ T cells, 4_Activated CD8⁺ T cells, and 5_Activated CD8⁺ T cells, was conducted on the basis of their temporal trajectories. The marker genes associated with these subsets include CD8A, X-C motif chemokine ligand 1, and TRB variable 12-4 (TRBV12-4). The CD8 antigen, a cell surface glycoprotein predominantly expressed on cytotoxic T lymphocytes, functions as a coreceptor alongside TCRs, facilitating the recognition of antigens presented by antigen-presenting cells in conjunction with class I major histocompatibility complex I molecules. This homodimeric molecule supports activated lymphocyte survival and differentiation into memory CD8⁺ T cells and modulates the immune response by influencing antigen presentation and lymphocyte activation and differentiation[25]. X-C motif chemokine ligand 1 encodes a class C chemokine that is produced primarily by activated CD8⁺ T cells and plays a significant role in inflammatory and immune responses. It is particularly notable for its ability to induce T-cell chemotaxis, thereby promoting leukocyte migration and activation[26]. In contrast, TRBV12-4 is a component of the TCR complex that is potentially active within the plasma membrane and is involved in antigen recognition within the variable region of the TRB chain[27]. The specific categorization of these three activated CD8⁺ T-cell subsets indicates that CD8⁺ T cells influence the immune response through distinct mechanisms at various stages of development.

CD8⁺ T cells and CD4⁺ T cells have two different directions of differentiation; therefore, to investigate the differences in the fates of the two branches, we performed an analysis using the branched expression analysis modeling method. CD8⁺ T cells in cell fate 1 expressed more exhaustion-related genes than those in cell fate 2, indicating the final depletion phenotype (Figure 3C). On the other hand, the expression of ribosomal protein-associated genes, which play an important role in the selective regulation of the activation state of Tregs, was increased in CD4⁺ T cells in cell fate 2 (Figure 3D). During the intermediate phase, genes such as T cell immunoreceptor with Ig and ITIM domains, C-C motif chemokine ligand 5, CD160, KLRD1, and tumor necrosis factor (TNF), demonstrated increased expression in this transitional phase (Figure 3C). Conversely, Figure 3D indicates a downregulation of ribosomal proteins. Pseudotime trajectory analysis is fundamentally a simulation and predictive approach with inherent limitations. However, it enables the categorization of cells along a virtual timeline by assessing the similarity of their gene expression profiles. This methodology facilitates the inference of the dynamic processes underlying cellular transitions from one state to another. Consequently, it helps in elucidating the molecular mechanisms involved in cell fate determination by identifying novel cell types and investigating critical regulatory factors during developmental processes. Taken together, the results of the pseudotime trajectory analysis revealed the developmental trajectories of cells in the same subtypes and identified differentiated intermediate populations. Further studies targeting these intermediate populations will greatly improve the precision and efficiency of immunotherapeutic strategies.

Information exchange between different cell subpopulations

In addition to vertically exploring the developmental trajectory of cells, horizontal cell-to-cell communication is also an important method of regulating signal transduction across different cells. Here, the “CellChat” package was employed to analyze the patterns of cellular communication across various cell types and disease states. Among the signaling pathways upregulated in tumors, migration inhibitory factor (MIF)-[CD74 + C-X-C motif chemokine receptor 4 (CXCR4)] is widely active in various intercellular interactions. In contrast, C-type lectin domain-containing 2C (CLEC2C)-KLRB1 exhibited high intensity only in the interaction of γδ T cells with two populations of CD8⁺ T cells (Figure 4A). Each cell population exhibits distinct patterns of incoming and outgoing signaling in both adjacent normal and tumor tissues. In both patients, the 8_Activated CD8⁺ T cells and 4_Resting CD8⁺ T cells were identified as the cell types receiving the highest levels of incoming signals (Supplementary Figure 2A). A similar pattern was observed in multiple cell types with respect to the outgoing signals. Among all the ligand-receptor pairs, major histocompatibility complex I, MIF, CLEC, and integrin subunit beta 2 were the most frequent incoming and outgoing molecules in different cell types (Figure 4B and C). Compared with those of the control, greater amounts of TNF, LAMININ, fibronectin 1, and COLLAGEN were observed in the outgoing signals (Supplementary Figure 2B and C). In conclusion, CD8⁺ T cells received the greatest number of afferent signals, indicating that the immune system in tumor tissues was in an abnormal activation state. Furthermore, the pronounced activation of two specific signaling pathways, namely, MIF-(CD74 + CXCR4) and CLEC2C-KLRB1, indicates their potential significance in advanced HCC.

Figure 4 Cell communication analysis of CD8⁺ T and CD4⁺ T-cell subsets. A: Bubble map showing signaling pathways upregulated among cell subpopulations in tumors. The abscissa is the ligand-acceptor pair; the ordinate is the cell type alignment group. The tissue type is noted in the parentheses; B and C: Strength mapping of outgoing signaling pathways in different tissues for each cell subpopulation in P0301. DC: Dendritic cell; NK: Natural killer; Treg: Regulatory T cells.

Ecological analysis of DC and NK cell subpopulations in HCC

In addition to T cells, the gene expression profiles of DCs and NK cells were also analyzed. The successful generation of antitumor immune responses depends on the strong capacity of DCs to launch and regulate adaptive immune responses[28]. Understanding the subpopulation composition and gene expression profile of DCs is important for understanding their potential role in liver tumor immune responses[29]. Four DC subpopulations were identified from the total population of cells, including plasmacytoid DCs (pDCs), activated DCs, chromosome 1 open reading frame 54⁺ (C1ORF54⁺) myeloid DCs (mDCs), and unspecified DCs (Figure 5A). The cells exhibiting comparatively elevated GZMB, interferon regulatory factor 4, and interferon regulatory factor 8 levels were classified as pDCs[30]. DC activation was defined according to CD80, CD86, and human leukocyte antigen expression. The C1ORF54⁺ mDC cluster was named because increased C1ORF54 expression was observed[31,32]. An analysis of gene expression differences between tumor and normal samples revealed that only a limited number of genes were upregulated in DCs activated by tumor tissues from the P0301 sample. Of these genes, C-X-C motif chemokine ligand 2 (CXCL2), CXCL3, and CXCL8 MT1G are notable (Figure 5B). In addition, we assessed the expression levels of cytokines and cytokine receptors across these four different cell clusters. IL 2 receptor subunit gamma (IL2RG) and IL3RA were identified in pDCs. Relatively high expression levels of interferon gamma (IFNG) receptor and TNF-α were detected in both DC-activated and C1ORF54⁺ mDCs, whereas DC-activated cells also expressed relatively high levels of IL13RA, IL10RB, IL1B and IL18 (Supplementary Figure 3).

Figure 5 Single-cell immune landscape of dendritic cells and natural killer cells. A: Uniform manifold approximation and projections plots showing dendritic cell subpopulations; B: Scatter plot of differentially expressed genes showing genes highly expressed in activated 0_DC tumors. Each dot represents a gene; C: Uniform manifold approximation and projections plots showing subpopulations of natural killer cells; D: Histogram showing the percentage of natural killer cells from different sources of tissues. UMAP: Uniform manifold approximation and projections; DC: Dendritic cell; C1ORF54: Chromosome 1 open reading frame 54; PDC: Plasmacytoid dendritic cell; NK: Natural killer.

For NK cells, three different cell populations were identified from the total number of NKs recognized from the major cell type analysis, including two NK cell subpopulations (CD69^highCD16^highCD56^low NK cells and CD69^lowCD16^highCD56^low NK cells). Notably, a group of CD8⁺ T cells was also mixed into the clustering plot of the NK cells, presumably because the expression profiles of the two subsets were relatively similar (Figure 5C). All the NK cell populations expressed higher levels of granulysin, GZMA, GZMB, NK cell granule protein 7 and perforin 1, whereas IFNG expression was greater in the CD69highCD16highCD56 Low NK cells than in the CD69 LowCD16highCD56low NK cells. Regarding cytokine levels, CD69highCD16highCD56 Low NK cells expressed higher levels of TNF, IL2RG and IFNG receptor 1, and a modestly greater amount of IL10RA was detected in CD69^highCD16^highCD56^low NK cells. On the basis of previous studies, it has been hypothesized that this population of NK cells functions as helper cells and can enhance antigen presentation by secreting TNF-α or altering T-cell activity by producing IL-10[33]. Next, we compared the proportions of these two NK subpopulations in these patients. In P0301, CD69^highCD16^highCD56^low NK cells were the major subtype, accounting for nearly 75% of the total NK cells (Figure 5D). However, in P0302, the proportions of both cell types were almost identical. With respect to gene expression differences, few significant genes between tumor and normal tissues were identified.

TCR diversities and pairing frequencies of V and J genes

Antigen recognition and clonal T-cell expansion of clone-specific TCRs are essential for adaptive immunity[34]. To determine the extent of clonal T-cell diversity in HCC patients, a clonal analysis of the TCRs was performed to identify the lineage composition. In terms of the usage of different V gene families, the major V genes for the beta chain of the TRB in P0301 were TRBV6 and TRBV9, whereas the most abundant in P0302 were TRBV20 and TRB28 (Supplementary Figure 4A and B). In terms of the TRA, preferential usage of TRAV 1-2 was noted in P0301, whereas P0302 exhibited a more uniform distribution of TRAV genes, although TRAV 1-2 remained the most frequently utilized. To increase the diversity of the whole variable region, different V genes and J genes can undergo a recombination process. Here, the composition of different V-J pairs was similar for each patient (Supplementary Figure 4C and D). Additionally, TCR length is another important determinant of the specificity of TCRs. After analysis of the distribution of the CDR3 length, the TRA and TRB lengths were between 6 and 21, and 5 and 21, respectively. As the function of the TCR requires both the α chain and β chain, we next analyzed the information for these two chains in each cell. Complete α and β double chains were detected in 60% of our sequencing data (Supplementary Figure 4E and F).

With the whole V(D)J sequence assembled for each cell, we also wanted to know whether the clonotypes existed for each sample and if only the lineages containing more than two different nucleic acid sequences were included. As a result, most clonotypes were composed of only two cells, and the nodes representing higher UMI values were detected mainly in the tumor tissues, suggesting that the TME stimulated the selective expansion of T cells. Notably, some lineages were shared by both the tumor and adjacent normal T-cell-derived cells, but the node size was relatively small. On the basis of the number of different sequences that existed in each lineage, we also grouped all the lineages identified into five major groups (Figure 6A). The hyperexpanded lineages constituted the predominant group within each sample. Subsequent functional analyses of these tumor-derived cell clones will be instrumental in elucidating their potential roles within the TME.

Figure 6 Analysis of T-cell receptor distribution and clonality in hepatocellular carcinoma patients. A: Sectoral graph showing the extent of T-cell receptor (TCR) clone amplification in hepatocellular carcinoma patients. The different colors correspond to each of the five TCR clone amplification groups; B: Uniform manifold approximation and projections clustering map of T cells generated by combining TCR expression data; C and D: The expansion levels of five different types of TCR amplification groups in each cell subtype of P0301 control and tumor.

Integration of the immune repertoire and mRNA analysis

With both immune repertoire data and mRNA sequencing data, we also wanted to determine the expression of TCRs in each integrated cell cluster. Using the cell labels annotated via scRNA-seq, we found that the expression of TCRs was observed in most T-cell subtypes. However, few TCRs were detected in Tregs, γδ T cells and C7_T cells (Figure 6B). Given the older age of the patients, we hypothesize that the phenomenon of T-cell senescence led to a decrease in the number and diversity of TCRs, which was particularly severe in the mentioned cell populations[35]. To gain insight into T-cell expression profiles, we investigated the degree of expansion within each cell subtype. In P0301, hyperexpanded cell lineages were predominantly observed in 3_CD8⁺ T cells, 4_resting CD8⁺ T cells and 8_activated CD8⁺ T cells. For 0_naive CD4⁺ T cells, a larger proportion of “single” cells was noted as the main type. Among the 1_ activated CD8⁺ T cells and 4_ activated CD8⁺ T cells, the numbers of these five different types were approximately the same (Figure 6C and D). In P0302, the “hyperexpanded” cell lineages were the dominant type in 4_resting CD8⁺ T cells, 1_activated CD8⁺ T cells and 8_activated CD8⁺ T cells, and the largest proportion of single types was observed in 0_naive CD4⁺ T cells. As the most expanded lineages may represent highly diverse cells that participate in antitumor reactions, we next analyzed the cellular composition of the most abundant hyperexpanded cell clones. As shown in Figure 6C, the top clones were generally derived from several different cell subtypes. In P0301, the highest cell lineage consisted of approximately 40 distinct V(D)J sequences, and the major cell type was CD8⁺ T cells. Taken together, our data suggest that the composition of TCR clonotypes varies greatly between tissues. Naive CD4⁺ T cells have the fewest TCR clones because they are in the initial state. In contrast, tumor-infiltrating CD8⁺ T cells have a high rate of clonal amplification. Furthermore, different T-cell populations are not completely isolated and may undergo extensive state transitions.

DISCUSSION

PLC is a highly heterogeneous tumor disease that is induced by many inner and outer stimuli. Given the variability in cancer development, cellular responses may differ among patients. In our study, the transcriptome and the TCR sequence in both tumor and adjacent normal tissues were detected at the single-cell level. Ultimately, we identified eight distinct T-cell subsets in two patients, including activated CD8⁺ T cells, TRM cells, NK-like CD8⁺ T cells, naive CD4⁺ T cells, Tregs, FCER1G⁺ T cells, MT1E⁺ T cells and γδ T cells.

Primary HCC is associated with a high mortality rate, and early detection remains a significant challenge because of its insidious onset and subtle early symptoms. Additionally, treatment options for HCC are limited, and the prognosis of patients with advanced-stage disease is poor. Various immune checkpoint inhibitors are currently used clinically to enhance antitumor T-cell activity by blocking checkpoint inhibitory receptors[36,37]. However, the efficacy of these agents in solid tumors, particularly in HCC, is inconsistent. Many individuals do not show significant improvement in their disease after treatment, with a low objective response rate[38]. These findings indicate a deficiency in our understanding of the TME in HCC, significantly hindering the advancement of drug classes based on immune checkpoint inhibitors.

Tumor-infiltrating lymphocytes are critical to the success of immunotherapy. However, our understanding of their quantity, type, distribution within tumors, and interactions with one another is limited. In addition, the mechanisms behind these varied responses still need to be clarified. The wide application of scRNA-seq has provided a convenient strategy to study the biological mechanisms of tumor immune cells, including their heterogeneity, kinetics, and potential role in disease progression and response to immune checkpoint inhibitors[39-42]. With the progressive development and refinement of scRNA-seq technologies[43], we can analyze single-cell-level sequences of highly diverse genes such as TCRs. TCR sequencing can be used to define the clonality of T cells and track the dynamic relationships of these cells[44]. Therefore, we integrated scRNA-seq data with TCR immunome library data from HCC, indicating T-cell heterogeneity between normal and tumor samples. We also characterized the dynamic relationship of each T-cell subpopulation in the TME.

Tumor-infiltrating lymphocytes are mainly T lymphocytes[45], and T lymphocytes are divided into CD4⁺ T cells and CD8⁺ T cells according to their cell phenotype. On the one hand, these immune cells act as killers of tumor cells (e.g., CD8⁺ T cells and NK cells) and play a role in promoting tumor development on the other hand[46]. Research has demonstrated that immune cell infiltration within tumors contribute to tumor development through various mechanisms, including the secretion of growth factors[47], enabling tumor cells to evade inhibitory signals, inhibiting apoptosis[48], promoting angiogenesis[49], enhancing tumor metastasis[50] and exerting immunosuppressive effects[51]. In our study, most of the cells detected were CD8⁺ T cells and CD4⁺ T cells. These cells accounted for approximately 80% of the total sequenced cells, indicating that T cells are in a very active state of proliferation and differentiation. We analyzed and characterized six typical T-cell clusters in cancer vs normal paracancerous tissues from two HCC patients. Based on previous experience, TRM cells occupy tissues without recirculation and differ from recirculating central and effector memory T cells in terms of transcription, phenotype, and function. These cells play a key role in the defense against infection and cancer[52,53]. Our method also identified two newly discovered T-cell clusters not detected in previous bulk RNA sequence HCC studies. FCER1G⁺ T cells share marker genes with NKT cells, such as FCER1G, tyrosine kinase binding protein, and T cell receptor delta constant, indicating their similar functional types. MT1E⁺ T cells proliferate significantly in tumors, and MT1E expression is also upregulated in HCC tumor tissues according to differential gene analysis, which differs from the findings of previous studies[54] and warrants further exploration of its mechanism in the TME. With respect to the T-cell population, activated CD8⁺ T cells severely under-infiltrated the tumor tissue. In contrast, the percentage of Tregs was increased, constituting a specific immunosuppressive microenvironment in HCC. However, we did not detect the corresponding subpopulation of exhausted T cells, probably due to the small number of samples and the different annotation methods we used.

The characteristics of T-cell infiltration in HCC are closely related to DEGs. Several genes significantly upregulated in tumor tissues were identified in our study and may serve as potential biomarkers. RGCC, an important complement response gene that directly binds to and stimulates the kinase activity of the cell cycle protein-dependent kinases cell division control protein 2 and Akt[55], plays a critical role in cell proliferation, differentiation, tumorigenesis, and metastasis[56]. In previous reports, RGCC genes were highly expressed in a variety of solid tumors, including colorectal breast cancer[56-58]. We observed that RGCC gene expression is commonly upregulated in HCC tumor tissues and may be involved in mediating tumor cell proliferation and invasion[59]. Previous studies have indicated that the CSRNP family is positively associated with acute inflammatory responses and humoral immune response pathways and is a potential prognostic biomarker of overall survival in patients with renal clear cell carcinoma[60]. Our findings of elevated expression of tumor CSRNP1 suggest that it may play an important role in immune infiltration through T-cell infiltration and thus influence the immune microenvironment. The TGF-β superfamily plays a key role in the regulation of liver fibrosis and HCC[61,62]. TGF-β family cytokines initiate signal transduction by binding to receptor complexes, with SMAD7 acting as an antagonist of TGF-β signaling. In HCC, we detected elevated SMAD7 expression in tumors, which is consistent with the findings of Park et al[63]. Notably, high SMAD7 expression did not promote the corresponding malignant process but supported the malignant proliferative, survival and invasive states of HCC. These findings suggest that the function of SMAD7 is largely dependent on the microenvironment of HCC and is a highly valuable biomarker indicative of HCC progression.

Longitudinal pseudotime trajectory analysis and horizontal intercellular communication analysis have provided new insights into the developmental relationships and dynamic communication of different cell subpopulations. In conclusion, the results of the pseudotime trajectory analysis revealed that even if the cell subpopulations are annotated as identical, they may be in different developmental locations and may be functionally differentiated. We also identified key genes at the CD8⁺ and CD4⁺ T-cell differentiation nodes, and further interpretation of these genes will greatly improve the accuracy and efficiency of immunotherapeutic strategies. Furthermore, analysis of intercellular communication enabled us to focus on the MIF-(CD74+CXCR4) and CLEC2C-KLRB1 signaling pathways, which are significantly activated in HCC. Further investigation of these pathways may uncover essential mechanisms underlying tumor progression.

T cells recognize a wide range of antigens through a large TCR population. TCRs vary from cell to cell and are unique identifiers for T cells. Linking the TCR clonotype to the scRNA-seq landscape provides a better understanding of the clonal expansion of the T-cell population. We inferred the clonal expansion characteristics of each T-cell cluster in HCC patients. C3_CD8⁺ T cells, C4_resting CD8⁺ T cells, and C8_activated CD8⁺ T cells exhibited a “hyperexpanded” degree of clonal expansion. In contrast, C0 naive CD4⁺ T cells presented the most “single” phenotype, confirming their inactivated state, which coincided with the results of the pseudotime trajectory analysis. The most amplified cell subpopulations usually represent greater diversity. Our data suggest that the different T-cell subpopulations are not completely isolated and may undergo a wide range of state transitions and share the same ancestry. Notably, some γδ T cells were also observed in these clones. It is possible that some CD8⁺ T cells were misclassified as γδ T cells during the clustering steps or that the γ/δ chains were captured by the single-cell V(D)J sequencing kit.

CONCLUSION

Overall, this study provides a detailed picture of HCC-infiltrating lymphocytes in terms of their degree of aggregation, differentiation, and gene expression characteristics. We also identified several DEGs that could serve as potential biomarkers. However, further investigation is warranted to elucidate their specific roles and mechanisms, particularly given the constraints imposed by the limited sample size. In addition, our characterization of T cells in HCC will help improve the development of TME-targeted immunotherapies.