Abstract Multiple myeloma (MM) is a distinguished hematologic malignancy, with existing studies elucidating its interaction with neutrophil extracellular traps (NETs), which may potentially facilitate tumor growth. However, systematic investigations into the role of NETs in MM remain limited. Utilizing the single-cell dataset [36]GSE223060, we discerned active NET cell subgroups, namely neutrophils, monocytes, and macrophages. A transcriptional trajectory was subsequently constructed to comprehend the progression of MM. Following this, an analysis of cellular communication in MM was conducted with a particular emphasis on neutrophils, revealing an augmentation in interactions albeit with diminished strength, alongside abnormal communication links between neutrophils and NK cells within MM samples. Through the intersection of differentially expressed genes (DEGs) between NET active/inactive cells and MM versus healthy samples, a total of 316 genes were identified. This led to the development of a 13-gene risk model for prognostic prediction based on overall survival, utilizing transcriptomics dataset [37]GSE136337. The high-risk group manifested altered immune infiltration and heightened sensitivity to chemotherapy. A constructed nomogram for predicting survival probabilities demonstrated encouraging AUCs for 1, 3, and 5-year survival predictions. Collectively, our findings unveil a novel NET-related prognostic signature for MM, thereby providing a potential avenue for therapeutic exploration. Subject terms: Cancer, Haematological cancer, Tumour biomarkers Introduction Multiple myeloma (MM) is the second most prevalent hematologic malignancy derived from plasma cells. Constituting 1–2% of all cancers, MM affects an estimated 34,920 individuals in the US and approximately 588,161 globally each year^[38]1. It is typified by an augmented count of plasma cells in the bone marrow (BM) and increased concentrations of monoclonal immunoglobulins (M-protein) in the serum. These changes lead to complications such as destructive bone lesions, renal impairment, anemia, and hypercalcemia. Notably, the initial symptoms of MM can be ambiguous and resemble other diseases, leading to potential diagnostic and therapeutic delays. For instance, many MM patients experience insidious disease progression with only mildly elevated M-protein concentrations and minimal bone changes, which is easy to overlook. Data from the Surveillance, Epidemiology, and End Results (SEER) program between 2010 and 2016 indicated a 5-year relative survival rate for MM of 53.9%^[39]2. To date, there have been different methods for MM staging and grading. The Durie-Salmon staging system, serving as one of the pioneering methods for MM staging, can effectively gauge the tumor burden but presents a certain limitation in the prognostic evaluation^[40]3. Subsequently, the International Staging System (ISS) was introduced for the preliminary risk stratification of MM. However, ISS lacks adequate considerations about cytogenetics, which play a crucial role in determining the disease's aggressiveness and response to therapy^[41]4. On the basis of the ISS, the Revised-ISS takes the cytogenetics and lactic dehydrogenase (LDH) levels into consideration^[42]5. Nonetheless, it addressed only a limited number of cytogenetic abnormalities with high reproducibility, neglecting the prognostic implications of several core genetic abnormalities and their associated phenotypes in the MM microenvironment. Hence, investigation into novel potential biomarkers is essential and meaningful for the prognostic improvement and therapeutic guidance of MM patients. The BM tumor microenvironment (TME) is pivotal in the pathogenesis and progression of MM. The intricate interaction between MM cells and the BM microenvironment underpins MM cell survival, proliferation, and drug resistance^[43]6. Neutrophils, the predominant cell population within the BM, engage in numerous interactions with myeloma cells. Research has suggested that neutrophils secrete an array of growth factors and cytokines, such as VEGF, TGF-β, and IL-6 to stimulate the growth and proliferation of myeloma cells^[44]7. Meanwhile, different enzymes altering the matrix composition can be released by neutrophils to enhance tumor cell migration^[45]8. Additionally, neutrophils foster angiogenesis^[46]9, facilitating nutrient supply for MM. In turn, MM cells cultivate a profoundly immunosuppressive BM microenvironment, within which several components including myeloid-derived suppressor cells (MDSCs) and N2 neutrophils are amplified^[47]10. Primarily found to neutralize harmful microorganisms, neutrophil extracellular traps (NETs) are extruded from dying neutrophils and present as web-like structures consisting of decondensed DNA chromatin scaffolds and assembled cytosolic and granule proteins. This decondensation of DNA chromatin occurs via citrullination, after which it is expelled from the cell in conjunction with citrullinated histones and neutrophilic cytoplasmic contents rich in granular enzymes–a process termed 'NETosis'^[48]11. NETs have been reported to present effects of a double-sided nature depending on the immune status and interaction with the TME^[49]12. From an antitumor immunity standpoint, NETs impede tumor growth by stimulating the immune system: they facilitate neutrophil interactions with T cells, thereby lowering the activation threshold and directly activating T cells. In contrast, NETs can offer a microenvironment suitable for the delivery of protumorigenic proteins to tumor cells^[50]13. Meanwhile, several studies have indicated that NETs promote the growth and development of tumors via the enhancement of mitochondrial function and induce the activation of corresponding signaling pathways^[51]12,[52]13. In hematological neoplasms, there are extremely intricate interactions between tumor cells and the immune system, making the formation of NETs a more universal and complex phenomenon. The importance of NETs has been described in a number of hematological malignancies, including their impact on various aspects of tumorigenesis^[53]14,[54]15, progression^[55]16, susceptibility to and severity of infection^[56]17, and thrombosis^[57]18. The immunomodulatory role of NETs in leukemia is also likely to be positive, with significant reductions in NETs found during childhood acute lymphoblastic leukemia (ALL) treatment, and increased NETs release as they recover from the disease^[58]19. Immature granulocytes usually persist in the blood of patients during treatment for hematological malignancies and do not release chromatin for use by NETs after activation^[59]20, which may partially responsible for the immunodeficiency. The ability of different hematologic tumors to form NETs appears to vary^[60]15,[61]19, as does their impact on disease pathogenesis. The significance of NETs in MM is not fully understood. It is known that myeloma cells can induce citrullination of histone H3 and prompt NET formation in neutrophils. Elevated concentrations of NETs or their components have been documented in MM patients and may correlate with disease severity and progression^[62]21. Citrullination is of great necessity for DNA chromatin decondensation which is one step of NET formation. In the process of citrullination, PAD4 is a key enzyme. One study showed that MM mice treated with BMS-P5, a specific PAD4 inhibitor, presented a noticeable delay in symptom onset and disease progression^[63]22. This underscores the potential significance of NETs in the pathogenesis and development of MM, suggesting that they could serve as a viable prognostic marker and therapeutic target for the disease^[64]23. The advent of single-cell RNA-sequencing (scRNA-seq) technology, coupled with advances in data analysis techniques, offers an unparalleled window into the molecular features of diverse immune cell populations within the TME^[65]24. Prior research suggests that harnessing gene expression signatures, grounded in the molecular attributes of immune cells extracted from scRNA-seq data, may robustly forecast the prognosis and immunotherapeutic responses of cancer patients^[66]25,[67]26. In the process of literature retrieval, there is not only no score model based on NET-related genes for MM prognostic evaluation but also no TME assessment and therapeutic guidance on the basis of these genes. Our study examined both single-cell and bulk RNA sequencing data from myeloma samples to pinpoint NET activity associated genes in MM. Leveraging a systematic analysis of these genes, we formulated a risk score model aimed at prognostic prediction for MM patients. Meanwhile, our findings further affirmed the model's stability and its effectiveness in predicting patient prognosis and provided a possible potential direction for MM therapy. Results The study's flow chart is illustrated in Fig. [68]1. Figure 1. [69]Figure 1 [70]Open in a new tab The workflow of the study. Single-cell sequencing analysis To identify the origins of highly expressed genes, we scrutinized the cell population of MM using the single-cell sequencing dataset [71]GSE223060. After quality control and removal of doublets, we derived single-cell transcriptomes from 166,757 cells. Figure [72]2A demonstrates that among the 60 samples included, the cell distribution was consistent, suggesting an absence of marked batch effects. This uniformity validates the data for subsequent analyses. Cells were then categorized into 22 distinct clusters, as depicted in Fig. [73]2B. Each cluster's genetic characteristics enabled us to annotate different cell types using cell type-specific biomarkers. Figure [74]2C reveals 11 distinct cell types, including T cells, plasma cells, monocytes, NK cells, B cells, neutrophils, macrophages, dendritic cells (DCs), MAST cells, platelets, and plasmacytoid dendritic cells (pDCs). The dot plots in Fig. [75]2E illustrate specific genes for each cell type, while Fig. [76]2D displays the proportions of these cell types across samples. Figure 2. [77]Figure 2 [78]Open in a new tab Identification of cell subgroups and expression of marker genes from scRNA-sequencing database. (A) UMAP map shows the distribution of MM and control group. (B) UMAP map shows the distribution of MM cell subgroups. (C) UMAP map shows annotation results of MM cell subgroups. (D) Cumulative histogram shows the distribution of cell types in patients with MM and control group. (E) Expression profiles of the marker genes in each cell type. Identification of neutrophil extracellular trap active cells We examined the expression patterns of NET-related genes at the single-cell level within active cell subgroups. Utilizing the optimal threshold to ascertain cell activity, we identified 11,259 cells exhibiting NET activity. Cell clusters with AUC values exceeding 0.17 displayed high NET activity, whereas those with AUC values below 0.17 exhibited low NET activity, as depicted in Fig. [79]3A. Figure [80]3B presents the UMAP diagram of these active cells, indicating that neutrophils, monocytes, and macrophages were predominantly active. Figure 3. [81]Figure 3 [82]Open in a new tab Identification of active cell subgroups. (A) AUC score of the NET-related marker genes, the threshold value was 0.17. (B) UMAP colorogram shows the score of cell activity. The brighter the color, the higher the activity. Pseudo-time trajectories analysis Utilizing the definitive NET-activated subgroups, we constructed a transcriptional trajectory to identify key gene expression programs governing MM progression. The trajectory's transcriptional states highlighted distinct paths, as depicted in Fig. [83]4A and B. Figure 4. [84]Figure 4 [85]Open in a new tab Transcriptional trajectory analysis revealed transcriptional patterns in NET activated subgroups. (A) The pseudo-time color gradient transitions from dark to light blue. (B) The pseudo-time trajectory is divided into three different states by Monocle 2. (C) The DEGs of different branches (different cell fates) shown in heatmap. The top GO BP pathways of different clusters in heatmap were listed nearby. To decipher the molecular underpinnings of this transformation, we investigated the genes influencing MM cell fate. Genes predominantly expressed in the pre-branch were chiefly associated with 'cell killing' and 'leukocyte activation involved in immune response' GO BP pathways. Meanwhile, genes related to 'regulation of hemopoiesis', 'response to virus', and 'regulation of leukocyte cell–cell adhesion' were predominant in cell fate 2. Conversely, cell fate 1 exhibited high expression of genes linked to 'cytoplasmic translation', 'energy-coupled proton transport', and 'ATP synthesis coupled proton transport', as illustrated in Fig. [86]4C and Table [87]S1. Cellular communication patterns in the MM microenvironment To delve deeper into the cellular interaction network within the MM microenvironment, we employed the 'CellChat' R package to discern variations in cell-to-cell communication between the MM and control groups. Comparison to normal tissues revealed an increase in the quantity of interactions among MM samples, accompanied by a decrease in the intensity of these interactions, as illustrated in Fig. [88]5A. Furthermore, Fig. [89]5B illustrates that, in the case of the majority of cell interactions, both the quantity and strength of these interactions exhibited an augmentation when contrasted with normal tissues. These findings underscore the intricate nature of the microenvironment. The outgoing and incoming signaling patterns for both normal and MM tissues are distinctly illustrated in Fig. [90]5. For instance, the macrophage migration inhibitory factor (MIF) signal targeting DCs was diminished in MM (Fig. [91]5C), and the lymphocyte specific protein tyrosine kinase (LCK) signal originating from T cells was reduced in MM (Fig. [92]5D). Figure 5. [93]Figure 5 [94]Open in a new tab Overall pattern of intercellular communication analysis. (A) Bar plot shows the interaction number and strength between MM and normal. (B) The network diagram displays the number and strength of interactions between cell types in the MM and control groups. The red bands represent an increase or enhancement in the number and strength of interactions, while the blue bands represent a decrease or weakening in the number and strength of interactions. (C) Heatmap depicting signals contributing the most to the outgoing signaling pathways in MM and normal. (D) Heatmap depicting signals contributing the most to the incoming signaling pathways in MM and normal. Additionally, we examined receptor ligands potentially mediating communication between neutrophils and other immune cells. Notably, neutrophils communicated with NK cells through the HLA-E-KLRC1 and HLA-E-CD94:NKG2A pathways, which were absent in normal tissues. This implies a role for neutrophil-derived HLA-E in MM progression, as depicted in Fig. [95]6A and B. Figure 6. [96]Figure 6 [97]Open in a new tab Comparison significant ligand-receptor pairs between neutrophil and other cells. (A) Significantly increased ligand receptor pairs in the MM group. (B) Significantly reduced ligand receptor pairs in the MM group. Subsequently, we delved into the expression of MIF and MHC-I pathway genes across different cells in both normal and MM tissues. In comparison to control tissues, the expression of the MIF ligand in neutrophils was notably reduced in MM, and the expression of its receptors in NK cells was also decreased (Fig. [98]7A and B). In contrast, the expression of the MHC-I ligand in neutrophils remained unchanged in MM, and a similar stability was observed in the expression of its receptors in NK cells (Fig. [99]7C and D). This provides insight into the diminished communication intensity of the MIF pathway between neutrophils and NK cells in MM tissues. Figure 7. [100]Figure 7 [101]Open in a new tab The expression of MIF and MHC-I signaling pathways as ligand receptors in tissues. (A) The expression distribution of MIF signaling ligand receptors in the control group. (B) The expression distribution of MIF signaling ligand receptors in the MM group. (C) The expression distribution of MHC-I signaling ligand receptors in the control group. (D) The expression distribution of MHC-I signaling ligand receptors in the MM group. Enrichment analysis of differentially expressed genes related to neutrophil extracellular traps in MM A total of 1,806 DEGs were discerned between NET active and inactive cells (Table [102]S2), exhibiting significant differences (adjusted p value < 0.05; | Log2-fold change |> 0.25). The heatmap in Fig. [103]8A displays the top 10 upregulated (CTSS, S100A9, S100A12, S100A8, MNDA, VCAN, FCN1, LYZ, RP11-1143G9.4, CST3) and downregulated genes (IGHG3, IGKC, IGHG1, JCHAIN, IGHA1, IGLC3, IGLC2, IGHG4, IL32, CCL5). Figure 8. [104]Figure 8 [105]Open in a new tab Enrichment analysis of DEGs related to NET activity in MM. (A) The heatmap shows the significantly DEGs in NET active cells of MM. (B) The heatmap shows the significantly DEGs between MM and controls in single-cell dataset. (C) The Venn diagram highlights the key genes. (D) The circle diagram shows GO and KEGG enrichment results of intersection genes. Differential expression analysis between MM samples and healthy controls was conducted separately for single-cell and bulk transcriptome datasets. From the single-cell dataset, 471 significant DEGs were identified (adjusted p value < 0.05; | Log2-fold change |> 0.25), as depicted in Table [106]S3. Figure [107]8B's heatmap showcases the top 10 upregulated (IGKC, IGHG3, IGLC3, IGHG1, IGHA1, JCHAIN, IGLC2, IGHG4, IGHG2, IGHGP) and downregulated genes (CD69, IGHM, TCL1A, CXCR4, IGHD, IGHD, LEPROTL1, RPS3A, CD79A, RPS12). An overlap between the two DEG sets revealed 316 intersection DEGs, as visualized in Fig. [108]8C and detailed in Table [109]S4. To elucidate the biological functions of these intersection DEGs, we undertook enrichment analysis for GO terms and KEGG pathways. GO analysis, detailed in Table [110]S5, indicated enrichment in biological processes such as cytoplasmic translation, ribosome assembly, and ribosome biogenesis. CC were dominated by features such as cytosolic ribosome and ribosomal subunit, while MF featured rRNA binding and antigen binding (Fig. [111]8D). Prominent KEGG pathways (Table [112]S6) included Coronavirus disease–COVID-19, Ribosome, and Protein processing in the endoplasmic reticulum (Fig. [113]8D). Construction and verification of prognostic risk model Through univariate Cox analysis of the 316 intersection DEGs (between 1806 DEGs form NET active vs. inactive cells and 471 DEGs from MM samples vs. healthy controls), 28 NET-related prognostic genes significantly associated with MM prognosis were identified (p < 0.05) (Table [114]S7). Using random sampling, 70% (n = 291) of the MM samples (n = 420) were allocated to the training set, and the remaining 30% (n = 129) formed the validation set. Genes that have no or little effect on the effectiveness of the predictive model are defined as redundant genes. Redundant genes within the training set were pruned via LASSO regression analysis, with the seed parameter set at 44. This yielded 13 NET-related hub genes significantly tied to MM patient prognosis (Table [115]S8), as illustrated in Fig. [116]9A and B, which were utilized to construct prognostic models. The other 15 genes out of the 28 genes were used as redundant genes and excluded from the model. Figure 9. [117]Figure 9 [118]Open in a new tab Cox and LASSO regression analysis of the MM dataset. (A) Change trajectory of LASSO regression independent variable, the abscissa represents the logarithm of the independent variable λ, and the ordinate represents the coefficient of the independent variable. (B) Confidence interval under each lambda in LASSO regression. (C) The survival curve of patients in high- and low-risk groups from training cohort, respectively. (D) The survival curve of patients in low- and high-risk groups from validation cohort, respectively. Red represents the high-risk group, and blue represents the low-risk group. (E) 1-, 3-, and 5-year time-dependent ROC curves of models for training cohorts. (F) 1-, 3-, and 5-year time-dependent ROC curves of models for validation cohorts. To validate the models crafted from these 13 gene signatures, samples were categorized into low- and high-risk cohorts, using the median risk value as a delimiter. Kaplan–Meier survival curves were plotted for both the training (Fig. [119]9C) and validation cohorts (Fig. [120]9D). Evidently, patients in the high-risk cohort exhibited notably poorer prognoses than their low-risk counterparts across both datasets. To gauge the predictive capacity of the model, ROC curves were generated (Fig. [121]9E and F). In the training set, the 1-, 3-, and 5-year survival AUC values were 0.762, 0.771, and 0.754, respectively (Fig. [122]9E). Correspondingly, in the validation set, these AUC values were 0.570, 0.639, and 0.727 (Fig. [123]9F). GSEA and GSVA To elucidate the potential mechanisms of the DEGs, we employed GSEA. Utilizing the MSigDB Collection, we identified the most significantly enriched signaling pathways according to their normalized enrichment score (NES) (Table [124]S9). Significantly enriched pathways in MM included DNA REPLICATION (NES = 2.1177, adjusted p = 0.0124, FDR = 0.0095, Fig. [125]10A), PARKINSON'S DISEASE (NES = 2.0419, adjusted p = 0.0124, FDR = 0.0095, Fig. [126]10B), SPLICEOSOME (NES = 2.0237, adjusted p = 0.0124, FDR = 0.0095, Fig. [127]10C), P53 SIGNALING PATHWAY (NES = 1.5615, adjusted p = 0.0291, FDR = 0.0223, Fig. [128]10D), RIBOSOME (NES = 1.5339, adjusted p = 0.0441, FDR = 0.0337, Fig. [129]10E), and ASTHMA (NES = -2.1028, adjusted p = 0.0396, FDR = 0.0303, Fig. [130]10F). Additionally, using the MSigDB Collection for GSVA, we highlighted the top 5 pathways with the most significant differential expression between low- and high-risk groups. These findings are visualized in the pathway activity heatmap (Fig. [131]10G and Table [132]S10). Figure 10. [133]Figure 10 [134]Open in a new tab GSEA and GSVA of significantly enriched pathways. DNA REPLICATION (A), PARKINSONS DISEASE (B), SPLICEOSOME (C), P53 SIGNALING PATHWAY (D), RIBOSOME (E), ASTHMA (F). (G) GSVA of significantly enriched pathways. Immune infiltration analysis We examined the infiltration levels of 28 immune cell types between the low- and high-risk groups using the ssGSEA method. The infiltration levels of different immune cells, including activated CD4^+ T cells, activated CD8^+ T cells, effector memory CD4^+ T cells, gamma delta T cells, macrophages, memory B cells, NK cells, NK T cells, regulatory T cells, and type 2 T helper cells, exhibited significant differences between the two groups (p < 0.05, Fig. [135]11A). While most immune cells exhibited positive correlations with each other, a subset demonstrated negative correlations. Specifically, MDSCs, effector memory CD4^+ T cells, type 1 T helper cells, memory B cells, CD56bright NK cells, pDCs, and NK cells infiltration levels were negatively correlated (Fig. [136]11B). Figure 11. [137]Figure 11 [138]Open in a new tab Distinction of immune infiltrations between the high and low risk groups. (A) Boxplot shows the estimated proportion of immune cells between low- and high-risk groups. (B) Correlation among immune cells. Asterisks represented p value (****p < 0.0001, ***p < 0.001, **p < 0.01, *p < 0.05). Additionally, we observed significant correlations between each hub gene and its corresponding immune cells (Fig. [139]12A–I). Notably, the genes ATF7IP2 (R = 0.2057, p < 0.001), MGAT4A (R = -0.1837, p < 0.001), and MEl1 (R = -0.1739, p < 0.001) were significantly associated with memory B cells (Fig. [140]12A–C). Genes ATF7IP2 (R = 0.2532, p < 0.001), RNF125 (R = 0.3258, p < 0.001), and C1orf56 (R = 0.217, p < 0.001) had a significant association with type 2 T helper cells (Fig. [141]12D–F), whereas ATF7IP2 (R = 0.2191, p < 0.001), C1orf56 (R = 0.1989, p < 0.001), and CPIP1 (R = 0.1739, p < 0.001) were significantly related to activated CD4^+ T cells (Fig. [142]12G–I). Figure 12. [143]Figure 12 [144]Open in a new tab Correlation between immune cells and genes. Correlation of gene ATF7IP2 (A), MGAT4A (B) and MEl1 (C) with Memory B cell; Correlation of gene ATF7IP2 (D), RNF125 (E) and C1orf56 (F) with Type2 T helper cell; Correlation of gene ATF7IP2 (G), C1orf56 (H) and CPIP1 (I) with Activated CD4^+ T cell. Construction and verification of the nomogram To ascertain the role of the risk score as an independent prognostic factor, we conducted both univariate and multivariate Cox regression analyses considering clinical characteristics such as age, sex, and risk. Our findings affirm that the risk score stands as an independent prognostic risk factor for patients, irrespective of the Cox regression analysis employed (Fig. [145]13A and B). Utilizing multivariate Cox regression analysis, we constructed a nomogram, demonstrating that the risk score can significantly forecast clinical outcomes (Fig. [146]13C). We employed the ROC curve to evaluate the predictive efficacy of the nomogram concerning patient prognosis. The AUC values for 1-, 3-, and 5-year survival were 0.735, 0.756, and 0.770, respectively (Fig. [147]13D). Figure 13. [148]Figure 13 [149]Open in a new tab Risk score is an independent prognostic factor for clinical characteristics. (A) Forest map shows the results of univariate Cox regression analysis performed on clinical characteristics. (B) Forest map shows the results of multivariate Cox regression analysis performed on clinical characteristics. (C) The nomogram of the prediction model. The line segment represents the contribution of the clinical factor to the outcome events, total points represent the total score of the sum of the corresponding individual scores of the value of all variables, and the bottom three lines represent the prognosis of 1-, 3-, and 5-year survival corresponding to each value point. (D) 1-, 3-, and 5-year time-dependent ROC curves of nomogram. Drug susceptibility analysis We investigated the ability of the risk score to predict chemotherapeutic sensitivity in MM patients. We assessed several agents, namely cyclophosphamide_1512, bortezomib_1191, cisplatin_1005, epirubicin_1511, AZD7762_1022, venetoclax_1909, doramapimod_1042, vincristine_1818, and KU-55933_1030, for their therapeutic effectiveness in MM treatments (Fig. [150]14, Table [151]S11). Our findings indicate that patients with a high-risk score might have increased sensitivity to conventional chemotherapy drugs, specifically cyclophosphamide_1512, cisplatin_1005, epirubicin_1511, vincristine_1818, and the proteasome inhibitor (PI) bortezomib_1191. This suggests that a combination of PI and multidrug chemotherapy may be optimal for this patient cohort. Figure 14. [152]Figure 14 [153]Open in a new tab Drug susceptibility between the low- and high-risk groups. (A) Difference in sensitivity to Cyclophosphamide_1512 between low- and high-risk groups. (B) Difference in drug sensitivity of Bortezomib_1191 between low- and high-risk groups. (C) Differences in drug susceptibility to Cisplatin_1005 between low- and high-risk groups. (D) Difference in Podophyllotoxin Epirubicin_1511 drug sensitivity between low- and high-risk groups. (E) Differences in AZD7762_1022 susceptibility between low- and high-risk groups. (F) Differences in drug sensitivity of Venetoclax_1909 between the low- and high-risk groups. (G) Differences in susceptibility to Doramapimod_1042 between low- and high-risk groups. (H) Difference in Vincristine_1818 susceptibility between low- and high-risk groups. (I) Difference in drug sensitivity to KU-55933_1030 between low- and high-risk groups. Discussion With immunotherapy advancing and gaining more attention, an increasing number of biomarkers have been explored for immunotherapy response prediction^[154]27. Due to the nonnegligible and vital effects of the TME on tumor occurrence and development, influence of the TME on cancer immunotherapeutic efficacy has been extensively examined, and TME-related biomarkers have attracted heightened attention^[155]28. However, there remains a paucity of reliable biomarkers and risk score models reflecting the roles of the tumorigenic TME in immunotherapeutic responses and prognosis in MM. The evolution of scRNA-seq technology offers a comprehensive lens into the molecular profiles of tumor-infiltrating immune cells within the TME^[156]29. In this study, 41 MM samples were analyzed via single-cell sequencing, and 11 distinct cell types were identified. We observed a higher proportion of neutrophils in myeloma samples than in normal controls, unlike lymphocytes which were suppressed by elevated tumorigenic plasma cells (Fig. [157]1D). The ability of neutrophils in MM to form NETs also appears to be like other hematological tumors and does not directly correlate with their count^[158]19, which needs to be further determined. We utilized NET-related genes sourced from the work of Zhang et al.^[159]30, for the subsequent evaluation of NET activity via the AUCell algorithm. Cell populations with AUC values exceeding 0.17 were characterized by high NET activity, suggesting that NETs might modulate these cells, influencing tumorigenesis and progression. Subsequently, the pivotal genes governing NET activity were further pinpointed by combining the results of the difference and enrichment analyses. The 13-gene signature effectively differentiates patients into low- and high-risk subpopulations. Across all training and validation sets, our signature demonstrated great consistency and stability. Specifically: (1) patients from distinct risk subpopulations were clearly distinguishable; (2) high-risk subpopulation patients exhibited a poor prognosis; (3) tumor immune microenvironments between low- and high-risk subpopulations exhibited significant differences; and (4) the signature's diagnostic values for 1-year, 3-year, and 5-year survival rates were commendable. We found that the 1, 3, and 5-year AUC values of our NET-related prognostic model were higher than of the UAMs GEP70^[160]31 model in the same dataset ([161]GSE136337), which suggested its good prognostic predictive accuracy (Supplementary Fig. [162]1A, B). Interestingly, we chose to externally validate the NET-related prognostic model with a dataset of screened myeloma cells ([163]GSE4581), which still resulted in a significant difference in prognosis between the differentiated low and high-risk groups (Supplementary Fig. [164]2A, B). Among the 13 signature-associated genes (RNF125, NPM1, CRIP1, HIST1H1C, C1orf56, S100A6, GAPDH, CCND1, RHOH, ANKRD28, MEI1, MGAT4A and ATF7IP2), RNF125, NPM1, CRIP1, HIST1H1C, C1orf56, S100A6 and GAPDH were risk genes, while CCND1, RHOH, ANKRD28, MEI1, MGAT4A and ATF7IP2 served as protective genes. Among those risk genes, RNF125 was significantly linked to a high score and poor prognosis. RNF125, an E3 ubiquitin ligase, is involved in tagging specific proteins, leading to their ubiquitination and subsequent degradation. In immune processes, RNF125 ubiquitinates key signaling molecules, influencing their stability and function. This mechanism is essential in MM and is regarded as the target of PIs such as bortezomib^[165]32. Additionally, RNF125 may function as a positive regulator in the T-cell receptor signaling pathway, potentially affecting T-cell infiltration^[166]33. Consistent with this conclusion, our study found that the high-risk cohort presented an increased infiltration of activated CD4^+/CD8^+ T cell. Nonetheless, the prognosis was unfavorable in this group, potentially due to T-cell anergy via mechanisms involving GRAIL^[167]34. NPM1, implicated in various cellular activities, shows significant overexpression in hyperdiploid MM due to chromosome 5 gains, suggesting its key role in the pathogenesis of hyperdiploid MM^[168]35. Additionally, CRIP might be linked to intestinal zinc transport and myeloma bone disease severity^[169]36. Our findings reflected that C1orf56 served as a risk gene in MM, which is supported by the fact that C1orf56 is a proto-oncogene repressed by DNMT3B methylation. S100A6 expression was notably higher in primary MM patients than in controls, associating it with MM progression and intramedullary metastasis^[170]37. MEI1, presumed to be involved in meiosis I, is linked to gestational trophoblastic neoplasms. In our framework, elevated MEI1 expression contributes favorably to MM prognosis, but the exact mechanism requires further study. MGAT4A, encoding a pivotal glycosyltransferase, was observed to function as a protective gene in MM and positively impact prognosis in our study. The protective effect of MGAT4A has been observed in breast cancer, in which diminished expression is related to drug resistance^[171]38. In a preliminary search of the literature, we found no overlap between the 13 genes in this model and those in existing MM prognostic models. This may be related to the fact that there are currently no prognostic models for MM based on the genes of NETs, as well as our focus on the potentially possible significant prognostic contribution of the microenvironment. However, the mechanism through which the different myeloma cells stimulate the activation of NET-related genes in microenvironmental cells and the occurrence of adverse prognostic expression profiles needs to be further explored. We noted that among the transcriptomic features of different cell types (malignant plasma cells vs. microenvironmental cells), six of the risk genes, C1orf56, CRIP1, GAPDH, HIST1H1C, RNF125, and S100A6, had a consistent trend of expression in the low- and high-risk groups differentiated by the NET-related prognostic model (Supplementary Fig. [172]3A, B). In order to clarify the difference in the contribution of these genes to prognosis in different cell types, we performed Quantitative Real-time PCR to detect the expression of the five genes (GAPDH as an internal reference gene) in myeloma and bone marrow stromal cells. The results (Supplementary Fig. [173]4A–E) suggest that three genes, CRIP1, HIST1H1C and RNF125, are significantly overexpressed in myeloma cell lines and may contribute to the poor prognosis in malignant plasma cell samples, whereas C1orf56 and S100A6 may contribute to the poor prognosis in microenvironmental cell samples. The GO and KEGG annotation results revealed primary enrichments in cytoplasmic translation, ribosome assembly, and ribosome biogenesis (BP). Furthermore, significant annotations were noted in the CC of the cytosolic ribosome, ribosomal subunit, and ribosome. MF predominantly involved the structural constituent of ribosome, rRNA binding, and antigen binding. Other enrichments were identified in pathways such as COVID-19, ribosome, and protein processing in the endoplasmic reticulum. Many of these factors have been previously linked to the pathogenesis of ribosome biogenesis. Notably, the expression of genes related to ribosome biogenesis correlates with disease progression and prognosis in MM patients, suggesting potential therapeutic targets, including BRD9^[174]39. GSEA facilitates the extraction of valuable insights from large-scale gene datasets, even with minimal fold changes. Utilizing GSEA of our gene datasets, we identified numerous gene sets significantly enriched within the MM group. Specifically, DNA replication gene overexpression in differentiated tissues could suggest a pathological state. In the MM dataset, this might imply an accelerated division of myeloma cells, a characteristic often observed in malignant tumors. The spliceosome, a complex molecular entity predominantly found in eukaryotic cell nuclei, primarily functions in excising introns and ligating exons during pre-mRNA processing, a critical step known as RNA splicing. This is essential for the optimal maturation of mRNA prior to its export for protein translation. Disruptions in spliceosome function can induce aberrant mRNA splicing, potentially producing dysfunctional proteins. Research by Hector H. Huang identified splicing interference as a novel aspect of the PI mechanism, unveiled further spliceosome modulation methods, and posited spliceosome targeting as a promising therapeutic approach for MM^[175]40. Given the likelihood that the DEGs we extracted came from cells from the BM microenvironment, common comorbidities of elderly such as Parkinson's disease and asthma may be confounding. We chose the MM dataset with screened myeloma cells ([176]GSE4581) for GSEA and GSVA analyses and obtained different results. Pathway co-enriched by malignant plasma cells and microenvironmental cells was DNA REPLICATION. Also, the difference was that malignant plasma cells were enriched to pathway properties of the tumor cell itself, such as the CELL CYCLE and PROTEASOME, whereas microenvironmental cells were enriched to the pathway of transcription and translation, such as SPLICEOSOME and RIBOSOME (Supplementary Fig. [177]5A–F). GSVA analyses were also suggestive of showing that myeloma cells and microenvironmental cells differ significantly in the pathways with the most significant differential expression in the low- and high-risk groups distinguished by our NET-related gene prognostic model (Supplementary Fig. [178]5G). The dynamic interaction between myeloma cells and the BM microenvironment plays a pivotal role in malignant transformation, treatment response, and disease progression. Our comprehensive investigation of the prognostic-signature-based immune distinctions revealed that the high-risk group exhibited augmented infiltration of cells linked to adaptive immunity. Conversely, the low-risk group had a pronounced infiltration of cells associated with the innate immune system. These variations in cellular infiltration suggest a heightened propensity in the high-risk group to develop immune evasion through mechanisms such as immune resistance, exhaustion, and suppression^[179]41. The high-risk cohort showed increased infiltration of type 2 T helper cells (HR = 5.7, p = 0.013) and decreased infiltration of NK cells (HR = 0.089, p = 0.046). Both have been previously identified as significant adverse prognostic factors in MM^[180]42. Type 1 T helper cells generate IFN-γ, bolstering the cell-mediated immune response, while type 2 T helper cells, which produce IL-4, counteract this type 1 T helper cells response. Research by Faqing Tian et al. uncovered that myeloma cells could serve as antigen presenting cells (APCs), showcasing microbial antigens to type 2 T helper cells, which spurs their proliferation and thereby aids tumor progression via intimate Th2-myeloma cell interactions^[181]43. NK cells possess a spectrum of antitumor and immunomodulatory functions. A direct correlation exists between NK cell activity and disease-free survival in MM patients. Reduced NK cell activity aligns with advanced clinical stages, elevated LDH, heightened BM plasma cell infiltration, and increased β2 microglobulin levels^[182]44. Those MM patients exhibiting long-term disease stability displayed an expansion of NK cells^[183]45. In our study, this finding is consistent with those previous conclusions made by others. Additionally, our findings suggest that neutrophils in myeloma patients engage in aberrant interactions with NK cells via receptors such as HLA-E, which could account for the diminished NK cell presence in the high-risk cohort, consequently affecting prognosis. Therefore, devising immune therapies targeting NK cells, such as BCMA CAR-NK^[184]46, seems promising for the high-risk group. Evidently, these immune cellular infiltration disparities could be the underlying factors for the adverse prognosis observed in the high-risk group, with this group exhibiting markedly lower survival rates than their low-risk counterparts. To elucidate the impact of immune cell infiltration in MM more profoundly, we employed ssGSEA for a thorough assessment of immune infiltration within MM contexts. Our analysis revealed that heightened infiltration of memory B cells, type 2 T helper cells, and activated CD4^+ T cells might be intrinsically linked to the onset and progression of MM. Notably, ATF7IP2 demonstrated a significant association with memory B cells, while ATF7IP2, RNF125, and C1orf56 showed substantial correlations with type 2 T helper cells. Furthermore, ATF7IP2, C1orf56, and CPIP1 exhibited strong associations with activated CD4^+ T cells. Conversely, MGAT4A and MEI1 displayed negative correlations with memory B cells. Within our model, ATF7IP2 emerges as the predominant gene contributing to a low-risk profile. This observation aligns with prior findings in studies on lung cancer^[185]47. Interestingly, ATF7IP2 bears positive correlations with diverse immune cell infiltrations, typically linked to unfavorable prognoses. This implies that ATF7IP2 might positively influence prognosis via alternative pathways. Such hypotheses necessitate further exploration to elucidate the intricate interplay between genes and immune cells, providing a direction for the in-depth exploration and supplementation of the potential molecular mechanism in MM. Drug sensitivity prediction analysis suggests that high-risk patients might have heightened sensitivity to bortezomib treatment. This raises the possibility that primary induction therapy, which includes the PI bortezomib such as the VCD (bortezomib, cyclophosphamide, dexamethasone) regimen, may not fully explain the observed prognostic disparities between groups. Moreover, venetoclax appears to be more effective in the high-risk cohort. Its use has been chiefly confined to patients with cytogenetic abnormalities characterized by the t(11:14) translocation, which results in the IgH/CCND1 fusion gene^[186]48. Notably, our prognostic model includes the CCND1 gene, which is believed to act as a protective factor. This is consistent with the favorable prognosis in the myeloma subgroup with the IGH/CCND1 fusion gene, implying the venetoclax might be beneficial for high-risk patients with elevated CCND1 expression. Additionally, Selinexor, a selective nuclear export inhibitor, has shown significant efficacy against relapsed and refractory MM^[187]49. Given its capacity to inhibit human NET formation in vitro^[188]50, selinexor may offer potential as a salvage therapy for high-risk patients. Furthermore, AZD7762, an ATP-competitive checkpoint kinase inhibitor, augments checkpoint termination and bolsters DNA-focused treatments. For refractory high-risk patients, its combined administration with cytotoxic agents or immune checkpoint inhibitors^[189]51 could be advantageous. In contrast, for low-risk patients, conventional treatments yield moderate drug sensitivity. For those unresponsive to standard treatments, small molecule drugs such as the MAPK inhibitor Doramapimod and the ATM inhibitor KU-55933 might be viable options. The results of drug sensitivity prediction analysis provide a possible direction for MM therapy, but all of speculations need in-depth experiments to verify. This study presents several limitations warranting acknowledgment.