Immunosuppressive landscape in hepatocellular carcinoma revealed by single-cell sequencing

Background/Aims: Hepatocellular carcinoma (HCC), accounting for 75-85% of primary liver cancer cases, is the third leading cause of cancer-related death worldwide. The purpose of this research was to examine the tumor immune microenvironment (TIME) in HCC.

Methods: We investigated the HCC TIME by integrated analysis of single-cell and bulk-tissue sequencing data to reveal the landscape of major immune cell types.

Results: Regulatory T(Treg) cells were found to be specifically distributed in the TIME of HCC. Several immune checkpoints, including TNFRSF4, TIGIT and CTLA4, were found to be uniquely overexpressed in Treg cells, and the glycolysis/gluconeogenesis pathway was enriched in Treg cells. We also discovered the presence of two NK-cell subsets with different cytotoxic capacities, one in an activated state with antitumor effects and another with an exhausted status. In addition, memory B cells in HCC were found to exist in a unique state, with high proliferation, low differentiation, and low activity, which was induced by overexpression of PRAP1 and activation of the MIF-CD74 axis.

Conclusions: We revealed the TIME landscape in HCC, highlighting the heterogeneity of major immune cell types and their potential mechanisms in the formation of an immunosuppressive environment. Hence, blocking the formation of the TIME could be a useful therapeutic strategy for HCC.

Introduction

Primary liver cancer is the sixth most commonly diagnosed cancer and the third leading cause of cancer death worldwide (1). Hepatocellular carcinoma (HCC), accounting for 75-85% of primary liver cancer cases, is the largest source of cancer-related disease burden. Specifically, approximately 70-80% of patients with HCC are diagnosed at an advanced stage (2), when curative therapies, including surgical resection, radiofrequency ablation, and liver transplantation, are not feasible. In addition, as 80~90% of HCC patients have concomitant cirrhosis, and the application of different therapeutic options might be limited because of the patient’s overall health status. Patients with advanced HCC typically receive multimodal therapy primarily comprising targeted therapy and immunotherapy (3).

Most tumors are complex and develop and evolve under robust selective pressure from their microenvironment, which includes nutrition-related, metabolic, immunological, and therapeutic components. Such pressure promotes the diversification of both malignant and nonmalignant compartments of the tumor microenvironment (TME), resulting in a degree of intratumoral heterogeneity (ITH) that enables aggressive disease progression and resistance to treatment (4). Multiple studies have proven that the immune component within the tumor, known as the tumor immune microenvironment (TIME), which includes immune cells, extracellular immune factors, and cell surface molecules, is closely associated with tumor development, recurrence, and metastasis (5, 6). The recognition of the importance of the TIME has given rise to another treatment option—immune checkpoint inhibitors (ICIs)—which have revolutionized cancer therapy. Blocking antibodies targeting immune inhibitory receptors such as CTLA-4, PD-1, and PD-L1 are by far the most extensively utilized immunotherapeutic drugs, with promising curative effects in a variety of cancers (7). Nivolumab, the first anti-PD1 monoclonal antibody approved as second-line therapy for HCC, has an objective response rate (ORR) of 15-20% and a disease control rate (DCR) of 58-64% (8). However, the ORR of tremelimumab (which targets CTLA-4) is only 7.2% (9). Thus, the intricate immune microenvironment of HCC needs to be explored, and doing so may improve the understanding of the effects of immunotherapy. Furthermore, the discovery of new therapeutic targets could pave the way for new treatments.

Different from traditional bulk RNA sequencing, single-cell RNA sequencing (scRNA-seq) has enabled researchers to investigate intratumoral heterogeneity at the single-cell level (10, 11). Thanks to the application of scRNA-seq, the field of single-cell genomics has massively expanded over the past few years, revealing surprising new insights into tissue biology and disease causes. Therefore, in this study, we performed a single-cell transcriptome analysis of seven HCC tumor tissue and paired normal tissue samples, thoroughly exploring the components and effector functions of major immune cell types and discovering novel biomarkers that might be employed as new therapeutic targets.

Materials and methods

Data acquisition and preparation

The Cancer Genome Atlas (TCGA) database was used to obtain the bulk gene expression profiles of HCC samples as well as the corresponding clinical data. There were 374 tumor samples and 50 paracancerous samples in the TCGA liver hepatocellular carcinoma (LIHC) cohort. Three fibrous lamellar carcinomas and seven mixed carcinomas in the 374 tumor samples were removed. In addition, samples from the same patients were removed. As a result, we ended up with 360 tumor samples. The scRNA-seq count matrix was downloaded from the GSE149614 dataset, which contains 21 samples, including 10 primary tumor samples, 8 corresponding peritumor liver samples, 1 metastatic lymph node sample and 2 portal vein tumor thrombus samples. Because our study focused on the heterogeneity between tumor and normal tissues, we selected corresponding tumor and normal samples. Since sample HCC07T only had data for 510 cells, HCC07T and the equivalent HCC07N sample were removed. Ultimately, a total of seven tumor samples and corresponding peritumor tissues were included in our study.

Single-cell analysis

The Seurat R package (12) was used for quality control (QC). Cells were removed when (a) RNA counts were fewer than 50; and (b) mitochondrial gene expression percentages were more than 5%. Data normalization was performed using the NormalizeData function in Seurat. The top 14 principal components and the top 1500 variable genes were chosen for the subsequent analysis. Then, cell clusters were detected by using Seurat’s FindClusters function (resolution = 0.5) and displayed via 2D t-distributed stochastic neighbor embedding (tSNE) (13). Cells from each cluster were compared to the annotated reference dataset using the SingleR package (14). According to the comparison results and the identified cell markers, each cluster was annotated. To reveal the differentiation of cell groups, we used Monocle2 (15), an R package built for pseudotime analysis. Furthermore, we estimated cell differentiation with online analysis via CytoTRACE (CytoTRACE (stanford.edu)), which is a robust approach for predicting cell differentiation with scRNA-seq data.

Functional enrichment analysis

Differentially expressed genes (DEGs) were selected by the FindMarker function in Seurat, and genes with adjusted p values < 0.05 were considered for subsequent analysis. Subsequently, Gene Ontology (GO) pathway enrichment analysis was performed utilizing clusterProfiler R packages (16). In this study, the gene set enrichment of each cell cluster was determined by gene set enrichment analysis (GSEA). Gene sets and KEGG pathways were obtained from MSigDB. Then, gene set variation analysis (GSVA) with 50 hallmark gene sets was performed among cell clusters. Furthermore, we used the AUCell R package to estimate the score of customized gene sets across cell groups.

InferCNV analysis

The InferCNV R package was developed to infer the copy number variation (CNV) of tumor cells based on the scRNA-seq matrix. Human genetic coordinate information was downloaded from https://data.broadinstitute.org/Trinity/CTAT/cnv/. Hepatocytes from normal tissue were considered the reference set, and the CNV of hepatoma cells was calculated. The final heatmap was generated after denoising.

SCENIC analysis

The SCENIC R package (17) is designed to assess gene regulatory networks and can guide the identification of transcription factors (TFs) and cell states. The GENIC3 method was used to extrapolate coexpression modules between transcription factors and candidate target genes. Each module, i.e., each regulon, contains a TF and its target genes. The cisTarget human motif database was utilized to enrich the gene signature, and targets in this signature were pruned according to the default set of cis-regulatory cues. Then, we assessed the activity of each regulon in each cell using the AUCell algorithm. The ComplexHeatmap and heatmap R packages were used for visualization.

Cell–cell communication

To study the cell group-to-cell group interactions, we employed CellPhoneDB (18), which can be used to analyze cell-to-cell communication at the molecular level. The receptor–ligand pairs with p <0.05 between each cell group were preserved, and we removed the receptor–ligand pairs within the same cell type. Finally, the number of receptor–ligand pairs was counted and visualized by Cytoscape.

Construction of a hepatoma cell differentiation trajectory

To verify the significance of genes related to cancer differentiation trajectory, we generated a hepatoma cell differentiation trajectory using the Monocle2 R package. Using the ConsensusClusterPlus package in R, 360 HCC samples in the TCGA LIHC cohort were divided into 3 clusters based on their expression of genes related to hepatoma cell evolution. Then, Kaplan–Meier (K-M) survival curves for each cluster were generated. Correlation analysis of clinical characteristics was performed for the 3 clusters.

Statistical analysis

Analyses between two groups were performed utilizing the Wilcoxon test. One-way ANOVA was used to compare three or more groups. Generally, statistical analyses were conducted by R studio (version 4.1.1) and GraphPad Software (version 8.0). The significance level was set at P<0.05.

Results

Single-cell transcriptome analysis identified cell compositions

Seven HCC samples and corresponding normal tissues were included in our study. A total of 39,667 cells passed QC, of which 21,121 were from tumor tissues and the rest were from normal tissues. These cells were then grouped into 36 clusters (Figure 1A), with the top five most significant genes highlighted and plotted on a heatmap (Figure 1B). Nine major cell types were identified in HCC. In addition to hepatocytes, smooth muscle cells, endothelial cells, and tissue stem cells, there were many immune cells (PRPTC⁺), including T cells, monocytes, natural killer (NK) cells, macrophages, and B cells (Figures 1D, E). When investigating the distributions of all 9 cell lineages, we noticed that the main cell type in tumor tissue was hepatoma cells, while in normal tissue, it was T cells (Figure 1C). Interestingly, there was an extreme difference in T cells and NK cells between intratumor and peritumor tissues. We believe that the depletion of T cells and NK cells in tumor tissue is associated with the formation of immunosuppression and is inextricably linked to tumor progression. In contrast to those of T and NK cells, the proportions of macrophages and monocytes were increased in HCC tumor tissues (Figure 1F). Neoantigens are tumor-specific antigens (TSAs) derived from the expression of mutated genes in cancer cells and are not present in normal tissues. These neoantigens may attract macrophages as well as DCs differentiated from monocytes that can engulf and present them, leading to the aggregation of macrophages and monocytes in HCC tumor tissues.

FIGURE 1

Figure 1 Cell-type classification in HCC. (A) t-SNE plot of 36 cell clusters. (B) Heatmap of the top five marker genes in each cell cluster. (C) Cell distribution in tumor and normal tissues. (D) Cell distribution in immune cells (PTPRC⁺) or nonimmune cells (PTPRC^-). (E) t-SNE plot exhibiting the cell types in HCC. The labels were automatically annotated by SingleR. (F) Bar plot of the distribution of all identified cell clusters: the proportion of each cell in tumor and normal samples (left) and in each patient (right). (G) The top five DEGs with the highest expression in tumor and normal tissues. (H) t-SNE plot of the top five DEGs.

Subsequently, we analyzed DEGs between tumor and normal tissues (Table S1) and selected the top five DEGs with the highest expression, which were SPP1, AKR1C2, AKR1B10, AGT, and EPHX1 in tumor tissues and NKG7, KLRD1, KLRB1, CCL5, and CST7 in normal tissues. The expression of the above genes was then marked in each cell lineage (Figures 1G, H). The highly expressed genes in normal tissues were mainly distributed in T and NK cells. NKG7, KLRD1, and KLRB1 are marker genes of T cells and NK cells, which is in line with the fact that the proportions of T and NK cells in tumor samples were significantly smaller than those in paracancerous tissues. In addition, EXPH1, AKR1C2 and AKRB10 were predominantly expressed in hepatoma cells. EPHX1 was first purified from the liver and is engaged in a variety of physiological activities, including lipid metabolism and detoxification of heterologous substances (19). Furthermore, it has been demonstrated that EPHX1 causes resistance to 5-fluorouracil in hepatoma cells and promotes chemoresistance in leukemia (20, 21). AKR1C1 and AKRB10 are both members of the human aldo-keto reductase family. Multiple studies have revealed that AKR1C1 is upregulated in various cancers, such as lung, breast and cervical cancers, and is associated with cancer metastasis and chemotherapy resistance (22, 23). Moreover, emerging experiments have also identified the role of AKR1B10 in HCC invasion and drug resistance (24). We validated that AKR1C1, AKRB10 and EXPH1 are promising biomarkers of HCC and potential therapeutic targets by scRNA-seq analysis.

Treg cells were enriched in HCC and had distinct metabolic characteristics

T cells have been found to be associated with prognosis in a variety of malignancies. During the past two decades, immunotherapies targeting T cells, including chimeric antigen receptor (CAR) T-cell therapy, adoptively transferred tumor-infiltrating lymphocytes, and immune checkpoint inhibitor (ICI) therapies, have been found to be effective in suppressing cancer growth (25). However, T cells have high heterogeneity, and different subtypes play varying effector functions within the TME. To investigate the heterogeneity of T cells, we clustered 12884 T cells into 15 subgroups, mainly four cell subtypes: CD8⁺Tem, CD4⁺Tem, CD8⁺Tcm, and Treg cells (Figures 2A, B). Among them, CD8⁺ Tem cells, which showed the highest proportion, were found in normal tissues, while CD4⁺ Tem, CD8⁺ Tcm, Treg cells were mainly found in HCC tissues (Figures 2C, D). This difference in distribution is likely related to the distinct functions of these T-cell subtypes: CD8⁺ Tem cells mostly reside in peripheral tissues and lymphoid tissues and play an effector role (26); CD4⁺ Tem cells travel to the infection site to exert killing functions (27), while Treg cells reside in tumor tissue and play an immunosuppressive role (28). The top three marker genes of the four T-cell subtypes are shown in Figure 2E.

FIGURE 2

Figure 2 Treg cells were enriched in HCC and had distinct metabolic characteristics. (A) t-SNE plot of 15 subclusters of T cells and (B) four major T-cell subtypes: CD8⁺ Tem, CD4⁺ Tem, CD8⁺ Tcm, and Treg cells. (C) Distribution of the four T-cell subtypes in tumor and normal samples and (D) in each patient. (E) Top five marker genes of the four T-cell subtypes. (F) Heatmap of immune checkpoint expression in the four T-cell subtypes. The row Z score represents the expression level. A differentiation trajectory of T cells, colored based on pseudotime (G) and cell type (H), is shown. (I) GSEA revealed that metabolic pathways were significantly enriched (FDR < 0.05) in Treg cells.

Next, we determined the expression levels of immune checkpoints in the four T-cell subtypes (Figure 2F). Several inhibitory checkpoints, such as TIGIT, CTLA4, TNFRSF4 and TNFRSF9, were only overexpressed in Treg cells. Despite being a minor component of immune cells, Treg cells can play a critical role in the TIME network. Treg cells suppress anticancer immunity, which hinders protective immunosurveillance of cancer and prevents the formation of an effective antitumor immune response, thereby promoting tumor development and progression. The immune checkpoints mentioned above that are only highly expressed in Treg cells can be considered potential therapeutic targets to restore immunity against cancer cells in HCC. Furthermore, we mapped the T cell differentiation trajectory by utilizing the Monocle2 R package (Figures 2G, H). Notably, Treg cells from HCC tissues were mainly found at the beginning and end of the differentiation trajectory, while CD8⁺ Tem cells were mainly found in the middle of the differentiation trajectory. Recent research on Treg cells has shown that targeting Treg cells in the end-stage of differentiation is an effective strategy to stimulate antitumor immunity (28). Strategies to deplete and control Treg cells by modulating Treg cell differentiation and development deserves further exploration.

In addition, the GSEA results showed that compared with other T-cell subtypes, Treg cells showed enrichment of metabolic pathways, such as glycolysis, gluconeogenesis, glutathione metabolism, and starch and sucrose metabolism (Figure 2I), indicating that Treg cells in HCC are not metabolically inhibited and have metabolic flexibility. Emerging studies have highlighted the critical role of metabolism in immune cells. Treg cells can proliferate with tumor development, maintaining their immunosuppressive effects because of increased metabolism in the TME. Tumors can also support Treg cell immunity by modulating metabolites, resulting in tumor immune escape (29, 30). In recent years, targeting metabolism has become a strategy to inhibit tumor development, and strategies employing DCs and macrophages have made progress (31). However, there is still much unknown about the connection between Treg cell metabolic pathways and the epigenetic control of gene expression. Our research will contribute to a better understanding of the tumor-promoting mechanism of aberrant Treg cell metabolism, as well as its effects on the differentiation and development processes.

Macrophages are related to the immunosuppressive environment in HCC

To study the interaction network in the HCC microenvironment, the Python-based cell–cell communication molecular analysis tool CellphoneDB was used to identify ligand–receptor pairs in the TME. Cytoscape was used to visualize the resulting network. Unexpectedly, except for epithelial cells and smooth muscle cells in HCC tissues, macrophages had the richest communication with other cell types (Figure 3A), revealing that macrophages play a significant role in the TME.

FIGURE 3

Figure 3 Macrophages are related to the immunosuppressive environment of the TME. (A) Interaction network constructed by CellPhoneDB. The size of the circles represents the interaction count. A larger size means more interaction with other cell types. (B) t-SNE plot of the eight subtypes of macrophages. (C) Distribution of the eight subtypes of macrophages in tumor and normal samples and (D) in each patient. (E) Violin plots of the expression of TREM2 in the eight subtypes of macrophages. (F) Kaplan–Meier survival curve for patients in the TCGA LIHC cohort. A log rank p value <0.05 was considered statistically significant. (G) Heatmap of immune checkpoint expression in macrophages. The row Z score represents the expression level. (H) Differences in hallmark pathway activities scored with GSVA. The t values calculated by a linear model are shown. (I) Heatmap of the AUC scores of the expression of transcription factors identified by SCENIC.

To assess the heterogeneity of macrophages, we clustered all macrophages into eight subclusters (Figure 3B). The phenotypes of clusters 1 to 8 are depicted in Figure 3B: macrophages: M−CSF/IFNγ/Pam3Cy; DCs: anti-CD40/VAF347; macrophages: M−CSF/IFNγ; macrophages: IL−4/Dex/cntrl; monocytes: leukotriene D4; monocytes: CD16⁺, monocytes; and monocytes: CD16⁻. By annotating the cell sources, we discovered that clusters 1, 2, 4, and 7 were derived from only cancer tissues, while clusters 3 and 8 were derived from only normal tissues. In addition, clusters 5 and 6 were derived from both cancer and normal tissues (Figures 3C, D). Therefore, cluster 3 was identified as a macrophage phenotype from normal tissues and thus was employed as a control. According to the expression of IFNγ and IL-4, cluster 1 was considered closer to M1 macrophages, and cluster 4 was considered closer to M2 macrophages. Although we could not accurately annotate the mononuclear and DC subsets, their functions should not be ignored.

As accumulating evidence has revealed a role of TREM2 in tumor-associated macrophages (TAMs) and myeloid-derived suppressor cells (MDSCs), we analyzed the expression of TREM2 in the eight clusters and found it to be high in cluster 1, cluster 2 and cluster 4 (Figure 3E), which were all derived from cancer tissues. TREM2 has been linked to poor prognosis and is a vital in inducing immunosuppression in the TME (32–36). To further study the clinical value of TREM2, we divided 360 HCC patients from the TCGA LIHC cohort into two groups according to the expression level of TREM2. Interestingly, the TREM2-overexpression group had a significantly worse prognosis (Figure 3F), consistent with previous research results, suggesting that TREM2 may be important for immunosuppression in the TME.

Subsequently, the expression of immune checkpoints among the 8 clusters was determined. Clusters 1 and 4 had higher expression levels of LAIR1 than the other clusters (Figure 3G). LAIR1 has been found to be associated with tumor immunosuppression and has recently been reported to block the LAIR1 and TGF-β signaling pathways to remodel the TME, making PD-L1-mediated tumor eradication possible (37).Our results also suggest that PD-L1-mediated tumor eradication could be a potential therapeutic strategy.

GSVA showed that the adipogenesis, fatty acid metabolism and bile acid metabolism pathways were enriched in cells from cluster 4, while the inflammatory response and the complement and interferon (IFN) pathways, which are the hallmarks of M2-like TAMs (as described in a previous study), were significantly suppressed (Figure 3H). It is worth mentioning that DCs were more significantly enriched in processes closely related to tumor metabolism, such as epithelial-mesenchymal transition, hypoxia, and angiogenesis.

Finally, we analyzed the expression of TFs in each subgroup. We found that compared to other clusters, cluster 1 and cluster 4 had higher expression of KLF6, PRDM1, STAT1, JUNB, and HES1 (Figure 3I). Researchers have reported that KLF6 plays an essential role in immunosuppression and modulates neutrophil maturation (38–42), while PRDM1 is a crucial epigenetic gene associated with T-cell terminal differentiation (43, 44). In other words, targeting PRDM1 would enable the generation of superior antitumor T cells. Furthermore, the IRF3/STAT1 pathway was reported to be associated with M2 polarization in a murine model of sarcoma (45). All of these findings indicate the vital role of macrophages in immunosuppression in the TME, as well as the potential mechanisms and TFs involved in this process.

Two distinct NK-cell subsets with different cytotoxic capacities in HCC

NK cells have recently been a hotspot in cancer immunology research. We detected 3151 and 187 NK cells in peritumor and tumor tissues, respectively, and clustered them into 11 subsets (Figures 4A, B). As previously described, the proportion of NK cells was dramatically decreased in tumor tissues compared to normal tissues. NK cells can kill cancer cells without MHC restriction and do not require preactivation by DCs or antibodies. NK cells play a critical role in the immunosurveillance protecting against cancer (46). The obvious decrease in NK cells in HCC suggests that cancer cells shape an environment that inhibits NK-cell proliferation, possibly by hypoxia or metabolic inhibition, which directly leads to resistance to attacks from NK cells.

FIGURE 4

Figure 4 Two distinct NK-cell subsets with different cytotoxic abilities in HCC. (A) t-SNE plot of the eleven subsets of NK cells. (B) Distribution of the eleven NK-cell subtypes in tumor and normal samples. (C) Volcano plot of the differentially expressed genes (DEGs) between tumor-derived NK cells and normal tissue-derived NK cells. The upregulated genes (log2(fold change) >1) are colored red, while the downregulated genes (log2(fold change) less than -1) are colored blue. Upregulated and downregulated genes are annotated. (D) Violin plots of the expression of several activating receptors, including CD160, NCR3, IFNG (IFNγ) and FASLG, in cluster 9, cluster 10 and normal tissues. (E) Differentiation trajectory of NK cells, colored for pseudotime (right) and cell subset (left). (F) Heatmap of the AUC scores of the expression of transcription factors identified by SCENIC.

DEGs of NK cells were identified between tumor and normal tissues (Figure 4C). Notably, HSP70 (including HSPA1A and HSPA1AB), a ubiquitous molecular chaperone, was highly expressed in tumor-derived NK cells. The major action of HSP70 is to maintain protein homeostasis and to mediate cytoprotective effects. Under adverse stress conditions, HSP70 can improve cell resistance to the environment and acts as a protector against stress (47). Tumor-resident NK cells are in a state of stress in the TME, and the antitumor effect of these NK cells can be supported by maintenance of hemostasis induced by upregulated HSP70. In addition, SERPINA1, which regulates hydrolase activity, and CLU, which stabilizes proteins, were upregulated (48, 49).

Interestingly, tumor-derived NK cells were distributed in clusters 9 and 10, while NK cells of other clusters came from normal tissues (Figures 4A, B). To investigate the heterogeneity of cluster 9 and cluster 10, we detected the expression of activating receptors of NK cells, including CD160, NCR3, IFNG (IFNγ), and FASLG, in cluster 9, cluster 10, and normal tissues (Figure 4D). CD160 potentiates NK-cell activation and cytotoxicity and induces the secretion of the cytokines IFN-γ, IL-6, IL-8 and TNF-α. The FasL protein is encoded by FASLG and is secreted by activated NK cells (50). After the ligand binds with the corresponding receptor on cancer cells, the apoptotic system of cancer cells is initiated. Notably, these genes were upregulated in cells from cluster 10 and normal tissue. Thus, compared to the cells of cluster 9, the cells of cluster 10 has more potent cytotoxicity, and the cells from cluster 9 exhibited an exhausted status. Compelling evidence has suggested that NK cells become exhausted in the presence of tumors and chronic infections, displaying low cytotoxicity and effector function (51). NK-cell exhaustion is also observed in HCC (52), yet the exact regulatory mechanisms have been poorly explored. Hence, we performed a pseudotime analysis (Figure 4E). It was found that cluster 10 cells were present at the beginning of the differentiation trajectory, while cluster 9 cells, which represented an exhausted NK-cell subset, existed at the end of the trajectory. We believe that cluster 10 cells gradually transform into cluster 9 cells with tumor infiltration and TME formation.

Subsequently, SCENIC analysis showed that the genes regulated by CEBPD, FOS, and JUN were significantly activated in cluster 10 (Figure 4F). JUN/FOS can regulate NK-cell immune activity through IF-2 (53), and CEBPD can promote NK-cell development and immunity (54). All of these findings confirmed the presence of two distinct NK-cell clusters with different cytotoxic capacities in HCC and suggest that the TME may gradually develop several mechanisms to suppress the immunotoxicity of NK cells. We revealed candidate TFs involved in this process, and the mechanisms of these TFs in NK cells are worth exploring. Furthermore, these TFs may serve as therapeutic targets.

Low activity state of memory B cells in HCC

A total of 1395 B cells were ultimately obtained after data processing and screening, and they were categorized into six clusters and annotated into four subtypes: B cells, immature B cells, memory B cells and plasma cells (Figures S1A and 5A). Notably, we noticed that almost all immature B cells were derived from normal tissues, while memory B cells were derived from tumor tissues (Figures 5B, C).

FIGURE 5

Figure 5 Low activity state of memory B cells in HCC. (A) t-SNE plot of the four B-cell subtypes in HCC. (B) Distribution of the four B-cell subtypes in tumor and normal samples and (C) in each patient. (D) Volcano plot of the DEGs between tumor-derived memory B cells and normal tissue-derived memory B cells. The upregulated genes (log2(fold change) >1) are colored red, while the downregulated genes (log2(fold change) less than -1) are colored blue. Upregulated and downregulated genes are annotated. (E and F) Differentiation status of different tissue-derived memory B cells as determined using CytoTRACE. A higher score indicates a lower degree of differentiation. (G) Pie charts of the proportions of cells in each stage of the cell cycle in different tissue-derived memory B cells. (H) Bar plot of the activated memory B-cell signature scores of memory B cells in tumor and normal samples. The signature scores were calculated based on the expression of activated memory B-cell-related genes, including CD86, AICDA, DHR59, EBI3, TBX21, KLF10, ZEB2, TFEC, ZBTB32, YBX3, CXCR3, ITGAX, and SIGLEC6. (I) Receptor–ligand pairs between hepatocytes and memory B cells.

Previous investigations have shown a strong infiltration of memory B cells in HCC, suggesting a specific immune response to the tumor (55–57). The prognostic impacts of memory B cells in HCC remain controversial. Memory B cells have been shown to have tumor-killing potential, with their presence in the HCC microenvironment indicating a good prognosis (58). However, other studies have reported that B cells restrain the antitumor response in several ways (59). Therefore, we conducted further analysis of memory B-cell subsets. The DEGs of memory B-cell subsets between tumor and normal tissues were identified (Figure 5D). Strikingly, we observed that PRAP1 was significantly upregulated in memory B cells derived from cancer tissue. PRAP1 has been demonstrated to be a novel p53 target gene that promotes cancer cell resistance to chemotherapy drugs such as 5-fluorouracil (5-FU) by cell cycle arrest to protect cells from apoptosis and contribute to cancer cell survival (60). In addition, PRAP1 can downregulate mitotic arrest deficient 1 (MAD1), which is a key factor in mitotic checkpoint signaling, leading to chromosomal instability and promoting the occurrence of HCC (61).

Next, we compared the differentiation potential of different tissue-derived memory B cells. Using CytoTRACE, we predicted a higher differentiation potential for tumor-derived memory B cells (Figures 5E, F). Similarly, a comparison of the cell cycles of the different tissue-derived memory B cell subtypes found that a larger proportion of tumor-derived memory B cells were in the G2/M phase, suggesting a more robust proliferative capacity, consistent with the low differentiation status of tumor-derived memory B cells (Figure 5G). To investigate functional heterogeneity, we assessed the expression of a gene set related to memory B-cell activation to predict the functional status (62). The results showed that there were significantly fewer memory B cells in the activated state in tumor tissues than in peritumor tissue (p = 0.011), which may be linked to the occurrence of an antitumor response (Figure 5H).

Given the crucial role of the crosstalk between the TME and cancer cells in tumorigenesis and progression, we compared receptor–ligand pairs between cancer cells and memory B cells. In particular, we discovered that macrophage migration inhibitory factor (MIF), which interacts with CD74⁺CXCR4⁺ and CD74⁺CD44⁺, was specifically expressed in tumor samples (Figure 5I). MIF has been confirmed to contribute to a variety of facets of tumor growth, including cell proliferation, differentiation, and angiogenesis. MIF can bind to its receptor CD74 in the TME, which is present on TAMs, DCs, Treg cells, and MDSCs, facilitating immunological escape and cancer growth (63, 64). For the first time, our study found that MIF may inhibit memory B-cell activity in HCC and that inhibiting the MIF-CD74 axis may be a new treatment strategy.

Hypermetabolism and immunosuppression in hepatoma cells

Hepatoma cells are derived from hepatocytes, so we compared CNV between hepatoma cells and normal hepatocytes. InferCNV analysis showed obvious CNV in hepatoma cells (Figure S2A). Meanwhile, hepatoma cells expressed high levels of APOA2, APOA1, AMBP, TTR, APOH, and ASGR1, which is consistent with what was previously observed by Sun Y et al. in HCC (Figure S2B) (65). Furthermore, hepatoma cells originating from the same patient tended to cluster together (Figure 6A), which indicated that there was significant heterogeneity between cancer cells across different patients. Heterogeneity is likely to lead to different responses to the same treatment among patients. In addition, multiple cancer cell subsets were discovered in the same lesion, revealing the existence of different tumor cell subtypes within tumor tissue (Figure S2C). Different subtypes of tumor cells exhibit differences in immune characteristics, growth rate, and invasive ability, resulting in different sensitivities to antitumor drugs.

FIGURE 6

Figure 6 Hypermetabolism and immunosuppression in hepatoma Cells. (A) Cancer cell distribution across all tumor samples. (B) Volcano plot of the DEGs between tumor-derived hepatoma cells and normal tissue-derived hepatocytes. The upregulated genes (log2(fold change) >0.5) are colored red, while the downregulated genes (log2(fold change) less than 0.5) are colored blue. (C) Gene Ontology analysis of the DEGs. The upregulated and downregulated DEGs are annotated. FDR <0.05 was considered significantly enriched. (D) Differences in hallmark pathway activities scored with GSVA. The t values calculated by a linear model are shown. (E) Differentiation trajectory of cancer cells colored for pseudotime. (F) TCGA LIHC patients were clustered into 3 clusters by ConsensusClusterPlus based on the expression of genes related to the cancer cell evolution states. (G) Kaplan–Meier survival curves for the patients in the 3 clusters. A log rank p value < 0.05 was considered statistically significant. (H) Heatmap of the AUC scores of the expression of transcription factors identified by SCENIC.

The DEGs between hepatoma cells and hepatocytes were detected by scRNA-seq data analysis, which can avoid the interference of other cells in the TME (Table S2 and Figure 6B). Strikingly, members of the metallothionein family, including MT1M, MT1H, MT1G, MT1F, MT1E, MT1X, and MT1A, were significantly downregulated in cancer cells. Metallothionein, a low-molecular-weight metal-binding protein, plays key roles in a range of biological processes in the human body, including participating in metal ion homeostasis and detoxification, regulating cell growth and proliferation, modulating immune inflammatory responses, and protecting the body from DNA damage and oxidative stress (66). Accumulating studies have shown the vital functions of MT1 proteins in tumor growth, invasion, and immune escape in kidney, breast, lung, and ovarian cancers (67). Recent studies have shown decreased MT1 expression in HCC (68). Promoter methylation can lead to significant repression of MT1G and MT1M expression. Moreover, MT1G and MT1M promoter methylation was found to be associated with an increased incidence of vascular invasion or metastasis (69). Therefore, MT1 proteins can serve as biomarkers of HCC, and these findings may lead to the development of new and effective therapeutic modalities if the mechanism and function of MT1 proteins in the development of HCC can be elucidated.

GO enrichment analysis revealed that the upregulated DEGs were mainly enriched in metabolism-related terms, such as ribonucleotide metabolic processes, purine nucleotide metabolic processes, precursor metabolites, and energy production (Figure 6C). This result indicates the presence of a hypermetabolic state in hepatoma cells. Tumor initiation and progression require metabolic reprogramming of cancer cells to alter their metabolic pathways to meet the biological needs of metabolism and biosynthesis. Therefore, high metabolism is a distinctive feature of tumor cells and promotes proliferation and invasion. However, immune functions were found to be suppressed in cancer cells. This is in accordance with our previous analysis showing that the microenvironment of HCC is immune-suppressed. The GSVA results were generally in line with the GO enrichment analysis results (Figure 6D), and we observed that the inflammatory response, as well as the interferon α-response and interferon γ-response pathways, were significantly repressed in cancer cells. The Mtorc1 and MYC signaling pathways, which lead to high metabolism and proliferation of cancer cells, were considerably enriched. Enrichment of the G2/M checkpoint and DNA repair pathways indicated the presence of considerable DNA damage in cancer cells. SCENIC analysis showed that heat shock factor 1 (HSF1), which maintains proteostasis in response to stress environments by inducing the expression of heat shock proteins, was significantly upregulated in cancer cells (Figure 6H). Recent studies have demonstrated the roles of HSF1 in tumorigenesis, such as inhibiting apoptosis, reprogramming metabolism, and regulating the TME (70, 71). In addition, DDIT3, a stress-induced TF, controls genes involved in cell cycle arrest and/or apoptosis as its primary function. A recent study discovered the role of DDIT3 in balancing glycolysis and oxidative phosphorylation during glutamine deprivation in cancer cells (72). Both of the above TFs involved in stress regulation have been found to be overexpressed in hepatoma cells and to be associated with poor prognosis (Figure S2D). In the HCC TME, DDIT3 and HSP1 may together regulate apoptosis and metabolism to promote tumor cell proliferation and tumor progression. Therefore, HSF1 and DDIT3 can serve as biomarkers of clinical prognosis and promising drug targets.

Cancer cells have multiple subgroups with varying degrees of differentiation that can exhibit differences in various aspects. To stratify HCC patients based on the degree of differentiation of cells and provide precise treatment recommendations, we generated a cancer cell differentiation trajectory using the Monocle2 R package (Figure 6E). Markers of HCC stem cells, including CD44, EPCAM, and THYI (CD90), were expressed at the beginning of the cancer cell differentiation trajectory (Figure S2E). DEGs were selected based on the different cancer cell differentiation states, which influence tumor cell differentiation and progression, as well as patient prognosis (Table S3). Therefore, we divided TCGA LIHC patients into 3 clusters based on these DEGs using ConsensusClusterPlus (Figure 6F). The K–M analysis illustrated that the patients in the cluster 2 group had considerably shorter OS than those in the cluster 1 and cluster 3 groups (p<0.05) (Figure 6G). When we analyzed the clinical features of the three clusters, we found that grade, T stage, and stage were significantly different among the patients in the three groups (Figure S2F). The cluster 2 group had a larger number of patients with advanced stage disease and tumor hypofractionation, indicating that cancer cells with varying degrees of differentiation have diverse invasive capacities that influence cancer progression.

Discussion

There remains a lack of treatment for HCC patients in the terminal stage who have missed the chance for surgery. Only a few patients can benefit from immunotherapy targeting immune inhibitory receptors of T cells, and some of them may develop drug resistance and stop responding to therapy. The TIME is the most dominant component of the TME and is critical for cancer progression as well as drug resistance. Immunosuppression is a characteristic of cancer, and crosstalk between cancer cells and immune cells ultimately leads to an environment that leaves patients with a weakened defense and often a worse prognosis (73). Because the TIME consists of a variety of immune cells and cytokines that frequently interact with each other, targeting one type of immune cell may lead to a series of changes in other microenvironments and relevant pathways. Combination therapies targeting several immune cells may thus be particularly effective in cancer treatment. Recent research has shown that the anti-tumor effects of Regorafenib, a multi-targeted tyrosine kinase inhibitor, are highly depended on its anti-tumor angiogenesis and anti-immunosuppressive properties, such as decrease of the TAM infiltration and enhancement of NK-cell cytolytic activity (74). Current research on the HCC TIME mainly focuses on T cells and macrophages, with few studies assessing the TIME as a whole. In this study, we investigated the TIME of HCC using single-cell sequencing and revealed the landscape of major immune cell types. The distribution and pathways of major cell subsets were analyzed, and some potential TIME-regulating mechanisms were discovered. The results of this study improve our understanding of the mechanisms of TIME formation and provide new ideas for immunotherapy.

Since the HCC TIME includes multiple immune cells and cytokines, we analyzed the predominant immune cell types in this study. We discovered a large number of Treg cells enriched in cancer tissue. Treg cells, the key component of immune homeostasis, maintain self-tolerance and suppress anticancer immunity. The presence of Treg cells is linked to tumor progression, aggressiveness, and metastasis. Their regulatory functions involve a wide range of immune cells other than T cells, including macrophages, DCs, neutrophils, NK cells, and innate lymphocytes. In addition to the typical expression of CD25, CD4, and FOXP3, Treg cells express a number of chemokine receptors and surface molecules, such as CTLA4, PD-1, and TIGIT, which are linked to anti-tumor immunity and may make them a direct target for ICI therapy. T-cell exhaustion, which is connected to tumor immune evasion, can be brought on by the overexpression of co-inhibitory receptors such CTLA4 and PD-1 (75). It has been reported that anti-CTLA4 and anti-PD1 monoclonal antibodies can deplete Treg cells and achieve clinical benefit. However, no other Treg cell-targeting therapies have been clinically proven to be effective, largely due to the difficulties associated with selectively targeting Treg cells (76). Our study revealed that, in addition to CTLA4, Treg cells uniquely overexpressed several immune checkpoints. Targeting these checkpoints may selectively deplete Treg cells, enabling precision immunotherapy. In addition, the glycolysis/gluconeogenesis pathway was enriched in Treg cells. Tumor cell metabolism is mainly dependent on glycolysis, which results in a high-lactate environment. Cancer cells can reprogram the metabolic pathways of Treg cells by communicating with them to induce adaptation to a high-lactate environment. Effector T cells are suppressed, while Treg cells are capable of surviving and proliferating within the TME; these Treg cells in turn support the formation of the TME, promoting immune escape of tumor cells and triggering tumor progression. Thus, altering the glycolytic pathway of Treg cells may induce activation of the entire immune system.

We also identified the presence of two NK-cell subsets with different activation statuses, one in an activated state with cancer-killing activity and another with exhaustion. The differentiation trajectory of NK cells in HCC was plotted using pseudotime analysis. The results showed a gradual depletion of NK-cell activity during TME formation. Recently, individuals with metastatic melanoma, NSCLC, and other tumor types have benefited significantly from ICIs. Strong T-cell anti-tumor immune responses are the basis of ICIs. Blocking NK cell-specific checkpoint receptors to reverse TME-induced NK cell exhaustion can revive NK cells’ direct cytotoxic activity against tumors and further initiate and enhance T cell-mediated adaptive anti-tumor immunity, according to a growing body of evidence (77). NK cell-based therapies are an effective complement to T cell therapy. Studies have found that combinations of anti-NKG2A monoclonal antibodies with anti-PD-1 monoclonal antibodies have showed promising outcomes in treating patients with advanced solid tumors (78). In this study, we observed the gradual functional depletion of NK-cell during the formation of TME by scRNA-seq analysis, and hindering this process is a new strategy for immunotherapy. In addition, the key TFs involved in this process was revealed, reversing the decreased expression of these TFs may restore the activity of NK cells.

Furthermore, memory B cells are enriched in tumor tissue, as observed by Tang B et al (56). However, the prognostic impacts of memory B cells in HCC remain controversial (58, 59). For the first time, our research found that memory B cells from HCC are in a unique state, with high proliferation, low differentiation, and low activity. Memory B cells are driven into this state by not only their own high expression of PRAP1 but also the influences of cancer cells via the MIF-CD74 axis. Currently, a detailed assessment of memory B cells in HCC has not been carried out, and our study suggests that the unique functional state and precise regulation of memory B cells may produce a good prognosis for patients.

Macrophages are at the center of cell–cell communication in the HCC TIME, and the prevailing perspective is that TAMs can be classified into M1-like (proinflammatory) and M2-like (anti-inflammatory) phenotypes. Our study suggests that the simple classification of M1-like and M2-like TAMs does not account for the complexity of TAMs. Some TAM subsets feature both M1-like and M2-like characteristics, and these subsets cannot be distinguished using conventional M1 and/or M2 macrophage markers. However, we discovered that TREM2 is uniquely expressed in tumor-derived DCs and macrophage subsets. Mounting evidence has revealed an immunosuppressive role of TREM2 in cancer (79, 80). TREM2 overexpression was associated with worse prognosis in the TCGA LIHC cohort. In addition, we revealed altered pathways and immune checkpoints in TAMs, which may provide new ideas for TAM-related immunotherapy.

Finally, using pseudotime analysis, we generated a hepatoma cell differentiation trajectory. For the first time, patients in the TCGA LIHC cohort were clustered into three distinct clusters based on their expression of differentiation-related genes. Notably, cluster 2 patients showed worse prognosis. Cancer cells in different differentiation stages exhibit differences in growth rate and invasive ability, resulting in different sensitivities to antitumor drugs. Therefore, individualized and precise treatment based on genetic testing will achieve increased clinical benefits.

In this study, we revealed the landscape of the TIME in HCC, highlighting the heterogeneity of major immune cell types and their potential mechanisms in the formation of an immunosuppressive environment. In summary, our research provides a novel theoretical basis for modulating the TIME in HCC and will aid the development of new immunotherapies.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

Author contributions

YB designed the research scheme and perform the bioinformatics analyzes. SHT downloaded and organized the gene expression matrix, clinical information in HCC. DPC performed the statistical analyzes. YB and DPC wrote the manuscript. QZ and YMZ critically revised the article for essential intellectual content and provided administrative support. All authors read and approved the final version of the manuscript. All authors reviewed and revised the manuscript. QZ and YMZ were the guarantors for this study.

Funding

This work was supported by the Tianjin Natural Science Foundation (20JCYBJC01310 and 21JCYBJC00320), the Tianjin Science and technology project (19ZXDBSY00010), and the Tianjin Health Science and technology project (TJWJ2021ZD002 and ZC20218).

Acknowledgments

We thank Yanlong Zhang for assistance with the statistical methods.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2022.950536/full#supplementary-material

Supplementary Figure 1 | t-SNE plot of the six subclusters of B cells

Supplementary Figure 2 | (A) CNV heatmaps with hierarchical clustering from the InferCNV analysis. The reference cells are normal tissue-derived hepatocytes, and the test cells are tumor-derived hepatoma cells. (B) t-SNE plot of the six markers of HCC, APOA2, APOA1, AMBP, TTR, APOH, and ASGR1. (C) t-SNE plot of the twenty subclusters of hepatoma cells. (D) Kaplan–Meier survival analysis of patients in the TCGA cohort. (E) Differentiation trajectory of cancer cells colored based on gene expression. (F) The proportions of patients with different T stages, stages, and grades in each cluster. A p value <0.05 was considered statistically significant.

Supplementary Table 1 | DEGs between tumor and normal tissues identified by Seurat’s FindClusters function

Supplementary Table 2 | DEGs between hepatoma cells and normal hepatocytes

Supplementary Table 3 | Genes related to hepatoma cell evolution based on pseudotime analysis

Abbreviations

CNV, copy number variation; DEG, differentially expressed genes; DCR, disease control rate; GSEA, gene set enrichment analysis; GSVA, gene set variation analysis; HCC, Hepatocellular carcinoma; ICIs, immune checkpoint inhibitors; ITH, intratumoral heterogeneity; MIF, macrophage migration inhibitory factor; MDSCs, myeloid-derived suppressor cells; ORR, objective response rate; QC, quality control; tSNE, t-Distributed Stochastic Neighbor Embedding; TIME, tumor immune microenvironment; TAMs, tumor-associated macrophages; TSAs, tumor-specific antigens.

References

1. Sung H, Ferlay J, Siegel RL, et al. Global cancer statistics 2020: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin (2021) 71(3):209–49. doi: 10.3322/caac.21660

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Hepatocellular carcinoma. Nat Rev Dis Primers (2021) 7(1):7. doi: 10.1038/s41572-021-00245-6

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Zhang T, Merle P, Wang H, Zhao H, Kudo M. Combination therapy for advanced hepatocellular carcinoma: do we see the light at the end of the tunnel. Hepatobiliary Surg Nutr (2021) 10(2):180–92. doi: 10.21037/hbsn-2021-7

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Vitale I, Shema E, Loi S, Galluzzi L. Intratumoral heterogeneity in cancer progression and response to immunotherapy. Nat Med (2021) 27(2):212–24. doi: 10.1038/s41591-021-01233-9

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Ciner AT, Jones K, Muschel RJ, Brodt P. The unique immune microenvironment of liver metastases: Challenges and opportunities. Semin Cancer Biol (2021) 71:143–56. doi: 10.1016/j.semcancer.2020.06.003

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Fu T, Dai LJ, Wu SY, Xiao Y, Ma D, Jiang YZ, et al. Spatial architecture of the immune microenvironment orchestrates tumor immunity and therapeutic response. J Hematol Oncol (2021) 14(1):98. doi: 10.1186/s13045-021-01103-4

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Bagchi S, Yuan R, Engleman EG. Immune checkpoint inhibitors for the treatment of cancer: Clinical impact and mechanisms of response and resistance. Annu Rev Pathol (2021) 16:223–49. doi: 10.1146/annurev-pathol-042020-042741

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Kudo M, Matilla A, Santoro A, Melero I, Gracián AC, Acosta-Rivera M, et al. CheckMate 040 cohort 5: A phase I/II study of nivolumab in patients with advanced hepatocellular carcinoma and child-pugh b cirrhosis. J Hepatol (2021) 75(3):600–9. doi: 10.1016/j.jhep.2021.04.047

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Kelley RK, Sangro B, Harris W, Ikeda M, Okusaka T, Kang YK, et al. Safety, efficacy, and pharmacodynamics of tremelimumab plus durvalumab for patients with unresectable hepatocellular carcinoma: Randomized expansion of a phase I/II study. J Clin Oncol (2021) 39:2991–3001. doi: 10.1200/JCO.20.03555

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Papalexi E, Satija R. Single-cell RNA sequencing to explore immune cell heterogeneity. Nat Rev Immunol (2018) 18:35–45. doi: 10.1038/nri.2017.76

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Ramachandran P, Matchett KP, Dobie R, Wilson-Kanamori JR, Henderson NC. Single-cell technologies in hepatology: new insights into liver biology and disease pathogenesis. Nat Rev Gastroenterol Hepatol (2020) 17:457–72. doi: 10.1038/s41575-020-0304-x

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Butler A, Hoffman P, Smibert P, Papalexi E, Satija R. Integrating single-cell transcriptomic data across different conditions, technologies, and species. Nat Biotechnol (2018) 36:411–20. doi: 10.1038/nbt.4096

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Narayan A, Berger B, Cho H. Assessing single-cell transcriptomic variability through density-preserving data visualization. Nat Biotechnol (2021) 39:765–74. doi: 10.1038/s41587-020-00801-7

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Aran D, Looney AP, Liu L, Wu E, Fong V, Hsu A, et al. Reference-based analysis of lung single-cell sequencing reveals a transitional profibrotic macrophage. Nat Immunol (2019) 20:163–72. doi: 10.1038/s41590-018-0276-y

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Qiu X, Mao Q, Tang Y, Wang L, Chawla R, Pliner HA, et al. Reversed graph embedding resolves complex single-cell trajectories. Nat Methods (2017) 14:979–82. doi: 10.1038/nmeth.4402

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Yu G, Wang LG, Han Y, He QY. clusterProfiler: an r package for comparing biological themes among gene clusters. OMICS (2012) 16:284–7. doi: 10.1089/omi.2011.0118

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Aibar S, González-Blas CB, Moerman T, Huynh-Thu VA, Imrichova H, Hulselmans G, et al. SCENIC: single-cell regulatory network inference and clustering. Nat Methods (2017) 14:1083–6. doi: 10.1038/nmeth.4463

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Vento-Tormo R, Efremova M, Botting RA, Turco MY, Vento-Tormo M, Meyer KB, et al. Single-cell reconstruction of the early maternal-fetal interface in humans. Nature (2018) 563:347–53. doi: 10.1038/s41586-018-0698-6

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Václavíková R, Hughes DJ, Souček P. Microsomal epoxide hydrolase 1 (EPHX1): Gene, structure, function, and role in human disease. Gene (2015) 571:1–8. doi: 10.1016/j.gene.2015.07.071

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Cheng H, Huang C, Tang G, Qiu H, Gao L, Zhang W, et al. Emerging role of EPHX1 in chemoresistance of acute myeloid leukemia by regurlating drug-metabolizing enzymes and apoptotic signaling. Mol Carcinog (2019) 58:808–19. doi: 10.1002/mc.22973

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Sun R, Dong C, Li R, Chu H, Liu J, Hao D, et al. Proteomic analysis reveals that EPHX1 contributes to 5-fluorouracil resistance in a human hepatocellular carcinoma cell line. Proteomics Clin Appl (2020) 14:e1900080. doi: 10.1002/prca.201900080

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Chang WM, Chang YC, Yang YC, Lin SK, Chang PM, Hsiao M. AKR1C1 controls cisplatin-resistance in head and neck squamous cell carcinoma through cross-talk with the STAT1/3 signaling pathway. J Exp Clin Cancer Res (2019) 38:245. doi: 10.1186/s13046-019-1256-2

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Zhu H, Chang LL, Yan FJ, Hu Y, Zeng CM, Zhou TY, et al. AKR1C1 activates STAT3 to promote the metastasis of non-small cell lung cancer. Theranostics (2018) 8:676–92. doi: 10.7150/thno.21463

PubMed Abstract | CrossRef Full Text | Google Scholar

24. DiStefano JK, Davis B. Diagnostic and prognostic potential of AKR1B10 in human hepatocellular carcinoma. Cancers (Basel) (2019) 11(4). doi: 10.3390/cancers11040486

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Kishton RJ, Sukumar M, Restifo NP. Metabolic regulation of T cell longevity and function in tumor immunotherapy. Cell Metab (2017) 26:94–109. doi: 10.1016/j.cmet.2017.06.016

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Kaech SM, Cui W. Transcriptional control of effector and memory CD8+ T cell differentiation. Nat Rev Immunol (2012) 12:749–61. doi: 10.1038/nri3307

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Luckheeram RV, Zhou R, Verma AD, Xia B. CD4⁺T cells: differentiation and functions. Clin Dev Immunol (2012) 2012:925135. doi: 10.1155/2012/925135

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Tanaka A, Sakaguchi S. Regulatory T cells in cancer immunotherapy. Cell Res (2017) 27:109–18. doi: 10.1038/cr.2016.151

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Menk AV, Scharping NE, Moreci RS, Zeng X, Guy C, Salvatore S, et al. Early TCR signaling induces rapid aerobic glycolysis enabling distinct acute T cell effector functions. Cell Rep (2018) 22:1509–21. doi: 10.1016/j.celrep.2018.01.040

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Watson MJ, Vignali P, Mullett SJ, Overacre-Delgoffe AE, Peralta RM, Grebinoski S, et al. Metabolic support of tumour-infiltrating regulatory T cells by lactic acid. Nature (2021) 591:645–51. doi: 10.1038/s41586-020-03045-2

PubMed Abstract | CrossRef Full Text | Google Scholar

31. O’Neill LA, Pearce EJ. Immunometabolism governs dendritic cell and macrophage function. J Exp Med (2016) 213:15–23. doi: 10.1084/jem.20151570

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Mulder K, Patel AA, Kong WT, Piot C, Halitzki E, Dunsmore G, et al. Cross-tissue single-cell landscape of human monocytes and macrophages in health and disease. Immunity (2021) 54:1883–900.e5. doi: 10.1016/j.immuni.2021.07.007

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Qiu H, Shao Z, Wen X, Jiang J, Ma Q, Wang Y, et al. TREM2: Keeping pace with immune checkpoint inhibitors in cancer immunotherapy. Front Immunol (2021) 12:716710. doi: 10.3389/fimmu.2021.716710

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Sayed FA, Kodama L, Fan L, Carling GK, Udeochu JC, Le D, et al. AD-linked R47H-TREM2 mutation induces disease-enhancing microglial states via AKT hyperactivation. Sci Transl Med (2021) 13:eabe3947. doi: 10.1126/scitranslmed.abe3947

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Shi D, Si Z, Xu Z, Cheng Y, Lin Q, Fu Z, et al. Synthesis and evaluation of (68)Ga-NOTA-COG1410 targeting to TREM2 of TAMs as a specific PET probe for digestive tumor diagnosis. Anal Chem (2022) 94:3819–30. doi: 10.1021/acs.analchem.1c04701

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Zhang H, Liu Z, Wen H, Guo Y, Xu F, Zhu Q, et al. Immunosuppressive TREM2(+) macrophages are associated with undesirable prognosis and responses to anti-PD-1 immunotherapy in non-small cell lung cancer. Cancer Immunol Immunother (2022) 1–12. doi: 10.1007/s00262-022-03173-w

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Horn LA, Chariou PL, Gameiro SR, Qin H, Iida M, Fousek K, et al. Remodeling the tumor microenvironment via blockade of LAIR-1 and TGF-β signaling enables PD-L1-mediated tumor eradication. J Clin Invest (2022) 132(8). doi: 10.1172/JCI155148

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Chen Q, Yin X, Jiang G, Huang J. Effects of circ-ABCB10/KLF6 on the progression of laryngocarcinoma. Panminerva Med (2021). doi: 10.23736/S0031-0808.21.04429-3

CrossRef Full Text | Google Scholar

39. Guo F, Du J, Liu L, Gou Y, Zhang M, Sun W, et al. lncRNA OR3A4 promotes the proliferation and metastasis of ovarian cancer through KLF6 pathway. Front Pharmacol (2021) 12:727876. doi: 10.3389/fphar.2021.727876

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Khoyratty TE, Ai Z, Ballesteros I, Eames HL, Mathie S, Martín-Salamanca S, et al. Distinct transcription factor networks control neutrophil-driven inflammation. Nat Immunol (2021) 22:1093–106. doi: 10.1038/s41590-021-00968-4

PubMed Abstract | CrossRef Full Text | Google Scholar

41. León X, Venegas M, Pujol A, Bulboa C, Llansana A, Casasayas M, et al. Predictive value of transcriptional expression of krüppel-like factor-6 (KLF6) in head and neck carcinoma patients treated with radiotherapy. Clin Transl Oncol (2021) 23:2507–12. doi: 10.1007/s12094-021-02651-4

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Wang Y, Lu K, Li W, Wang Z, Ding J, Zhu Z, et al. MiR-200c-3p aggravates gastric cell carcinoma via KLF6. Genes Genomics (2021) 43:1307–16. doi: 10.1007/s13258-021-01160-6

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Yoshikawa T, Wu Z, Inoue S, Kasuya H, Matsushita H, Takahashi Y, et al. Genetic ablation of PRDM1 in antitumor T cells enhances therapeutic efficacy of adoptive immunotherapy. Blood (2022) 139(14):2156–72. doi: 10.1182/blood.2021012714

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Zeng F, Luo L, Li D, Guo J, Guo M. KPNA2 interaction with CBX8 contributes to the development and progression of bladder cancer by mediating the PRDM1/c-FOS pathway. J Transl Med (2021) 19:112. doi: 10.1186/s12967-021-02709-5

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Hu J, Chen Z, Bao L, Zhou L, Hou Y, Liu L, et al. Single-cell transcriptome analysis reveals intratumoral heterogeneity in ccRCC, which results in different clinical outcomes. Mol Ther (2020) 28:1658–72. doi: 10.1016/j.ymthe.2020.04.023

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Shimasaki N, Jain A, Campana D. NK cells for cancer immunotherapy. Nat Rev Drug Discovery (2020) 19:200–18. doi: 10.1038/s41573-019-0052-1

CrossRef Full Text | Google Scholar

47. Multhoff G. Heat shock protein 70 (Hsp70): membrane location, export and immunological relevance. Methods (2007) 43:229–37. doi: 10.1016/j.ymeth.2007.06.006

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Boëlle PY, Debray D, Guillot L, Corvol H. SERPINA1 z allele is associated with cystic fibrosis liver disease. Genet Med (2019) 21:2151–5. doi: 10.1038/s41436-019-0449-6

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Wilson MR, Zoubeidi A. Clusterin as a therapeutic target. Expert Opin Ther Targets (2017) 21:201–13. doi: 10.1080/14728222.2017.1267142

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Myers JA, Miller JS. Exploring the NK cell platform for cancer immunotherapy. Nat Rev Clin Oncol (2021) 18:85–100. doi: 10.1038/s41571-020-0426-7

PubMed Abstract | CrossRef Full Text | Google Scholar

51. Zhang Q, Bi J, Zheng X, Chen Y, Wang H, Wu W, et al. Blockade of the checkpoint receptor TIGIT prevents NK cell exhaustion and elicits potent anti-tumor immunity. Nat Immunol (2018) 19(7):723–32. doi: 10.1038/s41590-018-0132-0

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Sun H, Huang Q, Huang M, Wen H, Lin R, Zheng M, et al. Human CD96 correlates to natural killer cell exhaustion and predicts the prognosis of human hepatocellular carcinoma. Hepatology (2019) 70(1):168–83. doi: 10.1002/hep.30347

PubMed Abstract | CrossRef Full Text | Google Scholar

53. Rakhshandehroo M, Gijzel SM, Siersbæk R, Broekema MF, de Haar C, Schipper HS, et al. CD1d-mediated presentation of endogenous lipid antigens by adipocytes requires microsomal triglyceride transfer protein. J Biol Chem (2014) 289:22128–39. doi: 10.1074/jbc.M114.551242

PubMed Abstract | CrossRef Full Text | Google Scholar

54. Ponti C, Gibellini D, Boin F, Melloni E, Manzoli FA, Cocco L, et al. Role of CREB transcription factor in c-fos activation in natural killer cells. Eur J Immunol (2002) 32:3358–65. doi: 10.1002/1521-4141(200212)32:12<3358::AID-IMMU3358>3.0.CO;2-Q

PubMed Abstract | CrossRef Full Text | Google Scholar

55. Rohr-Udilova N, Klinglmüller F, Schulte-Hermann R, Stift J, Herac M, Salzmann M, et al. Deviations of the immune cell landscape between healthy liver and hepatocellular carcinoma. Sci Rep (2018) 8:6220. doi: 10.1038/s41598-018-24437-5

PubMed Abstract | CrossRef Full Text | Google Scholar

56. Tang B, Zhu J, Li J, Fan K, Gao Y, Cheng S, et al. The ferroptosis and iron-metabolism signature robustly predicts clinical diagnosis, prognosis and immune microenvironment for hepatocellular carcinoma. Cell Commun Signal (2020) 18:174. doi: 10.1186/s12964-020-00663-1

PubMed Abstract | CrossRef Full Text | Google Scholar

57. Tang B, Zhu J, Zhao Z, Lu C, Liu S, Fang S, et al. Diagnosis and prognosis models for hepatocellular carcinoma patient’s management based on tumor mutation burden. J Adv Res (2021) 33:153–65. doi: 10.1016/j.jare.2021.01.018

PubMed Abstract | CrossRef Full Text | Google Scholar

58. Zhang Z, Ma L, Goswami S, Ma J, Zheng B, Duan M, et al. Landscape of infiltrating b cells and their clinical significance in human hepatocellular carcinoma. Oncoimmunology (2019) 8:e1571388. doi: 10.1080/2162402X.2019.1571388

PubMed Abstract | CrossRef Full Text | Google Scholar

59. Liu RX, Wei Y, Zeng QH, Chan KW, Xiao X, Zhao XY, 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. doi: 10.1002/hep.28020

PubMed Abstract | CrossRef Full Text | Google Scholar

60. Huang BH, Zhuo JL, Leung CH, Lu GD, Liu JJ, Yap CT, et al. PRAP1 is a novel executor of p53-dependent mechanisms in cell survival after DNA damage. Cell Death Dis (2012) 3:e442. doi: 10.1038/cddis.2012.180

PubMed Abstract | CrossRef Full Text | Google Scholar

61. Sze KM, Chu GK, Mak QH, Lee JM, Ng IO. Proline-rich acidic protein 1 (PRAP1) is a novel interacting partner of MAD1 and has a suppressive role in mitotic checkpoint signalling in hepatocellular carcinoma. J Pathol (2014) 233:51–60. doi: 10.1002/path.4319

PubMed Abstract | CrossRef Full Text | Google Scholar

62. Horns F, Dekker CL, Quake SR. Memory b cell activation, broad anti-influenza antibodies, and bystander activation revealed by single-cell transcriptomics. Cell Rep (2020) 30:905–13.e6. doi: 10.1016/j.celrep.2019.12.063

PubMed Abstract | CrossRef Full Text | Google Scholar

63. Tanese K, Hashimoto Y, Berkova Z, Wang Y, Samaniego F, Lee JE, et al. Cell surface CD74-MIF interactions drive melanoma survival in response to interferon-γ. J Invest Dermatol (2015) 135:2775–84. doi: 10.1038/jid.2015.204

PubMed Abstract | CrossRef Full Text | Google Scholar

64. Wirtz TH, Saal A, Bergmann I, Fischer P, Heinrichs D, Brandt EF, et al. Macrophage migration inhibitory factor exerts pro-proliferative and anti-apoptotic effects via CD74 in murine hepatocellular carcinoma. Br J Pharmacol (2021) 178:4452–67. doi: 10.1111/bph.15622

PubMed Abstract | CrossRef Full Text | Google Scholar

65. Sun Y, Wu L, Zhong Y, Zhou K, Hou Y, Wang Z, et al. Single-cell landscape of the ecosystem in early-relapse hepatocellular carcinoma. Cell (2021) 184:404–21.e16. doi: 10.1016/j.cell.2020.11.041

PubMed Abstract | CrossRef Full Text | Google Scholar

66. Dai H, Wang L, Li L, Huang Z, Ye L. Metallothionein 1: A new spotlight on inflammatory diseases. Front Immunol (2021) 12:739918. doi: 10.3389/fimmu.2021.739918

PubMed Abstract | CrossRef Full Text | Google Scholar

67. Si M, Lang J. The roles of metallothioneins in carcinogenesis. J Hematol Oncol (2018) 11:107. doi: 10.1186/s13045-018-0645-x

PubMed Abstract | CrossRef Full Text | Google Scholar

68. Krizkova S, Kepinska M, Emri G, Eckschlager T, Stiborova M, Pokorna P, et al. An insight into the complex roles of metallothioneins in malignant diseases with emphasis on (sub)isoforms/isoforms and epigenetics phenomena. Pharmacol Ther (2018) 183:90–117. doi: 10.1016/j.pharmthera.2017.10.004

PubMed Abstract | CrossRef Full Text | Google Scholar

69. Zhang S, Huang Z, Zhou S, Wang B, Ding Y, Chu JZ, et al. The effect and mechanism of metallothionein MT1M on hepatocellular carcinoma cell. Eur Rev Med Pharmacol Sci (2018) 22(3):695–701. doi: 10.26355/eurrev_201802_14295

PubMed Abstract | CrossRef Full Text | Google Scholar

70. Wang G, Cao P, Fan Y, Tan K. Emerging roles of HSF1 in cancer: Cellular and molecular episodes. Biochim Biophys Acta Rev Cancer (2020) 1874:188390. doi: 10.1016/j.bbcan.2020.188390

PubMed Abstract | CrossRef Full Text | Google Scholar

71. Kang H, Oh T, Bahk YY, Kim GH, Kan SY, Shin DH, et al. HSF1 regulates mevalonate and cholesterol biosynthesis pathways. Cancers (Basel) (2019) 11(9). doi: 10.3390/cancers11091363

CrossRef Full Text | Google Scholar

72. Li M, Thorne RF, Shi R, Zhang XD, Li J, Li J, et al. DDIT3 directs a dual mechanism to balance glycolysis and oxidative phosphorylation during glutamine deprivation. Adv Sci (Weinh) (2021) 8(11):e2003732. doi: 10.1002/advs.202003732

PubMed Abstract | CrossRef Full Text | Google Scholar

73. Jin MZ, Jin WL. The updated landscape of tumor microenvironment and drug repurposing. Signal Transduct Target Ther (2020) 5:166. doi: 10.1038/s41392-020-00280-x

PubMed Abstract | CrossRef Full Text | Google Scholar

74. Granito A, Forgione A, Marinelli S, Renzulli M, Ielasi L, Sansone V, et al. Experience with regorafenib in the treatment of hepatocellular carcinoma. Therap Adv Gastroenterol (2021) 14:17562848211016959. doi: 10.1177/17562848211016959

PubMed Abstract | CrossRef Full Text | Google Scholar

75. Granito A, Muratori L, Lalanne C, Quarneti C, Ferri S, Guidi M, et al. Hepatocellular carcinoma in viral and autoimmune liver diseases: Role of CD4+ CD25+ Foxp3+ regulatory T cells in the immune microenvironment. World J Gastroenterol (2021) 27(22):2994–3009. doi: 10.3748/wjg.v27.i22.2994

PubMed Abstract | CrossRef Full Text | Google Scholar

76. Togashi Y, Shitara K, Nishikawa H. Regulatory T cells in cancer immunosuppression - implications for anticancer therapy. Nat Rev Clin Oncol (2019) 16:356–71. doi: 10.1038/s41571-019-0175-7

PubMed Abstract | CrossRef Full Text | Google Scholar

77. Zhang C, Liu Y. Targeting NK cell checkpoint receptors or molecules for cancer immunotherapy. Front Immunol (2020) 11:1295. doi: 10.3389/fimmu.2020.01295

PubMed Abstract | CrossRef Full Text | Google Scholar

78. André P, Denis C, Soulas C, Bourbon-Caillet C, Lopez J, Arnoux T, et al. Anti-NKG2A mAb is a checkpoint inhibitor that promotes anti-tumor immunity by unleashing both T and NK cells. Cell (2018) 175(7):1731–43.e13. doi: 10.1016/j.cell.2018.10.014

PubMed Abstract | CrossRef Full Text | Google Scholar

79. Katzenelenbogen Y, Sheban F, Yalin A, Yofe I, Svetlichnyy D, Jaitin DA, et al. Coupled scRNA-seq and intracellular protein activity reveal an immunosuppressive role of TREM2 in cancer. Cell (2020) 182:872–85.e19. doi: 10.1016/j.cell.2020.06.032

PubMed Abstract | CrossRef Full Text | Google Scholar

80. Molgora M, Esaulova E, Vermi W, Hou J, Chen Y, Luo J, et al. TREM2 modulation remodels the tumor myeloid landscape enhancing anti-PD-1 immunotherapy. Cell (2020) 182:886–900.e17. doi: 10.1016/j.cell.2020.07.013

PubMed Abstract | CrossRef Full Text | Google Scholar