In this study, we acquired gene expression matrix (RNA-seq) and clinical information of HCC patients from The Cancer Genome Atlas (TCGA) (https://portal.gdc.cancer.gov/; TCGA_LIHC; 362 patient samples), International Cancer Genome Consortium (ICGC) (https://dcc.icgc.org/; ICGC_JP; 230 patient samples), and Gene Expression Omnibus (GEO) (GSE14520; 225 patient samples) ( Table 1 ). The tumor samples for RNA-seq included 374 cases for TCGA_LIHC, 243 cases for ICGC_JP, and 225 cases for GSE14520. Accordingly, the adjacent normal samples of RNA-seq included 50 cases for TCGA_LIHC, 202 cases for ICGC_JP, and 220 cases for GSE14520. TCGA_LIHC dataset was used as the internal training cohort while ICGC_JP and GSE14520 datasets were used as the external testing cohorts. GSE14520 dataset was downloaded through “GEOquery” R package, and clinical data was acquired from the website. All TCGA_LIHC and ICGC_JP data, gene-expressed profile and clinical details were manually downloaded from the website.

Meanwhile, 99 genes that are implicated in the CCL18 signaling pathway were obtained from WikiPathways (https://www.wikipathways.org/) ( Table S1 ). Immune-related genes were obtained from the ImmPort database (https://www.immport.org/shared/home).

We utilized the “survival” and “survminer” R packages to perform univariate Cox regression analysis in the TCGA_LIHC training cohort with the aim of identifying genes that have a substantial impact on survival. Consequently, we observed significant variations in the expression of 39 genes. Subsequently, the “glmnet” R package was employed to conduct LASSO regression, enabling the selection of the most crucial variables from this gene set. As a result, a subset of 6 genes, specifically BMI1, CCR3, CDC25C, CFL1, LDHA, and RAC1, were identified as prominent hub genes.

The multivariate Cox regression analysis was employed using the 6 hub genes. The risk score for individual patients in TCGA_LIHC cohort was calculated using the following formula.

In this formula, coefficients were acquired from multivariate Cox regression, where gene(i) represents mRNA expression. The patients in both the training and test cohorts were classified into two groups, namely the high- and low-risk groups, based on the median value of the risk score. Additionally, we examined whether the risk score independently served as a prognostic factor. The clinical characteristics within TCGA_LIHC cohort including age, gender, M, N, T stage, stage and risk score were analyzed through univariate and multivariate Cox regression analysis. Analogous analyses were performed using the ICGC_JP and GSE14520 cohorts.

A predictive nomogram was developed using the 6 hub genes to estimate the survival probability for individual HCC patients. The predictive capability of the nomogram was evaluated in both the training and test cohorts through the utilization of receiver operating characteristic (ROC) curves and calibration curves. The “rms”, “timeROC”, and “pROC” R packages were utilized for these assessments.

To evaluate the influence of the six hub genes on the overall survival rate of TCGA_LIHC cohort, survival curves were generated using the Kaplan-Meier method from the “survival” R package, and the median value was taken as the best cut-off. The Kaplan-Meier Plotter (http://kmplot.com/analysis/) platform was employed to examine the association between the expression levels of prognostic genes and various clinical endpoints, including overall survival (OS), recurrence-free survival (RFS), progression-free survival (PFS), and disease-specific survival (DSS). Additionally, the genes encompassed within pathways were collected and subjected to analysis using the “GSVA” R package, with the parameter method = ‘ssgsea’ being specifically chosen. Subsequently, the correlation between the prognostic genes and the pathways was assessed using Spearman correlation analysis (16).

The “limma” R package was utilized to identify differentially expressed genes (DEGs) in the TCGA_LIHC cohort. This included comparisons between HCC tumor and normal samples, as well as between high- and low-risk groups (Log2 fold change > 1, p value< 0.05). A total of 308 genes were identified as the intersection between the genes exhibiting high expression in tumors and the genes within the high-risk score group. The findings were visually represented through the utilization of volcano and Venn diagrams using “ggvenn”, “tidyverse”, and “ggrepel” R packages, effectively illustrating the intersection of 308 genes exhibiting high expression in tumors and the genes within the high-risk score group.

The functional enrichment analysis of DEGs was systematically performed using the “clusterProfiler” and “org.Hs.eg.db” R packages, which facilitated Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses for a comprehensive understanding of the biological functions involved. The protein-protein interaction (PPI) network was meticulously examined utilizing the STRING database (https://cn.string-db.org/), and a significant molecular cluster was identified through the application of MODE, a Cytoscape plugin (Cytoscape software version 3.9.1).

The estimation of stromal, immune, and ESTIMATE scores was performed using the “ESTIMATE” R package, which provides a computational approach to calculate and quantify the stromal and immune components within the tumor microenvironment. Meanwhile, the “CIBERSORT” R package was utilized to estimate the impact of the risk score on the proportions of 22 immune cell subtypes in the TCGA_LIHC training cohort.

The “OncoPredict” R package was utilized to predict drug sensitivity based on gene expression profiles. This approach enabled the calculation of drug sensitivity values for each sample, with lower values indicating higher efficacy of the drug. Comprehensive drug sensitivity analysis was conducted on all samples to determine the drug sensitivity values for commonly used tumor drugs.

CAMP (https://clue.io/) is an extensive database and analysis platform that offers valuable resources and tools to delve into and comprehend gene expression profiles and drug perturbations. In order to identify potential small-molecule drugs for the treatment of high-risk patients, we utilized the CAMP online database. Specifically, we selected downregulated genes from the pool of significant genes, as well as upregulated genes from the top 150 genes. By inputting this gene set into the CAMP database, we conducted a meticulous screening to identify promising small-molecule drugs capable of modulating the dysregulated gene expression patterns associated with high-risk patients.

L02, Huh7, HepG2, LM3 and MHCC-97H cells were purchased from the Cell Bank of the Chinese Academy of Sciences (Shanghai, China). Cells were incubated in DMEM medium with 10% fetal bovine serum (FBS) and maintained in penicillin (100 IU/mL) and streptomycin (100 mg/mL) in 5% CO2 at 37°C. The plasmids encoding the BMI1, CCR3, CDC25C, CFL1, LDHA, and RAC1 genes were constructed by cloning the sequence of the coding region using the appropriate primers ( Table S2 ) and inserting the fragment into the pcDNA3.1 (+) plasmid. The cells were transfected with the plasmids using Lipofectamine 2000 Transfection Reagent (Invitrogen) and then the medium was changed 6h after transfection.

Total RNA was isolated using the TRIzol (Invitrogen) RNA extraction method, following the manufacturers’ instructions. qRT-PCR measurements were performed as described previously (17) with the appropriate primers listed in Table S3 . Glyceraldehyde-3-phosphate dehydrogenase (GAPDH) was regarded as the internal reference, and the 2^−ΔΔCt method was applied to express the ratio of the target gene expression in the experimental group compared to the control group.

Protein extraction and western blotting analysis were performed using previously standard procedures (17). The following antibodies were used for western blotting: anti-E-cadherin antibody(#20874-1-AP Proteintech, China), anti-N-cadherin antibody(#22018-1-AP Proteintech, China), anti-ZEB1 antibody(#66279-1-Ig Proteintech, China), anti-Vimentin antibody(10366-1-AP Proteintech, China), anti-SOX2 antibody(#11064-1-AP Proteintech, China), anti-GAPDH antibody(#10494-1-AP Proteintech, China), and anti-PD-L1 antibody(#66248-1-Ig Proteintech, China).

Cell proliferation was detected using the Cell Counting Kit-8 (CCK-8, #K1018, APExBIO, USA) and the EdU cell proliferation assay kit (#C0071S, Beyotime Biotech, China), following the manufacturer’s instructions. The colony-formation assay was also performed to assess cell proliferation. 5000 cells were plated per well in triplicate in 6-well plates. The culture medium was changed every 3 days. Once visible clones were observed, each well was washed with PBS three times, fixed with methanol for 30 minutes at room temperature, and then stained with 0.05% crystal violet for 30 minutes. After washing, the colonies were counted and imaged.

According to the published method (18), transwell migration (without Matrigel) and Matrigel (Matrigel, Corning, China) invasion assays were performed to evaluate cell migration and invasion abilities, respectively. Additionally, cell migration was measured using wound healing assays as previously described (18).

A sphere formation assay was performed to assess the stemness properties of LM3 and MHCC-97H cells. 1×10⁵ cells were seeded into the 6-well ultra-low attachment plates (Corning, China) in sphere formation medium (#CCM0012, Minneapolis, USA). The cells were incubated in a CO₂ incubator for two weeks, and the number of spheres was counted under a stereomicroscope (Olympus).

Statistical analysis was performed using R Studio (R version 4.2.3) and GraphPad Prism 9.0.2. For genes with multiple probes, the maximum expression was selected. Cox regression analysis was conducted using the “survival” package to assess the association between variables and prognosis, including hazard ratios and 95% confidence intervals (CI). Lasso analysis was employed as a variable selection method to refine the scope of variables, with the lambda value chosen for optimal regularization. All data were presented as the mean and standard error of the mean (mean ± SD, n = 3). GraphPad Prism was used to create bar graphs. A two-tailed Student’s t-test was used to compare the means between two groups, and an ANOVA test was used to assess significant differences among various experimental groups. The p-values in multiple comparisons were adjusted to control the false discovery rate (FDR) using the Benjamini-Hochberg method. The OS was evaluated using a Kaplan-Meier (K-M) curve, with statistical significance assessed using a log-rank test. The correlation between two variables conforming to a normal distribution was calculated employing the Pearson method. P< 0.05 was considered as statistically significant.

HCC raw datasets were obtained from the TCGA, ICGC, and GEO databases. Prior to analysis, these datasets were normalized using the log2(TPM + 1) transformation. 99 genes within CCL18 signaling pathway were downloaded from WikiPathways ( Table S1 ). The expression patterns of these genes in the TCGA_LIHC cohort (n = 424) were illustrated using a heatmap ( Figure S1A ). Next, 39 genes exhibiting significant differences in hazard ration (HR) were identified from the TCGA_LIHC cohort by univariate Cox regression analysis ( Figure 1A ). Six hub genes, namely BMI1, CCR3, CDC25C, CFL1, LDHA, and RAC1, were identified using Lasso regression analysis ( Figures 1B, C ). Then, we compared the expression of these six hub genes in TCGA_LIHC, ICGC_JP (n = 445), and GSE14520 (n = 445) cohorts between normal and tumor samples ( Figures 1D–I ). The high expression of CCR3 and LDHA in tumor tissue is not prominent. This study shows that BMI1, CCR3, CDC25C, CFL1, and RAC1 consistently have higher expression levels, while LDHA only shows increased expression in ICGC_JP, not in TCGA_LIHC and GSE14520 cohorts. Additionally, the OS analysis revealed that patients with high expression of 6 hub genes exhibited a shorter survival time ( Figures 1J–O ) in the TCGA_LIHC cohort. Similarly, the Kaplan-Meier Plotter analyzed the correlation between the expression of six hub genes and survival, including OS, RFS, PFS, DSS ( Figures S1B–E ), and the results basically indicated that hub gene expression is associated with a poor prognosis in patients with HCC. On the other hand, the RT-qPCR assay was applied to detect the expression of hub genes in L02 normal hepatocytes and Huh7, HepG2, LM3, and MHCC-97H hepatoma cells ( Figure 1P ). Six hub genes were significantly up-regulated in LM3 and MHCC-97H cell lines, while only BMI1 and CDC25C were up-regulated in hepG2 cells, and CFL1 and LDHA did not show up-regulation in Huh7 cells. Additionally, the expression of hub genes was significantly upregulated upon stimulation with CCL18 in most HCC cells ( Figure 1Q ).

We calculated the risk score of individual HCC patients in all cohorts using the following formula: Risk score = [BMI1 expression × (0.1971976)] + [CCR3 expression × (0.2586062)] + [CDC25C expression × (0.181479)] + [CFL1 expression × (0.2796217)] + [LDHA expression × (0.3388622)] + [RAC1 expression × (0.1223471)]. Applying the median score as the best cut-off value, patients were divided into two groups: the high- and low-risk group. The K-M survival curve analysis revealed that the risk score served as a robust prognostic indicator for HCC patients. Notably, patients with higher risk scores exhibited significantly worse prognosis compared to those in the low-risk group (TCGA_LIHC p<0.001; ICGC_JP p<0.01; GSE14520 p<0.01) ( Figures 2A–C ).

To assess the potential of the risk score as an independent prognostic factor, we conducted univariate and multivariate Cox regression analyses to examine its association with other clinical characteristics, such as age, gender, and stage. We found that age and gender were not identified as independent predictors of prognosis in HCC patients, and the risk score served as an independent risk factor for OS in all cohorts (p<0.01) ( Figures 2D–F ). Risk score triptychs show the corresponding risk score, gene expression, and survival status ( Figures 2G-I ).

Additionally, a predictive nomogram was constructed using the TCGA_LIHC cohort, incorporating the expression levels of the six hub genes. This nomogram provided a quantitative assessment of the 1-, 3-, and 5-year survival rates for each HCC patient ( Figure 2J ), thereby offering potential clinical utility. The ROC curves showed the excellent predictive performance of the hub genes ( Figures 2K–M ), with area under the curve (AUC) values of 0.772 at 1-year, 0.735 at 3-year, and 0.725 at 5-year in the training cohort ( Figure 2K ). Moreover, the calibration curves were drawn to evaluate the consistency between the predictive survival possibility and the actual probability in TCGA_LIHC, ICGC_JP and GSE14520 cohorts ( Figures S2A–H ). These results highlight the remarkable precision and accuracy of the constructed nomogram.

CCL18 exerts distinct effects in various cancer types, demonstrating its ability to enhance the proliferation of specific malignancies, such as ovarian cancer and osteosarcoma. Furthermore, CCL18 acts as an inducer, promoting tumor metastasis¹³. However, the effects of CCL18 on HCC cells have not been deeply studied. Hence, our study aimed to investigate the effects of exogenous CCL18 on HCC cells.

To investigate the biological functions of CCL18 in LM3 and MHCC-97H cells, we stimulated the cells with or without 40ng/mL CCL18 for 48h observed its impact on cell proliferation, migration, invasion and stemness properties. CCK-8 and EdU assays revealed that CCL18 promoted the proliferation abilities of LM3 and MHCC-97H cells ( Figures 3A, B ). Transwell migration and invasion assays showed that CCL18 significantly promoted the migration and invasion abilities of LM3 and MHCC-97H cells ( Figure 3C ). Alterations in the expression of Epithelial-Mesenchymal Transition (EMT)-associated proteins were detected by western blotting, indicating that CCL18 facilitated the EMT process ( Figure 3D ). Cell stemness was determined by the sphere formation assay ( Figure 3E ) and the detection of stemness gene expressions ( Figure 3F ). The results suggest that CCL18 promotes stemness of LM3 and MHCC-97H cells. Taken together, these data collectively indicate that CCL18 may play a promoting role in HCC cell proliferation, migration, invasion, and stemness in vitro.

To investigate the genomic composition of the high-risk group in tumors, we reanalyzed the TCGA_LIHC cohort. Applying the thresholds of p< 0.05 and log2FoldChange > 1, we identified a total of 2722 DEGs in HCC tumor samples ( Figure 4A ) and 475 DEGs in the high-risk group ( Figure 4B ). By intersecting these two gene sets, we identified a total of 308 genes ( Figure 4C ).

Next, to better understand the function and specific mechanism of these genes, we utilized the “clusterProfiler” R package to conduct GO and KEGG enrichment analysis. The analysis of biological processes indicated a significant enrichment of genes specifically expressed in the high-risk group of tumors in processes such as “mitotic nuclear division” and “sister chromatid segregation” ( Figure 4D ). Furthermore, the analysis of cellular components highlighted the enrichment of these genes in cellular locations such as the “chromosomal region” and “spindle” ( Figure 4E ). Moreover, the molecular function analysis demonstrated a significant enrichment of these genes in functions such as “DNA replication origin binding” and “catalytic activity, acting on DNA” ( Figure 4F ). Additionally, the KEGG pathway analysis indicated a significant enrichment of these genes in pathways such as “cell cycle” and “DNA replication” ( Figure 4G ). Furthermore, the interaction network of the 308 proteins was analyzed using the STRING database. Subsequently, the MCODE plugin identified a highly significant cluster consisting of 104 proteins ( Figure 4H ). This observation suggests that these 104 proteins play a crucial role as the main regulatory agents within the high-risk group of HCC.

The direct chemotactic effect of CCL18 on Treg cells and its role in modulating the immunosuppressive tumor microenvironment are well-established (19). Using the ESTIMATE R package, we performed calculations of the relevant indicators for both the high- and low-risk groups. Our results revealed no significant differences in the Stromal Score and ESTIMATE Score values between the high- and low-risk groups. However, a notable disparity was observed in the Immune Score values, indicating a significantly higher level of immune cell infiltration in the tumor samples from the high-risk group ( Figure 5A ). This finding provides evidence that the high-risk group of tumor samples exhibits a greater extent of immune cell infiltration compared to the low-risk group.

Additionally, the CIBERSORT R package was employed to investigate the composition of immune cells and explore the correlation between hub genes and immune cells. The box plots illustrate the estimated proportions of 22 immune cell types in both the high- and low-risk groups ( Figure 5B ). The findings suggest that immune suppression may indeed occur in the high-risk group of HCC. This conclusion is supported by the higher relative abundance of Memory B cells, resting dendritic cells, eosinophils, M0 macrophages, neutrophils, T cells CD4 memory resting, T cells follicular helper, and Tregs observed in this group. Conversely, the low-risk group exhibits a higher relative abundance of M2 macrophages, resting mast cells, activated/resting NK cells, resting T cells CD4 memory, and T cells gamma delta. The scatter plot shows the proportions of the 22 immune cells at the individual patient level ( Figure 5C ). Additionally, the Pearson correlation coefficient was calculated to assess the correlation between hub genes and immune cell types ( Figure 5D ). The Venn diagram illustrates the high expression of 21 immune-related genes in the high-risk group of tumor samples ( Figure 5E ). We computed the correlation between these genes and immune cell types, and the results were similar to previous findings ( Figure 5F ).

Eventually, up-regulation in PD-L1 protein expression was observed in HCC cells after stimulation of CCL18 ( Figure 5G ) or transfection with CCR3, CDC25C, CFL1, LDHA or RAC1 plasmid ( Figure 5H ). Together, these results imply that CCL18 signaling pathway is associated with immune cell infiltration and immune escape in the HCC microenvironment.

We performed an analysis to examine the correlation between six hub genes and various signaling pathways. By applying a significance threshold of p< 0.05 and |cor| > 0.3, we successfully identified specific hub genes associated with different pathways. Within the DNA replication pathway, we found BMI1, CDC25C, and CFL1 to be the hub genes of interest. Similarly, the G2M pathway revealed the presence of BMI1, CDC25C, CFL1, and RAC1 as significant hub genes. Moving to the PI3K pathway, our analysis highlighted BMI1, CCR3, CFL1, and LDHA as the hub genes involved. Lastly, within the EMT pathway, the hub genes CCR3, CFL1, and RAC1 were found to play crucial roles ( Figure 6A ).

In order to elucidate the biological functions of the hub genes, various analyses including proliferation, migration, invasion, and stemness were conducted on HCC cells. Overexpression plasmids of BMI1, CCR3, CDC25C, CFL1, LDHA, and RAC1 were constructed and subsequently transfected into LM3 and MHCC-97H cells. The transfection efficiency was confirmed through qPCR analysis ( Figure 6B ). The CCK-8 and EdU assays demonstrated that all the aforementioned genes played a role in promoting HCC cell proliferation ( Figures 6C, D ). Transwell migration and invasion assays suggested that overexpression of CCR3, LDHA, and RAC1 enhanced the migration and invasion abilities of both LM3 and MHCC-97H cells ( Figures 6E, F ). The results of the sphere formation assay revealed that all of the above genes supported the maintenance of stemness properties in HCC cells ( Figure 6G ), and the detection of stemness marker genes suggested that BMI1, CCR3, CDC25C, and RAC1 maintained the stemness of LM3 and MHCC-97H cells by upregulating stemness transcription factors ( Figure 6H ). Western blot experiments showed that overexpression of BMI1, CCR3, CDC25C, CFL1, LDHA, and RAC1 could induce EMT ( Figure 6I ).

Together, these results imply that different hub genes are involved in different processes that contribute to HCC progression.

To provide better clinical recommendations, we calculated the sensitivity score of common drugs in TCGA_LIHC patients using the “OncoPredict” R package. We then calculated the Pearson correlation between the drug sensitivity score and risk score. Based on a significance level of p<0.05 and a correlation coefficient of cor>0.3 or cor<-0.6, we selected 14 anti-tumor drugs for horizontal lollipop mapping ( Figures 7A–C ). Meanwhile, we utilized the CAMP database to further predict potential small molecule compounds for the treatment of high-risk group patients. Based on the lowest scores, we identified the top 10 ranking compounds as PD-198306, fenretinide, MK-2206, wortmannin, vemurafenib, WYE-125132, BMS-754807, selumetinib, BGT-226 and GSK-269962 ( Figure 7D ). However, it is important to note that further experimental validation and in-depth studies are required to confirm the therapeutic potential and efficacy of these compounds.