Patient cohort

To delineate the genomic and immune factors contributing to resistance against Dara-Rd, we analyzed data from WGS and flow cytometry in a longitudinal cohort of 32 patients with RRMM enrolled in the NCT03848676 observational study (Fig. 1A). Clonal evolution analysis also included two patients: one treated with daratumumab single agent, and the other with daratumumab, bortezomib, and dexamethasone (Dara-Vd). The median age was 68 (range 46–80); 18 (56.2%) were women. The median number of prior treatments was 1 (range 1–3), with 23 (71.8%) patients previously exposed to high-dose melphalan chemotherapy followed by autologous stem cell transplantation (HDM-ASCT; Fig. 1B). Before starting Dara-Rd, 97% (N = 31) and 50% (N = 16) of patients were previously exposed to PI and IMIDs, respectively. After a median follow-up of 18.2 months (range 4.2–43.4), 40.6% (N = 13) patients experienced progression (i.e., progressors). Patients that did not experience progression were defined as “durable responders”. Overall, 8 (25%) patients achieved minimal residual disease (MRD) negativity by next-generation flow cytometry (Fig. 1B) [17]. The median event-free survival (EFS) was 38.4 months (Fig. 1C). Two patients (AIRC112 and AIRC113) died due pulmonary sarcoma and COVID pneumonia without relapse. Considering all the relapsed patients 9 out of 13 (69.2%) died due disease progression. Combining WGS and FISH data, the presence of at least one high-risk cytogenetic alteration (HRCA) [i.e., del17p and/or t(4;14) and/or t(14;16)] was observed in 8 out of 30 patients (26.7%; FISH and WGS data was missing in 2 patients) and it was associated with poor EFS (P = 0.00044; log-rank test; Fig. 1D and Supplemental Tables 1 and 2) [14]. We did not see an influence of other key clinical and serological features on the clinical outcome (Supplemental Table 2). Demographics, disease characteristics, and response to Dara-Rd for all patients are summarized in Table 1 and Supplemental Table 1.

Genomic determinant of resistance and progression to Dara-Rd

To investigate the intrinsic genomic drivers associated with resistance to Dara-Rd, we performed WGS to profile the genomic landscape in BM samples at baseline and at progression. Four samples out of 32 failed the sequencing due to low cellularity or cancer cell fraction and therefore were excluded. The median coverage for the 28 cases that successfully completed WGS was 70× (Supplemental Table 3). The median single base substitution (SBS) burden was 5800 (range 1808–17,411). There was no difference in total mutational burden between progressors and durable responders. To assess the clinical implications of major somatic drivers in MM, we utilized an extensive catalog of established myeloma genomic drivers derived from large cohorts of NDMM patients, that encompasses mutations, mutational signatures, recurrent aneuploidies, and canonical translocations (Supplemental Methods and Supplemental Table 4) [18, 19]. Among the catalogs of 90 nonsynonymous MM driver mutations [18,19,20], we found 31 mutated genes in the patients with RRMM treated with Dara-Rd (Supplemental Fig. 1A). Genes in MAPK pathway (NRAS, KRAS and BRAF) were the most frequently mutated, with 61% of RRMM cases (treated with Dara-Rd) carrying at least one mutation in one of these genes (Supplemental Fig. 1A). Nevertheless, neither MAPK pathway nor any other mutations in MM driver genes predicted poor anti-CD38 MoAb outcomes (Supplemental Table 5). Finally, we did not find any mutation involving CD38 gene either before or after treatment, but only a large deletion in one patient [21].

To investigate underlying mutational processes involved in these cases, we performed a mutational signature analysis (Supplemental Methods) [22, 23]. We identified seven SBS signatures involved in our cohort of RRMM: the clock-like aging signatures (SBS1 and SBS5), APOBEC mutational activity (SBS2 and SBS13), germline center associated polymerase eta (SBS9), radical oxygen damage (SBS18), SBS8, and SBS-MM1 (i.e., SBS99; melphalan mutagenesis; Fig. 2A and Supplemental Table 6) [19, 23,24,25,26,27,28,29,30]. SBS-MM1/SBS99 mutational signature was identified in 14 out of 19 (74%) patients who received HDM-ASCT. The lack of SBS-MM1/SBS99 in 5 patients can be explained through the single-cell expansion model as previously described [31]. Specifically, it has been shown that distinct chemotherapy agents promote their mutational activity by introducing a unique catalog of mutations in each exposed single cell [25, 31, 32]. These mutations can be detected by bulk WGS only if a single tumor cell exposed to chemotherapy expands, taking clonal dominance. In contrast, if the cancer progression is driven by multiple clones originating from different single cells exposed to chemotherapy, the chemotherapy-induced mutational signature will not be detectable because each clone harbors different catalogs of unique chemotherapy-related variants. Interestingly, two patients (AIRC119 and AIRC128) exhibited two distinct set of melphalan-related variants detected in both subclonal and clonal variants. This scenario was explained by the patients’ exposure to tandem HDM-ASCT. Among the mutational signatures identified in our cohort, high APOBEC contribution was associated to higher rate of progression (P = 0.03 estimated using Fisher exact test) and short EFS (P = 0.047; log-rank test; Fig. 2B and Supplemental Table 6). APOBEC has been shown to be one of the strongest prognostic markers for poor outcomes in NDMM treated with and without daratumumab-based regimens [18, 24, 26, 33]. These findings suggest that Dara-Rd treatment cannot fully counteract the negative prognostic impact of high APOBEC mutational activity. Among the other SBS signatures, the presence of SBS18, known to be caused by reactive oxygen species damage and associated with resistance to CART in diffuse large B-cell lymphoma [34], was associated with shorter EFS, with 4 out of 5 patients (80%) progressing after anti-CD38 MoAb combination treatment (P = 0.03 estimated using log-rank test; Fig. 2C and Supplemental Table 6). Finally, high SBS9 showed favorable outcomes, most likely due to the known inverse correlation between APOBEC and SBS9 (P = 0.041; log-rank test; Supplemental Table 6) [25, 33].

Fig. 2: Impact of genomic alterations on clinical outcome in RRMM treated with Dara-Rd.

A A heatmap showing all the genomic alterations associated with progression after Dara-Rd among patients with available WGS data (N = 28); at the top, the bar plot with the relative contribution per each SBS mutational signature. Minimal residual disease (MRD) was tested using Euroflow [17]. WT wild type, HDM-ASCT high-dose melphalan and autologous stem cell transplant, ROS radical oxygen stress. B–I Kaplan–Meier curves showing the impact of APOBEC (B), SBS18 (C), deletion 1p22.1 (RPL5); D), deletion 10p15.3 (E), deletion 16p13.3 (CREBBP); F), deletion 17p13.1 (TP53); G), structural variants on MYC (MYC, PTV1, and NSMCE2); H) and presence of chromothripsis (I) on event-free survival (EFS); P value is calculated with log-rank test.

Next, we investigated the clinical impact of recurrent MM aneuploidies and structural variants (Supplemental Fig. 1B and Supplemental Methods) [18, 35]. Different deletions were associated with shorter EFS. Specifically, deletion of 1p22.1 (RPL5) was associated with shorter EFS (P = 0.031; log-rank test; Fig. 2D). Mono and bi-allelic loss of RPL5 have been associated with poor response in both Dara and non-Dara-based regimens [33, 36]. In addition to RPL5, patients with loss of 10p15.3, 16p13.3 (CREBBP) and TP53 were characterized by shorter EFS (P = 0.0061, P = 0.034 and P = 0.015; log-rank test, respectively; Fig. 2E–G). Since the identified aneuploidies 16p13.3 and 10p15.3 are genomic event not extensively characterized in the literature, we investigated their biological impact using paired WGS and RNA sequencing data from 705 patient enrolled in the CoMMpass study (https://themmrf.org). Running a gene set enrichment analysis (GSEA), we found that patients harboring one of these two alterations (versus wild-type) were enriched for E2F, MYC targets and G2M checkpoint pathways (Supplemental Fig. 2). These pathways were also found to be enriched in patients with RRMM that progressed after Daratumumab, Carfilzomib, Lenalidomide and dexamethasone at first relapse (Kydar study) [13, 33]. These data suggest that 16p13.3 and 10p15.3 deletions can contribute to the deregulation of distinct high-risk pathways associated with resistance to Dara-Rd treatment.

Thanks to WGS resolution, we were able to investigate the clinical impact of structural variants (SVs) and complex events in patients treated with Dara-Rd (Supplemental Fig. 1B) [35]. SVs causing MYC gain of function were associated with shorter EFS, similar to what was observed in NDMM treated with Dara-KRd (i.e., NSMCE2 deletions, MYC and PVT1 events; P = 0.00016; log-rank test; Fig. 2H) [33]. Among complex events [18, 35], the presence of chromothripsis, a well-recognized adverse prognostic factor marker in NDMM, and it was found in 36% (10/28) of patients, with 7 experiencing progression (P = 0.05; log-rank test; Fig. 2I) [35]. Genomic events associated with shorted EFS were not differentially distributed between early (EFS < 12 months) and late progressors. Despite the limitations in our sample size, we conducted a multivariate analysis using Cox proportional hazard model with the genomic drivers associated with poor outcomes, revealing that the mutational process associated with reactive oxygen species damage (SBS18) and SVs involving MYC were significantly independent of HRCA and ISS classification (P = 0.01217 and P = 0.00581, respectively; Supplemental Table 7 and Supplemental Methods). Overall, these data suggest that distinct genomic events and patterns of genomic complexity are associated with resistance to Dara-Rd treatment.

Clonal evolution post Dara-Rd treatment

To reconstruct the genomic evolution from baseline to progression after Dara-based treatment, we analyzed WGS data from four patients with longitudinal paired samples collected before treatment and at progression. In addition, we included two patients treated with daratumumab as a single agent or in combination with bortezomib and dexamethasone (Dara-Vd), with available WGS data before and after treatment. Within these six patients, we identified patterns of branching evolution in all cases, with 50% experiencing a complete clonal shift (Fig. 3A–F). The presence of clonal shift was independent of the patients’ EFS. Although we did not observe any recurrent genomic events, we did detect distinct and intriguing events that were undetectable at baseline but became dominant upon progression. In the patient AIRC104 the progression was driven by a clone with a previously undetectable chromothripsis event involving chromosome 4, but not CD38 (Fig. 3A). In AIRC114 and AIRC115 the progression was driven by clone carrying translocations between chromosome 4 and 17 involving IRF4 and IKZF3 genes, and a large deletion of TP53, respectively (Fig. 3C, D). In AIRC117 at progression, the dominant clone had a deletion of chromosome 4p involving CD38 undetectable at baseline (Fig. 3E) [21].

Fig. 3: Clonal evolution of patients with RRMM treated with daratumumab-based treatment.

Top of each panel: the phylogenetic tree with the trunk in violet and the branches in salmon, light blue, and light green. The y axis is the number of single nucleotide variants. At the bottom left side, the cancer cell fraction (CCF) scatter-plot per each cluster at baseline (BAS) and progression (PD). In the bottom right side, the bar plot with the mutational signature relative contribution for each cluster. The mutational signature color legend is the same of the Fig. 2. HY hyperdiploid.

Association between immune profiles at baseline and clinical outcomes

To explore the immune composition and its impact on clinical outcomes, BM and PB samples collected at baseline from 32 patients were investigated by flow cytometry (Supplemental Methods; Supplemental Table 8). In line with what was previously shown, both the relative and absolute number of immune cells were highly correlated between BM and PB, with 77% of cell types having Spearman coefficient equal or greater than 0.50 (P < 0.001; Supplemental Fig. 3 and Supplemental Table 9) [37, 38]. Next, we investigated the association between distinct immune patterns and MM response to Dara-Rd. First, we characterized the immune cell landscape with particular attention to CD38 pos immune cells and pattern of immune exhaustion. This analysis revealed that most myeloid-derived suppressor cells (MDSC, CD33 + CD11b + ), monocytes (CD14 + ), and natural killer (NK, CD56 + ) cells were CD38 pos in both BM and PB samples (Supplemental Fig. 4A). Noticeably lower counts and proportions of CD38 pos cells were observed in conventional CD4 + T-helper cells and CD8 + cytotoxic T cells, Vγ9 Vδ2 T cells, B-reg (CD38hi CD19 + CD24 + ) [10], and T-reg (identified by CD4 + CD25 + + CD127lo and CD4 + CD25 + + FOXP3 + phenotypes) (Supplemental Fig. 4A). Despite monocytic-MDSCs (Mo-MDSCs, CD33 + CD14 + CD11b + ) exhibited a higher CD38 pos expression compared to granulocytic-MDSCs (Gr-MDSCs, CD33 + CD14 - CD11b + ) at the baseline in both BM and PB (Supplemental Fig. 4B), they were less prevalent than Gr-MDSCs (Supplemental Fig. 4C). Furthermore, at baseline all patients showed high prevalence of CD38 pos classical (CD14 + + CD16 -) and intermediate (CD14 + + CD16 + ) monocytes, and the classical subgroup was the most predominant compared with intermediate and nonclassical monocytes (Supplemental Fig. 4D, E). In addition, at baseline in BM samples the fraction of CD38 pos was higher among cytotoxic NK cells (CD56 + CD16 + ) than proliferative NK cells (CD56 + CD16 -) (Supplemental Fig. 4F). Across the whole BM immune populations, the proportion of cytotoxic NK cells was less abundant than the proliferative NK cells (Supplemental Fig. 4G).

Next, we went to investigate the impact of the BM and PB immune composition on clinical outcomes. At baseline, patients who experienced progression after Dara-Rd treatment showed an enrichment of CD38 pos NK cells (P = 0.04 in the BM and P = 0.05 in the PB estimated by Wilcoxon test; Fig. 4A and Supplemental Fig. 5A, B). This was also confirmed by measuring the mean fluorescence intensity (MFI); the CD38 expression on NK cells was higher at baseline in progressors compared with durable responders (P = 0.024 in the BM estimated by Wilcoxon test; Supplemental Fig. 5C). Furthermore, at baseline, progressor after Dara-Rd exhibited an enrichment of CD38 pos Mo-MDSC and Gr-MDSC cells compared to durable responders (P = 0.033 and P = 0.01 estimated by Wilcoxon test, respectively in the BM only; Supplemental Fig. 5A). Interestingly, at baseline, patients who progressed were enriched in higher levels of exhausted T cells, such as TIM3 + cytotoxic T cells and TIM3 + helper T cells (P = 0.0089 and P = 0.0031 respectively; Supplemental Fig. 5D–G) in PB, and lower presence of helper T cells (P = 0.015; Supplemental Fig. 5H, I) in BM. Moreover, at baseline, a significant enrichment of classical and intermediate monocytes was identified in progressed patient PB samples (P = 0.039 both; Supplemental Fig. 5J, K). These data suggest Dara-Rd efficacy in patients with RRMM, is influenced by distinct patterns of immune composition at baseline, particularly high CD38 pos immune cells, exhausted T cells, and reduced helper T cells in the BM and/or PB.

Fig. 4: Daratumumab immunomodulation on natural killer cells.

A Boxplot showing the percentage of CD38 pos among NK cells (CD56 pos). Light blue indicates durable responders, orange the progressors. P value was calculated with Wilcoxon test two-sided. B Heatmap showing CD38 pos in CD56 pos NK cells z-scores per each patient at seven different time points. Dots and black lines indicate that the sample was not collected at that time point, as the patient has previously progressed; one patient (AIRC124) was excluded because several samples were not collected at different time points. C Boxplots displaying the percentage of CD38 pos among CD56 pos NK cells over time. The red line connects the median value at each time point. D Boxplots showing the absolute number of cytotoxic NK cells over time. E Smooth plot presenting different NK cells patterns over time. F Boxplots showing the number of proliferative NK cells over time. C–F The P values were estimated using Wilcoxon test two-sided. The full list of all P values is reported in Supplemental Tables 11 and 12.

Immune modulation induced by Dara-Rd treatment over time

Given the concordance observed between BM and PB samples, we longitudinally assessed how Dara-Rd treatment influenced the immune system’s composition over time. This involved analyzing 202 PB samples collected at various time points from 31 study-enrolled patients. (Supplemental Table 8). One sample was excluded (AIRC124) due to lack of longitudinal samples after Dara-Rd treatment. As expected, both CD38 pos tumor and normal cells were rapidly depleted after exposure to daratumumab (Supplemental Fig. 6A and Supplemental Table 10) [15, 33, 39]. Among the main immune cell populations, a significant depletion was observed among cytotoxic and CD38 pos NK cells immediately after Dara-Rd treatment (P < 0.0001 estimated by Wilcoxon test at 3 month; Fig. 4B–D). In contrast, we observed an expansion of proliferative NK cells (Fig. 4E and Supplemental Table 11). This expansion was driven by CD38 neg proliferative NK cells that were not depleted by the daratumumab. Interestingly, the expansion of proliferative NK cells was significantly higher in durable responders compared to progressors (P = 0.015 estimated by Wilcoxon test; Fig. 4F and Supplemental Table 12). Durable responders were also characterized by a significant expansion of cytotoxic NK cells at 6 months (P = 0.043 estimated by Wilcoxon test, Fig. 4D and Supplemental Table 12). Consistent with earlier research, this data underscores the critical role of NK cells in the efficacy of anti-CD38 MoAb-based treatments [15, 33, 39]. Furthermore, these findings suggest a scenario in which daratumumab facilitates the fratricide depletion of CD38 pos NK cells, potentially constraining the subsequent availability of immune effector cells in the fight against myeloma cells. Because our flow cytometry panel utilized a CD38 multiepitope monoclonal antibody, the reduction in CD38-positive NK cells over time is likely due to depletion rather than CD38 epitope masking (see Supplemental Methods). This antibody effectively detected CD38, even when bound to daratumumab, using a distinct binding site, thereby avoiding false negative identification of CD38-masked NK cells. To further validate our observations, we analyzed CITEseq data from 20 NDMM patients undergoing Dara-KRd treatment [33]. We assessed NK cell count, CD38 expression, and protein levels in NK cells pre and post dara-based therapy. Our analysis revealed a significant reduction in NK cell count following Dara-KRd treatment (P = 0.0068; Supplemental Fig. 6B). Moreover, CD38 mRNA levels in NK cells exhibited a decrease post-treatment (P = 0.002; Supplemental Fig. 6C), which corresponded to a decrease in protein concentration (P = 0.0266; Supplemental Fig. 6D). Overall, these findings suggest that NK cells after Dara-based therapy exhibit decreased CD38 expression, rather than CD38 presence with daratumumab attached.

Next, we investigated the immune dynamics of monocyte cells and the classical (CD14 + + CD16 -), intermediate (CD14 + + CD16 + ), and nonclassical (CD14 + CD16 + +) subgroups. The longitudinal analysis showed a significant reduction of CD38 pos monocytes at 3 months after Dara-Rd treatment (P = 0.0001 estimated by Wilcoxon test), and the reduction was particularly relevant among the classical and intermediate subgroups (P = 0.0005 and P = 0.03 estimated by Wilcoxon test, respectively; Supplemental Fig. 7A and Supplemental Table 11). Interestingly, after being depleted at 3 months, the CD38 pos intermediate monocytes experienced a significant expansion after 12 months of therapy (P = 0.0042 estimated by Wilcoxon test; Supplemental Fig. 7A and Supplemental Table 11). In addition, independently from the CD38 status, reduction of intermediate and nonclassical monocyte subgroups was observed at 3 months after Dara-Rd treatment (P = 0.0036 and P = 0.0001 estimated by Wilcoxon test; Supplemental Fig. 7B, C and Supplemental Table 11).

Finally, we investigated the composition and clinical implications of T-reg and of patterns of T-cell exhaustion. While the proportion of T-reg (i.e., CD4 + CD25 + + FOXP3 + ) over both the total number of immune cells and the T-CD4 experienced a significant expansion over time, in particular after 9 months across all patients (Fig. 5A and Supplemental Table 11), the CD38 pos T-reg experienced a significant and sustained depletion immediately after the first cycles of Dara-Rd (P < 0.0001 estimated by Wilcoxon test, Fig. 5A, Supplemental Fig. 7D, and Supplemental Table 11). While the proportion of CD38 pos T-reg was not significantly different between durable responders and progressors, the latter group had a higher proportion of all T-reg (i.e., CD38 pos and CD38 neg) over the total number of T-CD4 at 3 months compared to durable responders (P = 0.05 estimated by Wilcoxon test; Fig. 5B; Supplemental Table 12). In terms of T-cell exhaustion, we found an increase in LAG3 + T cells following Dara-Rd treatment, particularly LAG3+ helper T cells (P = 0.014 estimated by Wilcoxon test; Fig. 5C; Supplemental Table 12). Notably, this increase was observed consistently over time. Furthermore, progressors showed higher expansion of LAG3 + helper T cells at 15 months compared to durable responders (P = 0.013 estimated by Wilcoxon test; Fig. 5C; Supplemental Table 12). In addition, progressors had a higher proportion of Gr-MDSC at 15 months (P = 0.0025 estimated by Wilcoxon test; Fig. 5B; Supplemental Table 12). Conversely, the durable responders showed enrichment of helper T cells (P = 0.046 at 6 months) and Vγ9 Vδ2 T cells at different time points (i.e., P = 0.04