A glycolysis-related gene signatures in diffuse large B-Cell lymphoma predicts prognosis and tumor immune microenvironment

1 Introduction

Diffuse large B-cell lymphoma (DLBCL) is the most common histological subtype of lymphoma worldwide, accounting for approximately 25%–30% of all non-Hodgkin lymphoma (Sehn and Salles, 2021). Currently, RCHOP (rituximab, cyclophosphamide, adriamycin, vincristine, and prednisone) regimen is considered the standard first-line chemotherapy (Coiffier et al., 2010). The prognosis of DLBCL patients is often prospectively predicted using the International Prognostic Index (IPI) based on the clinical features before treatment (International Non-Hodgkin’s Lymphoma Prognostic Factors, P, 1993). However, approximately 40% of patients with DLBCL have a poor response to treatment and develop relapsed/refractory disease after first-line therapy (Coiffier et al., 2010). Furthermore, even patients with DLBCL in the same IPI risk category exhibit different clinical outcomes owing to the high tumor heterogeneity (Ruppert et al., 2020). Hence, there is an urgent need to explore novel molecular biomarkers for accurate prediction of prognosis and reveal ideas for developing specific therapeutic targets.

Metabolic reprogramming is one of the most important features of malignant tumors. By altering metabolism and energy production ways, neoplastic cells satisfy the needs of sustainable proliferation, invasion, and metastasis (Faubert et al., 2020). Aerobic glycolysis is among the key features of metabolic reprogramming and is common among cancer cells. In 1924, Otto Heinrich Warburg first demonstrated that tumor cells preferentially converted glucose into lactate via glycolysis even with sufficient oxygen, which differed from normal cells in which glucose is metabolized via oxidative phosphorylation in the mitochondria. This phenomenon is referred to Warburg effect (DeBerardinis and Chandel, 2020). The glucose will not be metabolized to lactate through glycolysis in normal tissues when oxygen is available. Normal tissues convert glucose into lactic acid only under oxygen deprivation. By contrast, tumor tissues metabolize and convert about 66% of ingested glucose to lactate even in the presence of oxygen (Kalyanaraman, 2017; San-Millan and Brooks, 2017). Therefore, through glycolysis, lactate levels in cancer tissue are remarkably increased by 40-fold and are highly related to cancer aggressiveness and poor prognosis. By upregulating VEGF and TGF-β2 expression, inhibiting monocyte migration, suppressing T-cell activation, and promoting the release of IL-6, glycolysis promotes angiogenesis, cell migration, cell metastasis, and immune escape in the tumor microenvironment (San-Millan and Brooks, 2017). Glycolysis contributes about 50%–60% of the total ATP in cancer cells at most (Zheng, 2012). Although glycolysis is less efficient in producing ATP compared to mitochondrial oxidative phosphorylation (OXPHOS), the rate of ATP generation by glycolysis is approximately 100 times faster than that of OXPHOS, which sufficiently meets the energy demand of the rapid tumor proliferation (Ganapathy-Kanniappan and Geschwind, 2013). This demonstrates the advantages of glycolysis over OXPHOS in cancer metabolism. Besides, accumulating evidence has revealed that the tumor immune microenvironment (TIME) influences the efficacy of immunotherapy and prognosis of cancer patients (Karube, 2021; Watanabe, 2021). Apart from the metabolic function, glycolysis-induced hypoxia and lactate accumulation affect the TIME and the function of immune cells. It regulates the antitumor immune response by suppressing lymphocyte proliferation and cytotoxic activity, thereby indirectly promoting cancer progression (Wu et al., 2020). For example, metabolites of glycolysis were reported to attenuate the infiltration level of CD8⁺ T cells and natural killer cells and facilitate infiltration of immunosuppressive T cells.

By providing energy and synthesizing key molecules and enzymes, glycolysis actively regulates the proliferation, metastasis, gene expression regulation, and biosynthesis processes in cancer cells (Faubert et al., 2020). Oslund et al. found that one of the biological functions of bisphosphoglycerate mutase (BPGM) is to regulate serine synthesis by controlling the levels of intermediate glycolytic products (Oslund et al., 2017). Upregulated BPGM has been associated with the poor prognosis of hepatocellular carcinoma patients (Cai et al., 2020). Peptidylglycine alpha-amidating monooxygenase (PAM), a copper- and ascorbate-dependent enzyme, coordinates tissue-specific responses to hypoxia (Rao et al., 2021). Increased expression of PAM has been proven to be an independent predictor of favorable prognosis in patients with neuroendocrine neoplasms (Horton et al., 2020). However, there are no reliable glycolytic biomarkers for predicting the prognosis of DLBCL patients.

In this study, we developed a glycolysis-related prognostic gene signature in DLBCL patients and validated its prediction performance in different cohorts. We then comprehensively analyzed the TIME and immune cell infiltration at different risk levels. This study presents an effective prognostic model and reveals new therapeutic targets for DLBCL that are expected to promote individualized treatment of patients.

2 Methods

2.1 Data acquisition and preprocessing

The gene expression profile and clinicopathological characteristics of patients with DLBCL of GSE181063 (n = 1310), GSE10846 (n = 420), and GSE53786 (n = 119) datasets were downloaded from Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/). The GSE181063 profile was used as the training cohort and based on the Illumina HumanHT-12 WG-DASL V4.0 R2 expression beadchip platform (Illumina, Bethesda, MD, United States). The independent external validation cohorts GSE10846 and GSE53786 were based on the Affymetrix Human Genome U133 Plus 2.0 Array platform (Affymetrix, Santa Clara, CA, United States). The inclusion criteria of DLBCL patients were as follows: 1) histologically confirmed DLBCL; 2) had detailed corresponding clinicopathological information; 3) initially treated with RCHOP or CHOP regimen; 4) overall survival (OS) time ≥30 days. The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

Glycolysis-related gene sets (REACTOME_GLYCOLYSIS, KEGG_GLYCOLYSIS_GLUCONEOGENESIS, HALLMARK_GLYCOLYSIS, BIOCARTA_GLYCOLYSIS_PATHWAY, BIOCARTA_FEEDER_PATHWAY, GOBP_POSITIVE_REGULATION_OF_GLYCOLYTIC_PROCESS, GOBP_NEGATIVE_REGULATION_OF_GLYCOLYTIC_PROCESS, GOBP_GLYCOLYTIC_PROCESS_THROUGH_FRUCTOSE_6_PHOSPHATE) were obtained from the Molecular Signatures Database (MSigDB v7.5.1, http://www.broad.mit.edu/gsea/msigdb/).

2.2 Glycolytic risk model construction and validation

Cox and least absolute shrinkage and selection operator (LASSO) regression analyses were used in the training cohort to identify the glycolysis-related prognostic model. Firstly, we performed the univariate Cox regression analysis with the “survival” package in R to identify survival-associated glycolytic genes. Genes with p < 0.01 were selected. Then, LASSO regression was carried out to further select key glycolytic genes by using the “glmnet” package in R. Finally, the multivariate Cox analysis was performed to establish the glycolytic risk model. The risk model was constructed by each glycolytic gene expression value (G) and the corresponding coefficient (β). The $\mathrm{r}\mathrm{i}\mathrm{s}\mathrm{k}\mathrm{s}\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{e}={\displaystyle \sum }_{n=1}^j\left(\beta \mathrm{j}*\mathrm{G}\mathrm{j}\right)$. Patients with DLBCL were classified into low- and high-risk groups according to the median value of the risk score. Heatmap of the risk gene expression, risk score distribution plots, and survival status scatter plots of the low- and high-risk groups were established to assess the risk model. The predictive accuracy of the risk model was further evaluated by Kaplan-Meier (K-M) analyses and time-dependent receiver operating characteristic (ROC) curves by using the “survival”, “survminer”, and “timeROC” packages in R.

GSE181063 dataset was selected as the training cohort. The other two independent datasets, GSE510846 and GSE53786 were used as external validation cohorts. Additionally, samples in the GSE181063 dataset were randomly divided into two internal validation cohorts at a 1:1 ratio by using the “caret” package in R. Clinical features of the training and validation cohorts are shown in Table 1. The comparison of clinical characteristics between two internal validation cohorts is summarized in Supplementary Table S1. The performance of the risk model was assessed in all validation cohorts. p < 0.05 was considered as the significance threshold.

TABLE 1. Clinicopathological characteristics of patients with DLBCL enrolled from the GEO database.

2.3 Association between the clinical features and glycolytic risk model

To explore whether the risk model could accurately predict prognosis in the subgroup of clinical parameters including age, gender, Ann Arbor stage, ECOG-PS, lactate dehydrogenase (LDH) level, COO (cell of origin), number of extranodal sites, and B symptoms, we compared the survival between different risk groups based on the stratification of clinical parameters. To further figure out the relationship between the risk model and clinical factors, we compared the risk score distribution in subgroups of the clinical parameters by using the “ggpubr” package in R.

2.4 Independent prognostic factors and nomogram construction

Univariate and multivariate Cox regression analyses were applied in the training and validation (GSE10846) cohorts to determine whether the risk score and clinical parameters were independent prognostic factors in DLBCL patients. Next, we developed a prognostic nomogram to integrate independent prognostic factors with “rms”, “foreign”, and “survival” packages in R. According to the risk contribution to OS, prognostic factors were assigned score points. The total score points were calculated to predict the 1-, 3-, and 5-year prognosis of DLBCL patients. The Harrell’s concordance index (C-index) and calibration curves were performed to evaluate the consistency of the actual and predicted survival probabilities of patients by using the R packages “rms”, “foreign”, and “survival”. ROC curves were drawn and the area under the ROC curve (AUC) was calculated to assess the prognostic performance of nomogram and IPI by using the “survivalROC” package in R.

2.5 Tumor microenvironment (TME) and immune cell infiltration

TME consists of tumor cells, immune cells, stromal cells, and non-cellular components (Bader et al., 2020). Immune scores and stromal scores were calculated in each sample using the “estimate” package in R software (Yoshihara et al., 2013). They represented the proportion of stromal and immune cells in TME and were incorporated as ESTIMATE scores, to further infer the tumor purity in each DLBCL patient. The survival was compared between different score levels with a t.test.

To further explore immune cells (encompass activated CD8 T cells, activated CD4 T cells, activated B cells, activated dendritic cells, natural killer cells, and macrophages) infiltration profile in DLBCL patients, the single-sample gene set enrichment analysis (ssGSEA) was performed based on gene expressions of immune cell-specific markers to quantify the enrichment level of each immune cell type with “GSVA” package in R (Bindea et al., 2013). The correlation between the immune cells infiltration and risk model was evaluated with Pearson correlation coefficient and displayed by using the “ggplot2” package in R.

2.6 Association of the immune modulator and glycolytic risk model

Immune checkpoint inhibitors (ICIs) prevent checkpoint proteins from binding with the partner proteins and prompt T cells to kill cancer cells (Galluzzi et al., 2020). The immune checkpoint molecules, including PD-1, PD-L1, PD-L2, CTLA-4, LAG3, TIM-3, TIGIT, B7-H3, B7-H4, and CD47, were considered promising targets in B-cell lymphoma in previous studies (Ansell et al., 2009; Lee et al., 2017; Advani et al., 2018; Xu-Monette et al., 2018; Jiang et al., 2019; Hatic et al., 2021; Shimada et al., 2021; Lee et al., 2022). The association between immune checkpoint molecule expression and risk score were assessed with the Pearson correlation coefficient and presented using the “corrplot” package in R to explore the potential effect of our glycolytic risk model in tumor immunotherapy.

2.7 Measuring the glycolytic risk genes expression at the protein level

The immunohistochemical staining of the glycolytic risk genes in lymphoma patients and normal controls was obtained from the Human Protein Atlas (HPA) database (http://www.proteinatlas.org/) (Liu et al., 2022).

2.8 Analysis of the difference of metabolites

Details of DLBCL cell lines information were downloaded from Cancer Cell Line Encyclopedia (CCLE, https://portals.broadinstitute.org/ccle/data), including cell lines annotations, mRNA expression, and metabolites. The correlation analysis between gene expression and metabolites were performed. Metabolites with p-value < 0.05 and |cor| > 0.5 were selected in this study and shown in the heat map using “ggplot2” package in R software (Sui et al., 2022).

2.9 Evaluation of the sensitivity of chemotherapeutic agents

To predict the drug sensitivity based on the heterogenicity of tumor samples, the half-maximal inhibitory concentration (IC50) of DLBCL patients in high- and low-risk groups was estimated by using the “pRRophetic” package in R. The potential anti-cancer compounds were screened based on the drug response data in Genomics of Drug Sensitivity in Cancer (GDSC) database (Wang et al., 2021).

3 Gene set variation analysis (GSVA)

GSVA between high- and low-risk groups were performed by using the “GSVA package” in R to identify commonly activated/inhibited signaling pathways. The Hallmark gene set obtained from MSigDB was used as the reference gene set. The cut-off criteria were set as the p-value < 0.05 and the t value > 2 (Zhang et al., 2021).

3.1 Statistical analysis

All data analyses were performed in R version 4.0.3. p < 0.05 was considered statistically significant.

4 Results

4.1 Patient characteristics

The study included 947 eligible patients from the training cohort (GSE181063, n = 559) and two external validation cohorts (GSE10846, n = 305; GSE53786, n = 83) (Lenz et al., 2008; Scott et al., 2014; Lacy et al., 2020). Patient baseline characteristics are shown in Table 1. The training cohort was randomly allocated into two internal validation cohorts at a 1:1 ratio. There were no significant differences in clinical factors between the two groups (Chi-square test, p > 0.05). Results are shown in Supplementary Table S1.

4.2 Construction of the glycolytic risk model

We performed the univariate Cox regression analysis with the gene expression profile and survival information of DLBCL patients in the training cohort and identified 14 prognostic genes related to OS (p < 0.01), including eight protective genes (hazard ratio (HR) < 1.0) and six risk genes (HR > 1.0, Figure 1A). Lasso regression (Figures 1B, C) and multivariate analyses (Figure 1D) were further applied to screen prognostic genes. We established a prognostic signature composed of eight glycolysis-related genes in DLBCL, eventually. The glycolytic risk score for OS = (ADH1B × 0.17355739) + (ALDH2 × −0.121597937) + (ANGPTL4 × 0.301463871) + (BPGM × 0.396173449) + (CTH × 0.196458723) + (NUP98 × 0.655960453) + (PAM × −0.14314632) + (PLOD2 × −0.144167963). Among these prognostic genes, ADH1B (p < 0.001), ANGPTL4 (p = 0.030), BPGM (p < 0.001), NUP98 (p < 0.001), and PAM (p = 0.046) were independent risk factors for survival. Based on the risk score, 559 patients were assigned into high-risk (n = 279) and low-risk groups (n = 280) using the median risk score as the cutoff threshold. As revealed by the heatmap, ADH1B, ANGPTL4, BPGM, CTH, and NUP98 were significantly up-regulated while ALDH2, PAM, and PLOD2 were down-regulated in the high-risk group (all p < 0.001, Figure 2A). The risk score distribution and survival status of the individual patient indicated that the poor prognosis was tightly related to the high-risk score (Figure 2A). Patients in the high-risk group had a significantly worse survival than that in the low-risk group (p < 0.001, Figure 2B). The corresponding OS rates at 5-year were 57.2% and 75.7%, respectively. AUCs of the glycolytic risk model were 0.718, 0.695, and 0.688 at 1-, 3-, and 5-year, demonstrating its excellent predictive value (Figure 2C).

FIGURE 1. Identification of glycolysis-related prognostic genes in patients with DLBCL. (A)Fourteen glycolysis-related candidate genes were selected by univariate Cox regression (p < 0.01). (B, C) Further selection of the optimal parameter (lambda) by LASSO Cox regression. (D) Establishment of the prognostic signature with 8 glycolysis-related genes based on multivariate Cox regression analysis. DLBCL, diffuse large B-cell lymphoma; LASSO, least absolute shrinkage and selection operator.

FIGURE 2. Construction of the glycolytic risk model for patients with DLBCL in the training cohort. (A) Glycolysis-related risk gene expression, risk score distribution, and survival status of patients in high- and low-risk groups. (B) Kaplan-Meier analysis of overall survival in different risk groups. (C) Time-dependent ROC analysis of the glycolytic risk model. DLBCL, diffuse large B-cell lymphoma; ROC, receiver operating characteristic.

4.3 Validation of the risk model

To validate the prediction performance of the glycolytic risk model in DLBCL patients, we assessed the gene signature in two external validation cohorts and two internal validation cohorts. Similarly, high-risk patients exhibited a significantly unfavorable prognosis, compared to low-risk patients (Figures 3A, C, E, G; p < 0.001 for internal validation cohorts and GSE10846; p < 0.05 for GSE53786). AUCs at 1-, 3-, 5-year ranged from 0.669 to 0.76 in two internal validation cohorts (Figures 3B, D). The range of AUC values in two external validation cohorts was between 0.62 and 0.70 (Figures 3F, H), which further confirmed the predictive efficacy of the risk model.

FIGURE 3. Validation of the glycolytic risk model. Kaplan-Meier analyses of overall survival for patients in high- and low-risk groups in internal validation cohorts (A, C), GSE10846 (E), and GSE53786 (G). Time-dependent ROC analyses of the risk model in internal validation cohorts (B, D), GSE10846 (F), and GSE53786 (H). ROC, receiver operating characteristic.

4.4 Predictive performance of the glycolytic risk model in clinical subgroups

To assess the predictive performance of the glycolytic risk model in clinical subgroups, we conducted the clinical parameters stratification based on age (≤60 vs. >60), gender (male vs. female), ECOG-PS (<2 vs. ≥2), AnnAnbor Stage (I-II vs. III-IV), LDH (normal vs. elevated), COO (GCB vs. nonGCB), extranodal sites (<2 vs. ≥2), and B symptoms (no vs. yes) (Figures 4A–P). Overall, high-risk patients displayed poorer prognosis compared with low-risk patients regardless of age, gender, ECOG-PS, stage, LDH level, COO, number of extranodal sites, and with or without B symptoms, indicating our risk model was effective and stable.

FIGURE 4. Performance of glycolytic risk model in subgroups of clinical factors. Kaplan-Meier analyses of overall survival for high- and low-risk patients in the subgroups of age (A, B), gender (C, D), ECOG-PS (E, F), stage (G, H), LDH level (I, J), COO (K, L), number of extranodal sites (M, N), and B symptoms (O, P). ECOG-PS, Eastern Cooperative Oncology Group performance status; LDH, lactate dehydrogenase; COO, cell of origin.

4.5 Distribution of glycolytic risk scores in various clinical factors

The distribution of risk scores was analyzed among various clinical factors in the training cohort (Figures 5A–I). Age>60 (p = 0.0017), ECOG-PS≥2 (p < 0.001), advanced stage (p = 0.0038), elevated LDH (p < 0.001), nonGCB (p < 0.001), and with B symptoms (p = 0.019) subgroups contained more high-risk patients compared with age≤60, ECOG-PS<2, early stage, normal LDH, GCB, and without B symptoms subgroups.

FIGURE 5. Distribution of glycolytic risk score in different clinical characteristics. Comparison of the glycolytic risk score distribution between different subgroups of age (A), gender (B), ECOG-PS (C), stage (D), LDH (E), COO (F), rituximab (G), B symptoms (H), and extranodal sites (I). ECOG-PS, Eastern Cooperative Oncology Group performance status; LDH, lactate dehydrogenase; GCB, COO, cell of origin.

4.6 Development and validation of a glycolysis-clinical nomogram for prognostic prediction

By integrating the clinical factors and risk score, we constructed a nomogram to provide a quantitative method for predicting the 1-, 3-, and 5-year individual survival probability in DLBCL patients. Univariate Cox analysis showed that age, ECOG-PS, LDH, stage, rituximab, COO, and glycolytic risk score were OS-related factors in both training and validation cohorts (GSE10846) (Table 2, all p < 0.05). Multivariate Cox analysis revealed that the glycolytic risk score was an independent prognostic indicator in both cohorts (Table 2, all p ≤ 0.002). A higher total score in the nomogram indicated worse survival (Figure 6A). Calibration curves of the nomogram exhibited good agreement between the predicted value and the actual value (C-index = 0.73 for training cohort, Figure 6B; C- index = 0.73 for validation cohort, Figure 6C). In addition, AUCs of the nomogram at 1-, 3-, and 5-year showed its better predictive performance than IPI (Figures 6D–F for training cohort; Figures 6G–I for validation cohort), suggesting the glycolytic risk model improved the predictive value upon existing clinical predictors basis in DLBCL.

TABLE 2. Univariate and multivariate Cox regression analyses of clinicopathological parameters in training and validation cohorts.

FIGURE 6. Construction of prognostic nomogram for overall survival in patients with DLBCL. (A) The nomogram integrated age, ECOG-PS, use of rituximab, COO, and glycolytic risk score to predict 1-, 3- and 5-year survival. (B, C) Calibration plots of 5-year survival probability in training cohort (GSE181063) and validation cohort (GSE10846). ROC curves of the nomogram and IPI in training cohort (GSE181063, (D–F) for AUCs at 1-, 3-, and 5-year, respectively) and validation cohort (GSE10846, (G–I) for AUCs at 1-, 3-, and 5-year, respectively). DLBCL, diffuse large B-cell lymphoma; ECOG-PS, Eastern Cooperative Oncology Group performance status; COO, cell of origin; ROC, receiver operating characteristic; IPI, international prognostic index; AUC, area under the curve.

4.7 Glycolytic risk model and TME

TME was represented by the stromal and immune scores (Liang et al., 2021). Patients in the high-risk group exhibited low immune scores, stromal scores, and ESTIMATE scores (Figures 7A–C; all p < 0.01) and high tumor purities (Figure 7D; p < 0.001), indicating the glycolytic risk model was associated with immunosuppression. The corresponding results of survival analysis (Figures 7E–G, all p ≤ 0.001) were coincident with the high-risk score-related poor prognosis. To explore the underlying mechanism, we analyzed the correlation between the risk model and infiltration levels of immune cells. As shown in Figure 7, the glycolytic risk score was negatively correlated with activated CD8 T cells (Figure 7I; r = −0.12, p = 0.006), activated dendritic cells (Figure 7L; r = −0.11, p = 0.007), natural killer cells (Figure 7M; r = −0.18, p < 0.001), and macrophages (Figure 7N; r = −0.09, p = 0.027), thus might explained why glycolytic risk score had a negative influence on clinical outcomes in DLBLC patients.

FIGURE 7. Association between glycolytic gene signature and tumor microenvironment. (A–D) Comparison of the immune score, stromal score, ESTIMATE score, and tumor purity between low- and high-risk patients with DLBCL. Kaplan-Meier analyses of high- and low-immune score (E), stromal score (F), ESTIMATE score (G), and tumor purity (H) patients. Correlation between activated CD8 T cells (I), activated CD4 T cells (J), activated B cells (K), activated dendritic cells (L), natural killer cells (M), macrophages (N), and risk score. DLBCL, diffuse large B-cell lymphoma.

4.8 Glycolytic risk model and immune checkpoint molecules

Encouraged by plenty of achievements in antitumor immunotherapy, we evaluated the correlations of our glycolytic risk model and the immune checkpoint molecules PD-1, PD-L1, PD-L2, CTLA-4, LAG3, TIM-3, TIGIT, B7-H3, B7-H4, and CD47 in DLBCL. As shown in Figure 8, the glycolytic gene signature was negatively associated with PD-L2 (r = −0.09, p = 0.039), CTLA4 (r = −0.18, p < 0.001), TIM-3 (r = −0.23, p < 0.001), TIGIT (r = −0.09, p = 0.028), B7-H3 (r = −0.09, p = 0.036).

FIGURE 8. Association between the glycolytic risk model and immune checkpoint modulators. (A, B) Correlation between the risk score and immune checkpoint modulators. (B) Correlation between the risk score and expression of PD-L2, CTLA4, TIM-3, TIGIT, and B7-H3.

4.9 The protein expression level of risk genes

To explore the protein expression of the eight glycolytic risk genes in patients with DLBLC, we compared the protein staining intensity of these genes by using the HPA database. ADH1B, CTH, and PLOD2 showed overexpression in the lymph node tissue of non-Hodgkin lymphoma patients compared to normal lymph node tissue. Unfortunately, the information of NUP98 was not available in the HPA database. No expression difference was found of ALDH2, ANGPTL4, BPGM, and PAM between lymphoma patients and control samples (Supplementary Figure S1).

5 Identified differential metabolites

Based on the expression of risk/protective genes in high- and low-risk groups, we explored the associated metabolites and presented the 66 metabolites in the supplementary figure, such as lactate, tryptophan, isoleucine, sarcosine, citrate etc. (Supplementary Figure S2).

5.1 Screening of potential anti-cancer compounds

Based on the GDSC database, we assessed the drug therapeutic responses between high- and low-risk groups. Of note, patients in the high-risk group were more sensitive to cytarabine, methotrexate, vinblastine, gemcitabine, lenalidomide, and rapamycin, while patients in the low-risk group were more sensitive to docetaxel, bortezomib, bleomycin and LFM.A13, which provided the reference directing clinical treatment (Supplementary Figure S3).

5.2 Differential biological behaviors between high- and low-risk groups

GSVA analysis was performed to explore the difference in biological signaling pathways activated/inhibited between high- and low-risk groups. Our results showed significant differences in biological processes between high- and low-risk patients. Notably, MYC targets V2, MYC targets V1, E2F targets, G2M checkpoint, and unfolded protein response signaling were the top enriched signatures in the high-risk group. In comparison, low-risk patients were more relevant to the apical surface, coagulation, epithelial-mesenchymal transition, angiogenesis, and estrogen response early signaling pathways (Supplementary Figure S4).

6 Discussion

Accurate identification of high-risk DLBCL patients is important for prognostic prediction and decision-making regarding treatment choices. The IPI has been used to predict prognosis of patients with aggressive non-Hodgkin lymphoma since 1993 (Liu and Barta, 2019). However, IPI system only includes clinical characteristics and cannot directly reflect the biological heterogeneity of patients with DLBCL (International Non-Hodgkin’s Lymphoma Prognostic Factors, P, 1993). Moreover, there is no conclusive information for comparing the effects of glycolysis and OXPHOS in B-cell lymphoma. Through biological function experiments, the effects of AZD3965 and IACS-010759, the inhibitors of glycolysis and OXPHOS, were tested in 8 cell lines of B cell lymphoma. It was found that AZD3965 decreased the growth of lymphoma cells from 60% to 98% in 4 cell lines and inhibition of oxidative phosphorylation resulted in only 5%–45% inhibition of lymphoma cell growth in all cell lines (Noble et al., 2022). This confirmed that glycolysis regulates the proliferation of B cell lymphoma. In this study, we constructed and validated a glycolytic risk model in patients with DLBCL for the first time. The glycolytic risk model showed good prediction performance for clinical outcomes and characteristics of TIME, metabolites, activated signaling pathways, and the sensitivity of patients to chemotherapy and targeted therapies in DLBCL.

Glycolysis affects tumorigenesis and cancer development via interacting with the TME (Takahashi et al., 2021). Indeed, we found that DLBCL patients with high glycolytic risk scores were characterized by high tumor purities and low immune scores in this study. Sukumar et al. (2013) Showed that increased glycolytic activity downregulated CD8⁺ T cell memory and antitumor function. Previous immune cell enrichment analysis in patients with head and neck squamous cell carcinoma showed that enhanced glycolysis downregulated the infiltration level of T cells, dendritic cells, and B cells, which impaired T cell activation, tumor antigen presentation, and antibody production (Martinez-Reyes and Chandel, 2021). Several glycolysis-related prognostic models have been established to predict cancer survival and immune cell infiltration in the TIME (Yu et al., 2020; Xia et al., 2021). Yu et al. (2020) constructed a glycolysis-based risk scoring model for patients with gastric cancer. They found that the level of immune-positive cells, such as the NK cells, in the low-risk group was significantly higher compared with that of the high-risk group. Particularly, the high-risk group showed high infiltration of immunosuppressive cells (Yu et al., 2020). It was postulated that high glycolytic gene expression induced impaired immunity and immunosuppressive TIME, increased risk of tumor formation, and poor prognosis (Yu et al., 2020). In hepatocellular carcinoma, Xia et al. (2021) confirmed that the glycolysis-related model was negatively correlated with CD8 T cells and NK cells, indicating that the poor prognosis of high-risk group may be driven by infiltration of various immune cells. In accord with them, we found that the glycolytic risk score was negatively correlated to the activated CD8 T cell, natural killer cell, and activated dendritic cell. CD8⁺ T cell is the preferred immune cell type in the process of targeting cancer cells (Farhood et al., 2019). Natural killer cells recognize and kill cancer cells via releasing cytolytic granules. Dendritic cells capture antigens and present them to the T cells during the antitumor process (Galli et al., 2020). These findings might explain the poor prognosis observed in patients of the high-risk group in the present study. In this study, the glycolytic risk score was positively correlated with activated B cells. According to cell-of-origin, DLBCL was pathologically divided into germinal center B cell like, activated B cell like, and unclassifiable subtypes. The activated B cell like subtype exhibited a poor prognosis compared with the germinal center B cell like subtype (Li et al., 2018). The correlation between glycolytic risk score and activated B cells may be explained by the pathological features of DLBCL.

The glycolytic risk model developed in our study showed excellent prediction performance and its efficacy was verified from multiple aspects. High-risk patients had significantly worse clinical outcomes than low-risk patients in all cohorts. Besides, the developed risk model still had good performance after clinical factors stratification. High- and low-risk patients showed similar prognoses in ECOG-PS>2 group as other groups, however, without significant difference, which might be caused by the small number of cases in this group. An important finding is that some specific groups, including age>60, ECOG-PS≥2, advanced stage, elevated LDH, nonGCB, and with B symptoms, encompassed more high-risk patients. Warburg found that cancer cells preferred aerobic glycolysis to oxidative phosphorylation for glucose metabolism (DeBerardinis and Chandel, 2020). However, aerobic glycolysis is inefficient compared to oxidative phosphorylation in terms of generating adenosine 5′-triphosphate. Therefore, cancer cells competitively absorb more glucose than normal cells to meet their proliferation needs (Vander Heiden et al., 2009), and this forms the basis for clinical application of FDG-PET/CT (El-Galaly et al., 2018). It has been reported that aggressive lymphomas have higher maximum standardized uptake value than indolent lymphomas, indicating that the invasion activity of lymphoma is dependent on glucose uptake (El-Galaly et al., 2018). A study by Lee et al. (2021) showed that high-grade renal cell carcinoma had higher F-18 FDG uptake than low-grade tumors and F-18 FDG uptake was positively associated with high Fuhrman grades. These findings suggest that advanced stage and invasive disease status tend to correlate with increased glucose metabolism than early-stage disease and this might explain our results. However, further experiments are needed to verify our hypothesis.

In the past decade, ICIs has revolutionized cancer therapy (Robert, 2020). Several clinical trials have demonstrated that ICIs can effectively treat malignant lymphoma (Li et al., 2021). However, only a proportion of patients show satisfactory response to ICIs. To find effective biomarkers for predicting ICI response, Hu et al. performed a comprehensive analysis at a pan-cancer level and found that low expression of immune checkpoint gene indicated worse response to ICI immunotherapy, low immune infiltration, and poor prognosis (Hu et al., 2021). Ren et al. (2018) reported that higher PD-1 expression was correlated with favorable disease-free survival and OS in breast cancer patients. In our study, the glycolytic gene signature was negatively correlated with immune checkpoint modulator gene expression, including PD-L2, CTLA-4, TIM-3, TIGIT, and B7-H3, which might indicate that high-risk patients with DLBCL have a poor response to the immunotherapy of these targets. However, biomarkers for predicting the response to ICIs are not only based on these target molecules but also on other factors like antigen presentation, tumor mutational burden, and interferon-gamma release (Hatic et al., 2021).

In this retrospective study, we comprehensively analyzed the prognosis and immune environments of patients with DLBCL based on glycolysis. However, our findings are based on the analysis of the online databases and are short of further experimental validation. Further large-scale and prospective cohort studies are thus required to confirm our findings. Moreover, in-depth in vivo and in vitro experiments are necessary to be done to elucidate the biological function of glycolysis-related risk genes and to explore the underlying mechanism.

7 Conclusion

In conclusion, the glycolytic risk model and the nomogram developed in our study provided promising and stable prognostic predictors in patients with DLBCL. Our results indicated that the glycolytic risk model might provide potential therapeutic targets in clinical practice and reflected immune status, giving insights into the underlying mechanism of DLBCL.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

Author contributions

YC and CL conceptualized and designed the study. YC and CL collected the data and carried out the data analysis. YC created the figures and tables. YC drafted the initial manuscript. YC and CL reviewed and revised the manuscript. All authors approved the final submission of the manuscript.

Funding

This work was supported by the National Natural Science Foundation of China (Grant numbers 82200565 and 82200566). Funders were not involved in the study design, data collection and interpretation, manuscript preparation, and decisions of publication.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fcell.2023.1070777/full#supplementary-material

References