Identification of an Autophagy-Related Signature Based on Whole Bone Marrow Sequencing for the Prognosis and Immune Microenvironment Characterization of Multiple Myeloma

Li, Licheng; Chen, Ting; Wang, Jishi; Li, Mengxing; Li, Qinshan

doi:https://doi.org/10.1155/2022/3922739

Journal of Immunology Research

On this page

Abstract Introduction Materials and Methods Results Discussion Conclusions Data Availability Conflicts of Interest Authors’ Contributions Acknowledgments Supplementary Materials References Copyright Related Articles

Special Issue

Exploration of Novel Potential Biomarkers for Cancer Immunotherapy

View this Special Issue

Research Article | Open Access

Volume 2022 | Article ID 3922739 | https://doi.org/10.1155/2022/3922739

Identification of an Autophagy-Related Signature Based on Whole Bone Marrow Sequencing for the Prognosis and Immune Microenvironment Characterization of Multiple Myeloma

Licheng Li,^1,2Ting Chen,^1,2Jishi Wang,^1,2Mengxing Li ,^1,2and Qinshan Li^3,4

Academic Editor: Fu Wang

Received11 Mar 2022

Revised21 Apr 2022

Accepted04 May 2022

Published29 May 2022

Abstract

Myeloma (MM) is a malignant plasma cell disorder, which is incurable owing to its drug resistance. Autophagy performs an integral function in homeostasis, survival, and drug resistance in multiple myeloma (MM). Therefore, the purpose of the present research was to identify potential autophagy-related genes (ARGs) in patients with MM. We downloaded the transcriptomic data (GSE136400) of patients with MM, as well as the corresponding clinical data from the Gene Expression Omnibus (GEO); the patients were classified at random into two groups in a ratio of 6: 4, with 212 samples in the training dataset and 142 samples in the test dataset. Both multivariate and univariate Cox regression analyses were performed to identify autophagy-related genes. The univariate Cox regression analysis demonstrated that 26 ARGs had a significant correlation with overall survival (OS). We constructed an autophagy-related risk prognostic model based on six ARGs: EIF2AK2 (ENSG00000055332), KIF5B (ENSG00000170759), MYC (ENSG00000136997), NRG2 (ENSG00000158458), PINK1 (ENSG00000158828), and VEGFA (ENSG00000112715) using LASSO-Cox regression analysis to predict risk outcomes, which revealed substantially shortened OS duration in the high-risk cohort in contrast with that in the low-risk cohort. Therefore, the ARG-based model significantly predicted the MM patients’ prognoses and was verified in an internal test set. Differentially expressed genes were found to be predominantly enriched in pathways associated with inflammation and immune regulation. Immune infiltration of tumor cells resulted in the formation of a strong immunosuppressive microenvironment in high-risk patients. The potential therapeutic targets of ARGs were subsequently analyzed via protein–drug network analysis. Therefore, a prognostic model for MM was established via a comprehensive analysis of ARGs, through using the clinical models; we have further revealed the molecular landscape features of multiple myeloma.

1. Introduction

Multiple myeloma (MM) is a clonal proliferative disease of plasma cells whose characteristics include abnormal immunoglobulin production proliferation of malignant plasma cells into the bone marrow, and bone destruction [1]. Its annual incidence is 4/100,000 people, accounting for approximately 10% of all patients with hematological diseases [2]. Although various treatments including protease inhibitors and immune targeted therapy have been used in clinical practice, however, high recurrence rates and drug resistance rates still affect the prognosis of patients, and the pathogenesis has not been elucidated. Therefore, it is essential to identify new biomarkers for prognosis and understand the pathogenesis of MM.

Autophagy is type II programmed cell death, which is a process of self-digestion and catabolism in cells [3]. Recent evidence indicates that autophagy is a key mechanism for leukemia development and chemotherapy resistance [4]. In addition, autophagy is an essential pro-survival mechanism for medication resistance in MM. The combination therapy of autophagy blockers and traditional anti-MM therapy has enhanced their effects on drug-resistant MM plasma cells, thus making autophagy a novel treatment target [5]. Studies on the function of autophagy in MM are mainly focused on the regulatory mechanisms of several autophagy-related pathways [6]. As far as we know, no research has probed into the potential association between the expression of autophagy-related genes (ARGs) and the MM patients’ prognoses.

In the present research, we examined the relationship between ARGs and the prognosis of patients with MM and established a prognostic model for MM. The transcriptome sequencing data of MM patients were obtained from the GEO database, and the patients were classified into the training and test sets. LASSO regression was performed to screen for the key risk factors of patients with MM, construct an ARG-based prognostic risk model, and validate the accuracy of the model in an external verification dataset. The risk model was found to be an independent and effective predictor of prognosis and helped to characterize the MM-specific microenvironment landscape for tumor-infiltrating immune cells and identify potential drug-targeted ARGs. Therefore, the risk model helped to understand the prognostic characteristics and molecular mechanisms of ARGs in MM.

2. Materials and Methods

2.1. Data Acquisition and Processing

The transcriptome data of 1424 patients were downloaded from a public database, including 354 whole bone marrow samples and pretreated sequencing samples. The raw dataset GSE136400 [7] was acquired from GEO (https://www.ncbi.nlm.nih.gov/geo/), and arrays were processed using the robust multi-array averaging (RMA) algorithm that was performed with the aid of the affy R package [8]. Batch effects were equalized via ComBat analysis using the sva package [9]. The sequencing datasets of patients with MM were obtained from the UCSC Xena database (https://xenabrowser.net/datapages/).

2.2. Screening for ARGs

In the Human Autophagy Database (HADb, http://autophagy.lu/clustering/index. html), 209 ARGs were annotated, and 10 genes with expression values of 0 were excluded. The expression profile generated from GSE136400 contained a sum of 199 ARGs.

2.3. Screening for ARGs Related to Overall Survival

A total of 354 pretreated sequencing samples of whole bone marrow were randomly divided into two groups in a ratio of 6 : 4, with the training set comprising of 212 patients and the test set comprising of 142 patients. Patients in the training set were screened to elucidate the underlying risk characteristics of the identified ARGs in MM. Eventually, overall survival– (OS–) related ARGs with were selected utilizing a univariate Cox hazard regression analysis.

2.4. Enrichment Analyses of OS-Related ARGs

The Gene Ontology (GO) analysis and the Kyoto Encyclopaedia of Genes and Genomes (KEGG) enrichment analysis were conducted utilizing the ClusterProfiler package [10], which revealed the particular roles of ARGs associated with the OS in MM. The Benjamini–Hochberg with adjusted was set for the purpose of determining statistical significance.

2.5. Characteristics of Molecular Interactions of OS-Related ARGs

To explore the relationship among OS-related ARGs, we employed the STRING database to establish a protein-protein interaction (PPI) network [11] and visualized it using Cytoscape (version: 3.8.0) [12] to calculate the topology level for each module.

TRRUST, an online tool for exploring the relationship between human and mouse transcriptional regulatory networks [13], was used to screen for transcription factors (TFs) associated with key OS-related ARGs.

2.6. Establishment and Verification of an ARG-Related Prognostic Model for MM

To prevent these performances to be over-fitted, we used LASSO-Cox regression [14] analysis in the training set to determine key OS-related ARGs for establishing a prognostic model. Multivariable Cox proportional hazard regression analysis was conducted on these ARGs, and variables were gradually selected according to the Akaike information criterion (AIC) [15]. The following is the equation for calculating the final risk score based on the optimized prognostic features:

In this equation, denotes the regression coefficient, denotes the ARGs used to construct the signature, denotes the relative expression value of each ARG in the signature, whereas signifies the sum of genes within the signature. Based on the median risk score, patients were classified into 2 cohorts: low- and high-risk cohorts. For the purpose of comparing the two cohorts, the log-rank test and the Kaplan–Meier analysis were performed. The prediction accuracy of the signature was determined using the time-dependent receiver operating characteristic (ROC) curves [16].

To verify the prediction accuracy of the prognostic model, 354 whole bone marrow pretreatment sequencing samples were randomly divided into two groups in a ratio of 6 : 4. The verification group consisted of a sum of 142 patients. The abovementioned equation was employed to derive each patient’s risk score, and Kaplan–Meier curves were drawn to reassess differences in prognosis and survival across the two cohorts.

2.7. Enrichment Analyses and the Determination of Differentially Expressed Genes (DEGs)

The limma software package [17] was used to identify DEGs in the high- and low-risk groups. GO and KEGG enrichment analyses were conducted with the aid of the ClusterProfiler software package [10], which provides information regarding the GO terms cell component (CC), molecular function (MF), biological process (BP), and KEGG pathways, thus understanding the potential functions of DEGs in MM. In addition, ClustVis [18] was used to cluster the DEGs and generate heat maps.

2.8. Gene Set Variation Analysis

As mentioned above, the low- and high-risk cohorts were established based on the median risk score. Gene set variation analysis (GSVA) [19] was used to further confirm the main enrichment pathways of the low- and high-risk cohorts. Nominal and indicated significant enrichment.

2.9. Validation of Prognostic Accuracy of OS-Related ARGs in Clinicopathological Nomograms

The rms software package [20] was used to establish a quantitative prediction technique for the prognosis of patients with MM, and the prognostic features were included in the clinicopathological analyses of the training group. Subsequently, the final model was selected according to the AIC, and calibration curves were charted for the purpose of examining the predictive accuracy regarding the nomogram [21].

2.10. Immune Microenvironment Landscape of MM and Possible Immunotherapy Targets for Predicting Prognosis

CIBERSORT was employed to analyze the infiltration levels of 22 distinct immune cells in the low- and high-risk cohorts [22], such as NK cells, plasma cells, B cells, and T cells. The criterion of was defined as having statistical significance and was used for subsequent analyses as well.

Feasible therapeutic methods for people with MM, immunotherapy, as well as molecular targeted therapy and immunotherapy, have gained increasing attention over the years [23]. We used Spearman correlation analysis to analyze the relationship between the prognostic risk score and treatment goals in clinical practice, and the treatment-related targets were as follows: CTLA-4(ENSG00000163599), VTCN1 (ENSG00000134258), TNFSF18 (ENSG00000120337), TNFSF15 (ENSG00000181634), TNFRSF25 (ENSG00000215788), TNFRSF9 (ENSG00000049249), TNFRSF8 (ENSG00000120949), TMIGD2 (ENSG00000167664), CD276 (ENSG00000103855), PDCD1LG2 (PD-L2)(ENSG00000197646), PDCD1 (PDL1)(ENSSSAG00000045082), LAIR1 (ENSG00000167613), LAG3 (ENSG00000089692), and HAVCR2 (TIM3) (ENSG00000135077).

To screen for potential drug targets of these OS-related ARGs, we employed the NetworkAnalyst 3.0 (https://www.networkanalyst.ca/) to analyze the protein-drug interaction in OS-related ARGs. Information regarding the target and drug content was obtained from the DrugBank (version 5.0, https://go.drugbank.com/) [24].

3. Results

3.1. Screening of OS-Related ARGs in MM

ARGs whose P value was less than 0.05 were chosen for the purpose of conducting further investigations by means of a univariate Cox proportional hazard regression analysis so as to determine the relationship between each identified ARG and the MM patients’ prognoses. We found that 26 ARGs were correlated with OS (Table 1), and GO analysis of these ARGs revealed that they were mainly enriched in macrophage polarization and regulation of autophagy (Figure 1(a)). In addition, these ARGs were involved in autophagy (animals), apoptosis, PI3K–AKT signaling pathway, and other pathways (Figure 1(b)).

(a)

(b)

(c)

To examine the interaction among OS-related ARGs, a PPI network was established, which identified two key modules, namely, CASP3 and TP53, with the highest topology levels (Figure 1(c)). The TP53 module comprised 13 nodes and 26 edges, whereas the CASP3 module consisted of PINK1 and WIPI2 nodes. The regulatory mechanisms of these ARGs may be critical in the pathogenesis of MM.

Furthermore, the TRRUST database (Table S1) was used to identify 42 TFs that regulated these OS-related ARGs and included nuclear TFs (NF-κB1, E2F1, and SP1), histone deacetylase family (HDAC2, HDAC3, and SIRT1), TP53, and NF-κB signal transduction core factors (NF-κB1 and RELA). In addition, important ARG node factors, such as CASP3, MYC, and TP53, were significantly regulated by these TFs.

3.2. Screening and Verification of OS-Related ARGs for Prognosis

To avoid potential over-fitting, hub OS-related ARGs were selected for establishing an interaction network using LASSO-Cox regression analysis (Figures S1A and S1B). Six ARGs were identified via the stepwise multivariate Cox proportional hazard regression analysis and were applied to develop a prognostic model for patients with MM (Figure 2(a)). The following is the equation for calculating the risk score of each individual patient:

(a)

(b)

(c)

(d)

According to their median risk score, all patients were classified into two cohorts, namely, the high-risk cohort and the low-risk cohort. The death rate of patients increased with the increasing risk score (Figures 2(b) and 2(c)). In addition, six OS-related ARGs, namely, EIF2AK2, KIF5B, MYC, NRG2, PINK1, and VEGFA, were found to have an elevated expression in the high-risk cohort (Figure 2(d) and Figure S1C).

3.3. Evaluation of the Prognostic Features of OS-Related ARGs in Patients with MM

An analysis of Kaplan–Meier survival data was performed for the purpose of verifying the prediction accuracy of the risk signature. The results indicated that the patients’ overall survival (OS) in the high-risk cohort was remarkably shorter in contrast with that in the low-risk cohort (, Figure 3(a)). Furthermore, the AUC values for anticipating OS over one, three, and five years were shown to be 0.610, 0.704, and 0.752, indicating that the signature had a high predictive power (Figure 3(b)).

(a)

(b)

(c)

(d)

(e)

(f)

Figure 3

Evaluation of the prognostic features of OS-related ARGs in multiple myeloma patients. (a) Survival status of the low- and high-risk cohorts. (b) ROC curves for anticipating OS over one, three, and five years. (c) Univariate Cox regression analysis of clinical-pathological variables and risk score. (d) Multivariable Cox regression analysis of clinical characteristics and risk score. (e) Nomograms to predict the MM patients’ OS over one, three, and five years by combining ISS stage, age, and risk score. (f) The accuracy of the predicted survival rates over one, three, and five years was validated using calibration curves of nomograms. The dashed line indicates an ideal nomogram, whereas the green, blue, and red solid lines refer to the actual utility of the nomogram.

A univariate Cox regression analysis was performed for the purpose of examining if the risk signature independently served as a prognostic marker. It was demonstrated that age, β₂-microglobulin level, ISS stage, chromosomal karyotype abnormality, and risk scores were substantially associated with OS (Figure 3(c)). In addition, multivariable Cox regression analysis revealed that after adjusting these clinicopathological factors, the risk signature remained to be a prognostic indicator of MM in an independent manner (Figure 3(d)). Moreover, age and ISS stage were also identified as independently serving as prognostic predictors. Furthermore, the predictive ability of the risk signature was compared among patients with MM according to the risk scores and clinical characteristics. The AUC values for predicting the OS over three, five, and ten years were higher for the risk score than for clinical variables (Figures S2B, S2C, and S2D). However, the AUC value for predicting 1-year OS was lower for the risk score than for some clinical variables (Figure S2A). This finding indicated that the established risk signature is very important for long-term survival prognosis.

To evaluate the signature in a highly accurate manner, a nomogram that combined the risk score, age, and ISS stage was plotted (Figure 3(e)), which demonstrated excellent concordance in predicting OS over one, three, and five years as revealed by calibration plots (Figure 3(f)). The results suggested that the risk score, combined with clinical variables, could anticipate the MM patients’ OS in a highly accurate manner.

3.4. Enrichment Analysis of ARGs

After observing differences in OS between the low- and high-risk cohorts, we performed GSVA to detect biological functions and related signaling pathways. Humoral immune response and cell cycle pathways exhibited a remarkable enrichment in the low-risk cohort (Figures S3C and S3D), whereas the ubiquitin-mediated proteolysis and endoplasmic reticulum-related degradation (ERAD) pathways were shown to exhibit remarkable enrichment in the high-risk cohort (Figures S3A and S3B). Previous studies have shown that ERAD interacts with the ubiquitin-proteasome pathway and autophagy to reduce protein misfolding or its consequences in a coordinated manner [25]. Therefore, the results indicated that autophagy-related events were associated with high-risk patients with MM.

3.5. Validation of the Prognostic Characteristics of the Model

The risk score formula was used to calculate the risk scores of patients belonging to the low- and high-risk cohorts. The findings revealed that the OS of patients was shorter and their prognosis was worse in the high-risk cohort as opposed to the low-risk cohort (Figure 4(a)). In the validation set, the AUC values for predicting the OS one, two, and three years were 0.648, 0.629, and 0.533, correspondingly (Figure 4(b)). The above data indicated that the risk score independently served as a factor for predicting the MM patients’ OS profiles.

(a)

(b)

3.6. Enrichment and Identification of DEGs

DEGs were screened in the training and test groups, and their potential functions were assessed. We identified a sum of 56 DEGs in the training set using the limma software package. These DEGs had distinct gene expression patterns in the low- and high-risk cohorts. Of these, 50 genes were upregulated, whereas 6 genes were downregulated (Figures 5(a) and 5(b)). Further, GO analysis illustrated that the ARGs were involved in several biological functions such as regulation of epithelial morphogenesis, endoplasm reticulum lumen, and receptor-ligand activity (Figure 5(c)). KEGG analysis illustrated that the DEGs were primarily involved in the Hippo, Wnt, and TGF-β signaling pathways and intentional immune network for IgA production (Figure 5(d)). Considering these crucial biological implications, we inferred that the DEGs played a role in the progression and immunoregulation of MM.

(a)

(b)

(c)

(d)

3.7. Correlation between Prognostic ARGs and the Immune Microenvironment of MM

The interaction between tumor-infiltrating immune cells (TILs) and cancer cells in the tumor microenvironment (TME) is important for cancer progression and drug resistance. Therefore, we assessed the correlation between the OS-related ARG-based risk signature and infiltration of immune cells in TME in the training set and found that the infiltration levels of monocytes, resting mast cells, and neutrophils were remarkably elevated in the low-risk cohort, whereas those of memory B cells, eosinophils, activated mast cells, activated dendritic cells, and plasma cells were increased in the high-risk cohort (Figures 6(a) and 6(b)). These findings illustrated that the risk signature exhibited a significant correlation with TILs present in the microenvironment of MM.

(a)

(b)

3.8. Relationship between the Risk Signature and Immune Checkpoint Therapy

Immune checkpoint inhibitors are a new type of targeted immunotherapy, which has recently been used in the preclinical treatment and clinical trials of MM [26]. We used Spearman correlation analysis to examine the relationship between the risk score and various immune checkpoints, including CD274 (PD-L1), CTLA-4, LAG3, and HAVCR2 (TIM3). In addition, differential and correlation analyses were performed for 48 possible checkpoints for tumor-targeted therapy. In addition, the expression of 21 immune factors was found to differ significantly among the risk subgroups (Figures 7 and S4). These results illustrated that the risk score was strongly associated with immune checkpoint genes and tumor immune infiltration. Therefore, patients are more likely to benefit from treatment with immune checkpoints including CTLA-4, PD-L1, LAG3, and HAVCR2.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

(l)

(m)

(n)

(o)

3.9. Identification of Potential Drug Targets of OS-Related ARGs

To screen for potential drug targets of survival-related ARGs, data regarding the protein-drug interactions of these ARGs were obtained from DrugBank and analyzed using NetworkAnalyst. Among these ARGs, the HSP90AB1 gene was identified as a target.

(Table 2 shows drugs with potential targeting of the HSP90AB1 gene in DrugBank).

4. Discussion

MM has been identified as the second most prevalent hematological cancer in adults whose characteristics include the buildup of malignant plasma cells within the bone marrow [27]. Emerging therapies have been used in recent years; however, owing to high recurrence and drug resistance rates, their curative effects on patients and the prognosis of patients are poor. In the present research, an accurate prognosis prediction model was established, which may help clinicians in designing individualized treatment strategies for patients. We identified ARGs with a substantial association with the MM patients’ OS. Moreover, based on six prognosis-related ARGs, we developed a prognostic risk-score model to predict OS. Analysis of the validation set further suggested that the prognostic model is an independent, stable indicator of OS. Moreover, immune invasion analysis revealed the relationship between the risk score and TILs in the low- and high-risk cohorts, revealing a unique autophagic molecular landscape of the immune microenvironment in MM. Finally, the potential therapeutic targets of MM were identified.

Autophagy is a dynamic and continuous metabolic process. Dysregulation of autophagy is considered a mechanism of survival and drug resistance in many types of cancers. Cancer can upregulate autophagy to promote the survival of cancer cells in TME and increase growth and aggression. Autophagy augments the stemness by degrading ubiquitinated p53, thus relieving the tumor suppressor activity of p53 [28]. Several ATG proteins, as well as their corresponding core complexes, such as LC3 ubiquitin-like binding system, ATG9A transport system, ATG12, autophagy-specific class III PI3K complex, and ULK 1/2 kinase core complex, are the key regulatory factors of tumor occurrence and progression. These regulatory factors targeting autophagy-related key proteins may function as effective intervention approaches for cancer management [6]. In addition, regulators targeting autophagy may represent a promising treatment strategy for myeloma [29]. The use of OS-related ARGs and stratification of high- and low-risk patients help to identify potential prognostic markers for the recurrence and drug resistance of MM.

In this study, 26 ARGs were identified in patients with MM via univariate cox regression analysis, and CASP3 and TP53 were identified as the main modules of these ARGs via protein interaction analysis. Certain ARGs have been involved in the occurrence and progression of cancers and medication resistance to hematologic malignancies. For instance, CASP3 can control the development of AML1-ETO-driven leukemia in a ULK1-dependent manner [30]; TP53 is a well-recognized tumor suppressor gene and mutated TP53 serves as a regulatory factor and target of tumor autophagy [31]. In MM, autophagy is regulated via activation of the p53 pathway and inhibition of CK1 [32]. Furthermore, myc is another well-known oncogene and an autophagy regulator, which can promote the progression and drug resistance of MM [33]. Some other OS-related ARGs, such as DAPK2, EEF2, ITPR1, and VEGFA, have not yet been reported to regulate autophagy in MM; however, they have been confirmed to regulate the progression of other tumors through autophagy [34–37]. However, their specific regulatory mechanisms should be investigated further. In this study, we found that the identified ARGs were significantly enriched in pathways related to the regulation of autophagy, apoptosis, and PI3K–AKT. In addition, PI3K–AKT–mTORC is considered a key factor for the negative regulation of autophagy initiation [38], indicating that these ARGs rely on the PI3K–AKT pathway to regulate autophagy and are involved in the progression of MM.

The establishment of tumor-related prognostic models based on the characteristics of transcriptome expression profiles shows a strong application prospect for risk assessment [39]. In the present research, key OS-related ARGs were identified via LASSO-Cox regression analysis for establishing a risk model. Six ARGs with the best prognostic characteristics were used to construct the model. In both multivariate and univariate Cox regression analyses, adopted for adjustment analysis of clinical variables, the risk score was found to independently serve as a prognostic indicator for MM patients, and the AUC values for predicting OS over one, three, and five years were 0.610, 0.704 and 0.752, correspondingly. Moreover, the AUC value for anticipating 1-year OS was higher for a clinical variable than for the risk score, which suggested that the risk score was extremely significant for long-term survival prognosis. A nomogram integrating ISS stage, age, and risk scores and the calibration curves demonstrated that the model has the potential for accurately anticipating the patients’ prognoses. The prediction accuracy of the six OS-related ARGs was similar to those of the five prognostic characteristic genes identified based on RNA sequencing in a study [40]. In addition, the characteristics of the immune microenvironment were explored, and the potential targets of immunotherapy for MM were identified, which may promote the development of drug-based immunotherapy.

Furthermore, DEGs were significantly enriched in the Hippo, Wnt, and TGF-β pathways, which are critical signaling pathways that regulate the infiltration of immune cells in TME. For example, the Hippo pathway effector Yap inhibits T cell function and invasion in TME [41], and the Wnt/β-catenin pathway is the best-characterized pathway that enhanced the progression of cancer by modulating the tumor immune cycle in the majority of the nodes, such as tumor cells, dendritic cells, and T cells [42]. In addition, blocking the TGF-β signal in CD4⁺T cells can reshape TME and inhibit cancer progression [43]. TILs are one of the primary constituents of TME, and the density and types of TILs have marked prognostic associations in multiple cancers [44]. Long-lived plasma cells as well as memory B cells need autophagy in order to survive in MM and perform an integral function in maintaining the survival and medication resistance of tumor cells [45]. In this study, immune infiltration analysis revealed that high-risk patients had high levels of infiltrating eosinophils, activated mast cells, activated dendritic cells, plasma cells, and memory B cells. This finding suggests that a persistently high level of autophagy may be responsible for the intense immunosuppressive microenvironment in MM.

There has been the successful use of immune checkpoint inhibitors in the treatment of hematological cancers and a variety of solid tumors. However, achieving durable response rates in some cancer patients remains challenging because of immune evasion and acquired resistance [46]. In a study, the expression levels of immunosuppressive molecules in eight common tumors exhibited a negative association with the signature-defined risk score, including the classical checkpoints such as PD-L1 and LAG3; however, inhibitors targeting these checkpoints may be less effective in treating patients. Blocking of LAG3 has been shown to significantly enhance T cell activity, kill MM cells, and improve the prognosis of patients with MM [47]. Therefore, continuous activation or inhibition of autophagy may demonstrate the immunosuppressive status of MM patients and their sensitivity to immune-targeted therapy; however, the potential relationship between these two regulatory mechanisms should be investigated further.

Potential strategies to address new drug development challenges and enhance the therapeutic capabilities of existing drugs aim to improve drug target selection and repurpose the existing approved drugs. In this study, the protein-drug interaction analysis of prognosis-related ARGs using the DrugBank database identified only the HSP90AB1 gene as a drug target. There are 14 types of drugs that can target HSP90AB1, and HSP90AB1 has been positively correlated with long non-coding RNA MALAT1 in MM; however, high MALAT1 levels have been correlated with shorter overall progression-free survival [48]. Preclinical studies have shown that HSP90 inhibitors combined with bortezomib enhance anti-MM activity [49]. Therefore, HSP90AB1 is a promising therapeutic target for MM, and further studies should be conducted to develop the next generation of HSP90 inhibitors with higher efficacy and lower toxicity.

The study constructed a characteristic risk prognostic model based on the analysis of autophagy and multiple myeloma patients. We have demonstrated that our models can accurately predict prognosis, but they do have some limitations. We conducted our study entirely using bioinformatics, without an independent external validation cohort. For our study, wet laboratory validation was also necessary to provide additional explanation details. Furthermore, retrospective, published datasets have confirmed that ARGs contribute to disease progression and prospective data are necessary to validate their clinical value and assess their potential clinical relevance.

5. Conclusions

In the present research, we established a prognostic model for MM consisting of six ARGs, which was found to independently function as a predictor of OS in both training and verification sets. In addition, we revealed the molecular landscape characteristics of MM, including the cross-correlation among regulatory pathways, TME, and potential drug targets, which were verified based on several aspects. The prognostic model was verified in publicly available retrospective datasets; however, the verification of its potential clinical value requires further investigation using more prospective data, and the specific molecular processes warrant thorough experimental research.

Data Availability

The manuscript contains the actual data provided in the present research; additional inquiries can be sent to the respective authors.

Conflicts of Interest

No potential conflict of interest was reported by the authors.

Authors’ Contributions

LL and ML were assigned in the conceptualization and design. QL contributed vital ideas. LL and TC were assigned in the data acquisition. LL, TC, and ML contributed in the methodology. LL and TC contributed in the data analysis and interpretation. LL was assigned in the writing (original draft). LL, ML, and JW were assigned in the writing (review and editing). All authors contributed to the article and approved the submitted version.

Acknowledgments

We hereby express our gratitude to everyone who contributed to the public datasets used in the present research.

Supplementary Materials

Supplementary 1. Supplementary Figure 1: LASSO regression analysis is performed to identify crucial variables. (a) 1000-fold cross-validation for variable selection in LASSO regression. (b) Coefficients of LASSO regression for significant autophagy-related genes. Each curve corresponds to a gene involved in autophagy. (c) Six model genes were expressed in low- and high-risk cohorts ().

Supplementary 2. Supplementary Figure 2: different multiple receivers operating characteristic (ROC) curves for risk score and clinical factors. (a) Multiple receivers operating characteristic curves for estimating one-year OS. (b) Multiple ROC curves for anticipating the OS over three years. (c) Multiple ROC curves for anticipating the OS over 5 years. (d) Multiple ROC curves for anticipating the OS over 10 years.

Supplementary 3. Supplementary Figure 3: analysis of the gene set enrichment regarding the autophagy-related signature in low- and high-risk cohorts. (a, b) Genes enriched in the high-risk cohort. (c, d) Genes enriched in the low-risk cohort.

Supplementary 4. Supplementary Figure 4: the tumor immune microenvironment analysis in the low- and high-risk cohorts. (a) Differential expression of 21 types of immune checkpoints in the low- and high-risk cohorts. (b) Analysis of the correlation between the risk score and multiple immune checkpoints.

Supplementary 5. Supplementary Table 1: identification of the transcriptional regulators of OS-associated ARGs.

Supplementary 6. Supplementary Table 2: characteristics of all patients.

References

A. J. Cowan, D. J. Green, M. Kwok et al., “Diagnosis and management of multiple myeloma,” Journal of the American Medical Association, vol. 327, no. 5, pp. 464–477, 2022.
View at: Publisher Site | Google Scholar
J. Liu, Union for China Lymphoma Investigators of the Chinese Society of Clinical Oncology, W. Liu et al., “Incidence and mortality of multiple myeloma in China, 2006–2016: an analysis of the Global Burden of Disease Study 2016,” Journal of Hematology & Oncology, vol. 12, no. 1, p. 136, 2019.
View at: Publisher Site | Google Scholar
I. Dikic and Z. Elazar, “Mechanism and medical implications of mammalian autophagy,” Nature Reviews. Molecular Cell Biology, vol. 19, no. 6, pp. 349–364, 2018.
View at: Publisher Site | Google Scholar
W. Du, A. Xu, Y. Huang et al., “The role of autophagy in targeted therapy for acute myeloid leukemia,” Autophagy, vol. 17, no. 10, pp. 2665–2679, 2021.
View at: Publisher Site | Google Scholar
J. Xu, Y. Su, A. Xu et al., “miR-221/222-mediated inhibition of autophagy promotes dexamethasone resistance in multiple myeloma,” Molecular Therapy, vol. 27, no. 3, pp. 559–570, 2019.
View at: Publisher Site | Google Scholar
X. Li, S. He, and B. Ma, “Autophagy and autophagy-related proteins in cancer,” Molecular Cancer, vol. 19, no. 1, p. 12, 2020.
View at: Publisher Site | Google Scholar
S. A. Danziger, M. McConnell, J. Gockley et al., “Bone marrow microenvironments that contribute to patient outcomes in newly diagnosed multiple myeloma: a cohort study of patients in the Total therapy clinical trials,” PLoS Medicine, vol. 17, no. 11, article e1003323, 2020.
View at: Publisher Site | Google Scholar
L. Gautier, L. Cope, B. M. Bolstad, and R. A. Irizarry, “Affy–analysis of Affymetrix GeneChip data at the probe level,” Bioinformatics, vol. 20, no. 3, pp. 307–315, 2004.
View at: Publisher Site | Google Scholar
J. T. Leek, W. E. Johnson, H. S. Parker, A. E. Jaffe, and J. D. Storey, “The Sva package for removing batch effects and other unwanted variation in high-throughput experiments,” Bioinformatics, vol. 28, no. 6, pp. 882-883, 2012.
View at: Publisher Site | Google Scholar
G. Yu, L. G. Wang, Y. Han, and Q. Y. He, “clusterProfiler: an R package for comparing biological themes among gene clusters,” OMICS, vol. 16, no. 5, pp. 284–287, 2012.
View at: Publisher Site | Google Scholar
D. Szklarczyk, A. L. Gable, D. Lyon et al., “STRING v11: protein-protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets,” Nucleic Acids Research, vol. 47, no. D1, pp. D607–D613, 2019.
View at: Publisher Site | Google Scholar
P. Shannon, A. Markiel, O. Ozier et al., “Cytoscape: a software environment for integrated models of biomolecular interaction networks,” Genome Research, vol. 13, no. 11, pp. 2498–2504, 2003.
View at: Publisher Site | Google Scholar
H. Han, J. W. Cho, S. Lee et al., “TRRUST v2: an expanded reference database of human and mouse transcriptional regulatory interactions,” Nucleic Acids Research, vol. 46, no. D1, pp. D380–D386, 2018.
View at: Publisher Site | Google Scholar
J. B. Fontanarosa and Y. Dai, “Using LASSO regression to detect predictive aggregate effects in genetic studies,” BMC Proceedings, vol. 5, Suppl 9, p. S69, 2011.
View at: Publisher Site | Google Scholar
S. I. Vrieze, “Model selection and psychological theory: a discussion of the differences between the Akaike information criterion (AIC) and the Bayesian information criterion (BIC),” Psychological Methods, vol. 17, no. 2, pp. 228–243, 2012.
View at: Publisher Site | Google Scholar
P. J. Heagerty, T. Lumley, and M. S. Pepe, “Time-dependent ROC curves for censored survival data and a diagnostic marker,” Biometrics, vol. 56, no. 2, pp. 337–344, 2000.
View at: Publisher Site | Google Scholar
M. E. Ritchie, B. Phipson, D. Wu et al., “Limma powers differential expression analyses for RNA-sequencing and microarray studies,” Nucleic Acids Research, vol. 43, no. 7, article e47, 2015.
View at: Publisher Site | Google Scholar
T. Metsalu and J. Vilo, “ClustVis: a web tool for visualizing clustering of multivariate data using principal component analysis and Heatmap,” Nucleic Acids Research, vol. 43, no. W1, pp. W566–W570, 2015.
View at: Publisher Site | Google Scholar
S. Hänzelmann, R. Castelo, and J. Guinney, “GSVA: gene set variation analysis for microarray and RNA-seq data,” BMC Bioinformatics, vol. 14, no. 1, p. 7, 2013.
View at: Publisher Site | Google Scholar
Z. Zhang and M. W. Kattan, “Drawing nomograms with R: applications to categorical outcome and survival data,” Annals of Translational Medicine, vol. 5, no. 10, p. 211, 2017.
View at: Publisher Site | Google Scholar
A. C. Alba, T. Agoritsas, M. Walsh et al., “Discrimination and calibration of clinical prediction models,” Journal of the American Medical Association, vol. 318, no. 14, pp. 1377–1384, 2017.
View at: Publisher Site | Google Scholar
A. M. Newman, C. L. Liu, M. R. Green et al., “Robust enumeration of cell subsets from tissue expression profiles,” Nature Methods, vol. 12, no. 5, pp. 453–457, 2015.
View at: Publisher Site | Google Scholar
S. A. Minnie and G. R. Hill, “Immunotherapy of multiple myeloma,” The Journal of Clinical Investigation, vol. 130, no. 4, pp. 1565–1575, 2020.
View at: Publisher Site | Google Scholar
D. S. Wishart, Y. D. Feunang, A. C. Guo et al., “DrugBank5.0: a major update to the DrugBank database for 2018,” Nucleic Acids Research, vol. 46, no. D1, pp. D1074–D1082, 2018.
View at: Publisher Site | Google Scholar
D. Kwon, S. M. Kim, and M. A. Correia, “Cytochrome P450 endoplasmic reticulum-associated degradation (ERAD): therapeutic and pathophysiological implications,” Acta Pharmaceutica Sinica B, vol. 10, no. 1, pp. 42–60, 2020.
View at: Publisher Site | Google Scholar
S. Zanwar, B. Nandakumar, and S. Kumar, “Immune-based therapies in the management of multiple myeloma,” Blood Cancer Journal, vol. 10, no. 8, p. 84, 2020.
View at: Publisher Site | Google Scholar
A. Perrot, V. Lauwers-Cances, J. Corre et al., “Minimal residual disease negativity using deep sequencing is a major prognostic factor in multiple myeloma,” Blood, vol. 132, no. 23, pp. 2456–2464, 2018.
View at: Publisher Site | Google Scholar
J. Wang, D. Liu, Z. Sun et al., “Autophagy augments the self-renewal of lung cancer stem cells by the degradation of ubiquitinated p53,” Cell Death & Disease, vol. 12, no. 1, p. 98, 2021.
View at: Publisher Site | Google Scholar
Z. Yun, J. Zhichao, Y. Hao et al., “Targeting autophagy in multiple myeloma,” Leukemia Research, vol. 59, pp. 97–104, 2017.
View at: Publisher Site | Google Scholar
N. Man, Y. Tan, X. J. Sun et al., “Caspase-3 controls AML1-ETO-driven leukemogenesis via autophagy modulation in a ULK1-dependent manner,” Blood, vol. 129, no. 20, pp. 2782–2792, 2017.
View at: Publisher Site | Google Scholar
Y. Shi, E. Norberg, and H. Vakifahmetoglu-Norberg, “Mutant p53 as a regulator and target of autophagy,” Front Oncologia, vol. 10, article 607149, 2020.
View at: Publisher Site | Google Scholar
P. Janovská, E. Normant, H. Miskin, and V. Bryja, “Targeting casein kinase 1 (CK1) in hematological cancers,” International Journal of Molecular Sciences, vol. 21, no. 23, p. 9026, 2020.
View at: Publisher Site | Google Scholar
R. M. Usman, F. Razzaq, A. Akbar et al., “Role and mechanism of autophagy-regulating factors in tumorigenesis and drug resistance,” Asia-Pacific Journal of Clinical Oncology, vol. 17, no. 3, pp. 193–208, 2021.
View at: Publisher Site | Google Scholar
Y. Xu, R. Li, X. Li et al., “An autophagy-related gene signature associated with clinical prognosis and immune microenvironment in gliomas,” Frontiers in Oncology, vol. 10, article 571189, 2020.
View at: Publisher Site | Google Scholar
R. V. Lenchine, S. R. Rao, X. Wang, D. M. Fang, and C. G. Proud, “Elongation factor eEF2 kinase and autophagy jointly promote survival of cancer cells,” The Biochemical Journal, vol. 478, no. 8, pp. 1547–1569, 2021.
View at: Publisher Site | Google Scholar
S. Xu, P. Wang, J. Zhang et al., “Ai-lncRNA EGOT enhancing autophagy sensitizes paclitaxel cytotoxicity via upregulation of ITPR1 expression by RNA-RNA and RNA-protein interactions in human cancer,” Molecular Cancer, vol. 18, no. 1, p. 89, 2019, Published 2019 Apr 18.
View at: Publisher Site | Google Scholar
J. Guo, M. Chen, G. Ai, W. Mao, H. Li, and J. Zhou, “Hsa_circ_0023404 enhances cervical cancer metastasis and chemoresistance through VEGFA and autophagy signaling by sponging miR-5047,” Biomedicine & Pharmacotherapy, vol. 115, article 108957, 2019.
View at: Publisher Site | Google Scholar
Y. Wu, Y. Zhang, X. Qin, H. Geng, D. Zuo, and Q. Zhao, “PI3K/AKT/mTOR pathway-related long non-coding RNAs: roles and mechanisms in hepatocellular carcinoma,” Pharmacological Research, vol. 160, article 105195, 2020.
View at: Publisher Site | Google Scholar
H. Qiu, Y. Li, S. Cheng, J. Li, C. He, and J. Li, “A Prognostic Microenvironment-Related Immune Signature via ESTIMATE (PROMISE model) predicts overall survival of patients with glioma,” Frontiers in Oncology, vol. 10, article 580263, 2020.
View at: Publisher Site | Google Scholar
M. Zamani-Ahmadmahmudi, S. M. Nassiri, and F. Soltaninezhad, “Development of an RNA sequencing-based prognostic gene signature in multiple myeloma,” British Journal of Haematology, vol. 192, no. 2, pp. 310–321, 2021.
View at: Publisher Site | Google Scholar
Y. Zhang, H. Zhang, and B. Zhao, “Hippo signaling in the immune system,” Trends in Biochemical Sciences, vol. 43, no. 2, pp. 77–80, 2018.
View at: Publisher Site | Google Scholar
S. G. Pai, B. A. Carneiro, J. M. Mota et al., “Wnt/beta-catenin pathway: modulating anticancer immune response,” Journal of Hematology & Oncology, vol. 10, no. 1, p. 101, 2017, Published 2017 May 5.
View at: Publisher Site | Google Scholar
E. Batlle and J. Massagué, “Transforming growth factor-β signaling in immunity and cancer,” Immunity, vol. 50, no. 4, pp. 924–940, 2019.
View at: Publisher Site | Google Scholar
M. F. Sanmamed, X. Nie, S. S. Desai et al., “A burned-out CD8+ T-cell subset expands in the tumor microenvironment and curbs cancer immunotherapy,” Cancer Discovery, vol. 11, no. 7, pp. 1700–1715, 2021.
View at: Publisher Site | Google Scholar
E. Milan, M. Fabbri, and S. Cenci, “Autophagy in plasma cell ontogeny and malignancy,” Journal of Clinical Immunology, vol. 36, Suppl 1, pp. 18–24, 2016.
View at: Publisher Site | Google Scholar
S. Bagchi, R. Yuan, and E. G. Engleman, “Immune checkpoint inhibitors for the treatment of cancer: clinical impact and mechanisms of response and resistance,” Annual Review of Pathology, vol. 16, no. 1, pp. 223–249, 2021.
View at: Publisher Site | Google Scholar
J. Bae, F. Accardi, T. Hideshima et al., “Targeting LAG3/GAL-3 to overcome immunosuppression and enhance anti-tumor immune responses in multiple myeloma,” Leukemia, vol. 36, no. 1, pp. 138–154, 2022.
View at: Publisher Site | Google Scholar
H. Handa, Y. Kuroda, K. Kimura et al., “Long non-coding RNA MALAT1 is an inducible stress response gene associated with extramedullary spread and poor prognosis of multiple myeloma,” British Journal of Haematology, vol. 179, no. 3, pp. 449–460, 2017.
View at: Publisher Site | Google Scholar
J. Cavenagh, H. Oakervee, P. Baetiong-Caguioa et al., “A phase I/II study of KW-2478, an Hsp90 inhibitor, in combination with bortezomib in patients with relapsed/refractory multiple myeloma,” British Journal of Cancer, vol. 117, no. 9, pp. 1295–1302, 2017.
View at: Publisher Site | Google Scholar

Copyright

Copyright © 2022 Licheng Li et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

PDF Download Citation

Download other formats

Order printed copies

Views

321

Downloads

527

Citations