From the Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo/), expression profiles were downloaded for 4 MM datasets, including GSE136337, GSE24080, GSE4204, and GSE6477. The datasets were log2 transformed after normalising with a microarray muite 5.0 (MAS5.0) method. Detailed clinicopathological data and survival data were also extracted from GEO (Table 1). The GSE136337 was used as the training, while the GSE24080 and GSE4204 datasets were used for validation. And we used the GSE6477 to reveal the evolution of GRGs in normal individuals, MGUS, SMM, and active MM groups. Glycolysis-related genomes “c2.cp.biocarta.v7.2.symbols”, “c2.cp.hallmark.v5.0.symbols”, “c2.cp.reactome.v7.5.symbols”, “c2.cp.kegg.v7.2.symbols”, “c2.cp.wikipathways.v7.2.symbols” in Gene Set Enrichment Analysis (GSEA) (http://www.gsea-msigdb.org/gsea/msigdb) were the source of the GRGs. After taking intersections with the three GEO datasets (GEO136337, GSE24080 and GSE4204), 293 genes were finally used for further study (Supplementary Table S1).

The prognostic glycolytic genes identified by univariable Cox regression analysis (p < 0.001) were then subjected to the Least absolute shrinkage and selection operator (LASSO) Cox regression analysis. Using the R package “glmnet”, the penalty parameter λ was set as 0.09 to determine the best weighting coefficient of glycolytic genes.

The co-expression correlation matrix of the 12 genes was constructed by “ggcorrplot” R package. The protein-protein interaction (PPI) network of the 12 genes was acquired from the STRING database (version 11.5) (https://www.string-db.org/) (Szklarczyk et al., 2011). The genetic mutation landscape of the GRGs in MM was explored with The Cancer Genome Atlas (TCGA) (https://portal.gdc.cancer.gov/) with the “maftools” package and cBioPortal for Cancer Genomics (http://www.cbioportal.org/). We also used the Cancer Cell Line Encyclopedia database (CCLE, https://portals.broadinstitute.org/ccle) to further validate these prognostic genes’ expression.

To further investigate the heterogeneity between subgroups, we compared clinical characteristics and drug sensitivity among the subtypes. As for cytogenetics, based on the existing data in the training dataset and the previous findings, we defined high-risk cytogenetic abnormalities (HRCAs) as follows: del17p, amp1q, t (4; 14), t (14; 20), t (14; 16) or MYC aberrations determined by fluorescent in situ hybridization (FISH) or conventional cytogenetics (Chng et al., 2014; Sonneveld et al., 2016; Abdallah et al., 2020; Schmidt et al., 2021; Wallington-Beddoe and Mynott, 2021). Patients with at least one HRCA were defined as the high-risk cohort. The low-risk subgroup included those with other abnormalities [del13q, del16q, del1p, del1q, t (11; 14), t (12; 14)]. The rest belonged to the non-mutation cohort.

The “pRRophetic” package was used to assess the drug susceptibility between the low- and high-scoring group (Geeleher et al., 2014a; Geeleher et al., 2014b). This method builds the statistical model based on the data from the Cancer Genome Project (CGP) (Garnett et al., 2012), consisting of baseline gene expression microarray data and sensitivity to 138 drugs in a panel of almost 700 cell lines.

To reveal the underlying mechanism, we employed the weighted gene co-expression network analysis (WGCNA) (Langfelder and Horvath, 2008, 2012). Also, we conducted an association study between modules and clinical aspects and identified critical genes linked to the model. Further use of the Metascape platform (https://metascape.org/gp/index.html) was made to implement the Gene Ontology (GO) analysis of the important genes. GSEA (GSEA v4.2.2 software, http://software.broadinstitute.org/gsea/login.jsp) was also conducted to uncover the differences in biological functions with the Kyoto Encyclopedia of Genes and Genomes (KEGG). Statistical significance was defined as p < 0.05 and q < 0.25.

We employed five techniques to measure the immune microenvironment of subgroups, including EPIC, MCPCounter (Becht et al., 2016), quanTIseq (Plattner et al., 2020), the single-sample gene set enrichment analysis (ssGSEA) and xCell (Aran, 2020), to remove variances across different algorithms. The correlation analysis of the GRGs with immune-related genes and functional status from Thorsson et al. was performed. (Thorsson et al., 2018). mRNAsi is a tool for evaluating how closely cancer cells resemble stem cells (Malta et al., 2018). The sensitivity to immune checkpoint inhibitor (ICI) was evaluated using the immunophenotype score (IPS) (Charoentong et al., 2017) and the tumor immune dysfunction and exclusion (TIDE) (Jiang et al., 2018).

A nomogram combing age, ISS stage and glycolytic risk score was constructed. For the nomogram’s self-verification, the calibration curve was conducted. Based on the merged scores and other clinical parameters, time-dependent receiver operating characteristic (time-ROC) curves were computed for 1-, 3-, and 5-year survival. With the “ggDCA” package, survival net benefits of each clinical feature and the risk score were estimated with decision curve analysis (DCA). DCA curves focus on evaluating the clinical benefit of the model, predicting which threshold probability range treatment for patients in the model will achieve a higher clinical benefit, meeting the practical needs of clinical decision-making. For the prediction model, the net benefit (NB) is a composite indicator that incorporates both true positives and false positives. In the DCA curve, there are two reference lines, one with no intervention for anyone, called treat_none, and the other with intervention for everyone, called treat_all. The model has real value only if its NB is higher than both treat_all and treat_none at some threshold probability (Vickers and Elkin, 2006).

MM cell lines including LP-1, RPMI8226, NCI-H929, MM1.R, I9.2, and U266 cells were collected from Fenghui Biotechnology Co., Ltd. (Hunan, China). Cells were grown in an incubator set at 37 °C with a humid environment that contained 5% CO2. In addition to 10% fetal bovine serum, the cells were cultured in RPMI-1640 medium (Gibco, Shanghai, China) with 100 IU/mL penicillin and 100 mg/mL streptomycin.

35 MM patients were subjected to our experiment, of whom four were relapsed MM (RMM). As controls, normal bone marrow samples were taken from 21 healthy volunteers for PCR on cell lines and patient samples. Table 2 displayed the clinical traits of the patients. The study was authorized by the First Affiliated Hospital of Wenzhou Medical University’s ethics committee, and all operations were conducted according to the informed consent and Helsinki Declaration.

Using the Righton DNA&RNA Blood and Tissue Kit (Righton Bio, Shanghai, China), total RNA was isolated from total bone marrow-derived mononuclear cells as directed by the manufacturer. The cDNA synthesis kit (Vazyme, Nanjing, China) was used to execute reverse transcription in order to produce cDNAs. With the help of the Taq Pro universal SYBR qPCR Master Mix (Vazyme, Nanjing, China), qRT-PCR was used to find the levels of glycolysis-related gene expression, with β-ACTIN acting as an internal reference. The primers used was provided in Table 3.

The statistical analysis was performed by SPSS version 24.0 (SPSS Inc., Chicago, IL, United States), R software version 4.1.1 (R Foundation for Statistical Computing, Vienna, Austria), and GraphPad Prism version 9.0.0 (GraphPad-Software Inc., San Diego, CA, United States). The Student’s t test was used to compare two groups for quantitative variables with a normal distribution. The Chi-square test was for the categorical variables. The Mann-Whitney U test was employed for the abnormal distributional variables. Multiple groups were compared using one-way analysis of variance (ANOVA) and the Kruskal-Wallis test. LSD test is used for multiple comparisons after ANOVA. Pearson’s correlation test was for the correlation evaluation between variables with normal distributions, and the Spearman’s correlation test was for variables with aberrant distributions. Statistics were judged significant at p < 0.05.

The glycolytic model was developed using the GSE136337. The GSE24080 and GSE4204 were used for model validation. Survival data were available for 1,514 subjects in the three datasets (GSE136337, n = 424; GSE24080, n = 556; GSE4204, n = 534). Clinicopathological features in sufficient samples allowed for subsequent Cox regression analysis (GSE136337, n = 415; GSE24080, n = 556) (Table 1). Moreover, the GSE6477 (n = 162) was used to reveal the evolution of GRGs in normal individuals, MGUS, SMM, and active MM groups.

In the GSE136337 (n = 424), 22 genes related to glycolysis were associated with prognosis (p < 0.001) (Figure 1A). Then, the glycolytic model was constructed using 12 genes (λ = 0.09) (Figure 1B). The formula was as follows: glycolytic risk score = (0.0803 × LDHB) + (0.0365 × SOD1) + (0.2087 × MDH2) - (0.1617 × PAM) + (0.1063 × MPC2) + (0.0182 × PRPS1) - (0.0262 × GMPPB) + (0.0217 × CYB5A) + (0.0593 × POLR3K) + (0.0721 × PAXIP1) + (0.1125 × AURKA) + (0.0507 × NSDHL). The median risk score was used to divide the subjects into high-risk and low-risk subgroups.

Kaplan–Meier curves revealed survival differences among the subtype in the training (Figure 1C) and the two validation datasets (Supplementary Figures S1A, B). The high-scoring population had worse survival (HR = 3.47, 95% CI = 2.51–4.80, p < 0.0001; HR = 2.46, 95% CI = 1.78–3.41, p < 0.0001; HR = 2.09, 95% CI = 1.35–3.22, p < 0.001). The heterogeneity of survival outcomes was further validated in Figure 1D; Supplementary Figures S1C, D. Heat maps were also plotted to compare the expression of the 12 glycolytic genes. In general, the high-risk groups had lower levels of GMPPB and PAM expression, whereas other genes exhibit the opposite pattern (Figure 1D; Supplementary Figures S1C, D). AUCs of the 1-, 3-, and 5-year survival were 0.671, 0.715, and 0.731 in the GSE136337 (n = 424) to assess the sensitivity and specificity (Figure 1E). The results of the validation datasets were exhibited in Supplementary Figure S1E (GSE24080, n = 556) and 1F (GSE4204, n = 534).

Depending on the TCGA data, we explored the mutation landscape of 293 genes in the candidate glycolytic gene set in MM. Figure 2A showed the top 20 genes with the highest mutation rate. Among them, missense mutations accounted for the major part and the single nucleotide polymorphisms (SNPs) occurred more frequently than insertions or deletions. Furthermore, we used cBioPortal to reveal the single nucleotide variants (SNVs) and copy number variations (CNVs) of the 12 genes (Figure 2B). The alterations of LDHB, SOD1, MDH2, MPC2, PRPS1, AURKA, PAXIP1 and NSDHL in cancer cell lines were 6%, 2.1%, 5%, 3%, 1.4%, 10%, 10%, and 4% separately. Amplification was the most common change. In contrast, both PAM and GMPPB were altered at a frequency of 5%, and their most common change was deletion (Figure 2B). These results were in general agreement with our data.

Moreover, the co-expression matrix and the PPI network showed a close relationship between the 12 genes (Figures 2C, D). In the CCLE database, LDHB, SOD1, MDH2, MPC2, PRPS1, CYB5A, POLR3K, and AURKA were over-expressed at the cellular level in MM, while PAM and GMPPB were under-expressed (Figure 2E), corresponding to the model equation above.

The glycolytic risk score was shown to be independently associated with survival (GSE136337, n = 415, HR = 3.191 (2.295, 4.436), p < 0.001; GSE24080, n = 556, HR = 1.662 (1.214, 2.275), p = 0.020) (Table 4). Subsequently, we analyzed the relationship between the 12-gene signature and clinicopathological traits in the GSE136337 (n = 415). The risk scores of the higher LDH and the high-risk cytogenetics groups were higher than those counterparts (p < 0.001) (Figure 3A). We observed the gradually increasing tendency of the risk score with the increasing ISS or R-ISS stage (p < 0.05) (Figure 3A). The distribution of the patient’s characteristics and survival outcomes were illustrated in Figure 3B.

According to the estimated IC50 values of the chemotherapeutic agents, the high-risk cohort showed the resistance to bortezomib, doxorubicin, and cytarabine (p < 0.001), but more sensitive to methotrexate and vorinostat (p < 0.0001), compared to the low-risk group (Figure 3C). Additionally, glycolysis-related gene expression changes were also present in the evolving disease spectrum of myeloma (GSE6477, n = 162) (Figure 3D). The expression levels of LDHB and SOD1 exhibited an increasing trend as the disease progressed. In addition, MDH2, MPC2, PAXIP1, and NSDHL were upregulated in myeloma patients. GMPPB, however, revealed a pattern of decline.

We conducted an investigation of biological functions among subgroups to uncover the mechanism underlying the glycolytic signature. WGCNA was implemented to build a weighted gene co-expression network. To ensure proper connectivity and sustain the network nearly scale-free, the soft threshold was set at four (Figure 4A). 19 gene modules were produced after grouping related modules (Figure 4B). The red module and the risk score had the strongest association (r = 0.58, p < 0.001) (Figure 4C). The red module’s 221 genes were analyzed for GO enrichment using the online database Metascape. Important processes in tumors, including the mitotic cell cycle, cell cycle checkpoints, DNA metabolism, and the retinoblastoma gene in cancer, were the main sites of gene module concentration (Figure 4D).

In GSEA, the high-scoring population was home to the considerably enriched pathways, which were mostly linked to glycolysis and included pyrimidine metabolism, one carbon pool via folate, and oxidative phosphorylation. Consistent with previous results, the pathways of cell cycle, DNA repair, and nucleotide excision repair can also be found (Figure 4E).

The immune microenvironment is incredibly important to the growth of MM. Based on several different algorithms, the heat map demonstrated the heterogeneity in the immune landscape among the subgroups (Figure 5A). Immune cell infiltration levels were higher in the low-scoring cohort (p < 0.05), such as active B cells, effector memory CD8⁺ T cells, immature B cells, NK cells, plasmacytoid dendritic cells, and Th2, while stromal cell abundance was higher in the high-scoring population (p < 0.01) (Figures 5B, C). Moreover, we performed a correlation analysis of the GRGs with immune-related genes and functional status. The expression of most prognostic genes was negatively correlated with the immune-related genes (Figure 5D). For example, LDHB, SOD1, MDH2, POLR3K, and PAXIP1 were negatively correlated with immune-related genes in terms of antigen presentation and cell adhesion, while GMPPB was positively correlated (Figure 5D).

Meanwhile, higher risk scores were observed to be strongly associated with the stemness index (Figure 5E). IPS characterizes the immune phenotype of the cells from four perspectives (“antigen-presenting, AP; effector cells, EC; suppressor cells, SC; checkpoints, CP”) (Charoentong et al., 2017). And a total score (z-score, AZ) is finally generated after normalization. A higher z-score corresponds to higher ICI responsiveness. Consistent to our previous results, GMPPB was positively correlated with antigen presentations (p < 0.01). Compared to GMPPB, MDH2, and POLR3K were linked to higher concentrations of immunosuppressive cells (p < 0.01). Furthermore, high expression of GMPPB was often accompanied by greater sensitivity to ICI (p < 0.01), while LDHB and PAXIP1 showed the opposite pattern (p < 0.05) (Figure 5F). The quantity and functional status of T cells is important mediators in immune-targeted therapy. The T-cell exclusion was more prone to occur in the TME of the high-scoring cohort (p < 0.0001), foreshadowing their possible resistance to ICI (Figure 5G). In conclusion, glycolysis-related markers may be a useful tool for determining the immune status of MM and for gauging the effectiveness of immunotherapy.

For improving the survival prediction efficiency of the model, a combined nomogram model was constructed with age, ISS stage, and glycolytic risk score (GSE136337, n = 415) (Figure 6A). Considering the clinical information of the validation cohort, we did not include cytogenetics. The nomogram showed a strong predictive value, as shown by a C-index of 0.728. An ideal concordance between the forecast and observation was visible on the survival probability calibrations plot (Figure 6B). The 1-, 3- and 5-year AUCs were 0.75, 0.75, and 0.78, respectively, higher than the AUCs of the ISS and R-ISS (Figure 6C). The same was in the validation set GSE24080 (n = 556) (Figure 6D).

Moreover, graphical representations provided by DCA showed that the glycolytic risk score outperformed other metrics, such as the R-ISS, in terms of the net benefit to survival (GSE136337, n = 415) (1-year, threshold:0.125–0.250; 3-year, threshold:0.125–0.563; 5-year, threshold:0.125–0.750) (Figure 6E).

MM cell lines and patient samples were subjected to PCR to further validate our predictive model. In these six cell lines (LP-1, RPMI8226, NCI-H929, MM1.R, I9.2, and U266), AURKA, CYB5A, LDHB, MDH2, MPC2, NSDHL, PAXIP1, POLR3K, PRPS1, and SOD1 exhibited primarily upregulation (p < 0.01). GMPPB and PAM, on the other hand, underwent downregulation (p < 0.05) (Figure 7). Correspondingly, the experimental results from the patients were consistent with the cell lines (p < 0.0001) (Figure 8A). We then further investigated the changes in gene expression profiles related to glycolysis among the RMM patients. Figure 8B showed that the expression of GMPPB and PAM was significantly lower than that of controls (p < 0.001), while the opposite was true for other genes (p < 0.01).