Multi-omics characterization of metabolic and immune interactions in prostate cancer

Introduction

Prostate cancer (PCa) is one of the most prevalent malignancies among men globally, accounting for approximately 14.2% of new cancer cases in males and representing a leading cause of cancer-related mortality (1,2). Its clinical significance stems from its highly variable biological behavior, ranging from indolent tumors that may not require immediate intervention to aggressive forms with a substantial risk of metastasis and death (3). Histologically, PCa is predominantly an adenocarcinoma graded by the Gleason/Grade Group system, which remains a key prognostic factor (2,4). Advances in prostate-specific antigen (PSA) screening, multiparametric magnetic resonance imaging (MRI), and targeted biopsy have improved early detection and risk stratification (5,6). Moreover, emerging molecular and genomic biomarkers are being explored to refine diagnosis and prognostication beyond conventional parameters, paving the way for precision diagnostics in PCa (3,7).

Despite these advances, the molecular mechanisms underlying PCa heterogeneity and metabolic reprogramming remain insufficiently defined. Recent progress in single-cell technologies, such as single-cell RNA sequencing (scRNA-seq) and single-cell assay for transposase-accessible chromatin sequencing (scATAC-seq), has facilitated a deeper understanding of the extensive heterogeneity observed in PCa cells, particularly regarding changes in metabolic pathways and the immune microenvironment (8,9). To fully capture the interconnected layers of tumor biology, multiomics integration—combining transcriptomic, epigenomic, and metabolic data—has emerged as powerful tools that provide a more holistic view than single-omics analyses alone (10).

Metabolic reprogramming, including modifications in glucose, lipid, and amino acid metabolism, has been recognized as a fundamental contributor to tumor cell proliferation, survival, and metastatic potential (11-13). In PCa, a prominent reliance on aerobic glycolysis (the Warburg effect) supports rapid biomass synthesis and proliferation, whereas dysregulation of oxidative phosphorylation is increasingly linked to hypoxia adaptation and resistance to androgen deprivation therapy, highlighting these two pathways as central axes of metabolic rewiring (10). Simultaneously, dynamic shifts in immune cell populations within the tumor microenvironment (TME), such as T lymphocytes, macrophages, and natural killer (NK) cells, play a key role in modulating immune evasion and therapeutic resistance (14-16). However, the intricate crosstalk between metabolic reprogramming, particularly in these core pathways, and immune microenvironmental dynamics in PCa has not been fully elucidated.

This study aims to bridge this gap using an integrative multiomics strategy combining scRNA-seq and scATAC-seq to delineate metabolic heterogeneity and immune interactions in PCa, thereby providing new mechanistic and therapeutic insights. We present this article in accordance with the TRIPOD reporting checklist (available at https://tau.amegroups.com/article/view/10.21037/tau-2025-502/rc).

Methods

Data acquisition

The Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/info/datasets.html), developed and maintained by the National Center for Biotechnology Information (NCBI), was used as a repository for gene expression datasets. The GSE168668 dataset, comprising of a complete scRNA-seq expression profile, was obtained from the GEO public database for the purpose of single-cell transcriptomic analysis. Similarly, the GSE168667 dataset, containing a complete scATAC-seq expression profile, was acquired from the same source for chromatin accessibility analysis. The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/), recognized as the largest repository for cancer genomics data, provides access to gene expression profiles, copy number variations, and single-nucleotide polymorphisms, among other genomic features. In this study, messenger RNA (mRNA) expression data from 554 patients with PCa [TCGA-Prostate Adenocarcinoma (PRAD)], including 52 normal and 502 tumor samples, were retrieved from the TCGA public database. The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.

Quality control of scRNA-seq data

The scRNA-seq expression matrix was imported using the Seurat R package (v4.3.0). Cell filtering was conducted based on the total number of unique molecular identifiers (UMIs), the number of detected genes, and the proportion of mitochondrial gene expression. The percentage of mitochondrial gene expression was defined as the proportion of total gene expression attributed to mitochondrial genes. Cells exhibiting elevated mitochondrial gene expression along with low overall RNA expression were considered indicative of apoptotic or dying states. Quality control was implemented using the median absolute deviation (MAD), and cells with values exceeding three MADs from the median were classified as outliers and excluded from further analysis.

Dimensionality reduction, clustering, and annotation of scRNA-seq data

Global normalization was performed using the LogNormalize method. The total expression of each cell was scaled to 10,000 through multiplication by a predefined scaling factor, followed by logarithmic transformation. Cell cycle scores were computed using the CellCycleScoring function. Highly variable genes were identified through the FindVariableFeatures function. To account for variations in UMI counts and mitochondrial gene expression across cells, the ScaleData function was applied. Principal component analysis was conducted on the normalized expression matrix, and selected principal components were retained for subsequent analyses. Nonlinear dimensionality reduction was carried out using Uniform Manifold Approximation and Projection (UMAP). Cell type classification and corresponding marker gene identification were performed using a combination of reference databases, including CellMarker and PanglaoDB, in conjunction with automated annotation provided by the SingleR software.

Processing of scATAC-seq data

The scATAC-seq data were processed using the Signac R package (v1.10.0). Low-quality cells were excluded based on the following criteria: peak region fragments between 3,000< peak region fragments <20,000, reads in peaks >15%, blacklist ratio <0.05, nucleosome signal <4, and transcription start site (TSS) enrichment >2. Peak region fragments were defined as the total number of fragments within identified peak regions, while reads in peaks represented the proportion of fragments located within ATAC-seq peaks. The blacklist ratio corresponded to the proportion of reads mapping to genomic regions identified by ENCODE as blacklist regions. Nucleosome signal referred to the intensity of nucleosome-associated signals, and TSS enrichment denoted the enrichment score at TSSs, as defined by ENCODE. Data normalization was performed using term frequency-inverse document frequency, and all detected peaks were used as input features for dimensionality reduction via Latent Semantic Indexing (LSI). UMAP was subsequently applied for nonlinear dimensionality reduction. Gene activity was inferred using the GeneActivity function, and the resulting gene activity matrices were incorporated into the Seurat object using the CreateAssayObject function.

Integrated analysis of scRNA-seq and scATAC-seq data

To integrate scRNA-seq and scATAC-seq data, the Seurat integration framework was utilized to identify corresponding cell pairs across the two modalities. Shared patterns of gene activity from scATAC-seq and gene expression from scRNA-seq were identified using the FindTransferAnchors function (reduction = ‘cca’). Cell type labels for each cell in the scATAC-seq dataset were predicted using the TransferData function (weight.reduction = ‘lsi’, dims =2:30).

Peak annotation

Peak annotation was conducted to determine the genomic regions associated with the identified peaks, including promoter regions, introns, and exons. The annotatePeak function from the ChIPSeeker R package was employed for this purpose, with promoter regions defined as spanning from −3,000 to +3,000 base pairs when compared to the TSS. Genomic region mapping was conducted using the TxDb.Hsapiens.UCSC.hg38.knownGene transcript object.

Analysis of transcription factor (TF) activity in single cells

To quantify TF motif enrichment, the chromVAR package from Signac (v1.10.0) was used. Cell type-specific peak sets derived from scATAC-seq data were first merged, and motif information from the JASPAR2020 core vertebrate collection was incorporated. Following motif annotation, chromVAR deviation scores were calculated using the RunChromVAR function. Differentially enriched chromVAR motifs were identified using the FindMarkers function.

Pseudotime analysis

Single-cell research has enabled the characterization of complex physiological processes and highly heterogeneous cell populations. Such analyses facilitate the identification of genes specific to particular cell subtypes, as well as those associated with intermediate states and transitions between distinct cellular fates. In many single-cell studies, gene expression is captured asynchronously across individual cells, with each cell representing a snapshot of the transcriptional dynamics under investigation. The Monocle algorithm applies a pseudotime-based strategy to order individual cells along trajectories that reflect underlying biological processes, such as cell differentiation.

Immune infiltration

Cell Type Identification by Estimating Relative Subsets of RNA Transcripts (CIBERSORT) is a widely used method for assessing the composition of immune cell types within the TME. This approach performs deconvolution analysis on immune cell subtype expression matrices by using support vector regression. The method incorporates 547 biomarkers that distinguish 22 human immune cell phenotypes, including various subsets of T cells, B cells, plasma cells, and myeloid cells. In this study, the CIBERSORT algorithm was applied to patient data to estimate the relative proportions of the 22 immune cell types. Correlation analyses were subsequently conducted between gene expression levels and the inferred immune cell composition.

Gene set enrichment analysis (GSEA)

Patients were stratified into high- and low-expression groups according to the expression levels of key genes. GSEA was subsequently conducted to assess pathway-level differences between the two groups. The background gene set was obtained from version 7.0 of the Molecular Signatures Database (MSigDB), which provides annotated gene sets for a wide range of biological pathways. Differential pathway expression between the groups was analyzed, and gene sets demonstrating significant enrichment (adjusted P value <0.05) were ranked based on their consistency scores.

Gene set variation analysis (GSVA)

GSVA is a non-parametric, unsupervised method used to assess gene set enrichment within transcriptomic datasets. This approach assigns enrichment scores to predefined gene sets, thereby converting gene-level variation into pathway-level changes to assess the biological functions of samples. Gene sets were obtained from the MSigDB, and the GSVA algorithm was applied to calculate enrichment scores for each gene set, enabling the assessment of potential functional changes across samples.

Non-coding RNA network related to key genes

The microRNAs (miRNAs), obtained from the mircode database, are small non-coding RNAs that regulate gene expression by facilitating mRNA degradation or inhibiting translation. To explore the potential regulatory roles of miRNAs in the transcriptional or post-transcriptional modulation of key genes, relevant miRNAs were identified using the mircode database. The miRNA-gene interaction network was subsequently visualized using Cytoscape software.

Nomogram model construction and analysis of ENO1/CKB expression in different N stages

To integrate key gene expression with clinical features, we used the TCGA-PRAD cohort (502 tumor samples with complete survival and clinical data). ENO1 and CKB expression were included as continuous variables in the model. Clinical features, including age, tumor stage (T stage), and node stage (N stage), were combined as independent variables, with overall survival as the dependent variable.

A multivariate Cox proportional hazards regression model was constructed. Regression coefficients were converted to a scoring scale, and a nomogram was plotted using the rms R package, with each variable corresponding to a separate tick mark. Total scores were calculated as the sum of individual scores, and 3- and 5-year survival probabilities were derived from the nomogram.

For analysis across N stages, samples were stratified into N0 (no regional lymph node metastasis) and N1 (regional lymph node metastasis) groups. ENO1 and CKB expression data were extracted from TCGA-PRAD RNA-seq and normalized by log₂[fragments per kilobase per million mapped fragments (FPKM) +1]. Differences in gene expression between N0 and N1 groups were assessed using the Wilcoxon rank-sum test, and results were visualized using ggplot2.

Statistical analysis

All statistical analyses were conducted using R software (version 4.3.0). A value of P<0.05 was considered statistically significant.

Results

Single-cell data quality control

From the sample data quality, cells with fewer than 200 detected genes were excluded using the following filtering criteria: (nFeature_RNA >200 & percent.mt ≤ median + 3MAD & nFeature_RNA ≤ median + 3MAD & nCount_RNA ≤ median + 3MAD). In this context, nFeature_RNA denotes the number of detected genes, nCount_RNA represents the total UMIs per cell, and percent.mt indicates the proportion of mitochondrial reads. Following filtration, a total of 2,190 cells were retained. The violin plot and scatter plot illustrating the quality control results are presented in Figure S1A,S1B.

Dimensionality reduction and clustering of single-cell data

A total of 2,000 highly variable genes were identified (Figure S1C). The data were subsequently subjected to normalization, scaling, PCA, and UMAP analysis (Figure 1A, Figure S1D,S1E). Cell subtypes were annotated, and all clusters were categorized into six distinct cell types: LNCaP_1, LNCaP_2, LNCaP_3, LNCaP_4, LNCaP_5, and LNCaP_6 (Figure 1B). A bubble plot of differentially expressed genes (DEGs) across the six cell types is presented in Figure 1C, along with the distribution of cell type proportions presented in Figure 1D.

Figure 1 scRNA-seq landscape of PCa. (A) Cells are clustered into six groups using the UMAP algorithm based on significant principal components from PCA. (B) Cell annotation for the six clusters, with the clusters labeled as LNCaP_1 to LNCaP_6. (C) Dot plot of cell markers for the six cell types. (D) Proportional representation of the six cell types. PCa, prostate cancer; PCA, principal component analysis; scRNA-seq, single-cell RNA sequencing; UMAP, Uniform Manifold Approximation and Projection.

scATAC-seq data processing

Following quality control, a total of 2,892 cells were retained for downstream analysis in the scATAC-seq dataset (Figure 2A). Six clusters were identified using a resolution parameter of 0.8 and dimensions 2 to 30 (Figure 2B). Label transfer from the scRNA-seq dataset was conducted using Seurat’s transfer learning algorithm, resulting in the identification of six cell types in the scATAC-seq dataset that corresponded to those defined in the scRNA-seq data (Figure 2C).

Figure 2 scATAC-seq data processing and integration with scRNA-seq. (A) Quality control of scATAC-seq data, showing the number of cells, genes, and sequencing depth for each sample. (B) Cells are clustered into six groups using the UMAP algorithm based on significant principal components from PCA. (C) Cell annotation for the six clusters, labeled as LNCaP_1 to LNCaP_6. PCA, principal component analysis; scATAC-seq, single-cell assay for transposase-accessible chromatin sequencing; scRNA-seq, single-cell RNA sequencing; TSS, transcription start site; UMAP, Uniform Manifold Approximation and Projection.

Quantification of hotspot mechanisms and differential analysis

A total of 355 glycolysis-related genes were obtained from the literature, and AUCell was applied to quantify their activity in both scRNA-seq and scATAC-seq datasets (17). The results indicated that the LNCaP_4 cell type exhibited the highest glycolysis scores in both datasets (Figure 3A,3B), indicating a strong association with glycolytic metabolism. From the median glycolysis score, LNCaP_4 cells were stratified into two subgroups: Hexp.LNCaP_4 and Lexp.LNCaP_4. DEGs and differential accessibility regions (DARs) were subsequently analyzed between these subgroups. In the scRNA-seq dataset, 37 DEGs were significantly upregulated in Hexp.LNCaP_4 compared to Lexp.LNCaP_4 (logFC >0.25, P_val <0.05, pct.1 >0.1) (Figure 3C, table available at https://cdn.amegroups.cn/static/public/tau-2025-502-1.csv).

Figure 3 Multi-omics identification of key metabolic genes ENO1 and CKB. (A) Glycosylation scores and their differences across the six cell types in scRNA-seq. (B) Glycosylation scores and their differences across the six cell types in scATAC-seq. (C) Volcano plot of differential genes, with blue indicating downregulated genes and red indicating upregulated genes. (D) Intersection of DARs and DEGs, revealing 658 specific DARs and 35 DEGs, with two intersecting genes identified. (E) Correlation between key gene expression and cell developmental trajectories. DAR, differential accessibility region; DEG, differentially expressed gene; FC, fold change; scATAC-seq, single-cell assay for transposase-accessible chromatin sequencing; scRNA-seq, single-cell RNA sequencing; UMAP, Uniform Manifold Approximation and Projection.

Genomic annotation and pathway enrichment of DARs

The 673 DARs were annotated using the ChIPSeeker R package. Most of these regions were located within promoter and intronic regions (Figure S2A). Pathway enrichment analysis indicated that the annotated DARs were predominantly associated with the HIF-1 signaling pathway and fructose and mannose metabolism (Figure S2B). Footprint analysis further demonstrated that Tn5 insertion was enriched at the TF binding sites, with comparable distribution patterns observed between groups (Figure S2C). In the scATAC-seq dataset, 673 DARs were significantly upregulated in Hexp.LNCaP_4 when compared to Lexp.LNCaP_4 (Figure S2D, table available at https://cdn.amegroups.cn/static/public/tau-2025-502-2.csv).

Key gene selection and pseudotime analysis

An intersection analysis of the 37 DEGs and 673 DAR-associated genes identified two overlapping genes: CKB and ENO1 (Figure 3D). Pseudotime analysis was conducted on the scRNA-seq dataset using Monocle2 to compute cell similarity and construct a differentiation trajectory. The resulting visualization depicts the cellular developmental process and facilitates the examination of differentiation stages and gene expression dynamics over pseudotime. Pseudotime-based, state-based, and cluster-based color-coded visualizations are depicted in Figure S2E-S2G. Additionally, the expression patterns of CKB and ENO1 throughout the differentiation process were visualized. Both genes exhibited increased expression during the progression of differentiation, followed by a decline at later stages (Figure 3E).

TF regulation analysis

TF activity analysis identified six TFs—HOXC13, HOXA13, HOXD13, CDX2, CDX4, and CDX1—that exhibited specific binding to the DARs. Sequence logo plots depicting the preferred binding motifs of these TFs are presented in Figure 4A.

Figure 4 Functional mechanisms and immune modulation by ENO1 and CKB. (A) TF motif analysis of differentially accessible chromatin regions, showing the motif sites for HOXC13, HOXA13, HOXD13, CDX2, CDX4, and CDX1. (B) Correlation analysis between key gene expression and the levels of immunoinhibitory factors in the TME. (C) KEGG signaling pathways associated with key genes, along with pathway regulation and involved genes. (D) Correlation between key genes and immune cells. KEGG, Kyoto Encyclopedia of Genes and Genomes; NK, natural killer; TF, transcription factor; TME, tumor microenvironment.

Immune infiltration and immune regulatory factors

The TME is composed of fibroblasts, immune cells, extracellular matrix components, various growth factors, inflammatory cytokines, and distinct physicochemical properties, all of which exert significant influence on disease diagnosis, survival outcomes, and sensitivity to clinical treatments. Immune infiltration levels and intercellular immune correlations are illustrated in Figure S3A,S3B. When compared to the control group, the disease group demonstrated significantly elevated levels of M1 macrophages, plasma cells, and activated CD4 memory T cells, while exhibiting reduced levels of activated dendritic cells, resting mast cells, and monocytes (Figure S3C).

In addition, correlations between the expression of key genes and various immune components, including immunosuppressive factors, immunostimulatory factors, chemokines, and their receptors, were examined. These findings suggest a strong association between key genes and levels of immune cell infiltration, underscoring their potential roles in modulating the immune microenvironment (Figure 4B, Figure S3D-S3G).

GSEA indicated that ENO1 was significantly enriched in pathways such as nucleotide metabolism, glycolysis/gluconeogenesis, and the pentose phosphate pathway (Figure 4C). CKB was enriched in pathways related to oxidative phosphorylation, arginine and proline metabolism, and glutathione metabolism (Figure 4C).

Further analysis indicated that ENO1 expression was positively correlated with CD8⁺ T cells, regulatory T cells (Tregs), and M0 macrophages, and negatively correlated with resting CD4 memory T cells, monocytes, activated dendritic cells, mast cells, and neutrophils. Expression of CKB was positively associated with CD8 T cells, activated NK cells, M0 macrophages, and resting mast cells, while negative correlations were observed with naïve B cells, memory CD4 T cells, activated CD4 T cells, M1 macrophages, activated dendritic cells, and mast cells (Figure 4D).

Pathway enrichment analysis of key genes

Specific signaling pathways associated with the key genes were further investigated to elucidate their potential molecular roles in disease progression. GSVA indicated enrichment of ENO1 in the MYC_TARGETS_V2 and MYC_TARGETS_V1 pathways (Figure S4A), while CKB was associated with xenobiotic metabolism and estrogen response pathways (Figure S4B). These results suggest that ENO1 and CKB may contribute to disease progression through the regulation of these signaling pathways.

Nomogram model construction and analysis of ENO1/CKB expression in different N stages

The nomogram predicted patient survival based on clinical features and key gene expression (Figure 5A). Higher total scores corresponded to worse outcomes, and CKB contributed substantially to the total points, indicating strong prognostic relevance.

Figure 5 Clinical prognostic value of ENO1 and CKB. (A) Nomogram for predicting 3- and 5-year overall survival. The model integrates age, T stage, N stage, and the expression levels of ENO1 and CKB. The C-index of the model is 0.82. The bottom axes show the predicted probability of survival at 1,825 days (5-year survival) and 1,095 days (3-year survival). (B) Expression analysis of ENO1 and CKB across N stages (*, P<0.05; **, P<0.01; Wilcoxon test). C-index, concordance index; N, node; T, tumor.

Expression analysis across N stages revealed that ENO1 was significantly higher in N1 patients (P<0.01), whereas CKB was higher in N0 patients (P<0.05) (Figure 5B). These results suggest that ENO1 is associated with tumor progression and lymph node metastasis.

Non-coding RNA network and correlation analysis of disease-regulatory genes

Reverse prediction of the key genes was conducted using the mircode database, resulting in the identification of 30 miRNAs and 39 mRNA-miRNA interaction pairs, which were visualized using Cytoscape (Figure S4C). Disease-regulatory genes were retrieved from the GeneCards database (https://www.genecards.org/). Expression differences of the top 20 disease-related genes, selected based on Relevance Score and gene expression data, were analyzed. Significant differential expression between the two patient groups was observed for HOXB13, PTEN, CHEK2, TP53, SPOP, PCA3, RNASEL, ELAC2, PCAT1, PRNCR1, SCHLAP1, PCAT7, TMPRSS2, BRAF, and PCAT14. Correlation analysis between key genes and disease-regulatory genes presented several significant associations, including a positive correlation between ENO1 and CHEK2 (r=0.380), and a negative correlation between CKB and ATM (r=−0.452) (Figure S5).

Disease gene co-expression network and immune-metabolic pathways

Disease-related genes were selected from the GeneCards database (https://www.genecards.org/), and co-expression networks between key genes and disease-related genes were analyzed using correlation-based methods (Figures S6A-S6F,S7A-S7F). AUCell was applied to quantify the activity of immune-metabolic pathway genes at the single-cell level. Bubble plots illustrated variations in the activity of key genes within immune-metabolic pathways. ENO1 and CKB demonstrated higher activity in pathways such as oxidative phosphorylation (Figure S8).

Discussion

PCa is the second most frequently diagnosed malignancy among men worldwide, affecting millions globally (18). Despite the availability of multiple therapeutic modalities such as surgical resection, radiation therapy, androgen deprivation therapy, chemotherapy, immunotherapy, and targeted therapy, the prognosis for patients with advanced or metastatic PCa remains poor (19-21). Clinically, both disease progression and therapeutic response are strongly influenced by molecular-level variations (22). Metabolic reprogramming and changes within the immune microenvironment have been recognized as key contributors to tumor initiation, progression, and immune evasion (23). In this study, scRNA-seq and scATAC-seq were applied to comprehensively analyze the transcriptomic and chromatin accessibility profiles of PCa cells, thereby presenting significant cellular heterogeneity. In addition, key regulatory genes implicated in metabolic reprogramming were identified, along with their associations with immune microenvironmental characteristics.

PCa displays considerable clinical heterogeneity, which is reflected in the diversity of molecular changes and variability in clinical progression. Data from TCGA indicated that PCa can be categorized into seven molecular subtypes based on specific gene fusions (e.g., ERG, ETV1/4, and FLI1) or mutations (e.g., SPOP, FOXA1, and IDH1), accounting for approximately 74% of all cases (22). This heterogeneity is further evidenced by the genomic complexity of PCa, including localized hypermutation events and their association with specific genomic features (24). Previous studies have demonstrated that molecular heterogeneity significantly influences tumor biology and therapeutic response, particularly in multifocal PCa, where both intra- and inter-tumoral heterogeneity present substantial challenges to the clinical application of molecular classification systems (25). Such heterogeneity complicates the prediction of clinical treatment outcomes. In this study, six distinct PCa cell subtypes were identified, each exhibiting considerable variability in gene expression profiles and chromatin accessibility. Notably, substantial differences were observed in the expression of metabolism-related genes across these subtypes.

Metabolic reprogramming is recognized as a hallmark of cancer, enabling tumor cells to modify their metabolic pathways in order to support survival and proliferation (11). These findings further underscore the significance of metabolic reprogramming as a defining characteristic of PCa. In PCa, changes in glucose metabolism frequently manifest as the Warburg effect, wherein cancer cells preferentially use aerobic glycolysis for energy production rather than relying on the more efficient oxidative phosphorylation pathway (26). This metabolic shift facilitates not only the production of adenosine triphosphate (ATP) through enhanced glucose utilization but also the generation of biosynthetic intermediates required for cellular growth and division (27).

Beyond glycolysis, metabolic reprogramming in PCa encompasses other metabolic processes, including lipid and creatine metabolism. Enhanced creatine metabolism has been identified as a key contributor to PCa progression. Experimental evidence indicates that creatine treatment is associated with increased basal respiration in vitro and elevated tumor cell proliferation in vivo (28). These metabolic adaptations not only promote rapid tumor cell proliferation but also modulate immune cell function by altering nutrient availability and signaling pathways within the TME (29,30).

In this study, chromatin accessibility was significantly increased for genes related to glucose metabolism and oxidative phosphorylation in the LNCaP_4 cell subtype, indicating a potentially pivotal role for this subtype in metabolic reprogramming. Notably, high expression levels of ENO1 (a glycolytic enzyme) and CKB (a creatine metabolism-related enzyme) were observed in this subtype, indicating their central roles as key regulators of metabolic processes.

ENO1, a critical enzyme within the glycolytic pathway, catalyzes the conversion of 3-phosphoglycerate to 2-phosphoglycerate and contributes to the regulation of glucose metabolism, thereby supporting tumor cell proliferation and survival (31). This analysis demonstrated that ENO1 expression was significantly elevated in the LNCaP_4 subtype, corresponding to the heightened energy demand and proliferative capacity of tumor cells. The upstream drivers of this expression pattern may involve both hypoxia-induced HIF-1 signaling and transcriptional activation, as suggested by increased chromatin accessibility at its promoter/enhancer regions and the enrichment of the HIF-1 pathway in our differential accessibility analysis (32). Beyond its metabolic role, ENO1 appears to shape a pro-tumorigenic microenvironment by enhancing immune suppression—its strong positive correlation with Tregs supports a role in immune evasion. Moreover, its positive correlation with CHEK2 (r=0.380) implies potential crosstalk between metabolic regulation and DNA damage response, possibly contributing to genomic instability and therapy resistance.

In contrast, CKB plays an essential role in maintaining cellular energy homeostasis. Through its involvement in creatine metabolism, CKB contributes to energy regeneration via the creatine-phosphate shuttle, ensuring a continuous energy supply under conditions of high metabolic stress (33). Upregulation of CKB expression was observed in the LNCaP_4 subtype, proposing a key function in supporting tumor cell energetics. However, its downregulation in advanced disease may result from alterations in oncogenic pathways such as AR signaling or MYC activity, which broadly orchestrate metabolic reprogramming in PCa (34). Functionally, the role of CKB extends beyond metabolism. Its positive correlation with NK cell infiltration suggests a capacity to potentiate anti-tumor immunity, positioning it as a favorable regulator. Conversely, its strong negative correlation with ATM (r=−0.452) indicates a potential inverse relationship with DNA repair activation in advanced disease.

Dynamic changes within the tumor immune microenvironment play a key role in tumor initiation, progression, and immune evasion—consistent with the immune-modulatory roles of ENO1 and CKB observed above. Immune cell infiltration analysis presented significant changes in immune cell distribution in PCa samples, most notably an increase in M1 macrophages, plasma cells, and activated CD4⁺ memory T cells. These findings indicate that the immune microenvironment may shift toward a state that facilitates tumor progression. In contrast, reduced proportions of dendritic cells, resting mast cells, and monocytes were observed, indicative of an immunosuppressive microenvironment. Such changes may be driven by metabolic reprogramming: previous studies have demonstrated that changes in energy metabolism pathways can impact immune cell differentiation and function, thereby contributing to the immunosuppressive state in the TME (35,36). This further supports a model where ENO1 and CKB not only regulate metabolism but also indirectly shape immune responses by influencing the TME landscape.

To translate these molecular and immunological findings into clinical practice, we constructed a Nomogram prognostic model to quantify the association between ENO1/CKB and patient outcomes. Using the TCGA-PRAD cohort, the model integrated ENO1/CKB expression (included as continuous variables) with core clinical parameters—age, T stage (reflecting local tumor invasion), and N stage (indicating lymph node metastasis status). It exhibited robust predictive performance: higher total scores (sum of individual variable scores) correlated significantly with worse 3- and 5-year survival outcomes, and CKB emerged as one of the most influential predictors, underscoring its strong prognostic relevance to disease progression. Complementary stratified analysis across N stages further linked these genes to clinical aggressiveness: ENO1 was significantly upregulated in N1 tumors (with lymph node metastasis, P<0.01), while CKB was higher in N0 tumors (without metastasis, P<0.05). These results not only validate ENO1 and CKB as clinically meaningful markers but also clarify their distinct roles—ENO1 as a specific indicator of metastatic potential, and CKB as a prognostic anchor for overall disease trajectory—effectively bridging molecular mechanisms with key clinical endpoints of PCa.

Collectively, these findings indicate that ENO1 acts as a driver of tumor progression and immunosuppression, while CKB contributes to energy stability and may support immune surveillance. This duality suggests a complex role for CKB: while supporting tumor energetics, its higher expression in less aggressive disease may indicate a context-dependent function related to maintaining energy balance and potentially modulating immune activity. The complementary yet contrasting functions of these two genes highlight their potential as biomarkers for risk stratification and therapeutic targeting in PCa.

Regarding clinical translation, ENO1 is a promising therapeutic target. Recent preclinical work demonstrates that ENO1-specific siRNA delivered via a proteolysis-triggered nanotherapeutic system can suppress tumor growth and remodel the immunosuppressive microenvironment (37), providing proof-of-concept for targeted therapy. This approach aligns with emerging strategies that target metabolic vulnerabilities and microenvironmental dependencies within the complex tumor ecosystem (38). Moreover, recent studies have highlighted the interplay between metabolic reprogramming and epigenetic regulation in PCa, suggesting that combination therapies targeting both ENO1 and key epigenetic regulators may represent a promising avenue for future intervention (39). In contrast, CKB currently lacks direct therapeutic strategies, but its prognostic value could guide patient stratification and treatment decisions, including identifying patients who might benefit from immunotherapy or less intensive interventions. Together, these insights underscore the molecular and clinical heterogeneity of PCa and emphasize the importance of integrating metabolic, immune, and epigenetic perspectives into personalized therapeutic design (38,39).

Additionally, correlations between ENO1-CHEK2 and CKB-ATM suggest potential crosstalk between metabolic regulation and DNA repair pathways, offering avenues for future combination therapies that simultaneously target tumor metabolism and genomic stability.

Taken together, these findings consolidate the dual role of ENO1 and CKB in coordinating tumor metabolism and immune regulation, validate their clinical utility as prognostic and metastatic markers, and highlight their potential to inform personalized, metabolism-informed therapeutic strategies for PCa. This work not only advances our understanding of PCa’s pathobiology but also provides actionable insights for translating molecular findings into clinical practice.

Several limitations should be acknowledged. First, the study primarily relied on data derived from publicly available databases, and future work will validate ENO1 and CKB expression patterns in independent clinical cohorts. Second, the intracellular mechanisms by which these metabolic genes influence immune cell function were not extensively investigated. Therefore, functional experiments will be essential to elucidate their mechanistic roles in regulating tumor metabolism and immune interactions, which could be explored through cell-based assays and in vivo models.

Conclusions

In summary, using single-cell technologies (scRNA-seq and scATAC-seq), we identified ENO1 and CKB as key metabolic regulators in PCa and elucidated their dual roles in linking metabolic reprogramming to immune microenvironment modulation. By integrating these genes into a clinical nomogram (with age, T stage, and N stage) and conducting N stage-stratified analysis, we further confirmed two critical clinical associations: ENO1’s specific link to lymph node metastasis and CKB’s strong prognostic value for patient survival. These findings bridge molecular mechanisms with clinical outcomes, clarifying how ENO1 and CKB influence PCa progression through metabolic and immune pathways, and validate their potential for personalized risk stratification and metabolism-targeted therapy. This work provides a new perspective for understanding PCa’s complex pathogenesis and advancing precision oncology approaches for this disease.

Acknowledgments

We are particularly grateful to all the people who have given us help on our article.

Footnote

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

Peer Review File: Available at https://tau.amegroups.com/article/view/10.21037/tau-2025-502/prf

Funding: This study was supported by the Fundamental Research Program of Shanxi Province (No. 202203021221244).

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://tau.amegroups.com/article/view/10.21037/tau-2025-502/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. The 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. Bray F, Laversanne M, Sung H, et al. Global cancer statistics 2022: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin 2024;74:229-63. [Crossref] [PubMed]
2. Vlajnic T, Bubendorf L. Molecular pathology of prostate cancer: a practical approach. Pathology 2021;53:36-43. [Crossref] [PubMed]
3. Kiełb P, Kowalczyk K, Gurwin A, et al. Novel Histopathological Biomarkers in Prostate Cancer: Implications and Perspectives. Biomedicines 2023;11:1552. [Crossref] [PubMed]
4. Pepe P, Pepe L, Fiorentino V, et al. PSMA PET/CT Accuracy in Diagnosing Prostate Cancer Nodes Metastases. In Vivo 2024;38:2880-5. [Crossref] [PubMed]
5. Pepe P, Pepe L, Fiorentino V, et al. Multiparametric MRI targeted prostate biopsy: When omit systematic biopsy? Arch Ital Urol Androl 2024;96:12992. [Crossref] [PubMed]
6. Fiorentino V, Martini M, Dell'Aquila M, et al. Histopathological Ratios to Predict Gleason Score Agreement between Biopsy and Radical Prostatectomy. Diagnostics (Basel) 2020;11:10. [Crossref] [PubMed]
7. Pepe P, Pepe L, Tamburo M, et al. 68Ga-PSMA PET/CT evaluation in men enrolled in prostate cancer Active Surveillance. Arch Ital Urol Androl 2023;95:11322. [Crossref] [PubMed]
8. Brady L, Nelson PS. RISING STARS: Heterogeneity and the tumor microenvironment in neuroendocrine prostate cancer. J Endocrinol 2023;256:e220211. [PubMed]
9. Bian X, Wang W, Abudurexiti M, et al. Integration Analysis of Single-Cell Multi-Omics Reveals Prostate Cancer Heterogeneity. Adv Sci (Weinh) 2024;11:e2305724. [Crossref] [PubMed]
10. Li X, Han Z, Ai J. Synergistic targeting strategies for prostate cancer. Nat Rev Urol 2025;22:645-71. [Crossref] [PubMed]
11. Faubert B, Solmonson A, DeBerardinis RJ. Metabolic reprogramming and cancer progression. Science 2020;368:eaaw5473. [Crossref] [PubMed]
12. Beloribi-Djefaflia S, Vasseur S, Guillaumond F. Lipid metabolic reprogramming in cancer cells. Oncogenesis 2016;5:e189. [Crossref] [PubMed]
13. Wang Z, Wu X, Chen HN, et al. Amino acid metabolic reprogramming in tumor metastatic colonization. Front Oncol 2023;13:1123192. [Crossref] [PubMed]
14. Jiang W, He Y, He W, et al. Exhausted CD8+T Cells in the Tumor Immune Microenvironment: New Pathways to Therapy. Front Immunol 2020;11:622509. [Crossref] [PubMed]
15. Nakamura K, Smyth MJ. Myeloid immunosuppression and immune checkpoints in the tumor microenvironment. Cell Mol Immunol 2020;17:1-12. [Crossref] [PubMed]
16. Chambers AM, Lupo KB, Matosevic S. Tumor Microenvironment-Induced Immunometabolic Reprogramming of Natural Killer Cells. Front Immunol 2018;9:2517. [Crossref] [PubMed]
17. Yang Y, Yang Y, Liu J, et al. Establishment and validation of a carbohydrate metabolism-related gene signature for prognostic model and immune response in acute myeloid leukemia. Front Immunol 2022;13:1038570. [Crossref] [PubMed]
18. Torre LA, Bray F, Siegel RL, et al. Global cancer statistics, 2012. CA Cancer J Clin 2015;65:87-108. [Crossref] [PubMed]
19. Larsson R, Mongan NP, Johansson M, et al. Clinical trial update and novel therapeutic approaches for metastatic prostate cancer. Curr Med Chem 2011;18:4440-53. [Crossref] [PubMed]
20. Valenca LB, Sweeney CJ, Pomerantz MM. Sequencing current therapies in the treatment of metastatic prostate cancer. Cancer Treat Rev 2015;41:332-40. [Crossref] [PubMed]
21. Turco F, Buttigliero C, Delcuratolo MD, et al. Hormonal Agents in Localized and Advanced Prostate Cancer: Current Use and Future Perspectives. Clin Genitourin Cancer 2024;22:102138. [Crossref] [PubMed]
22. The Molecular Taxonomy of Primary Prostate Cancer. Cell 2015;163:1011-25. [Crossref] [PubMed]
23. Perri F, Della Vittoria Scarpati G, Pontone M, et al. Cancer Cell Metabolism Reprogramming and Its Potential Implications on Therapy in Squamous Cell Carcinoma of the Head and Neck: A Review. Cancers (Basel) 2022;14:3560. [Crossref] [PubMed]
24. Fraser M, Sabelnykova VY, Yamaguchi TN, et al. Genomic hallmarks of localized, non-indolent prostate cancer. Nature 2017;541:359-64. [Crossref] [PubMed]
25. Carm KT, Hoff AM, Bakken AC, et al. Interfocal heterogeneity challenges the clinical usefulness of molecular classification of primary prostate cancer. Sci Rep 2019;9:13579. [Crossref] [PubMed]
26. Levine AJ, Puzio-Kuter AM. The control of the metabolic switch in cancers by oncogenes and tumor suppressor genes. Science 2010;330:1340-4. [Crossref] [PubMed]
27. Schulze A, Harris AL. How cancer metabolism is tuned for proliferation and vulnerable to disruption. Nature 2012;491:364-73. [Crossref] [PubMed]
28. Patel R, Ford CA, Rodgers L, et al. Cyclocreatine Suppresses Creatine Metabolism and Impairs Prostate Cancer Progression. Cancer Res 2022;82:2565-75. [Crossref] [PubMed]
29. Peitzsch C, Gorodetska I, Klusa D, et al. Metabolic regulation of prostate cancer heterogeneity and plasticity. Semin Cancer Biol 2022;82:94-119. [Crossref] [PubMed]
30. Shiao SL, Chu GC, Chung LW. Regulation of prostate cancer progression by the tumor microenvironment. Cancer Lett 2016;380:340-8. [Crossref] [PubMed]
31. Kumar AA. Prostate cancer genotyping for risk stratification and precision treatment. Curr Urol 2024;18:87-97. [Crossref] [PubMed]
32. Dong D, Wang Z, Liu M, et al. Combined SNPs sequencing and allele specific proteomics capture reveal functional causality underpinning the 2p25 prostate cancer susceptibility locus. Nat Commun 2025;16:8950. [Crossref] [PubMed]
33. Teishima J, Hayashi T, Nagamatsu H, et al. Fibroblast Growth Factor Family in the Progression of Prostate Cancer. J Clin Med 2019;8:183. [Crossref] [PubMed]
34. Guo H, Wu Y, Nouri M, et al. Androgen receptor and MYC equilibration centralizes on developmental super-enhancer. Nat Commun 2021;12:7308. [Crossref] [PubMed]
35. Comito G, Ippolito L, Chiarugi P, et al. Nutritional Exchanges Within Tumor Microenvironment: Impact for Cancer Aggressiveness. Front Oncol 2020;10:396. [Crossref] [PubMed]
36. Yang K, Wang X, Song C, et al. The role of lipid metabolic reprogramming in tumor microenvironment. Theranostics 2023;13:1774-808. [Crossref] [PubMed]
37. Zhang SM, Jin XK, Chen H, et al. Proteolysis-triggered RNA Interference for Mitochondrial Iron Dyshomeostasis to Activate Antitumor Immunity in Hepatic Carcinoma. Adv Mater 2025;37:e08457. [Crossref] [PubMed]
38. Pecci V, Troisi F, Aiello A, et al. Targeting of H19/cell adhesion molecules circuitry by GSK-J4 epidrug inhibits metastatic progression in prostate cancer. Cancer Cell Int 2024;24:56. [Crossref] [PubMed]
39. Fiorentino V, Pepe L, Pizzimenti C, et al. PD-L1 expression in prostate cancer and Gleason Grade Group: Is there any relationship? Findings from a multi-institutional cohort. Pathol Res Pract 2025;269:155916. [Crossref] [PubMed]