Waldenstrom macroglobulinemia (WM) and its precursor IgM gammopathy are distinct disorders characterized by clonal mature IgM-expressing B-cell outgrowth in the bone marrow. Here, we show by high-dimensional single-cell immunogenomic profiling of patient samples that these disorders originate in the setting of global B-cell compartment alterations, characterized by expansion of genomically aberrant extrafollicular B cells of the nonmalignant clonotype. Alterations in the immune microenvironment preceding malignant clonal expansion include myeloid inflammation and naïve B- and T-cell depletion. Host response to these early lesions involves clone-specific T-cell immunity that may include MYD88 mutation–specific responses. Hematopoietic progenitors carry the oncogenic MYD88 mutations characteristic of the malignant WM clone. These data support a model for WM pathogenesis wherein oncogenic alterations and signaling in progenitors, myeloid inflammation, and global alterations in extrafollicular B cells create the milieu promoting extranodal pattern of growth in differentiated malignant cells.
These data provide evidence that growth of the malignant clone in WM is preceded by expansion of extrafollicular B cells, myeloid inflammation, and immune dysfunction in the preneoplastic phase. These changes may be related in part to MYD88 oncogenic signaling in pre–B progenitor cells and suggest a novel model for WM pathogenesis.
This article is highlighted in the In This Issue feature, p. 549
Waldenstrom macroglobulinemia (WM) is a clinically distinct B-cell malignancy characterized by progressive growth and accumulation of an IgM-expressing lymphoplasmacytic clone (1, 2). Clinical malignancy in WM is preceded by a precursor state termed IgM monoclonal gammopathy of undetermined significance (MGUS; ref. 1). Genetic analyses have demonstrated recurrent mutations in the myeloid differentiation primary response-88 (MYD88) gene in the vast majority of WM patients, wherein they are found to involve the entire expanded B-cell clone (3–5). However, mutations in the MYD88 gene are neither specific nor sufficient for the pathogenesis of WM, and can be detected in IgM MGUS as well as in other B-cell lymphomas, including central nervous system or testicular lymphomas (4, 6–11). The clinical phenotype of both IgM MGUS and WM is distinct and characterized by the infiltration of the bone marrow as the dominant site of tumor growth. Although patients with IgM MGUS are asymptomatic and exhibit low clonal burden, onset of clinical malignancy is often characterized by the development of anemia and progressive tumor infiltration. In spite of recent advances in the genetics of WM tumor cells, the mechanisms underlying the pattern of tumor growth as well as the cell of origin of this malignancy remain poorly understood (12–14).
Current therapy for WM is based on monoclonal antibodies targeting B cells, strategies targeting B-cell receptor (BCR) signaling, as well as systemic chemotherapy (1). Although these therapies lead to high rates of tumor regression, they are not curative. The immune system has emerged as a promising strategy to target tumors, including those in the lymphoid system. Properties of the host immune response and the nature of the tumor immune microenvironment are important determinants of outcome in several malignancies and may affect the evolution of premalignant states. However, current information about the immune microenvironment in WM and IgM MGUS is limited to studies evaluating only a small number of parameters (12, 15–17), and data with application of newer high-content single-cell approaches in WM and MGUS to probe early changes in the immune microenvironment are limited. Importantly, evidence that the immune system can mediate tumor-specific immune recognition of WM tumor cells is lacking. In this study, we combine several single-cell proteomic and genomic tools with functional studies to gain insights into immunobiology and host response in WM and its precursor, IgM MGUS.
To better understand early events in the pathogenesis of WM, we analyzed the bone marrow mononuclear cells of patients with IgM MGUS (n = 8) or newly diagnosed/previously untreated WM (n = 8; Supplementary Table S1 for patient characteristics) with a combination of mass cytometry and cellular indexing of transcriptomes and epitopes by sequencing (CITE-seq) and compared it with those from age-matched healthy donors (HD; n = 5; Fig. 1A). Mass cytometry analyses revealed an increase in CXCR5neg B cells in WM, along with a relative decline in myeloid cells (Fig. 1B; Supplementary Fig. S1). Importantly, these changes were observed as early as the MGUS stage, when the clonal burden is low. Mass cytometry analysis also revealed an increase in bone marrow T cells in patients with MGUS (Fig. 1B; Supplementary Fig. S1). CITE-seq analyses on 84,128 single cells (HD 20,946 cells, MGUS 29,780 cells, and WM 33,402 cells) identified 42 distinct clusters, which were also broadly classified into B/lymphoplasmacytoid (LPC), T/natural killer (NK), monocyte/myeloid and precursor cell types (Fig. 1C and D; Supplementary Figs. S2 and S3), and again confirmed alterations in myeloid, B-, and T-cell lineages (Fig. 1E; Supplementary Fig. S4). Together, these data utilize complementary tools to show that bone marrow cells in both WM and IgM MGUS are characterized by alterations in several hematopoietic lineages compared with HD counterparts.
Changes in B Cells and LPCs
As WM is a B-cell malignancy, we first focused on analyses of CD19+ cells by mass cytometry. Changes in the B-cell compartment were evident as early as IgM MGUS, with an increase in CXCR5neg B cells and a decline in naïve B cells, although the clonal population at this stage is still small (Fig. 2A–D; Supplementary Fig. S5). CXCR5 is well recognized as a critical gene for entry of B cells into lymphoid follicles (18), and CXCR5negCD21lo B cells have been recently appreciated as cells in the extrafollicular (EF) pathway of B-cell differentiation (19). The CXCR5negCD19+ population consisted of two distinct populations, with one representing mature B cells (CXCR5neg population 1) and another expressing CD38 that includes LPCs (CXCR5neg population 2; Fig. 2B). The malignant clone in WM is characterized by the expression of clonotypic kappa or lambda light chain, which secretes a clonal Ig clinically documented as an M spike. Although the observed alterations in B cells were enriched in the malignant B cells expressing clonotypic light chain, they were also observed in the B cells expressing the opposite (or nonclonal) light chain (Fig. 2E). To validate these findings, we also analyzed CD19+ B cells by CITE-seq (overall strategy in Supplementary Fig. S2). Uniform Manifold Approximation and Projection (UMAP) clustering of B cells based on antibody staining identified 10 distinct clusters, which again revealed progressive loss of IgD+CD27− naïve B cells (cluster 2) from HDs to IgM MGUS and WM, along with an increase in clusters containing CXCR5neg B cells (cluster 3; Fig. 3A–C). Several of these clusters also expressed plasma cell markers such as BCMA and CD138, consistent with lymphoplasmacytoid differentiation of B cells in WM (Supplementary Fig. S6). Trajectory analysis of CD19+ B lineage cells identified two distinct trajectories originating from pre–B cells, with one enriched in HD and another in WM (Supplementary Fig. S7). UMAP clustering of CD19+ cells based on transcriptome identified 24 distinct clusters (Fig. 3D–F; Supplementary Fig. S8). Of note, B cells from WM patients as well as some MGUS patients formed distinct clusters, indicating that the malignant clone in each patient is transcriptionally distinct, analogous to prior studies in other malignancies (refs. 20, 21; Fig. 3E). Clonality of B cells in these patient-specific clusters was confirmed by BCR sequencing as well as computational imputation of BCRs from single-cell transcriptome data (Supplementary Fig. S9). Examples of top genes characterizing specific clones in individual patients include genes implicated in lymphomagenesis, including EZH2, CD79b, CXCR4, BCL2, and LT-β (Fig. 3F). Overall, B cells from WM and IgM MGUS patients were transcriptionally distinct from those in HDs. Top pathways increased in WM and IgM MGUS B cells included the NF-kB pathway and the IL6/STAT3 pathway, consistent with prior genomic studies with bulk tumor cells (Supplementary Figs. S10 and S11; ref. 7).
Changes in Nonclonal Light Chain Expressing B Cells and Pre–B Cells
Based on the findings from mass cytometry relating to B cells expressing nonclonal light chain (for example, kappa light chain in a patient with IgM lambda M spike), we further evaluated these cells by CITE-seq. The polyclonal nature of these B cells and lack of clonal BCR was verified by BCR sequencing (shown as an example in Supplementary Fig. S9C). Transcriptome of these nonclonal light chain–expressing B cells in WM as well as MGUS patients also revealed several changes compared with their HD counterparts, with increased expression of genes such as Pim1 and TCF4, and decrease in HLA-DR previously implicated in lymphoma biology (Supplementary Fig. S12).
In addition to mature B cells, transcriptomes of CD19+CD10+ pre–B cells from WM or IgM MGUS patients also revealed several alterations compared with HD pre–B cells (Supplementary Fig. S13). Interestingly, signatures derived from differentially expressed genes in WM B cells compared with HD were also found to be enriched in WM pre–B cells (Supplementary Fig. S14). Top pathways enriched in WM pre–B cells included the IL6–STAT3 pathway, as seen in WM B cells as well (Supplementary Fig. S15). Taken together, these data demonstrate that clonal populations in WM and IgM MGUS originate in the context of more global alterations in B as well as pre–B cells, with resultant expansion of EF B cells.
Changes in the Myeloid Compartment
Besides B cells, the most prominent changes in the immune microenvironment were in the myeloid compartment (Fig. 1B and E). UMAP clustering of CD3- and CD19-depleted cells was utilized to identify distinct clusters that broadly fell into myeloid, NK, and progenitor clusters and were analyzed separately (Fig. 4A; Supplementary Fig. S16). Analysis of the CD14+/CD11c+ myeloid compartment revealed that both IgM MGUS and WM were associated with a marked decline in clusters of classic monocytes (cluster 1) and a relative increase in a distinct cluster of myeloid cells (cluster 7), particularly in MGUS (Fig. 4A–D; Supplementary Fig. S17). This myeloid population was characterized by a distinct genomic signature of inflammation-associated genes such as IL1β, CCL4, CCL3, IL6, NLRP3, and CXCL3 (Fig. 4E). Pathway analysis was consistent with inflammation-associated signaling, particularly in MGUS (Fig. 4F). Together, these data suggest that activation of myeloid inflammation is an early feature of MGUS, occurring before the evolution of the malignant clone.
Changes in Innate Cells
Evaluation of CD3−CD56+ NK clusters identified alteration in NK cells, most prominently in the WM cohort (Fig. 4A and C; Supplementary Fig. S18). These changes were characterized by greater expression of lytic genes/exhaustion markers in WM NK cells, along with loss of interferon signature in these cells (Supplementary Fig. S18), consistent with the concept of immune exhaustion and dysfunction in NK cells with evolution of WM. Pathway analysis of differentially expressed genes in NK cells from both MGUS and WM revealed decreased interferon signaling consistent with NK dysfunction (Supplementary Fig. S19).
Changes in T Cells and Evidence of Tumor-Specific Immunity
Analysis of the T-cell compartment by mass cytometry revealed evidence of T-cell activation with a decrease in naïve CD8 T cells and increased CD8+ T cells expressing granzyme (Fig. 5A–D). These cells also expressed high levels of Eomes and T-bet and low levels of TCF1 (Fig. 5B). Notably, expression of PD1 was not a prominent feature of these cells, which expressed Tigit instead (Fig. 5B). Analysis of transcriptomes of T cells between WM or MGUS cohorts and HDs revealed that, as with mass cytometry, the most prominent changes were increases in lytic genes and markers associated with T-cell exhaustion in WM/MGUS CD8+ T cells, concurrent with loss of interferon-response genes (Fig. 5E and F). Changes in interferon-response genes were also seen in CD4+ T cells (Supplementary Fig. S20). Pathway analysis confirmed changes in T cells, most notably as decline in interferon response (Supplementary Figs. S21 and S22). Together, these data show that changes in the T-cell compartment begin early in MGUS, before the establishment of progressive malignant clone, and are characterized by progressive depletion of naïve T cells and enrichment of terminal effector T cells.
Although the T-cell phenotypes are consistent with evidence of T-cell activation in vivo, whether the immune system mediates specific recognition of WM cells in vivo is not known. To test this directly, we analyzed reactivity of freshly isolated blood/bone marrow T cells to autologous monocyte-derived dendritic cells (DC) loaded with autologous tumor cells. In all patients tested, freshly isolated T cells could mediate tumor-specific immune recognition (Fig. 6A). To further evaluate if these T cells still retained proliferative capacity and could be enriched further, we cultured them with autologous tumor-loaded DCs. DC-mediated stimulation led to expansion of tumor-specific T cells in culture (Fig. 6B). Tumor reactivity persisted following depletion of CD4+ T cells, indicating that it consisted predominantly of CD8+ T cells (Fig. 6C). MYD88 L265P is a common and recurrent mutation found in WM cells. To test whether antitumor T-cell reactivity includes MYD88-specific T cells, we analyzed reactivity to mutant and wild-type peptides in blood/marrow mononuclear cells, with peptide-specific stimulation. Data from representative responder and nonresponder patients are shown in Fig. 6D. Overall, reactivity specific for mutant MYD88 peptide was detected in 5 of 15 patients tested (Fig. 6E). Together, these data demonstrate that the immune system is capable of mediating tumor-specific recognition of WM cells, which includes reactivity to mutant MyD88 in a subset of patients.
Genetic Alterations in Precursor Cells
The findings discussed above demonstrate that although the expanded clone in WM has the phenotype of mature IgM+ B cells, the evolution of the malignant clone occurs in the context of more global genomic and phenotypic aberrancies in hematopoiesis, which precede the expansion of the clone. In order to test whether genomic alterations (including MYD88 mutations) previously characterized in the mature B-cell clone can also be detected in earlier progenitors, we flow-sorted CD19+CD10− mature B cells, as well as less mature progenitors (LinnegCD19+CD10+ pre–B cells, LinnegCD19negCD34neg progenitors, and LinnegCD19negCD34+ progenitors) from five patients with WM and two HDs and analyzed their genomes by exome sequencing (gating strategy in Supplementary Fig. S23). In all patients studied, the MYD88 L265P mutation could be detected in pre–B progenitor compartments and involved the entire mature B-cell clone (Fig. 7A–C). In four of the five patients studied, the MYD88 L265P mutation was detected with variant allele fraction (VAF) of >10% in at least one progenitor compartment (>10% VAF in CD19+CD10+ cells in three of five patients), indicating its detection in over 20% of progenitor cells (Fig. 7C). In contrast, MYD88 L265P mutations were not detected in any fraction from HD marrow (Supplementary Fig. S24). The finding that MyD88-mutant fractions show clear expansion (>10% VAF/20% of cells) in CD19+CD10+ pre-B fractions suggests that this is not a function of low-level contamination during flow sorting. This is also supported by the finding that several SNVs were not shared between sorted fractions. The detection of MYD88 L265P by exome sequencing was verified by allele-specific qPCR, which yielded similar VAF estimates (Supplementary Fig. S25A) and confirmed the detection of MYD88 L265P mutations in precursor compartments (Supplementary Fig. S25B). MYD88 L265P was also detected by PCR at low levels in the CD19+CD10− cells and at least one other progenitor compartment in MGUS marrow but not in any of the sorted populations in HD (Supplementary Fig. S25C and S25D). We also tested whether MYD88 mutations could be detected in circulating B cells expressing the nonclonal light chain. Analysis of flow-sorted B cells (based on expression of CXCR5 and Ig light chain) by PCR revealed that low levels of MYD88 mutation could be detected in the nonclonal IgL+ B cells in WM patients but not HDs (Supplementary Fig. S25E and S25F). Together, these data demonstrate that MYD88 L265P mutations can occur early in lymphopoiesis in WM, before the expansion of the malignant B-cell clone, and are not restricted to the expanded clonal population.
These data demonstrate that the preneoplastic phase of WM is characterized by a distinct immune landscape with marked alterations in the nonclonal B-cell compartment. Coupled with the detection of MYD88 L265P mutations in early pre–B progenitors, these data also suggest that although the expanded clonal population in WM has the phenotype of a mature B cell, the clone emerges in the context of alterations in earlier B progenitors that begin well before expansion of the malignant clone.
Changes in the B-cell compartment are manifest not only as a decline in normal naïve B cells but a marked increase in EF B cells with the CXCR5negCXCR4+ phenotype. Importantly, although clonal B cells in WM are also often CXCR5neg, changes in EF B cells are also observed in those expressing nonclonal immunoglobulin light chain. These nonclonal cells also overexpress several genes, such as Pim1 and TCF4, previously implicated in lymphomagenesis (22, 23). These data are, therefore, consistent with the concept that WM and IgM MGUS clones develop in the context of global alterations of B cells with expansion of aberrant nonclonal EF B cells. Recent studies have demonstrated enhanced autoreactivity in EF B cells, implicating them in several autoimmune states (19, 24). The observed expansion of EF B cells in WM and IgM MGUS may, therefore, underlie the increased risk of autoimmunity in these disorders (25). The observed decline in naïve B cells may lead to reduced response to pathogens and vaccines in these patients (26, 27).
The mechanism underlying the observed increase in EF B cells is not known but may be related in part to aberrant MYD88 signaling, based on our finding of MYD88 mutations in pre–B cells. MYD88 has been implicated in both B-cell and myeloid differentiation, with the latter being responsible for its name (28, 29). MYD88 plays a central role in Toll-like receptor (TLR) signaling, and this pathway was recently suggested to play a critical role in B-cell differentiation along the EF pathway (30, 31). Transformation of cells in the EF pathway may also help explain the predominantly extranodal/bone marrow–restricted pattern of tumor growth in WM. It is also notable that in addition to WM, MYD88 L265P mutations are also highly prevalent in some other extranodal lymphomas such as the central nervous system or testicular lymphoma (32).
Our studies also provide the first evidence of tumor-specific immune recognition in WM. T-cell compartment in WM and IgM MGUS also exhibited depletion of naïve T cells and concomitant increase in effector populations consistent with in vivo activation and similar to findings in other hematologic malignancies (21, 33). Importantly, evidence of immune dysfunction with attrition of TCF1+ stem-like memory T cells and increase in terminally differentiated T as well as NK cells begins early in the MGUS stage, before the malignant clone is established, similar to recent findings in IgG MGUS (33). Although these T cells express markers of exhausted T cells, tumor-specific T cells could be expanded in culture and retained capacity for tumor-reactive cytokine production. To measure tumor-specific immunity, we prioritized DC-based cross-presentation assays, which allow detection of immunity against multiple antigens. Our studies also show that at least in some patients, the MYD88 L265P mutation is immunogenic. Although WM-specific T cells were present in all patients tested, MYD88-specific reactivity was observed only in a subset, suggesting that other neoantigens may also be targets of WM-specific response in this setting. Immunogenicity of MYD88 L265P in some WM patients is consistent with recent studies priming mutation-specific T-cell receptors (TCR) in HDs (34); however, further studies are needed to better understand the functional properties of MYD88 L265P–specific T cells in WM tumor bed and identify specific MHC haplotypes wherein high-affinity TCRs could be harnessed for immunotherapy. Further studies are also needed to quantify the degree to which mutant MYD88 is a target of global WM-specific immunity. This analysis was hampered in the current study by limitations of sample and HLA restriction. Harnessing preexisting immune response in WM patents may, in principle, also be feasible via blockade of immune checkpoints, particularly Tigit, overexpressed in these cells.
Our data suggest a novel model for the pathogenesis of WM, wherein aberrant increase in polyclonal CXCR5negCXCR4+ EF B-cell compartment creates the milieu for the emergence of MGUS (Fig. 7D). MYD88 mutations are present not only in the malignant B-cell clone, but also in earlier pre–B progenitors. In some patients, the mutant pre–B compartment was expanded to >20% of the pre–B compartment. The presence of mutations in even earlier progenitors requires further study; however, in at least one patient (WM patient 1 in Fig. 7), the mutant fraction accounted for >20% of CD34+ CD19− cells, indicating interpatient heterogeneity relating to detection and expansion of the mutant fraction in earlier progenitors in some WM patients. Humanized models (35) that permit efficient growth of WM cells may be needed to better understand the biology of progenitor compartments in WM. Evolution of precursor states to WM may depend on both acquisition of additional genomic changes in tumor cells and alterations in nonmalignant cells. This model is similar to emerging data for polyclonal origins of non-IgM monoclonal gammopathies with the pre-MGUS phase beginning in early decades of life (36–38). It is also consistent with the view that MYD88 mutations are not sufficient for the development of WM or other MYD88-driven lymphomas, and additional alterations in tumor cells and/or host response are needed (8, 9). The immune microenvironment in the precursor state is also characterized by myeloid inflammation. The inflammatory environment may promote the malignant clone but also induce immune exhaustion. Our data that the immune system can mediate WM-specific immune responses support the potential for immune-mediated control, also suggested in MYD88-driven murine lymphoma models (39). The concept that alterations in hematopoietic progenitors may underlie the development of mature B-cell malignancies has also been previously suggested in the case of some other B-cell malignancies such as chronic lymphocytic leukemia and hairy cell leukemia (40–42).
These data have several implications for biology and therapy of WM. Homing properties of EF cells based on the lack of CXCR5 expression may explain the predominantly extranodal pattern of tumor involvement in WM. The presence of oncogenic mutations in earlier progenitors lacking clonal BCR suggests that oncogenic signaling may contribute to changes not just in tumor cells, but also in nonmalignant cells in the tumor microenvironment. Targeting MYD88-mutant progenitors as well as harnessing the ability of the immune system to target these tumors may be needed to improve therapy and curability of WM. With increasing appreciation of age-associated genomic alterations in hematopoiesis, this issue also needs broader evaluation across human tumors.
Bone marrow specimens from patients with a diagnosis of MGUS or WM (1) were obtained following written informed consent from patients. Studies were approved by the Emory institutional review board and were conducted in accordance with recognized ethical guidelines. Specimens from HDs were purchased from All Cells, Inc. Patient characteristics are shown in Supplementary Table S1.
Thawed bone marrow mononuclear cells (BMMNC) were stained with custom panels of metal-conjugated antibodies according to manufacturer-suggested concentrations (Fluidigm; antibodies as noted in Supplementary Table S2). Cells were fixed, permeabilized, and washed in accordance with the manufacturer's cell-surface and nuclear staining protocol as described. After antibody staining, cells were incubated with intercalation solution, mixed with EQ Four Element Calibration Beads (catalog number 201708), and acquired using a Helios mass cytometer (all from Fluidigm). Gating and data analysis were performed using Cytobank (https://www.cytobank.org). Viable cells and doublets were excluded using cisplatin intercalator and DNA content with iridium intercalator. Equal numbers of cells from each donor were utilized when data were concatenated prior to analysis.
After thawing bone BMMNCs, live cells were isolated using the EasySep Dead Cell Removal (Annexin V) Kit (STEMCELL Technologies) per the manufacturer's protocol. Cells were then stained using the TotalSeq-C antibody cocktail (Supplementary Table S3) following the 10x Genomics protocol for Chromium Single-Cell Immune Profiling with Feature Barcoding Technology (ver. 1.0). Gene expression and cell-surface libraries were prepared according to the protocol from 10x Genomics. The quality of the prepared libraries was assessed using the Agilent HS Bioanalyzer 2100, and sequencing was conducted with the Illumina NovaSeq 6000. To minimize variability, all patients and controls were sequenced together.
Data Analysis for CITE-seq
Reads were aligned to human reference sequence GRCh38/hg38, filtered, deduplicated, and converted into a feature barcode matrix using Cell Ranger 4.0. Samples with low fraction reads (<70%), high reads mapped to antisense strand (>10%), and problematic barcode rank plot were excluded. Final samples with feature barcode matrix containing data for genes and 55 surface proteins were analyzed using Seurat V3 (43). Cells with fewer than 200 unique sequenced genes, more than 10% mitochondrial genes, and more than 70,000 sequenced features were filtered to exclude dead cells and doublets. The data were batch corrected using canonical correlation analysis method using Seurat V3 data integration pipeline. Gene expression for each cell was normalized using R package SCTransform (44), which fits the regularized negative binomial (NB) regression model. For each gene, unique molecular identifier (UMI) counts are considered as the response and cellular sequencing depth as the predictor variable to obtain regularized parameter estimates. The estimated regularized regression parameter was used to transform the UMI counts for each gene into Pearson residuals serving as scaled gene-expression data. A regularized NB model in comparison with log-normalization technique has been suggested to prevent overestimation of true variance in low-abundance genes (44). The final data set used for analyses consisted of 84,128 cells from 18 samples (healthy, n = 4; IgM MGUS, n = 7; and WM, n = 7). Protein expression data were treated as compositional data and normalized independently of transcriptomic data using centered log-ratio normalization, where counts were divided by the geometric mean of the corresponding feature across cells, and log-transformed (45). Principal component analysis (PCA) was then performed using the top 3,000 variable genes ranked by residual variance. Elbow method was utilized to identify the number of principal components that explained maximum variability in data, which resulted in 20 as an optimal number of principal components. Local neighborhood for each cell was defined by taking the 20 nearest neighbors in the k-nearest neighbor graph calculated using Euclidean distance in the PCA space. Rare and nonconvex cell populations were further identified using shared nearest-neighbor networks by calculating neighborhood overlap using the Jaccard index. Finally, Louvain community detection algorithm with resolution of 0.5 was used to identify 42 clusters of cells with similar gene-expression profiles. These 42 clusters were further classified into four major categories—B, T/NK, myeloid, and precursor cell types—based on the protein expression of CD19, CD3, and CD14 and confirmed by reference transcriptomic data sets of immune cells with SingleR (46). Cells in each of these four major categories were analyzed separately, and subclustering within each category was performed based on transcriptomic and protein expression data. The subcluster identity within each category was determined by using significantly differentially expressed genes and proteins using Wilcoxon rank-sum test between each cluster and the rest. Differentially expressed genes were denoted as statistically significant for the Bonferroni adjusted P < 0.05 with a fold change exceeding 1.2× (i.e., ≥ 1.2 or ≤ 0.83). Pathway analysis for significantly differentially expressed genes was performed with the DAVID Gene Functional Classification Tool (47) using the overrepresentation analysis method. Gene set enrichment analysis (GSEA; ref. 48) was used for pathway analysis based on preranked genes. Reactome (49) and hallmark (50) gene sets were utilized, and pathways with q-value < 0.05 were retained. Data are uploaded to Gene Expression Omnibus with accession number GSE179221.
BCR Sequencing and BCR Imputation from Transcriptomic Data
BCR sequencing reads from FASTQ were aligned to GRCh38-alts-ensembl using Cellranger VDJ pipeline version 6.0.2 to generate single-cell V(D)J (variable, diversity and joining) segments for heavy and light chains. The unique combination of heavy- and light-chain genes is referred to as clonotype. The clonotype data for BCR sequencing were mapped to Seurat object containing CITE-seq data by matching cell barcodes. Clonotypic frequency obtained from Cellranger output was depicted using pie chart. We also imputed the BCR clonotype from transcriptomic data using Cellranger VDJ pipeline 6.0.2. Clonotypic frequency was generated similar to BCR sequencing data.
Frozen BMMNCs or peripheral blood mononuclear cells from HDs (n = 3) and patients with IgM MGUS (n = 2) and WM (n = 5) were thawed and stained prior to cell sorting using BD FACSAria. BMMNCs were stained with live dead dye (Thermo Scientific) as well as antibodies to detect CD19 (clone HIB19, BD Bioscience) and CD3 (HIT3a), CD14 (M5E2), CD56 (HCD56), CD10 (HI10a), and CD34 (581)—all from BioLegend. Cells were gated to obtain single, living, Lin−(CD3−, CD14−, CD56−) cells, and four different fractions were collected for exome sequencing: CD19+CD10−, CD19+CD10+, CD19−CD34−, and CD19−CD34+. DNA was extracted using the AllPrep DNA/RNA Micro Kit (Qiagen) and used for exome sequencing or qPCR analysis to detect MYD88 L265P mutation. Peripheral blood mononuclear cells were labeled with live dead dye, and antibodies to detect CD19 and CXCR5 (clone RF8B2, BD) as well as Igkappa light chain (TB28-2) and Iglambda light chain (MHL-38) from BioLegend. Cells were gated to obtain single, live cells, and four different fractions of CD19+ cells were sorted: Kappa+CXCR5+, Kappa+CXCR5−, Lambda+CXCR5+, and Lambda+CXCR5− cells. DNA was extracted using the AllPrep DNA/RNA Micro Kit (Qiagen) and used to perform qPCR for the detection of MyD88 L265P.
Exome Sequencing and Data Analysis
Exome capture was performed using the IDT xGen v1 capture kit (Integrated DNA Technologies) according to the manufacturer's instructions. Exome FASTQ sequencing files were quality and adapter trimmed using Trim Galore (v0.6.4) with Cutadapt (v2.8; https://www.bioinformatics.babraham.ac.uk/projects/trim_galore/) prior to mapping to the GRCh38 reference genome (GRCh38.d1.vd1.fa.tar.gz; https://gdc.cancer.gov/about-data/gdc-data-processing/gdc-reference-files) with Burrows–Wheeler alignment tool (v0.7.17-r1188; ref. 51). Mapped SAM files were converted to BAM files with SAMtools (v1.10; ref. 52), and read groups and duplicate reads were marked using the Genome Analysis Toolkit (v184.108.40.206; ref. 53). Variant alleles were called with Mutect2 (v220.127.116.11; ref. 53), with all BAM files for a given patient included and using the 1000 Genomes Project “Panel of Normals.” VCF files were converted to a matrix format in R (v3.6.3) with the bedr (v1.0.7) package, and annotation of variants was performed with the VariantAnnotation (v1.32.0) package (54). Variants with at least 40× coverage in all samples were analyzed for different allele frequencies between any two cell types using a two-sided Fisher exact test with a Benjamini–Hochberg FDR correction. VAFs from Mutect2 are shown in the figures, and Venn diagrams of mutations required a minimum of two reads to call a mutation present in a given cell type. Overall, exome sequencing was performed to a median depth of 162,486,965 reads with greater than 99% mappability, resulting in a median exome coverage of 213×.
MYD88 L265P Allele-Specific qPCR
Allele-specific qPCR for MYD88 L265P was performed as previously described (55). Briefly, the qBiomarker Somatic Mutation PCR Assay (Qiagen) MYD88_85940 was used to determine the total MYD88 fraction as compared with the MYD88 L265P alleles on a CFX96 Touch Real-Time PCR Detection System (Bio-Rad). Variant allele fraction was calculated as the ratio of variant allele to total MYD88 using the following formula as per the manufacturer's recommendations:
Evaluation of WM-Specific Immunity
The presence of WM-specific T cells in blood or bone marrow was quantified both directly ex vivo and after in vitro stimulation with autologous tumor-loaded DCs using a 16-hour Elispot assay as described previously (56). This assay allows for detection of global tumor-specific immunity against multiple tumor-derived antigens. Briefly, autologous T cells were cultured with unpulsed DCs (as controls) or CD19+ tumor-loaded DCs at a T:DC ratio of 1:10. After 16 hours, the presence of tumor-specific IFNγ spot-forming cells was quantified as described (56). For in vitro expansion, T cells were cultured with autologous tumor-loaded DCs for 10 to 14 days at a T:DC ratio of 1:30 in the presence of IL2 (Chiron) at 10 U/mL on days 4 and 7. In this assay, DCs are needed both for expansion of tumor-specific T cells and for readout of T-cell immunity.
Detection of MYD88-Specific Immunity
Thawed mononuclear cells were cultured overnight with peptides containing wild-type or mutant (L265P) MYD88 mutation sequences. After overnight culture, culture supernatant was harvested and analyzed for the presence of IP-10 as previously described (57).
Statistical analysis of mass cytometry data was performed using 2D graphing and statistics software GraphPad Prism. Nonparametric Mann–Whitney (for comparing two groups) and Kruskal–Wallis (for comparing three groups) tests with a significance threshold of P < 0.05 were used to compare different cohorts. Wilcoxon rank-sum test with a significance threshold of P < 0.05 after Bonferroni correction was used to identify differentially expressed genes between clusters and disease states in the single-cell RNA-sequencing data. We used χ2 with a significance threshold of P < 0.005 to identify clusters with differential composition by disease state. Data in bar graphs are plotted as mean ± SEM.
A.K. Nooka reports personal fees from Amgen, Takeda, Janssen, GlaxoSmithKline, Bristol Myers Squibb, Sanofi, Oncopeptides, Karyopharm, and Adaptive outside the submitted work. B.G. Barwick reports grants from Multiple Myeloma Research Foundation, American Society of Hematology, and Paula and Roger Riney Foundation during the conduct of the study. J.L. Kaufman reports grants from NIH during the conduct of the study, as well as grants and personal fees from Janssen, Bristol Myers Squibb, Roche, Amgen, Takeda, and Sanofi, personal fees from TG Therapeutics, Incycte, and Pharmacyclics, and grants from Fortis, Sutro, AbbVie, and Novartis outside the submitted work. S.M. Ansell reports grants from International Waldenstrom Macroglobulinemia Foundation outside the submitted work. L.H. Boise reports grants and personal fees from AstraZeneca, and personal fees from Genentech and AbbVie outside the submitted work. S. Lonial reports personal fees from Amgen, Takeda, Janssen, Bristol Myers Squibb, GlaxoSmithKline, Celgene, AbbVie, and Novartis, and personal fees and other support from TG Therapeutics outside the submitted work. M.V. Dhodapkar reports other support from Bristol Myers Squibb/Celgene, Roche/Genentech, Janssen, and Lava Therapeutics outside the submitted work. No disclosures were reported by the other authors.
A. Kaushal: Formal analysis, visualization, writing–review and editing. A.K. Nooka: Project administration, writing–review and editing. A.R. Carr: Writing–review and editing, performed mass cytometry. K.E. Pendleton: Writing–review and editing, performed single-cell sequencing. B.G. Barwick: Writing–review and editing, exome sequencing analysis. J. Manalo: Writing–review and editing, sample processing. S.S. McCachren: Formal analysis, writing–review and editing. V.A. Gupta: Writing–review and editing. N.S. Joseph: Writing–review and editing. C.C. Hofmeister: Writing–review and editing. J.L. Kaufman: Writing–review and editing. L.T. Heffner: Writing–review and editing. S.M. Ansell: Writing–review and editing. L.H. Boise: Writing–review and editing. S. Lonial: Writing–review and editing. K.M. Dhodapkar: Conceptualization, resources, formal analysis, supervision, funding acquisition, writing–original draft, project administration. M.V. Dhodapkar: Conceptualization, resources, data curation, formal analysis, supervision, funding acquisition, writing–original draft, project administration.
This work was supported in part by funds from the International Waldenstrom Macroglobulinemia Foundation and NIH R35CA197603 (to M.V. Dhodapkar). K.M. Dhodapkar is supported in part by funds from NIH CA238471 and AR077926. M.V. Dhodapkar and K.M. Dhodapkar are also supported in part by funds from a Specialized Center for Research (SCOR) award from the Leukemia & Lymphoma Society and from the Riney Foundation. The authors acknowledge the help of Deon Doxie and Roman Radzievski in Winship mass cytometry resource and colleagues in the Emory Integrated Genomics Core (both supported by P30CA138292 to Winship shared resources).