Single-cell transcriptome analysis identifies neutrophil-related prognostic signatures in hepatocellular carcinoma

Introduction

Hepatocellular carcinoma (HCC) is the most prevalent form of primary liver cancer and ranks as the third leading cause of cancer-related mortality worldwide (1). The incidence of HCC has been rising globally, particularly in Western countries, due to increasing rates of metabolic dysfunction-associated steatotic liver disease (MASLD) and nonalcoholic steatohepatitis (NASH) (2). In the United States, the annual incidence of HCC in patients with NASH ranges from 0.7% to 2.6%, with a notable increase in cases among individuals without cirrhosis (3).

Neutrophils, as key components of the innate immune system, play a dual role in tumor progression. In the context of HCC, neutrophils can exert antitumor effects through direct cytotoxic actions and the release of antimicrobial peptides. Conversely, tumor-associated neutrophils (TANs) can promote tumor progression and resistance to therapies by releasing cytokines and forming neutrophil extracellular traps (NETs) (4). The presence of NETs in the inflammatory microenvironment associated with liver cancer has been shown to promote tumor growth (5). These complex interactions underscore the importance of understanding neutrophil dynamics within the HCC microenvironment.

Advancements in bioinformatics, particularly single-cell RNA sequencing and gene expression profiling, have revolutionized our understanding of tumor biology. These technologies enable the identification of cell-type-specific gene signatures and the dissection of tumor microenvironments at unprecedented resolution. Integrating single-cell data with bulk transcriptomic analyses allows for the development of robust prognostic models that can inform personalized therapeutic strategies (6,7).

Despite these advances, current prognostic tools for HCC, including the Barcelona Clinic Liver Cancer (BCLC) staging system and alpha-fetoprotein (AFP) levels, primarily rely on clinicopathological parameters and fail to capture the molecular heterogeneity of the tumor immune microenvironment (8). Neutrophil-related molecular signatures may bridge this gap by providing cell-type-specific prognostic information that reflects the immunological landscape of HCC. Specifically, characterizing neutrophil-derived transcriptional programs could improve prognostic stratification by identifying patients with immunosuppressive microenvironments who are more likely to experience treatment resistance and poor outcomes. Furthermore, such signatures could guide treatment selection by identifying patients who may benefit from immune-modulating therapies or specific targeted agents that counteract neutrophil-mediated immunosuppression (9,10).

In this study, we aim to comprehensively investigate the role of neutrophils in HCC progression and develop a novel prognostic model based on neutrophil-related gene signatures. By integrating single-cell and bulk transcriptomic data, we seek to elucidate the functional contributions of neutrophils to the HCC microenvironment and identify potential therapeutic targets to improve patient outcomes. We present this article in accordance with the TRIPOD reporting checklist (available at https://jgo.amegroups.com/article/view/10.21037/jgo-2025-1-1010/rc).

Methods

Data acquisition

This study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments. Single-cell data (GSE282701) containing 6 samples with complete expression profiles (3 disease, 3 control) were downloaded from the National Center for Biotechnology Information (NCBI) Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/info/datasets.html). GSE14520 Series Matrix Files with GPL3921 annotation including 221 patients with survival information were obtained from GEO. The Cancer Genome Atlas-Liver Hepatocellular Carcinoma (TCGA-LIHC) processed expression data encompassing 424 patients were downloaded from the TCGA portal (https://portal.gdc.cancer.gov/).

Single-cell quality control

Expression profiles were processed using the Seurat package (v4.3.0) in R (v4.3.0). Cells were filtered based on unique molecular identifier (UMI) counts, expressed gene numbers, and mitochondrial gene expression percentages. Median absolute deviation (MAD) quality control removed cells deviating >3 MAD from median values. DoubletFinder package (v2.0.4) filtered doublets from each sample.

Single-cell dimensionality reduction and cell annotation

Global normalization used the LogNormalize method to scale total expression to 10,000 per cell followed by log transformation. CellCycleScoring function calculated cell cycle scores. FindVariableFeatures function identified highly variable genes. ScaleData function removed expression variation from mitochondrial genes, ribosomal genes, and cell cycle differences. RunPCA function performed linear dimensionality reduction on the expression matrix. Harmony package (v1.0) removed batch effects. RunUMAP function performed nonlinear dimensionality reduction using Uniform Manifold Approximation and Projection (UMAP). Cell annotation utilized Cell Marker and Panglao DB databases combined with SingleR package (v1.6.1) for automated annotation.

Ligand-receptor interaction analysis

CellChat package (v1.4.0) quantitatively inferred intercellular communication networks from single-cell data using normalized expression profiles and cell subtype information. Interaction weights and counts quantified interaction relationship closeness to assess cellular activity and influence in disease states.

Functional enrichment analysis

The clusterProfiler R package (v4.0.5) performed Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) functional annotation of neutrophil differential genes to explore functional correlations. Pathways with p-values and q-values <0.05 were considered significant categories.

Model construction and prognosis

Cox univariate regression selected prognosis-related genes from neutrophil markers. LASSO regression using the glmnet package (v4.1-4) constructed prognostic models with 10-fold cross-validation. Risk scores incorporated gene expression values weighted by estimated regression coefficients from LASSO regression analysis. Median risk scores stratified patients into high-risk and low-risk groups. The survival package (v3.3-1) performed Kaplan-Meier (KM) analysis to assess survival differences using log-rank tests. The survivalROC package (v1.0.3) generated receiver operating characteristic (ROC) curves to evaluate model predictive accuracy.

Immune cell infiltration analysis

The CIBERSORT method using the CIBERSORT package estimated relative proportions of 22 immune cell types from RNA-seq data. This method applies support vector regression (SVR) principles for deconvolution analysis using 547 biomarkers distinguishing immune cell phenotypes including T cells, B cells, plasma cells, and myeloid cell subsets.

Drug sensitivity analysis

The oncoPredict R package predicted chemotherapy sensitivity for tumor samples using the Genomics of Drug Sensitivity in Cancer (GDSC) database with regression methods and 10-fold cross-validation testing regression and prediction accuracy.

Gene set variation analysis (GSVA) and gene set enrichment analysis (GSEA) analysis

GSVA used the GSVA package (v1.40.1) to provide comprehensive scoring of gene sets downloaded from the Molecular Signatures Database (MSigDB) (v7.0). GSEA using the clusterProfiler package compared pathway differences between high- and low-risk groups, with significantly enriched gene sets (adjusted P value <0.05) ranked by enrichment scores.

Nomogram construction

Nomogram models using the rms package (v6.3-0) integrated risk scores and clinical characteristics through multivariable regression analysis, with variable contributions determined by regression coefficient magnitudes for scoring calculation.

Statistical analysis

All statistical analyses used R software (version 4.3.0). The survival package generated Kaplan-Meier survival curves compared with log-rank tests. Cox proportional hazard models using the survival package performed multivariable analyses. Statistical significance was set at P<0.05.

Results

Single-cell data processing and quality control

Following a rigorous and comprehensive quality control process, cells expressing fewer than 200 captured genes were systematically filtered out across multiple samples. The DoubletFinder algorithm was leveraged to accurately identify and effectively remove doublets. This meticulous quality control and doublet-removal procedure ultimately retained 52,841 high-quality cells for subsequent downstream analyses (Figure S1A,S1B). Principal component analysis (PCA) was then conducted. The optimal number of principal components was determined to capture the major sources of variation within the dataset. To facilitate dimensionality reduction, variance analysis was performed. This analysis successfully identified 2,000 highly variable genes, selected based on their standardized variance and average expression levels (Figure S1C). These highly variable genes are crucial for capturing the biological heterogeneity within the data. The integrated data processing workflow, encompassing normalization, scaling, PCA, and Harmony batch correction, was executed with accuracy. This sophisticated workflow demonstrated exceptional efficacy in eliminating batch-specific clustering patterns while preserving essential biological variation (Figure S1D-S1F).

Single-cell clustering and neutrophil identification

UMAP dimensionality reduction revealed distinct cell populations with clear separation boundaries, identifying 18 clusters based on transcriptional similarity (Figure 1A). Using canonical marker gene expression patterns and reference database annotations, these 18 clusters were classified into 10 major cell types: T cells, endothelial cells, hepatocytes, cancer-associated fibroblasts (CAFs), neutrophils, plasma cells, monocytes, fibroblasts, macrophages, and B cells (Figure 1B). To verify the accuracy of these cell type classifications, we analyzed the expression of marker genes across the 10 cell types (Figure 1C). The bubble plot in Figure 1C displays the expression levels of specific marker genes within different cell types, where the size of the bubbles corresponds to the percentage of cells expressing the gene within the respective cell type, and the color gradient represents the average expression level of the gene. The expression patterns of these marker genes were consistent with the known characteristics of the cell types, thereby providing robust evidence for the classification of the cell types.

Figure 1 Single-cell transcriptomic analysis identifies neutrophils and their interactions in the HCC microenvironment. (A) UMAP plot showing 18 distinct clusters identified from 52,841 high-quality cells after dimensionality reduction and clustering analysis. Each dot represents a single cell, colored by cluster assignment. (B) Cell type annotation of the 18 clusters into 10 major cell populations based on canonical marker gene expression and reference database annotation. Colors indicate different cell types: T cells (red), endothelial cells (blue), hepatocytes (green), CAFs (purple), neutrophils (orange), plasma cells (brown), monocytes (pink), fibroblasts (cyan), macrophages (yellow), and B cells (gray). (C) Bubble plot displaying marker gene expression patterns across the 10 identified cell types. Bubble size represents the percentage of cells expressing each gene within the respective cell type, while color intensity indicates average expression level. (D,E) Cell-cell communication analysis using CellChat showing ligand-receptor interaction networks. (D) Network plot illustrating communication patterns, with neutrophils showing active participation in signaling networks with various immune and stromal cell populations. (E) Heatmap of interaction strength between different cell types, with the number of interactions and interaction weights/strength indicated. (F) Heatmap of differentially expressed genes in neutrophils compared to other cell types. Red dots represent significantly upregulated genes (average log₂ fold change >1 and adjusted P value <0.05), establishing these as neutrophil-specific transcriptional signatures for prognostic model development. CAFs, cancer-associated fibroblasts; UMAP, Uniform Manifold Approximation and Projection.

CellChat analysis revealed extensive ligand-receptor interaction networks among different cell populations, demonstrating that neutrophils actively participated in signaling networks with various immune and stromal cell populations within the HCC microenvironment (Figure 1D,1E). Differential expression analysis using the FindAllMarkers function identified 304 genes specifically upregulated in neutrophils (average log2 fold change >1 and adjusted P value <0.05), establishing these as neutrophil-specific transcriptional signatures for prognostic model development (Figure 1F).

Functional enrichment analysis of neutrophil-specific genes

Gene Ontology enrichment analysis of the 304 neutrophil-specific genes revealed significant enrichment in biological processes fundamental to neutrophil function. The most significantly enriched GO terms included cell chemotaxis (P<0.001), neutrophil chemotaxis (P<0.001), and cell-substrate junction organization (P<0.01), confirming the biological relevance of these genes to neutrophil biology (Figure 2A). KEGG pathway analysis demonstrated enrichment in disease-relevant pathways including leishmaniasis (P<0.01), chemokine signaling pathway (P<0.001), and nuclear factor kappa-light-chain-enhancer of activated B cells (NF-κB) signaling pathway (P<0.01), indicating involvement in inflammatory and immune response processes critical to HCC pathogenesis (Figure 2B).

Figure 2 Functional enrichment analysis of neutrophil-specific genes. (A) GO biological process enrichment analysis of 304 neutrophil-specific differentially expressed genes. Bar plot shows the top significantly enriched GO terms ranked by −log₁₀(P value). The most enriched terms include cell chemotaxis (P<0.001), neutrophil chemotaxis (P<0.001), and cell-substrate junction organization (P<0.01). The color gradient represents the adjusted P value, indicating the significance of enrichment. (B) KEGG pathway enrichment analysis displaying disease-relevant pathways. Bar plot presents significantly enriched pathways (P<0.05) including leishmaniasis (P<0.01), chemokine signaling pathway (P<0.001), and NF-κB signaling pathway (P<0.01), indicating involvement in inflammatory and immune response processes critical to HCC pathogenesis. The color gradient represents the adjusted P value, indicating the significance of enrichment. BP, biological process; CC, cell composition; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; MF, molecular function; NF-κb, nuclear factor kappaB.

Prognostic model construction and internal validation

Cox univariate regression analysis identified 75 genes significantly associated with patient overall survival (OS) from the 304 neutrophil differential genes (P<0.05) (Figure 3A). LASSO regression with 10-fold cross-validation selected the optimal λ parameter to prevent overfitting, ultimately identifying 14 genes for the final prognostic model (Figure 3B-3D). The resulting risk score formula incorporated each gene’s expression weighted by its LASSO regression coefficient.

Figure 3 Construction and internal validation of neutrophil-related prognostic signature. (A) Forest plot showing hazard ratios and 95% confidence intervals from Cox univariate regression analysis of 75 prognosis-related genes selected from 304 neutrophil differential genes (P<0.05). (B) LASSO regression coefficient profiles for candidate genes. The vertical line indicates the optimal λ value selected by cross-validation. (C) Cross-validation curve for LASSO regression showing mean cross-validated error versus log(λ). The optimal λ corresponds to the minimum cross-validated error. (D) Coefficient plot of the 14 genes selected for the final prognostic model after LASSO regression. (E,F) Kaplan-Meier survival curves comparing high-risk and low-risk groups in training cohort (E) and testing cohort (F). P values were calculated using log-rank test. High-risk patients (red) showed significantly shorter overall survival compared to low-risk patients (blue) in both cohorts (P<0.001). (G,H) Time-dependent ROC curves evaluating the predictive performance of the prognostic model in training cohort (G) and testing cohort (H). AUC values are shown for 1-, 2-, 3-, 4-, and 5-year survival predictions, with AUC >0.70 indicating good discriminatory ability. AUC, area under the curve; CI, confidence interval; LASSO, least absolute shrinkage and selection operator; ROC, receiver operating characteristic; TCGA, The Cancer Genome Atlas.

Patients were stratified into high-risk and low-risk groups using the median risk score as the cutoff. Kaplan-Meier survival analysis demonstrated statistically significant differences in overall survival between risk groups in both the training cohort (log-rank test, P<0.001) and testing cohort (log-rank test, P<0.001) (Figure 3E,3F). Time-dependent ROC curve analysis evaluated the model’s predictive performance, with area under the curve (AUC) values exceeding 0.70 at multiple time points, indicating good discriminatory ability according to established criteria for prognostic models (Figure 3G,3H).

External validation and model generalizability

The neutrophil-related prognostic signature was applied to the independent GSE14520 external validation cohort (n=221) to assess generalizability across different patient populations and experimental platforms. The external cohort demonstrated significantly different overall survival between high-risk and low-risk groups (log-rank test, P<0.001) (Figure 4A), confirming the model’s predictive capability beyond the discovery dataset. Time-dependent ROC analysis in the external dataset yielded AUC values consistently above 0.60 across multiple time points, indicating moderate discriminatory performance. While these AUC values are lower than those observed in the training cohort, they remain within the acceptable range for externally validated prognostic models and support the model’s potential for further clinical evaluation (Figure 4B).

Figure 4 External validation of the neutrophil-related prognostic signature. Kaplan-Meier survival curve for the external validation cohort GSE14520 (n=221). High-risk patients (red line) demonstrated significantly shorter overall survival compared to low-risk patients (blue line) (log-rank test, P<0.001), confirming the model’s predictive capability across different patient populations and experimental platforms. (B) Time-dependent ROC curves for the external validation cohort showing AUC values consistently above 0.60 across multiple time points (1-, 3-, and 5-year), indicating moderate discriminatory performance that supports the model’s potential for further clinical evaluation. AUC, area under the curve; FPR, false positive rate; ROC, receiver operating characteristic; TPR, true positive rate.

Clinical correlation analysis and independent prognostic validation

Risk score distributions demonstrated significant associations with established clinical prognostic indicators across multiple clinical variables. The risk scores showed statistically significant differences among survival status categories (alive vs. deceased, P<0.001), pathological stages (Stage I–II vs. Stage III–IV, P<0.01), and T classifications (T1–2 vs. T3–4, P<0.01), indicating that the molecular signature effectively captured clinically relevant biological differences in tumor aggressiveness (Figure 5A-5D). Additional clinical variables including age, gender, and other pathological parameters showed varying degrees of association with risk stratification (Figure 5E-5H).

Figure 5 Clinical correlation analysis and independent prognostic validation. (A-H) Box plots showing risk score distributions across different clinical variables. (A) Survival status (alive vs. deceased, P<0.001), (B) Pathological stage (Stage I–II vs. Stage III–IV, P<0.01), (C) T classification (T1–2 vs. T3–4, P<0.01), (D-H) additional clinical variables including age, gender, and other pathological parameters. Box plots display median, quartiles, and outliers. Statistical significance was determined using Wilcoxon rank-sum test or Kruskal-Wallis test as appropriate. (I) Forest plot of univariate Cox regression analysis showing hazard ratios and 95% confidence intervals for clinical variables and risk score. The risk score demonstrated strong prognostic significance (P<0.001). (J) Forest plot of multivariable Cox regression analysis adjusting for conventional clinical variables. The risk score remained an independent prognostic factor (P<0.01) after adjustment for age, gender, pathological stage, and T classification. CI, confidence interval.

Univariate Cox regression analysis confirmed the prognostic significance of the risk score as a continuous variable (P<0.001), demonstrating strong association with overall survival (Figure 5I). Multivariable Cox regression analysis, adjusting for conventional clinical variables including age, gender, pathological stage, and T classification, established the risk score as an independent prognostic factor (P<0.01) (Figure 5J). These findings validated the additive predictive value of the neutrophil-related signature beyond conventional clinical prognostic factors.

Immune microenvironment characterization and drug sensitivity analysis

CIBERSORT deconvolution analysis revealed distinct immune cell infiltration profiles between risk groups, with statistically significant differences observed across multiple immune populations. High-risk patients demonstrated altered proportions of memory B cells, resting dendritic cells, M0 macrophages, neutrophils, activated and resting natural killer (NK) cells, memory CD4⁺ T cell subsets, and regulatory T cells (Tregs) compared to low-risk patients (all P<0.05) (Figure 6A,6B). These findings indicated that the neutrophil-related signature captured broader immune dysfunction patterns extending beyond neutrophil biology alone.

Figure 6 Analysis of immune cell infiltration and functional pathways associated with risk stratification in HCC patients. (A) Stacked bar chart showing the distribution of immune cell infiltration profiles across high-risk (Hrisk) and low-risk (Lrisk) groups. The analysis was performed using CIBERSORT deconvolution, revealing distinct immune cell profiles between the two groups. (B) Box plots illustrating the expression levels of various immune cell types, including memory B cells, resting dendritic cells, M0 macrophages, neutrophils, activated and resting NK cells, memory CD4+ T cell subsets, and regulatory T cells, across high-risk and low-risk groups. (C) Violin plots displaying the expression of six compounds (Linsitinib_1510, Cyclophosphamide_1512, Sapitinib_1549, Uprosertib_1553, LCL161_1557, and Lapatinib_1558) across high-risk and low-risk groups. Differential therapeutic responses were identified (P<0.05), suggesting potential utility for guiding personalized therapeutic selection. (D) Bar chart from GSEA showing significantly enriched pathways in high-risk patients, including PI3K/AKT/mTOR signaling (P<0.001), mTORC1 signaling (P<0.01), DNA repair pathways (P<0.01), and cell cycle regulation (normalized enrichment score >1.5, P<0.01). (E) Heatmap of pathway activity scores derived from GSVA comparing high-risk and low-risk groups, highlighting the differential enrichment of cellular pathways. (F) Pie chart showing the distribution of activities in the cell cycle, DNA replication, and Notch signaling pathways between high-risk and low-risk groups. (G) Stacked bar charts showing the proportion of patients with predicted responses (True) and non-responses (False) to immune checkpoint inhibitors, stratified by high-risk (n=182) and low-risk (n=183) groups. The overall chi-square test indicated a significant difference between the two groups (χ²=11.16, P=8.36×10⁻⁴). ns, no significance; *, P<0.05; **, P<0.01; ***, P<0.001. GSEA, gene set enrichment analysis; GSVA, gene set variation analysis; HCC, hepatocellular carcinoma; NK, natural killer.

Correlation analysis using the TISIDB database revealed associations between risk stratification and various immunoregulatory factors, including immunosuppressive factors, immunostimulatory factors, chemokines, and major histocompatibility complex (MHC) molecules (Figure S2). Drug sensitivity prediction using the GDSC database identified differential therapeutic responses between risk groups for six compounds: Linsitinib_1510, Cyclophosphamide_1512, Sapitinib_1549, Uprosertib_1553, LCL161_1557, and Lapatinib_1558 (all P<0.05) (Figure 6C). These differential sensitivity patterns provide preliminary evidence suggesting potential utility for guiding personalized therapeutic selection in precision medicine approaches for HCC treatment.

Pathway enrichment analysis revealed that high-risk patients demonstrated significant enrichment in PI3K/AKT/mTOR signaling (P<0.001), mammalian target of rapamycin complex 1 (mTORC1) signaling (P<0.01), DNA repair pathways (P<0.01), and cell cycle regulation (normalized enrichment score (NES) >1.5, P<0.01) (Figure 6D,6E). Immunotherapy response prediction indicated that high-risk patients exhibited significantly lower predicted responses to immune checkpoint inhibitors (ICIs) compared to low-risk patients (χ² test, P<0.001) (Figure 6F,6G), suggesting potential therapeutic resistance mechanisms.

Integrated clinical model and mechanistic validation

A comprehensive nomogram model was developed by integrating the molecular risk score with clinical variables to enhance prognostic accuracy for clinical application. The nomogram demonstrated that the risk score provided significant independent contributions to survival prediction when combined with conventional clinical parameters (Figure 7A). The model was extended to provide clinically relevant time-specific survival predictions, generating 3- and 5-year survival probability estimates for individual patients (Figure 7B). Performance evaluation through ROC analysis and decision curve analysis (DCA) confirmed superior predictive accuracy of the integrated nomogram compared to individual predictors, with concordance index (C-index) values exceeding 0.75 and positive net benefit across relevant threshold probabilities (Figure 7C,7D).

Figure 7 Development and validation of a prognostic nomogram and single-cell expression analysis. (A) Nomogram integrating clinical variables and molecular risk score for predicting patient survival. The points assigned to each variable correspond to the total points, which are used to estimate survival probabilities. (B) Calibration plots showing the agreement between the nomogram-predicted survival probabilities and observed outcomes for 1-, 3-, and 5-year survival. (C) Time-dependent ROC curves comparing the predictive performance of the integrated nomogram with individual clinical predictors, demonstrating superior accuracy with C-index values above 0.75. (D) Decision curve analysis illustrating the net benefit of the nomogram across various threshold probabilities, indicating a positive clinical utility. (E) Heatmap depicting the expression levels of 14 prognostic signature genes across 10 identified cell types, highlighting neutrophil-specific expression patterns. (F) UMAP plots showing the distribution of 14 signature genes across different immune cell populations, with color gradients representing expression levels, confirming higher expression in neutrophils. ROC, receiver operating characteristic; UMAP, Uniform Manifold Approximation and Projection.

Single-cell expression validation confirmed the neutrophil-specific nature of the prognostic signature genes. Examination of expression patterns across all 10 identified cell types demonstrated that the 14 model genes exhibited significantly higher expression levels in neutrophils compared to other cell populations (P<0.001 for all genes), providing strong biological validation for the neutrophil-specific approach (Figure 7E,7F). AUCell pathway activity analysis revealed that signature genes participated in distinct functional networks: VNN2 in complement and inflammatory responses; UBE2D1, TREM1, SLC11A1, and S100A9 in interferon-gamma (IFN-γ) and tumor necrosis factor alpha (TNF-α)/NF-κB signaling; RTN3 and APLP2 in coagulation and oxidative phosphorylation; ETS2 in angiogenesis; and ETF1, DHX34, colony stimulating factor 3 receptor (CSF3R), and CD58 in inflammatory responses (Figure S3). This validation supported the biological rationale underlying the prognostic model and confirmed the cellular specificity and functional relevance of the identified signature genes.

Discussion

HCC remains a major global health burden, with immune microenvironment dynamics playing a pivotal role in its progression and therapeutic resistance (8). In particular, neutrophils have dual roles in tumor biology. While they can exert anti-tumor effects, their infiltration into tumors often correlates with immunosuppression and poor prognosis in HCC patients (9). Our study provides important insights into the neutrophil-related prognostic signatures in HCC and their implications for both prognosis and treatment.

In this study, we developed a neutrophil-related prognostic signature based on single-cell transcriptomic data, which effectively stratified HCC patients into high- and low-risk groups. The resulting signature was validated across multiple independent cohorts, confirming its robustness and predictive accuracy. High-risk patients, characterized by increased neutrophil infiltration, demonstrated poorer overall survival, supporting the established link between neutrophil infiltration and worse prognosis in HCC. These findings are consistent with previous studies, such as those by Wu et al. (11) and Li et al. (12), which also showed that immune-related signatures are crucial for predicting HCC outcomes. Specifically, Li et al. (12) reported that immune infiltration, particularly neutrophils, correlates with the progression and recurrence of HCC. Furthermore, the role of neutrophils in promoting tumor growth and metastasis through the secretion of pro-inflammatory cytokines and matrix metalloproteinases has been well documented (13,14).

Our immune microenvironment analysis revealed significant differences in the immune cell infiltration profiles between the high-risk and low-risk groups. High-risk patients exhibited higher proportions of immunosuppressive cell populations, including regulatory T cells (Tregs) and M0 macrophages, while anti-tumor immune cells such as cytotoxic T cells were significantly reduced (8,15). These immune infiltration differences have direct prognostic implications. Increased Treg infiltration has been consistently associated with poor OS in HCC, as Tregs suppress cytotoxic CD8⁺ T cell function through the secretion of immunosuppressive cytokines such as interleukin-10 (IL-10) and transforming growth factor-beta (TGF-β), thereby facilitating immune escape (15). Similarly, elevated M0 macrophage proportions may reflect an undifferentiated macrophage pool that is readily polarized toward the immunosuppressive M2 phenotype within the tumor microenvironment, further contributing to immune tolerance and tumor progression (8). In contrast, the depletion of cytotoxic CD8⁺ T lymphocytes in high-risk patients likely diminishes the capacity for effective anti-tumor immune surveillance, which has been reported to be associated with poor survival in HCC (16). Collectively, these immune infiltration patterns suggest that the immunosuppressive microenvironment in high-risk patients creates a vicious cycle of immune evasion and tumor progression, ultimately leading to the observed survival differences between risk groups. This finding aligns with previous studies showing that high neutrophil counts in HCC are associated with an immunosuppressive environment and resistance to immunotherapy (9).

A significant observation in our study was the differential drug sensitivity between the high-risk and low-risk groups. Specifically, high-risk patients were more sensitive to drugs like Linsitinib and Cyclophosphamide, while low-risk patients showed differential sensitivity profiles to other compounds. Mechanistically, Linsitinib targets insulin-like growth factor 1 receptor (IGF-1R) signaling, which has been implicated in hepatic tumor immune microenvironment modulation and the regulation of pro-tumorigenic immune responses (17). Cyclophosphamide, when used at metronomic doses, has been shown to selectively deplete Tregs and restore T cell and NK cell effector functions (18), which may counteract the immunosuppressive environment enriched in high-risk patients. Sapitinib and Lapatinib target epidermal growth factor receptor (EGFR)/human epidermal growth factor receptor 2 (HER2) family members, and EGFR signaling has been shown to mediate tumor immune evasion by upregulating programmed death-ligand 1 (PD-L1) expression and suppressing anti-tumor immune responses (19). Uprosertib inhibits AKT signaling, directly targeting the PI3K/AKT/mTOR pathway that was significantly enriched in high-risk patients. LCL161, a second mitochondria-derived activator of caspases (SMAC) mimetic, promotes tumor cell apoptosis and has demonstrated preliminary clinical activity in solid tumors, with its mechanism of action involving modulation of TNF-dependent signaling pathways relevant to inflammation and immune regulation (20). This finding is critical for precision medicine approaches, as it suggests that neutrophil-related signatures could serve as biomarkers to guide chemotherapy and targeted therapy selection in HCC patients. However, it should be noted that these drug sensitivity predictions are derived from computational analyses using the GDSC database and should be considered hypothesis-generating findings that require prospective clinical validation before therapeutic recommendations can be made. These findings are further corroborated by a series of studies conducted by Geh et al. (10) and Wang et al. (4). Their research has demonstrated that diverse immune microenvironments substantially influence the sensitivity of cancer cells to therapeutic interventions, encompassing both chemotherapy and targeted therapies (4,21). Moreover, differential drug sensitivity patterns, such as those observed with sapitinib and lapatinib, underscore the potential for using immune microenvironment profiling to predict individual responses to treatment and optimize therapeutic strategies.

Additionally, the identification of distinct drug sensitivities between the risk groups supports the notion that targeting specific immune pathways, such as those mediated by neutrophils, could improve therapeutic outcomes in HCC. Yang et al. (22) highlighted the role of the PI3K/AKT/mTOR pathway in HCC progression and its potential as a therapeutic target. In our study, high-risk patients demonstrated increased activity in the PI3K/AKT/mTOR pathway, suggesting the mTOR pathway’s role in conferring chemotherapy resistance and promoting tumor survival (23,24). Therefore, our results could support the design of clinical trials testing combinations of neutrophil-targeting therapies with chemotherapy or mTOR inhibitors.

Importantly, the enriched pathways identified in high-risk patients—PI3K/AKT/mTOR, NF-κB, DNA repair, and inflammatory signaling—while commonly observed in cancer prognosis models, have well-established mechanistic connections to neutrophil biology that strengthen the biological rationale of our model. The PI3K/AKT/mTOR pathway is a master regulator of neutrophil survival, migration, reactive oxygen species (ROS) production, and degranulation, and inhibition of PI3Kγ has been shown to impair neutrophil chemotaxis and inflammatory functions (25). NF-κB signaling is central to neutrophil activation, cytokine production, and NET formation (26), which are key mechanisms by which TANs promote HCC progression (5). Furthermore, the enrichment of DNA repair pathways in high-risk patients may reflect a compensatory cellular response to genotoxic stress imposed by neutrophil-derived ROS and reactive nitrogen species within the tumor microenvironment, as myeloid cell-derived oxidants have been demonstrated to induce epithelial mutagenesis (27). These mechanistic connections suggest several testable hypotheses for future investigation. For instance, pharmacological inhibition of PI3K signaling in TANs may reduce their pro-tumorigenic functions and restore anti-tumor immunity in high-risk HCC patients. Similarly, targeting NF-κB-dependent NET formation may synergize with ICIs to overcome immunotherapy resistance in neutrophil-enriched tumors. Additionally, the neutrophil-derived oxidative stress–DNA repair axis may represent a druggable vulnerability in high-risk HCC, where combined inhibition of DNA repair mechanisms and neutrophil recruitment could enhance therapeutic efficacy. These hypotheses provide a framework for future mechanistic studies and clinical trials to establish the causal role of neutrophils in the pathway activation patterns observed in our high-risk group.

The integration of molecular risk scores with clinical data through a nomogram enhances the prognostic accuracy for HCC patients. Our model demonstrated strong performance with high C-index values and positive net benefit in DCA, indicating that the combination of molecular and clinical data offers significant advantages over conventional prognostic methods. The addition of neutrophil-related signatures to clinical parameters can provide a more personalized approach to patient stratification, which is essential for treatment decision-making.

To contextualize the performance and clinical relevance of our model, we compared it with recently published HCC prognostic signatures. Several immune-related prognostic models have been reported in HCC, including tumor microenvironment-related signatures (11) and immune cell-related signatures (12), with AUC values generally ranging from 0.65 to 0.80 in internal validation cohorts. Ferroptosis-related (28) and metabolism-related (29) prognostic models have also reported comparable predictive performance. Our model achieved AUC values exceeding 0.70 in internal validation and >0.60 in external validation, placing it within a broadly comparable performance range. However, the distinctive advantage of our signature lies in its neutrophil-specific biological basis. Unlike general immune-related models that incorporate mixed immune cell markers or cell death-related models that focus on non-immune pathways, our model specifically captures the transcriptional landscape of tumor-associated neutrophils—a dominant innate immune population in the liver that plays critical roles in chronic inflammation-driven hepatocarcinogenesis. This neutrophil-specific perspective provides complementary prognostic information that may be particularly relevant for patients with inflammation-associated HCC and those being considered for immunotherapy. Future studies should directly compare our model with existing signatures in the same patient cohorts to rigorously evaluate its incremental prognostic value and determine whether the neutrophil-specific information provides non-redundant predictive power.

Previous studies have also explored the clinical utility of nomograms in HCC. For instance, Shen et al. (30) constructed a nomogram using immune-related biomarkers and clinical variables, which outperformed traditional models in predicting patient survival and guiding therapy decisions. Our study adds to this body of work by demonstrating that the neutrophil-related signature can serve as an independent prognostic factor, providing valuable insights into tumor aggressiveness and patient outcomes. The utility of our model in predicting survival probabilities and guiding clinical decisions, particularly for high-risk patients, supports its potential integration into clinical practice.

Pathway enrichment analysis of the signature genes revealed significant involvement of inflammatory pathways, particularly TNF-α/NF-κB and PI3K/AKT/mTOR signaling, in high-risk patients. These pathways are crucial for tumor progression and immune evasion in HCC. Our findings align with recent studies that have shown the activation of the PI3K/AKT/mTOR pathway in HCC and its role in resistance to therapy (31,32). Furthermore, the enrichment of TNF-α signaling highlights the central role of chronic inflammation in liver carcinogenesis, which has been well documented in the literature (33).

Our study also identified specific genes involved in neutrophil-mediated inflammatory responses, such as S100A9 and TREM1. These genes are known to enhance inflammation and promote immune evasion, suggesting that their inhibition could improve therapeutic outcomes (34,35). Targeting such genes may offer new therapeutic opportunities for overcoming immune suppression and enhancing the efficacy of existing treatments.

Our findings also have important implications for the use of immunotherapy in HCC. The lower predicted response to ICIs observed in high-risk patients supports the notion that neutrophil-mediated immune suppression contributes to immunotherapy resistance. This finding is consistent with studies showing that high neutrophil infiltration correlates with resistance to ICIs in HCC. Targeting neutrophils or their signaling pathways may enhance the efficacy of immunotherapy, particularly in patients with high-risk profiles.

Regarding the clinical translational potential of the 14-gene prognostic model, we acknowledge that a large gene panel may present challenges for routine clinical testing. However, multi-gene prognostic panels have been successfully translated into clinical practice in other cancer types, such as OncotypeDX (21 genes) and MammaPrint (70 genes) for breast cancer risk assessment (36). With advances in multiplex quantitative polymerase chain reaction (qPCR) and targeted next-generation sequencing (NGS) technologies, a 14-gene panel is technically feasible for clinical implementation. Nevertheless, future studies should explore the possibility of model streamlining through stepwise gene reduction analysis to identify a minimal gene set that retains adequate predictive performance, thereby enhancing clinical applicability, reducing costs, and facilitating broader adoption in resource-limited settings.

Several limitations of this study warrant consideration. First, the single-cell dataset used for neutrophil gene discovery comprised only 6 samples (3 HCC and 3 controls), which may limit the generalizability of the identified neutrophil-specific gene signatures. However, the single-cell data primarily served as a discovery tool for candidate gene identification, and the prognostic model was subsequently constructed and validated using large-scale bulk transcriptomic cohorts (TCGA-LIHC, n=424; GSE14520, n=221), which strengthens the robustness of the prognostic findings. Second, the external validation AUC values (>0.60) indicate moderate discriminatory performance, and further prospective validation in larger, multi-center cohorts is needed to confirm the model’s clinical utility. Third, all drug sensitivity predictions are based on computational analyses using the GDSC database and require experimental and clinical validation before informing treatment decisions. Fourth, the retrospective nature of the study design limits causal inferences, and the relationship between neutrophil-related gene signatures and clinical outcomes should be further investigated through prospective cohort studies and functional experiments. Fifth, the study did not incorporate proteomic or functional validation of the identified signature genes, which would strengthen the mechanistic conclusions. Future studies integrating multi-omics data and experimental validation in independent prospective cohorts are warranted to address these limitations.

Further research is needed to validate these findings in clinical trials, particularly those investigating combinations of neutrophil-targeted therapies with immune checkpoint inhibitors. The identification of neutrophil-related biomarkers could help stratify patients for more personalized treatment approaches, potentially overcoming resistance to immunotherapy and improving patient outcomes.

Conclusions

In conclusion, our study provides a comprehensive neutrophil-related prognostic signature for HCC, offering new insights into the tumor microenvironment and potential therapeutic targets. The integration of this signature with clinical variables through a nomogram enhances prognostic accuracy and provides a personalized approach to treatment decision-making. While the current evidence is based on retrospective cohort validation, our model represents a promising candidate for clinical risk stratification that warrants further prospective validation in independent cohorts. Our findings underscore the critical role of neutrophils in HCC progression and highlight their potential as therapeutic targets, particularly in the context of immunotherapy. These results pave the way for future clinical applications of immune microenvironment profiling in HCC.

Acknowledgments

None.

Footnote

Reporting Checklist: The authors have completed the TRIPOD reporting checklist. Available at https://jgo.amegroups.com/article/view/10.21037/jgo-2025-1-1010/rc

Peer Review File: Available at https://jgo.amegroups.com/article/view/10.21037/jgo-2025-1-1010/prf

Funding: None.

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://jgo.amegroups.com/article/view/10.21037/jgo-2025-1-1010/coif). The authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. This study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Hwang SY, Danpanichkul P, Agopian V, et al. Hepatocellular carcinoma: updates on epidemiology, surveillance, diagnosis and treatment. Clin Mol Hepatol 2025;31:S228-54. [Crossref] [PubMed]
2. Singal AG, Kanwal F, Llovet JM. Global trends in hepatocellular carcinoma epidemiology: implications for screening, prevention and therapy. Nat Rev Clin Oncol 2023;20:864-84. [Crossref] [PubMed]
3. Kim DY. Changing etiology and epidemiology of hepatocellular carcinoma: Asia and worldwide. J Liver Cancer 2024;24:62-70. [Crossref] [PubMed]
4. Wang J, Li H, Liu Y, et al. CD24 recruits tumor-associated neutrophils to promote the progression of hepatocellular carcinoma. J Immunother Cancer 2025;13:e012118. [Crossref] [PubMed]
5. Zhu W, Fan C, Dong S, et al. Neutrophil extracellular traps regulating tumorimmunity in hepatocellular carcinoma. Front Immunol 2023;14:1253964. [Crossref] [PubMed]
6. Huang D, Ma N, Li X, et al. Advances in single-cell RNA sequencing and its applications in cancer research. J Hematol Oncol 2023;16:98. [Crossref] [PubMed]
7. Li Y, Jin J, Bai F. Cancer biology deciphered by single-cell transcriptomic sequencing. Protein Cell 2022;13:167-79. [Crossref] [PubMed]
8. Yin Y, Feng W, Chen J, et al. Immunosuppressive tumor microenvironment in the progression, metastasis, and therapy of hepatocellular carcinoma: from bench to bedside. Exp Hematol Oncol 2024;13:72. [Crossref] [PubMed]
9. Liu S, Zhou Y, Li G, et al. PLAUR(+) Neutrophils Drive Anti-PD-1 Therapy Resistance in Patients with Hepatocellular Carcinoma by Shaping an Immunosuppressive Microenvironment. Adv Sci (Weinh) 2025;12:e07167. [Crossref] [PubMed]
10. Geh D, Leslie J, Rumney R, et al. Neutrophils as potential therapeutic targets in hepatocellular carcinoma. Nat Rev Gastroenterol Hepatol 2022;19:257-73. [Crossref] [PubMed]
11. Wu Y, Zhai Y, Ding Z, et al. Single-cell transcriptomics reveals tumor microenvironment changes and prognostic gene signatures in hepatocellular carcinoma. Int Immunopharmacol 2024;143:113317. [Crossref] [PubMed]
12. Li MT, Zheng KF, Qiu YE. Identification of immune cell-related prognostic genes characterized by a distinct microenvironment in hepatocellular carcinoma. World J Clin Oncol 2024;15:243-70. [Crossref] [PubMed]
13. Yang L, Liu Q, Zhang X, et al. DNA of neutrophil extracellular traps promotes cancer metastasis via CCDC25. Nature 2020;583:133-8. [Crossref] [PubMed]
14. Deryugina EI, Zajac E, Juncker-Jensen A, et al. Tissue-infiltrating neutrophils constitute the major in vivo source of angiogenesis-inducing MMP-9 in the tumor microenvironment. Neoplasia 2014;16:771-88. [Crossref] [PubMed]
15. Shen KY, Zhu Y, Xie SZ, et al. Immunosuppressive tumor microenvironment and immunotherapy of hepatocellular carcinoma: current status and prospectives. J Hematol Oncol 2024;17:25. [Crossref] [PubMed]
16. Zheng C, Zheng L, Yoo JK, et al. Landscape of Infiltrating T Cells in Liver Cancer Revealed by Single-Cell Sequencing. Cell 2017;169:1342-1356.e16. [Crossref] [PubMed]
17. Osher E, Macaulay VM. Therapeutic Targeting of the IGF Axis. Cells 2019;8:895. [Crossref] [PubMed]
18. Ghiringhelli F, Menard C, Puig PE, et al. Metronomic cyclophosphamide regimen selectively depletes CD4+CD25+ regulatory T cells and restores T and NK effector functions in end stage cancer patients. Cancer Immunol Immunother 2007;56:641-8. [Crossref] [PubMed]
19. Concha-Benavente F, Srivastava RM, Trivedi S, et al. Identification of the Cell-Intrinsic and -Extrinsic Pathways Downstream of EGFR and IFNγ That Induce PD-L1 Expression in Head and Neck Cancer. Cancer Res 2016;76:1031-43. [Crossref] [PubMed]
20. Infante JR, Dees EC, Olszanski AJ, et al. Phase I dose-escalation study of LCL161, an oral inhibitor of apoptosis proteins inhibitor, in patients with advanced solid tumors. J Clin Oncol 2014;32:3103-10. [Crossref] [PubMed]
21. Yang C, Yang H, Mao Z, et al. Autophagy inhibitor-sensitized artificially activated neutrophils against hepatocellular carcinoma. Theranostics 2025;15:7197-218. [Crossref] [PubMed]
22. Yang YM, Kim SY, Seki E. Inflammation and Liver Cancer: Molecular Mechanisms and Therapeutic Targets. Semin Liver Dis 2019;39:26-42. [Crossref] [PubMed]
23. Liu Y, Azizian NG, Sullivan DK, et al. mTOR inhibition attenuates chemosensitivity through the induction of chemotherapy resistant persisters. Nat Commun 2022;13:7047. [Crossref] [PubMed]
24. Zhou Y, Pang J, Liu H, et al. Fibronectin type III domain-containing protein 5 promotes autophagy via the AMPK/mTOR signaling pathway in hepatocellular carcinoma cells, contributing to nab-paclitaxel chemoresistance. Med Oncol 2022;40:53. [Crossref] [PubMed]
25. Hawkins PT, Stephens LR, Suire S, et al. PI3K signaling in neutrophils. Curr Top Microbiol Immunol 2010;346:183-202. [Crossref] [PubMed]
26. Papayannopoulos V. Neutrophil extracellular traps in immunity and disease. Nat Rev Immunol 2018;18:134-47. [Crossref] [PubMed]
27. Canli Ö, Nicolas AM, Gupta J, et al. Myeloid Cell-Derived Reactive Oxygen Species Induce Epithelial Mutagenesis. Cancer Cell 2017;32:869-883.e5. [Crossref] [PubMed]
28. Liang JY, Wang DS, Lin HC, et al. A Novel Ferroptosis-related Gene Signature for Overall Survival Prediction in Patients with Hepatocellular Carcinoma. Int J Biol Sci 2020;16:2430-41. [Crossref] [PubMed]
29. Zhao Y, Zhang J, Wang S, et al. Identification and Validation of a Nine-Gene Amino Acid Metabolism-Related Risk Signature in HCC. Front Cell Dev Biol 2021;9:731790. [Crossref] [PubMed]
30. Shen XT, Xie SZ, Zheng X, et al. Cirrhotic-extracellular matrix attenuates aPD-1 treatment response by initiating immunosuppressive neutrophil extracellular traps formation in hepatocellular carcinoma. Exp Hematol Oncol 2024;13:20. [Crossref] [PubMed]
31. Ferrín G, Guerrero M, Amado V, et al. Activation of mTOR Signaling Pathway in Hepatocellular Carcinoma. Int J Mol Sci 2020;21:1266. [Crossref] [PubMed]
32. Kahraman DC, Kahraman T, Cetin-Atalay R. Targeting PI3K/Akt/mTOR Pathway Identifies Differential Expression and Functional Role of IL8 in Liver Cancer Stem Cell Enrichment. Mol Cancer Ther 2019;18:2146-57. [Crossref] [PubMed]
33. Tiegs G, Horst AK. TNF in the liver: targeting a central player in inflammation. Semin Immunopathol 2022;44:445-59. [Crossref] [PubMed]
34. Le Y, Zhang J, Lin Y, et al. S100A9 Exacerbates the Inflammation in Rosacea through Toll-Like Receptor 4/MyD88/NF-κB Signaling Pathway. J Invest Dermatol 2024;144:1985-1993.e1. [Crossref] [PubMed]
35. Colonna M. The biology of TREM receptors. Nat Rev Immunol 2023;23:580-94. [Crossref] [PubMed]
36. Sparano JA, Gray RJ, Makower DF, et al. Adjuvant Chemotherapy Guided by a 21-Gene Expression Assay in Breast Cancer. N Engl J Med 2018;379:111-21. [Crossref] [PubMed]