- Open Access
Single-cell transcriptomics reveals the role of Macrophage-Naïve CD4 + T cell interaction in the immunosuppressive microenvironment of primary liver carcinoma
Journal of Translational Medicine volume 20, Article number: 466 (2022)
Liver carcinoma generally presents as an immunosuppressive microenvironment that promotes tumor evasion. The intercellular crosstalk of immune cells significantly influences the construction of an immunosuppressive microenvironment. This study aimed to investigate the important interactions between immune cells and their targeting drugs in liver carcinoma, by using single-cell and bulk transcriptomic data.
Single-cell and bulk transcriptomic data were retrieved from Gene Expression Omnibus (GSE159977, GSE136103, and GSE125449) and The Cancer Genome Atlas (TGCA-LIHC), respectively. Quality control, dimension reduction, clustering, and annotation were performed according to the Scanpy workflow based on Python. Cell–cell interactions were explored using the CellPhone database and CellChat. Trajectory analysis was executed using a partition-based graph abstraction method. The transcriptomic factors (TFs) were predicted using single-cell regulatory network inference and clustering (SCENIC). The target genes from TFs were used to establish a related score based on the TCGA cohort; this score was subsequently validated by survival, gene set enrichment, and immune cell infiltration analyses. Drug prediction was performed based on the Cancer Therapeutics Response Portal and PRISM Repurposing datasets.
Thirty-one patients at four different states, including health, hepatitis, cirrhosis, and cancer, were enrolled in this study. After dimension reduction and clustering, twenty-two clusters were identified. Cell–cell interaction analyses indicated that macrophage-naive CD4 + T cell interaction significantly affect cancerous state. In brief, macrophages interact with naive CD4 + T cells via different pathways in different states. The results of SCENIC indicated that macrophages present in cancer cells were similar to those present during cirrhosis. A macrophage-naive CD4 + T cell (MNT) score was generated by the SCENIC-derived target genes. Based on the MNT score, five relevant drugs (inhibitor of polo-like kinase 1, inhibitor of kinesin family member 11, dabrafenib, ispinesib, and epothilone-b) were predicted.
This study reveals the crucial role of macrophage-naive CD4 + T cell interaction in the immunosuppressive microenvironment of liver carcinoma. Tumor-associated macrophages may be derived from cirrhosis and can initiate liver carcinoma. Predictive drugs that target the macrophage-naive CD4 + T cell interaction may help to improve the immunosuppressive microenvironment and prevent immune evasion. The relevant mechanisms need to be further validated in experiments and cohort studies.
Liver cancer remains the fourth most common cause of tumor-mediated death and ranks sixth in terms of incidence worldwide . As the first-line treatment for liver cancer, sorafenib prolongs survival modestly compared with placebo, conferring very limited survival benefits . Several immune checkpoint inhibitors have been approved for treating various cancers. Nivolumab, pembrolizumab, and atezolizumab, immune checkpoint inhibitors of the PDL1–PD1 pathway, have shown impressive responses against tumors and have a manageable safety profile in liver cancer [3,4,5]. However, it has been reported that patients with nonalcoholic steatohepatitis (NASH)-mediated liver cancer may not benefit from anti-PD1 treatment; besides, NASH-related aberrant T-cell activation is a potential cause . The immune context should be recognized as an important determinant of immunotherapy efficacy. To identify patients who would benefit most from immunotherapeutic interventions, it is important to elucidate the composition and properties of cells in the tumor microenvironment (TME).
Most patients with liver cancer have cirrhosis or chronic hepatitis . A complex balance exists between immunity and tolerance in healthy livers, and with progression from hepatitis to cirrhosis and liver cancer, the liver immune microenvironment gets changed. With functions of maintaining tolerance, immune suppression, and oncogenesis, immune cells such as macrophages, T regulatory (Treg) cells, and myeloid-derived suppressor cells (MDSCs) are constantly activated. Cells that contribute to protective antimicrobial or antitumor immunity, such as natural killer (NK) cells and T effector (Teff) cells, are mostly inhibited . Multiple suppressive cells and molecules function in the TME, thus impairing antitumor responses . Nevertheless, how these cells and molecules change in the dynamic and highly complex microenvironment of a healthy liver that converts to a cancerous liver, remains poorly understood.
Here, we integrated 31 cases of single-cell transcriptomic data from healthy patients and those with hepatitis, cirrhosis, and cancer. Interestingly, we found that macrophage-naïve CD4 + T cell interaction significantly influenced TME, which presented various communication pathways in different states. On the one hand, this interaction became active in cirrhosis and cancer; on the other hand, the differentiation of Treg cells from naïve CD4 + T cell occurred first in cirrhosis and then became more obvious in cancer. Consequently, we hypothesized that naïve CD4 + T cell induced by Macrophage differentiated into Treg in cirrhosis and they finally contributed to immunosuppressive TME. Furthermore, we generated a score to quantify the macrophage-naïve CD4 + T cell interaction and predicted the drugs targeting this interaction. Our results enable to shed light on the key interactions and sources of immunosuppressive TME in the liver, potentially assisting to prevent tumor evasion, and further guide the development of rational immunotherapy for primary liver carcinoma.
Materials and methods
Single-cell transcriptomic data were retrieved from Gene Expression Omnibus (GEO; https://www.ncbi.nlm.nih.gov/geo/) database. Three datasets were employed in our research: GSE159977 , GSE136103 , and GSE125449 . We selected the cases among these three datasets, finally settling on 31 cases (hepatitis = 3, cirrhosis = 5, health = 5, and cancer = 18) for this study. Besides, bulk transcriptomic data were retrieved from The Cancer Genome Atlas (TCGA-LIHC; https://portal.gdc.cancer.gov/), GEO databases (GSE25097; GSE33814; GSE54236) and ICGC database (LIRI cohort; https://dcc.icgc.org/). After data collection, we summarized and presented the workflow of this study (Additional file 1: Fig. S1). Simultaneously, we collected some available clinical features of employed datasets and samples from single-cell transcriptomic data in Additional file 1: Table S1.
Single-cell RNA-seq data processing
We downloaded the three files (barcodes, features, and matrix) for each case and employed them to generate the object of Scanpy . Considering the immune cells employed in one dataset (GSE159977), we only extracted the immune cells for further analyses. Genes expressed in fewer than three cells in a sample were excluded, as were cells that expressed fewer than 200 genes. Further quality control was performed on cells based on the number of genes expressed in the count matrix and the percentage of mitochondrial gene counts. Cells with > 2500 genes were filtered, as well as cells with > 5% of mitochondrial gene count (Additional file 1: Fig. S2A–E). We then applied the library-size correction method to normalize the data matrix using the “scanpy.pp.normalize_total” function in Scanpy. A logarithmized normalized data matrix was employed for downstream analysis.
Dimension reduction and unsupervised clustering
Dimension reduction and unsupervised clustering were performed based on the workflow in Scanpy. We selected highly variable genes for downstream analysis by using the “scanpy.pp.highly_variable_genes” function with parameter “highly_variable_nbatches ≥ 4”. A total of 2303 highly variable genes were identified. To investigate the effect of the cell cycle, we calculated the “S_score” and “G2M_score” by using “sc.tl.score_genes_cell_cycle”. Then, the effects of total counts per cell, the percentage of mitochondrial genes expressed, the “S_score” and “G2M_score” were regressed out using the “scanpy.pp.regress_out” function. We also scaled each gene to unit variance using “scanpy.pp.scale” with parameter “max_value = 10”. After data preprocessing, we reduced the dimensionality of the data by performing a principal component analysis (PCA). A PCA matrix was calculated to reveal the main axes of variation and denoise the data through the “scanpy.tl.pca” function with parameter “svd_solver = ‘arpack’, n_pcs = 11” (Additional file 1: Fig. S2F). To remove the batch effects from different datasets, a BBKNN (Batch Balanced K Nearest Neighbors) algorithm with parameter “n_pcs = 35, metric = ‘euclidean” was performed. BBKNN is a fast and lightweight batch alignment method, which is written in Python and compatible with Scanpy, and its output can be immediately used for dimensionality reduction, clustering and pseudotime inference . Previous study  has compared the different algorithms for batch correction and found that BBKNN showed the good performance in high number of batches and RAM usage. Furthermore, the dimensionality of merged datasets was reduced using uniform manifold approximation and projection (UMAP) implemented by the “scanpy.tl.umap” function. Then, to cluster the neighborhood graph of the cells, we employed the Leiden graph-clustering method. The marker genes of each cluster were identified using the “scanpy.tl.rank_genes_groups” function.
Cell cluster annotation
We employed the marker genes that were validated experimentally and calculated from Scanpy to annotate cell clusters. Genes that were highly and specifically expressed in the cell cluster were considered as marker genes. CellMarker (http://biocc.hrbmu.edu.cn/CellMarker/)  and Panglao DB (https://panglaodb.se/)  websites were used to identify the cell types representing the marker genes.
Cell–cell interaction analyses
We used CellPhone database  and CellChat  to infer cell–cell interactions between immune cell subsets. The potential interaction strength was predicted based on the expression of immune-associated receptors and ligands. The enriched ligand-receptor interactions were calculated based on a permutation test. p < 0.01 were considered significant. To further complement the cell–cell interaction analysis using CellPhoneDB, we employed the CellChat package for analysis. The Cellchat package includes a comprehensive signaling molecule interaction database that considers the known structural composition of receptor-ligand interactions, such as multimeric receptor-ligand complexes, soluble agonists and antagonists, as well as stimulatory and inhibitory membrane-bound coreceptors . The inference of cell–cell interactions includes identification of differentially expressed signaling genes, calculation of ensemble average expression and intercellular communication probability, and identification of statistically significant intercellular communications. We divided the single-cell transcriptomic data into four parts according to the states (health, hepatitis, cirrhosis, and cancer) and explored the differences in ligand and receptor interactions and signaling pathways, across these four states.
Gene regulatory network
To identify cell type-specific gene regulatory networks, single-cell regulatory network inference and clustering (SCENIC) analysis were performed. Genes that were expressed in 3% of samples and cells that expressed > 0 unique molecular identifier (UMI) and normalized by log2(filtered expr + 1) were filtered according to the SCENIC workflow . Genes present in RcisTarget’s human feather databases (hg19-500 bp-upstream-7species.mc9nr.feather and hg19-tss-centered-10 kb-7species.mc9nr.feather) were utilized. These genes are predicted in a region of 500 to 10,000 base pairs, upstream of the hg19 human reference genome. The SCENIC workflow consisted of three steps: (1) identification of co-expression modules between transcription factor (TF) and potential target genes; (2) inference of direct target genes based on those potential targets for which the motif of the corresponding TF is significantly enriched (The regulon is defined as a TF and its direct target genes); (3) calculation of the regulon activity score in every single cell using the area under the recovery curve. To further quantify the cell-type specificity of a regulon, we adapted an entropy-based strategy that was previously used for gene expression data analysis, to generate a regulon-specific score (RSS) . The algorithm was comprehensively described in a previous study .
To characterize the developmental trajectory between naïve CD4 + T cells and Treg cells, trajectory analysis was performed by partition-based graph abstraction (PAGA)  in this research. We employed the PAGA method of Scanpy to assess the most likely trajectories between naïve CD4 + T and Treg cells. The computations were performed using the default parameters.
TCGA data analyses
Score generation and validation
To further validate the role of macrophage-naïve CD4 + T cell interaction in bulk transcriptomic data, we established a score based on the target genes of top regulons from these two cell types. First, differential expression analysis was performed to select the target genes. Genes with both p < 0.05, and |log2fold change|> 1 were considered significantly differentially expressed. Subsequently, we performed univariate Cox regression analysis to further explore the prognostic genes. Genes in the univariate Cox analysis were eligible for further selection if they exhibited a p < 0.01. The least absolute shrinkage and selection operator (LASSO) regression analysis was performed to establish the score. Here, a lasso penalty was used to account for shrinkage and variable selection. The optimal value of the lambda penalty parameter was defined by performing 10 cross-validations. The formula for calculating scores was as follows:
Furthermore, to investigate the correlation between macrophage–naïve CD4 + T cell (MNT) score and overall survival, we performed survival analysis between the high and low-score groups using the “survival” package. We also performed gene set enrichment analysis (GSEA) to explore signaling pathways in the high-score group. Considering the IL16 related pathways activated between macrophages and naïve CD4 + T cells in cancer, we selected the potential signaling pathways in Molecular Signatures Database as input in GSEA. Regarding the immune cell infiltration in bulk transcriptomic data, we employed four algorithms (CIBERSORT, QUANTISEQ, MCP counter, and TIMER) to estimate and compare the difference between the high and low-score groups.
Drug sensitivity data of human cancer cell lines were obtained from the Cancer Therapeutics Response Portal (CTPR, https://portals.broadinstitude.org/ctrp) and PRISM Repurposing dataset (https://depmap.org/portal/prism/). The CTRP and PRISM contain the sensitivity data for 481 compounds over 835 CCLs and 1448 compounds over 482 CCLs, respectively. Both datasets employed the area under the dose–response curve (area under the curve, AUC) values to estimate drug sensitivity. Lower AUC values indicated higher sensitivity to treatment. For drug response prediction, we employed the ridge regression model located in the “pRRophetic” package to measure drug response and compared the difference in drug sensitivity between high and low-score groups. The algorithm of drug prediction was comprehensively described in previous study .
Identification of Macrophage subtypes
Considering the various functions of macrophage subtypes, we investigated the correlation between identified macrophages (single-cell transcriptomic data) and macrophage infiltration (bulk transcriptomic data). We employed the differential expression genes of each state of macrophage to define the macrophage by ssGSEA. Considering the different states of macrophage, three datasets involved health (GSE33814), hepatitis (GSE33814), cirrhosis (GSE25097), and cancer (TCGA-LIHC) were using in this research. In terms of the subtypes of macrophage, we calculated the M1 polarization and M2 polarization score using ssGSEA. The defined genes were presented in Additional file 1: Table S2. Spearman method was used to calculate the correlation between different states of macrophage and polarization score.
Spatial transcriptomic analyses
To further investigate the spatial location of single cell clusters in HCC TME, we employed “CARD” package  to deconvolute the spatial transcriptomic data based on our single cell data of cancer state. The spatial transcriptomic data were retrieved from previous study . We selected four HCC tissue section for deconvolution. Using “CARD” package, we created “CARD” object through “CreateCARDObject” function. And we employed the “CARD_deconvolution” function with default parameters to calculate the results.
To explore the underlying relationship between single cell clusters and clinical characteristics, our study used “Scissor” package  to analyze. Based on bulk transcriptomic data and phenotype information, “Scissor” can automatically identify cell subpopulations from single-cell data that are most responsible for the differences of phenotypes . We employed TCGA-LIHC cohort and select “tumor size” and “stage” to analyze the single cell data of cancer. The “Scissor” function with “alpha = 0.01, family = binomial” parameters was executed.
A multi-sector single-cell atlas of liver
In total, 53,184 cells from 31 cases were analyzed in our study. Leiden clustering of these cells identified 22 distinct major clusters representing epithelial, immune, endothelial, and fibroblast populations (Fig. 1A). Although we selected immune cells in data preprocessing, a fraction of stromal cells was found; this may be attributed to the different clustering methods and annotations. As illustrated in Additional file 1: Fig. S2G–J, the batch effects were removed. Figure 1B shows the distribution of different states in the UMAP. Next, we employed the specific markers to annotate leiden clusters into 21 cell types of the liver (cluster 17 was undefined), comprising T helper cells (CD40LG, CCR6), naïve CD4 + T cells (MAL, CCR7, LEF1), CD8 + exhausted T cells (CD8A, TIGIT, PDCD1), ZNF683 + cytotoxic T lymphocytes (ZNF683, IFNG, TNF), FGFBP2 + cytotoxic T lymphocytes (GZMK, NKG7, FGFBP2), XCL1 + natural killer cells (XCL1, IL2RB, XCL2), GNLY + natural killer cells (GNLY, GZMB, FCGR3A), mucosal-associated invariant T cells (NCR3, SLC4A10), macrophage (TREM2, FCGR1A,GPR34), CD1C + _A dendritic cells (CLEC10A, CD1C, FCER1A), CD1C + _B dendritic cells (S100A9, S100A8, FCN1), CD141-CD1C-dendritic cells (LST1, CFD), B cells (CD79A, MS4A1, CD79B), Plasma cells (IGKC, IGHG3, IGHG1), endothelial cells (GNG11, ID1, ID3), Treg cells (TNFRSF4, CTLA4, FOXP3), AXL + dendritic cells (LGMN, AXL, DAB2), plasmacytoid dendritic cells (IRF7, SERPINF1, LILRA4), CD141 + CLEC9A + dendritic cells (CLEC9A, IDO1, BATF3), Hepatic stellate cell (ALB, IGFBP6, COL1A1), and mast cells (HPGDS, VWA5A, LTC4S) (Fig. 1C and Additional file 1: Fig. S3). To characterize the change in each cell type among the four states, we calculated the proportion of cell types during each state. Interestingly, we observed that some immune cells have a relatively higher proportion in cancer, including T helper cells, naïve CD4 + T cells, exhausted CD8 + T cells, macrophages, B cells, plasma cells, and Treg cells (Fig. 1D). The statistical significance of each cell type in four states were presented in Additional file 1: Fig. S4.
To decipher the ligand-receptor interactions in different states, we performed cell–cell interaction analyses. The results from CellPhoneDB indicated that different states exhibited various immune cell interactions (Fig. 2A–D). Furthermore, we observed that macrophage-naïve CD4 + T cell interaction had an important effect on the cirrhosis and cancer but not the hepatitis. Simultaneously, the results from CellChat revealed that macrophages exhibited higher outgoing interaction strength in cirrhosis and cancer, but lower strength in hepatitis (Fig. 3A). Meanwhile, naïve CD4 + T cells showed a lower outgoing interaction strength in health, hepatitis, and cancer states, but their outgoing strength increased in the cirrhosis state. We speculated that naïve CD4 + T cells may be regulated by macrophages through paracrine mechanisms, mainly in hepatitis, cirrhosis, and cancer. Consequently, we investigated the signaling pathways by which macrophages act as senders and naïve CD4 + T cells as receivers, in subsequent analyses. From the health to hepatitis, the MIF signaling network and related ligand-receptor interactions, including MIF-(CD74 + CXCR4), MIF-(CD74 + CXCR2), and MIF-(CD74 + CD44) showed a significant effect on macrophage–naïve CD4 + T cell interaction (Fig. 3B and E). From hepatitis to cirrhosis, the CCL signaling network and related ligand-receptor interactions, including CCL5-CCR5 and CCL5-CCR1, presented the highest communication probability in macrophage-naïve CD4 + T cell interaction (Fig. 3C and F). From the cirrhosis to cancer, IL16 and related ligand-receptor interaction (IL16-CD4) exhibited a significant difference in macrophage-naïve CD4 + T cell interaction (Fig. 3D and G).
Furthermore, we visualized the network and interactions of the above-mentioned signaling pathways (MIF, CCL, and IL16). In the healthy state, the communication probability exhibited by the MIF signaling pathway was insignificant, similar to that of IL16 in the cirrhosis state. Therefore, the only difference observed was in the CCL signaling pathway, between hepatitis and cirrhosis in our study. As illustrated in Figs. 4A and B, the interactions based on the MIF signaling pathway between macrophages and other immune cells, were mediated mainly by paracrine signaling. From hepatitis to cirrhosis, the total interactions of the CCL signaling pathway increased significantly (Fig. 4C and E), and both the strength of autocrine and paracrine cells from macrophages and naïve CD4 + T cells, increased in the CCL signaling pathway (Fig. 4D and F). From cirrhosis to cancer, macrophages and naïve CD4 + T cells employ autocrine and paracrine mechanisms to communicate in the IL16 signaling pathway (Fig. 4G and H).
Gene regulatory network in macrophages and naïve CD4 + T cells
Considering the importance of macrophage-naïve CD4 + T cell interaction in the TME, we employed SCENIC to decipher the gene regulatory network (GRN) of macrophages and naïve CD4 + T cells, followed by regulon-activity-based hierarchical clustering. The results indicated that cirrhosis and cancer presented similar GRN of macrophages, while the difference was observed between hepatitis and the other three states (Fig. 5A). Naïve CD4 + T cells exhibited distinct GRN in the cancer (Fig. 5F). Furthermore, we employed the RSS to evaluate regulons and ranked them in each state (Fig. 5B–E,G–J). We found that the transcriptomic factor (XBP1) was the top-ranked in both macrophages and naïve CD4 + T cells in cancer (Fig. 5E and J).
Treg cells may derive from naïve CD4 + T cells in cancer
A previous study  revealed that Treg cells in human breast cancer are mainly derived from naïve CD4 + T cells recruited by tumor-associated macrophage-derived CCL18. Our results also demonstrated that Treg cells accounted for a relatively high proportion in the tumor microenvironment (Fig. 1D). We speculated that Treg cells may also be converted from naïve CD4 + T cells in cancer and thus, we performed a trajectory analysis. The results from PAGA validated our speculation (Fig. 5K). Furthermore, we found that naïve CD4 + T cells convert into Treg in specific states (cirrhosis and cancer), which may be related to macrophage.
Validation from bulk transcriptomic data
Macrophage-naïve CD4 + T cell interaction score generation
To further elucidate the relationship between macrophage-naïve CD4 + T cell interaction and the clinic, we established a macrophage-naïve CD4 + T cell (MNT) score. First, we selected the target genes of the top regulon from macrophages (cancer) and naïve CD4 + T cells (cancer), respectively, and visualized them in network (Fig. 6A). A total of 125 target genes were utilized in our study, comprising 78 in macrophages (cancer) and 47 in naïve CD4 + T cells (cancer) (Additional file 1: Table S3). We then performed differential expression analysis of selected genes between the normal and tumor tissue based on TCGA-LIHC dataset. Forty-seven genes were eligible for Cox regression analysis (Fig. 6B). The results from the Cox regression analysis showed that 18 genes were significantly correlated with prognosis (p < 0.01) (Fig. 6C). Subsequently, we employed a LASSO regression model to establish the MNT score. The results from Fig. 6D revealed that the MNT score was generated by seven genes (BOD1, SEC61A1, RHEB, CFL1, PTMA, C1orf109, and E2F5). We divided the cases into high-score and low-score groups based on the median. A significant survival difference (p < 0.01) was found between the high and low-score groups (Fig. 7A). Another two HCC datasets (ICGC-LIRI and GSE54236) also demonstrated that the patients with higher MNT score presented the worst prognosis (Additional file 1: Fig. S5). The cases in the high-score group were significantly enriched in leukocyte-and cytokine-related pathways (Fig. 7B). We also found that the cases in the high-score groups presented higher infiltration of Treg cells and macrophages (Fig. 7C–F).
Drug prediction based on CTPR and PRISM
We predicted the potential drugs based on the CTPR and PRISM databases (Fig. 7G–J) and found that the cases in the high-score group may be more sensitive to GSK461364 (inhibitor of polo-like kinase 1), SB-743921 (inhibitor of kinesin family member 11), dabrafenib, ispinesib, and epothilone-b.
Macrophages in TME tend to exhibit the M2 phenotype
To further explore the function of macrophages in this study, we investigated the correlation between differential expression genes of macrophage and polarization score (M1, and M2). The results demonstrated that macrophage tended to present the M1 phenotype in health while M2 in hepatitis, cirrhosis, and cancer (Fig. 7K–R).
Translational level validation from HPA database
To further validate the association between macrophage, naïve CD4 + T cell and Treg cell in liver carcinoma, we estimated the translational level of their marker genes using immunohistochemistry. Three maker genes were found in liver carcinoma of HPA database (TREM2 and GPR34 from macrophage, PDCD1 from Treg cell), which indicated that macrophages and Treg cells existed in liver carcinoma (Additional file 1: Fig. S6).
Spatial co-location was observed between Macrophage and Treg
The deconvolution results of HCC spatial transcriptomic data indicated that macrophage and Treg existed in HCC tissues (Additional file 1: Fig. S7). Furthermore, we observed that both macrophage and Treg located in the similar position (Additional file 1: Fig.S7), which indicated that there was higher communication probability between macrophage and Treg in cancer.
Macrophage, naïve CD4 + T cell and Treg associated with HCC progression
To further explore the association between cell types and clinical characteristics, we combined the bulk transcriptomic data (TCGA-LIHC) and single-cell transcriptomic data (cancer) to performed scissor analyses. The results showed that 80.4% cells in cancer associated with T3/T4, in which naïve CD4 + T cell accounts for 21.5%, macrophage accounts for 11.5% and Treg accounts for 5% (Additional file 1: Fig. S8). Regarding the stage, 77% cells in cancer associated with stage III/stage IV, including 22.8% of naïve CD4 + T cell, 13.9% of macrophage and 4% of Treg (Additional file 1: Fig. S8).
Primary liver carcinoma commonly arises in damaged liver and is characterized by extensive inflammation and fibrosis. The TME, including immune cells, reacts to liver injury by producing cytokines and components of the extracellular matrix, which promotes angiogenesis and survival of cancer stem cells . Generally, primary liver carcinoma presents an immunosuppressive microenvironment. Nevertheless, the correlation between the immunosuppressive microenvironment and liver injury (hepatitis and cirrhosis) remains unclarified. Here, we generated a comprehensive single-cell atlas of the liver to understand the immune cell interactions among health, hepatitis, cirrhosis, and cancer. We first focused on the proportion of each immune cell type in the four states and found that macrophages account for a relatively high proportion of tumor microenvironments. Simultaneously, macrophages and naïve CD4 + T cells exhibit stronger interactions in cancer, which indicated that macrophage-naïve CD4 + T cell interaction essentially affects the cancerous state. These results led us to investigate the differences of macrophage-naïve CD4 + T cell interaction in the four states. Notably, macrophages present higher outgoing and incoming strength in cirrhosis and cancer but not in hepatitis, possibly due to the decreased number of macrophages in hepatitis. Interestingly, we found that MIF-related ligand-receptor interactions were highly activated in hepatitis. MIF exhibits a protective effect against steatosis, resulting in a decreased quantity of macrophages . The interaction between CCL5 and CCR5 significantly influences macrophage numbers from the state of hepatitis to cirrhosis. The receptor CCR5 induces the recruitment of macrophages ; thus, more macrophages can be observed in cirrhosis. A previous study found that the ligand CCL5 promotes steatosis, inflammation , and early cirrhosis . The ligand-receptor pair of IL16 and CD4 is activated in the interaction of macrophages and naïve CD4 + T cells, and in conversion from cirrhosis to the cancer. IL16 induces the recruitment of macrophages in breast cancer . The above-mentioned results elucidated the differences of macrophage-naïve CD4 + T cell interaction under the four states, which may be helpful in exploring the source of immunosuppressive microenvironment of liver carcinoma.
Next, we predicted the GRNs of macrophages in the four states. Importantly, the GRNs of macrophages between cirrhosis and cancer were similar, indicating that tumor-associated macrophages may stem from a cirrhotic state. Macrophages in hepatitis presented a distinct GRN compared to the other three states, revealing that macrophages may perform different functions in hepatitis. Besides, the identified macrophages tended to exhibit the M2 phenotype, which demonstrates that macrophages have an immunosuppressive effect in cancer. Furthermore, a previous study reported that naïve CD4 + T cells are recruited and converted to Treg cells by macrophages-derived CCL18 in breast cancer . In our trajectory analyses of four states, the phenomenon that naïve CD4 + T cell differentiated into Treg occurred first in cirrhosis and then became more obvious in cancer, which was unobserved in health and hepatitis. Interestingly, the interaction between macrophage and naïve CD4 + T cell exhibited in CellphoneDB is consistent with the findings of trajectory analyses. The consistency further confirms the connection between the macrophage-naive CD4 + T cell interaction and Treg population. In summary, macrophages produced in response to cirrhosis play a crucial role in constructing and enhancing the immunosuppressive microenvironment of primary liver carcinoma.
Considering the limited therapy available for primary liver carcinoma, we predicted potential sensitive drugs based on macrophage-naïve CD4 + T cell interaction. First, we generated the MNT score to quantify macrophage-naïve CD4 + T cell interaction. For all cases, high MNT scores led to a worse prognosis, validating the importance of macrophage-naïve CD4 + T cell interaction in bulk transcriptomic data. Meanwhile, the cases with high MNT scores presented higher infiltration of Treg cells and macrophages, which corresponds with the results of single-cell analyses. Regarding drug prediction, inhibitors of polo-like kinase 1 and kinesin 11 based on CTPR, were significantly correlated with high MNT scores. Polo-like kinase 1 has been shown to regulate liver tumor cell death by phosphorylation of Tap63 . Meanwhile, kinesin family member 11 has been shown to correlate with tumor size and prognosis in liver cancer . Three other drugs from PRISM also cause concern. Dabrafenib is a single-agent treatment for patients with BRAF V600E mutation-positive advanced melanoma . Epothilone-b is an anti-cancer drug for treating aggressive metastatic or locally advanced breast cancer, and prevents cancer cells from dividing by interfering with tubulin . Interestingly, ispinesib is also an inhibitor of KIF11, which has been used to treat patients with locally advanced, recurrent, or metastatic liver cancer (NCT00095992). We reasoned that employing the predictive drugs may target the interaction of macrophages and CD4 naïve T cells, which results in reprogramming of the TME and strengthening of immunosurveillance in primary liver carcinoma.
Compared with the previous studies [6, 10, 11], our research integrated the single cell transcriptomic data of four states from healthy liver to primary liver carcinoma, which was comprehensive and creative. By comparing the TME of four states, we found that different states presented the vary proportion and communication among immune cells. The key point we found in this research was the interaction of macrophage and naïve CD4 + T cell. Macrophage induced naïve CD4 + T cell to differentiate into Treg in cirrhosis and cancer. Simultaneously, we observed that Macrophage and Treg co-located in HCC tissue. These findings may be helpful in decoding the immunosuppressive TME of HCC under single cell resolution. Furthermore, we linked the single cell transcriptomic and clinical phenotypes and verified the clinical importance of Macrophage. At present, macrophage-related therapy in HCC mainly divided into three approaches, including cutting off the source and eliminating the production of M2 macrophages, remodeling M2 macrophages to M1 macrophages, and blocking communication between M2 macrophages and liver cancer cells . Blocking the interaction of macrophage and naïve CD4 + T cell may be the novel and interesting therapeutic strategy in HCC immunotherapy.
In conclusion, this study reveals the crucial role of macrophage-naïve CD4 + T cell interaction in the immunosuppressive microenvironment of primary liver carcinoma. Tumor-associated macrophages may derive from cirrhosis. Naïve CD4 + T cell induced by macrophage may differentiate into Treg in cirrhosis and they finally contributed to immunosuppressive TME. Predictive drugs that target the macrophage-naïve CD4 + T cell interaction may help to improve the immunosuppressive microenvironment and prevent immune evasion. However, due to the bioinformatics used in our research, the relevant mechanisms need to further validations by experiments and cohort studies.
Availability of data and materials
The datasets generated during the current study are available in the TCGA database (https://portal.gdc.cancer.gov/), GEO database (https://www.ncbi.nlm.nih.gov/geo/) and ICGC database (https://dcc.icgc.org/).
Villanueva A. Hepatocellular carcinoma. New Engl J Med. 2019;380:1450–62. https://doi.org/10.1056/NEJMra1713263.
Llovet JM, Ricci S, Mazzaferro V, Hilgard P, Gane E, Blanc JF, et al. Sorafenib in advanced hepatocellular carcinoma. New Engl J Med. 2008;359:378–90. https://doi.org/10.1056/NEJMoa0708857.
El-Khoueiry AB, Sangro B, Yau T, Crocenzi TS, Kudo M, Hsu C, et al. Nivolumab in patients with advanced hepatocellular carcinoma (CheckMate 040): an open-label, non-comparative, phase 1/2 dose escalation and expansion trial. Lancet. 2017;389:2492–502. https://doi.org/10.1016/S0140-6736(17)31046-2.
Finn RS, Qin S, Ikeda M, Galle PR, Ducreux M, Kim TY, et al. Atezolizumab plus bevacizumab in unresectable hepatocellular carcinoma. New Engl J Med. 2020;382:1894–905. https://doi.org/10.1056/NEJMoa1915745.
Zhu AX, Finn RS, Edeline J, Cattan S, Ogasawara S, Palmer D, et al. Pembrolizumab in patients with advanced hepatocellular carcinoma previously treated with sorafenib (KEYNOTE-224): a non-randomised, open-label phase 2 trial. Lancet Oncol. 2018;19:940–52. https://doi.org/10.1016/S1470-2045(18)30351-6.
Pfister D, Núñez NG, Pinyol R, Govaere O, Pinter M, Szydlowska M, et al. NASH limits anti-tumour surveillance in immunotherapy-treated HCC. Nature. 2021;592:450–6. https://doi.org/10.1038/s41586-021-03362-0.
Fattovich G, Stroffolini T, Zagni I, Donato F. Hepatocellular carcinoma in cirrhosis: incidence and risk factors. Gastroenterology. 2004;127:S35-50. https://doi.org/10.1053/j.gastro.2004.09.014.
Keenan BP, Fong L, Kelley RK. Immunotherapy in hepatocellular carcinoma: the complex interface between inflammation, fibrosis, and the immune response. J Immunother Cancer. 2019;7:267. https://doi.org/10.1186/s40425-019-0749-z.
Palucka AK, Coussens LM. The Basis of oncoimmunology. Cell. 2016;164:1233–47. https://doi.org/10.1016/j.cell.2016.01.049.
Ramachandran P, Dobie R, Wilson-Kanamori JR, Dora EF, Henderson BEP, Luu NT, et al. Resolving the fibrotic niche of human liver cirrhosis at single-cell level. Nature. 2019;575:512–8. https://doi.org/10.1038/s41586-019-1631-3.
Ma L, Hernandez MO, Zhao Y, Mehta M, Tran B, Kelly M, et al. Tumor cell biodiversity drives microenvironmental reprogramming in liver cancer. Cancer Cell. 2019;36(4):418-430.e6. https://doi.org/10.1016/j.ccell.2019.08.007.
Wolf FA, Angerer P, Theis FJ. SCANPY: large-scale single-cell gene expression data analysis. Genome Biol. 2018;19:15. https://doi.org/10.1186/s13059-017-1382-0.
Polański K, Young MD, Miao Z, Meyer KB, Teichmann SA, Park JE. BBKNN: fast batch alignment of single cell transcriptomes. Bioinformatics. 2020;36:964–5. https://doi.org/10.1093/bioinformatics/btz625.
Chazarra-Gil R, van Dongen S, Kiselev VY, Hemberg M. Flexible comparison of batch correction methods for single-cell RNA-seq using BatchBench. Nucleic Acids Res. 2021;49: e42. https://doi.org/10.1093/nar/gkab004.
Zhang X, Lan Y, Xu J, Quan F, Zhao E, Deng C, et al. Cell Marker: a manually curated resource of cell markers in human and mouse. Nucleic Acids Res. 2019;47:D721–8. https://doi.org/10.1093/nar/gky900.
Franzén O, Gan LM, Björkegren JLM. PanglaoDB: a web server for exploration of mouse and human single-cell RNA sequencing data. Database (Oxford). 2019. https://doi.org/10.1093/database/baz046.
Efremova M, Vento-Tormo M, Teichmann SA, Vento-Tormo R. Cell PhoneDB: inferring cell–cell communication from combined expression of multi-subunit ligand–receptor complexes. Nat Protoc. 2020;15:1484–506. https://doi.org/10.1038/s41596-020-0292-x.
Jin S, Guerrero-Juarez CF, Zhang L, Chang I, Ramos R, Kuan CH, et al. Inference and analysis of cell-cell communication using Cell Chat. Nat Commun. 2021;12:1088. https://doi.org/10.1038/s41467-021-21246-9.
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. https://doi.org/10.1038/nmeth.4463.
Cabili MN, Trapnell C, Goff L, Koziol M, Tazon-Vega B, Regev A, et al. Integrative annotation of human large intergenic noncoding RNAs reveals global properties and specific subclasses. Gene Dev. 2011;25:1915–27. https://doi.org/10.1101/gad.17446611.
Suo S, Zhu Q, Saadatpour A, Fei L, Guo G, Yuan GC. Revealing the critical regulators of cell identity in the mouse cell atlas. Cell Rep. 2018;25:1436-45.e3. https://doi.org/10.1016/j.celrep.2018.10.045.
Wolf FA, Hamey FK, Plass M, Solana J, Dahlin JS, Göttgens B, Rajewsky N, Simon L, Theis FJ. PAGA: graph abstraction reconciles clustering with trajectory inference through a topology preserving map of single cells. Genome Biol. 2019;20:59. https://doi.org/10.1186/s13059-019-1663-x.
Yang C, Huang X, Li Y, Chen J, Lv Y, Dai S. Prognosis and personalized treatment prediction in TP53-mutant hepatocellular carcinoma: an in silico strategy towards precision oncology. Brief Bioinform. 2020. https://doi.org/10.1093/bib/bbaa164.
Ma Y, Zhou X. Spatially informed cell-type deconvolution for spatial transcriptomics. Nat Biotechnol. 2022;40:1349–59. https://doi.org/10.1038/s41587-022-01273-7.
Wu R, Guo W, Qiu X, Wang S, Sui C, Lian Q, et al. Comprehensive analysis of spatial architecture in primary liver cancer. Sci Adv. 2021. https://doi.org/10.1126/sciadv.abg3750.
Sun D, Guan X, Moran AE, Wu LY, Qian DZ, Schedin P, et al. Identifying phenotype-associated subpopulations by integrating bulk and single-cell sequencing data. Nat Biotechnol. 2022;40:527–38. https://doi.org/10.1038/s41587-021-01091-3.
Su S, Liao J, Liu J, Huang D, He C, Chen F, et al. Blocking the recruitment of naive CD4(+) T cells reverses immunosuppression in breast cancer. Cell Res. 2017;27:461–82. https://doi.org/10.1038/cr.2017.34.
Hernandez-Gea V, Toffanin S, Friedman SL, Llovet JM. Role of the microenvironment in the pathogenesis and treatment of hepatocellular carcinoma. Gastroenterology. 2013;144:512–27. https://doi.org/10.1053/j.gastro.2013.01.002.
Heinrichs D, Berres M-L, Coeuru M, Knauel M, Nellen A, Fischer P, et al. Protective role of macrophage migration inhibitory factor in nonalcoholic steatohepatitis. FASEB J. 2014;28:5136–47. https://doi.org/10.1096/fj.14-256776.
Barashi N, Weiss ID, Wald O, Wald H, Beider K, Abraham M, et al. Inflammation-induced hepatocellular carcinoma is dependent on CCR5 in mice. Hepatology. 2013;58:1021–30. https://doi.org/10.1002/hep.26403.
Kim B-M, Abdelfattah AM, Vasan R, Fuchs BC, Choi MY. Hepatic stellate cells secrete Ccl5 to induce hepatocyte steatosis. Sci Rep. 2018;8:7499. https://doi.org/10.1038/s41598-018-25699-9.
Li B, He F, Yang X, Chen YW, Fan JG. Steatosis induced CCL5 contributes to early-stage liver fibrosis in nonalcoholic fatty liver disease progress. Transl Res. 2017;180:103–17. https://doi.org/10.1016/j.trsl.2016.08.006.
Richmond J, Tuzova M, Cruikshank W, Center D. Regulation of cellular processes by interleukin-16 in homeostasis and cancer. J Cell Physiol. 2014;229:139–47. https://doi.org/10.1002/jcp.24441.
Komatsu S, Takenobu H, Ozaki T, Ando K, Koida N, Suenaga Y, et al. Plk1 regulates liver tumor cell death by phosphorylation of TAp63. Oncogene. 2009;28:3631–41. https://doi.org/10.1038/onc.2009.216.
Hu ZD, Jiang Y, Sun HM, Wang JW, Zhai LL, Yin ZQ, et al. KIF11 Promotes proliferation of hepatocellular carcinoma among patients with liver cancers. Biomed Res Int. 2021;2021:2676745. https://doi.org/10.1155/2021/2676745.
Puszkiel A, Noé G, Bellesoeur A, Kramkimel N, Paludetto MN, Thomas-Schoemann A, et al. Clinical pharmacokinetics and pharmacodynamics of dabrafenib. Clin Pharmacokinet. 2019;58:451–67. https://doi.org/10.1007/s40262-018-0703-0.
Cheng H, Huang G. Synthesis & antitumor activity of epothilones B and D and their analogs. Future Med Chem. 2018;10:1483–96. https://doi.org/10.4155/fmc-2017-0320.
Cheng K, Cai N, Zhu J, Yang X, Liang H, Zhang W. Tumor-associated macrophages in liver cancer: from mechanisms to therapy. Cancer Commun (Lond). 2022. https://doi.org/10.1002/cac2.12345.
This study was supported by the National Natural Science Foundations of China (NO. 81873248, NO. 82174173, NO. 81972785, NO. 81773162, NO. 81903967, NO. 82104962).
Ethics approval and consent to participate
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Additional file 1: Figure S1.
The workflow of study. Figure S2. Three violin plots of the computed quality measures: (A) the number of genes expressed in the count matrix, (B) the total counts per cell, and (C) the percentage of counts in mitochondrial genes. Two scatter plots: (D) the number of genes expressed in the count matrix and the total counts per cell, (E) the percentage of counts in mitochondrial genes and the total counts per cell. (F) The result of principal component analysis. This gives us information about how many PCs we should consider in order to compute the neighborhood relations of cells. The umap of four states before (G) and after (H) batch effect removing. The umap of three datasets before (I) and after (J) batch effect removing. Figure S3. The marker genes presented in umap. Figure S4. Statistical significance between four states in each cell type. Note: each circle represents the sample. Figure S5. Validation of survival analyses based on MNT score from ICGC-LIRI and GSE54236 cohorts. Figure S6. The translational level of TREM2, GPR34, and PDCD1 in immunohistochemistry. Figure S7. Deconvolution of spatial transcriptomic data. Note: The higher value represents the higher probability of the target cell type location. Figure S8. The association between cells from cancer and clinical characteristics (tumor size and stage). The left umaps were the results of Scissor analyses. The right bar plots were the proportion of each cell type in Scissor+ cell group. Table S1. The details of three public datasets. Table S2. The representative genes of M1 and M2 polarization. Table S3. The transcription factors of top 1 regulons from macrophage and naïve CD4+ T cells in cancer.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.
About this article
Cite this article
Mo, Z., Liu, D., Chen, Y. et al. Single-cell transcriptomics reveals the role of Macrophage-Naïve CD4 + T cell interaction in the immunosuppressive microenvironment of primary liver carcinoma. J Transl Med 20, 466 (2022). https://doi.org/10.1186/s12967-022-03675-2
- Single cell transcriptomic
- Naïve CD4 + T cell
- Tumor microenvironment