Transcriptome sequencing reveals a profile that corresponds to genomic variants in Waldenström macroglobulinemia

Whole-genome sequencing identified several highly recurring somatic mutations in patients with Waldenström macroglobulinemia (WM).^1,2  In over 90% of WM patients, a single point mutation at NM_002468:c.978T>C (rs387907272) in MYD88 is found, resulting in a p.Leu265Pro (L265P) amino acid change.^1,3  MYD88 is an adaptor for Toll-like (TLR) and interleukin 1 (IL1) receptors, and the MYD88^L265P mutation triggers constitutive activation of NF-κB through IRAK and BTK.^1,4,5  Mutated MYD88 WM patients show greater overall survival and clinical responses to the BTK inhibitor ibrutinib.^6,7 

Activating CXCR4 frameshift or nonsense mutations in the C-terminal tail are found in 30% to 40% of WM patients, are primarily subclonal, and almost always associated with MYD88^L265P.^2,3,8  These somatic mutations are similar to the causal germ line variants that underlie WHIM (warts, hypogammaglobulinemia, infection, and myelokathexis) syndrome.^2,9  In WM, somatic CXCR4 mutations (CXCR4^WHIM) are determinants of disease presentation, as well as resistance to ibrutinib.^3,6,10  Somatic mutations in ARID1A, a member of the SWItch/sucrose nonfermentable family and epigenetic regulator, are also found in 20% of WM patients. Gene losses affecting NF-κB signaling (ie, HIVEP2, TNFAIP3), as well as the ARID1A homolog ARID1B are present in most WM patients.^2,11  Deletions in chromosome 6q with and without concurrent 6p amplifications, trisomy 3, and amplifications of 3q, as well as trisomy 4, are also commonly found in WM.^12-14  Previous array-based gene expression studies of WM were largely conducted prior to these genomic discoveries and therefore the effects of recurrent somatic events on transcriptional regulation remain to be clarified.^15-17  We therefore performed next-generation RNA sequencing in 57 WM patients and compared findings to sorted healthy donor-derived nonmemory (CD19⁺CD27⁻) and memory (CD19⁺CD27⁺) B cells. The latter represent the B-cell population from where most cases of WM are thought to be derived.^18,19 

Bone marrow (BM) aspirates were collected from 57 patients with the WM consensus diagnosis.²⁰  Participants provided informed consent for sample collection per the Dana-Farber/Harvard Cancer Center Institutional Review Board. WM cells were isolated by CD19⁺ magnetic-activated cell sorting (MACS) microbead selection (Miltenyi Biotec, Auburn, CA) from Ficoll-Paque (Amersham-Pharmacia Biotech, Piscataway, NJ) separated BM mononuclear cells. Peripheral blood mononuclear cells from nine healthy donors (HDs) were sorted for nonmemory (CD19⁺CD27⁻) and memory B cells (CD19⁺CD27⁺) using a memory B-cell isolation kit (Miltenyi Biotec). RNA and DNA were purified using the AllPrep mini kit (Qiagen, Valencia, CA). Most samples were previously characterized by whole-genome sequencing and all samples were screened for MYD88 and CXCR4 gene mutations by Sanger sequencing.²  MYD88^L265P and CXCR4 c.1013C>G and c.1013C>A mutations were analyzed by allele-specific polymerase chain reaction (PCR) as previously described.^8,21 

Transcriptome profiling was conducted by the Center for Cancer Computational Biology at the Dana-Farber Cancer Institute (Boston, MA) using the NEBNext Ultra RNA library prep kit (New England BioLabs, Ipswich, MA). The paired-end samples were run 2 per lane for 50 cycles on an Illumina HiSeq (Illumina, San Diego, CA). Read-level data are available through dbGAP accession (applied). Reads were aligned to KnownGene HG19/GRCh37 reference using STAR (Spliced Transcripts Alignment to a Reference).²²  Genes with mean raw read counts of <10 were not analyzed, leaving 16 888 expressed genes for analysis. Read counts per gene were obtained using featureCounts from Rsubread, and analyzed using voom from the edgeR/limma Bioconductor packages in R (R Foundation for Statistical Computing, Vienna, Austria).^23-27  Differential expression models accounted for sex, prior treatment, as well as MYD88 and CXCR4 mutation status. A false discovery rate (FDR) cutoff of 10% was used to determine significant differentially expressed genes. Functional enrichment analysis was conducted using Ingenuity Pathway Analysis (Qiagen). Clustering and correlation analysis was conducted using the variance stabilizing transformation of the count data from the Bioconductor DESeq2 package.²⁸  In all other cases, estimates of gene expression levels are represented in transcripts per million (TpM). Log-fold change (LFC) listed in text and tables are derived from the limma analysis.

Previously described BCWM.1 and MWCL-1 cell lines transduced to express CXCR4 with or without activating mutations observed in WM patients were used to model CXCR4-stimulated gene expression changes in WM.¹⁰  Briefly, CXCR4 complementary DNA (cDNA) transcripts were subcloned into plenti-internal ribosomal entry site (IRES)–green fluorescent protein (GFP) vectors, and stably transduced using a lentiviral system.⁴  In addition to wild-type (WT) CXCR4, vectors with the mutations c.1013C>G (p.Ser338*), c.932_933insT (p.Thr311fs), and c.1030_1041delinsGT (p.Ser344fs) were generated. Cell lines were stimulated with 50 nM of the CXCR4 ligand CXCL12 (SDF1A) (R&D Systems, Minneapolis MN) for 2 hours. RNA from each sample was extracted at baseline and at 2 hours. Gene expression was assessed using Affymetrix Human Gene ST 2.0 arrays (Affymetrix, Santa Clara, CA). Array data are archived in the Gene Expression Omnibus database (accession number GSE83150). Array data were analyzed using the Oligo and limma R bioconductor packages.^24,29 

Random forest analysis was implemented using the cforest function in the party package in R.³⁰  Models were constructed using the cforest_unbiased function with 20 predictors per tree over 10 000 trees. As variable importance scores of nonsignificant genes are distributed around zero, genes were further filtered for those with permutation importance scores greater than the absolute value of the lowest scoring predictor. For reproducibility, the random seed for the analysis was set to 265. Principal component analysis (PCA) was conducted using R’s prcomp, Pearson’s product-moment coefficient was used for correlation studies, and cohort mean values were compared using the Wilcoxon rank sum test.

Gene expression results were validated by using TaqMan gene expression assays Hs01027785_m1 (DUSP4), hs00244839_m1 (DUSP5), Hs00243182_m1 (RGS13), Hs00892674_m1 (RGS16), Hs00698249_m1 (RASSF6), Hs00396602_m1 (WNK2), Hs00924602_m1 (PRDM5) (Thermo Fisher Scientific, Waltham MA). Methylation-specific PCR assays for RASSF6, WNK2, and PRDM5 were conducted on bisulfite-converted DNA using previously established protocols.^31-33 

The clinical characteristics for the 57 WM patients are presented in Table 1. The distribution of somatic mutations and cytogenetic findings in these patients was similar to those previously described.³⁴  The top 500 differentially expressed genes for each comparison discussed below are available as a part of the supplemental Data (available on the Blood Web site).

Comparison of WM and HD samples generated a list of 13 571 differentially expressed genes. To gain further insight into these gene expression differences, this gene list was sorted by absolute log₂ fold change, revealing a marked upregulation of DNTT, RAG1, RAG2, and IGLL1 in WM samples. These genes have important functional roles in VDJ recombination and B-cell receptor (BCR) signaling. Other upregulated genes relevant to WM biology and pathogenesis included CXCL12, VCAM1, CD5L, BMP3, and IGF1. Among the top 40 differentially expressed genes sorted by absolute expression in WM as measured in TpM were CXCR4, the WM prognostic marker B2M, and 2 genes related to BCR signaling: CD79A and CD79B. To examine transcriptional similarities between the samples in an unbiased manner, multidimensional scaling of the 500 genes with the highest variance was used to best represent the data in 2 dimensions (Figure 1A). Although peripheral B cells (PBs) and HD memory B cells (MBs) clustered together, separation between MYD88^L265PCXCR4^WHIM and MYD88^L265PCXCR4^WT genotyped WM patient samples was observed. To examine whether the patterns of expression were similar between HD samples and patient mutational groups, a correlation coefficient based on overall gene expression was generated for each sample-to-sample comparison. These values were averaged together by HD PB/MB cell type and WM MYD88/CXCR4 genotype to show the extent of correlation between groups (Figure 1B). Strong correlations were observed across all subsets, though mean correlation values between HD B cells and MYD88^L265PCXCR4^WHIM patients were significantly higher vs those between HD B cells and MYD88^L265PCXCR4^WT patients (P = .025 and P = .036 for PB and MB cells, respectively).

No preferential association was noted between WM samples and HD B-cell type using PCA of the 596 genes differentially expressed between HD PB and MB cells (Figure 1C; supplemental Figure 1). To better understand the role of B-cell differentiation in WM, we analyzed 19 genes linked to the transition from MB to plasmablasts and plasma cells (Figure 1D).^35,36  Many of these genes were not only significantly differentially expressed between HD and WM samples, but between WM patients based on MYD88 and CXCR4 mutation status suggesting reduced differentiation in both MYD88^L265PCXCR4^WHIM and, in particular, MYD88^WTCXCR4^WT samples. Separation of patient samples by genotype can be seen in the PCA of these 19 genes in Figure 1E. BCL2 and BCL2L1 were upregulated in all WM samples, however, analysis of the larger BCL2 family revealed stratification of samples by HD cell type and MYD88/CXCR4 mutation status (supplemental Figure 2). The proapoptotic BAX gene was downregulated in WM samples and PMAIP1 was uniquely upregulated in MYD88^L265PCXCR4^WT WM.

Restricting the analysis to WM patients, PCA of the top 500 high variance genes accounted for 41% of the overall variability within the first 2 components (Figure 2A). Tumor suppressors including WNK2, TP53I11, PRDM5, and CABLES1 influence the second principal component that separates MYD88^L265PCXCR4^WT from MYD88^L265PCXCR4^WHIM and MYD88^WTCXCR4^WT samples (supplemental Figure 1). Supervised clustering of the 3103 genes differentially expressed between MYD88^L265PCXCR4^WT and MYD88^L265PCXCR4^WHIM revealed that most expression differences in MYD88^L265PCXCR4^WHIM patients followed a pattern resembling HD samples, despite carrying the MYD88^L265P mutation (Figure 2B). The top predicted upstream regulator for this gene signature was the inhibition of lipopolysaccharide signaling (Figure 2C). To further explore how CXCR4 mutations modulate MYD88 signaling in WM, the gene signature was filtered for genes involved in TLR signaling. MYD88^L265PCXCR4^WHIM patients showed downregulation of the lipopolysaccharide receptors TLR4 and NOD2, and upregulation of TLR7 and IRAK3 (supplemental Figure 3A). The most significant results associated with CXCR4 mutation status were silencing of genes in patients with MYD88^L265PCXCR4^WHIM that were elevated in MYD88^L265PCXCR4^WT patients (supplemental Figure 3B).

To investigate the role of methylation in genes differentially expressed based on CXCR4 and MYD88 mutation status, 3 such genes were selected for screening based on their regulation by promoter methylation in other malignancies: PRDM5, WNK2, and RASSF6.^32,33,37  Gene expression was validated by PCR in HD B-cell samples as well as in MYD88 and CXCR4 genotyped WM patient samples (supplemental Figure 4A). Mirroring the transcriptome data, RASSF6 expression was elevated in MYD88^L265PCXCR4^WT and MYD88^L265PCXCR4^WHIM samples but not in HD or MYD88^WT samples whereas WNK2 and PRDM5 were elevated only in the MYD88^L265PCXCR4^WT patients. To determine relevance of these methylation sites in hematologic malignancies, WM, B-cell lymphoma, and myeloma cell lines were assessed for promoter methylation using methylation-specific PCR assays (supplemental Figure 4B). Sixteen MYD88 and CXCR4 genotyped primary WM patient samples and 4 healthy donor B-cell samples were then tested (supplemental Figure 4C). Robust promoter methylation was observed for WNK2 and PRDM5 in the MYD88^WTCXCR4^WT patients but was absent in HD samples. Partial promoter methylation of WNK2 and PRDM5 was also observed in some MYD88^L265P samples. This correlated to a relative reduction in gene expression only for WNK2 in the MYD88^L265PCXCR4^WHIM genotyped patient sample with the strongest methylation signal.

To better understand interactions between CXCR4^WHIM and MYD88^L265P in WM, lentiviral transduction was performed to overexpress CXCR4^WT, or 1 of 3 documented CXCR4^WHIM transcripts (c.1013C>G, c.932_933insT, or c.1030_1041delinsGT) in MYD88^L265P expressing BCWM.1 WM cells. Transduced cells were stimulated in triplicate with CXCL12 for 2 hours and gene expression profiling performed. The relative change in gene expression from baseline was then compared between CXCR4^WT and CXCR4^WHIM transduced cells revealing 61 differentially expressed genes. Of these genes, 16 of 61 (26.2%) were encompassed in the gene expression signature for MYD88^L265PCXCR4^WHIM patient samples (supplemental Table 1). In 8 of 16 (50%) of these genes, the direction of the change relative to MYD88^L265PCXCR4^WT was same as that observed in the patient samples. However, it was the genes whose direction of change was opposed that revealed the most about the underlying biology. These discordant genes were enriched for suppressors of MAPK and G-protein signaling downstream of CXCR4. They were preferentially induced in CXCR4^WHIM cell lines by CXCL12 stimulation and suppressed in MYD88^L265PCXCR4^WHIM patient samples.^10,38  Interrogation of these dual specificity phosphatases (DUSP) and regulator of G-protein signaling (RGS) genes revealed RGS1, RGS2, RGS13, DUSP1, DUSP2, DUSP4, DUSP5, DUSP10, DUSP16, and DUSP22 were all significantly suppressed in MYD88^L265PCXCR4^WHIM vs MYD88^L265PCXCR4^WT patients. Expression in patients of DUSP4, DUSP5, RGS13, and RGS16, which were selected for validation based on the cell line studies, can be seen in supplemental Figure 5A. The increased expression of these genes in CXCR4^WHIM compared with CXCR4^WT transduced lines was validated by real-time PCR (supplemental Figure 5B). Dose-response curves for these genes following CXCL12 stimulation was conducted in transduced BCWM.1 and MWCL-1 (supplemental Figure 5C).

For the 5 patients with ARID1A mutations, 16 differentially expressed genes were observed. The ARID1A-mutated population also demonstrated elevated BM disease infiltration (median, 90%; range, 70%-95%) when compared with non-ARID1A-mutated patients (median, 58%; range, 5%-95%; P = .0259). No statistically significant gene expression changes were observed in the 4 patients with CD79B mutations. Analysis based on somatic cytogenetic abnormalities included deletion of chromosome 6q (131 genes), amplification of chromosome 6p (65 genes), amplification of chr3q including trisomy 3 (11 genes), and trisomy 4 (776 genes). A distinct gene expression signature for 263 genes was also observed in the 3 patients who had first-degree relatives with WM. No statistically significant differences in gene expression were observed based on presence of extramedullary disease, elevated B₂M levels (>3 mg/dL), or International Waldenström Macroglobulinemia Prognostic Scoring System score.

To identify genes associated with relevant WM clinical parameters, all genes were tested for correlation with serum immunoglobulin M (IgM) levels, hemoglobin, and BM disease involvement. After filtering the results by significance (P < .05) and the absolute value of the correlation estimate >0.3, 312 genes were found to be associated with IgM, 965 with hemoglobin, and 3738 with BM involvement. These genes were combined with MYD88/CXCR4 genotype, age, sex, and treatment status in a random forest regression analysis for the first 37 study patients, leaving the remaining 20 patients for cross-validation. Filtering genes/predictors by variable importance score resulted in 40, 236, and 143 associated genes for IgM, BM, and hemoglobin, respectively (supplemental Figure 6). Functional analysis was conducted on each gene list and top genes from each group were selected based on biological significance. The 40 genes associated with IgM were predicted to be downstream of IL6 signaling (P = .0018) and included many genes of WM pathogenic interest including TNFRSF17, BMP3, IRS1, and CEPZOS (Figure 3A). Many genes relevant to WM biology, including CXCL13, TP53, CXCR4, MYD88, CDC23, and AKAP1, were associated with BM disease involvement (Figure 3B). Although MYD88 and CXCR4 were both correlated to BM disease involvement (ρ = −0.50; P < .0001 and ρ = 0.46; P = .0003, respectively), they were negatively correlated with each other (ρ = −0.46; P = .0004). As shown in supplemental Figure 7, this relationship may be influenced by MD88/CXCR4 mutation status. The gene list associated with hemoglobin included CXCL13, as well as PIK3AP1, PIK3CD, AJUBA, and OS9 (Figure 3C). Based on these findings, CXCL13 was included in an independent serum cytokine profiling study of 86 WM patients that validated the CXCL13 correlation between BM (Spearman ρ = 0.562; P = 2.192 × 10−⁸) and hemoglobin (Spearman ρ = −0.567; P = 1.53 × 10⁻⁸).³⁹ 

The percentage of reads supporting the mutant allele relative to WT allele was calculated for MYD88, CXCR4, ARID1A, and CD79B in the RNA sequencing data. Although the relevant mutations were observed in the RNA in all cases, many samples with CXCR4^WHIM mutations had mutant allele burdens in excess of 50%, an unexpected finding for a typically subclonal variant.⁸  As there can be differences in mapping efficiencies between indels and single-nucleotide variants, MYD88^L265P allele burden was then compared with the CXCR4^WHIM allele burden at the RNA level for the 11 patients with nonsense CXCR4 somatic mutations (Figure 4A). The median percentage of mutant reads was 37.4% (range, 18.4%-61.4%) and 54.5% (range, 18.2%-60.9%) for MYD88 and CXCR4, respectively (P = .0537) with the mutant allele burden of CXCR4 being greater than MYD88 in 9 of 11 patients (81.8%). Sanger sequencing of paired DNA and cDNA confirmed the mutant allele was underrepresented in the cDNA relative to the DNA for MYD88^L265P in many patients whereas the opposite was true for CXCR4 mutations (Figure 4B). This reduced expression of MYD88^L265P at the RNA level was also observed in some CXCR4^WT patients.

This study represents the first profile of the transcriptional landscape of WM by next-generation sequencing. Among the most striking findings was a >7.8 log-fold upregulation of VDJ recombination related genes including RAG1, RAG2, and DNTT in WM patients samples (Figures 5 and 6). The class switch recombination gene AICDA was not observed at meaningful levels consistent with the lack of immunoglobulin class switching in WM. The overexpression of RAG1 and RAG2 is notable because specific somatic mutation patterns such as deletions in BTG1 and ETV6 are associated with their aberrant expression, and are commonly observed in WM patients.^2,40,41  The universal expression of CXCR4 and the upregulation of its ligand CXCL12 are suggestive of autocrine/paracrine signaling. CXCR4 activation has been reported to increase cell adhesion to VCAM1 via VLA4, and VCAM1 was the eighth most upregulated gene relative vs HD samples.^42,43  Taken together, these findings suggest that VCAM1, CXCR4, and CXCL12 may facilitate the homotypic WM cell clustering observed in the BM of many WM patients. Likewise, IGF1, an inducer of AKT1 survival signaling in WM, was a top hit among the most upregulated genes. IGF1 and its receptor IGF1R may therefore be potential therapeutic targets in WM and warrant further investigation.⁴⁴ 

Despite previous studies suggesting that WM was derived from MB, we did not observe increased similarity at the transcription level for MB relative PB as shown in Figure 1B, C, and E. WM represents activated B cells in the processes of differentiating to plasma cells, evidenced by the high levels of XBP1, PRDM1, and SDC1, which may obscure signaling indicative of the true cell of origin. However, clustering by genes related to B-cell differentiation effectively sorted samples from patients with MYD88^L265PCXCR4^WT vs MYD88^L265PCXCR4^WHIM and MYD88^WTCXCR4^WT suggesting that the latter may represent an earlier stage of B-cell differentiation, consistent with the previously established lower rate of IgH somatic hypermutation in MYD88^WT WM.⁴⁵ 

The majority of transcriptional differences between MYD88^L265PCXCR4^WT and MYD88^L265PCXCR4^WHIM represented normalization of the latter for TLR4 signaling associated gene expression. This finding may help explain why CXCR4^WHIM mutations are found almost exclusively in MYD88-mutated patients. Others have also reported that CXCR4 activation can affect TLR4 signaling.^46-48  Regardless of mechanism, the direct downregulation of TLR4, and upregulation of IRAK3 in CXCR4^WHIM patients may be the most proximal explanation for this phenomenon. Although IRAK3 is a negative regulator of the IRAK4/1 signaling cascade, it has been shown that BTK activation downstream of MYD88 occurs independently of IRAK4/1 and is likely unaffected by IRAK3 upregulation.^4,49  In addition, TLR7 is upregulated in CXCR4^WHIM patients and can use the MYD88/IRAK4/IRAK3 complex to initiate regulatory signaling through MAP3K3-dependent activation of NF-κB.⁵⁰ 

Compared with the other WM genotypes and healthy donor B cells, WM patients with the MYD88^L265PCXCR4^WT genotype strongly upregulate a number of surface receptors and signaling molecules including IL17RB, GPER1, WNT5A, and IGF1. IL17RB signaling may provide an additional source of NF-κB activation whereas GPER1 and IGF1 activate AKT1 and MAPK signaling in marginal zone lymphoma and WM, respectively.^44,51,52  However, WNT5A is a suppressor of B-cell proliferation.⁵³  Many of the genes upregulated in the MYD88^L265PCXCR4^WT patient cohort such as PMAIP1, WNK2, PRDM5, CABLES1, CXXC4, and CDKN1C are tumor suppressors that either promote apoptosis, inhibit cell cycle, or block MAPK signaling.^54-58  Normalization of MYD88^L265P-associated tumor suppressor expression in CXCR4-mutated WM patients may provide the clearest explanation for the selective inhibition of MYD88 signaling observed in CXCR4^WHIM patients. This yin-and-yang relationship between MYD88 and CXCR4 is further supported by the modulation of allele specific transcription, decreasing and increasing the mutation allele burden of each gene, respectively. CXCR4 and MYD88 total transcription levels were negatively correlated and inversely predictive of BM disease involvement. Although aberrant promoter methylation of WNK2 and PRDM5 was prominent in MYD88^WTCXCR4^WT samples, it did not explain the low levels of WNK2 and PRDM5 observed in patients with CXCR4^WHIM mutations, thereby strengthening the case that these findings are related to CXCR4^WHIM signaling rather than unrelated epigenomic events.

The strongest gene markers for patients with MYD88^L265PCXCR4^WHIM was the upregulation of CXCR7, which like CXCR4 is activated by CXCL12, and TSPAN33 which is upregulated in activated B cells.^59,60  IL15 was uniquely suppressed in patients with MYD88^L265PCXCR4^WHIM and warrants functional studies to assess its role in WM microenvironment. Preclinical studies have shown that CXCR4^WHIM transduced cell lines are resistant to PIK3CD inhibitors such as idelalisib.^10,38  Both PIK3R5 and PIK3CG were significantly higher in CXCR4^WHIM patients. PIK3R5 interacts with the β/γ G-protein subunits downstream of G-protein–coupled receptors such as CXCR4 to activate PIK3CG which may explain the observed resistance.⁶¹  The G-protein and MAPK inhibitors RGS1, RGS2, RGS13, DUSP1, DUSP2, DUSP4, DUSP5, DUSP10, DUSP16, and DUSP22 as well as the MAPK-inducible tumor suppressor ERRFI1 were significantly downregulated in patients with MYD88^L265PCXCR4^WHIM vs MYD88^L265PCXCR4^WT.⁶²  The fact that DUSP2, DUSP4, DUSP5, RGS13, and ERRFI1 were all significantly induced by CXCL12 in CXCR4^WHIM vs CXCR4^WT-transduced lines supports a role for these genes as negative regulators of CXCR4 signaling and may have implications for MYD88 MAPK signaling as well.^63,64  This is similar to a number of secondary events such as MYBBP1A mutations and deletions of TRAF3, TNFAIP3, and HIVEP2, which are thought to represent loss of negative regulators for NF-κB signaling in WM.^2,11,65  The mechanism(s) by which negative regulators of CXCR4 signaling are lost in patients with CXCR4^WHIM mutations remain to be clarified warranting further investigation.

Of the 1155 genes that were differentially expressed between MYD88^L265P and MYD88^WT patients, 552 were not observed in the MYD88^L265PCXCR4^WT vs MYD88^L265PCXCR4^WHIM gene signature. MYD88^WT patient expression was quite heterogeneous indicating pathogenetic diversity in this population. Large-scale genomic exploration will invariably be needed to help clarify the genetic basis for MYD88^WT WM disease. Regardless, a clear downregulation of genes associated with NF-κB signaling including IL6, IRAK2, TNFAIP3, NFKBIZ, NFKB2, TIRAP, PIM1, and PIM2 was observed in these patients. Other genes of interest included the downregulation of CD40, and upregulation of PTBP3, CD86, and CXCR3. Notably, IGF1R, PIK3AP1, and AKT2 were all upregulated, suggesting a reliance on PIK3 signaling. The antiapoptotic gene BCL2, known to play an important role in WM, was overexpressed expressed in all patient samples while the proapoptotic BAX was underexpressed, suggesting a mechanism independent of the MYD88 and CXCR4 mutations.^65-67 

Somatic mutations in ARID1A were associated with increased CXCL13, as well as increased BM infiltration. CXCL13 has not been previously described in WM and was also a strong predictor of BM disease involvement and hemoglobin levels. CXCL13 may also play a role in mast cell recruitment.⁶⁸  Mast cells can be observed admixed with the WM B cells in patient BM histology and are thought to play an important supportive role in the tumor microenvironment of WM.⁶⁹  Deletions of chromosome 6q were associated with over 131 differentially expressed genes. Levels of previously identified target genes such as the NF-κB negative regulator HIVEP2, as well as BCLAF1, FOXO3 and ARID1B were all suppressed in the presence of 6q deletions.²  Although only 3 patients had strong familial predisposition to WM, these patients demonstrated significant dysregulation of several cancer-associated genes including RB1, STAT5B, ZNF300, MAPK9, and NFKB1.

These studies highlight the pivotal roles of MYD88 and CXCR4 signaling in WM. CXCR4 mutations appear to function primarily as dampeners of negative regulators that are upregulated in response to mutant MYD88 signaling. The upregulation of the ligand, adhesion targets, and CXCR4 itself in all WM patients provides evidence for uniform CXCR4 dysregulation in WM and supports the development of CXCR4 antagonists for WM therapy. The upregulation of PIK3 pathway members and the increased promoter methylation in MYD88^WT patients creates a strong rationale for preclinical studies of PIK3 inhibitors and demethylating agents in this population. The upregulation of BCL2 across all WM patient genotypes also supports the development of BCL2 antagonists. Finally, CXCL13 is highly expressed by WM cells, and its expression correlates with important clinical parameters. Further studies to delineate therapeutic targeting of this cytokine in WM are warranted.

The authors acknowledge the contributions of Yaoyu Wang and John Quackenbush at the Center for Cancer Computational Biology at the Dana-Farber Cancer Institute for RNA Sequencing, and Mathew Temple, Leutz Buon, Terry Haley, and the Dana-Farber Research Computing Center for their invaluable assistance. The authors acknowledge the generous support of the WM patients who provided their samples.

This work was supported by Peter S. Bing, a Translational Research grant from the Leukemia & Lymphoma Society, a grant from the International Waldenström’s Macroglobulinemia Foundation, a National Institutes of Health, National Cancer Institute Spore 5P50CA100707-12 Research Development Award, and the Orzag Family Fund for WM Research.

Contribution: Z.R.H. and S.P.T. designed the study and wrote the manuscript; Z.R.H. performed the data analysis and informatics studies; X.L. and N.T. performed the validation and methylation studies; G.Y., X.L., and J.C. transduced the cell lines and ran the stimulations; L.X., X.L., J.C., J.G.C., J.M.V., and N.T. prepared the study samples; S.P.T., J.J.C., and T.D. provided patient care and obtained consent and samples; R.J.M., P.B., J.G., and C.J.P. selected samples and provided clinical data analysis; and S.P.T., N.M.M., and K.C.A. reviewed the data and provided expert guidance.

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

Correspondence: Steven P. Treon, Bing Center for Waldenström’s Macroglobulinemia, Dana-Farber Cancer Institute, M547, 450 Brookline Ave, Boston, MA 02215; e-mail: steven_treon@dfci.harvard.edu.