Resolving therapy resistance mechanisms in multiple myeloma by multiomics subclone analysis

Intratumor heterogeneity arising from genetically distinct subclones as well as their adaptation and evolution during treatment is the main source of therapy resistance.¹ A paradigm disease for this process is multiple myeloma (MM), which is characterized by the clonal expansion of malignant plasma cells (PCs) in the bone marrow (BM).^2-5 MM intratumor heterogeneity as a clinical challenge becomes most evident during the course of several lines of treatment, when various resistance mechanisms accumulate.⁶ With the recent introduction of novel targeted immunotherapies⁷ for heavily pretreated relapsed or refractory MM (RRMM), the identification of common vulnerabilities of distinct multidrug-resistant subclones emerges as an urgent unmet need in the field. Toward this goal, MM subclones have been defined using single-cell RNA-sequencing (scRNA-seq) through the analysis of copy number aberrations (CNAs)⁸ and transcriptional clusters.⁹ Furthermore, mitochondrial DNA (mtDNA) mutations called from sc chromatin accessibility (scATAC-seq) data have been exploited to dissect and trace subclones in other hematologic malignancies.^10-12 At the same time, genomic,^13-17 transcriptomic,^9,18-20 or epigenetic¹⁹ sequencing readouts revealed molecular mechanisms that govern the treatment response in MM. Moreover, scRNA-seq methods provided information on the interplay of MM cells with their BM microenvironment (BME), which have been shown to be an essential part of the disease phenotype.^6,8,21 However, despite these advancements, it remains challenging to address the multifactorial nature of therapy resistance. Recent studies defined subclones based on a single readout and, therefore, focused only on individual aspects of molecular resistance. As such, it is not clear to what extend different resistance mechanisms can occur in combination in the same subclone and change during treatment.

In this study, we profiled serially collected BM samples from patients with RRMM using a multiomics approach that combines scRNA-seq, scATAC-seq, and bulk whole-genome sequencing (WGS) as well as bulk RNA-seq data. We show how copy number information from WGS and both sc modalities together with mtDNA mutations can be integrated to define subclones and their ancestral relationships at high resolution. Based on our matched transcriptome and epigenome analysis of individual subclones, we simultaneously address various layers of resistance, and specifically, we relate resistance to (1) subclones with preexisting epigenetic profiles that are associated with survival advantages, (2) converging phenotypic adaptation of genetically distinct subclones, and (3) subclone-specific interactions of MM and BME cells. From our analysis, it is apparent that multifactorial subclonal resistance mechanisms can occur in parallel, which has a number of implications for the development of strategies that target shared vulnerabilities across individual subclones.

Patient samples and related clinical information from both male (n = 7) and female (n = 8) patients with RRMM were obtained after written informed consent and in accordance with the Declaration of Helsinki. The study was approved by the ethics committee of the Medical Faculty at the University of Heidelberg (#S206/2011v13).

Viably frozen CD138⁺-enriched or CD138⁺-depleted BM mononuclear cells were thawed at 37°C, resuspended in RPMI with 10% fetal bovine serum, and washed in cold phosphate-buffered saline, with cells being collected via centrifugation at 500g for 5 minutes. scRNA-seq was performed according to the Chromium Single Cell 3ʹ Reagent Kit version 2 user guide (10x Genomics). Generated gene expression libraries were pooled and paired-end sequenced (26 bp and 74 bp) on an Illumina NovaSeq 6000 S2. For scATAC-seq, cell pellets were carefully resuspended in an ice-cold NP-40 lysis buffer and spun down immediately. Nuclei were resuspended in the nuclei buffer provided by 10x Genomics, counted and subjected to Tn5 tagmentation. The subsequent steps were performed according to the manufacturer's instructions for the 10x Genomics Single Cell ATAC version 1.0. Generated libraries were pooled and paired-end sequenced (26 bp and 74 bp) on an Illumina NovaSeq 6000 S2.

Preprocessing was performed using the DKFZ OTP WGS pipeline.²² Calling of CNAs and single nucleotide variants (SNVs) was performed as previously described.²³ For the supervised mutational signature analysis, we applied mmsig.²⁴ In order to call SNVs in the mitochondrial genome, we established a bulk sequencing analysis workflow.

Preprocessing, data normalization, and CNA prediction were performed as previously described.⁸ Differential and pseudobulk gene expression analysis was performed with Find(All)Markers and AverageExpression in Seurat.²⁵ Intercellular communication networks between tumor and BME cells were quantitatively inferred and analyzed with CellChat.²⁶ 

Preprocessing was performed with CellRanger ATAC and downstream analysis with ArchR.²⁷ All cells with a total number of <3000 fragments, a transcription start site enrichment score <7, a doublet enrichment score >6, and predicted doublet score >100 were excluded (supplemental Figure 1, available on the Blood website). PCs were distinguished from BME cells based on inferred gene-activity scores of specific marker genes, including TNFRSF17/CD138/SLAMF7/CD38, using MAGIC.²⁸ Nonmalignant PCs from different patients clustered together (supplemental Figure 1F). Myeloma cells from all patients were normalized, dimensionality reduced using iterative latent semantic indexing,^29,30 and clustered. Peak calling was performed with MACS2.³¹ Motif deviations were calculated based on the JASPAR database.^32,33 For each transcription factor (TF) motif, a mean value was calculated over all cells per sample. Coaccessibility among genomic regions was separately calculated for time points and subclones, adjusting the ArchR-framework to sc resolution without the aggregation of cells, as described previously.³⁴ We predicted CNAs by adopting an approach suggested by Lareau et al.¹⁰ Mitochondrial variants were called with mgatk.¹⁰ 

Subclones were identified separately in scRNA-seq and scATAC-seq data via supervised iterative hierarchical WGS-guided clustering of specific chromosomes or chromosomal regions using custom R-scripts (supplemental Figure 1K; supplemental Information; https://github.com/a-poos/MM_subclones).

For further details, see supplemental Information.

Previous studies have highlighted treatment-related subclonal dynamics in MM.¹³ However, how individual subclones respond to treatment and interact with the BME remains largely unknown, although this could provide crucial information of therapy resistance mechanisms. Accordingly, we performed droplet-based scATAC-seq and scRNA-seq in combination with WGS for CD138⁺-sorted BM aspirates of 15 patients with RRMM. Serial samples were collected for 12 of 15 patients before salvage therapy (T1) and at the time of subsequent relapse (T2; Figure 1A; supplemental Table 1), yielding 44 637 and 37 280 tumor cells after quality control, for scATAC-seq and scRNA-seq, respectively (Figure 1B; supplemental Figure 1; supplemental Tables 1 and 2). With clonal and subclonal CNAs occurring in virtually all patients with MM,³⁶ we implemented a WGS-guided clustering strategy to resolve the subclonal structure in scRNA-seq and scATAC-seq data (Figure 1C; supplemental Figures 1M and 2; supplemental Table 3). The results obtained for both sc-seq modalities were highly correlated (ρ = 0.972, linear model; average deviation, ∼5), enabling us to combine transcriptional and epigenetic profiles of individual subclones (Figure 1D; supplemental Figure 2Q).

Our initial approach assumed that MM cells with the same CNA profile belong to the same subclone. Next, we investigated whether mtDNA mutations further dissected CNA-based subclone definition in MM, because they have recently been proposed as an alternative approach for studying clonal dynamics in patients with chronic lymphocytic leukemia^10,11 and those with acute myeloid leukemia.¹² We called mtDNA mutations in the scATAC-seq data and confirmed the abundance of the respective variants in matched WGS (ρ = 0.98, linear model; supplemental Figure 3C-F). Overall, each patient with RRMM displayed a largely distinct mtDNA mutational profile (mean = 6 mutations; range 1-14 mutations; Figure 2A; supplemental Table 4). In line with previous reports, mutations followed a pattern equivalent to replication asymmetry and selective neutrality with C>T and T>C mutations^11,37,38 (supplemental Figure 3D).

For patient RRMM19, in line with WGS-based cancer clonal fractions for somatic nuclear mutations, 2 mtDNA mutations (15777G>A and 3749T>C) were distinctly enriched in single cells at T1 and T2, respectively. Thus, 2 subclones with the same CNA profile can be distinguished based on their mtDNA mutations (supplemental Figure 3A-B). We next compared subclones defined based on mtDNA mutations alone with those from our initial CNA-based assessment. Unique mtDNA mutations were detectable in 33% of the CNA-defined subclones (18 of 55; supplemental Figure 3G-L; supplemental Table 5), demonstrating that mtDNA mutations alone are not enough to resolve the subclonal architecture in RRMM.

By combining CNAs and mtDNA mutations, we were able to reveal further details of the subclonal architecture. Patient RRMM15 is the prime example for this finding (Figure 2B-G). In this patient, we initially identified 9 subclones using only the WGS-guided CNA approach (supplemental Figure 2P). By incorporating mtDNA mutation information, we identified 2 additional subclones in the scATAC-seq data (Figure 2B-D). These 2 subclones could be distinguished by the 2 lineage-defining mtDNA mutations, 10114T>C and 14769A>G (branch A) in scATAC-seq, and depicted unique chromosomal aberrations (supplemental Figure 3M). Importantly, these initially undetected CNAs could also be assigned to corresponding subclones in the scRNA-seq data (supplemental Figure 3N). In addition, we detected mtDNA mutations jointly enriched in multiple subclones that as such marked subclonal branches (eg, 15928G>A in subclones 1 and 2, 15356G>A and 12962G>A in subclones 7 and 8, and 13958G>A in subclones 9, 10, and 11; Figure 2B-C,E), which allowed us to delineate clonal relationships (Figure 2B-C,E; supplemental Figure 3G-L; supplemental Table 5). Cells with these jointly enriched mtDNA mutations segregated on the UMAP plot and marked distinct populations, including the branch-specific activity of POU5F1B, TCF3, and IRF4 TF motifs (Figure 2F-G), supporting the close relationship of those subclones.

Next, we assessed whether emerging relapse disease can be traced via mtDNA mutations. For this, sequencing data of tumor samples from 6 patients consecutively sampled over a time span of up to 3899 days were analyzed (Figure 2H; supplemental Figure 4A-C). Dynamic mtDNA mutations, including RRMM-enriched mutations, were detected in 4 of 6 patients (RRMM2/9/10/15; Figure 2H; supplemental Figure 4C). Remarkably, in these patients, RRMM-enriched mutations were seen up to 2464 days before T2 (Figure 2H). However, none of them were detectable at baseline (Figure 2H), suggesting clonal expansion of single tumor cells upon exposure to melphalan during first-line or relapse therapy (supplemental Table 1). Indeed, in our study, all patients with RRMM-enriched mtDNA mutations carried the single-base substitution MM1 mutational signature (supplemental Figure 4D), which is strongly linked with exposure to melphalan.^39-41 

Together, our results describe how mtDNA mutations can be best applied as lineage markers when combined with a CNA-based subclone assignment. This integrated approach yielded a comprehensive picture of the subclonal architecture with an average of 4 (range 1-11) subclones per patient (supplemental Figure 2Q) and matched transcriptome and epigenome profiles. In patients with paired samples, stable evolution, clonal selection, and branching evolution were seen in 5, 4, and 3 patients, respectively (supplemental Figure 2Q; supplemental Figure 5). It provides the foundation to investigate subclonal dynamics before and after treatment, as described in subsequent sections.

Having determined subclones with matching transcriptional and epigenetic profiles, we investigated the impact of treatment on individual subclones in all patients with multiple identical subclones that were present at both time points and defined based on our combined CNA- and mtDNA-based assignment (supplemental Figure 2Q). In patients with branching evolution, even at the sc level, we were unable to detect newly emerging subclones already at T1. Therefore, we focused on patients with stable evolution or clonal selection patterns to compare subclonal transcriptional and epigenetic changes over time (supplemental Figure 2Q). Although heterogeneous treatment regimens were applied (supplemental Table 1), a common pattern emerged in these patients: we found that the majority of therapy-induced changes in gene expression (median, 86%; range 64%-96%) and TF motif activity (median, 93%; range, 71%-95%) were shared across individual subclones for each patient (Figure 3A-G; supplemental Figure 6; supplemental Table 6). These converging adaptation mechanisms can be illustrated by known resistance mechanisms to 2 treatment approaches: proteasome inhibition via increased expression of heat-shock proteins (HSPs)^42-44 (Figure 3B-D) and MEK/BRAF inhibition via NF-κB pathway⁴⁵ upregulation (Figure 3E-G). In patient RRMM7, several HSPs known as resistance drivers, such as HSP90AA1,^42-44 were significantly upregulated upon carfilzomib treatment in both subclones (5 of 25 genes; P_adj ≤ .05; 1.5-fold enrichment both ways; Wilcoxon rank sum test; Figure 3C). In line with this finding, the motifs of TFs, including ZNF384, IRF1, and MEF2 family members, which can bind to the promoter of the upregulated HSPs,⁴⁶ displayed a significantly higher activity after therapy in both subclones (P_adj ≤ .05, Wilcoxon rank sum test; Figure 3D). In RRMM3, who was treated with MEK/BRAF inhibitors, both subclones demonstrated a significant upregulation of genes of the NF-κB pathway (18 of 37 upregulated genes; P_adj ≤ .05; 1.5-fold enrichment; Figure 3F).^47-50 This was corroborated by a significantly higher activity of NF-κB motifs upon MEK/BRAF inhibition in both subclones (P_adj ≤ .05; Figure 3G). Combined BRAF/MEK inhibition has previously been shown to induce high rates of clinical response.^23,51 Accordingly, we analyzed an additional patient belonging to this treatment group, who exhibited major changes in its clonal composition (RRMM9; Figure 3H). Similar to RRMM3, 11 of 37 (30%) genes upregulated upon MEK/BRAF inhibition are part of the NF-κB pathway⁴⁷ (P_adj ≤ .05; log2 fold–change > 0.25; Wilcoxon rank sum test; Figure 3I). These genes included BCMA (TNFRSF17),^51,CD79A,^48,49 and BTG2.⁵⁰ For RRMM9, only 1 of 4 subclones was still detectable after therapy (Figure 3H). This subclone depicted the highest NF-κB–family member motif activity of all subclones at T1 (Figure 3J). We conclude that, independent of subclonal composition, treatment resistance can be mediated either by converging phenotypic adaptation of all subclones or preexisting epigenetic profiles that do not necessarily have to be therapy-induced. This emphasizes the importance to assess MM on the level of individual subclones as well as the use of sc multiomics analysis for our understanding of the role of transcriptional and epigenetic features in treatment-induced resistance.

Based on the observation that subclonal properties can underlie differential treatment response, we next considered other patients in whom treatment-induced changes in the subclonal architecture were seen (supplemental Figure 2Q). We found that MCL-1 inhibition in patient RRMM6 resulted in the depletion of a subclone with a biallelic TP53 aberration (Figure 4A; supplemental Figures 4A and 5C). To address the molecular mechanism underlying this depletion, we first explored potential subclonal alterations of the target. We detected differences in the coaccessibility at the MCL-1 locus between the 2 main subclones (subclone 1 vs subclone 2) before initiation of MCL1 therapy. However, these differences did not result in significant MCL-1 expression changes (supplemental Figure 7). Next, we assessed whether TP53 aberrations directly affect the vulnerability of cells to MCL-1 induced cytotoxicity. Therefore, genetically engineered biallelic TP53 deleted + R175H mutated (del/mut)^52,53 and wild-type (wt/wt) AMO-1 cell lines were exposed to different doses of an MCL-1 inhibitor in vitro (Figure 4B). In contrast to what is suggested by the findings for RRMM6, the biallelic TP53 AMO-1 cell line was less sensitive compared with the wt cell line (P < .05, Welch two-sample t test, Figure 4B; supplemental Figure 7B). Thus, inactivation of TP53 does not appear to underlie the depletion of subclone 1 in RRMM6. We, next, turned toward transcriptional differences between the 2 main subclones. The top differentially expressed genes between the 2 main subclones at T1 covered several surface markers (Figure 4C; supplemental Table 6). Among those, CD44 has been reported as a key cell adhesion–mediated drug resistance (CAM-DR)⁵⁴ molecule and is a known mediator of tumor-BME interactions.^55,56 In line with the differential gene expression analysis, a coaccessibility analysis of the scATAC-seq data showed a strong correlation between the CD44 promoter and a distal putative enhancer element for subclone 1 that was absent for subclone 2 at time point T1 (Figure 4D). Thus, activating an enhancer in one but not the other subclone could lead to the observed differences in CD44 expression. To test whether a higher CD44 expression is related to biallelic TP53 inactivation, we performed bulk RNA-seq in a set of 46 patients with RRMM.^23,57 We found CD44 expression significantly increased in patients with TP53 biallelic mutation compared with that in patients with wt TP53 (P = .012, two-sided Wilcoxon rank sum test) as well as mutated TP53 (P = .051; Figure 4E). As a functional validation, we turned toward publicly available bulk RNA-seq data^52,53 from genetically engineered biallelic TP53 del/mut (R175H) and wt/wt AMO-1 cell lines. We found a significantly higher CD44 gene expression in the cell line with inactivated TP53 than in wt TP53 (Figure 4F; adjusted P = 7.103e−07; log₂ fold–change = 0.524, Wald test after Benjamini-Hochberg correction for multiple testing). This finding could also be confirmed using western blot analysis (Figure 4F). The high CD44 expression at T1 subsequently decreased upon MCL1-inhibition, suggesting a potential mechanism involving altered subclone-BM microenvironment interactions (Figure 4G). Indeed, CellChat²⁶ predicted a specific interaction between subclone 1 and cDC2/GMPs/CD14⁺ as well as CD16⁺ monocytes via CD44 at T1 (Figure 4H; supplemental Table 7). Importantly, this respective interaction was lost upon MCL1-inhibition (Figure 4H). Together these results suggest differences in tumor-BME interactions among subclones with differential treatment response.

The results described above prompted us to broaden our interaction predictions to all patients with at least 2 coexisting subclones per time point and matched scRNA-seq of the BME (n = 7). Across these patients, CellChat predicted, on average, 32 ligand-receptor interactions (range 23-45) between tumor and BME cells, of which ∼20% (median = 7; range 3-12) were not predicted in all subclones per patient (unshared interaction, Figure 5A-B; supplemental Table 7). The intercellular adhesion molecule (ICAM) signaling network was found to mediate most of the unshared interactions with the BM microenvironment across our patient cohort (Figure 5C-D; supplemental Figure 8A-B). In particular, the interaction between ICAM2 on the myeloma cells and several integrins (ITGB2/ITGAL/ITGAM/ITGAX) on the BME cells was predicted to be subclone-specific in 5 patients.

All patients with longitudinal samples showed interactions that changed between time points (temporal changes; Figure 5B-C; supplemental Table 7). For example, in patient RRMM15, all subclones at T2 showed a temporal interaction mediated by the ICAM1 pathway, which was not detected at T1 (Figure 5C; supplemental Figure 8B). ICAM1-targeted immunotherapies in the form of anti-ICAM1 monoclonal antibody bersanlimab (BI-505, BioInvent)⁵⁸ or CD38 × ICAM1 bispecific antibody⁵⁹ are already being tested in clinical trials (NCT01838369) or have been suggested for MM and B-cell non-Hodgkin lymphoma.^60-62 Accordingly, we further investigated ICAM1 temporal expression changes in an extended bulk RNA-seq data set of 16 patients with RRMM. We found that ICAM1 expression was significantly upregulated at T2 (P = .02, paired Wilcoxon signed rank test), particularly in patients RRMM6 and RRMM15 (Figure 5E). To further delineate the mechanism underlying increased ICAM1 expression, we assessed the chromatin accessibility at the ICAM1 locus through our scATAC-seq data. We observed a higher coaccessibility at the ICAM1 promoter region for cells at T2 than those at T1, suggesting a therapy-induced higher activity of ICAM1 (Figure 5F; supplemental Figure 8C). Moreover, RRMM cells also showed higher accessibility of TF motifs previously described in MM pathogenesis, including IRF,^63-72 NF-κB,^73-75 and STAT^76-78 family members at T2 (Figure 5G; supplemental Figure 8D). Within these families, IRF4, IRF1, RELA, STAT1, and STAT2 are known TFs that drive ICAM1 expression.⁴⁶ This upregulation mechanisms derived from our scATAC-seq data was in agreement with chromatin immunoprecipitation assays followed by sequencing (ChIP-seq) analysis of the MM cell line KMS12BM.⁷² We found an overlap of the ICAM1 promoter peak with IRF4 ChIP-seq peaks (Figure 5F; supplemental Figure 8C). These orthogonal data sets suggest a differential ICAM1 expression driven by TFs that include IRF4.

In summary, we observed subclone-specific interactions with BME cells that change upon treatment in all patients. Some interactions, such as those mediated by ICAM1, were shared among subclones at a given time point and, thus, represent potentially actionable immunotherapy targets across subclones.

The assessment of MM subclones and their molecular characterization during treatment are key for advancing patient stratification and therapy response prediction. In this study, we implemented an integrative multiomics approach that combines information from different bulk and sc readouts in a longitudinal analysis of RRMM to determine subclones at unprecedented resolution. Firstly, using WGS as a reference largely improved the accuracy of the CNA-based subclone definition from sc-seq data. Secondly, integrating 2 sc-seq approaches enabled us to further cross-validate our WGS-guided subclone definition and to analyze subclone-specific transcriptome and epigenome plasticity. Thirdly, including mtDNA mutations called from scATAC-seq resolved subclones with identical CNA profile. Our framework increased the ability to capture the subclonal diversity and their ancestry. For example, we identified 11 subclones for patient RRMM15 when using all readouts, although only 4 subclones were detected by clustering based on CNAs of all chromosomal regions. This was because CNA subclones from the same ancestral branch clustered together, suggesting that defining subclones based on transcriptional and epigenetic clusters alone may lead to an underestimation of the subclonal diversity. In addition, we conclude that mtDNA mutations on their own are insufficient for MM subclone definition, with only one-third of CNA-defined subclones in this study displaying unique mtDNA mutations. This observation could also be relevant for other hematologic malignancies, including chronic lymphocytic leukemia ^10,11 and acute myeloid leukemia,¹² in which mtDNA mutations have been used to define subclones.

Our approach provided a high-resolution map of the dynamic subclonal architecture with matched RNA and ATAC profiles for individual subclones before and after therapeutic interventions. We find that relapse can occur via different resistance mechanisms that include the expansion of subclones with preexisting epigenetic profiles associated with resistance, plasticity, and adaptation of gene expression programs in individual subclones as well as changes in MM-BME interactions. Plasticity in gene expression is in line with previously described resistance mechanisms to proteasome inhibitors via the upregulation of HSPs^42-44 and activation of the NF-κB pathway⁴⁵ upon MEK/BRAF inhibition. Thus, it will be important to monitor these pathways in relation to subclone composition for assessing treatment response.

In several patients, we found a branching evolution pattern as identified previously as the major evolution pattern in deep responding patients with MM.¹³ This pattern is characterized by depletion of previously dominant subclones and emergence of new subclones, which makes it challenging to study subclones over time in these patients. In this study, we were unable to detect the new subclones present at T2 at sc level at T1 in patients with branching evolution. To account for this issue, we envision our approach to be applied to tumor cells purified from patients receiving therapy, for example, as previously performed 48 hours after therapy start,⁷⁹ or from patients tested positive for minimal residual disease at various time points during the course of disease.⁸⁰ In patients with stable evolution and clonal selection, we were able to longitudinally track individual subclones over time. Future studies focusing on larger cohorts with uniform treatment will be crucial to understand common transcriptomic or epigenetic features among subclones that have disappeared or emerged.

Previous work has characterized the BME of different MM disease stages,^6,8,81 providing evidence for deregulation of cell surface marker expression as resistance mechanisms.^8,19 Through the combined analysis of tumor subclones and their BME, we show that individual subclones can have distinct interaction profiles with the BME, which further increases the complexity of MM heterogeneity. This includes subclone-specific expression and transcriptional plasticity for CD44, a key cell adhesion–mediated drug resistance⁵⁴ molecule and promising immunotherapy target.^55,56 In future studies, subclone-specific cellular interactions, such as CD44 on MM cells and LGALS9 on BME cells, could be validated by emerging spatially resolved techniques.^82,83 Thus, at least in some patients, subclones might display individual resistance mechanisms and respond differently to immunotherapy.^7,84,85 Consequently, to improve the outcome of available therapies, the identification of shared targets across subclones is essential and remains one of the biggest challenges to achieve long-term remissions. In this context, we show highly similar transcriptional and epigenetic changes upon treatment in patients with multiple distinct subclones, indicating converging adaptation mechanisms. The repertoire of shared treatment-induced changes included upregulation of ICAM1, another key mediator of cell-cell interactions.⁸⁶ As BME cells nurture and facilitate the proliferation of MM subclones,⁸⁷ such shared mediators could constitute actionable candidates for immunotherapeutic targeting of a heterogenous disease.² In summary, we anticipate that the integrative multiomics analysis, introduced in this article for RRMM, will prove useful for therapy response prediction based on the subclone composition as well as the identification of novel molecular targets against multidrug-resistant subclones.

The authors are grateful to the patients for participation in this study and thank the Sample Processing Lab, the High Throughput Sequencing unit of the Genomics & Proteomics Core Facility and the Omics IT and Data Management Core Facility of the German Cancer Research Center (DKFZ), the DKFZ–Heidelberg Center for Personalized Oncology (DKFZ-HIPO) office, the Biobank Multiple Myeloma UKHD and the Myeloma Registry for excellent services.

Support and funding of the project via the NCT Heidelberg Molecular Precision Oncology Program through projects HIPO H067, K08K and K43R (H.G., K.R., M.S.R.) is gratefully acknowledged. N.P. was supported by a European Union Horizon Europe Marie Skłodowska-Curie Postdoctoral Fellowship (101068158); S.H., H.G., M.S.R., the Biobank Multiple Myeloma UKHD, and the Myeloma Registry are supported by the Dietmar-Hopp Foundation; K.R. is supported by the German Research Foundation (DFG) grant for subproject Z1 within SFB1074. Data storage service via SDS@hd is supported by the Ministry of Science, Research and the Arts Baden-Württemberg and the DFG through grants INST 35/1314-1 FUGG and INST 35/1503-1 FUGG.

Contribution: A.M.P., N.P., M.S.R., K.R., and N.W. conceptualized the study; N.P., J.-P.M., S.M.T., K.B., A.B., and J.R. performed experiments; A.M.P., M.J.P., S.S., I.S., and N.W. analyzed data; A.M.P., N.P., M.J.P., K.R., and N.W. interpreted the data; N.P. drafted the figures; A.M.P. and N.P. wrote the original draft; I.S., M.J.P., M.S.R., K.R., and N.W. edited the manuscript; C.M.-T., H.G., M.S.R., and K.R. acquired funding; L.J., U.M., L.R., K.M.K., P.R., S.H., N.G., and M.S.R. acquired resources; O.S., M.S.R., K.R., and N.W. supervised the study; and all authors gave final approval.

Conflict-of-interest disclosure: The authors declare no competing financial interests.

Correspondence: Karsten Rippe, Division of Chromatin Networks, German Cancer Research Center and BioQuant, Heidelberg, Germany; e-mail: karsten.rippe@dkfz-heidelberg.de; and Niels Weinhold, Department of Internal Medicine V, University Hospital Heidelberg, Heidelberg, Germany; e-mail: niels.weinhold@med.uni-heidelberg.de.