Skip to main content

Integrative analysis of the prognostic value and immune microenvironment of mitophagy-related signature for multiple myeloma



Multiple myeloma (MM) is a fatal malignant tumor in hematology. Mitophagy plays vital roles in the pathogenesis and drug sensitivity of MM.


We acquired transcriptomic expression data and clinical index of MM patients from NCI public database, and 36 genes involved in mitophagy from the gene set enrichment analysis (GSEA) database. Least absolute shrinkage and selection operator (LASSO) Cox regression analysis was conducted to construct a risk score prognostic model. Kaplan–Meier survival analysis and receiver operation characteristic curves (ROC) were conducted to identify the efficiency of prognosis and diagnosis. ESTIMATE algorithm and immune-related single-sample gene set enrichment analysis (ssGSEA) was performed to uncover the level of immune infiltration. QRT-PCR was performed to verify gene expression in clinical samples of MM patients. The sensitivity to chemotherapy drugs was evaluated upon the database of the genomics of drug sensitivity in cancer (GDSC).


Fifty mitophagy-related genes were differently expressed in two independent cohorts. Ten out of these genes were identified to be related to MM overall survival (OS) rate. A prognostic risk signature model was built upon on these genes: VDAC1, PINK1, VPS13C, ATG13, and HUWE1, which predicted the survival of MM accurately and stably both in training and validation cohorts. MM patients suffered more adverse prognosis showed more higher risk core. In addition, the risk score was considered as an independent prognostic element for OS of MM patients by multivariate cox regression analysis. Functional pathway enrichment analysis of differentially expressed genes (DEGs) based on risk score showed terms of cell cycle, immune response, mTOR pathway, and MYC targets were obviously enriched. Furthermore, MM patients with higher risk score were observed lower immune scores and lower immune infiltration levels. The results of qRT-PCR verified VDAC1, PINK1, and HUWE1 were dysregulated in new diagnosed MM patients. Finally, further analysis indicated MM patients showed more susceptive to bortezomib, lenalidomide and rapamycin in high-risk group.


Our research provided a neoteric prognostic model of MM based on mitophagy genes. The immune infiltration level based on risk score paved a better understanding of the participation of mitophagy in MM.


The study reported a mitophagy-related genes signature in multiple myeloma.

The mitophagy-related genes signature had ideal prognostic independence in multiple myeloma.

A nomogram to predict the overall survival of multiple myeloma was built by combining the five-gene signature, LDH and ISS stage.

Immune infiltration was related to the mitophagy-related risk signature.

Peer Review reports


Multiple myeloma (MM) is a fatal hematologic cancer featured with abnormal propagation of monoclonal plasma cells in bone marrow [1]. Over the past decades, many effective treatments, such as bortezomib (injectable proteasome inhibitor), lenalidomide (oral immunomodulatory drug), chimeric antigen receptor-engineered T cells, and autologous stem cell transplantation [2,3,4,5], improve the outcome of MM patients. However, the medical need of MM remains unmet as a result of the significant heterogeneity [6]. Various factors are involved in MM progress, mainly including genetic abnormalities [7], changes in bone marrow microenvironment (BM-ME) [8], and epigenetic alterations [9]. Studies confirmed that immune cells in BM-ME and the dysregulation of genes involved in immune checkpoint are related to the immune infiltration level of BM-ME [10, 11]. The protection from BM-ME and high genetic instability guard MM cells against chemotherapies, or receptor-targeting drugs, which virtually leads to resistance and relapse [12].

Mitochondria is not only a crucial house of producing energy via oxidative phosphorylation, but also a center of producing cellular metabolites [13]. Mitophagy is a process of lysosome-dependent mitochondrial autophagy, which protects cells against proapoptotic proteins, poisonous reactive oxygen species (ROS), and the unavailing hydrolytic action of adenosine triphosphate (ATP), induced by the depolarization of mitochondrial membrane or the mitochondrial DNA (mtDNA) changes [14, 15]. Accumulating evidences have recently confirmed that mitophagy is a double-edged sword in cancer development. On one hand, the decrease of mitophagy promotes the cancer progression [16]. On the other hand, the increased mitophagy facilitates cancer cells proliferation and progression by defending cancer cells from apoptosis [17]. Nonetheless, roles of mitochondrial dysfunction in the immune microenvironment and prediction of outcome in MM remain indistinct.

In our study, we performed the least absolute shrinkage and selection operator (LASSO) Cox regression analysis to build a five-mitophagy-related-gene prognostic risk signature model. The risk model revealed a great predicted value of actual survival probabilities. In addition, enrichment analysis identified the alteration of immune checkpoint and immune microenvironment based on the risk score. Finally, drug sensitivity analysis predicted latent drugs for treating MM.

Materials and methods

Data acquisition

MM transcriptomic expression data and clinical features were acquired from the NCI Gene Expression Omnibus database (GEO). The raw data was normalized and log2 transformed. All detailed information of GEO dataset were shown in Supplementary Table 1. The overall design of this study was shown in Supplementary Fig. 1. Among them, GSE6477 and GSE13591 were applied to appraisal differentially expressed genes (DEGs). GSE9782 was used as the training set to establish the prognostic risk score model, and GSE24080 and GSE4204 were utilized as validation cohorts. GSE24080 was applied to perform univariate analysis and multivariate analysis for overall survival (OS) rate and construct the nomogram model. GSE6477 and GSE47552 were applied to evaluate the performance of risk score for MM diagnosis.

Establishment of the prognostic risk model

First, univariate cox regression analysis was applied to obtain the OS-related genes with p < 0.05, and ten genes (SLC25A4, VDAC1, RNF41, SLC25A5, PINK1, SQSTM1, VPS13C, ATG13, HUWE1, and OPTN) were significantly correlated with MM OS time. Next, we performed LASSO Cox regression to optimize the prognostic model by further compressing the genes and constructing the prognostic model by “glmnet” package (version 4.1-1), and five genes finally came into the risk score formula (VDAC1, PINK1, VPS13C, ATG13, and HUWE1). Moreover, MM patients were divided into two groups upon the optimal cutoff of the risk score with “Survminer” package (version 0.4.9). Receiver operating characteristic (ROC) curve was applied to estimate the prognostic value of risk score model in MM patients.

Construction of the nomogram

To assess the prediction value of risk score in MM, we performed univariate and multivariate cox regression analysis. Variables with p < 0.10 was supposed to the multivariate cox regression analysis, and p < 0.05 was considered as remarkable independent prognostic factors. Then, we used “rms” package (version 6.2-0) to construct the nomogram independent prognostic model. Finally, we assessed the predictive value of nomogram by ROC curve and calibration curve.

Functional pathway enrichment analysis

We used the online website Metascape ( to analyze gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway [18,19,20] of DEGs based on risk score group. We executed Gene set enrichment analysis (GSEA) by GSEA software (version 4.1.0). The absolute value of NES value > 1, p value < 0.05, and false discovery rate (FDR) value < 0.25 showed significance.

Estimation of immune infiltration

We computed the immune score, stromal score, tumor purity, and ESTIMATE score for each sample in GSE24080 by ESTIMATE algorithm [21]. According to the previous study [22], single sample GSEA (ssGSEA) was conducted to evaluate the abundance of infiltrating cells by “GSVA” package (version 1.38.2).

Prediction of drug sensitivity

To achieve precise treatment upon mitophagy-related signature and identified potential drugs for MM, we utilized the Genomics of drugs sensitivity in cancer (GDSC, to predict the chemotherapeutic response [23]. R package “pRROphetic” was implemented to evaluate the half-maximal inhibitory concentration (IC50) by ridge regression [24]. The forecast precision was determined through 10-fold cross-validation using the GDSC training set. Additional methods are described in supplementary methods.

Statistical analysis

At least three dependent experiments were performed in qRT-PCR and values were presented as mean ± SD. We used unpaired Student’s t-test and Mann-Whitney U test to judge the otherness in two groups. More than two groups, we used One-way ANOVA (for parametric data) and Kruskal-Wallis (for non-parametric data) test to contrast significance. ROC curve was performed to assess the diagnosis value of risk score and the prognostic value of nomogram model in MM. The Kaplan-Meier method with a two-sided log-rank test was performed to evaluate the OS of MM patients. SPSS 21 software (SPSS, Chicago, USA) and GraphPad Prism 8 were used for statistical analysis. P value lower than 0.05 was defined markedly different.


Identification of 15 mitophagy-related genes differentially expressed in MM

We acquired a total of 36 mitophagy-related genes from GSEA GOBP database (supplementary Table 2). The datasets of GSE13591 (including 5 normal plasma cells (NPC) and 133 MM patients) and GSE6477 (including 15 NPC and 73 MM patients) were used to analyze the genes involved in mitophagy which were differentially expressed in MM. We found 19 genes and 21 genes were differentially expressed in MM patients compared to human healthy donors in GSE6477 and GSE13591, respectively (Fig. 1A and B). We observed that SLC25A4, PHB2, CERS1, VPS13C, HUWE1, VDAC1, and SLC25A5 were significantly upregulated, while OGT, ATG13, PINK1, and OPTN were significantly downregulated in MM patients (supplementary Fig. 2A and 2B). Then, venn plot was used to identify 15 mitophagy-related genes overlapped in the two datasets (Fig. 1C). To ulteriorly probe the relationship among 15 mitophagy-related genes, we performed PPI network analysis (Fig. 1D). We found PINK1, SQSTM1 and VDAC1 were hub genes. Besides, we conducted correlation analysis to ascertain the relationship among these genes (Fig. 1E). The correlation between VDAC1 and SLC25A5 was significantly positive (r = 0.678), whereas the correlation between SQSTM1 and PHB2 was significantly negative (r = -0.416). Furthermore, we used cBioPortal, an online database, and found no frequent mutation of 15 mitophagy-related genes in 211 MM samples (Fig. 1F) [25].

Fig. 1
figure 1

Identification of 15 mitophagy-related genes differentially expressed in MM. (A). 19 mitophagy-related genes were differentially expressed in GSE6477. (B). 21 mitophagy-related genes were differentially expressed in GSE13591. (C). The veen plot of overlapped mitophagy-related genes. (D). The PPI network of the 15 mitophagy-related genes. (E). The correlation of 15 mitophagy-related genes in GSE6477. (F). The mutant frequency of 15 mitophagy-related genes in 211 multiple myeloma samples, including chromosome translocation and hyper-diploid

Establishment of prognostic risk model upon on mitophagy-related genes

We chose GSE9782, which includes 264 MM patients with survival time, to screen the OS-related mitophagy-related genes. We found the low expression of ATG13 (p < 0.0001), HUWE1 (p < 0.0001), OPTN (p < 0.0001), PINK1 (p < 0.0001), SQSTM1 (p = 0.001) and VPS13C (p = 0.029) was notably related to adverse OS rate, while the high expression of SLC25A5 (p = 0.004), VDAC1 (p < 0.0001), RNF41 (p = 0.039), and SLC25A4 (p = 0.009) was notably related to poor OS (Fig. 2A). Then, we executed LASSO cox regression analysis to establish a prognostic model in the training cohort GSE9782 (Fig. 2B C). Five mitophagy-related genes, including VDAC1, PINK1, VPS13C, ATG13 and HUWE1, were involved in the risk model. We evaluated the total risk score of every patient with formula as follows: (5.476831e-04*the expression of VDAC1) + (-4.209033e-03* the expression of PINK1) + (-1.629760e-03* the expression of VPS13C) + (-2.450343e-03* the expression of ATG13) + (-7.849962e-05* the expression of HUWE1). On the bias of the cutoff risk score, 264 MM patients were separated into low- and high-risk group. Kaplan-Meier survival analysis found that MM patients had poorer outcome observably in high-risk group than that in low-risk score (p < 0.0001) (Fig. 2D). ROC curve suggested that the AUC values were 0.6807 (95% confidence interval (CI): 0.6147 to 0.7467, p < 0.0001) and 0.7205 (95%CI: 0.6471 to 0.7940, p < 0.0001) for survival rate of 1- year and 2-year, respectively (Fig. 2E). Next, the total risk score was calculated in other two independent cohorts for validation. MM patients had poorer OS with high-risk score than those with low-risk score in the GSE24080 (p = 0.0002) and GSE4204 (p < 0.0001) (supplementary Fig. 3A). The AUC values were 0.6236 (95%CI: 0.5064 to 0.7408, p = 0.0484), 0.6432 (95%CI: 0.5589 to 0.7275, p = 0.0028), and 0.6085 (95%CI: 0.5355 to 0.6015, p = 0.0038) for survival rate of 1-, 2-, and 3-year in GSE24080, respectively (supplementary Fig. 3B). Parallel values were 0.5739 (95%CI: 0.5180 to 0.6297, p < 0.0001), 0.5894 (95%CI: 0.5416 to 0.6373, p = 0.0003), and 0.6195 (95%CI: 0.5680 to 0.6710, p < 0.0001) for survival rate of 1-, 2- and 3-year in GSE4204, respectively (supplementary Fig. 3B).

Fig. 2
figure 2

Construction of prognostic risk model based on mitophagy-related genes. (A). Forest plot of the univariate cox regression in GSE9782. (B). The LASSO Cox analysis identified mitophagy-related genes in the training cohort. (C). Partial likelihood deviance of different numbers of variables. One-thousand- fold cross-validation was applied for tuning penalty parameter selection. (D). Kaplan-Meier survival analysis in training cohort based on the risk score. (E). ROC curve for the prognostic risk signature. LASSO: least absolute shrinkage and selection operator; ROC: receiver operating characteristic

The mitophagy-related risk score was an independent prognostic element in MM

We applied GSE24080 database to performed univariate and multivariate cox regression analysis. We observed lactate dehydrogenase (LDH) (hazard ratio (HR), 2.134; 95% CI, 1.39–3.277; p = 0.001), international staging system (ISS) stage (HR, 1.83; 95% CI, 1.414–2.369; p < 0.0001), and risk score (HR, 1.601; 95% CI, 1.059–2.422; p = 0.026) were the independent prognostic elements for MM survival (Table 1). In order to visualize the prognosis prediction, we further constructed a nomogram model with the above independent prognostic features (Fig. 3A). Calibration curves are typically used to assess the agreement between predicted probabilities and observed event rates or frequencies in the real world, with the 45° line representing the best prediction scenario. Calibration ROC and curves were applied to estimate the precision and efficacy of the prognostic risk model. Decision curve analysis (DCA) is usually used to evaluate the clinical utility of different predictive models [26]. We observed that the precision of the total points in nomogram for 1- year survival was 0.8240 (95%CI: 0.7570 to 0.8911, P < 0.0001), 2- year survival 0.7729 (95%CI: 0.6912 to 0.8545, P < 0.0001), and 3- year survival 0.6653 (95%CI: 0.5933 to 0.7372, P < 0.0001), respectively (Fig. 3B). The calibration curves and DCA revealed a great predicted value of actual survival probabilities for 1-year, 3-year, and 5-year, respectively base on the nomogram model (Fig. 3C and D).

Table 1 Univariate analysis and Multivariate analysis for OS in GSE24080
Fig. 3
figure 3

The mitophagy-related risk score was an independent prognostic factor in MM. (A). Nomogram predicting 1-, 3-, and 5-year survival for MM patients based on mitophagy-related genes risk score. To use this nomogram, the specific point for each variable of the patient lies on each variable axis. Draw a vertical line upward to determine the point at which each variable accepts; the sum of these points is located on the Total Points axis, and draw a vertical line down to the survival axis to determine the probability of 1-, 3- and 5- year overall survival. LDH: lactate dehydrogenase; ISS: international staging system. (B). ROC curve for the nomogram prognosis system. (C). DCA curve for the nomogram prognosis system. The x-axis represents the threshold probability, while the y-axis stands for the net benefits. (D). Calibration plot of the nomogram for 1-year, 3-year, and 5-year OS. The risk score in normal donors and MM patients with different stages in GSE6477 (E) and GSE47552 (F). (G). ROC curve for MM with different myeloma stages. DCA: decision curve analysis * p < 0.05, ** p < 0.01, *** p < 0.001

The development of MM includes 4 steps: monoclonal gammopathy of unknown significance (MGUS), smoldering multiple myeloma (SMM), multiple myeloma, and plasma cell leukemia (PCL) or refractory and/or relapse MM (RRMM). The results found that MM patients and RRMM patients had more dramatically higher risk score than that in normal people (Fig. 3E F). Besides, ROC curve suggested that the value of risk score was able to distinguish MM from healthy donor (AUC = 0.8073, 95%CI: 0.6962 to 0.9184, P = 0.0002), MGUS (AUC = 0.8599, 95%CI: 0.7622 to 0.9516, P < 0.0001), and SMM (AUC = 0.7135, 95%CI: 0.5939 to 0.8330, P = 0.0018) (Fig. 3G).

The mitophagy-related risk signature was related to clinical characteristics of MM

Next, we used the database of GSE24080 to estimate the relation between the risk score and clinical characteristics (Table 2). The results showed that the value of risk score was related to LDH (p = 0.03), albumin (ALB) (p = 0.025), hemoglobin (HGB) (p = 0.02), beta-2 microglobulin (B2M) (p = 0.037), ISS stage (p = 0.022), the percentage of plasma cells in bone marrow biopsy (BMPC) (p = 0.018), and cytogenetic abnormalities (p < 0.0001). Cytogenetic abnormalities has been regarded as a poor prognostic factor for MM [27]. Therefore, we choose GSE136337 as the validation dataset, which consists with several cytogenetic information. We observed that MM patients with cytogenetic abnormalities, including del13q (p < 0.0001), del1p32 (p < 0.0001), del1p (p < 0.0001), del1q (p < 0.0001), myc_8q24 (p = 0.001), x1qplus (p < 0.0001), and hyperdiploid (p < 0.0001) had higher risk score (Table 3). But we don’t find the relationship between risk score and t (4;14), t (14;20), and t (14;16), because the number limitation of MM patients with these cytogenetic abnormalities. Further, we analyzed the gene levels of five mitophagy-related genes in the risk score formula. We observed MM patients with cytogenetic abnormalities had higher VDAC1 expression and lower PINK1 expression (supplementary Table 3). These results indicated mitophagy-related risk signature was associated with clinical characteristics of MM, especially cytogenetic abnormalities.

Table 2 The correlation of prognostic risk signature and clinical characteristics
Table 3 Impact of risk-score on the most common cytogenetic abnormalities of multiple myeloma (GSE136337)

To further assess the applicability of the risk signature model for different MM conditions, we performed Kaplan–Meier survival analysis. The risk score showed stabilized predictive value in term of ISS stage (ISS I/II, p = 0.0052, ISS III, p = 0.0079, Fig. 4A and B), gender (male, p = 0.0075; female, p = 0.0016, Fig. 4C and D), the condition of cytogenetic (no cytogenetic abnormalities, p = 0.022, with cytogenetic abnormalities, p = 0.027, Fig. 4E F). Among these, MM patients had the shorter OS in higher-risk score group. However, this phenomenon was found in young people (p = 0.00096, Fig. 4G), but not in older people (p = 0.06, Fig. 4H). In addition, higher risk score was related to more adverse OS in MM patients received various treatment, including bortezomib (p = 0.017, Fig. 4I), thalidomide (p = 0.0013, Fig. 4J), dexamethasone (DEX, p < 0.0001, Fig. 4K), and PS341 (proteasome inhibitor, p < 0.0001, Fig. 4L). Taken together, the prognostic model satisfied the personalized medical needs of MM.

Fig. 4
figure 4

Survival analysis in different subgroup. Kaplan–Meier survival analysis of MM patients grouped by ISS stage I/II/III (A, B), male and female (C, D), no cytogenetic abnormalities and with cytogenetic abnormalities (E, F), age < 65 and ≥ 65 (G, H), bortezomib treatment and thalidomide treatment (I, J) in GSE24080. Kaplan–Meier survival analysis of MM patients stratified by dexamethasone treatment and PS341 treatment (K, L) in GSE9782.

Gene expression and functional pathway enrichment analysis of mitophagy-related risk signature in MM

To further explore the underlying mechanisms of mitophagy-related genes in MM, we used limma package to get DEGs in GSE24080 based on the risk score. A total of 83 increased genes and 80 decreased genes were discerned (supplementary Fig. 4A and supplementary Table 4). Enrichment analysis of DEGs uncovered that upregulated genes were markedly enriched in mitotic cell cycle, ATPase activity, and ATP binding in GO term (supplementary Fig. 4B), and p53 signaling pathway, cell cycle, Hippo signaling pathway, FoxO signaling pathway, and TGF-beta signaling pathway in KEGG (supplementary Fig. 4C). Downregulated genes were overtly enriched in term of immune response, cytokine production, T cell receptor binding, and Toll-like receptor binding in GO term (Fig. 5A), and Hematopoietic cell lineage, cell adhesion molecules, and Th17 cell differentiation in KEGG (Fig. 5B). Moreover, GSEA revealed that high-risk group significantly enriched tumorigenesis and mitochondrial terms, including cell cycle (NES=-2.23, p < 0.0001, FDR = 0.012), p53 signaling pathway (NES=-1.62, p = 0.008, FDR = 0.076) in KEGG gene set, MYC targets (NES=-2.48, p < 0.0001, FDR = 0.0005), G2M checkpoint (NES=-2.25, p < 0.0001, FDR = 0.001), E2F target (NES=-2.17, p < 0.0001, FDR = 0.003), and mTORC1 signaling (NES=-1.92, p = 0.014, FDR = 0.013) in HALLMARK gene set, and ribosome biogenesis (NES=-2.51, p < 0.0001, FDR = 0.004), and mitochondrial gene expression (NES=-2.30, p = 0.002, FDR = 0.004) in GOBO gene set (Fig. 5C and supplementary Fig. 4D).

Fig. 5
figure 5

Enrichment analysis of mitophagy-related risk signature. Gene ontology and KEGG pathway enrichment analysis of DEGs in GSE24080 dataset grouped by risk score. (A). The significantly enriched gene ontology biological process (GOBP) of downregulated DEGs. (B). The significantly enriched KEGG pathway of downregulated DEGs. (C). GSEA analysis results of KEGG gene set, HALLMARK gene set, and GOBP gene set. NES, normalized enrichment score; FDR, false discovery rate

Immune infiltration was related to the mitophagy-related risk signature

Due to the significantly enrichment of immune response in downregulated DEGs in MM patients, ESTIMATE was performed to assess the different immune infiltration levels of MM patients in different risk group. First, Kaplan-Meier analysis revealed that low stromal score, immune score and ESTIMATE score were associated with adverse OS, while high tumor purity was associated with poor OS (Fig. 6A). In addition, MM patients had lower levels of stromal score, immune score and ESTIMATE score in the high-risk group, and higher levels of tumor purity in both training cohort and validation cohort (Fig. 6B C). Furthermore, we conducted ssGSEA to probe the disparate immune cell subsets based on the risk score. The results found 19 of 28 immune cells were notably dysregulated in high-risk group (Fig. 6D). Activated B cell, central memory CD4 T cell, activated dendritic cell, central memory CD8 T cell, CD56 bright killer cell, effector memory CD8 T cell, immature B cell, nature killer cell, nature killer T cell, and type 17 T helper cell, participated in anti-tumor immune response, were all decreased in high-risk group. Besides, the levels of immature dendritic cell, regulatory T cell, MDSC, macrophage and mast cell, belonging to immunosuppressive cells, were also decreased in MM patients with high-risk score.

Fig. 6
figure 6

Immune infiltration was related to the mitophagy-related risk signature. (A). Kaplan-Meier survival analysis upon stromal score, immune score, ESTIMATE score, and tumor purity. (B). The distribution of stromal score, immune score, ESTIMATE score, and tumor purity upon risk score in GSE24080. (C). The distribution of stromal score, immune score, ESTIMATE score, and tumor purity upon risk score in GSE9782. (D). The heatmap of the comparison in 28 immune-related gene sets upon risk score in GSE24080. * p < 0.05, ** p < 0.01, *** p < 0.001

Moreover, we explored the relation between risk score and immune checkpoint, which was regarded as potential therapy targets for MM [28]. Intriguingly, the results showed that risk score was connected with immune checkpoint negatively, including IDO1, CD276, CD86, PD-L1, and PD-L2, while positively related to CD279 (Fig. 7A). Taken together, the prognostic risk score was connected with tumor immune infiltration levels and genes involved in immune checkpoint. Therefore, our research verified the immune checkpoint inhibitors, including PD-L1, IDO1, CD276, CD86 had potential clinical values for MM patients.

Fig. 7
figure 7

Validations for mitophagy-related genes in clinical MM samples. (A). Correlation analysis between risk score and immune checkpoint for tumor-targeted treatment. The mRNA expression of VDAC1 (B), PINK1 (C), HUWE1 (D), VPS13C (E), and ATG13 (F) by qRT-PCR in primary bone marrow mononuclear cells from NDMM and normal donors. Data represent the mean ± SD. * p < 0.05, ** p < 0.01, *** p < 0.001

Validation for mitophagy-related genes in MM samples

To validate these results of bioinformatic analysis, we collected the mononuclear cells of NDMM patients and healthy donors in bone marrow. Then, the qRT-PCR was executed. We observed the expression levels of VDAC1 were markedly upregulated in MM patients than that in healthy donors (P = 0.0056, Fig. 7B), while the expression levels of PINK1 and HUWE1 were downregulated in MM patients (P = 0.001, P = 0.019, Fig. 7C and D). Although, the expression level of VPS13C was lower and ATG13 was higher in MM patients than these in controls, no evident differences were observed (P = 0.14, P = 0.838, Fig. 7E F), maybe for the reason of insufficient sample size.

Prediction of possible drugs for MM upon risk score signature

To further identified potential drugs for MM, we explored the estimated half-maximal inhibitory concentration (IC50) between the two group on the bias of the Genomics of Drug Sensitivity in Cancer (GDSC). We found decreased IC50 of 16 drugs in MM patients with high-risk score, including bortezomib and lenalidomide, which are widely used in MM therapy (Fig. 8). Besides, the estimated IC50 of ABT.263 (bcl-2 inhibitor), ABT. 888 (PARP inhibitor), AICAR (AMPK activator), ATRA (all-trans retinoic acid), dasatinib (tyrosine kinase inhibitor), AZD8055 (mTOR inhibitor), erlotinib (EGFR inhibitor), etoposide, MG.132 (proteasome inhibitor), parthenolide (HDAC inhibitor), rapamycin (mTORC1 complex inhibitor), and thapsigargin (ATPase inhibitor) were also lower in high-risk score group (Fig. 8), which might provide novel insights into MM treatment. To sum up, the risk score can be considered as an index to choose appreciate therapeutic targets for MM precisely.

Fig. 8
figure 8

Prediction of potential drugs for MM based on risk score signature. The distribution of the IC50 of various drugs upon risk score in GSE24080.


As the second most common hematological malignancy, MM remains incurable because of drug resistance and relapse [29]. There is no excellent survival prediction for MM so far. Mitophagy is a vital cellular progress, which results in degradation of dysfunctional mitochondria. Several proteins have been reported involved in mitophagy, including PINK1, Parkin, OPTN, ATG13, and p62/SQSTM1, and FUNDC1, etc. [30]. Up to now, increasing researches have revealed that mitophagy acts essential part in the development and drug sensitivity of various cancers, especially in MM. It has been reported that Parkin/PARK2 carried mutation in glioma [31], lung cancer [32], and breast cancer [33]. In addition, the expression of thioredoxin was notably increased in MM cells, which resistant to bortezomib. Increased thioredoxin increased resistance to bortezomib in MM via mitophagy inactivation [34]. However, as far as we known, there is no research focus on the potential relationship between mitophagy-related genes and MM prognosis.

Herein, we constructed a five-mitophagy-related risk signature for MM. The prognostic model worked well in both training and validation cohort via ROC curve and Kaplan-Meier survival analysis. Further enrichment analysis displayed that the DEGs between high- and low-risk groups were observably enriched in term of cell cycle, ATP-binding, and tumor-related pathways. Immune infiltration analysis had manifested that MM patients in high-risk score group suffered a notably decreased condition of antitumor immune activity. Finally, drug sensitivity analysis provided potential strategies for treating MM.

In our study, the mitophagy-related risk signature consists of VDAC1, PINK1, VPS13C, ATG13, and HUWE1, which can accurately predict MM prognosis. MM patients in high-risk score group suffered more adverse OS. Our findings are consistent with previous studies. VDAC1(voltage-dependent anion channel 1), a main factor in the outer mitochondrial membrane, was upregulated in multiple cancers, including lung cancer [35], breast tumor [36], cervical tumor [37], and liver cancer [38]. Overexpressed VDAC1 was a prognostic element and was related to immune infiltrates in breast cancer [36]. Besides, VDAC1 involved in bromodomain inhibitor resistance in breast cancer [39]. VDAC1 was upregulated in mesothelioma patients. The expression levels of VDAC1 were positive related to disease stage and high expression of VDAC1 had correlation with shorter survival rates. Loss of VDAC1 could inhibit cell propagation in mesothelioma cancer cells [40]. PINK1 (PTEN-induced kinase 1)/Parkin was a classical axis to regulate mitophagy in mammalian cells [41]. Previous researches uncovered that PINK1 and Parkin decreased in colorectal cancer [42]. Loss of PINK1 restrained mitophagy, facilitated the Warburg effect, and promoted macrophages polarization in gastric cancer [43]. Besides, another study found PINK1 was reduced in MM and regulated the MOB1B-mediated Hippo‐YAP/TAZ pathway leading to MM migration and homing [44]. Matrine, a natural alkaloid, accelerated the apoptosis of liver cancer cells by intercepting the PINK1/Parkin pathways and restraining mitophagy [45]. Consistent with our research, Li et al. identified an autophagy-related signature for MM prognosis, which contained PINK1, EIF2AK2, KIF5B, MYC, NRG2, and VEGFA [46]. These finding revealed the important role of PINK1 in MM. Recent study revealed that VPS13C (vacuolar protein sorting 13 homolog C) was descended and associated with poor prognosis in skin cutaneous melanoma [47]. ATG13 (autophagy related 13) is a member of autophagy initiation complex. ULK1 could active autophagy via phosphorylating ATG13 to inhibit the progress of breast cancer [48]. Zhou et al. reported that estrogen receptor α (ERα) could bind the promoter of ATG13 leading to the increasing of mitophagy. Erα inhibitor oxabicycloheptene sulfonate could reduce the expression of ATG13, leading to restrain the viability of breast cancer cells [49]. Niu et al. reveled that chemotherapy drug licochalcone A could active the upstream of autophagy ULK1/ATG13 complex, inducing hepatocellular carcinoma cells apoptosis by inducing autophagy [50]. Previous studied demonstrated HUWE1 played disparate role in different tumors, as an oncogene in liver cancer, lung cancer, colon adenocarcinoma, and stomach adenocarcinoma, whereas an antioncogene in skin tumors, and thyroid cancer [51]. The reason for the opposite function of HUWE1 in tumor may be the different tumor microenvironment and genetic backgrounds. The function and mechanism of these five mitophagy-related genes in various tumors had been overwhelmingly illustrated. However, the further study is urgently needed to be conducted in MM.

GSEA results indicated oncogene-related gene set, including MYC targets, E2F target, and mTORC1 signaling was significantly associated with high-risk group. Our findings were in line with the lower rapamycin IC50 in the high-risk group by subsequently drug sensitivity analysis, revealing the potential role of mTOR inhibitors in MM therapy. In fact, there is a cross-talk between mTOR signaling pathway and ubiquitin proteasome system [52]. The combination of these drugs offers multiple possibilities for treating MM. A clinical trial aims to RRMM (NCT00483262) uncovered the synergistic effect of bortezomib and mTOR inhibitor. In this clinical trial, 14 of 43 MM patients obtained partial response or better [53].

The results of ESTIMATE and ssGSEA disclosed that MM patients with high-risk score suffered lower levels of stromal score, immune score and ESTIMATE score, whereas higher levels of tumor purity. Besides, cell populations participated in anti-tumor immune response, such as, immature B cell, central memory CD4 T cell, type 17 T helper cell, central memory CD8 T cell, CD56 bright killer cell, effector memory CD8 T cell, nature killer cell, and nature killer T cell, were all markedly reduced in high-risk score group. Research has confirmed the essentials of mitochondrial dynamics in immune cells [54]. Paul et al. revealed the elevation of mitophagy facilitated anti-tumor immunity in intestinal epithelial cells by activated CD8 T cell through cross-dressing of dendritic cells [55]. Moreover, the activation of anti-tumor immune cells demands energy. The lessen antitumor activity of immune cells may be associated with mitophagy dysfunction [56]. In addition, the levels of immature dendritic cell, regulatory T cell, MDSC, macrophage and mast cell, belong to immunosuppressive cells, were also decreased in high-risk score group. A study revealed that macrophages was as a harmful prognostic index in innate immunity [57]. To sum up, the results indicated the tight correlation between mitophagy-related risk score and immunosuppression, which may explain the reason of poor prognosis in MM patients with high risk.

New strategies for the treatment of tumors focus on immune checkpoint inhibitors. We ulteriorly probed the correlation between risk score and immune checkpoint. These results disclosed that risk score had negative association with immune checkpoint, including IDO1, CD276, CD86, PD-L1, and PD-L2, which indicated inhibitors targeting checkpoints, such as PD-L1 may be less valid in MM patients in high-risk score group. However, a clinical trial has been reported that RRMM patients with extramedullary disease obtained benefit from the combination of PD-L1 inhibitor avelumab with radiotherapy [58]. Therefore, these findings suggested that mitophagy status and sensibility to immune checkpoints may be heterogeneous between primary and recurrent MM patients.

However, there are still some limitations in our study. First, the expression of five-mitophagy-genes were verified only in mRNA level, and the expression at protein level should be further clarified. Second, the nomogram model was not be applied in validation cohort, due to the lack of clinical information in the two independent datasets. Third, the function of the five-mitophagy-genes in MM has not yet to be cleared. Therefore, further experiments in vitro and in vivo need to be conducted.


In conclusion, we constructed a five-mitophagy-genes (VDAC1, PINK1, VPS13C, ATG13, and HUWE1) prognostic risk model, which as an independent element for MM OS, could estimate the survival of MM accurately and stably both in training and validation cohorts. The molecular landscape characteristics upon the risk score, including the regulatory pathways, immune infiltration level, and potential drug targets, improved our cognition for MM, which provide novel insights into MM treatment.

Data Availability

The public datasets used in this article, including GSE6477, GSE13591, GSE9782, GSE24080, GSE4204, and GSE47552 can be downloaded in the GEO database. The original contributions presented in the study are included in the article/ Supplementary Material. Further inquiries can be directed to the corresponding authors.



Multiple myeloma


Bone marrow microenvironment


Reactive oxygen species


Adenosine triphosphate


Mitochondrial DNA


Least absolute shrinkage and selection operator


Gene Expression Omnibus database


Differentially expressed genes


Overall survival


Protein-Protein interaction


Receiver operating characteristic


Gene ontology


Kyoto Encyclopedia of Genes and Genomes


Gene set enrichment analysis


False discovery rate


Single sample GSEA


Genomics of drugs sensitivity

IC50 :

Half-maximal inhibitory concentration


Real-time polymerase chain reaction


Bone marrow mononuclear cells


Lactate dehydrogenase


Hazard ratio


Confidence interval


Monoclonal gammopathy of unknown significance


Smoldering multiple myeloma


Refractory and/or relapse MM






Beta-2 microglobulin


Voltage-dependent anion channel 1


PTEN-induced kinase 1


Vacuolar protein sorting 13 homolog C


Autophagy related 13


Programmed cell death-ligand 1


Programmed cell death-ligand 2


  1. Kumar SK, Rajkumar V, Kyle RA, Duin MV, Anderson KC. Multiple myeloma. Nat Rev Dis Primers. 2017;3:17046.

    Article  PubMed  Google Scholar 

  2. Natural history of. Relapsed myeloma, refractory to immunomodulatory drugs and proteasome inhibitors: a multicenter IMWG study. Leukemia 2017.

  3. Rajkumar SV, Kumar S. Multiple myeloma current treatment algorithms. Blood Cancer J 2020, 10(9).

  4. Zhao WH, Wang BY, Chen LJ, Fu WJ, Xu J, Liu J, Jin SW, Chen YX, Cao XM, Yang Y et al. Four-year follow-up of LCAR-B38M in relapsed or refractory multiple myeloma: a phase 1, single-arm, open-label, multicenter study in China (LEGEND-2). J Hematol Oncol 2022, 15(1).

  5. Pinto V, Bergantim R, Caires HR, Seca H, Guimaraes JE, Vasconcelos MH. Multiple myeloma: available therapies and causes of Drug Resistance. Cancers 2020, 12(2).

  6. Abdallah N, Rajkumar SV, Greipp P, Kapoor P, Gertz MA, Dispenzieri A, Baughn LB, Lacy MQ, Hayman SR, Buadi FK et al. Cytogenetic abnormalities in multiple myeloma: association with disease characteristics and treatment response. Blood Cancer J 2020, 10(8).

  7. Corre J, Munshi N, Avet-Loiseau H. Genetics of multiple myeloma: another heterogeneity level? Blood 2015, 125(12):1870–6.

  8. Bianchi G, Munshi NC. Pathogenesis beyond the cancer clone(s) in multiple myeloma. Blood. 2015;125(20):3049–58.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  9. De Smedt E, Maes K, Verhulst S, Lui H, Kassambara A, Maes A, Robert N, Heirman C, Cakana A, Hose D, et al. Loss of RASSF4 expression in multiple myeloma promotes RAS-Driven malignant progression. Cancer Res. 2018;78(5):1155–68.

    Article  PubMed  Google Scholar 

  10. Guillerey C, Harjunpaa H, Carrie N, Kassem S, Teo T, Miles K, Krumeich S, Weulersse M, Cuisinier M, Stannard K, et al. TIGIT immune checkpoint blockade restores CD8(+) T-cell immunity against multiple myeloma. Blood. 2018;132(16):1689–94.

    Article  CAS  PubMed  Google Scholar 

  11. Di Marzo L, Desantis V, Solimando AG, Ruggieri S, Annese T, Nico B, Fumarulo R, Vacca A, Frassanito MA. Microenvironment drug resistance in multiple myeloma: emerging new players. Oncotarget. 2016;7(37):60698–711.

    Article  PubMed  PubMed Central  Google Scholar 

  12. Holthof LC, Mutis T. Challenges for Immunotherapy in multiple myeloma: bone marrow microenvironment-mediated Immune suppression and Immune Resistance. Cancers 2020, 12(4).

  13. Richardson DR, Lane DJ, Becker EM, Huang ML, Whitnall M, Suryo Rahmanto Y, Sheftel AD, Ponka P. Mitochondrial iron trafficking and the integration of iron metabolism between the mitochondrion and cytosol. Proc Natl Acad Sci U S A. 2010;107(24):10775–82.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  14. Lemasters JJ. Variants of mitochondrial autophagy: types 1 and 2 mitophagy and micromitophagy (type 3). Redox Biol. 2014;2:749–54.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  15. Chang JY, Yi HS, Kim HW, Shong M. Dysregulation of mitophagy in carcinogenesis and tumor progression. Bba-Bioenergetics. 2017;1858(8):633–40.

    Article  CAS  PubMed  Google Scholar 

  16. Yu MF, Nguyen ND, Huang YQ, Lin D, Fujimoto TN, Molkentine JM, Deorukhkar A, Kang Y, San Lucas FA, Fernandes CJ et al. Mitochondrial fusion exploits a therapeutic vulnerability of pancreatic cancer. Jci Insight 2019, 4(16).

  17. Maes H, Rubio N, Garg AD, Agostinis P. Autophagy: shaping the tumor microenvironment and therapeutic response. Trends Mol Med. 2013;19(7):428–46.

    Article  CAS  PubMed  Google Scholar 

  18. Kanehisa M, Furumichi M, Sato Y, Kawashima M, Ishiguro-Watanabe M. KEGG for taxonomy-based analysis of pathways and genomes. Nucleic Acids Res. 2023;51(D1):D587–92.

    Article  CAS  PubMed  Google Scholar 

  19. Kanehisa M. Toward understanding the origin and evolution of cellular organisms. Protein Sci. 2019;28(11):1947–51.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  20. Kanehisa M, Goto S. KEGG: Kyoto Encyclopedia of genes and genomes. Nucleic Acids Res. 2000;28(1):27–30.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  21. Yoshihara K, Shahmoradgoli M, Martínez E, Vegesna R, Kim H, Torres-Garcia W, Treviño V, Shen H, Laird PW, Levine DA, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun. 2013;4:2612.

    Article  PubMed  Google Scholar 

  22. Jia QZ, Wu W, Wang YQ, Alexander PB, Sun CD, Gong ZH, Cheng JN, Sun HB, Guan YF, Xia XF et al. Local mutational diversity drives intratumoral immune heterogeneity in non-small cell lung cancer. Nat Commun 2018, 9.

  23. Lu XF, Jiang LY, Mang LY, Zhu Y, Hu WJ, Wang JS, Ruan XJ, Xu ZB, Meng XW, Gao J, et al. Immune signature-based subtypes of cervical squamous cell Carcinoma tightly Associated with Human Papillomavirus Type 16 expression, molecular features, and clinical outcome. Neoplasia. 2019;21(6):591–601.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  24. Geeleher P, Cox N, Huang RS. pRRophetic: an R package for prediction of clinical chemotherapeutic response from tumor gene expression levels. (1932–6203 (Electronic)).

  25. Cerami E, Gao J, Dogrusoz U, Gross BE, Sumer SO, Aksoy BA, Jacobsen A, Byrne CJ, Heuer ML, Larsson E, et al. The cBio cancer genomics portal: an open platform for exploring multidimensional cancer genomics data. Cancer Discov. 2012;2(5):401–4.

    Article  PubMed  Google Scholar 

  26. Fitzgerald M, Saville BR, Lewis RJ. Decision curve analysis. Jama-J Am Med Assoc. 2015;313(4):409–10.

    Article  CAS  Google Scholar 

  27. Morgan GJ, Walker BA, Davies FE. The genetic architecture of multiple myeloma. Nat Rev Cancer. 2012;12(5):335–48.

    Article  CAS  PubMed  Google Scholar 

  28. Zanwar S, Nandakumar B, Kumar S. Immune-based therapies in the management of multiple myeloma. Blood Cancer J 2020, 10(8).

  29. Robak P, Drozdz I, Szemraj J, Robak T. Drug resistance in multiple myeloma. Cancer Treat Rev. 2018;70:199–208.

    Article  CAS  PubMed  Google Scholar 

  30. Mizushima N, Levine B. Autophagy in mammalian development and differentiation. Nat Cell Biol. 2010;12(9):823–30.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  31. Maugeri G, D’Amico AG, Magro G, Salvatorelli L, Barbagallo GMV, Saccone S, Drago F, Cavallaro S. D’Agata V: expression profile of parkin isoforms in human gliomas. Int J Oncol. 2015;47(4):1282–92.

    Article  CAS  PubMed  Google Scholar 

  32. D’Amico AG, Maugeri G, Magro G, Salvatorelli L, Drago F, D’Agata V. Expression pattern of parkin isoforms in lung adenocarcinomas. Tumor Biol. 2015;36(7):5133–41.

    Article  Google Scholar 

  33. Tay SP, Yeo CWS, Chai C, Chua PJ, Tan HM, Ang AXY, Yip DLH, Sung JX, Tan PH, Bay BH, et al. Parkin enhances the expression of cyclin-dependent kinase 6 and negatively regulates the proliferation of breast Cancer cells. J Biol Chem. 2010;285(38):29231–8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  34. Zheng ZH, Fan SJ, Zheng J, Huang W, Gasparetto C, Chao NJ, Hu JD, Kang YB. Inhibition of thioredoxin activates mitophagy and overcomes adaptive bortezomib resistance in multiple myeloma. J Hematol Oncol 2018, 11.

  35. Zhang GX, Jiang GX, Wang C, Zhong K, Zhang JJ, Xue Q, Li X, Jin H, Li BL. Decreased expression of microRNA-320a promotes proliferation and invasion of non-small cell lung cancer cells by increasing VDAC1 expression. Oncotarget. 2016;7(31):49470–80.

    Article  PubMed  PubMed Central  Google Scholar 

  36. Fang YT, Liu JP, Zhang QC, She CH, Zheng RJ, Zhang RD, Chen ZX, Chen CF, Wu JD. Overexpressed VDAC1 in breast cancer as a novel prognostic biomarker and correlates with immune infiltrates. World J Surg Oncol 2022, 20(1).

  37. Zhang CL, Ding WC, Liu Y, Hu Z, Zhu D, Wang XL, Yu L, Wang LM, Shen H, Zhang WC, et al. Proteomics-based identification of VDAC1 as a tumor promoter in cervical carcinoma. Oncotarget. 2016;7(32):52317–28.

    Article  PubMed  PubMed Central  Google Scholar 

  38. Pittala S, Krelin Y, Shoshan-Barmatz V. Targeting Liver Cancer and Associated Pathologies in mice with a mitochondrial VDAC1-Based peptide. Neoplasia. 2018;20(6):594–609.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  39. Yang GC, Zhou DW, Li J, Wang W, Zhong W, Fan W, Yu MC, Cheng HT. VDAC1 is regulated by BRD4 and contributes to JQ1 resistance in breast cancer. Oncol Lett. 2019;18(3):2340–7.

    CAS  PubMed  PubMed Central  Google Scholar 

  40. Pandey SK, Machlof-Cohen R, Santhanam M, Shteinfer-Kuzmine A, Shoshan-Barmatz V. Silencing VDAC1 to treat Mesothelioma Cancer: Tumor Reprograming and Altering Tumor Hallmarks. Biomolecules 2022, 12(7).

  41. Eiyama A, Okamoto K. PINK1/Parkin-mediated mitophagy in mammalian cells. Curr Opin Cell Biol. 2015;33:95–101.

    Article  CAS  PubMed  Google Scholar 

  42. Poulogiannis G, McIntyre RE, Dimitriadi M, Apps JR, Wilson CH, Ichimura K, Luo FJ, Cantley LC, Wyllie AH, Adams DJ, et al. PARK2 deletions occur frequently in sporadic colorectal cancer and accelerate adenoma development in apc mutant mice. P Natl Acad Sci USA. 2010;107(34):15145–50.

    Article  CAS  Google Scholar 

  43. Xu Y, Lu JW, Tang YB, Xie WJ, Zhang HT, Wang BB, Zhang SL, Hou WJ, Zou C, Jiang PC, et al. PINK1 deficiency in gastric cancer compromises mitophagy, promotes the Warburg effect, and facilitates M2 polarization of macrophages (retracted article. See vol. 549, 2022). Cancer Lett. 2022;529:19–36.

    Article  CAS  PubMed  Google Scholar 

  44. Fan SJ, Price T, Huang W, Plue M, Warren J, Sundaramoorthy P, Paul B, Feinberg D, MacIver N, Chao N et al. PINK1-Dependent Mitophagy regulates the Migration and Homing of multiple myeloma cells via the MOB1B-Mediated Hippo-YAP/TAZ pathway. Adv Sci 2020, 7(5).

  45. Wei R, Cao J, Yao S. Matrine promotes liver cancer cell apoptosis by inhibiting mitophagy and PINK1/Parkin pathways. (1466 – 1268 (Electronic)).

  46. Li LC, Chen T, Wang JS, Li MX, Li QS. Identification of an Autophagy-Related Signature Based on Whole Bone Marrow Sequencing for the Prognosis and Immune Microenvironment Characterization of Multiple Myeloma. J Immunol Res 2022, 2022.

  47. Wu XP, Zhao JM. Novel oxidative stress-related prognostic biomarkers for melanoma associated with tumor metastasis. Medicine 2021, 100(8).

  48. Karabi Z, Moradian F, Kheirabadi M. The effect of lactoferrin on ULK1 and ATG13 genes expression in breast cancer cell line MCF7 and bioinformatics studies of protein interaction between lactoferrin and the autophagy initiation complex. Cell Biochem Biophys. 2022;80(4):795–806.

    Article  CAS  PubMed  Google Scholar 

  49. Zhou J, Shen R, Zhou HB, Huang J. OBHS impairs the viability of breast cancer via decreasing ERα and Atg13. (1090–2104 (Electronic)).

  50. Niu Q, Zhao W, Wang J, Li CM, Yan T, Lv W, Wang GJ, Duan WH, Zhang T, Wang KN, et al. LicA induces autophagy through ULK1/Atg13 and ROS pathway in human hepatocellular carcinoma cells. Int J Mol Med. 2018;41(5):2601–8.

    CAS  PubMed  PubMed Central  Google Scholar 

  51. Gong XF, Du DY, Deng YR, Zhou YQ, Sun L, Yuan ST. The structure and regulation of the E3 ubiquitin ligase HUWE1 and its biological functions in cancer. Invest New Drug. 2020;38(2):515–24.

    Article  CAS  Google Scholar 

  52. Eichner R, Fernandez-Saiz V, Targosz BS, Bassermann F. Cross Talk Networks of mammalian target of Rapamycin Signaling with the Ubiquitin Proteasome System and their clinical implications in multiple myeloma. Int Rev Cel Mol Bio. 2019;343:219–97.

    Article  CAS  Google Scholar 

  53. Ghobrial IM, Weller E, Vij R, Munshi NC, Banwait R, Bagshaw M, Schlossman R, Leduc R, Chuma S, Kunsman J, et al. Weekly bortezomib in combination with temsirolimus in relapsed or relapsed and refractory multiple myeloma: a multicentre, phase 1/2, open-label, dose-escalation study. Lancet Oncol. 2011;12(3):263–72.

    Article  CAS  PubMed  Google Scholar 

  54. Xie JH, Li YY, Jin J. The essential functions of mitochondrial dynamics in immune cells. Cell Mol Immunol. 2020;17(7):712–21.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  55. Ziegler PK, Bollrath J, Pallangyo CK, Matsutani T, Canli O, De Oliveira T, Diamanti MA, Muller N, Gamrekelashvili J, Putoczki T, et al. Mitophagy in Intestinal epithelial cells triggers adaptive immunity during tumorigenesis. Cell. 2018;174(1):88–.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  56. O’Sullivan TE, Johnson LR, Kang HH, Sun JC. BNIP3-and BNIP3L-Mediated Mitophagy promotes the generation of natural killer cell memory. Immunity. 2015;43(2):331–42.

    Article  PubMed  PubMed Central  Google Scholar 

  57. Zhou SL, Zhou ZJ, Hu ZQ, Huang XW, Wang Z, Chen EB, Fan J, Cao Y, Dai Z, Zhou J. Tumor-Associated Neutrophils Recruit Macrophages and T-Regulatory cells to promote progression of Hepatocellular Carcinoma and Resistance to Sorafenib. Gastroenterology. 2016;150(7):1646–.

    Article  CAS  PubMed  Google Scholar 

  58. Kazandjian D, Dew A, Hill E, Ramirez EG, Morrison C, Mena E, Lindenberg L, Yuan C, Maric I, Wang HW, et al. Avelumab, a PD-L1 inhibitor, in combination with Hypofractionated Radiotherapy and the Abscopal Effect in Relapsed Refractory multiple myeloma. Oncologist. 2021;26(4):288–.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

Download references


We appreciate the generosity of GEO database for sharing the huge amount of data.


This study was supported by the National Key Research and Development Program of China grant nos. 2019YFA0801804; the National Natural Science Foundation of China, grants Nos. 82170140, 81922003 and 82100150; Natural Science Foundation of China Shaanxi Province, grant No. 2020JQ-558.

Author information

Authors and Affiliations



KGY and HAL conceived the project, supervised, and designed all research. JYC, LR and KGY wrote the manuscript. JYC, FYD and LR designed and conducted the bioinformatic analysis. JYC and SLY performed the experimental work. ZLL and GN revised the manuscript.

Corresponding authors

Correspondence to Aili He or Guangyao Kong.

Ethics declarations

Competing interests

The authors declare no competing interests.

Ethics approval and consent to participate

All samples from 16 healthy donors and 16 newly diagnosed MM (NDMM) patients were purified from bone marrow aspirates from 2021 to 2022. The study was conducted in accordance with the declaration of Helsinki, and the study protocol was approved by the Ethics Committee of the second affiliated hospital of Xi’an Jiaotong University, Shaanxi, China (case No. 2022128). Written informed consent was obtained from each study participant.

Consent for publication

Not applicable.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Electronic supplementary material

Below is the link to the electronic supplementary material.

Supplementary Material 1

Supplementary Material 2

Supplementary Material 3

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Jia, Y., Liu, R., Shi, L. et al. Integrative analysis of the prognostic value and immune microenvironment of mitophagy-related signature for multiple myeloma. BMC Cancer 23, 859 (2023).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: