Introduction
The advent of anti-CD38 monoclonal antibodies has transformed the treatment landscape of multiple myeloma (MM). In particular, the combination of daratumumab with the immunomodulatory drug lenalidomide (DR) has demonstrated remarkable efficacy, with 91% of overall response rate, 69% of very good partial response rate or better and a median of progression free survival of 53.4 months, becoming a standard treatment for newly diagnosed elderly frail patients with MM [ 1 – 3 ]. However, despite the substantial success of this regimen, some patients fail to achieve complete remission, and all eventually relapse. To improve patients selection for an optimal response to DR therapy, a deeper understanding of the immune microenvironment necessary is critical.
Daratumumab targets CD38-overexpressing MM cells for destruction through multiple mechanisms including antibody-dependent cell-mediated cytotoxicity (ADCC), antibody-dependent cellular phagocytosis (ADCP), and complement-dependent cytotoxicity (CDC) [ 4 , 5 ]. It has also been shown to deplete CD38 + regulatory T cells, potentially enhancing the activity of helper and cytotoxic T cells [ 6 ]. Lenalidomide acts directly on plasma cells by promoting the degradation of the lymphoid transcription factors IKZF1 and IKZF3, which are crucial for plasma cell differentiation and MM cell progression, via the CRBN-CRL4 ubiquitin ligase complex [ 7 , 8 ]. Lenalidomide also exerts an indirect immunomodulatory effect by enhancing the immune system, including the induction of interleukin (IL)-2 production in T cells and the increased expression of death effector molecules [ 8 , 9 ]. When co-administered to patients with MM, daratumumab and lenalidomide have been found to exhibit synergistic effects via activation of natural killer (NK) cells and expansion of effector memory T cells [ 9 ].
In MM, the bone marrow (BM) immune microenvironment undergoes significant alterations throughout disease progression [ 10 – 12 ]. For instance, an increase in NK cells is often observed in the early stages, accompanied by altered chemokine receptor expression [ 12 ]. The frequency of mature cytotoxic CD56dim NK cells decreases, while the frequency of CD16/CD226low NK cells with adhesion defects and impaired effector functions increases [ 13 ]. At the smoldering MM stage, the loss of granzyme K+ memory cytotoxic T cells and the dysregulation of CD14 + monocytes contribute to tumor progression [ 12 ]. These immune alterations may be driven by inflammatory mesenchymal stromal cells, which co-localize with tumor and immune cells and express genes involved in tumor survival and immune modulation [ 14 ].
The implications of these immune alterations in the context of MM treatment remain poorly understood, and their impact on patient outcomes is unclear. We hypothesized that patient-level differences in the BM immune microenvironment may associated with response, to DR treatment. To test this, we conducted single-cell RNA sequencing (scRNA-seq) of 21 BM samples from patients treated with DR as part of the IFM2017-03 clinical trial ( NCT03993912 ). Our findings aim to investigate whether an enrichment of inflammatory signatures in the BM of poor responders is associated with dysfunctional states of CD56dim NK cells and impaired phagocytic capacity of CD14 + monocytes.
Results
Study population
To assess the impact of the BM immune microenvironment on MM treatment, we performed scRNA-seq on BM samples from 21 patients with newly diagnosed MM (NDMM) (Fig. 1 a and Supp. Table 1 ). All patients were enrolled in the IFM2017-03 clinical trial and treated with the DR combination therapy. This phase 3 study evaluated a dexamethasone-sparing strategy for frail patients with NDMM. The 21 BM patients analyzed by scRNA-seq, were obtained from patients included in arm B of the clinical trial, which received DR until progression or unacceptable toxicity. The main clinical characteristics of the patients are summarized in Supplementary Table 1 . Patients were divided into two groups based on treatment response: 11 complete responders (CR) and 10 non-responders or poor responders (NR).
Inflammatory bone marrow microenvironment at diagnosis is associated with poor responses
A total of 153,913 single cells were sequenced by paired transcriptome and a panel of 35 antibodies, with a median of 6,505 cells (range: 3,646 − 11,360) per sample. BM scRNA-seq data from 20 healthy donors (HD) from Oetjen et al. [ 15 ] were selected as the control group (Supp. Table 2 ). After quality control and normalization, we analyzed a total of 149,365 cells from MM patients and HD, encompassing all the major cell types of the BM: hematopoietic stem and progenitor cells (HSPC), erythroid progenitor cells (EPC), stromal cells, T & NK cells, B cells, myeloid cells, normal plasma cells (nPC, from HD) and MM plasma cells (MMPC, from MM patients) (Fig. 1 b and Supp. Figure 1 a).
Matched gene expression levels and antibody expression level analysis effectively identified key clusters, including CD34 for HSPC, PTPRC and KLRB1 for T & NK cells, CD19 and CD79A for B cells, and FCN1 and ITGAM for myeloid cells. nPC clustered together and exhibited high expression of MZB1 and low expression of SDC1 , while MMPC had high expression of SDC1 (Fig. 1 c-d and Supp. Figure 1 b-c). There was no significant difference in the proportions of major immune cell populations between CR and NR (Fig. 1 e). However, when compared to HD, MM patients had higher proportions of myeloid cells and an enrichment of T and NK cell ratios, consistent with previous findings [ 16 ] (Fig. 1 e). Consistent with these findings, flow cytometry analysis of BM samples from CR and NR patients showed no significant changes in the percentage of B, T, and NK cells, however we found increased myeloid cell proportions in NR patients (Fig. 1 f and Supp. Figure 2 ).
Next, to begin exploring specific factors that may influence patient responses to DR treatment, we compared the characteristics of the different response groups. Heatmap visualization and clustering of cell subset proportions across samples identified four modules (Fig. 1 g). Notably, Module 4 was highly enriched in the NR group and was characterized by enrichment of inflammatory pathways, including hallmark_inflammatory_response and hallmark_TNFa_signaling_via_NFkB pathways (Fig. 1 g-h). Furthermore, we compared three inflammatory, cytokine and interferon (IFN) signature scores, associated with hallmark pathways, among CR and NR groups at the BM microenvironment (BME) level. These signatures were predominantly enriched in myeloid cells and subsets of T & NK cells (Fig. 1 i). A clear distinction emerged, with significantly higher scores for all three signatures (inflammation, cytokine, and IFN) in NR patients compared to CR patients (Fig. 1 j), with a consistent trend of higher combined scores in NR at the patient level (Fig. 1 k. This consistent single-cell and patient-level trend supports the biological relevance of the coordinated elevation in inflammatory, cytokine, and IFN signaling and suggests that an inflamed BME may be correlated with poor response to DR therapy.
Inflammation-associated classical monocytes exhibit lower phagocytosis scores in NR
We next sought to the specific myeloid subsets underlying the inflammatory BME and their capacity to activate NF-κB signaling. We integrated myeloid cells from all 21 patients and 20 HD and annotated them into 13 clusters based on expression levels of marker genes (Fig. 2 a). The largest cluster consisted of classical monocytes (cMono) characterized by high expression of FCN1 and S100A12 , indicative of their origin from blood-derived monocytes. Intermediate monocytes had relatively weak expression of CD14 and CD68 , while classical monocyte/macrophages highly expressed CD14 with the macrophage marker CD163 (Fig. 2 b and Supp. Figure 3 a-b), with both clusters representing intermediate states of differentiation (Supp. Figure 3 c). Nonclassical monocytes exhibited high expression of FCGR3A but lacked CD14 (Fig. 2 b). TOP2A and MKI67 were highly expressed in the proliferating monocyte (prolMono) cluster (Fig. 2 b). Two clusters exhibited ELANE monocyte marker, with one cluster defined as monocyte precursors (preMono, high in RETN expression) and the other as granulocyte-monocyte progenitors (GMP). We also identified four dendritic cell (DC) clusters: plasmacytoid DCs (pDC) exhibited high IL3RA expression, while conventional DCs of types 1, 2 and 3 (cDC1, cDC2, cDC3) expressed high levels of ITGAX and were characterized by specific expression of CLEC9A , CD1C and CD22 , respectively (Fig. 2 b).
Single-cell analysis enhanced heightened inflammatory activity, cytokine signaling and NF-κB pathway activation signatures in classical monocytes (Fig. 2 c), while phagocytosis-related genes and pathways were highly enriched in both classical and nonclassical monocyte subsets in the CR group (Fig. 2 d-e). Given prior evidence that ADCP contributes to daratumumab efficacy [ 4 , 5 ], the significantly higher ADCP-related gene signature score observed in cMono from CR group further suggests that enhanced ADCP-related transcriptional programs may be aligned with favorable therapeutic responses (Fig. 2 f).
Consistent with previous reports describing interferon-responsive monocyte states in MM [ 12 ], we identified an ISG+ cMono subset characterized by high expression of IFN-inducible genes (e.g., IFI6 , ISG15 , and IFIT2 ) and elevated IFN scores (Fig. 2 c and g-h ). Notably, differential expression analysis further revealed that ISG+ cMono cells in the NR group exhibited an inflammatory phenotype, with upregulation of the NF-κB signaling genes FKBP5 and TNFAIP3 and high expression of proinflammatory factors such as CXCR4 , IFITM3 , IFI6 and IRF1 (Fig. 2 i). Pathway analysis of these differentially expressed genes revealed distinct transcriptional programs corresponding to treatment response: phagocytosis-related and T cell s activation pathways were enriched in the CR group, whereas genes involved in the negative regulation of T cell responses and the production of key immunoregulatory factors such as IL12 were enriched in the NR group (Fig. 2 j). Together, these findings indicate that ISG+ cMono cells in the NR group display a more inflammatory and immunoregulatory transcriptional profile than those in the CR group.
Although myeloid cells are significantly increased in NR patients (Fig. 1 f), a deeper flow cytometry analysis of such population did not reveal any significant change in the proportions of cMono and nonclassical monocytes (Supp. Figure 4 ). Interestingly, we observed a consistent increase in several inflammatory markers within the cMono population. We found that cMono from NR patients were characterized by a significantly increased expression of CD27, a member of the TNF-receptor superfamily, that induces NF-κB signalling upon engagement [ 17 ]. Together with an increased, although not significant, expression of HLA-DR, KLF4 and IRF1 [ 18 , 19 ], these data suggest a more inflammatory profile of cMono.
All together, the single-cell transcriptomic and flow cytometric analyses indicate that inflammatory reprogramming of ISG+ cMono in NR patients was characterized by reduced ADCP-related transcriptional programs, altered T cell activation signatures, and enhanced NF-κB-driven inflammatory profiles. This suggests that ISG + cMono states may contribute to the inflammatory features underlying poor response to DR treatment.
Inflammation-driven NK cells exhibit impaired cytotoxicity-related transcriptional programs in NR patients
We next explored the role of NK cells in predicting the response to DR. NK cells can be functionally categorized into cytotoxic CD56dim and cytokine-producing CD56bright (NKbright) subsets. We defined clusters for these subsets based on their transcriptional features (Fig. 3 a). The NKdim was characterized by elevated expression of FCGR3A (encoding CD16), GZMB and PRF1, whereas the NKbright subset exhibited high expression of GZMK and NCAM1 (encoding CD56) (Fig. 3 b). Differential analysis of these two NK clusters revealed that FKBP5, which promotes NF-κB signaling, was significantly enriched in the NR group (Fig. 3 b). In contrast, the CR group exhibited enrichment of granzymes (GZMA and GZMK) which are known to induce cancer cell apoptosis, along with several chemokines (CX3CR1, CCL3, CCL4 and CCL5) (Fig. 3 c). Functional enrichment analysis of the NKdim group indicated reduced ADCC-related gene signatures in NR compared with CR patients (Fig. 3 d-e). Moreover, the enrichment of the cytotoxic pathway in both NKdim and NKbright revealed higher cytotoxicity-related transcriptional programs in CR patients, whereas these signatures were reduced in NR patients, providing transcriptomic evidence for distinct cytotoxicity-associated states of NK cells between the two response groups.
CellChat analysis further predicted potential interactions between MMPC and NK cells. Specifically, MMPC in CR patients were predicted to interact with both NKdim and NKbright cells, whereas NK cells in the NR group showed few predicted ligand-receptor interactions with MMPC ( Fig. 3 f-g). Additionally, we observed significantly lower CD16 expression in NKdim cells at both the gene and protein levels in the NR group (Fig. 3 h). Flow cytometry on matched samples validated this finding (Fig. 3 i) and similar to cMono populations showed a trend in increased HLA-DR expression in NK cells derived from NR patients (Supp. Figure 5 ). These data suggest that reduced CD16 expression on NKdim cells in NR patients may be associated with impaired NK cell-tumor cell interactions and reduced cytotoxicity-related transcriptional profiles.
NR and CR patients display distinct T cell profiles
T cells play a crucial role in adaptive immunity and are essential for the elimination of tumor cells by the immune system [ 20 ]. To dissect the diversity of T cells in MM, we further extracted and clustered T cells into three conventional CD4 + clusters: CD4 naïve T cells (CD4 TN), CD4 central memory T cells (CD4 TCM), and CD4 effector memory T cells (CD4 TEM); and three CD8⁺ subsets: CD8 naïve T cells (CD8 TN), CD8 memory T cells (CD8 TM), and CD8 cytotoxic T cells (CD8 tox). Additional clusters included regulatory T cells (Treg), T helper 17 cells (Th17), ISG + T cells, γδ T cells (gdT), NKT cells, and two NK clusters (NKbright and NKdim) (Fig. 4 a-c and Supp. Figure 6 a-d). We next compared the transcriptional profiles associated with functional activity of these subsets across response groups, reflecting its relevance to antitumor immunity. Consistently, effector T cell populations exhibited higher effector scores in the CR group, with significant increases in CD8 tox cells, NKT cells and both NKbright and NKdim subsets (Fig. 4 d). These results suggest that transcriptional signatures of enhanced effector activity in CR patients may contribute to more effective antitumor immune responses.
We found that the ISG + T cell cluster, defined by expression of TNF , IFIT1 and IFIT2 , had the higher median scores for inflammation, cytokine and IFN signaling, consistent with the overall significant difference by Kruskal-Wallis test (Supp. Figure 6 e). Interestingly, ISG + T cells were more abundant in NR (Fig. 4 e). These cells exhibited high expression of proinflammatory genes (e.g., LTA , CCL3L1 , CCL5 , KLF6 and IFNG ) and were enriched in inflammatory pathways (e.g., IL-1β, IFN-γ and NF-κB) that are key regulators of immune and inflammatory responses (Fig. 4 e-g). Additionally, ISG + T cells in the NR group showed higher expression of inflammatory factors and greater enrichment in inflammation-related pathways (Fig. 4 h-i). These results indicate that distinct T cell functional states are associated with differential responses to DR therapy in MM, with enriched effector T cells signatures observed in CR patients and proinflammatory ISG⁺ T cells contributing to impaired responses in NR patients.
Flow cytometry analysis of MM infiltrating T cells confirms not significant differences in the distribution of CD4 + ad CD8 + T cells (Supp. Figure 7 ), however we found in naïve/central memory (CCR7 + CD27+) CD4 + T cells from NR patients a significant overexpression of HLA-DR (Fig. 4 j), further confirming the presence of an inflammatory environment in NR subjects.
Inflammatory microenvironment correlates with NF-κB activation in MM plasma cells
We performed bulk RNA sequencing (RNA-seq) of MM tumor cells from 11 patients (five CR and six NR) to examine the clinical and biological properties of CR and NR and to assess consistency with our single-cell findings (Fig. 5 a). Principal component analysis (PCA) clearly separated CR and NR patients along the first two principal components (Fig. 5 a). The top 100 differentially expressed genes were identified for both groups, with NR highly expressing NF-κB-related genes including CLEC7A , CCR7 , TRIM8 , CARD9 , SLC44A2 , MAP3K3 , FLOT1 and TLR2 (Fig. 5 b).
We also analyzed the MMPC present in the single-cell sequencing data from 21 patients. To better characterize the tumor cells, we used 694 nPC from HD as a reference for calculating the copy number variation (CNV) scores. All MMPC had higher CNV scores than the healthy references (Supp. Figure 8 a), and the CNV score was significantly higher in NR than in CR patients, indicating a greater level of genomic variants (Supp. Figure 8 a). A total of 12,928 MMPC were identified in the single-cell sequencing data.
To explore tumor cell-specific alterations, we extracted cells with high tumor gene scores for re-clustering, which showed clear patient-specific segregation in uniform manifold approximation and projection (UMAP) space (Supp. Figure 8 b-c). Consistent with the bulk RNA-seq results, the top two principal components from the scRNA-seq data revealed a separation trend between cancer cells from CR and NR patients, reflecting underlying transcriptional differences (Fig. 5 c). The overall variance appeared more related to treatment response than to inter-patient variability. Importantly, when clinical factors such as gender, age and clinical stage were considered, patient distribution along PC1 remained independent of these variables (Fig. 5 d). This observation indicates that the variance observed between CR and NR patients primarily relates to drug response.
We then calculated the genes highly expressed in each response group at the PC1 level, identifying the top 200 genes in NR patients, which were associated with cytokines, immune inflammatory response and the NF-κB signaling pathway (Fig. 5 d). The NF-κB pathway regulates inflammatory genes that influence cancer development and pathogenesis of MM [ 20 ]. Further, differential gene expression analysis at the scRNA-seq level showed enrichment of important transcription factors such as KLF6 and KLF2 in the NR group (Fig. 5 e). Gene ontology (GO) analysis confirmed that MM tumor cells from NR patients were characterized by increased immune inflammatory response, NF-κB signaling and cytokine production (Fig. 5 f). Pseudo-bulk analysis showed consistent results for NR (Supp. Figure 8 d-f).
We also identified a set of 34 genes associated with NF-κB signaling by integrating single-cell and bulk transcriptomic analyses, which were preferentially expressed in NR patients (Supp. Figure 8 g). Among these, 23 genes were significantly associated with poor prognosis after Benjamini–Hochberg false discovery rate (FDR) correction in the MMRF CoMMpass cohort (FDR < 0.05 & hazard ratio [HR] > 1), and 16 genes were further predictive of poor prognosis in 177 cases treated with lenalidomide and/or daratumumab (Supp. Figure 8 h-i). These findings suggest that elevated NF-κB signaling is associated with poor response to the combination daratumumab-lenalidomide therapy.
To validate these observations, we employed flow cytometry to assess in MM cells the expression of inflammation-related markers, such as KLF4, HLA-DR, and IRF1 [ 18 , 19 ]. Importantly, these markers were robustly upregulated in MM cells from NR patients (Fig. 5 g). This observation, supports scRNA-seq data indicating that BM of NR patients are characterized by a more inflammatory TME, compared to CR.
We further observed a correlation between the inflammation signature of the immune microenvironment and NF-κB activation in tumor cells (Fig. 5 h). Text mining and the co-expression databases STRING revealed that the cytokine genes TNF and IL1B were co-regulated with NF-κB signaling genes in the tumor cells of NR patients (Fig. 5 i). Moreover, we found that TNF and IL1B were enriched in the BME, especially in T, NK and myeloid cells of MM patients (Fig. 5 j). These findings illustrate that the enrichment of NF-κB signaling within tumor cells cooccurs with an inflammatory microenvironment, and is associated with poor response in NR patients.
Predictive models of drug response
Having identified 34 genes that were enriched in the NR group and positively regulate NF-κB signaling (Supp. Figure 8 g), we first examined their expression at the single-cell level, where NR patients consistently showed higher than CR patients (Fig. 6 a). In the Multiple Myeloma Research Foundation (MMRF) CoMMpass cohort ( n = 853), patients with above the median expression of these 34 genes exhibited significantly worse survival (Fig. 6 b). This finding supports the prognostic significance of the 34-gene NF-κB signature in MM, indicating that elevated expression of these NF-κB signaling genes correlates with adverse clinical outcomes.
To further test the ability of the model for the 34 genes in single-cell data, we conducted a random forest (RF) classifier. Using 80% tumor cells (10,342 tumor cells) as the training set and the remaining 20% tumor cells (2,586 tumor cells) as the test set, the RF classifier achieved an AUC (area under curve) of 1.0 on the training set. The test set obtained an AUC value of 0.96, demonstrating high accuracy in distinguishing CR from NR patients (Fig. 6 c). To address potential patient-level information leakage, we further performed leave-one-patient-out cross-validation, in which cells from each patient were independently held out for testing. This analysis achieved an AUC of 0.75 across 14 patients (Supp. Figure 9 a). These analyses were used to evaluate the single-cell discriminative capacity of the signature at the single-cell level and was not used for feature selection in bulk-level modeling.
At the bulk level, gene importance was re-estimated using only the 11 baseline bone marrow RNA-seq samples (Fig. 6 d), followed by forward feature selection to construct candidate models. To avoid overfitting given the limited sample size, we applied leave-one-out cross-validation (LOOCV) during model training. Candidate models containing two or more genes achieved LOOCV AUCs of 1.0 in the discovery cohort and were subsequently evaluated in 14 pseudobulk samples as an independent assessment, achieving an AUC of approximately 0.8 (Fig. 6 d). Guided by the principle of minimal feature complexity, we refined the model to eight genes ( VAPA , NOD2 , CD74 , CCR7 , BCL10 , SLC44A2 , RHOA , LAPTM5 ) based solely on the discovery cohort. The resulting signature showed clear separation between CR and NR groups (Fig. 6 e-f) and was subsequently evaluated for prognostic association in the CoMMpass cohort and the subgroup of 177 cases receiving lenalidomide and/or daratumumab treatment (Fig. 6 g, Supp. Figure 9 b).
This eight-gene classifier achieved an AUC of 0.92 in 14 pseudobulk samples (Fig. 6 h) and demonstrated an AUC of 0.81 when applied to an external cohort of 28 MM patients profiled by Affymetrix microarray (Fig. 6 i). Kaplan–Meier analysis of this external cohort demonstrated a significant survival difference between high- and low-signature groups ( P = 0.001; Fig. 6 j), further supporting the prognostic relevance of the NF-κB signature.
Finally, to explore a more clinically accessible approach, a RF model was trained using the same eight-gene NF-κB signature on 40 peripheral blood mononuclear cell (PBMC) samples from newly diagnosed MM patients in the IFM2017-03 trial using 5-fold cross-validation, achieved an AUC of 0.77 (Fig. 6 k). These results indicate that the NF-κB–based gene signature retains predictive potential across different sample types, providing a framework for anticipating treatment response in MM.
Discussion
In this study, we conducted a comprehensive single-cell analysis of the BME in elderly, frail NDMM patients treated with DR in the IFM2017-03 trial. We identified distinct immune and molecular differences between the CR and NR groups. NR patients exhibited an inflamed microenvironment characterized by elevated inflammatory cytokines and IFN signatures, predominantly within myeloid and NK cell populations. NF-κB signaling was also significantly enriched in NR patients and correlated with altered immune states, including decreased NK cell cytotoxicity-related gene signatures and reduced ADCP-related gene signatures in classical monocytes. Using these findings, we developed a predictive model based on NF-κB-related genes that accurately classified patients into CR or NR groups, highlighting a potential role of NF-κB-associated inflammatory states in shaping response to DR therapy.
Previous single-cell studies have described IFN-responsive immune states in MM, including ISG-high transcriptional programs in T cells and cMono [ 12 ]. Our study extends these observations by showing that these IFN-responsive immune states are associated with NR to DR therapy, with ISG+ cMono displaying enhanced inflammatory features. The upregulation of NF-κB signaling that we observed in NR patients is also supported by existing literature, where this pathway has been implicated in promoting cancer cell proliferation, survival and resistance to apoptosis. Additionally, the reduced functionality of NK cells in NR patients – particularly the decreased expression of CD16 on cytotoxic NKdim cells – is in line with previous reports of impaired NK cell activity in MM patients compared to healthy donors [ 13 ].
One critical question is whether the activation of NF-κB in MM cells precedes or follows the development of an inflammatory microenvironment. It remains unclear whether NF-κB activation in tumor cells initiates the release of proinflammatory cytokines such as TNF and IL-1β, which subsequently shape the immune landscape, or whether the inflamed microenvironment itself triggers NF-κB activation within the tumor cells. Some studies have suggested that tumor cells, through NF-κB activation, can induce immune cells towards a pro-tumorigenic phenotype [ 21 , 22 ], while others propose that an already inflamed microenvironment can drive NF-κB activation in tumor cells, promoting resistance to treatment [ 23 , 24 ]. Future studies in precursor stages could help to resolve this question. Given the observed inflammatory microenvironment in NR patients, the well-established anti-inflammatory properties of dexamethasone – particularly its inhibition of NF-κB activation and proinflammatory cytokine production – may function as a critical modulator of treatment efficacy [ 25 ]. This suggests that inflammatory profiling could inform optimal corticosteroid dosing strategies within DR-containing regimens to maximize therapeutic outcomes.
Our predictive model, which incorporated key NF-κB signaling genes, achieved high accuracy in identifying NR patients, paving the way for developing tools to enable personalized treatment strategies. Validation in the CoMMpass cohort further supported the prognostic relevance of the NF-κB signature across a clinically heterogeneous MM population, although additional uniformly treated DR cohorts will be needed to evaluate its predictive utility for DR response. Beyond predicting responses, another approach could involve targeting inflammation pathways – through cytokines or NF-κB activation – may improve the efficacy of MM treatment. Future studies are needed to determine whether inhibition of proinflammatory cytokines such as TNF or IL-1β, or modulation of NF-κB signaling, can enhance responses to DR-based therapies.
While our study primarily focused on multiple myeloma, the role of IFN signaling in driving resistance to immunotherapy has been extensively studied in various cancer types, particularly in solid tumors. Several studies have shown that chronic type I IFN signaling can promote resistance to checkpoint blockade immunotherapy [ 26 , 27 ] and that tumor-intrinsic type I IFN signaling can induce a T cell-low state associated with resistance to immunotherapy [ 28 , 29 ]. Persistent STAT1 activation or sustained type I IFN signaling can inhibit NK cell proliferation and functions including cytolytic capacity and IFN-γ production, which could further compromise NK cell-mediated anti-tumor immune responses [ 30 ]. Although these studies were primarily conducted in preclinical solid tumor models, the high levels of IFN signaling we observed in myeloid cells of NR patients suggest that similar mechanisms may be at play in MM, underscoring the importance of the inflamed microenvironment in driving resistance to immunotherapy across different cancer types [ 30 , 31 ].
Despite the valuable insights provided by our study, several limitations should be acknowledged. Although our cohort was sufficient to identify statistically significant transcriptomic differences associated with treatment response, it remains relatively small given the biological and clinical heterogeneity of multiple myeloma. We evaluated key features associated with DR activity, including plasma cell CD38 expression, CD16-related immune effector programs, and ADCP-related signatures, larger cohorts with comprehensive clinical and functional characterization will be required to further define the specificity of these features in DR response. Expanding the cohort in future studies will be essential to further validate our findings and refine the predictive model. Second, our study focused on the pretreatment BM context; therefore, it does not fully reflect the dynamic interactions between immune and tumor cells throughout the course of treatment. To gain a more comprehensive understanding of immune responses to DR therapy, longitudinal sampling of the BME would be required. Finally, future studies should integrate additional approaches, such as genomic or epigenetic profiling, to uncover the broader regulatory mechanisms underlying NF-κB-mediated drug resistance and identify novel therapeutic targets.
In conclusion, our study provides novel insights into the role of the inflammatory microenvironment and NF-κB signaling in shaping inflammatory states associated with poor response to DR therapy. By identifying key immune cell dysfunctions and developing a predictive model based on NF-κB-related genes, we have laid the foundation for the development of more personalized therapeutic strategies aimed at improving treatment stratification and outcomes for patients with poor responses.
Methods and materials
Clinical trial design and patient selection
The IFM 2017-03 trial was conceived as a multicenter, international, open-label, randomized, parallel-group controlled study conducted in so-called “frail” patients aged ≥ 65 years with previously untreated MM who are not eligible for autologous hematopoietic stem cell transplantation following high-dose chemotherapy conditioning. The study included two arms: arm A (Rd), representing the standard treatment with lenalidomide and dexamethasone in 28-day cycles, and arm B (DR), representing the investigational combination of subcutaneous daratumumab and lenalidomide in 28-day cycles with dexamethasone administered only during the first two cycles.
According to the anonymized response data from the R-DaraSC arm, 21 samples were selected based on patients’ response to the treatment, as typically assessed using the criteria defined by the International Myeloma Working Group (IMWG). Patients referred to as complete responders (CR) showed at least a Very Good Partial Response (VGPR) by the start of the 3rd treatment cycle, thus including patients in VGPR or in CR, whether stringent or not. The non-responder (NR) group of patients showed a < 75% reduction in the monoclonal component by the start of the 7th treatment cycle, thus including patients with progressive disease, stable disease, or a partial response. Importantly, the two response groups showed comparable baseline clinical and cytogenetic characteristics, including age, sex, ISS stage, and detected high-risk cytogenetic abnormalities, with no statistically significant differences observed (Supplementary Table 1 ).
Sample preparation for single-cell RNA sequencing
Immune cells were enriched from fresh BM samples by Ficoll separation and then frozen in FBS supplemented with 10% DMSO until further processing. Frozen samples were thawed at 37 °C and diluted in RPMI supplemented with 20% FBS and 10 mg/ml DNAse I (StemCell). Cells were pelleted for 10 min at 300 g and resuspended in PBS/0.5% BSA. Dead cells were removed using the Dead Cell Removal Kit (Miltenyi) following manufacturer instructions. The live cell fraction was subsequently incubated with 20 ml anti-CD138 and 20 ml anti-CD235a MicroBeads (Miltenyi). CD138 + CD235 + and CD138 − CD235a − fractions were obtained through magnetic isolation using MS columns (Miltenyi) following manufacturer instructions. Cell viability was checked by Calcein/Draq7 staining and flow cytometry, and the CD138 − CD235a − fraction was then labelled with antibodies for analysis on a BD Rhapsody system following manufacturer instructions. Briefly, cells were incubated with Human BD Fc Block and then co-labelled with a Sample Tag and 35 AbSeq antibodies for surface protein expression profiling according to the manufacturer’s recommendations (Single Cell Labelling with the BD™ Single-Cell Multiplexing Kit and BD™ AbSeq Ab-Oligos Doc ID: 214419 Rev. 1.0). The AbSeq antibodies specific for the following markers were used (all from BD Biosciences): CD127 (Cat# 940012, RRID: AB_2875903), CD19 (Cat# 940247, RRID: AB_2876128), CD45 (Cat# 940002, RRID: AB_2875893), CD38 (Cat# 940013, RRID: AB_2875904), CD11c (Cat# 940265, RRID: AB_2876144), IgD (Cat# 940026, RRID: AB_2875917), CD123 (Cat# 940020, RRID: AB_2875911), CD16 (Cat# 940006, RRID: AB_2875897), CD194 (Cat# 940047, RRID: AB_2875938), CD3 (Cat# 940000, RRID: AB_2875891), CD45RA (Cat# 940011, RRID: AB_2875902), CD27 (Cat# 940018, RRID: AB_2875909), CD28 (Cat# 940017, RRID: AB_2875908), CD183 (Cat# 940030, RRID: AB_2875921), CD161 (Cat# 940283, RRID: AB_2876160), CD185 (Cat# 940042, RRID: AB_2875933), CD24 (Cat# 940028, RRID: AB_2875919), CD45RO (Cat# 940022, RRID: AB_2875913), CD197 (Cat# 940014, RRID: AB_2875905), CD8 (Cat# 940003, RRID: AB_2875894), CD25 (Cat# 940463, RRID: AB_2876314), CD20 (Cat# 940016, RRID: AB_2875907), HLA-DR (Cat# 940010, RRID: AB_2875901), CD4 (Cat# 940001, RRID: AB_2875892), CD56 (Cat# 940007, RRID: AB_2875898), CD196 (Cat# 940033, RRID: AB_2875924), CD14 (Cat# 940005, RRID: AB_2875896), CD138 (Cat# 940213, RRID: AB_2876095), CD34 (Cat# 940264, RRID: AB_2876143), CD279 (Cat# 940015, RRID: AB_2875906), TIM-3 (Cat# 940066, RRID: AB_2875957), LAG-3 (Cat# 940080, RRID: AB_2875971), CD152 (Cat# 940034, RRID: AB_2875925), CD11b (Cat# 940008, RRID: AB_2875899), TCRgd (Cat# 940057, RRID: AB_2875948). After labelling and washing steps, the cells were stained with viability markers (Calcein AM - Thermo Fisher Scientific Cat# C1430; Draq7 - BD Biosciences Cat# 564904) to estimate cell concentration and viability using the BD Rhapsody Scanner (BD Biosciences, La Jolla, CA, USA). The viability was estimated to be between 64 and 93% (average of 83%).
Single-cell capture and cDNA synthesis
The single-cell isolation using the BD Rhapsody platform and cDNA preparation were performed according to the manufacturer’s instructions (Single Cell Capture and cDNA Synthesis with the BD Rhapsody Single-Cell Analysis System Doc ID: 210966 Rev. 1.0). The use of a unique sample Tag allows combining and loading multiple samples on a single BD Rhapsody cartridge. In this study, three samples were loaded per cartridge to reach the target of 8,000 cells per sample, and seven cartridges were used to analyze the 21 samples selected. After quantity and quality controls, labelled samples were pooled in equal proportions in cold Sample Buffer (BD Biosciences Cat# 650000062) for a total of 35,000 cells. Single cells were isolated in cartridge microwells, with magnetic beads for polyA-based mRNA capture. At the loading step, the Rhapsody scanner detected (on average for the seven cartridges) 28,000 wells with viable cells and 23,000 wells with viable cells and a magnetic bead, which theoretically corresponds to 7,700 cells per sample, a value close to our target. After cell lysis, beads were magnetically retrieved from the microwells and pooled into a single tube before reversing transcription of the polyA-captured mRNA, AbSeq and Sample Tag.
Preparation of the whole transcriptome mRNA, AbSeq and Sample Tag libraries
After cDNA synthesis, two parallel workflows were performed following the manufacturer’s protocol (mRNA Whole Transcriptome Analysis [WTA], AbSeq, and Sample Tag Library Preparation Protocol 23-21752-00 10/2019): one for the whole transcriptome libraries using a random priming approach and another for the Abseq and Sample Tag libraries based on PCR amplifications. Final libraries were quantified with the Qubit dsDNA HS assay kit (Invitrogen Cat# Q32851 ) and the quality was analyzed with a BioAnalyzer instrument (Agilent Technologies, RRID SCR_018043) using a High Sensitivity DNA chip (Agilent Technologies Cat# 5067 − 4626). Before pooling, library concentration was normalized to 4nM.
Library pooling and sequencing
WTA, AbSeq and Sample Tag libraries were pooled and sequenced on a NovaSeq™ 6000 sequencer (Illumina, San Diego, California, United States) using the NovaSeq 6000 S4 Reagent Kit v1.5 (200 cycles) (Illumina, Cat# 20028313) with a loading concentration of 300pM and the following configuration: 75 bp for Read 1, 75 bp for Read 2 and 8 bp for index. For pooling, the library ratio was calculated to target a sequencing depth of 50,000 reads/cell for WTA, 17,500 reads/cell for AbSeq and 600 reads/cell for Sample Tag.
Sample preparation for bulk RNAseq
CD138 + cells were isolated from BM aspirates using anti-CD138 MicroBeads (Miltenyi), after which RNA was extracted. Library preparation and sequencing were performed by Novogene.
Flow cytometry analysis of BM aspirates
Samples (12 CR and 11 NR) were thawed at 37 °C, and 2 × 10 6 cells were stained with 1:1,000 with Zombie Aqua Fixable Viability Kit (BioLegend, Cat# 423101) for 20 min at 4 °C. Fc receptors were blocked using human TruStain FcX (BioLegend, Cat# 422301) following manufacturer’s instructions. Samples were then fixed with 0.2% formaldehyde for 20 min at 4 °C and stained overnight in Permeabilization buffer (eBioscience, Cat# 00-8333-56) as described by Whyte et al. [ 31 ]. The following antibodies were used for staining at 1:500 final concentration: anti-human CD45 APC-Cy7 (clone 2D1, BioLegend, Cat# 368515, RRID: AB_2566375), anti-human CD3 BUV395 (clone UCHT1, BD Biosciences, Cat# 563546, RRID: AB_2744387), anti-human CD56 FITC (clone HCD56, BioLegend, Cat# 318304, RRID: AB_604100), and anti-human CD16 BV605 (clone 3G8, BioLegend, Cat# 302039, RRID: AB_2561354). Alternatively, viability staining was performed with 1:1000 LIVE/DEAD fixable Blue viability kit (ThermoFisher Scientific, Cat# L23105 ), and after blocking the Fc receptors surface staining with anti-human HLA-DR PE-Fire 810 (clone L243, BioLegend, Cat# 307683, RRID: AB_ 2904336), anti-human CD138 PE-Fire 700 (clone MI15, BioLegend, Cat# 356553, RRID: AB_3083343), and the anti-human CD16 APC-Fire 810 (clone 3G8, BioLegend, Cat# 302073, RRID: AB_2910378), anti-human CD27 KIRAVIA blue 520 (clone M-T271, BioLegend, Cat# 356435, RRID: AB_2832700), anti-human CCR7 BV421 (clone G043H7, BioLegend, Cat# 353207, RRID: AB_10915137). Samples were then fixed and permeabilized using the eBioscience Foxp3/Transcription Factor Staining Buffer Set (ThermoFisher Scientific, Cat# 00-5523-00), and stained with anti-human CD45 Spark PLUS UV 395 (clone HI30, BioLegend, Cat# 304095, RRID: AB_3097561), anti-human CD3 BV510 (BioLegend, Cat# 300447, RRID: AB_2563467), anti-human CD19 Spark NIR685 (clone HIB19, BioLegend, Cat# 302269, RRID: AB_2860769), anti-human CD8 BV605 (clone HIT8a, BioLegend, Cat# 300935, RRID: AB_2832574), anti-human CD4 APC-Fire 750 (clone A161A1, BioLegend, Cat# 357425, RRID: AB_2716182), anti-human CD56 PE-Dazzle 594 (clone 5.1H11, BioLegend, Cat# 362543, PRRID: AB_2565921), anti-human CD11b StarBright UltraViolet 740 (clone M1/70, BioLegend, Cat# 101274), anti-human CD38 BV785 (clone HIT2, BioLegend, Cat# 303529, RRID: AB_2561368), anti-human CD14 BV711 (clone 63D3, BioLegend, Cat# 367139, RRIO: AB_2814322), anti-human KLF4 APC (R&D Systems, Cat# IC3640A, RRID: AB_2044690), and anti-human IRF1 PE (clone D5E4, Cell Signaling Technologies, Cat# 12732, RRID: AB_2798008). Samples were acquired using an LSRFortessa X-20 (BD Biosciences) or the Aurora CS (Cytek Biosciences). Data unmixing was performed using AutoSpectral R package (drcytometer.github.io/AutoSpectral/) and data analysis was performed using OMIQ (Dotmatics).
Bulk RNA-seq data analyses
The bulk RNA-seq data sequencing counts of each sample were greater than 20 million. Raw counts were utilized for normalization, principal component analysis (PCA) and enrichment analysis. PCA was based on a normalized matrix of counts per million (CPM) displaying a two-dimensional plane of PC1 and PC2 coordinates. Heatmaps were presented using the Complex Heatmap package (version 2.12.1) with specific loading genes in order of PC1 coordinates for each patient, including 100 CR group genes and 100 NR group genes. The specific genes were then enriched separately by the cluster Profiler package (version 4.4.4), and the highly enriched entries that were presented met a significance level of q -value < 0.05 including at least four genes.
scRNA-seq analyses
Raw fastq data were mapped to the human reference genome (GRCh38) via the standard BD Rhapsody workflow according to the manufacturer’s recommendations. Unique molecular identifier (UMI) count matrices were imported into Seurat (R package, versions 4.2.1 and 4.1.1) [ 32 ] for quality control that included selecting cells with a library complexity of more than 500 features, less than 50,000 counts, filtering out cells with high percentages of mitochondrial genes (> 25%) and removing multiple and undetermined cells. To account for differences in sequencing depth across cells, RNA UMI counts were log-normalized, while antibody UMI counts were applied to a centered log ratio transformation to account for unspecific binding background signal.
Healthy BM data
To provide additional comparisons for the patients collected in this study including the responder and non-responder groups, scRNA-seq data including 20 healthy individuals published by Oetjen et al. were collected and analyzed, and the original sequencing files were downloaded from GEO ( GSE120221 ) [ 15 ].
scRNA-seq data integration
The scRNA-seq dataset of 21 patients and 20 normal BM samples was merged by integration, normalized and analyzed by PCA across all cells using the Seurat (v4.1.1, RRID: SCR_016341) R package [ 32 ]. We controlled batch effects using the FindIntegrationAnchors function to identify anchors and the IntegrateData function for integrating the dataset base on anchors. The top 2,000 variable genes from each dataset and the first 30 PCs were used for integration. Based on the integration results, the top 20 PCs were clustered unsupervised, using a shared nearest neighbor modular optimized clustering algorithm integrated in Seurat (resolutions 1–2), and cells were projected in two-dimensional coordinates of UMAP.
Cell type annotation
Cell types were annotated based on manual review of marker genes (as above). The subgroup clustering results were decisive for the identification of cell type identity. Mixed cell populations containing two apparent subgroups in a cluster due to clustering problems were further subdivided and defined. Cell populations expressing mixed marker genes (e.g., highly expressing both erythroid and immune cell features) were considered mixed cells and were excluded from subsequent analysis. Cluster-specific highly expressed gene analyses were performed using the FindAllMarkers function of the Seurat package. Only genes expressed in at least 25% of cells in either populations were tested. Adjusted p -values < 0.05 were considered significant using the default Wilcoxon test.
Module calculation based on cell proportions
First, the proportions of different cell subpopulations were calculated across MM patient response groups and the HD group. The table function was used to determine the proportion of each fine subpopulation within each Response group, followed by normalization. A heatmap was then generated using pheatmap to visualize the distribution of cell subpopulations under different response conditions. Additionally, hierarchical clustering was performed on both samples and cell subpopulations, and the column order of each module was extracted based on the clustering results.
For differentially expressed gene (DEG) analysis of Module4, we applied FindMarkers to identify genes that were differentially expressed between the CR and NR groups, using the thresholds min.pct = 0.25 and logFC.threshold = 0.25. Genes with significant differential expressions (p_val_adj < 0.05 and |avg_log2FC| ≥ 0.25) were further selected. ggplot2 was used to visualize gene expression changes, highlighting genes with higher expression in either the CR or NR groups.
GSEA-based functional enrichment analysis of Module4
Following the identification of Module4-associated DEGs, we performed gene set enrichment analysis (GSEA) to assess functional enrichment. First, genes were ranked in descending order based on avg_log2FC, and bitr was used to convert gene symbols (SYMBOL) to Entrez IDs to meet the input requirements for GSEA. Next, the Hallmark gene sets were retrieved using msigdbr (RRID: SCR_022870), and GSEA was conducted using the GSEA function with a significance threshold of pvalueCutoff = 0.05 to identify significantly enriched gene sets. Finally, gseaplot2 was employed to visualize the enrichment curves of the most significantly enriched Hallmark pathways, providing insights into the inflammation-related characteristics of Module4 in the NR group.
Gene set enrichment scores
To assess the corresponding properties of cell clusters or different response groups, we evaluated NF-κB signaling, cytokine and inflammation scores for each cell with the ‘AddModuleScore’ function in Seurat v4.1.1. Inflammatory signature genes were obtained from the gene set ‘HALLMARK_INFLAMMATORY_RESPONSE’ downloaded from the Molecular Signatures Database (MSigDB, RRID: SCR_016863) [ 33 ], cytokine genes were obtained from the GO database called ‘positive regulation of cytokine production’ (GO: 0001819), and the NF-κB signaling gene set included ‘I kappaB kinase/NF-kappaB signaling’ (GO: 0007249) and ‘positive regulation of I-kappaB kinase/NF-kappaB signaling’ (GO: 0043123). Additionally, the effector gene set includes PRF1 , IFNG , NKG7 , GZMB , GZMA , GZMK , GZMD , KLRK1 , KLRB1 , KLRD1 , CTSW , and CST7 .
Gene ontology and gene set enrichment analysis
The clusterProfiler package (RRID: SCR_016884) was used to assess the enrichment of the top ranked genes in Gene Ontology biological processes categories. Significantly enriched categories were identified with adjusted p -values < 0.05. Group-specific genes for CR and NR in single-cell clusters were calculated according to the FindMarkers function in Seurat v4.1.1, with thresholds of expression in at least 25% of cells, log 2 fold change ≥ 0.25 and adjusted p -value < 0.05.
RNA velocity and pseudo-time analysis
RNA velocity analysis was performed with scVelo (version 0.1.25, RRID: SCR_018168) [ 34 ] in Python. Bam files obtained from BD flow-based analysis were used to generate loom files. Briefly, after gene selection and normalization, first- and second-order moments were calculated using the ‘scv.pp.ments’ function. The complete splicing dynamics was recovered with the ‘scv.tl.recover_dynamics’ function and the velocity values in dynamic mode were obtained with the ‘scv.tl.velocity’ function. The velocity was projected onto the diffusion map with the ‘scv.pl.velocity_embedding_stream’ function and displayed as a streamline. The spliced and unspliced phase maps of individual genes were visualized with the ‘scv.pl.velocity’ function. The pseudotime of the cell was obtained with the ‘scv.tl.recover_latent_time’ function.
scRNA-seq correlation analysis
For each type of defined myeloid cells, we calculated the average expression values of all genes for each subpopulation using the AverageExpression function, and the top 2,000 genes with the largest standard deviation were prioritized to calculate the Pearson correlation coefficients.
Detection of CNVs in cancer cells
We inferred CNVs of the 21 patients by InferCNV using scRNA-seq data [ 35 ]. Genes with a mean UMI > 0.1 were selected for calculation using plasma cells from normal samples as a reference and all plasma cells from patients as a test set, as described in InferCNV (RRID: SCR_021140). In brief, the heatmap shows the variation score for each chromosome based on genomic position according to regions with a window of every 101 genes. Plasma cells from patients all exhibited large CNVs and clarified the significantly high variant status of NRs.
Pseudobulk data analyses
Pseudobulk analysis was performed on a patient-by-patient basis by aggregating gene expression counts across all malignant plasma cells from each patient to construct a new expression matrix. First, MMPCs were extracted based on cell type annotations. To ensure robust analysis, patients with < 10 detected MMPCs were excluded. For each retained patient, raw gene counts were summed across all MMPCs within that patient to generate a patient-level expression profile. To filter lowly expressed genes, only genes with CPM values > 1 in at least two patients were retained. The resulting count matrix was then log-transformed using log 2 (CPM + 1) normalization. Subsequent PCA and enrichment analyses were consistent with the bulk data analysis process described above.
Protein-protein interaction network analysis
To investigate the interaction between NF-κB-related genes and the key inflammatory cytokines IL1B and TNF , we constructed a protein-protein interaction (PPI) network. The STRING database ( https://string-db.org/ , RRID: SCR_005223) was utilized to retrieve known and predicted protein-protein interactions. The interaction network was generated by inputting the 34 NF-κB-related genes, which were identified through a series of computational analyses, along with IL1B and TNF . A confidence score threshold of 0.4 was applied to filter for high-confidence interactions. The resulting interaction network was then imported into Cytoscape (RRID: SCR_003032) for visualization and further analysis.
Single-cell random forest model construction (bone marrow)
We used scRNA-seq data derived from BM plasma cells of MM patients to construct a classification model based on 34 NF-κB-related genes. Data preprocessing included gene length normalization and TPM transformation. The dataset was split into a training set (80%) and a test set (20%), with response labels defined as CR and NR. An RF classifier was trained on the training set, and its performance was evaluated using ROC curves and AUC. The model demonstrated high predictive power, achieving an AUC of 1.0 on the training set and 0.96 on the test set, indicating strong generalizability.
Bulk RNA-seq random forest model construction (bone marrow)
Like the single-cell analysis, we applied gene length normalization and TPM transformation before performing feature selection using the RF algorithm. Eight key genes were identified as the most informative for model construction. A predictive model was developed using bulk RNA-seq data and validated through leave-one-out cross-validation, achieving an AUC of 1.0. To further assess its robustness, the model was tested on pseudobulk samples. This yielded an AUC of 0.92, demonstrating strong predictive performance.
Independent validation and preprocessing of Affymetrix microarray data (bone marrow)
To evaluate the predictive performance of the NF-κB–related gene signature in an independent external cohort, we applied the previously trained eight-gene random forest model to a set of 28 multiple myeloma patients treated at Montpellier University Hospital. Among these patients, 57% (16/28) achieved CR, while 43% (12/28) were NR (Supp. Table 3 ). Prior to initiation of daratumumab treatment, multiparametric flow cytometry of tumor plasma cells was performed for all patients. Gene expression was profiled using Affymetrix U133 Plus 2.0 microarrays, with RNA extracted following standard protocols (RNeasy Kit, Qiagen). Expression values were log2-transformed and z-score standardized across samples, and probes with missing values were excluded.
PBMC-based model construction for clinical applicability
To enhance the clinical applicability of our findings, we extended the model to PBMC-derived patient samples. We selected the same eight key genes identified in the BM bulk RNA-seq analysis. Using data from 40 PBMC samples, we trained an RF model and evaluated its performance through five-fold cross-validation. The model achieved an AUC of 0.77, suggesting its potential utility in peripheral blood-based patient stratification.
Across all model analyses, we employed a consistent data preprocessing workflow including gene length normalization, TPM transformation, and RF-based feature selection, ensuring methodological coherence across BM scRNA-seq data, PBMC scRNA-seq data and BM bulk RNA-seq data.
Survival analysis
Survival data were downloaded from MMRF CoMMpass Study (RRID: SCR_003721) and included a total of 853 patients with MM and 172 patients treated with lenalidomide or daratumumab. The R packages survivor and survminer (RRID: SCR_026512) were used for prognostic assessment; specifically, the obtained NF-κB signaling gene and patient survival information were used as inputs to generate patient survival curves and statistical tests for overall survival. Adjusted p -values < 0.05 were considered statistically significant.
Supplementary Information
Below is the link to the electronic supplementary material.
Supplementary Material 1
Acknowledgements
We thank the facilities UMS2014-US41-PLBS BioImaging Center Lille, Flow Cytometry Core Facility, and Génomique Fonctionelle at Structurale (GFS), F-59000 Lille, France, for technical support. Figure 1 A was created with BioRender.com. We thank Dr. Daniel Ackerman of Insight Editing London for editing the manuscript prior to submission.
Author contributions
Conceptualization: SMa, SMi. Methodology: WC, SG, AO, CV. Investigation: WC, SG. Visualization: WC. Resources: LHBI, AO, GE, AC, YTN, NC, EW, GD, RT, BP, ML, MS, AK, SGD, JM, BA, JC, CV, LM, XL, MF, TF. Funding acquisition: SMa. Project administration: SMa, SMi, WJ. Supervision: SMa, SMi, WJ. Writing – original draft: SMa, SMi, WJ.
Funding
This study was financed by Johnson and Johnson and the Intergroup Francophone du Myélome (IFM).
Data availability
All single-cell RNA, surface protein, and bulk sequencing data generated in this study have been deposited in the BioStudies under the reference S-BSST2098 and in the European Genome-Phenome Archive (EGA) under the reference EGAS00001008201. The publicly available scRNA-seq dataset from 20 healthy individuals ( GSE120221 ) was used for comparison. All R code used for the analysis is available on GitHub at https://github.com/Chengwenwen1/MM-France/tree/main .
Declarations
Ethical approval
The French National Agency for the Safety of Medicines and Health (ANSM) and the Ethical Committee “Sud-Mediterranee II” (reference 219 A10) approved the protocol (ID# 19.02.28.43641). The IFM2017-03 trial (registered as NCT03993912 on 2019-10-17, https://www.clinicaltrials.gov/ ) and the correlative studies were conducted in accordance with the principles of the Declaration of Helsinki and the International Conference on Harmonization Good Clinical Practice guidelines. All the patients provided written informed consent. The trial was promoted by Lille University Hospital (CHU Lille) in collaboration with the IFM and with the support of Johnson and Johnson.
Consent to participate
Informed consent was obtained from all participants in this study.
Competing interests
The authors declare no competing interests.
Footnotes
Publisher’s note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Wenwen Cheng and Silvia Gaggero contributed equally to this work.
Wenfei Jin, Suman Mitra and Salomon Manier are co-senior authors.
Contributor Information
Suman Mitra, Email: suman.mitra@inserm.fr.
Salomon Manier, Email: salomon.manier@inserm.fr.
References
1. Dimopoulos MA, Oriol A, Nahi H, San-Miguel J, Bahlis NJ, Usmani SZ, et al. Daratumumab, Lenalidomide, and Dexamethasone for Multiple Myeloma. N Engl J Med. 2016;375(14):1319–31. 10.1056/NEJMoa1607751.
2. Facon T, Kumar S, Plesner T, Orlowski RZ, Moreau P, Bahlis N, et al. Daratumumab plus Lenalidomide and Dexamethasone for Untreated Myeloma. N Engl J Med. 2019;380(22):2104–15. 10.1056/NEJMoa1817249.
3. Facon T, Leleu X, Manier S. How I treat multiple myeloma in geriatric patients. Blood. 2024;143(3):224–32. 10.1182/blood.2022017635.
4. Sanchez L, Wang Y, Siegel DS, Wang ML. Daratumumab: a first-in-class CD38 monoclonal antibody for the treatment of multiple myeloma. J Hematol Oncol. 2016;9(1):51. 10.1186/s13045-016-0283-0.
5. Liu J, Xing L, Li J, Wen K, Liu N, Liu Y, et al. Epigenetic regulation of CD38/CD48 by KDM6A mediates NK cell response in multiple myeloma. Nat Commun. 2024;15(1):1367. 10.1038/s41467-024-45561-z.
6. Krejcik J, Casneuf T, Nijhof IS, Verbist B, Bald J, Plesner T, et al. Daratumumab depletes CD38 + immune regulatory cells, promotes T-cell expansion, and skews T-cell repertoire in multiple myeloma. Blood. 2016;128(3):384–94. 10.1182/blood-2015-12-687749.
7. Fink EC, Ebert BL. The novel mechanism of lenalidomide activity. Blood. 2015;126(21):2366–9. 10.1182/blood-2015-07-567958.
8. Kronke J, Udeshi ND, Narla A, Grauman P, Hurst SN, McConkey M, et al. Lenalidomide causes selective degradation of IKZF1 and IKZF3 in multiple myeloma cells. Science. 2014;343(6168):301–5. 10.1126/science.1244851.
9. Casneuf T, Adams HC 3rd, van de Donk N, Abraham Y, Bald J, Vanhoof G, et al. Deep immune profiling of patients treated with lenalidomide and dexamethasone with or without daratumumab. Leukemia. 2021;35(2):573–84. 10.1038/s41375-020-0855-4.
10. Schinke C, Poos AM, Bauer M, John L, Johnson S, Deshpande S, et al. Characterizing the role of the immune microenvironment in multiple myeloma progression at a single-cell level. Blood Adv. 2022;6(22):5873–83. 10.1182/bloodadvances.2022007217.
11. Cheng Y, Sun F, Alapat DV, Wanchai V, Mery D, Siegel ER, et al. Multi-omics reveal immune microenvironment alterations in multiple myeloma and its precursor stages. Blood Cancer J. 2024;14(1):194. 10.1038/s41408-024-01172-x.
12. Zavidij O, Haradhvala NJ, Mouhieddine TH, Sklavenitis-Pistofidis R, Cai S, Reidy M, et al. Single-cell RNA sequencing reveals compromised immune microenvironment in precursor stages of multiple myeloma. Nat Cancer. 2020;1(5):493–506. 10.1038/s43018-020-0053-3.
13. Blanquart E, Ekren R, Rigaud B, Joubert MV, Baylot V, Daunes H, et al. NK cells with adhesion defects and reduced cytotoxic functions are associated with a poor prognosis in multiple myeloma. Blood. 2024;144(12):1271–83. 10.1182/blood.2023023529.
14. de Jong MME, Kellermayer Z, Papazian N, Tahri S, Hofste Op Bruinink D, Hoogenboezem R, et al. The multiple myeloma microenvironment is defined by an inflammatory stromal cell landscape. Nat Immunol. 2021;22(6):769–80. 10.1038/s41590-021-00931-3.
15. Oetjen KA, Lindblad KE, Goswami M, Gui G, Dagur PK, Lai C, et al. Human bone marrow assessment by single-cell RNA sequencing, mass cytometry, and flow cytometry. JCI Insight. 2018;3(23). 10.1172/jci.insight.124928.
16. Tirier SM, Mallm JP, Steiger S, Poos AM, Awwad MHS, Giesen N, et al. Subclone-specific microenvironmental impact and drug response in refractory multiple myeloma revealed by single-cell transcriptomics. Nat Commun. 2021;12(1):6960. 10.1038/s41467-021-26951-z.
17. Akiba H, Nakano H, Nishinaka S, Shindo M, Kobata T, Atsuta M, et al. CD27, a member of the tumor necrosis factor receptor superfamily, activates NF-kappaB and stress-activated protein kinase/c-Jun N-terminal kinase via TRAF2, TRAF5, and NF-kappaB-inducing kinase. J Biol Chem. 1998;273(21):13353–8. 10.1074/jbc.273.21.13353.
18. Feinberg MW, Cao Z, Wara AK, Lebedeva MA, Senbanerjee S, Jain MK. Kruppel-like factor 4 is a mediator of proinflammatory signaling in macrophages. J Biol Chem. 2005;280(46):38247–58. 10.1074/jbc.M509378200.
19. Ivashkiv LB. IFNgamma: signalling, epigenetics and roles in immunity, metabolism, disease and cancer immunotherapy. Nat Rev Immunol. 2018;18(9):545–58. 10.1038/s41577-018-0029-z.
20. Vrabel D, Pour L, Sevcikova S. The impact of NF-kappaB signaling on pathogenesis and current treatment strategies in multiple myeloma. Blood Rev. 2019;34:56–66. 10.1016/j.blre.2018.11.003.
21. Grivennikov SI, Greten FR, Karin M. Immunity, inflammation, and cancer. Cell. 2010;140(6):883–99. 10.1016/j.cell.2010.01.025.
22. Saccani A, Schioppa T, Porta C, Biswas SK, Nebuloni M, Vago L, et al. p50 nuclear factor-kappaB overexpression in tumor-associated macrophages inhibits M1 inflammatory responses and antitumor resistance. Cancer Res. 2006;66(23):11432–40. 10.1158/0008-5472.CAN-06-1867.
23. Luo JL, Maeda S, Hsu LC, Yagita H, Karin M. Inhibition of NF-kappaB in cancer cells converts inflammation- induced tumor growth mediated by TNFalpha to TRAIL-mediated tumor regression. Cancer Cell. 2004;6(3):297–305. 10.1016/j.ccr.2004.08.012.
24. Hagemann T, Lawrence T, McNeish I, Charles KA, Kulbe H, Thompson RG, et al. Re-educating tumor-associated macrophages by targeting NF-kappaB. J Exp Med. 2008;205(6):1261–8. 10.1084/jem.20080108.
25. Auphan N, DiDonato JA, Rosette C, Helmberg A, Karin M. Immunosuppression by glucocorticoids: inhibition of NF-kappa B activity through induction of I kappa B synthesis. Science. 1995;270(5234):286–90. 10.1126/science.270.5234.286.
26. Jacquelot N, Yamazaki T, Roberti MP, Duong CPM, Andrews MC, Verlingue L, et al. Sustained Type I interferon signaling as a mechanism of resistance to PD-1 blockade. Cell Res. 2019;29(10):846–61. 10.1038/s41422-019-0224-x.
27. Zak J, Pratumchai I, Marro BS, Marquardt KL, Zavareh RB, Lairson LL, et al. JAK inhibition enhances checkpoint blockade immunotherapy in patients with Hodgkin lymphoma. Science. 2024;384(6702):eade8520. 10.1126/science.ade8520.
28. Mathew D, Marmarelis ME, Foley C, Bauml JM, Ye D, Ghinnagow R, et al. Combined JAK inhibition and PD-1 immunotherapy for non-small cell lung cancer patients. Science. 2024;384(6702):eadf1329. 10.1126/science.adf1329.
29. Qiu J, Xu B, Ye D, Ren D, Wang S, Benci JL, et al. Cancer cells resistant to immune checkpoint blockade acquire interferon-associated epigenetic memory to sustain T cell dysfunction. Nat Cancer. 2023;4(1):43–61. 10.1038/s43018-022-00490-y.
30. Tabellini G, Vairo D, Scomodon O, Tamassia N, Ferraro RM, Patrizi O, et al. Impaired natural killer cell functions in patients with signal transducer and activator of transcription 1 (STAT1) gain-of-function mutations. J Allergy Clin Immunol. 2017;140(2):553–e644. 10.1016/j.jaci.2016.10.051.
31. Whyte CE, Tumes DJ, Liston A, Burton OT. Do more with Less: Improving High Parameter Cytometry Through Overnight Staining. Curr Protoc. 2022;2(11):e589. 10.1002/cpz1.589.
32. Argelaguet R, Clark SJ, Mohammed H, Stapel LC, Krueger C, Kapourani CA, et al. Multi-omics profiling of mouse gastrulation at single-cell resolution. Nature. 2019;576(7787):487–91. 10.1038/s41586-019-1825-8.
33. Liberzon A, Subramanian A, Pinchback R, Thorvaldsdottir H, Tamayo P, Mesirov JP. Molecular signatures database (MSigDB) 3.0. Bioinformatics. 2011;27(12):1739–40. 10.1093/bioinformatics/btr260.
34. Bergen V, Lange M, Peidli S, Wolf FA, Theis FJ. Generalizing RNA velocity to transient cell states through dynamical modeling. Nat Biotechnol. 2020;38(12):1408–14. 10.1038/s41587-020-0591-3.
35. Tirosh I, Izar B, Prakadan SM, Wadsworth MH 2nd, Treacy D, Trombetta JJ, et al. Dissecting the multicellular ecosystem of metastatic melanoma by single-cell RNA-seq. Science. 2016;352(6282):189–96. 10.1126/science.aad0501.
Associated Data
Supplementary Materials
Supplementary Material 1
Data Availability Statement
All single-cell RNA, surface protein, and bulk sequencing data generated in this study have been deposited in the BioStudies under the reference S-BSST2098 and in the European Genome-Phenome Archive (EGA) under the reference EGAS00001008201. The publicly available scRNA-seq dataset from 20 healthy individuals ( GSE120221 ) was used for comparison. All R code used for the analysis is available on GitHub at https://github.com/Chengwenwen1/MM-France/tree/main .