Abstract
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
Introduction
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.
Results
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 the bone marrow microenvironment comparing HDs, IgM MGUS, and WM. A, Overall strategy. Bone marrow mononuclear cells (BMMNC) were obtained from patients with IgM MGUS (n = 8) and WM (n = 8), as well as age-matched HDs (n = 5). CITE-seq, single-cell mass cytometry, BCR sequencing, and exome sequencing were performed on the samples. BMMNCs were also used for functional assays to test T-cell reactivity to tumor. B, BMMNCs from HDs (n = 5), MGUS (n = 8), and WM (n = 8) were stained with metal-conjugated antibodies. Data were analyzed using Cytobank software. Figure shows t-distributed Stochastic Neighbor Embedding (t-SNE) plots of concatenated live-cell gated data from mass cytometry analysis. Concatenations were done with equal cell numbers of cells from each donor. The overlay plots show differences in different immune cell subsets including differences in B-cell subsets [naïve B cells (brown) and CXCR5neg B cells in MGUS and WM (orange)], myeloid/monocyte population (green), and T cells (pink) in patients with MGUS (pink). NK, natural killer. C, BMMNCs from HDs (n = 4), MGUS (n = 7), and WM (n = 7) were labeled with TotalSeq-C antibodies and processed using the 10x DropSeq platform. Figure shows Uniform Manifold Approximation and Projection (UMAP) clustering of 84,128 single BMMNCs based on transcriptome. Forty-two distinct clusters could be identified, including B/lymphoplasmacytoid (LPC) cells (clusters 3, 4, 10, 14, 18, 22, 29, 34, 32, 37, 39), T/NK cells (0, 1, 6, 7, 9, 11, 12, 19, 20, 27, 31, 33, 36, 41), myeloid cells including monocytes and dendritic cell subsets (clusters 2, 5, 23, 25, 30), as well as progenitors/precursor cells (clusters 8, 13, 15, 16, 21, 24, 26, 28, 40). D, Feature plots showing surface expression of lineage antibodies CD3, CD56, CD19, and CD14 on clusters in C. E, UMAP clustering of BMMNCs cells by cohort. Figure shows differences in distinct myeloid (clusters 2, 5) and B-cell populations (cluster 3) in MGUS and WM, as well as in B/LPC populations in MGUS and WM (e.g., clusters 4, 10, 17, 22).
Changes in the bone marrow microenvironment comparing HDs, IgM MGUS, and WM. A, Overall strategy. Bone marrow mononuclear cells (BMMNC) were obtained from patients with IgM MGUS (n = 8) and WM (n = 8), as well as age-matched HDs (n = 5). CITE-seq, single-cell mass cytometry, BCR sequencing, and exome sequencing were performed on the samples. BMMNCs were also used for functional assays to test T-cell reactivity to tumor. B, BMMNCs from HDs (n = 5), MGUS (n = 8), and WM (n = 8) were stained with metal-conjugated antibodies. Data were analyzed using Cytobank software. Figure shows t-distributed Stochastic Neighbor Embedding (t-SNE) plots of concatenated live-cell gated data from mass cytometry analysis. Concatenations were done with equal cell numbers of cells from each donor. The overlay plots show differences in different immune cell subsets including differences in B-cell subsets [naïve B cells (brown) and CXCR5neg B cells in MGUS and WM (orange)], myeloid/monocyte population (green), and T cells (pink) in patients with MGUS (pink). NK, natural killer. C, BMMNCs from HDs (n = 4), MGUS (n = 7), and WM (n = 7) were labeled with TotalSeq-C antibodies and processed using the 10x DropSeq platform. Figure shows Uniform Manifold Approximation and Projection (UMAP) clustering of 84,128 single BMMNCs based on transcriptome. Forty-two distinct clusters could be identified, including B/lymphoplasmacytoid (LPC) cells (clusters 3, 4, 10, 14, 18, 22, 29, 34, 32, 37, 39), T/NK cells (0, 1, 6, 7, 9, 11, 12, 19, 20, 27, 31, 33, 36, 41), myeloid cells including monocytes and dendritic cell subsets (clusters 2, 5, 23, 25, 30), as well as progenitors/precursor cells (clusters 8, 13, 15, 16, 21, 24, 26, 28, 40). D, Feature plots showing surface expression of lineage antibodies CD3, CD56, CD19, and CD14 on clusters in C. E, UMAP clustering of BMMNCs cells by cohort. Figure shows differences in distinct myeloid (clusters 2, 5) and B-cell populations (cluster 3) in MGUS and WM, as well as in B/LPC populations in MGUS and WM (e.g., clusters 4, 10, 17, 22).
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).
Mass cytometry analysis of changes in B cells. Single-cell mass cytometry was performed on bone marrow mononuclear cells from HDs (n = 5) or patients with IgM MGUS (n = 8) and WM (n = 8). A, T-distributed Stochastic Neighbor Embedding (t-SNE) plot gated on CD19+ B cells by cohort. Figure shows differences in different B-cell subsets, including differences in naïve B cells (green) and switched memory B cells in MGUS and WM (brown) and in CXCR5neg B-cell subset in MGUS and WM patients compared with HDs (purple). B, Phenotypes of major subpopulations of CD19+ cells in A. Heatmap shows expression of cell-surface markers (CXCR5, CD27, CD38, CD21, CD22, HLA-DR, CD40, CXCR4, C-KIT, Ki67, CD79b, and CD20) and immunoglobulins (IgD, IgM) in the different CD19+ cell populations shown in A. Histogram shows expression of CD19 on the various B-cell populations and T cells as control. C, Differences in CXCR5-positive B cells. Graphs show the percentage of CXCR5-positive cells in HDs, MGUS, and WM. Each dot represents a unique patient. D, Differences in naïve B cells. Figure shows naïve B cells as a percentage of total B cells in HDs, MGUS, and WM. Each dot represents a unique patient. E, Differences in B cells based on expression of clonal light chain. CD19+ B cells from MGUS and WM were gated based on the presence of clonal Ig light chain (kappa/lambda, depending on the light chain of the M spike). Figure shows differences in B-cell populations including in naïve B cells (right) and CXCR5+ B cells (left) in both clonal light chain positive (clonal LC+) as well as B cells expressing the opposite light chain (termed as nonclonal LC+) in MGUS and WM. #, P = 0.06; *, P < 0.05; **, P < 0.01, Kruskal–Wallis.
Mass cytometry analysis of changes in B cells. Single-cell mass cytometry was performed on bone marrow mononuclear cells from HDs (n = 5) or patients with IgM MGUS (n = 8) and WM (n = 8). A, T-distributed Stochastic Neighbor Embedding (t-SNE) plot gated on CD19+ B cells by cohort. Figure shows differences in different B-cell subsets, including differences in naïve B cells (green) and switched memory B cells in MGUS and WM (brown) and in CXCR5neg B-cell subset in MGUS and WM patients compared with HDs (purple). B, Phenotypes of major subpopulations of CD19+ cells in A. Heatmap shows expression of cell-surface markers (CXCR5, CD27, CD38, CD21, CD22, HLA-DR, CD40, CXCR4, C-KIT, Ki67, CD79b, and CD20) and immunoglobulins (IgD, IgM) in the different CD19+ cell populations shown in A. Histogram shows expression of CD19 on the various B-cell populations and T cells as control. C, Differences in CXCR5-positive B cells. Graphs show the percentage of CXCR5-positive cells in HDs, MGUS, and WM. Each dot represents a unique patient. D, Differences in naïve B cells. Figure shows naïve B cells as a percentage of total B cells in HDs, MGUS, and WM. Each dot represents a unique patient. E, Differences in B cells based on expression of clonal light chain. CD19+ B cells from MGUS and WM were gated based on the presence of clonal Ig light chain (kappa/lambda, depending on the light chain of the M spike). Figure shows differences in B-cell populations including in naïve B cells (right) and CXCR5+ B cells (left) in both clonal light chain positive (clonal LC+) as well as B cells expressing the opposite light chain (termed as nonclonal LC+) in MGUS and WM. #, P = 0.06; *, P < 0.05; **, P < 0.01, Kruskal–Wallis.
CITE-seq analysis of differences in CD19+ cells. B cells were identified based on their surface binding of anti-CD19 antibody, and antibody-based UMAP clustering analysis was performed on the B-cell–associated antibodies used in CITE-seq (IgD, CD27, CXCR5, CD138, PD-1, CD10, CD16, CD21, CD79b, CD27, CD38, and CD20). B cells from patients with IgM MGUS (4,039 cells) and WM (13,077 cells) were compared with those from HDs (3,403 cells). A, UMAP of B cells clustered based on binding to antibodies. Figure shows distribution of B cells in 10 distinct clusters within the three cohorts, showing differences in IgD+CD27− naïve B cells (cluster 2) and in CXCR5neg B cells (cluster 3) in MGUS and WM. B, Feature plot showing binding of antibodies against IgD, CD27, CD10, and CXCR5 in B cells within the different clusters. C, Proportion of cell clusters by cohort. Bar graph shows distribution of B cells from HDs, MGUS, and WM within different clusters shown in A. MGUS and WM cohorts had decreased proportion of cells in cluster 2 and a relative increase in cluster 3 when compared with HD B cells. Graph shows mean ± SEM. *, P < 0.05; **, P < 0.01, Kruskal–Wallis. D, UMAP of CD19+ B cells based on transcriptome. UMAP clustering analysis was performed on all B cells (n = 20,519) from HDs, MGUS, and WM. This clustering analysis revealed 24 transcriptionally distinct B-cell populations. E, UMAP of CD19+ B cells based on transcriptome by cohort. Figure shows decline in B cells in cluster 1 in MGUS and WM, whereas this is the dominant cluster in HDs. Transcriptome-based UMAP also revealed unique (patient-specific) B-cell clusters (circled) in WM. F, Heatmap of patient-specific clusters of clonal B cells. Comparison of transcriptomes of unique B-cell clusters identified in WM revealed distinct patient-specific patterns of gene expression. Data shown in the heatmap include top 10 differentially expressed genes from each B-cell cluster.
CITE-seq analysis of differences in CD19+ cells. B cells were identified based on their surface binding of anti-CD19 antibody, and antibody-based UMAP clustering analysis was performed on the B-cell–associated antibodies used in CITE-seq (IgD, CD27, CXCR5, CD138, PD-1, CD10, CD16, CD21, CD79b, CD27, CD38, and CD20). B cells from patients with IgM MGUS (4,039 cells) and WM (13,077 cells) were compared with those from HDs (3,403 cells). A, UMAP of B cells clustered based on binding to antibodies. Figure shows distribution of B cells in 10 distinct clusters within the three cohorts, showing differences in IgD+CD27− naïve B cells (cluster 2) and in CXCR5neg B cells (cluster 3) in MGUS and WM. B, Feature plot showing binding of antibodies against IgD, CD27, CD10, and CXCR5 in B cells within the different clusters. C, Proportion of cell clusters by cohort. Bar graph shows distribution of B cells from HDs, MGUS, and WM within different clusters shown in A. MGUS and WM cohorts had decreased proportion of cells in cluster 2 and a relative increase in cluster 3 when compared with HD B cells. Graph shows mean ± SEM. *, P < 0.05; **, P < 0.01, Kruskal–Wallis. D, UMAP of CD19+ B cells based on transcriptome. UMAP clustering analysis was performed on all B cells (n = 20,519) from HDs, MGUS, and WM. This clustering analysis revealed 24 transcriptionally distinct B-cell populations. E, UMAP of CD19+ B cells based on transcriptome by cohort. Figure shows decline in B cells in cluster 1 in MGUS and WM, whereas this is the dominant cluster in HDs. Transcriptome-based UMAP also revealed unique (patient-specific) B-cell clusters (circled) in WM. F, Heatmap of patient-specific clusters of clonal B cells. Comparison of transcriptomes of unique B-cell clusters identified in WM revealed distinct patient-specific patterns of gene expression. Data shown in the heatmap include top 10 differentially expressed genes from each B-cell cluster.
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 myeloid cells. Bone marrow mononuclear cells that did not bind to either anti-CD3 or anti-CD19 antibody from HDs (11,160 cells), MGUS (10,376 cells), and WM (11,262 cells) were analyzed using transcriptome-based UMAP clustering to evaluate changes in myeloid cells, NK cells, and progenitors. A, UMAP of CD3−CD19− cells based on transcriptome expression revealed 24 different clusters including myeloid/monocyte clusters (0, 1, 7, 13, 14, 17), NK cell clusters (5, 9, 19, 22), and progenitor/precursor cell clusters (2, 3, 4, 6, 8, 10, 12, 15, 16, 18, 21, 23). B, Feature plot showing antibody binding for CD14, CD11c, CD16, CD56, and C-kit as well as hemoglobin expression (by transcript). C, UMAP of CD3−CD19− cells by cohort showing decreased proportion of cells in myeloid clusters 0 and 1 in WM. D, Proportion of myeloid clusters by cohort. Bar graph showing the proportions of all myeloid clusters in HDs, MGUS, or WM, with decline in cluster 1 and increase in cluster 7 in MGUS and WM compared with HD. E, Volcano plot of genes differentially expressed in myeloid cluster 7. Many of the genes overexpressed in this cluster include those associated with myeloid inflammation. F, Pathway analysis of differentially regulated genes enriched in myeloid cluster 7. ER, endoplasmic reticulum; TCR, T-cell receptor.
Changes in myeloid cells. Bone marrow mononuclear cells that did not bind to either anti-CD3 or anti-CD19 antibody from HDs (11,160 cells), MGUS (10,376 cells), and WM (11,262 cells) were analyzed using transcriptome-based UMAP clustering to evaluate changes in myeloid cells, NK cells, and progenitors. A, UMAP of CD3−CD19− cells based on transcriptome expression revealed 24 different clusters including myeloid/monocyte clusters (0, 1, 7, 13, 14, 17), NK cell clusters (5, 9, 19, 22), and progenitor/precursor cell clusters (2, 3, 4, 6, 8, 10, 12, 15, 16, 18, 21, 23). B, Feature plot showing antibody binding for CD14, CD11c, CD16, CD56, and C-kit as well as hemoglobin expression (by transcript). C, UMAP of CD3−CD19− cells by cohort showing decreased proportion of cells in myeloid clusters 0 and 1 in WM. D, Proportion of myeloid clusters by cohort. Bar graph showing the proportions of all myeloid clusters in HDs, MGUS, or WM, with decline in cluster 1 and increase in cluster 7 in MGUS and WM compared with HD. E, Volcano plot of genes differentially expressed in myeloid cluster 7. Many of the genes overexpressed in this cluster include those associated with myeloid inflammation. F, Pathway analysis of differentially regulated genes enriched in myeloid cluster 7. ER, endoplasmic reticulum; TCR, T-cell receptor.
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.
Changes in T cells. A, T-distributed Stochastic Neighbor Embedding (t-SNE) plot of changes in T cells by cohort. Bone marrow mononuclear cells from HDs, MGUS, and WM were analyzed using mass cytometry. CD3 T cells from HDs, MGUS, and WM cohorts were concatenated using equal numbers of cells from individual patients and analyzed using t-SNE analysis. Overlay plot shows distribution of T cells across different cohorts with progressive increase in CD8 granzyme+ T cells (red) in MGUS/WM and decline in naïve CD8+ T cells (brown) compared with HDs. TM, memory T cell; TN, naïve T cell. B, Phenotype of major CD4/CD8+ T-cell subsets. Heat map shows protein expression of cell-surface markers as well as transcription factors/lytic molecules in the different T-cell subsets shown in A. The CD8+ granzyme+ T cells enriched in WM have a phenotype of TCF1loTigit+ KLRG1+ T cells. C, Proportions of CD8 T-cell subsets in HDs and WM. Supplementary figure shows naïve (TN; CCR7, CD45RO−), effector memory (TEM; CCR7−CD45RO+), central memory (TCM; CCR7−, CD45RO+CD45RA+), and terminal effector (TERM EFF; CCR7−CD45RO−) cells as a percentage of total CD8 T cells in HDs and WM patients. Each dot represents a unique patient. *, P < 0.05, Kruskal–Wallis. D, Proportion of granzyme-positive CD8 T cells as a proportion of total CD8 T cells in HDs as well as MGUS and WM patients. Each dot represents a unique patient. *, P < 0.05, Kruskal–Wallis. E, Volcano plot of differentially expressed genes in CD8+ T cells in MGUS versus HDs based on CITE-seq. F, Volcano plot of differentially expressed genes in CD8+ T cells WM versus HD based on CITE-seq.
Changes in T cells. A, T-distributed Stochastic Neighbor Embedding (t-SNE) plot of changes in T cells by cohort. Bone marrow mononuclear cells from HDs, MGUS, and WM were analyzed using mass cytometry. CD3 T cells from HDs, MGUS, and WM cohorts were concatenated using equal numbers of cells from individual patients and analyzed using t-SNE analysis. Overlay plot shows distribution of T cells across different cohorts with progressive increase in CD8 granzyme+ T cells (red) in MGUS/WM and decline in naïve CD8+ T cells (brown) compared with HDs. TM, memory T cell; TN, naïve T cell. B, Phenotype of major CD4/CD8+ T-cell subsets. Heat map shows protein expression of cell-surface markers as well as transcription factors/lytic molecules in the different T-cell subsets shown in A. The CD8+ granzyme+ T cells enriched in WM have a phenotype of TCF1loTigit+ KLRG1+ T cells. C, Proportions of CD8 T-cell subsets in HDs and WM. Supplementary figure shows naïve (TN; CCR7, CD45RO−), effector memory (TEM; CCR7−CD45RO+), central memory (TCM; CCR7−, CD45RO+CD45RA+), and terminal effector (TERM EFF; CCR7−CD45RO−) cells as a percentage of total CD8 T cells in HDs and WM patients. Each dot represents a unique patient. *, P < 0.05, Kruskal–Wallis. D, Proportion of granzyme-positive CD8 T cells as a proportion of total CD8 T cells in HDs as well as MGUS and WM patients. Each dot represents a unique patient. *, P < 0.05, Kruskal–Wallis. E, Volcano plot of differentially expressed genes in CD8+ T cells in MGUS versus HDs based on CITE-seq. F, Volcano plot of differentially expressed genes in CD8+ T cells WM versus HD based on CITE-seq.
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.
Tumor-specific immunity in WM. A, Detection of WM-specific immunity in freshly isolated T cells from blood/bone marrow. Freshly isolated blood or marrow-derived T cells were cultured overnight in the presence of autologous monocyte-derived mature dendritic cells (DC) loaded with autologous tumor cells (or unpulsed mature DCs as control). Tumor-specific IFNγ-producing T cells were quantified with an Elispot assay. Each dot represents mean of replicates from an individual patient. PBMC, peripheral blood mononuclear cell. DC+ Tum, DCs loaded with tumor cells. B, Ex vivo expansion of WM-specific T cells by autologous tumor-loaded DCs. Blood or bone marrow T cells were stimulated for 10 to 14 days in the presence of autologous unloaded mature DCs [DC(-)], DCs loaded with CD19− cells (nontumor, DC NT), or tumor-loaded DCs (DC Tum). The presence of tumor-specific T cells was quantified using IFNγ Elispot following overnight culture with autologous unpulsed or tumor-loaded DCs. Each dot represents mean of replicates from an individual patient. APC, antigen-presenting cell. C, CD4 T-cell depletion: Responder T cells from the experiment described in B were depleted of CD4+ T cells utilizing magnetic beads prior to testing for reactivity using tumor-loaded DCs. The presence of tumor-specific T cells was quantified using IFNγ Elispot. Each dot represents mean of replicates from an individual patient. D and E, Detection of MYD88 L265P–specific reactivity in WM patients. Mononuclear cells were cultured with peptides spanning wild-type (WT) or mutant (Mut) MyD88 L265P sequences. The presence of peptide reactivity was assayed based on the detection of IP-10/CXCL-10 in the supernatant. D, Data from representative patients with or without MYD88 reactivity. E, Summary of MyD88 L265P reactivity in all patients and HDs studied.
Tumor-specific immunity in WM. A, Detection of WM-specific immunity in freshly isolated T cells from blood/bone marrow. Freshly isolated blood or marrow-derived T cells were cultured overnight in the presence of autologous monocyte-derived mature dendritic cells (DC) loaded with autologous tumor cells (or unpulsed mature DCs as control). Tumor-specific IFNγ-producing T cells were quantified with an Elispot assay. Each dot represents mean of replicates from an individual patient. PBMC, peripheral blood mononuclear cell. DC+ Tum, DCs loaded with tumor cells. B, Ex vivo expansion of WM-specific T cells by autologous tumor-loaded DCs. Blood or bone marrow T cells were stimulated for 10 to 14 days in the presence of autologous unloaded mature DCs [DC(-)], DCs loaded with CD19− cells (nontumor, DC NT), or tumor-loaded DCs (DC Tum). The presence of tumor-specific T cells was quantified using IFNγ Elispot following overnight culture with autologous unpulsed or tumor-loaded DCs. Each dot represents mean of replicates from an individual patient. APC, antigen-presenting cell. C, CD4 T-cell depletion: Responder T cells from the experiment described in B were depleted of CD4+ T cells utilizing magnetic beads prior to testing for reactivity using tumor-loaded DCs. The presence of tumor-specific T cells was quantified using IFNγ Elispot. Each dot represents mean of replicates from an individual patient. D and E, Detection of MYD88 L265P–specific reactivity in WM patients. Mononuclear cells were cultured with peptides spanning wild-type (WT) or mutant (Mut) MyD88 L265P sequences. The presence of peptide reactivity was assayed based on the detection of IP-10/CXCL-10 in the supernatant. D, Data from representative patients with or without MYD88 reactivity. E, Summary of MyD88 L265P reactivity in all patients and HDs studied.
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.
Exome sequencing of early B-cell progenitors and proposed model for WM evolution. A–C, MYD88 mutations are present throughout the hematopoietic system. A, VAFs in CD19+CD10− cells (y-axis) as compared with CD19+CD10+ pre–B cells (left), CD19−CD34− cells (middle), and CD19−CD34+ cells (right; x-axis). All variant alleles are shown, and those with significantly different allelic fractions (FDR ≤ 0.05) between any two cell fractions are shown in red, with select genes labeled in colored triangles. B, Venn diagrams of variants shared and distinct between cell fractions. C, VAF of MYD88L265P across CD19−CD34+, CD19−CD34−, CD19+CD10+, and CD19+CD10− cell fractions. D, Proposed model for WM development: Acquisition of MYD88 mutation in hematopoietic progenitors is an early event in the origin of WM and associated with several changes in the immune microenvironment, including increase in extrafollicular B cells, myeloid inflammation, and alteration in immune function that begin as early as the MGUS stage. The B-cell clone emerges in this milieu and undergoes progressive growth and evolution in the WM stage. Host immune system mediates tumor-specific recognition of the clone but undergoes immune exhaustion over time.
Exome sequencing of early B-cell progenitors and proposed model for WM evolution. A–C, MYD88 mutations are present throughout the hematopoietic system. A, VAFs in CD19+CD10− cells (y-axis) as compared with CD19+CD10+ pre–B cells (left), CD19−CD34− cells (middle), and CD19−CD34+ cells (right; x-axis). All variant alleles are shown, and those with significantly different allelic fractions (FDR ≤ 0.05) between any two cell fractions are shown in red, with select genes labeled in colored triangles. B, Venn diagrams of variants shared and distinct between cell fractions. C, VAF of MYD88L265P across CD19−CD34+, CD19−CD34−, CD19+CD10+, and CD19+CD10− cell fractions. D, Proposed model for WM development: Acquisition of MYD88 mutation in hematopoietic progenitors is an early event in the origin of WM and associated with several changes in the immune microenvironment, including increase in extrafollicular B cells, myeloid inflammation, and alteration in immune function that begin as early as the MGUS stage. The B-cell clone emerges in this milieu and undergoes progressive growth and evolution in the WM stage. Host immune system mediates tumor-specific recognition of the clone but undergoes immune exhaustion over time.
Discussion
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.
Methods
Patients
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.
Mass Cytometry
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.
CITE-seq
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.
Cell Sorting
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 (v4.1.9.0; ref. 53). Variant alleles were called with Mutect2 (v4.1.9.0; 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
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.
Authors' Disclosures
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.
Authors' Contributions
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.
Acknowledgments
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).