Zongbo Dai1†
Chengshuo Zhang- 1Hepabobiliary Surgery Department, First Hospital of China Medical University, Shenyang, China
- 2Department of General Surgery, Anshan Central Hospital, Anshan, China
Background: The heterogeneity of the tumor microenvironment significantly influences the prognosis of hepatocellular carcinoma (HCC) patients, with cell communication through ligand-receptor complexes playing a central role.
Methods: We conducted single-cell transcriptomic analysis on ten HCC tissues to identify ligand-receptor genes involved in malignant HCC cell communication using CellChat. Leveraging RNA-Seq data from the TCGA Liver Cancer (TCGA-LIHC) and Liver Cancer - RIKEN, JP (LIRI-JP) cohorts, we employed Cox regression analysis to screen for prognosis-related genes. Prognostic risk models were constructed through unsupervised clustering and differential gene expression analysis. Subsequently, a co-culture system involving tumor cells and macrophages was established. A series of experiments, including Transwell assays, immunofluorescence staining, immunoprecipitation, flow cytometry, and immunohistochemistry, were conducted to elucidate the mechanism through which HCC cells recruit macrophages via the CCL16-CCR1 axis.
Results: Single-cell analysis unveiled significant interactions between malignant HCC cells and macrophages, identifying 76 related ligand-receptor genes. Patients were classified into three subtypes based on the expression patterns of eight prognosis-related ligand-receptor genes. The subtype with the worst prognosis exhibited reduced infiltration of T cell-related immune cells, downregulation of immune checkpoint genes, and increased M2-like tumor-associated macrophage scores. In vitro experiments confirmed the pivotal role of the CCL16-CCR1 axis in the recruitment and M2 polarization of tumor-associated macrophages. Clinical samples demonstrated a significant association between CCL16 protein expression levels and advanced stage, lymph node metastasis, and distant metastasis. Immunohistochemistry and immunofluorescence staining further confirmed the correlation between CCL16 and CCR1, CD68, and CD206, as well as CD68+CCR1+ macrophage infiltration.
Conclusions: Our study identified molecular subtypes, a prognostic model, and immune microenvironment features based on ligand-receptor interactions in malignant HCC cell communication. Moreover, we revealed the pro-tumorigenic role of HCC cells in recruiting M2-like tumor-associated macrophages through the CCL16-CCR1 axis.
Graphical Abstract
Introduction
Hepatocellular carcinoma (HCC), a prevalent and aggressive cancer, poses a substantial global burden (1). Despite advancements in HCC diagnosis and treatment, patients with this condition continue to face an unfavorable prognosis. This underscores the urgent need for a more comprehensive understanding of the fundamental mechanisms driving tumor progression and the identification of novel treatment targets (2). The progression of HCC is profoundly influenced by the tumor microenvironment (3). This ecosystem is intricate, encompassing various cell types, including cancer cells, immune cells, stromal cells, and components of the extracellular matrix. Recognizing the diversity of the tumor microenvironment has gained increasing importance as a critical determinant of HCC prognosis (4). Unraveling the molecular mechanisms underpinning this heterogeneity and its association with patient subtyping and prognosis is essential for advancing HCC management.
Macrophages, among the diverse cell populations within the tumor microenvironment, play a pivotal role in regulating both the immune response and tumor behavior (5). In HCC, the interaction between macrophages and tumor cells holds particular significance, influencing tumor growth, invasion, metastasis, and immune evasion (6, 7). Notably, Weng et al. identified that patients with a significant presence of PPT1+ macrophages experienced heightened infiltration of CD8+ T cells and increased expression of programmed cell death-1 (PD-1) (8). Additionally, Zhao et al. demonstrated that CD168+ M2-like macrophages regulate TOP2A+CENPF+ liver cancer stem-like cells, contributing to the excessive growth of HCC (9). Furthermore, Wang et al. reported that M2 tumor-associated macrophages (TAMs) enhance stemness levels and reduce sorafenib-induced cell apoptosis through paracrine secretion of CXCL1 and CXCL2 (10). A comprehensive exploration of the intricate interaction between macrophages and cancer cells provides valuable insights into HCC development and guides the development of innovative treatment approaches (11). In this dynamic microenvironment, cell communication mediated by ligand-receptor complexes orchestrates a broad spectrum of cellular responses, significantly influencing tumor behavior (12). However, the specific ligand-receptor interactions operating within the HCC microenvironment and their functional implications remain largely unexplored.
CCL16, a chemokine renowned for its chemotactic properties, has attracted attention due to its potential involvement in HCC progression and immune modulation. While Shen et al. demonstrated its role in promoting the stemness of breast cancer cells (13), its specific implications in HCC remain unexplored. Notably, CCL16 has exhibited the ability to engage with various chemokine receptors, specifically CCR1, CCR2, CCR5, and CCR8 (14–17), expressed on immune cells, including macrophages. This interaction suggests that CCL16 could influence the tumor microenvironment by recruiting immune cells, such as macrophages, through the activation of these receptors. The presence of CCL16-related interactions may further modulate the polarization and functional properties of macrophages, potentially influencing their pro- or anti-tumorigenic activities within the HCC microenvironment. Despite these intriguing observations, the specific roles and underlying mechanisms of CCL16 in the context of HCC remain incompletely understood, warranting further dedicated investigation.
This research aimed to investigate the intricate landscape of ligand-receptor interactions within the HCC microenvironment and elucidate their functional consequences. Utilizing cutting-edge methodologies, including single-cell transcriptomic analysis and comprehensive cohort data, our objective was to characterize the complex network of ligand-receptor interactions in HCC. By employing a combination of computational analyses, experimental validations, and clinical correlations, we sought to uncover valuable insights into the role of these interactions in tumor progression and immune modulation. This comprehensive exploration aimed to contribute to our understanding of the pathogenesis of HCC, with the ultimate goal of identifying potential therapeutic targets for precision medicine approaches in the treatment of HCC.
Methods
Source of data and preprocessing
The count matrix for the dataset GSE149614 based on the single-cell RNA-seq technology was extracted from the GEO database (https://ncbi.nlm.nih.gov/geo/). The TCGA Liver Cancer (TCGA-LIHC) dataset’s Level 3 Bulk RNA-Seq expression data were acquired from UCSC Xena (https://xenabrowser.net/datapages/) (18), whereas the ICGC database (https://dcc.icgc.org/) provided the Level 3 data for the Liver Cancer - RIKEN, JP (LIRI-JP) dataset. For further analysis, the Bulk RNA-seq data at Level 3 underwent a log2(x+1) transformation.
Analysis of data at the single-cell level
The single-cell data obtained from 10 liver cancer samples in GSE149614 underwent quality control and filtering with the Seurat R package (19, 20). The specific requirements were established as follows: nFeature_RNA must be between 500 and 4000, percent.mt should be less than 5, and percent.HB should be less than 1.Next, the data underwent normalization using the LogNormalize technique via the NormalizeData function. The FindVariableFeatures function was used to identify 2000 genes that exhibit high variability in each sample. Afterwards, the CCA method was utilized to eliminate batch discrepancies and merge the information. Principal Component Analysis was performed using the RunPCA function, with the parameters dims set to 20 and resolution set to 0.5. The identification of cell clusters was done using the FindNeighbors and FindClusters functions. The FindMarkers function was used to differentiate cell populations and identify marker genes for each cluster. Using default parameters, the CopyKAT R package was utilized to estimate Copy Number Variations (CNV) for every cell, enabling the detection of aneuploid cells (21). In the end, we utilized the CellChat R package (22) to deduce and examine the cellular communication network and discover ligand-receptor pairs linked to aneuploid cells.
Identification of ligand-receptor-related subtypes
To further select prognostic-related variables, univariate Cox regression analysis was conducted in the TCGA-LIHC cohort, based on the expression profiles of ligand-receptor genes linked to aneuploid cells. Afterwards, the ConsensusClusterPlus R package (23) was applied to perform separate unsupervised clustering in the TCGA-LIH (N = 365) and LIRI-JP (N = 232) cohorts. The clusterAlg parameter utilized was pam, with the Euclidean distance, 500 reps, a pItem of 0.8, and pFeature of 1. The log-rank test from the survival package was utilized to conduct survival analysis among various subtypes, resulting in the generation of Kaplan-Meier curves.
Analysis of the tumor microenvironment using quantitative methods
Using the GSVA package, we assessed the levels of 28 immune cell types’ infiltration in the TCGA-LIHC cohort through the single-sample Gene Set Enrichment Analysis (ssGSEA) technique. The background gene set was referenced from the previous study (24). The Stromal score and Immune Score for each sample were calculated using the Estimate R package. We examined the variation in expression of immune checkpoint genes among different groups by utilizing a gene list obtained from a prior study (25). Afterwards, the gene expression patterns were normalized using the scale method, and the Tumor Immune Dysfunction and Exclusion (TIDE) score for each specimen was calculated using the TIDE online tool (http://tide.dfci.harvard.edu/) (26, 27). The ssGSEA method mentioned earlier was used for quantification of all gene sets in this study, and Wilcoxon test was applied for differential analysis.
Detection of variably expressed genes across subtypes
Using the eBayes test method from the limma R package, we performed a comparative analysis of gene expression profiles across various subtypes. For each subtype, genes that exhibited a |FoldChange| > 1.5 and fdr < 0.05 were identified as differentially expressed. Afterwards, the clusterProfiler R package was utilized to perform functional enrichment analysis (28).
Construction of the prognostic model and nomogram
We initially conducted univariate Cox analysis to select initial variables with a threshold of P < 0.05, based on the distinct genes among various molecular subtypes. Afterwards, we employed the glmnet R software package to conduct Lasso regression analysis, specifically in the TCGA-LIHC group. To identify variables corresponding to lambda.min, Cox regression was performed using 10-fold cross-validation. Next, a stepwise regression analysis was utilized to construct the model for computing the risk score of every individual. Validation of the identical model was conducted in the LIRI-JP group. Patients were classified into high- and low-risk groups based on the cohort’s median risk score as the threshold. The survival package was utilized for conducting survival analysis, while the timeROC R package was employed to generate receiver operating characteristic (ROC) curves for assessing the predictive performance based on the Area Under Curve (AUC). Furthermore, the Wilcoxon test was employed to examine the disparities in the immune microenvironment among the two cohorts, while the Pearson correlation test was utilized to explore the association between the risk score and both the immune score and TIDE score. Prognostic impact of the risk score and clinical variables was observed in the TCGA-LIHC cohort through univariate and multivariate Cox analyses. The rms R package was utilized to construct a nomogram by integrating important variables. Afterwards, AUC curves were generated by utilizing the plotAUCcurve function from the timeROC R package, calibration curves were produced using the calibrate function, and the decision curve was plotted with the assistance of the rmda R package.
Cell culture
The Chinese Academy of Sciences Cell Bank provided HEK293T cells and liver cancer cell lines HEPG2/JHH7/HUH7/SNU761. DMEM medium supplemented with 10% fetal bovine serum (FBS) was used to culture HEPG2 and HUH7 cells, whereas RPMI1640 medium supplemented with 10% FBS was used to culture JHH7 and SNU761 cells. The THP1 cells, acquired from our lab, were grown in RPMI1640 solution with 10% fetal bovine serum that had been heat-inactivated, along with 10 U/mL penicillin and 10μg/mL streptomycin. The cells were cultured in a flask with a cell density of around 1×108, in an incubator at a temperature of 37°C and with a CO2 concentration of 5%. Cultured THP-1 cells were thinned to a concentration of 1 million cells per milliliter and then placed in culture dishes with a diameter of 35 millimeters. Afterwards, the cells were prompted to undergo differentiation by being incubated in RPMI1640 solution with a concentration of 100 ng/ml of phorbol 12-myristate 13-acetate (PMA) for a duration of 48 hours. Summary of relevant reagent information in the experimental section of this study can be found in Table 1.
Table 1 Reagent information used in this study.
Cell line knockdown and overexpression
To generate HEPG2 cells with CCL16 knockdown and THP1 cells with CCR1 knockdown, validated shRNA sequences targeting the respective proteins were searched online at the Sigma-Aldrich website (https://www.sigmainformatics.com/Informatics_tools/batch-search.php). Primers were designed based on the selected sequences and cloned into the plKO.1 vector. Subsequently, the plKO.1 vector was co-transfected with the lentiviral packaging vectors psPAX2 and pMD2 into HEK293T cells. After 48 hours, the viral supernatant was collected and used to infect the target cell lines. The infected cells were then selected using 4μg/mL puromycin, and the knockdown efficiency of the target proteins was assessed by qPCR. The shRNA sequences for CCL16 knockdown were as follows: shCCL16-1: TCCTTATCATTACTTCGGCTT, shCCL16-2: AGGAGAAGTATTTCGAATATT. The shRNA sequences for CCR1 knockdown were as follows: shCCR1-1: CCCTACAATTTGACTATACTT, shCCR1-2: CCCTGGTAGAAAGAAGATGAA. To generate SNU761 cells overexpressing CCL16, a Plvx-IRES-ZsGreen1-flag-CCL16 overexpression plasmid was constructed. Subsequently, the Plvx-IRES-ZsGreen1 vector was co-transfected with the lentiviral packaging vectors psPAX2 and pMD2 into HEK293T cells. After 48 hours, the viral supernatant was collected and used to infect SNU761 cells. The infected cells were then selected using 4μg/mL puromycin, and the efficiency of CCL16 overexpression was assessed by qPCR.
Macrophages migration and co-culture assays using Transwell
For the migration assay of macrophages, a 24-well Transwell chamber system with an 8μm membrane was used. PMA-induced THP1 cells (5×104 cells/well) were placed in the upper chamber, while the conditioned medium from the control, CCL16-knockdown HEPG2 cells/CCL16-overexpressing SNU761 cells was placed in the lower chamber. Following a 24-hour period, the THP1 cells that moved towards the bottom chamber were dyed using 0.1% crystal violet. To calculate the average, three fields were randomly chosen under a microscope to determine the number of migrated cells. In the co-culture assay of macrophages and liver cancer cells, a 6-well Transwell chamber system with a 0.4μm membrane was used to separate the upper and lower chambers. Control, CCL16-knockdown HEPG2 cells/CCL16-overexpressing SNU761 cells (1×105 cells/well) were placed in the lower chamber, while PMA-induced THP1 cells (1×105 cells/well) were placed in the upper chamber. After 48 hours of co-culture, the THP1 cells in the upper chamber were subjected to flow cytometry. Three independent replicates were conducted.
Elisa
CCL16 expression in the control, CCL16-knockdown HEPG2 cells, and CCL16-overexpressing SNU761 cells were measured using human CCL16 ELISA kits (Abcam). The particular experimental protocols were carried out in accordance with the provided guidelines. Three independent replicates were conducted.
Immunofluorescence
The SNU761 cells overexpressing CCL16 were co-cultured with PMA-treated THP1 cells. Afterward, the cells were fixed with 4% paraformaldehyde for 10 minutes, blocked with PBS containing 2% BSA for 30 minutes, and then incubated overnight at 4°C with antibodies against Flag and CCR1. On the following day, the cells were stained with corresponding fluorescent secondary antibodies for 1 hour at room temperature. The cell nuclei were stained with DAPI. Afterward, the slides were mounted using an anti-fade mounting medium and observed under a ZEISS microscope for fluorescence imaging. Liver cancer tissue samples were embedded in paraffin and prepared as sections. The sections underwent a series of steps including deparaffinization, hydration, and antigen retrieval. After blocking with 5% goat serum, the sections were incubated overnight at 4°C with appropriately diluted CD68 and CCR1 antibodies. On the next day, the sections were stained with corresponding fluorescent secondary antibodies for 1 hour at room temperature. Three distinct fields were randomly selected from each slide, and the number of CD68+CCR1+ cells was quantified. Three independent replicates were conducted. The average value of the cell counts from the three fields was used for statistical analysis.
Coimmunoprecipitation
Flag-CCL16 recombinant protein was added into THP1 cells. Following a 24-hour period, the cells were disrupted by employing RIPA buffer that consisted of proteinase inhibitors. Ten percent of the entire cell lysate was retained as the control sample, whereas the rest of the lysate was subjected to an overnight incubation at 4°C with Flag and IgG antibody. On the following day, beads containing protein A/G were introduced and left to incubate for an extra 4 hours. Afterwards, the beads were rinsed thrice using RIPA solution and heated in 20μL of 1x loading buffer. Flag/CCR1/CCR2/CCR5/CCR8 antibodies were utilized for immunoblotting analysis to identify the interacting proteins, in addition to the input protein samples. Three independent replicates were conducted.
Western blot
After undergoing SDS-PAGE gel electrophoresis, the proteins were then transferred onto a nitrocellulose (NC) membrane. Specific primary antibodies were incubated with the membrane overnight at 4°C. Following this, the sample was exposed to secondary antibodies for a duration of 1 hour at room temperature. The Immobilon Western chemiluminescent HRP substrate (Millipore) was utilized to visualize protein bands on a luminescent imaging workstation (Tanon). ImageJ software was used to measure the band intensities of various proteins, with ACTIN serving as the chosen loading control.
qPCR
TRIzol reagent was used to extract total RNA according to the instructions provided by the manufacturer. The PrimeScript RT Reagent Kit (Takara) was utilized for the synthesis of cDNA. The SYBR Green Master mix (Roche) was utilized to perform quantitative real-time PCR (QPCR). The 2-ΔΔCt method was utilized to determine the comparative levels of gene expression, which were then normalized to the expression of ACTIN. Table 2 contains a list of the specific primers utilized for gene-specific PCR.
Table 2 Primer sequences for qPCR.
Flow cytometry
After 48 hours of co-culture between macrophages and liver cancer cells in a 0.4μm pore size Transwell system, the PMA induced THP1 cells from the upper chamber were subjected to flow cytometry analysis. Cells were incubated with CD80-FITC and CD206-PE Cyanine7 flow cytometry antibodies at 4°C for 30 minutes. Following the addition of 1 mL PBS and centrifugation at 1500 rpm for 5 minutes, the cell pellet was resuspended in 200μL PBS. Flow cytometry analysis was performed using the Beckman Gallios flow cytometer, and the experimental results were analyzed using Flow Jo. Three independent replicates were conducted.
Immunohistochemistry
Liver cancer tissues were embedded in paraffin and then cut into sections. After deparaffinization and hydration, the sections were subjected to antigen retrieval. The activity of endogenous peroxidase was inhibited using 3% H2O2, followed by blocking the sections with 5% goat serum. Antibodies against CCL16, CD68, and CD206 were diluted, added, and incubated overnight at 4°C. Following the PBST (PBS with Tween 20) wash, the sections were exposed to biotinylated secondary antibodies and left at ambient temperature for a duration of 20 minutes. After being washed three times with PBST, the sections were then exposed to streptavidin-horseradish peroxidase at room temperature for a duration of 15 minutes. After three washes with PBST, fresh DAB working solution was added for color development. To halt the reaction, the sections were rinsed with water for a duration of 5 minutes, and subsequently examined under a ZEISS microscope. The Image Pro Plus software was used to calculate the Mean Optical Density (MOD) value.
Ethics statement
The research was carried out with the endorsement of the Ethics Committee at the First Hospital of China Medical University. All participants provided written consent prior to the study.
Statistical analysis
GraphPad Prism 8 was utilized for conducting statistical analysis. Means ± standard deviations (SD) were used to present quantitative data, whereas numbers were used to present qualitative values. For the analysis of suitable quantitative data, either the Unpaired Student’s t-test or two-way ANOVA was utilized. Statistical significance was determined by a two-sided P-value less than 0.05, indicated as * for P < 0.05, ** for P < 0.01, and *** for P < 0.001.
Results
The microenvironment of malignant liver cancer cells exhibits cellular communication characteristics, with macrophages showing the most notable interaction
To examine the communication properties of cancerous liver cells within the tumor microenvironment, we conducted an analysis of single-cell transcriptomic data using tissue samples obtained from 10 individuals diagnosed with malignant liver cancer. After quality control and clustering, we obtained a total of 13 cell types, comprising 44,320 cells (Figure 1A), including 3,648 aneuploid cells. In addition, we showcased the leading 5 marker genes for every cell cluster (Figure 1B). Among all cell types, T/NK cells had the highest proportion, followed by Macrophages, aneuploid cells, and Monocytes. On the other hand, there was significant heterogeneity in the composition of cell types between different patients (Figure 1C). The tSNE distribution of aneuploid cells identified by CopyKAT primarily overlapped with Hepatocytes and Cycling cells, demonstrating high accuracy (Figure 1D). The interactions between cell clusters identified by Cellchat revealed that the most prominent interacting cell types with aneuploid cells were macrophages and monocytes, followed by fibroblasts and myeloid DC cells (Figure 2A). Cell-cell interactions are often mediated by receptor-ligand pairs, with aneuploid cells acting as either source or receptor cells. Through Cellchat enrichment analysis, we identified 96 pairs of important receptor-ligand pairs, with a significant enrichment of aneuploid cell-derived interactions with other cell types compared to their interactions as receptor cells (Supplementary File 1). Particularly, the potential receptor-ligand pairs involving aneuploid cells acting on macrophages were found to be at a higher level compared to other cell types, including the CCL16-CCR1 axis (Figure 2B). These findings suggest that cell-cell interactions initiated by aneuploid cells as source cells are crucial role in the immune microenvironment of liver cancer.
Figure 1 Single-cell analysis identified cellular subpopulation categories and compositional features of liver cancer tissue samples. (A) The t-SNE clustering plot was generated after quality control and integration using the canonical correlation analysis (CCA) method. The left panel shows the t-SNE clustering of 10 HCC patients. The middle panel represents cell subgroups identified using the Seurat R package’s FindClusters function with a resolution of 0.5. The right panel displays the t-SNE clustering of cell subpopulations after manual annotation. (B) Based on the ascending order of Wilcoxon test values, the top 5 marker genes for each cell type were identified using the FindMarkers function. (C) The number of cells for each cell type and the cellular composition of each sample. (D) The t-SNE distribution of malignant aneuploid cells identified by CopyKAT.
Figure 2 Cellular communication features of malignant liver cancer cells in the microenvironment. (A) The number (top) and weight (bottom) of interactions between different cell types as identified by CellChat analysis. (B) Enrichment results of receptor-ligand pairs, with aneuploid cells as the source (left) and receptor cells (right), as determined by CellChat.
Molecular subtyping and immune cell infiltration characteristics based on prognostic receptor-ligand genes
As is well-established, tumor cells exhibit pronounced heterogeneity, and similarly, the intercellular interactions based on aneuploid cells are also expected to display such complexity (29). Therefore, leveraging the previously identified core set of receptor-ligand genes (comprising a total of 76 genes) associated with aneuploid cells, we conducted an extensive investigation of prognosis and molecular subtyping in two large-scale cohorts, namely TCGA-LIHC (training cohort) and LIRI-JP (testing cohort). Using a strict univariate Cox regression analysis with a p-value threshold of less than 0.1, we were able to identify 8 receptor-ligand genes that have significant prognostic value in the TCGA-LIHC cohort. Significantly, CCL16 stood out as the only factor posing a risk among these candidates (Figure 3A). Using the expression patterns of these eight genes, patients from both cohorts were successfully classified into three distinct subtypes through unsupervised clustering analysis (Figure 3B). Significantly, all three subcategories demonstrated comparable prognostic disparities in both datasets, underscoring the strength and uniformity of the categorization (Figure 3C). The analysis of immune cell infiltration indicated that the C1 subtype, which had the worst prognosis, displayed notably reduced scores for the majority of immune cell categories, specifically T cells, including Activated CD4 T, Central memory CD4 T, and Central memory CD8 T, in contrast to the C2 and C3 subtypes (Supplementary Figure 1A). Consistently, the results of the Estimate analysis also indicated a notable decrease in the Stromal score and Immune score within the C1 subtype (Supplementary Figure 1B). Moreover, the levels of gene expression for most immune checkpoint genes were notably decreased in the C1 subtype when compared to the other two subtypes, indicating a possible diminished immune response in the C1 subtype (Supplementary Figure 1C). According to the TIDE analysis, it was found that the C1 subtype exhibited notably reduced scores in relation to the IFNG pathway, which is linked to the innate immune responses (30). Remarkably, in circumstances of generally limited immune cell penetration, the C1 category, associated with the most unfavorable prognosis, exhibited notably increased ratings for M2 macrophages linked to tumors (Supplementary Figure 1D). Therefore, we speculate that the polarization of macrophages towards M2 phenotype might have a crucial impact on the progression of liver cancer.
Figure 3 Molecular subtyping based on prognostic receptor-ligand genes in TCGA-LIHC and LIRI-JP cohorts. (A) Univariate Cox regression analysis results of receptor-ligand genes in the TCGA-LIHC cohort. (B) Unsupervised clustering of the two cohorts based on the expression pattern of prognosis-associated receptor-ligand genes. (C) Kaplan-Meier survival curves of the three subtypes in both cohorts. Log-rank test.
Construction of the four-gene risk model and prognostic nomogram
Afterwards, we conducted differential gene expression analysis on the three subtypes. The genes that showed differential expression were analyzed for functional enrichment (Supplementary Figure 2A, B; Supplementary File 2) and used to construct a prognostic model. Initially, we conducted univariate Cox analysis on these genes that were expressed differentially (DEGs) in the TCGA-LIHC cohort, employing a threshold of P < 0.05. This yielded 274 prognostic-related DEGs. Following that, an analysis using Lasso regression was performed, leading to the identification of 13 genes (Figure 4A). By employing a stepwise regression approach, we eventually pinpointed four genes and developed a risk model to calculate the risk score with the following formula: risk score = -0.349 * CSF2RA + 0.245 * LRRC3 + 0.095 * UGT3A1 + 0.115 * EFHD1. Risk scores were computed for every patient in both the TCGA-LIHC and LIRI-JP cohorts (Supplementary File 3). Patients were categorized into high and low-risk groups based on the threshold set by the median value. The findings indicated that the high-risk category exhibited notably worse prognosis. In addition, the AUC plots demonstrated good predictive precision (Figures 4B, C). In order to improve the accuracy of predictions, we took into account the inclusion of clinical factors in the model. In the multivariable Cox analysis (Figure 5A), it was found that the risk score and stage continued to be significant prognostic factors, demonstrating their independence. Afterwards, we built a comprehensive calibration graph (Figure 5B). According to Figure 5C, the AUC curve indicated that the nomogram outperformed other clinical factors like the risk score and stage in terms of predictive accuracy. Moreover, the calibration plot demonstrated a strong correlation between the estimated and observed occurrence frequencies (Figure 5D), while the decision curve analysis (DCA) graph displayed the greatest overall net advantage for the nomogram (Figure 5E). In addition, a substantial inverse relationship was noted between the risk score and immune score (R = -0.23, P = 6.4e-6) (Supplementary Figure 3A). The group with minimal risk displayed notably elevated levels of Tumor-Reactive T Cell Signature (TRS score), which measures the responsiveness of T cells towards tumors (31), and increased levels of cytolytic activity score (CYT score), which evaluates the effectiveness of T cells in targeting tumor cells within the tumor microenvironment (32). Nevertheless, no notable distinction was detected between the two cohorts regarding the Th1/IFNγ genetic pattern (Supplementary Figure 3B). In contrast to the high-risk category, the low-risk category demonstrated elevated PD-1 and CTLA4 expression, along with a notable rise in the infiltration of immune cells related to T cells (Supplementary Figure 3C, D).
Figure 4 Construction of the 4-gene risk model. (A) Lasso regression analysis was performed to assess model fit. The analysis included two components. Left: The curve of partial-likelihood deviance was plotted against Log(λ), where a smaller value indicated a better fit of the model. Right: The curve of regression coefficients for each variable was plotted against the change of Log(λ). (B) Kaplan-Meier survival curves and 1, 3, 5-year ROC curves were generated for the high-risk and low-risk groups in the TCGA-LIHC cohort. (C) Kaplan-Meier survival curves and 1, 3, 5-year ROC curves were generated for the high-risk and low-risk groups in the LIRI-JP cohort. The Log-Rank test was used for statistical analysis.
Figure 5 Integrated nomogram with improved predictive performance. (A) Results of univariate and multivariate analyses of risk score and clinical variables in the TCGA-LIHC cohort. P < 0.05 indicates statistical significance. (B) Integrated prognostic nomogram combining risk score and stage for predicting 1, 3, and 5-year patient outcomes. ***: P < 0.001 in the multivariate Cox regression analysis. (C) AUC curves of the nomogram and other variables, showing the highest accuracy of the nomogram across different survival time points. (D) Calibration curves of the nomogram, with closer proximity to the diagonal line indicating a closer match between predicted and actual event rates. (E) Decision curve analysis (DCA) curves for different variables, demonstrating the highest standard net benefit of the nomogram.
Tumor cell-derived CCL16 mediates the recruitment and M2 polarization of macrophages in the liver cancer microenvironment
Through single-cell data analysis, we identified specific expression of the chemokine CCL16 in liver cancer cells, while its receptor CCR1 was specifically expressed in macrophages (Figure 6A). Our hypothesis suggests that the heightened occurrence of M2 macrophages in the microenvironment of liver cancer could potentially be attributed to CCL16-CCR1. In order to examine this hypothesis, we initially examined the CCL16 expression in various liver cancer cell lines by utilizing the CCLE database (Figure 6B). Subsequently, we validated the mRNA expressions of CCL16 in several liver cancer cell lines available in our laboratory through qPCR, which were consistent with the database results, with HEPG2 cells showing the highest expression and SNU761 cells showing the lowest expression (Figure 6C). Next, to confirm whether CCL16 expressed by liver cancer cells is associated with macrophage recruitment, we generated HEPG2 cells with CCL16 knockdown and SNU761 cells with CCL16 overexpression (Figures 6D, E). Transwell experiments revealed that knocking down CCL16 expression in HEPG2 cells significantly reduced their recruitment of THP1 cells, while overexpressing CCL16 promoted the recruitment of THP1 cells by SNU761 cells (Figures 6F, G). Given the widespread existence of tumor-associated macrophages (TAMs) in the tumor microenvironment, known for their M2-polarized characteristics linked to tumor advancement, we utilized a flow cytometry technique (33) that was previously explained to assess the polarization condition of macrophages in the co-culture setup. The schematic diagram illustrating the cell co-culture system and the macrophage migration assay is presented (Figure 6H). M1 and M2 macrophages were distinguished using CD80 as a common marker, whereas CD206 was exclusively used as a marker for M2 macrophages (34). The findings indicated that the combination of THP1 cells with HEPG2 cells lacking CCL16 expression notably decreased the percentage of CD80+CD206+ cells in THP1.In contrast, co-culturing THP1 cells with SNU761 cells overexpressing CCL16 significantly increased the proportion of CD80+CD206+ cells in THP1 (Figure 6I). The findings present convincing proof that cancer cells attract macrophages and support their M2 polarization by releasing CCL16.
Figure 6 Tumor cell-derived CCL16 mediates the recruitment and M2 polarization of macrophages in the liver cancer microenvironment. (A) Expression of CCL16 and CCR1 in different cell types at the single-cell level. (B) Expression levels of CCL16 in different liver cancer cell lines from the CCLE database. (C) mRNA expression of CCL16 in different liver cancer cell lines detected by qPCR. (D) Validation of CCL16 knockdown in HEPG2 cell line. (E) Validation of CCL16 overexpression in SNU761 cell line. (F) Transwell assay to evaluate the recruitment of THP1 cells by control and CCL16 knockdown HEPG2 cells, and ELISA assay to measure CCL16 concentration in the culture supernatant of HEPG2 cells. (G) Transwell assay to evaluate the recruitment of THP1 cells by control and CCL16 overexpressing SNU761 cells, and ELISA assay to measure CCL16 concentration in the culture supernatant of SNU761 cells. (H) Schematic diagram illustrating the cell co-culture system and the macrophage migration assay. (I) Flow cytometry analysis to detect the proportion of M2-polarized cells after co-culture of THP1 cells with CCL16 knockdown or overexpressing tumor cells. Transwell Scale Bar = 100μm. Three independent replicates were conducted. Statistical data are presented as mean ± SD, and each data point represents an independent measurement. Unpaired Student’s t-test. ns, not significant; **: P < 0.01; ***: P < 0.001; ****: P < 0.0001.
The recruitment of tumor-associated macrophages mediated by CCL16 depends on the macrophage receptor CCR1
Next, we aimed to identify the receptor through which CCL16 acts on macrophages. There have been studies reporting that CCL16 binds to known receptors including CCR1, CCR2, CCR5, and CCR8 (14–17). We added synthetic Flag-CCL16 protein into THP1 cells cultured in vitro and performed co-immunoprecipitation experiments. The results revealed that CCL16 predominantly interacts with the CCR1 receptor on macrophages, while the binding affinity to other receptors was minimal (Figure 7A). This interaction was further confirmed by immunofluorescence, showing colocalization of CCL16 and CCR1 on the cell membrane of THP1 cells (Figure 7B). To further investigate the role of CCR1, we performed CCR1 knockdown in THP1 cells (Figure 7C) and co-cultured them with HEPG2 cells to assess the migration ability of macrophages. The results demonstrated that CCR1 knockdown in THP1 cells significantly inhibited the recruitment of macrophages by HEPG2 cells (Figure 7D). Interestingly, after treatment with the CCR1 inhibitor BX471, the overexpression of CCL16 no longer promoted THP1 cell recruitment (Figure 7E). Consistently, when CCR1 was knocked down in THP1 cells, co-culture with tumor cells overexpressing CCL16 no longer facilitated macrophage recruitment (Figure 7F). These results demonstrate that the secretion of CCL16 by liver cancer cells promotes macrophage recruitment through binding to the CCR1 receptor on macrophages.
Figure 7 The recruitment of tumor-associated macrophages mediated by CCL16 depends on the macrophage receptor CCR1. (A) Immunoprecipitation assay to detect the interaction between Flag-CCL16 and CCR1, CCR2, CCR5, CCR8 in THP1 cell culture medium. (B) Immunofluorescence assay to detect the co-localization of Flag-CCL16 and CCR1 in THP1 cells after the addition of Flag-CCL16. Scale bar = 20 μm. (C) qPCR validation of CCR1 knockdown in THP1 cells. (D) Transwell assay to evaluate the cell migration ability of CCR1 knockdown THP1 cells co-cultured with HEPG2 cells. (E) Recruitment of THP1 cells by control or 5μM BX471-treated CCL16 overexpressing SNU761 cells after 24 hours. (F) Recruitment of THP1 cells by CCL16 overexpressing SNU761 cells after CCR1 knockdown. Transwell Scale Bar = 100μm. Three independent replicates were conducted. Statistical data are presented as mean ± SD, and each data point represents an independent measurement. Unpaired Student’s t-test or two-way ANOVA was used for statistical analysis. ns, not significant; ***: P < 0.001; ****: P < 0.0001.
M2 macrophage infiltration and CCR1 expression in clinical tissues are positively correlated with CCL16
In order to examine the involvement of the CCL16-CCR1 axis in clinical samples, a total of 42 liver cancer specimens were gathered. The specimens were subjected to immunohistochemical analysis to assess the expression of CCL16. The specimens were categorized into groups of high and low expression based on the median MOD value. Table 3 provides a summary of the variations in clinical characteristics observed in both groups. The results of the chi-squared test indicated a notable correlation between the levels of CCL16 protein expression and advanced stage, as well as the presence of lymph node and distant metastasis. Additionally, immunohistochemical staining was performed to evaluate the presence of CCR1, CD68, and CD206 expression. Significant positive associations between CCL16 and the other three proteins were validated (Figure 8A). In addition, the proportions of CD68+CCR1+ macrophages were determined using immunofluorescence staining, revealing a robust positive association between CCL16 expression and infiltration of CD68+CCR1+ macrophages (Figure 8B). The clinical pathological analyses confirm that the CCL16-CCR1 axis has the capability to enhance the infiltration of M2 macrophages in the microenvironment of liver cancer, indicating that targeting the CCL16-CCR1 axis could be a promising therapeutic approach for liver cancer.
Table 3 Clinical characteristics of patients with high and low expression of CCL16.
Figure 8 High expression of CCL16 is positively correlated with the infiltration of M2 macrophages and the expression of CCR1 in clinical tissues. (A) Immunohistochemistry representative images and Pearson correlation analysis of CCL16 with CCR1, CD68, and CD206 in liver cancer patient samples, as well as the Mean Optical Density (MOD) values. Scale Bar = 100μm. (B) Immunofluorescence detection of CCR1+ macrophage infiltration. Scale Bar = 20μm. Statistical analysis of the difference in CD68+CCR1+ cell numbers between high and low CCL16 expression groups using unpaired Student’s t-test. **: P < 0.01. Pearson correlation analysis of CCL16 expression and CCR1+ macrophage infiltration, r = 0.5743, P < 0.001.
Discussion
Hepatocellular carcinoma (HCC), a prevalent and aggressive cancer, has a significant global impact and a high fatality rate (1). The interaction among cells facilitated by ligand-receptor complexes, which contributes to the diversity of the tumor microenvironment, has been acknowledged as a vital element impacting the prognosis of HCC patients and propelling tumor advancement (35). Song et al. provide a comprehensive perspective on analyzing the immune microenvironment in HCC and emphasize the presence of a unique subset of macrophages characterized by CCL18 and CREM expression, which is notably enriched in advanced HCC patients. This specific macrophage subset plays a significant role in driving tumor progression and holds promising potential for future immunotherapeutic strategies (36). Chen et al. examined the role of cancer-associated fibroblast (CAF)-induced M2-polarized macrophages in promoting the progression of HCC through the plasminogen activator inhibitor-1 pathway. Additionally, they assessed macrophage polarization and identified key paracrine factors involved in their interactions with both CAFs and cancer cells in the context of HCC (37). Li et al. made a significant finding demonstrating that targeting the CCL2/CCR2 axis therapeutically could effectively impede the recruitment of inflammatory monocytes, inhibit tumor-associated macrophage (TAM) infiltration, and reverse M2 polarization. Consequently, this intervention effectively counteracted the immune-suppressive conditions within the tumor microenvironment, subsequently activating anti-tumor CD8 T cell responses (38). Yang et al. found that upregulated CD36 in metastasis-associated macrophages (MAMs) promoted M2 polarization, facilitates liver cancer metastasis through interactions with tumor cells, and its loss in MAMs can attenuate liver metastasis in mice (39).
Chemokines derived from tumor cells play a crucial role in the complex tumor microenvironment, with dysregulated cytokine production in the tumor microenvironment influencing all stages of carcinogenesis and affecting cancer initiation, progression, and responses to therapy (40). Zha et al. discovered tumor cells utilized complement-derived C3 to inhibit antitumor immunity by regulating tumor-associated macrophages through the C3a-C3aR-PI3Kγ pathway (41). Park et al. employed quantitative proteomics to unveil that exosomes originating from tumor cells under hypoxic conditions exhibit a significant enrichment of immunomodulatory proteins and chemokines. This enrichment includes notable factors such as CSF-1, CCL2, FTH, FTL, and TGFβ, all of which contribute to the promotion of macrophage M2 polarization (42). In both mouse models and immunotherapy-treated patients, House et al. found that macrophages, as the main source of CXCL9, significantly upregulated the ligand of CXCR3, which, upon CXCL9 knockout, led to reduced CD8 T cell infiltration and compromised therapeutic efficacy of PD-1/CTLA-4 immune checkpoint blockade therapy (43). Korbecki et al. emphasized the significant role of human CC motif chemokine ligands and their corresponding receptors in mediating the chemotaxis of immune cells to the tumor microenvironment (44). However, the role of CCL16 in the hepatocellular carcinoma microenvironment remains to be further investigated. To uncover the underlying mechanisms, we conducted a series of experiments, including Transwell co-culture and migration assays, immunofluorescence, immunoprecipitation, flow cytometry, and immunohistochemistry. Specifically, in co-culture experiments, CCL16-overexpressing hepatocellular carcinoma cells promoted migration and M2 polarization of macrophages towards tumor cells. This effect was mediated by the interaction between CCL16 secreted by tumor cells and the CCR1 receptor on macrophages. In clinical tissue samples, we observed a significant positive correlation between CCL16 and CCR1, as well as CD206 and CD68 at the protein level. Furthermore, CCL16 high-expressing patient tissues showed a significant increase in CD68+CCR1+ cells, further validating the recruitment role of CCL16 in CCR1-positive macrophages. These experiments allowed us to elucidate the role of the CCL16-CCR1 axis in recruiting tumor-associated macrophages and promoting M2 polarization in HCC.
Using the CellChat algorithm, we conducted an analysis of single-cell transcriptomic data obtained from HCC tumor samples. Our primary objective was to investigate the interactions between different cell populations within the HCC microenvironment, with a specific focus on ligand-receptor interactions involving malignant HCC cells. In this context, we identified key genes, such as CCL16 and CCR1, that are instrumental in these interactions. Furthermore, we leveraged RNA-Seq data from comprehensive datasets, including TCGA-LIHC and LIRI-JP, and employed Cox regression analysis to uncover predictive genes. We also employed unsupervised clustering to categorize individuals into molecular subgroups. Additionally, we developed a risk model comprising four genes by analyzing differential gene expression among subtypes. This model can be utilized for predicting the prognosis of HCC patients. Moreover, through bioinformatics analysis of the immune microenvironment, we made the noteworthy observation that the subtype associated with the poorest prognosis exhibited a significant increase in the infiltration score of M2 tumor-associated macrophages.
Our study provides valuable insights into the molecular subtyping, prognostic modeling, and the immune microenvironment of HCC. We achieved this by characterizing interactions between ligands and receptors and uncovering the pro-tumorigenic role of the CCL16-CCR1 axis within the HCC microenvironment. The findings from this research could have implications for the development of novel treatment strategies targeting specific ligand-receptor interactions or the modulation of the immune microenvironment in HCC. It’s worth noting that our study did not involve animal trials due to the absence of a homologous gene in mice, which presents a significant limitation when investigating the role of CCL16 in the immune microenvironment. Consistent with this, a study published in the journal Cell demonstrated that CCL16, produced by hepatocytes, binds to CCR1 expressed by human Kupffer cells (KCs) but not murine KCs. In this context, KCs refer to hepatic macrophages. Therefore, human CCL16 cannot interact with murine macrophages through CCR1 (45). On the other hand, as the first-discovered C-C chemokine receptor, CCR1 is overexpressed in several types of cancers and is associated with increased immune-suppressive cell infiltration and tumor metastasis (46, 47). Through literature and patent searches, several CCR1 antagonists are currently in development, including AZD-4818, BI-638683, BL-5923, BX-471, C-6448, C-4462, CCX9588, CCX354, CCX721, CP-481715, MLN-3701, MLN-3897, PS-031291/PS-375179, and UCB-35625. However, as of now, none of them have entered clinical trials in the field of oncology (48). The selective CCR1 antagonist CCX721 has demonstrated efficacy in reducing tumor burden and osteolytic lesions in a murine model of multiple myeloma (MM) by blocking osteoclasts (49). Additionally, studies have reported that inhibiting CCR1 with the receptor antagonist BL5923 suppresses the recruitment of immature myeloid cells, leading to a reduction in metastatic colon cancer and a significant prolongation of survival in mice with colon cancer liver metastasis (50). The combination of the CCR1 antagonist CCX9588 with an anti-PD-L1 antibody has shown promise as a therapeutic approach, synergistically inhibiting primary tumor growth and lung metastasis in an in situ breast cancer mouse model (51). Recently, in a mouse model of ovarian cancer, the small-molecule CCR1 inhibitor UCB35625 demonstrated the ability to reduce cell migration to the greater omentum, a preferential metastatic site for such cancers (52). Overall, these findings suggest that targeting CCR1 is a viable therapeutic strategy capable of limiting metastasis and delaying disease progression.
Data availability statement
The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding author.
Ethics statement
The studies involving humans were approved by the Ethics Committee at the First Hospital of China Medical University. The studies were conducted in accordance with the local legislation and institutional requirements. The participants provided their written informed consent to participate in this study.
Author contributions
ZD: Conceptualization, Formal Analysis, Investigation, Writing – original draft. YW: Conceptualization, Formal Analysis, Resources, Writing – review & editing. NS: Conceptualization, Investigation, Validation, Writing – review & editing. CZ: Conceptualization, Resources, Supervision, Validation, Writing – review & editing.
Funding
The author(s) declare financial support was received for the research, authorship, and/or publication of this article. This work is supported by grants from the Natural Science Foundation of Liaoning Province [No. 2021-MS-187].
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.2023.1299953/full#supplementary-material
Supplementary Figure 1 | Analysis of immune microenvironment differences among the three molecular subtypes. (A) Comparison of differences in immune cell infiltration gene set scores obtained through ssGSEA method. (B) Comparison of differences in Stromal score and Immune Score obtained through Estimate algorithm. (C) Differential expression of immune checkpoint genes among the three subtypes. (D) Comparison of differences in immune feature scores calculated by TIDE tool. Wilcoxon test. ns, not significant; *: P<0.05; **: P<0.01; ***: P<0.001; ****: P<0.0001.
Supplementary Figure 2 | Identification of differentially expressed genes among three subtypes. (A) Volcano plots showing differential expression analysis using the Limma R package to identify genes that are differentially expressed between each group and the other two groups. (B) Functional enrichment analysis results of the commonly differentially expressed genes.
Supplementary Figure 3 | Differential analysis of the immune microenvironment between high and low-risk groups. (A) There is a significant negative correlation between the risk score and immune score (R = -0.23, P = 6.4e-6). (B) Differences in Tumor-Reactive T Cell Signature (TRS score), cytolytic activity score (CYT score), and Th1/IFNγ gene signature score between high and low-risk groups. (C) Differential analysis of the scores of 28 immune cell gene sets between the two groups. Statistical analysis was performed using Wilcoxon test. ns, not significant; *: P < 0.05; **: P < 0.01; ***: P < 0.001; ****: P < 0.0001.
Supplementary file 1 | The receptor-ligand pairs and pathway enrichment results for the interactions between different cell types, as identified by CellChat.
Supplementary file 2 | The differential expression genes and functional enrichment results among the three molecular subtypes.
Supplementary file 3 | The results of univariate Cox analysis in the TCGA-LIHC cohort and risk scores of patients in both cohorts.
References
1. Vogel A, Meyer T, Sapisochin G, Salem R, Saborowski A. Hepatocellular carcinoma. Lancet (2022) 400(10360):1345–62. doi: 10.1016/s0140-6736(22)01200-4
2. Yang JD, Hainaut P, Gores GJ, Amadou A, Plymoth A, Roberts LR. A global view of hepatocellular carcinoma: trends, risk, prevention and management. Nat Rev Gastroenterol Hepatol (2019) 16(10):589–604. doi: 10.1038/s41575-019-0186-y
3. Chen Y, Zhou Y, Yan Z, Tong P, Xia Q, He K. Effect of infiltrating immune cells in tumor microenvironment on metastasis of hepatocellular carcinoma. Cell Oncol (Dordr) (2023) 46(6):1595–604. doi: 10.1007/s13402-023-00841-6
4. Cariani E, Missale G. Immune landscape of hepatocellular carcinoma microenvironment: Implications for prognosis and therapeutic applications. Liver Int (2019) 39(9):1608–21. doi: 10.1111/liv.14192
5. Filiberti S, Russo M, Lonardi S, Bugatti M, Vermi W, Tournier C, et al. Self-renewal of macrophages: tumor-released factors and signaling pathways. Biomedicines (2022) 10(11):2709. doi: 10.3390/biomedicines10112709
6. Zheng H, Peng X, Yang S, Li X, Huang M, Wei S, et al. Targeting tumor-associated macrophages in hepatocellular carcinoma: biology, strategy, and immunotherapy. Cell Death Discovery (2023) 9(1):65. doi: 10.1038/s41420-023-01356-7
7. Xu W, Cheng Y, Guo Y, Yao W, Qian H. Targeting tumor associated macrophages in hepatocellular carcinoma. Biochem Pharmacol (2022) 199:114990. doi: 10.1016/j.bcp.2022.114990
8. Weng J, Liu S, Zhou Q, Xu W, Xu M, Gao D, et al. Intratumoral PPT1-positive macrophages determine immunosuppressive contexture and immunotherapy response in hepatocellular carcinoma. J Immunother Cancer (2023) 11(6):e006655. doi: 10.1136/jitc-2022-006655
9. Zhao HC, Chen CZ, Tian YZ, Song HQ, Wang XX, Li YJ, et al. CD168(+) macrophages promote hepatocellular carcinoma tumor stemness and progression through TOP2A/β-catenin/YAP1 axis. iScience (2023) 26(6):106862. doi: 10.1016/j.isci.2023.106862
10. Wang HC, Haung LY, Wang CJ, Chao YJ, Hou YC, Yen CJ, et al. Tumor-associated macrophages promote resistance of hepatocellular carcinoma cells against sorafenib by activating CXCR2 signaling. J BioMed Sci (2022) 29(1):99. doi: 10.1186/s12929-022-00881-4
11. Arvanitakis K, Koletsa T, Mitroulis I, Germanidis G. Tumor-associated macrophages in hepatocellular carcinoma pathogenesis, prognosis and therapy. Cancers (Basel) (2022) 14(1):226. doi: 10.3390/cancers14010226
12. Sung PS. Crosstalk between tumor-associated macrophages and neighboring cells in hepatocellular carcinoma. Clin Mol Hepatol (2022) 28(3):333–50. doi: 10.3350/cmh.2021.0308
13. Shen W, Zhang X, Tang J, Zhang Z, Du R, Luo D, et al. CCL16 maintains stem cell-like properties in breast cancer by activating CCR2/GSK3β/β-catenin/OCT4 axis. Theranostics (2021) 11(5):2297–317. doi: 10.7150/thno.51000
14. Nomiyama H, Hieshima K, Nakayama T, Sakaguchi T, Fujisawa R, Tanase S, et al. Human CC chemokine liver-expressed chemokine/CCL16 is a functional ligand for CCR1, CCR2 and CCR5, and constitutively expressed by hepatocytes. Int Immunol (2001) 13(8):1021–9. doi: 10.1093/intimm/13.8.1021
15. Howard OM, Dong HF, Shirakawa AK, Oppenheim JJ. LEC induces chemotaxis and adhesion by interacting with CCR1 and CCR8. Blood (2000) 96(3):840–5. doi: 10.1182/blood.V96.3.840
16. Strasly M, Doronzo G, Cappello P, Valdembri D, Arese M, Mitola S, et al. CCL16 activates an angiogenic program in vascular endothelial cells. Blood (2004) 103(1):40–9. doi: 10.1182/blood-2003-05-1387
17. Kim IS, Jang SW, Sung HJ, Lee JS, Ko J. Differential CCR1-mediated chemotaxis signaling induced by human CC chemokine HCC-4/CCL16 in HOS cells. FEBS Lett (2005) 579(27):6044–8. doi: 10.1016/j.febslet.2005.09.064
18. Goldman MJ, Craft B, Hastie M, Repečka K, McDade F, Kamath A, et al. Visualizing and interpreting cancer genomics data via the Xena platform. Nat Biotechnol (2020) 38(6):675–8. doi: 10.1038/s41587-020-0546-8
19. Lu Y, Yang A, Quan C, Pan Y, Zhang H, Li Y, et al. A single-cell atlas of the multicellular ecosystem of primary and metastatic hepatocellular carcinoma. Nat Commun (2022) 13(1):4594. doi: 10.1038/s41467-022-32283-3
20. Hao Y, Hao S, Andersen-Nissen E, Mauck WM 3rd, Zheng S, Butler A, et al. Integrated analysis of multimodal single-cell data. Cell (2021) 184(13):3573–87.e29. doi: 10.1016/j.cell.2021.04.048
21. Gao R, Bai S, Henderson YC, Lin Y, Schalck A, Yan Y, et al. Delineating copy number and clonal substructure in human tumors from single-cell transcriptomes. Nat Biotechnol (2021) 39(5):599–608. doi: 10.1038/s41587-020-00795-2
22. Jin S, Guerrero-Juarez CF, Zhang L, Chang I, Ramos R, Kuan CH, et al. Inference and analysis of cell-cell communication using CellChat. Nat Commun (2021) 12(1):1088. doi: 10.1038/s41467-021-21246-9
23. Wilkerson MD, Hayes DN. ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking. Bioinformatics (2010) 26(12):1572–3. doi: 10.1093/bioinformatics/btq170
24. Charoentong P, Finotello F, Angelova M, Mayer C, Efremova M, Rieder D, et al. Pan-cancer immunogenomic analyses reveal genotype-immunophenotype relationships and predictors of response to checkpoint blockade. Cell Rep (2017) 18(1):248–62. doi: 10.1016/j.celrep.2016.12.019
25. Danilova L, Ho WJ, Zhu Q, Vithayathil T, De Jesus-Acosta A, Azad NS, et al. Programmed cell death ligand-1 (PD-L1) and CD8 expression profiling identify an immunologic subtype of pancreatic ductal adenocarcinomas with favorable survival. Cancer Immunol Res (2019) 7(6):886–95. doi: 10.1158/2326-6066.Cir-18-0822
26. Jiang P, Gu S, Pan D, Fu J, Sahu A, Hu X, et al. Signatures of T cell dysfunction and exclusion predict cancer immunotherapy response. Nat Med (2018) 24(10):1550–8. doi: 10.1038/s41591-018-0136-1
27. Fu J, Li K, Zhang W, Wan C, Zhang J, Jiang P, et al. Large-scale public data reuse to model immunotherapy response and resistance. Genome Med (2020) 12(1):21. doi: 10.1186/s13073-020-0721-z
28. Yu G, Wang LG, Han Y, He QY. clusterProfiler: an R package for comparing biological themes among gene clusters. Omics (2012) 16(5):284–7. doi: 10.1089/omi.2011.0118
29. Lu LC, Hsu CH, Hsu C, Cheng AL. Tumor heterogeneity in hepatocellular carcinoma: facing the challenges. Liver Cancer (2016) 5(2):128–38. doi: 10.1159/000367754
30. Qin Y, Zhang C. The regulatory role of IFN-γ on the proliferation and differentiation of hematopoietic stem and progenitor cells. Stem Cell Rev Rep (2017) 13(6):705–12. doi: 10.1007/s12015-017-9761-1
31. Yan M, Hu J, Ping Y, Xu L, Liao G, Jiang Z, et al. Single-cell transcriptomic analysis reveals a tumor-reactive T cell signature associated with clinical outcome and immunotherapy response in melanoma. Front Immunol (2021) 12:758288. doi: 10.3389/fimmu.2021.758288
32. Rooney MS, Shukla SA, Wu CJ, Getz G, Hacohen N. Molecular and genetic properties of tumors associated with local immune cytolytic activity. Cell (2015) 160(1-2):48–61. doi: 10.1016/j.cell.2014.12.033
33. Raggi F, Pelassa S, Pierobon D, Penco F, Gattorno M, Novelli F, et al. Regulation of human macrophage M1-M2 polarization balance by hypoxia and the triggering receptor expressed on myeloid cells-1. Front Immunol (2017) 8:1097. doi: 10.3389/fimmu.2017.01097
34. Cossarizza A, Chang HD, Radbruch A, Abrignani S, Addo R, Akdis M, et al. Guidelines for the use of flow cytometry and cell sorting in immunological studies (third edition). Eur J Immunol (2021) 51(12):2708–3145. doi: 10.1002/eji.202170126
35. Barcena-Varela M, Lujambio A. The endless sources of hepatocellular carcinoma heterogeneity. Cancers (Basel) (2021) 13(11):2621. doi: 10.3390/cancers13112621
36. Song G, Shi Y, Zhang M, Goswami S, Afridi S, Meng L, et al. Global immune characterization of HBV/HCV-related hepatocellular carcinoma identifies macrophage and T-cell subsets associated with disease progression. Cell Discovery (2020) 6(1):90. doi: 10.1038/s41421-020-00214-5
37. Chen S, Morine Y, Tokuda K, Yamada S, Saito Y, Nishi M, et al. Cancer−associated fibroblast−induced M2−polarized macrophages promote hepatocellular carcinoma progression via the plasminogen activator inhibitor−1 pathway. Int J Oncol (2021) 59(2):59. doi: 10.3892/ijo.2021.5239
38. Li X, Yao W, Yuan Y, Chen P, Li B, Li J, et al. Targeting of tumour-infiltrating macrophages via CCL2/CCR2 signalling as a therapeutic strategy against hepatocellular carcinoma. Gut (2017) 66(1):157–67. doi: 10.1136/gutjnl-2015-310514
39. Yang P, Qin H, Li Y, Xiao A, Zheng E, Zeng H, et al. CD36-mediated metabolic crosstalk between tumor cells and macrophages affects liver metastasis. Nat Commun (2022) 13(1):5782. doi: 10.1038/s41467-022-33349-y
40. Saxena S, Singh RK. Chemokines orchestrate tumor cells and the microenvironment to achieve metastatic heterogeneity. Cancer Metastasis Rev (2021) 40(2):447–76. doi: 10.1007/s10555-021-09970-6
41. Zha H, Wang X, Zhu Y, Chen D, Han X, Yang F, et al. Intracellular activation of complement C3 leads to PD-L1 antibody treatment resistance by modulating tumor-associated macrophages. Cancer Immunol Res (2019) 7(2):193–207. doi: 10.1158/2326-6066.Cir-18-0272
42. Park JE, Dutta B, Tse SW, Gupta N, Tan CF, Low JK, et al. Hypoxia-induced tumor exosomes promote M2-like macrophage polarization of infiltrating myeloid cells and microRNA-mediated metabolic shift. Oncogene (2019) 38(26):5158–73. doi: 10.1038/s41388-019-0782-x
43. House IG, Savas P, Lai J, Chen AXY, Oliver AJ, Teo ZL, et al. Macrophage-derived CXCL9 and CXCL10 are required for antitumor immune responses following immune checkpoint blockade. Clin Cancer Res (2020) 26(2):487–504. doi: 10.1158/1078-0432.Ccr-19-1868
44. Korbecki J, Grochans S, Gutowska I, Barczak K, Baranowska-Bosiacka I. CC chemokines in a tumor: A review of pro-cancer and anti-cancer properties of receptors CCR5, CCR6, CCR7, CCR8, CCR9, and CCR10 ligands. Int J Mol Sci (2020) 21(20):7619. doi: 10.3390/ijms21207619
45. Guilliams M, Bonnardel J, Haest B, Vanderborght B, Wagner C, Remmerie A, et al. Spatial proteogenomics reveals distinct and evolutionarily conserved hepatic macrophage niches. Cell (2022) 185(2):379–96.e38. doi: 10.1016/j.cell.2021.12.018
46. Korbecki J, Kojder K, Simińska D, Bohatyrewicz R, Gutowska I, Chlubek D, et al. CC chemokines in a tumor: A review of pro-cancer and anti-cancer properties of the ligands of receptors CCR1, CCR2, CCR3, and CCR4. Int J Mol Sci (2020) 21(21):8412. doi: 10.3390/ijms21218412
47. Bule P, Aguiar SI, Aires-Da-Silva F, Dias JNR. Chemokine-directed tumor microenvironment modulation in cancer immunotherapy. Int J Mol Sci (2021) 22(18):9804. doi: 10.3390/ijms22189804
48. Gilchrist A, Echeverria SL. Targeting chemokine receptor CCR1 as a potential therapeutic approach for multiple myeloma. Front Endocrinol (Lausanne) (2022) 13:846310. doi: 10.3389/fendo.2022.846310
49. Dairaghi DJ, Oyajobi BO, Gupta A, McCluskey B, Miao S, Powers JP, et al. CCR1 blockade reduces tumor burden and osteolysis in vivo in a mouse model of myeloma bone disease. Blood (2012) 120(7):1449–57. doi: 10.1182/blood-2011-10-384784
50. Kitamura T, Fujishita T, Loetscher P, Revesz L, Hashida H, Kizaka-Kondoh S, et al. Inactivation of chemokine (C-C motif) receptor 1 (CCR1) suppresses colon cancer liver metastasis by blocking accumulation of immature myeloid cells in a mouse model. Proc Natl Acad Sci U.S.A. (2010) 107(29):13063–8. doi: 10.1073/pnas.1002372107
51. Jung H, Bischof A, Ebsworth K, Ertl L, Schall T, Charo I. Combination therapy of chemokine receptor inhibition plus PDL-1 blockade potentiates anti-tumor effects in a murine model of breast cancer. J ImmunoTher Cancer (2015) 3(2):P227. doi: 10.1186/2051-1426-3-S2-P227
Keywords: ligand-receptor interaction, macrophages, CCL16, CCR1, hepatocellular carcinoma
Citation: Dai Z, Wang Y, Sun N and Zhang C (2024) Characterizing ligand-receptor interactions and unveiling the pro-tumorigenic role of CCL16-CCR1 axis in the microenvironment of hepatocellular carcinoma. Front. Immunol. 14:1299953. doi: 10.3389/fimmu.2023.1299953
Received: 23 September 2023; Accepted: 26 December 2023;
Published: 11 January 2024.
Edited by:
Dongxiao Li, Henan Provincial People’s Hospital, ChinaCopyright © 2024 Dai, Wang, Sun and Zhang. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Chengshuo Zhang, cszhang@cmu.edu.cn
†These authors have contributed equally to this work