Molecular Clusters and Tumor-Immune Drivers of IgM Monoclonal Gammopathies

IgM monoclonal gammopathy of undetermined significance (MGUS) and Waldenström macroglobulinemia (WM) are two clinicopathologic entities that fall under the umbrella of IgM monoclonal gammopathies. They represent a disease spectrum which extends from an asymptomatic condition with a mild increase in serum monoclonal IgM (IgM MGUS) to a symptomatic disease characterized by organ infiltration by a lymphoplasmacytic process and hypergammaglobulinemia causing organ disfunction and potential hyperviscosity (WM). The management of these diseases may be very different, ranging from observation to treatment with Bruton's tyrosine kinase (BTK) inhibition or chemoimmunotherapy (1, 2). Despite the recent scientific advances in the field, the current disease classification still relies solely on clinical features that do not explain the heterogeneity that exists within these conditions. This is in part due to the lack of a comprehensive evaluation of the mechanisms that drive the IgM monoclonal gammopathies (3, 4).

While previous studies have explored the pathobiology of this disease by investigating their genetic and transcriptional profiles, all of them employed a narrow focus, either by virtue of using predefined assays and/or a single-omics evaluation, or by limiting the population using the same arbitrary diagnostic criteria defined almost two decades ago (5–9). To shed light on the biology that may account for the clinical differences in IgM monoclonal gammopathies and overcome some of the constraints from previous studies, we performed the first comprehensive multi-omics study including whole-exome sequencing (WES), RNA sequencing (RNA-seq), proteomics, metabolomics, and mass cytometry (CyTOF) of clonal bone marrow (BM) CD19⁺ and/or CD138⁺ sorted cells, matched BM supernatant and peripheral blood (PB) serum from patients with IgM MGUS and WM to tackle this open question. Here, we identify three molecular clusters of patients with distinct transcriptional signatures, and these clusters demonstrate concordant proteomic, metabolomics, and immune features that have the potential to inform therapeutic decision-making.

A total of 32 patients with WM (n = 25) or IgM MGUS (n = 7; Supplementary Table S1), and S5 control subjects who underwent surgery for other reasons than cancer (e.g., hip replacement) were included in our study. IgM MGUS was defined as < 10% BM involvement by lymphoplasmacytic lymphoma (LPL), serum IgM monoclonal protein < 3 g/dL, and no signs or symptoms of disease. WM was defined as ≥ 10% BM involvement by LPL and signs or symptoms secondary to the lymphoplasmacytic infiltrative process (10). Only subjects with (i) BM clonal CD19⁺ and/or CD138⁺ sorted cells, (ii) BM supernatant, and (iii) PB serum collected at the same time point were included. Samples were collected at diagnosis in 27 patients and after one line of therapy in 5 patients with WM (median time from last therapy to sample collection 44 months, range 6–91 months), but not in a consecutive fashion (pre- and posttreatment). None of the patients with WM received first-line BTK inhibitor therapy because most of them were diagnosed before BTK inhibitor FDA approval in 2015, while a further 8 were treated with bendamustine plus rituximab after 2015 as per Mayo Clinic preference (11). The study was reviewed and approved by the human subjects Institutional Review Board at the Mayo Clinic, and written informed consent was obtained from all participants. The study was conducted according to the Declaration of Helsinki.

WES was performed on fresh-frozen CD138⁺ sorted BM cells and matched normal PB nucleotide cells from 27 patients (24 patients with BM infiltration and 3 MGUS without BM involvement) and executed as previously described by Hartert and colleagues (12). In summary, DNA was extracted using the QIAamp DNA mini kit (Qiagen GmbH, Hilden, Germany). WES was carried out at the Mayo Clinic Genome Analysis Core Facility using the Illumina HiSeq 4000 with an intended coverage of 80X. Library preparation was performed with the Low Input SureSelect kit. Prior to alignment to human reference, Agilent AGeNT LocatIt (version 4.0.1) was used to remove duplicate reads using unique molecular barcodes per fragment information. Non-duplicate reads were mapped to human genome reference build 38 (GRCh38) using BWA aligner version 0.7.10, followed by realignment and recalibration using GATK version 3.4–46 (13). Somatic mutations were called using consensus calls from three somatic callers, SomaticSniper (v1.0.4.2), JointSNVMix (v0.8-b2), and MuTect (v1.1.7; refs. 14–16), all of which use a tumor normal pair for variant calling. Mutations were annotated using an in-house annotation tool, BioR (17). Post annotation, mutations present at a frequency above one percent in public germline datasets ExAC and The 1000 Genome project were filtered out. In addition, any mutation with a depth of less than 10 supporting reads and/or with a variant allele frequency of less than one percent or with less than 2 alternate supporting reads were filtered out. Mutations in WM and lymphoma genes which have been previously published were analyzed (18–21). Generated segmentation files by PatternCNV (22) were analyzed using the GISTIC 2.0 tool (23) and these results were used to determine copy-number loss and gain on previous reported regions (9). For the 6q deletion, multiple single regions were lumped together. All WES data, including mutations and copy number were analyzed and visualized by using the R tools maftools (24) and complex heatmaps R program tool (25).

Total RNA was extracted using the Qiagen miRNeasy Micro Kit (217004). The library was prepared using the RNA access kit TruSeq Stranded Total RNA v2 from Illumina. The raw RNA-seq reads were processed through the Mayo Clinic RNA-seq bioinformatics pipeline, MAP-RSeq version 3.0, as previously described (26). Briefly, MAP-RSeq used the STAR tool (27) to align reads to the reference human genome build hg38. Gene and exon expression quantification were performed using the Subread (28) package to obtain raw counts. We further removed genes with a median count below 30. We next computed transcripts per million (TPM) for the expression data. For our gene expression clustering, we used the nonnegative matrix factorization (NMF) unsupervised clustering approach that ensure robustness of clustering (29). Only protein coding genes, with the exclusion of mitochondrial genes and sex chromosome genes, were used for the NMF and principal component analysis (PCA). In addition, transformed values to log2(TPM+1) were used for the clustering and the algorithm was run with the Brunet method (30) iterating 100 times to provide results from k = 2 until k = 6. Cluster number was selected on the basis of cophenetic value (Supplementary Fig. S1) and methods described by Brunet and colleagues (30). For differential expression analysis, the R package edgeR was used (31, 32). Only genes with FDR value < 0.05 were considered as statistically significant.

Gene set enrichment analysis was performed using the fgsea R tool (33). For the enrichment analysis all genes from the differently expression analysis were used and a pre-rank gene list was generated using the following formula: fold change * −log₁₀FDR. Genset considered in the analysis were the Hallmark, C2, CP:REACTOME, CP:BIOCART, CP:KEGG, CGP (chemical and genetic perturbations), CP:PID and CP:WIKIPATHWAYS. Gensets with a P value < 0.05 were considered significant. Pathways containing less than 40 and more than 500 genes were not tested.

BM supernatant from 32 patients (WM, n = 25; IgM MGUS, n = 7) and PB serum from16 patients (WM, n = 11; IgM MGUS, n = 5) were available for proteomic analysis using the SOMAscan Platform (SomaLogic Inc.) v3.1 and v4.1, respectively. SOMAscan is an aptamer-based assay allowing for the simultaneous measurement and quantification proteins through chemically modified nucleotides to transform a protein signal into a nucleotide signal that can be quantified by standard DNA technique using microarray. The assay used in this study measures more than 1,300 protein analytes quantifying proteins that span over 10 orders of magnitude (34). Raw protein intensities observed in each sample were adjusted using hybridization control probes. Transformed intensity distributions of samples were normalized using median centering. Normalized intensities of each protein were compared between groups using a generalized linear model configured to use Gaussian distribution and ANOVA contrasts. Resulting differential expression P values were adjusted using Benjamini–Hochberg method. Proteins have an adjusted P value < 0.05 were considered as significantly differentially expressed between groups.

A functional pathway analysis between the three cohorts identified by NMF was performed using the Ingenuity Pathway Analysis (IPA) tool (QIAGEN Inc.,) using the differentially expressed proteins. The IPA Core Analysis function was used to estimate canonical pathway dysregulation and upstream regulators based on protein expression differences between the cohorts using the Ingenuity Knowledge Base and tested by the Fisher exact test with a P value < 0.05. An activation Z-score (activation/inhibition) was generated on the basis of consistency between the experimental protein expression and the expected canonical pathways patterns considering the activation state of each molecule and the relationship to each with the combination of other two clusters.

PB serum of 25 patients (WM, n = 20; IgM MGUS, n = 5) and BM supernatant of 12 patients (WM, n = 12) were available for metabolomic analysis through Capillary Electrophoresis Time-of-Flight Mass Spectrometry (CE-TOFMS) and Liquid Chromatography Time-of-Flight Mass Spectrometry (LC-TOFMS) in two modes for cationic and anionic metabolites as described by Jalali and colleagues (35).

CyTOF was performed in 11 patients as previously described by Mondello and colleagues (36). In brief, all the samples were normalized and analyzed simultaneously to account for variability in signal across long acquisition times. Linage markers (Supplementary Table S2) were used to separate T cells into regulatory, naïve, follicular helper, memory, effector, or exhausted subsets, monocytes and macrophages. All the antibodies were purchased from Fluidigm. We used the SPADE method (37) to first perform hierarchical cluster analysis on the immune cell subsets.

Comparison of quantitative data between groups was done by Student t test or one-way ANOVA test (assuming normal distribution). Statistical analysis was done in JMP (Pro 14.1, SAS, Cary, NC). All reported P values were two-sided, and P values < 0.05 were considered statistically significant (*, P < 0.05; **, P < 0.01; ***, P < 0.001). The correlation between the processed intensities of each protein observed in the BM samples versus those in the PB samples was assessed using Spearman's rank correlation coefficient (38). The correlation between the overall protein profile of BM and PB was calculated using RV coefficient (39). Similar methods were applied to assess the relationship between PB and BM metabolites.

The data generated in this study are available upon request from the corresponding author.

To explore the molecular characteristics of IgM monoclonal gammopathies, we performed RNA-seq of sorted BM clonal CD19⁺ and/or CD138⁺ B cells from patients with IgM MGUS (n = 7), WM (n = 25), and healthy controls (n = 5). As the first step, we ran a PCA to provide insight into differences in gene expression profiling. While most WM samples clustered together, a subset of them grouped separately in close proximity to those from MGUS, suggesting the existence of variability among patients with WM. Conversely, MGUS and controls remained relatively localized to their respective groups except for one MGUS sample that grouped with WM (Fig. 1A). Of note, this MGUS patient presented with more advanced clinical features compared with the others as indicated by a greater extent of BM infiltration, higher IgM levels and mild cytopenia. To further explain the transcriptional relationship between these samples, we employed the unsupervised NMF plus consensus clustering for identifying a stable clustering solution between 2 and 6 clusters. We determined three clusters to be the most stable and robust result (Fig. 1B; Supplementary Fig. S1) based on established methods (30). Cluster 1 (C1) included only patients with WM, cluster 2 (C2) included patients with both WM and IgM MGUS, and cluster 3 (C3) included all controls as well as a small number of WM and IgM MGUS patients (Fig. 1C). Next, we investigated the transcriptional differences among these newly identified clusters. Consistent with our separation between WM and mixed MGUS/WM, C1, and C2 clustered apart. Similarly, C3 formed an intermediate group between C1 and C2 that was reminiscence of the control group. (Fig. 1D). Accordingly, these clusters showed distinct gene expression profiles (Fig. 1E; Supplementary Table S3), which in turn translated into different clinical features (Table 1). To further exclude potential confounding factors associated with residual nucleated cells after sorting, we investigated whether the expression signatures of immune cell subtypes (e.g., T cells, dendritic cells, macrophages) were differentially represented in these 3 clusters. As expected, there was no significant enrichment of immune cells, which might have altered our results (Supplementary Fig. S2).

To investigate whether genetic mutations and/or structural alterations could explain the differences among the clusters, we performed WES of the same matched patients analyzed by RNA-seq. On the basis of feasibility of DNA extraction, WES data were available in 14/16 (87%) patients in C1 and 7/11 (64%) in C2. For this analysis we excluded healthy controls from C3 and focused on the subgroup comprising only diseased patients (C3d, n = 5). Historically, somatic mutations of the MYD88 and CXCR4 genes are most frequently found in patients with WM (86%–97% and 24% respectively) compared with IgM MGUS (42%–54% and 7% respectively; refs. 28–30). Similarly, chromosome 6q deletion is the most common structural abnormality reported in WM (approximately 55% of the cases), but rarely described in IgM MGUS (31). In our study, the overall incidence of MYD88 L265P and CXCR4 mutations and 6q deletion were 78%, 26%, and 30%, respectively (Fig. 2). While MYD88 L265P was confirmed as the most common mutation in all patients, CXCR4 mutation and deletion of 6q showed a lower frequency rate compared with historical controls (7–9). Interestingly, there was no difference in mutation burden of MYD88 L265P and CXCR4 among the three clusters, likely due to the spread of patients with WM among the groups. Conversely, alterations in TNFAIP3 including mostly 6q deletion (which encodes for TNFAIP3) and one somatic mutation were significantly enriched in C1 compared with C2 and C3d (P = 0.019 and P = 0.04 respectively; Supplementary Table S4). TNFAIP3 copy-number loss was associated with aberrant expression of TNFAIP3 (Supplementary Fig. S3). In addition, deletion of 11p15.4 and deletion of 17q25.3 were more frequent in C2 compared with others (P = 0.004 and P = 0.011, respectively; Supplementary Table S5). No additional significant differences in genetic mutations or DNA abnormalities were found among the three clusters. For example, the incidence of 1p36.33 deletion [which involves TP73 (ref. 32), a tumor suppressor with close homology to the tumor suppressor gene TP53] was comparable among the three clusters (P = 0.661). Similarly, the frequency of MYC amplification was equivalent in all groups (P = 0.605; Supplementary Table S4–6). While these DNA abnormalities may play a role in the tumorigenesis of IgM monoclonal gammopathies, they do not provide a clear explanation to the different transcriptional signatures of these newly identified clusters.

Next, we took a deeper look at the molecular differences of each cluster compared with the others. RNA-seq revealed inverse transcriptional profiles between C1 and C2, with 1,948 and 1,942 differentially expressed genes respectively (FDR < 0.05) that were markedly skewed toward transcriptional suppression in C1 (n = 1,724 vs. n = 224; Fig. 3A,i) while toward transcriptional activation in C2 (n = 1,600 vs. n = 342; Fig. 3A,ii). C3d showed a transcriptional profile intermediate between C1 and C2 with 475 significantly upregulated and 175 significantly downregulated genes (Fig. 3A,iii). In line with this, we found that genes involved in cell cycle (e.g., EIF4E2, EIF2B1, ATM), cell proliferation (e.g., PIK3IP1, AKT3, MAP4K1, IKBKB), apoptosis (e.g., BCL2, BAX, PMAIP1), immune response (e.g., CD40, BATF, CD80, BTLA), antigen presentation (e.g., HLA-DRA, HLA-DOA, HLA-DQA1, HLA-DPB1), and metabolism (ENO3, ALDOC) were significantly upregulated in C1 while downregulated in C2 (Fig. 3B). Similarly, C3d showed significant increase of genes related to cell cycle (e.g., E2F2, CDK1, CHEK1, CDK6), cell division (e.g., AURKA, AURKB, PLK1, PLK4), and DNA repair (e.g., RAD51) as well as downregulation of cell proliferation (e.g., STAT3, TRAF3, IRAK2) and immune response genes (e.g., IL10RA, IL6R, BATF, IFNGR2; Fig. 3B). Gene set analysis revealed a negative enrichment in C1 for (i) signatures associated with stem cells and senescence; (ii) direct targets of cell cycle such as E2F and 4EBP1, which may promote tumor growth; while there was a positive enrichment for (iii) mRNA splicing and processing, which are major regulators in initiation of translation and posttranslational modification. The opposite was observed in C2. C3d confirmed intermediate features with combined upregulation of proliferation and stem cell functionality along with mRNA processing (Fig. 3C). Importantly, we observed significant negative enrichment of transcriptional signatures linked to immune response including antigen response and chemokine signaling in C1 (FDR = 0.002) and C3d (FDR = 0.001), while the same were positively enriched in C2 (FDR = 0.0004; Fig. 3C and D). Along the same lines, we observed a divergent metabolic reprograming in C1 and C2 with transcription activation of mitochondrial metabolism but suppression of glycolysis, lipid and carbohydrate metabolism in C1 while the opposite was true in C2. Like C1, C3d showed enrichment for mitochondrial metabolism, however no significant difference was observed for the other pathways (Fig. 3C).

Taken together, these data show that each of the molecular clusters has discrete and different transcriptional profiles. The signatures enriched in C1 suggest increased proliferation due to the combined suppression of cell cycle, stem cell, and senescence signatures, which may be further exacerbated by concomitant immune suppressive mechanisms. In contrast, C2 displays predominantly inflammatory features linked to upregulation of immune response and stemness signatures along with concordant suppression of genes involved in cell cycle. Finally, C3d presents an intermediate profile between C1 and C2, which may suggest a transitional stage between the two clusters (Supplementary Fig. S4).

To investigate whether the transcriptional signatures associated with the three molecular clusters translate into concordant proteomic abnormalities, we performed a SomaScan assay of BM supernatant obtained from patients with matched RNA-seq data. In accordance with the transcriptomic analysis, C1 showed a protein expression profile that was markedly skewed toward protein downregulation (n = 298 down vs. n = 65 up; Fig. 4Ai) while C2 toward protein upregulation (n = 88 up vs. n = 44 down; Fig. 4Aii). In line with its intermediate feature, C3d showed only few proteins differently expressed compared with the other two clusters (n = 4 vs. n = 5; Fig. 4Aiii). C1 was characterized by decreased expression of proteins involved in inflammation, including IL6, IFNγ, and TGFβ, along with increase of IL1R2, which promotes immunosuppression. As these cytokines play a critical rule in the immunosurveillance against tumor, their dysregulation highlights a potential immune escape mechanism linked to C1. In support of the putative immune evasion, we found lower concentrations of TNFRSF1A (regulator of immune response), IRF1 and STAT1 (essential transcription factors for MHC class I), FOXO1 (critical regulator of M2 macrophages) and of stimulatory immune receptors such as CD40, CD28, and ICOSL. There was also downregulation of proteins regulating the inflammatory milieu, including JAK2, STAT3, and BTK, as well as upregulation of SOCS1, which is negative regulator of JAK/STAT pathway. Notably, we found downregulation of CREBBP and EP300 with concomitant increase in BCL6, which is reminiscent of recent findings indicating that CREBBP loss-of-function induces malignant transformation in part by silencing of antigen presentation signatures (40). In addition, we confirmed downregulation of cell-cycle checkpoints including aurora kinase B, cyclin B1, CDKN1B, which in turn favor cell proliferation and growth. In contrast, the hallmark of C2 was upregulation of proteins involved in inflammation, with IL6, IL2, and IFNγ being the most dysregulated cytokines. Additional highly upregulated proteins associated with inflammation included TNFRSF18, CCL13, TGFβ, and CXCL5. Similarly, the expression of proteins linked with immune signaling (e.g., CD40, CD28, BTK) was significantly increased. Among the few differently expressed proteins in C3d, we found downregulation of the cell-cycle protein aurora kinase A and pro-inflammatory CCL24 (Fig. 4B and C). Pathway analysis confirmed a significant enrichment of proteins involved in the immune response (including CD40, BCR, ICOS/ICOSL, and PD1/PDL1 signaling), cytokine signaling (IL6, IL2, TGFβ), cell cycle (such as P53 and P38 signaling), senescence and metabolism (including GP6 and ID1 signaling) with opposite trend in C1 compared with C2 (Fig. 4D). We further investigated whether the altered proteins in the BM translated in similar changes in the PB. As expected, there was no significant correlation among the proteins of the two compartments (RV = 0.19, P = 0.26; Supplementary Fig. S5A–S5B), supporting the relevance of a tumor-specific milieu.

Using LC-TOFMS, we found a total of 25 of 44 serum metabolites with a significantly different concentration between at least 2 of the clusters (Fig. 5A). The most notable differences were detected in the levels of 3-Phosphoglyceric acid (3-PG) and 3-hydroxybutyric acid. The 3-PG concentration was undetectable in C3d, moderate in C2 (106.4 μmol/mL) while increased in C1 (188.5 μmol/mL). 3-PG is an intermediate metabolite of glycolysis whose conversion to 2-PG is dependent on phosphoglycerate mutase 1 (PGAM1), which is a mTOR-dependent enzyme with a critical role in promoting the Warburg effect (aerobic glycolysis; ref. 33). Hence, activation of mTOR pathway increases expression of PGAM1, which promotes glycolysis and in turn lowers 3-PG levels. While transcriptional activation of mTOR signaling in C3d and C2 explains the undetectable/low levels of 3-PG in these clusters, its transcriptional silencing in C1 may suppress glycolysis with a concordant increase of 3-PG due to failure to process this metabolite. 3-hydroxybutyric acid was undetectable in C1, low in C3d (60.6 μmol/mL), while elevated in C2 (316.2 μmol/mL). Normally, 3-hydroxybutyric acid increases during ketosis as an end product of fatty acid oxidation (41). This is in line with the preferential mitochondrial metabolism identified in C1 and C3d. In addition, we observed a higher level of lactic and malic acid in C2 compared with C1 and C3d, as result of preferential activation of glycolysis in C2 (Fig. 5B). We also found an increase of multiple amino acids, which however was similar among the three groups (Supplementary Fig. S6). We then explored whether the altered metabolites in the PB reflected those in the BM. However, this association did not seem significant (RV = 0.17, P = 0.30; Supplementary Fig. S7A–S7B). Collectively, these data confirmed the metabolic rewiring linked with the three molecular clusters, supporting the prevalence of mitochondrial metabolism in C1, glycolysis in C2 and their combination in C3d.

To confirm a different immune response mainly between C1 and C2, we performed CyTOF on fresh-frozen single-cell suspension (36) from matched BM biopsies of 11 patients in our cohort (C1, n = 7; C2, n = 4). We first evaluated whether C1 and C2 displayed a distinct T-cell profile. To this end, we performed t-distribution stochastic neighbor embedding (tSNE) analysis to characterize the immune profile in the tumor microenvironment. To decrease variability among the samples, we concatenated all files of CD3⁺ T cells. Overall, the number of CD4⁺ T cells was more abundant in C2 compared with C1 while the opposite was observed for CD8⁺ T cells (Fig. 6A). Among CD4⁺ T cells, activated T follicular helper cells (T_FH, PD1⁺TIGIT⁺ICOS⁺; P = 0.02) were predominantly expressed in C2. In support of T-cell activation, there was an increased expression of the costimulatory receptors CD69 (P = 0.02; Supplementary Fig. S8A) and BTLA (P = 0.01; Supplementary Fig. S8B) in C2 compared with C1. Also, T regulatory cells (T_reg, PD1⁺CD25⁺CTLA4⁺; P = 0.008) were significantly increased in C2 (P = 0.04; Supplementary Fig. S7C). In contrast, C1 was characterized by a higher expression of T effector memory cells (T_EM, CD45RA⁺CCR7⁻CD27⁻; P = 0.001). The concurrent presence of CD57 and KLRG1 suggested that these are senescent T_EM cells with reduced replicative capacity and decreased survival (Fig. 6A). To confirm our findings, we performed SPADE clustering analysis. The three main identified clusters were confirmed to be T_FH, T_reg, and T_EM cells (Fig. 6,B and C; Supplementary Fig. S9).

Patient characteristics based on their clinicopathologic diagnosis at the time of sample collection are summarized in Supplementary Table S1. Besides the more extensive BM involvement and elevated serum IgM monoclonal protein in patients with WM, the remaining clinical characteristics were similar in both diseases and matched the expected findings in these populations. The median time from diagnosis to treatment was 2.6 months (range, 0–55 months) in patients with WM while not reached in those with IgM MGUS (range of follow up, 26–281 months). Over the first 6 months, most of the patients with WM were either observed (n = 10, 41%) or underwent immunochemotherapy (n = 12, 50%), and only few received single agent rituximab (n = 2, 9%). In contrast, none of IgM MGUS required treatment during the follow up. We next compared the clinical characteristics of patients based on molecular clustering (Table 1). All patients assigned to C1 had a clinicopathologic diagnosis of WM, half of whom required initiation of treatment within 6 months (n = 7, 44%). Notably, all the patients requiring early treatment harbored MYD88 mutations (n = 7), with two of them co-expressing CXCR4 mutation. C2 included patients with both WM and IgM MGUS, but all had a lower proportion of BM involvement and serum IgM monoclonal gammopathy levels compared with C1. To better evaluate the C2 intra-cluster characteristics, we compared C2 IgM MGUS and C2 WM patients and noted similar clinical findings, including overall low median serum IgM levels (983 mg/dL in WM vs. 484 mg/dL in IgM MGUS; P = 0.54). The only difference within C2 subgroups was an expected higher BM involvement in patients with WM compared with IgM MGUS (20% in WM vs. 1% IgM MGUS; P = 0.01). Most of the patients in C2 were observed (64%). Even when treatment was required, rituximab monotherapy was the preferred initial regimen after a median of 12 months (range, 0.5–55 months) from diagnosis. Lastly, the clinical feature of patients in C3d was consistent with mild to moderate disease. However, while the BM involvement was minimal only in the patient with IgM MGUS, it was of at least 20% in those with WM. The majority of C3d patients received immunochemotherapy (60%) after a median of 7.8 months (range, 0.2–29 months), while the rest were observed (40%). Of note, the relapses were all patients with WM who were previously treated (range, 6–50 months), suggesting the possible persistence of cancer stem cells that may drive tumor growth and recurrence. The notion of such tumor niche is plausible given the transcriptional activation of stem cell and senescence genes in tandem with signaling pathways that control cell proliferation and antigen response. Loss of this balance may result in neoplastic transformation. Collectively, the data indicate that these clusters may improve identification of patients with molecular features reflecting different treatment need.

In this multi-omics study of IgM monoclonal gammopathies, agnostic to the clinicopathologic diagnosis, we identified 3 molecular clusters based on distinctive transcriptomic, proteomic, metabolomics, and immune features. The difference among them is linked to aberrant proliferation and/or inflammatory signaling which may drive their clinical behavior. While cell proliferation prevailed in C1, the inflammatory milieu defined C2, and intermediate features with combined proliferation and inflammation were observed in C3d, which may represent an intermediate state of this disease. While it is possible that C3d might be a posttreatment status, this probability is less likely as all the included patients had not been treated for at least 6 months before sampling.

Previous studies suggest that the pathobiology behind the progression from IgM MGUS to WM involves early genomic events followed by a stepwise accumulation of additional abnormalities (42). On the basis of this model, the proposed early event is the MYD88 L265P somatic mutation, which has been detected in 60% of IgM MGUS and 90% of patients with WM (43). However, in our study MYD88 L265P mutation was found in 80% of patients with WM. While this difference might be imputed to the lower limit of detection of variant allele frequency in WES compared with standard PCR (44), a similar mutation rate was also identified by Abeykoon and colleagues using AS-PCR (37), therefore supporting our findings. Another likely early event is mutations of CXCR4, which have been reported in 9% of IgM MGUS and 23% of WM cases (9). However, CXCR4 mutations are primarily subclonal, which might suggest a secondary acquisition after the MYD88 L265P mutation (7, 8, 45). Recently, a study using targeted next-generation sequencing of 11 selected genes (including MYD88 and CXCR4) confirmed an higher number of somatic mutations in patients with WM compared with those with IgM MGUS (8). Similarly, the frequency of copy-number alterations has also been shown to progressively increase from IgM MGUS (36%) to WM (73%–82%) with some of the most common abnormalities in WM (del 6q, +18q, +4, +12) rarely seen in IgM MGUS (46). Surprisingly, our WES data did not identify significant differences in the frequency of MYD88 and CXCR4 mutations among the three clusters with an overall similar genomic landscape. These contrasting results might stem from the reclassification of patients based on biologic rather than clinical criteria, or alternatively they might be due to the small number of patients included in our analysis. On the other hand, we found a significant enrichment in alterations of TNFAIP3 (including gene mutations and 6q deletion) in C1 compared with the other groups, suggesting a plausible role in malignant transformation. TNFAIP3 gene encodes the tumor suppressor A20, which terminates NF-kB activation (47). Loss of function of A20 might form an oncogenic cooperation with mutated MYD88 that converges on NF-kB, ultimately promoting disease progression. It should be noted that 6q deletion may affect also other genes involved in cell growth (e.g., BCLAF1, IBTK, HIVEP2; ref. 48) and antigen presentation (e.g., HLA-DRA, HLADQAQ, HLA-DOA; ref. 49). In addition, our data supports the role of nongenetic events in the evolution to WM. For example, spliceosome signaling was significantly enriched in C1 and moderately activated in C3d as opposed to C2. The oncogenic properties of spliceosome have been shown in several hematologic and solid malignancies as they lead to proliferative dysregulation in tandem with immune escape (50–52). It is therefore possible that posttranslational events might favor transformation from the indolent C2 to the more aggressive C1. However, mechanistic studies are warranted to elucidate if and how aberrations in this pathway may contribute to disease progression and aggressive behavior in WM.

On the basis of our findings, we propose that C2 represents an early stage of disease. In support of this, patients in C2 presented with low BM involvement and serum IgM levels, which are both suggestive of low tumor burden. The characteristic increase in the inflammatory milieu and prolonged cell survival are also reminiscent of the early antigenic stimulation that seems to represent the foundation of disease (53). Over time, this stemness property combined with inflammatory stimulation might trigger additional causative events that drive transformation to C1. During this process, malignant cells may transition through the intermediate stage C3d, which shares characteristics of C2 and C1, including the coexistence of inflammation and cell proliferation. This inflammatory phenotype may explain the co-clustering of controls with an underlying reactive state. Alternatively, C3d could also represent an intermediate stage after treatment response in which normal BM coexists with cancer stem cell niches. Loss of this balance enables aberrant silencing of genes, which facilitate lymphoid transformation. It is possible that this process may be accelerated by dysregulation of cell-cycle checkpoints, which promotes tumor growth and in turn provides survival advantage to one or more cellular subclones. Other preneoplastic events such as clonal hematopoiesis of indetermined potential (CHIP) are unlike as these commonly present in elderly people (>70 years) and involve myeloid progenitors harboring mutations of chromatin modifier genes (e.g., TET2, DNMT3A, ASXL1; ref. 54), as opposed to our study, in a relative younger population, which identified a lymphoid phenotype with none of the typical CHIP mutations. Finally, C1 depicts a more advanced disease stage, in which the inflammatory response subsidies in favor of uncontrolled cell proliferation, which may be further enhanced by failing immunosurveillance. Accordingly, patients in C3d had a higher BM involvement, but similar serum IgM levels compared with C2, while those in C1 presented a high tumor burden as demonstrated by a more severe BM involvement and serum IgM levels compared with the other two groups. In line with the evolving spectrum of disease, we identified concordant changes in the cellular metabolism, which initially favor glycolysis in C2, transitioning through the intermediate C3d state that combines glycolysis and mitochondrial metabolism, to then preferentially depend on mitochondria in C1. This metabolic rewiring can be explained by an emerging notion that posits mitochondrial metabolism as instrumental in promoting malignant transformation via both cell-intrinsic and cell-extrinsic mechanisms (55).

To the best of our knowledge, this is the first report that identifies distinct biological clusters of IgM monoclonal gammopathies. While these diseases have been previously evaluated with genomic and transcriptomic analysis, to date no study has reported a comprehensive multi-omics investigation. Our finding of a substantial subset of patients with WM with a multi-omic phenotype similar to IgM MGUS uncovers the limitations of the current diagnostic criteria and suggests the need for improvements in selecting patients who require treatment. Furthermore, the presence of a cohort with intermediate features may open the opportunity to preventive intervention to limit transformation to more aggressive disease.

From a therapeutic perspective, our data reinforce the importance of precision medicine and pave the way for agents that restore immunosurveillance. While the current therapeutic paradigm depicts single agent BTK inhibitor as mainstay of targeted therapy, most responses remain partial (56, 57), raising the question as to whether a combined therapy that might more effectively inhibit the NF-kB pathway may enhance clinical responses (58, 59). Learning from other B-cell lymphomas, additional effort may be needed to focus on identifying rational combination targeted therapy based on genetic (60, 61) and epigenetic (62) profiles. Hence, dual inhibition of NF-kB pathway might be required in our therapeutic armory to maximize treatment efficacy and prolong survival. Furthermore, additional consideration should be paid to novel immunotherapy strategies given the striking immunosuppression in C1, which might account for resistant disease. While no conclusions can be drawn regarding immune checkpoint blockade in WM due to severe toxicity and small patient numbers in the initial investigation (63), CAR T-cell therapy has recently shown promise in WM, warranting further studies (64).

In conclusion, we have performed the first comprehensive analysis of IgM monoclonal gammopathies, which identifies 3 molecular clusters with distinct multi-omic characteristics. Even though our findings will need independent validation, it is important to emphasize the fact that the multilayered analysis of each subject with different omics corroborating the same findings in this study provides a strong internal validation. These clusters have the potential to change the current definitions of IgM MGUS and WM and lay the foundation for a biological approach to this spectrum of diseases.

J. Paludo reports grants from Karyopharm and Biofourmis outside the submitted work. E. Braggio reports personal fees from DASA outside the submitted work. S. Ailawadhi reports grants and other support from GSK, BMS, Janssen, Amgen, and Pharmacyclics; other support from Sanofi, Takeda, Beigene, AstraZeneca, and Regeneron; and grants from Cellectar, AbbVie, Xencor and Ascentage outside the submitted work. M.A. Gertz reports personal fees from Ionis/Akcea, Prothena, Sanofi, Janssen, Aptitude Health, Physicians Education Resource, AbbVie, Johnson & Johnson, Celgene, Research to Practice, Ionis/Akcea, Juno, Sorrento, and i3Health; honorarium from Alnylym; and grants and personal fees from Ashfield outside the submitted work. A.J. Novak reports other support from Bristol-Myers Squibb outside the submitted work. No disclosures were reported by the other authors.

P. Mondello: Conceptualization, data curation, formal analysis, methodology, writing–original draft. J. Paludo: Formal analysis, investigation. J.P. Novak: Software, formal analysis. K. Wenzl: Data curation, formal analysis. Z.-Z. Yang: Writing–review and editing. S. Jalali: Investigation. J.E. Krull: Software, formal analysis. E. Braggio: Writing–review and editing. S. Dasari: Software, formal analysis, writing–review and editing. M.K. Manske: Investigation. J.A. Abeykoon: Investigation. V. Sarangi: Investigation. P. Kapoor: Investigation, writing–review and editing. A. Paulus: Investigation, writing–review and editing. C.B. Reeder: Investigation. S. Ailawadhi: Investigation. A.A. Chanan-Khan: Investigation. R.A. Kyle: Investigation, writing–review and editing. M.A. Gertz: Investigation, writing–review and editing. A.J. Novak: Software, formal analysis, investigation, writing–review and editing. S.M. Ansell: Conceptualization, resources, data curation, supervision, funding acquisition, methodology, project administration, writing–review and editing.

This work was supported by grants from the Binks Foundation, the International Waldenström Macroglobulinemia Foundation, and the Predolin Foundation.

The publication costs of this article were defrayed in part by the payment of publication fees. Therefore, and solely to indicate this fact, this article is hereby marked “advertisement” in accordance with 18 USC section 1734.