GPD2 & Glycolysis in Cholangiocarcinoma: Biomarkers & Immunity

The expression of GPD2 in various cancer types and healthy tissues

GPD2 expression was analyzed across a range of cancer types compared to healthy tissues. It was found that most cancers exhibited increased GPD2 expression (Fig. 1A). However, specific cancer types such as KIRP, KIPAN, WT, READ, TGCT, and ACC demonstrated decreased expression (Fig. 1A). The UALCAN database revealed elevated GPD2 protein levels in multiple cancers including PAAD, COAD, BRCA, UCEC, OV, lung cancer, and LIPH (Fig. 1B). IHC results from the HPA database confirmed high GPD2 protein expression in samples from BRCA, CHOL, and lung cancer (Fig. 1C). Additionally, IF staining localized GPD2 in the mitochondria of human osteosarcoma U2OS cell lines (Fig. 1D), supporting the critical role of GPD2 across various neoplasms.

Association of GPD2 with prognosis and genomic instability

Analysis of TCGA data revealed that GPD2 holds significant prognostic value for numerous cancers (Fig. 2A). Kaplan-Meier curves displayed that increased GPD2 expression were associated with reduced overall survival (OS) in cancers such as LIPH, LUAD, PAAD, STES, and LGG. Conversely, higher GPD2 expression predicted a better prognosis in KIPAN patients (Fig. 2A). Heterogeneity within tumor cells, measured by TMB, MSI, HRD, and LOH, related positively or negatively with GPD2 expression in various contexts. Specifically, positive associations were noted between GPD2 and MSI in STES, STAD, and COAD, but negative correlations were found in BLCA, PRAD, KIPR, and BRCA (Fig. 2B). TMB analysis showed correlations in multiple tumors such as BLCA, COAD, and LUAD among others, while negative correlations appeared in BRCA, GBM, and KIRP (Fig. 2C). The HRD scores correlated positively in LGG, LUAD, SARC, and other cancers (Fig. 2D), and LOH exhibited significant positive correlation across multiple tumors including LGG, LUAD, SARC, and LIHC (Fig. 2E). Stemness, indicated by RNAss and EREG.EXPss tumor stemness scores, positively correlated with GPD2 expression across most tumors, suggesting a role in maintaining tumor cell dedifferentiation status (Fig. 2F). Genetic analysis using cBioPortal mapped GPD2 mutations, identifying higher frequencies particularly in UCEC, SKCM, PRAD, LUSC, and LUAD, primarily alterations in specific mutation types (Fig. 2G). The distribution of these mutations was systematically depicted (Fig. 2H).

Prognostic and Genomic Characteristics. (A) Survival analysis of high- and low-GPD2 in LIPH, LUAD, PAAD, STES, and LGG, (B-F) Associations between GPD2 and MSI, TMB, HRD, LOH and stemness score, (G-H) Genomic alterations and mutation landscape of GPD2 in TCGA pan-cancer

Association of GPD2 and tumor immune infiltration

ESTIMATE scores could assess the abundance of stromal and immune cells, thereby aiding in the evaluation of TME characteristics [22]. Our analysis revealed that GPD2 exhibited correlations with ESTIMATE scores of 17 different tumor types, showing positive correlations with ACC, DLBC, COAD, and COADREAD, while demonstrating negative correlations with other cancers (Fig. 3A). Utilizing the TIMER2.0 database, we investigated the association between GPD2 and immune cell infiltration, which indicated predominantly inverse associations with CD8 + T cells and positive correlations with M2 macrophages across various tumors, as determined by algorithms such as CIBERSORT, CIBERSORT-ABS, and QUANTISEQ (Fig. 3B). Notably, GPD2 displayed positive correlations with CD8 + T cells specifically in COAD, PCPG, DLBC, HNSC, LIHC, TGCT, and UVM. Additionally, single-cell analysis from the TISCH database revealed that GPD2 was primarily expressed in malignant epithelial cells and macrophages (Fig. 3C). Across pan-cancer tissues, GPD2 showed significant positive correlations with numerous immune checkpoint genes (Fig. 3D) and was linked to various immune-related pathways (Fig. 3E).

Association of GPD2 and tumor immune infiltration. (A) The correlation between GPD2 and StromalScore, ImmuneScore, and ESTIMATEScore, (B) Association between GPD2 expression and CD8 + T cells and M2 macrophages across multiple tumors by using CIBERSORT, CIBERSORT-ABS, and QUANTISEQ algorithms, (C) GPD2 expression at the single-cell level in various tumors,(D, E) The correlation of GPD2 expression and immune checkpoint genes and immune-related pathways

Function and drug sensitivity analysis of GPD2

To investigate the function of GPD2 and its relationship with drug sensitivity, we performed enrichment analysis and identified candidate agents using the BEST web application across ESCA, CHOL, and LIPH (Fig. 4A, B, C). Elevated GPD2 expression in ESCA and CHOL was linked to the metabolism-related pathways, while in LIPH, these pathways were predominantly enriched in low GPD2 expression group. Furthermore, we predict the top ten drugs that exhibited either resistance or sensitivity to high GPD2 expression in ESCA, CHOL, and LIPH. We further examined the relationship between increased GPD2 expression and resistance to small molecule drugs across pan-cancer tissues using the CTRP datasets (Fig. 4D).

Function and drug sensitivity analysis of GPD2. (A.B.C) Enrichment analysis and candidate agent analysis by using BEST in ESCA, CHOL, and LIPH, (D) Drug sensitivity analysis of GPD2 in CTRP datasets

The role of GPD2 in suppressing proliferation and glycolytic activity in CCA

To investigate the molecular role of GPD2 in CCA, several in vitro functional assays were performed. Initial observations indicated that GPD2 levels were reduced in Hibepic cell line, in contrast to the higher expression observed in CCA cell lines (Fig. 5A, B). We successfully silenced the GPD2 gene in CCA cells using a lentiviral vector-mediated approach. The effectiveness of GPD2 knockout was subsequently validated by WB analysis (Fig. 5C). Among the shRNA constructs, sh2 and sh3 displayed the most effective knockdown efficiency and were thus used for subsequent experiments. The CCK-8 assay indicated a reduction in CCA cell proliferation following GPD2 knockdown (Fig. 5D). Similarly, the colony formation assay indicated a significant decrease in colony formation and proliferative abilities in CCA cells after GPD2 silencing (Fig. 5E). In addition, to assess the migratory abilities of CCA cells, wound healing and Transwell assays were performed. The results showed a substantial reduction in migration in HUCCT1 and RBE cells following GPD2 knockdown (Fig. 5F, G). Abnormal glucose metabolism provides essential energy for tumor cell growth and proliferation, with tumor cells often generating ATP from glucose through glycolysis, during which pyruvate is converted to lactate. Our findings indicated that silencing GPD2 led to reduced glucose uptake and decreased production of lactate (Fig. 5H). These results suggest that GPD2 knockdown not only inhibits the proliferation but also the glycolytic activity of CCA cells.

Functional validation of GPD2 in CCA cell lines. (A-B) Expression levels of GPD2 mRNA and protein in Hibepic and CCA cell lines, (C) Assessment of GPD2 protein levels in HUCCT1 and RBE cells after GPD2 knockdown using shRNA, (D) Evaluation of cell proliferation in RBE and HUCCT1 post-GPD2 knockdown using the CCK-8 assay, (E) Clonal formation assays were performed in RBE and HUCCT1 cells after GPD2 knockdown, (F-G) Transwell and wound healing assays were conducted to assess the migration abilities in RBE and HUCCT1 cells, (H) Glycolytic activities in RBE and HUCCT1 cells were measured by analyzing glucose consumption and lactate production. Data are presented as means ± SD (*P < 0.05; **P < 0.01; ***P < 0.001; ****P < 0.0001)

The clinical significance of GPD2 in CCA and construction of glycolysis‑related prognostic signature

Our functional studies demonstrated that the potential role of GPD2 in promoting CCA progression through metabolic reprogramming. We next investigated GPD2’s clinical relevance in patient cohorts. Consistent with our experimental findings, GPD2 was significantly overexpressed in tumor tissues across three independent cohorts (Figure S1A), with elevated expression correlating with reduced patient survival (Figure S1B, S1C). Additionally, GPD2 was independently associated with poor prognosis in CCA patients. These clinical observations strongly support our experimental evidence of GPD2’s oncogenic role in CCA.

The connection between GPD2’s functional impact and clinical outcomes was further strengthened by our finding that GPD2 expression positively correlated with glycolysis scores (Fig. 6A). Patients with high GPD2 expression exhibited significantly elevated glycolysis activity (Fig. 6B). This metabolic reprogramming translated to clinical outcomes, as patients with high glycolytic scores showed markedly worse survival (Fig. 6C), mirroring our in vitro results showing GPD2’s critical role in promoting glycolytic metabolism and CCA progression.

Construction of the glycolysis‑related prognostic signature. (A) Correlation between GPD2 expression and glycolysis scores, (B) Comparison of glycolysis scores across groups with low and high GPD2 expression, (C) Survival analysis comparing high and low glycolysis scores, (D-E) Visualization of LASSO regression results, (F.)Survival analysis for low and high expression levels of GPRC5A, HN1, REEP3, and PTBP3, (G) Examination of the association between PTBP3, REEP3, GPRC5A, HN1, GPD2, and glycolysis scores, (H)Multivariate Cox analyses of clinical factors and risk groups

To translate these mechanistic and clinical insights into a practical prognostic tool, we identified a four-gene signature through integrated bioinformatics analysis. Initial screening identified glycolysis-associated genes correlated with GPD2 expression, which were subsequently refined via LASSO regression to yield the optimal prognostic markers: PTBP3, REEP3, GPRC5A, and HN1 (Fig. 6D-E). Notably, elevated expression of these genes impacted patient prognosis adversely (Fig. 6F). The correlations between these four genes, GPD2, and glycolysis scores are presented in Fig. 6G. We then calculated glycolysis-related risk scores using the Cox risk model and classified individuals into low-risk and high-risk categories. Multivariate Cox analyses showed that high-risk group, with the highest hazard ratio (HR) value, was a significant prognostic factor for patients (Fig. 6H). Additionally, vascular invasion and lymphatic metastases were independently associated with poor prognosis.

Validation of the glycolysis‑related prognostic signature

To further assess the prognostic significance of the glycolysis-related risk model, survival analyses were conducted in both training and validation cohorts (Fig. 7A, B). Patients in high-risk category exhibited worse outcomes, with higher scores positively associated with poor prognosis. Analysis using ROC curves identified the model’s predictive ability. The risk model achieved AUC values of 0.78, 0.77, and 0.74 for 1-year, 2-year, and 3-year survival in the training dataset, respectively (Fig. 7C). The AUC values in the validation cohort were 0.72, 0.76, and 0.78, respectively (Fig. 7D). A nomogram incorporating both risk group and clinical features was then developed (Fig. 7E). The model’s predictive performance was further evaluated through calibration curves and C-index calculations. The calibration curve revealed a close alignment between observed and predicted OS (Fig. 7F). Notably, the nomogram exhibited the highest C-index compared to other clinical factors (Fig. 7G).

Validation of the glycolysis‑related prognostic signature. (A-B) Survival analysis of risk model in training and validation cohorts,(C-D) ROC curves demonstrating the accuracy of the risk model in predicting 1-, 2-, and 3-year survival in training and validation cohorts, (E) Development of a nomogram that integrates clinical factors and risk groups, (F) Calibration curves of the nomogram for predicting 1-, 2-, and 3-year OS, (G) Comparison of the C index for the nomogram against various other clinical characteristics

The relationship between risk score and tumor immune infiltration

CIBERSORT was used to evaluate the proportion of infiltrating immune cells, aiming to explore differences in immune infiltration between the high- and low-risk categories (Fig. 8A). The high-risk group showed decreased counts of CD8 + T cells, resting dendritic cells, and resting CD4 memory T cells. Conversely, it exhibited increased levels of M0 macrophages, neutrophils, memory B cells, and activated dendritic cells. Similar results were obtained using the ssGSEA algorithm (Fig. 8B). The correlation analysis also revealed that higher risk scores were linked to elevated proportions of M0 macrophages, neutrophils, activated dendritic cells, and memory B cells, while showing an inverse relationship with CD8 + T cells (Fig. 8C). Additionally, immune infiltration proportions may be associated with OS. For instance, individuals with high proportions of M0 macrophages and neutrophils demonstrated shorter OS (Fig. 8D, E). Multiple immune-related pathways exhibited significant activation in high-risk group (Fig. 8F). Moreover, in the IMvigor210 cohort, our four-gene signature risk model predicted patients’ responses to immunotherapy. The high-risk group displayed an increased rate of complete/partial response (CR/PR), while the low-risk group showed a higher incidence of stable/progressive disease (SD/PD) (Fig. 8G). The risk scores were higher in patients with CR/PR than those with SD/PD (Fig. 8H), implying that individuals in high-risk category may respond more favorably to immunotherapy.

Relationship between risk scores and tumor immune infiltration. (A-B). Comparison of immune cell infiltration between low-risk and high-risk categories as identified by CIBERSORT and ssGSE,(A.C) Correlation between risk scores and immune cell infiltration, (D-E) Survival analysis was conducted by comparing high and low expression levels of M0 macrophages and neutrophils.F. Correlation between risk scores and immune-related pathways, (G) Percentage of patients achieving CR/PR or SD/PD following immunotherapy in low- and high-risk groups within the IMvigor210 cohort. (H) Comparison of risk scores between low- and high-risk groups among CR/PR and SD/PD patients

Function enrichment and drug sensitivity analysis of risk gene set

Utilizing GO and KEGG functional enrichment analyses, we explored key biological functions of the glycolysis-related gene set. There were several main categories in GO-enriched functions, including metabolic processes, hemostasis, response, transport, adhesion, migration, proliferation, and others (Fig. 9A). KEGG analysis revealed that the DEGs between low- and high-risk groups were primarily involved in PI3K-AKT signaling pathway, Glycolysis/Gluconeogenesis, Complement and coagulation cascades, Central carbon metabolism in cancer, Cholesterol metabolism pathways and ECM-receptor interaction (Fig. 9B). We subsequently investigated possible therapeutic targets and compounds for patients with high-risk. IC50 values for various drugs from GDSC database were determined for each sample. The top ten drug candidates with negative correlations between IC50 values and risk scores were identified (Fig. 9C). Figure 9D illustrates the correlation between drug targets and signaling pathways, indicating that EGFR signaling is a potential drug target. In CTRP and PRISM databases, we investigated the correlation between risk scores and the AUC values of various drugs, identifying the top five drug candidates with negative correlations (Fig. 9E, G). The high-risk group displayed reduced AUC values (Fig. 9F, H), implying greater sensitivity to these drug candidates in comparison to the low-risk group.

Function enrichment and drug sensitivity analysis of risk gene set.(A) Clustering of GO terms based on key themes, (B) Enriched pathways identified through KEGG analysis,(C) Candidate compounds identified from the GDSC database,(D) Signaling pathways associated with candidate compounds sourced from GDSC, (E, g) Top five compounds showing the highest negative relationship between risk scores and drug AUC values in the CTPR and PRISM databases, (F, H) Comparison of AUC values between low-risk and high-risk groups

Analysis of risk gene set in single‑cell level

The t-SNE plot illustrated the arrangement of various cell populations. (Fig. 10A). Expression of the four risk genes was examined in single-cell level, revealing their predominant expression in malignant epithelial cells, macrophages, and hepatocytes (Fig. 10B). Using single-cell transcriptomic data, we determined the risk scores for each cell based on the expression of four genes, employing the AUCell method. The risk gene set showed an upregulation trend in malignant epithelial cells and a downregulation trend in plasma cells, T cells, B cells, endothelial cells, and macrophages (Fig. 10C). Among all cell types, the risk gene set was predominantly expressed and activated in malignant epithelial cells, which also had the highest scores (Fig. 10D, E, F). According to the risk scores of malignant cells, they were divided into high- and low-risk categories, followed by enrichment analysis of the DEGs. This revealed notable enrichment primarily in ATP biosynthetic process, ATP metabolic process, oxidative phosphorylation, aerobic respiration, and other pathways (Fig. 10G). We also explored how these malignant cells interacted with various cell types within TME. Our findings showed that high-risk and low-risk malignant cells displayed distinct communication profiles with other cells (Fig. 10H). Malignant cells in high-risk group could impact the output of signaling behaviors such as Midkine (MK), growth differentiation factor (GDF), and macrophage migration inhibitory factor (MIF) pathways, while also influencing input signaling pathways like MK and insulin-like growth factor (IGF) (Fig. 10I). Low-risk malignant cells could disrupt signaling pathways such as MK and vascular endothelial growth factor (VEGF). These signaling pathways influence cellular behaviors such as growth, migration, adhesion, and response to environmental cues, thereby affecting tumor progression.

Analysis of risk gene set in single‑cell level. A. Distribution of various cell types. (B) Proportions of PTBP3, REEP3, GPRC5A, and HN1 expression across different cell types. (C) Trends of AUCell risk scores among different cell types. (D-F) Distribution and expression patterns of AUCell risk scores in various cell types. (G) Functional enrichment analysis for high- and low- risk groups. (H) Ligand-receptor interactions originating from malignant cells. (I) Interactions between various signaling pathways and different cell types with varying risk scores.