Bone marrow microenvironments that contribute to patient outcomes in newly diagnosed multiple myeloma: A cohort study of patients in the Total Therapy clinical trials

Background The tumor microenvironment (TME) is increasingly appreciated as an important determinant of cancer outcome, including in multiple myeloma (MM). However, most myeloma microenvironment studies have been based on bone marrow (BM) aspirates, which often do not fully reflect the cellular content of BM tissue itself. To address this limitation in myeloma research, we systematically characterized the whole bone marrow (WBM) microenvironment during premalignant, baseline, on treatment, and post-treatment phases. Methods and findings Between 2004 and 2019, 998 BM samples were taken from 436 patients with newly diagnosed MM (NDMM) at the University of Arkansas for Medical Sciences in Little Rock, Arkansas, United States of America. These patients were 61% male and 39% female, 89% White, 8% Black, and 3% other/refused, with a mean age of 58 years. Using WBM and matched cluster of differentiation (CD)138-selected tumor gene expression to control for tumor burden, we identified a subgroup of patients with an adverse TME associated with 17 fewer months of progression-free survival (PFS) (95% confidence interval [CI] 5–29, 49–69 versus 70–82 months, χ2 p = 0.001) and 15 fewer months of overall survival (OS; 95% CI –1 to 31, 92–120 versus 113–129 months, χ2 p = 0.036). Using immunohistochemistry-validated computational tools that identify distinct cell types from bulk gene expression, we showed that the adverse outcome was correlated with elevated CD8+ T cell and reduced granulocytic cell proportions. This microenvironment develops during the progression of premalignant to malignant disease and becomes less prevalent after therapy, in which it is associated with improved outcomes. In patients with quantified International Staging System (ISS) stage and 70-gene Prognostic Risk Score (GEP-70) scores, taking the microenvironment into consideration would have identified an additional 40 out of 290 patients (14%, premutation p = 0.001) with significantly worse outcomes (PFS, 95% CI 6–36, 49–73 versus 74–90 months) who were not identified by existing clinical (ISS stage III) and tumor (GEP-70) criteria as high risk. The main limitations of this study are that it relies on computationally identified cell types and that patients were treated with thalidomide rather than current therapies. Conclusions In this study, we observe that granulocyte signatures in the MM TME contribute to a more accurate prognosis. This implies that future researchers and clinicians treating patients should quantify TME components, in particular monocytes and granulocytes, which are often ignored in microenvironment studies.

months of progression-free survival (PFS) (95% confidence interval [CI] 5-29, 49-69 versus 70-82 months, χ 2 p = 0.001) and 15 fewer months of overall survival (OS; 95% CI -1 to 31, 92-120 versus 113-129 months, χ 2 p = 0.036). Using immunohistochemistry-validated computational tools that identify distinct cell types from bulk gene expression, we showed that the adverse outcome was correlated with elevated CD8 + T cell and reduced granulocytic cell proportions. This microenvironment develops during the progression of premalignant to malignant disease and becomes less prevalent after therapy, in which it is associated with improved outcomes. In patients with quantified International Staging System (ISS) stage and 70-gene Prognostic Risk Score (GEP-70) scores, taking the microenvironment into consideration would have identified an additional 40 out of 290 patients (14%, premutation p = 0.001) with significantly worse outcomes (PFS, 95% CI 6-36, 49-73 versus 74-90 months) who were not identified by existing clinical (ISS stage III) and tumor (GEP-70) criteria as high risk. The main limitations of this study are that it relies on computationally identified cell types and that patients were treated with thalidomide rather than current therapies.

Conclusions
In this study, we observe that granulocyte signatures in the MM TME contribute to a more accurate prognosis. This implies that future researchers and clinicians treating patients should quantify TME components, in particular monocytes and granulocytes, which are often ignored in microenvironment studies.

Author summary
Why was this study done?
• The cells around a tumor, also known as the tumor microenvironment (TME), can help a tumor grow by suppressing the immune system or fight a tumor by mounting an immune response.
• Most studies of multiple myeloma (MM) have focused on the tumor itself, rather than the bone marrow (BM) TME in which the tumor is growing.
• We hypothesized that the MM TME held clues that could help us better treat patients.

What did the researchers do and find?
• We used a gene-expression-based computational technique to determine which cell types were present in patient BM.
• Patients with BM lacking a family of innate immune cells called granulocytes presented with worse outcomes compared to other patients.
• As MM progresses from a predisease to a cancerous state, the percentage of granulocytes decreases; the patients with the fewest granulocytes had more serious diseases.

Introduction
A major aim of modern cancer therapy is to develop precision immuno-oncology in which treatment is targeted to the biology of distinct disease segments [1]. Multiple myeloma (MM) is a bone marrow (BM) cancer for which modern therapy has increased 5-year patient survival from approximately 35% in 2000 to nearly 70% in 2017 [2]; however, more than 10% of patients do not respond well [3], necessitating novel treatment modalities. To date, work on disease stratification in MM has focused on the clonal plasma cells (the tumor) rather than the environment (the BM) in which the tumor proliferates. This work has led directly to the translocation cyclin D classification [4][5][6]; however, tumor-mutation-based stratification does not fully explain patient outcomes. Furthermore, these approaches have failed to identify a consistent mutational pattern driving disease progression from abnormal, but premalignant, plasma cells to cancer. MM often starts as monoclonal gammopathy of undetermined significance (MGUS), in which a clonal plasma cell has multiplied to constitute <10% of the plasma cells in the BM but the patient is otherwise asymptomatic. Each year, about 1% of patients with MGUS will progress to MM. Should that clonal plasma cell proliferate to account for >10% of plasma cells, the patient is considered to have asymptomatic or smoldering multiple myeloma (SMM). About 10% of patients with SMM will progress to MM each year [7]. Attempts to define the nature of the plasma cell mutations driving disease progression from MGUS through SMM to MM requiring therapy have largely failed and instead focused attention on changes in subclonal structure driven potentially by the tumor microenvironment (TME) [8].
The importance of microenvironmental interactions between clonal plasma cell and stromal cells mediated by autocrine and paracrine loops has been appreciated for some time [9][10][11]. However, little information is available on quantifying the global pathologic impact of the clone on the microenvironment [12,13]. Previous work has addressed the impact of clonal tumor cells on T cell subsets, natural killer (NK) cells, macrophages, dendritic cells, and highdensity neutrophils [4][5][6][8][9][10], but none of these studies took into account the contextual composition of nontumor cells within the whole bone marrow (WBM). We hypothesized that understanding the TME will improve our ability to identify patients who will have worse outcomes with current therapies and point toward potential new treatment avenues for those patients. This hypothesis was validated using the experimental design shown in Fig 1.

Patient data
Patients considered for analysis were enrolled into the Total Therapy (TT) 3-5 trial series; details of these protocols have been described in previous publications [14][15][16]. These patients The putatively pure tumor CD138 + samples were combined with the LM22 leukocyte signature matrix and 5 additional BM-specific cell types to determine the cell-type mix in the WBM samples. Illustration based on [18]. (C) The proportion of gene expression from the tumor was computationally removed from the WBM samples, which were then assigned to one of the 5 unique microenvironment clusters shown on the left side of the panel. One of these clusters (termed "low-granulocyte" Cluster 5) was associated with lower patient PFS/OS as shown on the were enrolled between February 25, 2004 and January 11, 2014 with 25-5,241 days of followup (mean 2,628 days, standard deviation 1,387 days, S1-S3 Tables). They included 266 males and 170 females aged 32.4-75.2 years at baseline (mean = 58.48, median = 59.85); 387 identified as White, 34 as Black, 5 as Hispanic, 5 as Asian, 3 as Native American, and 2 refused to identify a race. TT 3-5 generally includes treatment with thalidomide, dexamethasone, and autologous BM transplant with personalized treatment regimens informed by tumor gene expression and cytogenetics. In this study, samples were grouped by treatment stage into time points: pretreatment, postinduction, post-transplantation, postconsolidation, and postmaintenance ( Fig 1A). All TT protocols were approved by the Institutional Review Board of the University of Arkansas for Medical Sciences. All patients provided written informed consent prior to enrollment in keeping with institutional, federal, and international guidelines. All patients enrolled in TT 3-5 with at least one WBM biopsy sample taken prior to the initiation of therapy were considered for analysis. All patients provided written informed consent to approve the use of their BM biopsy samples for research in accordance with the Declaration of Helsinki of 2013. Additional patient and sample data are provided in S1-S3 Tables.

Ethics statement
Previously untreated patients with newly diagnosed multiple myeloma (NDMM) recruited to the TT clinical trials between February 2004 and August 2017 were included in this investigation after written informed consent. The sample collection protocol was approved by the University of Arkansas for Medical Sciences Institutional Review Board (protocol #2012-12). This study was performed in accordance with the 1964 Declaration of Helsinki. Ethics approval was not required for this analysis. This study has been reported according to the STROBE guideline available in S1 Text, STROBE checklist.

Methods overview
Analysis of RNA extracted from WBM biopsies potentially contains information necessary to provide a comprehensive analysis of the microenvironment [17]. To extract this information, we developed a custom algorithm able to define and quantitate the cellular subcomponents of the BM microenvironment. We applied it to WBM trephine biopsies taken from the aforementioned TT patients. These biopsies were paired with pathologist-estimated tumor fractions and cluster of differentiation (CD)138 + -selected plasma cells from those patients, allowing the signal from the malignant plasma cells to be removed and the focus put on the residual signal derived from the microenvironment.

Generation of myeloma genome signature matrix 27 (MGSM27)
To deconvolve the cellular content of the BM of plasma cell disorders, we upgraded a pre-existing 22-leukocyte signature matrix [18] by adding tumor and 4 BM cell types (S1A Fig). Deconvolution of samples purified to contain only single cell types revealed cell types that are likely to be mis-deconvolved as other cell types (S1B Fig). This spillover, combined with biological considerations, led to post hoc combination of several deconvolution estimates, resulting in the 18 cell types presented in this analysis. More details about signature matrix construction, including a usable version of MGSM27, are available in the Automated Deconvolution Augmentation of Profiles for Tissue Specific cells [19] package available on the Comprehensive R Archive Network (https://cran.r-project.org/web/packages/ADAPTS/index.html).

Testing deconvolution cell-type predictions
The digital cell quantifier algorithm [20] was modified to deconvolve cell types from bulk mRNA expression data using our custom MGSM27 signature matrix. DCQ was chosen because it outperformed other deconvolution methods in our internal trials (S1H and S1I Fig).
To ensure that the algorithm could accurately estimate malignant plasma cells and other cell types from bulk gene expression, we calibrated CD138 + cell-type estimates (normal and malignant plasma cells) against clinical estimates of myeloma plasma cell percentage in WBM (immunohistochemistry and flow cytometry) and purified CD138 + aspirates. The deconvolution method was then validated against publicly available PBMC data [21] and pathologist estimates of eosinophils (S1E Fig, S4 Table), neutrophils (S1F Fig, S5 Table), and mast cells (S6 Table) in subsets of our patient samples.

Clustering of BM based on its microenvironment
To address whether there were distinct clusters within the microenvironment of NDMM, we constructed pseudo-CD138gene-expression profiles from the WBM by removing the signature of CD138 + cells. Patient samples from all time points (N = 867) were clustered using unsupervised machine learning to identify groups of patients with similar BM microenvironments.

Statistical analyses
For supervised machine-learning survival analysis, cell-type estimates for cells with similar functions and gene-expression signatures were added together to create 19 cell types for assessment (S1B Fig). To determine which combinations of cell types and cytogenetic characteristics best separated patients with relative long and short overall survival (OS) and progression-free survival (PFS), we used conditional inference trees [22] with the Cox proportional hazards model [23], which conveniently included permutation tests to help determine generalizability.
Kaplan-Meier curves were generated using the R Survival package [24], and p-values were estimated with a χ 2 test. Estimated survival times were calculated by the survfit function and generally reported as restricted mean ± 1 standard error. Regression analysis (including lasso regression) to evaluate feature performance was performed using Cox regression in the R package GLMnet (Bioconductor) [24,25]. Decision trees were generated using the conditional inference trees algorithm, which accounts for censoring using a Cox proportional hazards model and only adds variables when they pass a p-value cutoff as calculated by a permutation test [22,26]. Stepwise Cox regression was performed by first selecting the variable with the lowest univariate p-value. At each subsequent step, variables were evaluated with an F test to determine which had the lowest p-value. If that p-value was < 0.1, then that variable was added to the model. A t test was then performed on the regression coefficients, and any variable (except age and sex) was removed if its coefficient p-value was >0.05. The procedure was continued until no more variables could be added to the model. To prevent bias in statistical analysis, results were calculated in parallel by 2 teams. One focused on lasso regression, which features predominantly in machine-learning literature; the other used stepwise regression, a technique rooted in conservative clinical and statistical analysis. Differences in these analyses (e.g., quartiles versus continuous regression) arose as a consequence of this independent verification.
Additional methods are available in S2 Text, supplementary methods.

Validation of MGSM27 cell-type predictions
The MGSM27 signature matrix correctly identified that, post-CD138 + selection, the samples were almost entirely pure CD138 + cells (N = 423; Wilcoxon p < 0.001; S1C Fig). In WBM samples in which the global RNA contribution of the CD138 + cells to the admixture was heterogeneous (S1D Fig Morphologically mature CD117 + mast cells were unexpectedly rare (much less than the previously reported 4.8% of cells in active myeloma) [27,28], and these CD117 + cells did not correlate with deconvolved mast cell estimates.
The MGSM 27 signature matrix was also applied to publicly available peripheral blood mononuclear gene-expression data (E-MTAB-3732) [21]. This test of the algorithm's ability to reconstruct all components of mixed cell-type samples revealed that the estimated cell percentages fell within the expected clinical ranges (S1G Fig).

BM microenvironment clusters
Clustering identified 5 main microenvironment clusters (Fig 2A), with patient demographics outlined in Table 1. Of these clusters, Cluster 5 was unique both in terms of immune composition and patient outcomes (Fig 3 and S2 Fig).
Cell-type estimates from the WBM revealed that Cluster 5 was also noteworthy for its distinct microenvironment. It was significantly enriched (false discovery rate [FDR] � 0.05) for estimated T cells, dendritic cells, M1 macrophages, tumor plasma cells, and other B cells ( Fig  3C and 3D) and was significantly depleted for adipocytes, osteoclasts, eosinophils, mast cells, point. Gray lines indicate transitions between time points. Portions of colored lines with no attached gray lines had no samples in the immediately preceding or following time point. � Proportions significantly different to pretreatment time point (i.e., binomial false discovery rate < 0.05). CD, cluster of differentiation; TME, tumor microenvironment.
The "high-granulocyte" Cluster 4 was the most different from Cluster 5 (Fig 3C and 3D): it had the highest number of eosinophils, mast cells, neutrophils, monocytes, M0 macrophages, and osteoclasts, and it was also enriched for ISS stage I samples. Cluster 1 was enriched for amplification of 1q, a high-risk factor [29], and showed a lack of BM resident cells, such as adipocytes, osteoblasts, osteoclasts, and inflammatory cells such as M1 macrophages and dendritic cells. However, this cluster had standard PFS/OS times and was effectively eradicated post-transplantation. Clusters 2 and 3 showed similar levels of adipocytes, innate cells, and lymphocytes, but they had inverse ratios of osteoblasts and osteoclasts. Whereas Cluster 2 had an M2 macrophage signal, Cluster 3 had an increase in the NK cell population with modest positive outcomes (estimated mean pretreatment PFS was 95 ± 9 months, 10 events, χ

Cluster 5 showed a negative impact on survival
Microenvironment Cluster 5 had significantly impaired mean PFS of 59 ± 5 months (66 events, n = 106) compared with 76 ± 3 months (124 events, n = 248) for the other clusters combined (χ 2 p = 0.001; Fig 3A); the mean OS was 106 ± 7 months (50 events, n = 106) compared with 121 ± 4 months (103 events, n = 248) for the other clusters (χ 2 P = 0.036; Fig 3B). As might be expected, Cluster 5 was significantly enriched for the level of malignant plasma cell infiltration ( Fig 3D); however, the cluster persisted even when the tumor load was reduced by therapy ( Fig 2B) and so was not driven by it. The adverse outcome of Cluster 5 was an independent prognostic feature in multivariate analysis including percentage myeloma infiltration, GEP-70, ISS stage III status, high-risk cytogenetics, age, sex, and Cluster 5 membership to  predict OS and PFS using Cox lasso regression [30]. Even in the presence of these variables, Cluster 5 was still significantly associated with adverse PFS, but not OS (Fig 3E and 3F).

Microenvironment cluster status was modified by treatment
Microenvironment clusters revealed patient dynamics in 4 settings related to treatment: postinduction (N = 245 individual patient samples available), post-transplantation (n = 83), postconsolidation (n = 51), and postmaintenance (n = 134) (Fig 2B). Although the fraction of patients in each cluster changed during treatment, no new clusters emerged (Fig 2B). The "high-granulocyte" Cluster 4 was the second most frequent subtype at diagnosis after "lowgranulocyte" Cluster 5 and became the most common TME cluster after transplantation, constituting 41%-58% of patients post-transplantation ( Fig 2B). Other smaller clusters showed different dynamics; e.g., "low-osteoblast" Cluster 1 decreased from 16% in the pretreatment patient population to 2% after transplantation ( Fig 2B). Finally, "low-granulocyte" Cluster 5 decreased from 30% pretreatment to 28% postinduction to 12% after transplantation. Fig 2C shows  When patient samples taken at relapse (within 180 days prior to clinically defined relapse or up to 30 days after relapse, n = 17) were compared with samples taken during complete remission (n = 108), all 8 of the significantly different (p � 0.05) TME components (adipocytes, eosinophils, healthy plasma cells, mast cells, monocytes, neutrophils, osteoclasts, and tumor plasma cells; S3 Fig) were significantly different in and concordant with the "low-granulocyte" Cluster 5. In particular, estimated granulocytes were significantly lower for patients who were about to relapse. Taken together, this implies a diverse array of initial MM microenvironment states that tend to respond to treatment and converge on a "high-granulocyte" Cluster 4 remission state. In contrast, the subset of patients who retain a "low-granulocyte" Cluster 5 microenvironmental state had impaired outcomes.

The cellular content of the high-risk microenvironment was associated with outcome
The granulocytes and other depleted cell types typical of Cluster 5 were strong univariate predictors of both poor OS and PFS across all 354 patients with NDMM (Fig 4A and 4B); e.g., eosinophils (OS, log-rank p = 0.02 PFS, log-rank p = 0.05), mast cells (OS, log-rank p = 0.001; PFS, log-rank p = 0.003), monocytes (OS, log-rank p = 0.006; PFS, log-rank p = 0.006), and variables in pretreatment samples. Those cell types, ratios, and clinical characteristics marked with an asterisk remain after multivariate Cox proportional hazards lasso regression, implying that they contain independent information about patient outcomes. CD, cluster of differentiation; del, deletion; FDR, false discovery rate; GEP-70, 70-gene Prognostic Risk Score; HR, hazard ratio; ISS, International Staging System; NK, natural killer; OS, overall survival; PFS, progression-free survival; t, translocation.  [35,36]. Biclusters significantly enriched for granulocyteannotated genes were also positively correlated with PFS.
Cox lasso regression determined that TME features had not been confounded by tumor burden or other known high-risk features. Regression models including CD138 + percentage, ISS stage, GEP-70 risk status, age � 65 years, high-risk cytogenetics, or sex tested whether cellular subtypes would be removed from the model. These models, for both PFS and OS, retained tumor burden, mast cells, osteoblast/osteoclast ratio, GEP-70, ISS stage III, del17p, and age � 65 years as significant predictors of patient outcomes (Fig 4A and 4B).
Stepwise Cox regression similarly demonstrated that granulocytes retained predictive power after correcting for known prognostic markers. The stepwise regression selected GEP-70 high risk, age � 65 years, low mast cells (quartile 1 versus other), and ISS stage III as the most predictive factors for PFS and OS (Fig 4C and 4D). Thus, the deconvolved mast cell counts add significant information beyond that provided by GEP-70 and ISS stage.
Immunohistochemistry-estimated granulocyte percentages were also predictive of patient outcomes. Conditional inference decision trees [22,26] built for the 290 patients with immunohistochemistry data (S6A- S6D Fig) revealed that low-granulocyte percentages were associated with adverse outcomes; e.g., low eosinophil counts as determined by immunohistochemistry were a significant predictor of outcome (permutation p < 0.001 for PFS; permutation p = 0.045 for OS) with a mean PFS of 67 ± 4 months (115 events, N = 195) and OS of 112 ± 5 months (90 events, N = 195) compared with a mean PFS of 97 ± 10 months (13 events, N = 31) and OS of 131 ± 6 months (34 events, N = 95) for patients with high counts.
Granulocytes were important for stratifying the 286 pretreatment samples with fully quantified cytogenetics, GEP-70 high-risk status, ISS stage, and tumor burden. Conditional inference trees showed that the pretreatment samples can be best divided by GEP-70 (estimated mean PFS of 42 ± 6 months, 27 events, n = 34), ISS stage III (estimated mean PFS of 54 ± 6 months, 40 events, n = 68), and low mast cells (estimated mean PFS of 61 ± 6 months, 26 events, n = 40) compared with an estimated mean PFS of 82 ± 4 months (60 events, n = 148) for the remainder (Fig 5A). Further, patients with tumor-specific high-risk status, as defined by ISS stage or GEP-70 status, had even worse outcomes if they had low mast cell estimates (Fig 5B).

Tumor gene-expression profiles were associated with depleted granulocytes
To address what tumor-side molecular features may be driving the identification of low mast cells and thus perhaps the poor outcomes associated with depleted granulocytes, we compared the pretransplantation gene expression of the 328 patients for whom we had paired purified CD138 + -selected samples and WBM samples. The WBM samples were used to determine the estimated percentage of mast cells and were divided into 2 groups: those with the lowest quartile of mast cells in the WBM and the remainder with higher mast cells. The top 32 genes in the CD138 + RNA-expression data that were differentially expressed between these 2 groups, as measured by a Wilcoxon signed-rank test (FDR � 0.05) and with an absolute log 2 fold change > 0.66, are shown in S7 Fig.

The development of a high-risk microenvironment during MM progression
To understand how the patterns seen in NDMM develop during the multistep progression from MGUS through SMM to NDMM, we applied the same deconvolution techniques to investigate WBM samples from 131 premalignant patients (MGUS, n = 55; SMM, n = 76). This revealed that patients with more advanced SMM have significantly increased tumor plasma cells (Wilcoxon signed-rank p = 0.004) and plasma cells (p = 0.01) and significantly decreased monocytes (p = 0.01) and eosinophils (p = 0.018), consistent with the patterns observed as patients moved toward worse outcomes in the MM cohort. Only 6 of the patients with MGUS/SMM were known to progress, so no significantly different cell types were observed. Using our previously defined 5 microenvironment clusters, 97% (127/131) of these patients (54/55 of MGUS and 73/76 of SMM) were assigned to "high-granulocyte" Cluster 4, the dominant post-transplantation remission microenvironment cluster. , and no other attribute can significantly split the patient population. The algorithm then selects the attribute that best divides the remaining patient population, which is ISS stage III status (p < 0.001). These 68 patients (40 of whom had an event) without GEP-70 high-risk tumors but who are ISS stage III cannot be significantly divided by any other attribute and have relatively poor outcomes (estimated mean PFS = 54 months). The remaining 188 patients (86 of whom had an event) are best divided by deconvolution-based mast cell estimates (p = 0.001). Those 40 patients (26 of whom had an event) with predicted mast cell levels of 8.97% or lower also have relatively short PFS (estimated mean PFS = 61 months). The remaining 148 patients had a relatively long PFS of 82 months. Importantly, the model does not select any more attributes that can pass a cutoff of p = 0.01, implying that of these deconvolved immune cell estimates and clinical characteristics, GEP-70, ISS stage III, and mast cell estimates are the strongest predictors. Clinically estimated tumor burden, which is inversely correlated with deconvolved mast cell percentages, was not selected by the tree and is, therefore, less informative. Unsupervised clustering of those premalignant patients, using the top 10% of genes ranked by variance, resulted in 3 premalignant microenvironment clusters (pClusters). pCluster 1 comprised 38 patients and 29% (11/38) had SMM; pCluster 2 comprised 64 patients and 63% (40/64) had SMM; pCluster 3 comprised 29 patients and 86% (25/29) had SMM. As these pClusters became enriched for SMM, plasma cell percentages increased while the estimated population of granulocytes and monocytes decreased (Fig 6A). As the disease progressed, the TME patterns became more like those of the MM patients with the worst outcomes, i.e., "lowgranulocyte" Cluster 5. Notably, patients with NDMM outside Cluster 5 had profiles similar to patients in pCluster 3, including similar levels of estimated plasma cells, implying that the Cluster 5 microenvironment is the final step in disease progression (Fig 6B).
The mast cell signature remained relevant in post-treatment patients with relapsed/refractory MM (RRMM), so we explored the same patterns in 295 patients with RRMM treated with lenalidomide in the STRATUS (Evaluation of Safety of Pomalidomide in Combination With Dexamethasone (Low Dose) in Patients With Refractory or Relapsed and Refractory Multiple Myeloma) trial [37]. As may be expected for patients with poor outcomes, the patients with RRMM were enriched for a depleted mast cell signature; 70% (207/295) had low mast cell estimates compared with 27% (87/320) of baseline pretreatment patients, as shown in Fig 5B  (Fisher exact p < 0.001).

Discussion
Here, we present an analysis of 867 WBM samples taken during the course of treatment of 436 patients with MM (Fig 1A and 1B) and a follow-up study using samples from predisease patients (55 with MGUS and 76 with SMM, S1-S3 Tables). Purified tumor samples were used to isolate the nontumor portion of WBM gene-expression data, which clustered patients by their microenvironment signatures (Fig 1C). One of the clusters (Cluster 5) stood out both in terms of patient survival and cellular composition. The "low-granulocyte" Cluster 5 microenvironment signature was shown to add significant additional prognostic information to standard tumor-specific risk-stratification tools, including the ISS stage; GEP-70 high-risk status; cytogenetic features including amp1q, del17p, adverse translocation t(4;14), and t(14;16); and "double-hit" disease [38] (Fig 1D). Follow-up studies on predisease SMM and MGUS patients revealed that this microenvironment developed as the disease progressed ( Fig 1E). Machinelearning techniques identified WBM components that independently impacted survival beyond known high-risk tumor features, with mast cell percentage being the single most predictive factor (Fig 1F). This signature was validated in the STRATUS trial, in which patients with RRMM were significantly enriched for the low-granulocyte phenotype (Fig 1G). The results of this analysis show that the nontumor composition of the MM microenvironment independently impacts survival and can be used to refine risk-stratification approaches.

What this study adds to existing research
Existing research in MM has tended to focus on tumor mutations and patient clinical characteristics. This is the first large-scale study of the MM TME (i.e., the WBM). It reveals a lowgranulocyte microenvironment that is associated with negative patient outcomes and shows how that microenvironment develops during predisease states. This study demonstrates that a WBM granulocyte (i.e., mast cell) signature may be used to stratify patients and identify an elevated risk population that would be missed by existing high-risk criteria.
It is important that this study implicates mast cells and other innate immune cells because most TME research focuses on the adaptive immune system (e.g., T cells). This may be because the deconvolution techniques do a relatively bad job of detecting T cell populations based on The bottom section shows a plausible mechanism of interaction between the tumor and microenvironment involving VCAM1 and other factors shown our internal validation data set, indicating room for further improvements in deconvolution technology (S1I Fig). The WBM samples in this study contain many fewer mature mast cells than the 4.8% reported in the myeloma literature [27,287]. Although this may be due to differences in staining methods, it emphasizes the need for additional studies of innate immune cells in myeloma patients.

Strengths and limitations of this study
This study's great strength is that it has up to 15 years of follow-up data for a high number of MM patients (436) sampled at different time points during treatment (867 samples). The patterns observed in this cohort have been validated in 131 high-risk predisease samples and 295 post-treatment relapse samples in the STRATUS cohort.
This study is limited in that the patients were primarily treated with thalidomide rather than the current front-line drug, lenalidomide (although lenalidomide was used in the STRA-TUS cohort). This study is also limited by the computational deconvolution used to enumerate cell types. Deconvolution demonstrates poor ability to differentiate cell types in the same lineage (S1D Fig). This challenges the identification of mast cells as the cell type best able to predict patient outcomes. Highly correlated mast cell, eosinophil, neutrophil, and monocyte counts make it difficult to tell whether mast cell estimates reflect another cell type or even a common progenitor. Implications of this limitation are presented in S3 Text.

Implications and next steps for research
The importance of nonmacrophage myeloid lineage cells in this study implies several plausible partial mechanisms that should be tested experimentally in the context of MM. Mast cells have been shown to play a role in the initiation of antitumor immune response [39], and differential gene-expression analysis supports this observation. Comparison of tumor gene-expression profiles from patients with low mast cell signatures with the other patients yields only 2 genes with significantly lower expression; both are associated with mast cell recruitment and activation (Fig 6B and S7 Fig). Vascular cell adhesion protein 1 is an adhesion protein, utilized by circulating premast cells in order to move into a tissue [40]; thus, tumors lacking this protein likely attract fewer mast cells. Proteoglycan 2 is a mast cell activation protein, the absence of which renders recruited mast cells less likely to be activated and thereby recruited by other immune cells [41].
https://doi.org/10.1371/journal.pmed.1003323.g006 activate mast cells and neutrophils [41]. Thus, it might be profitable for future researchers to purify these cell types (both individually and combined) from MM patient samples. This could lead to a more precise biomarker than the mast cell signature.
Beyond biomarker discovery, in vitro investigations of stromal cells highlighted in this study may be particularly profitable. By illuminating the interactions between these cell types and the tumor, new findings may suggest treatments for myeloma that complement existing therapies. One outstanding question is why patients who enter microenvironment Cluster 5 and leave do not have better outcomes than those who stay in Cluster 5. Has the microenvironment suffered some damage that is undetectable by the assays used here? One way to answer this sort of question would involve reconstructed BM 3D bioassays [45][46][47]. Models of these kind have shown improved ability to model tumor cell survival in drug-resistant MM. Therefore, they might provide a controlled environment to measure signaling as tumor and microenvironment components interact to determine tumor proliferation and drug resistance.

Implications and next steps for clinical practice
Mast cell estimates derived from WBM gene expression are a potentially useful biomarker to determine which patients have elevated risk status. As shown in Fig 5 (and S8 Fig), patients with low mast cell estimates (<8.97%) have worse outcomes than patients with higher mast cell estimates, whether or not the patient has been identified as having elevated risk by ISS stage or GEP-70 criteria. Patients who relapse show signs of the same "low-granulocyte" Cluster 5 that develops as MGUS and progresses to MM. If the microenvironment change precedes the malignant cell changes, then it is possibility that the "low-granulocyte" microenvironment may indicate patients who are about to relapse. Taken together, this implies that clinicians should be mindful of the MM WBM tumor microenvironment. WBM gene-expression data can flag elevated risk patients. Enumerations of the cell types present in the WBM of patients with MM, particularly of myeloid lineage cells such as mast cells, eosinophils, neutrophils, and monocytes, could test the conclusions of this study and potentially reveal new ways to treat patients who would otherwise suffer poor outcomes.

Conclusion
Taken together, the observations in this study provide evidence that the microenvironment is a significant contributor to tumor behavior and treatment response. Further investigation of the microenvironment cellular content may have significant impacts on outcome after immune therapies and, in addition, the role of atypical myeloid cells including granulocytes, myeloid derived suppressor cells, myeloid precursors, and other myeloid lineage cell types should be explored further in this context. To summarize, this article should encourage myeloma clinicians and researchers to study the TME and pay particular attention to granulocytes; their presence (or absence) may contribute to patient survival.   Also shown are the cell types that were combined to make the 18 cell types used throughout the study. Note that cell types with very small distances were collapsed unless they were biologically very different cell types (e.g., macrophages and dendritic cells). (C) Shown is the comparison between estimated deconvolution and pathologist-determined tumor purity for 423 patients with gene-expression data from purified CD138 + samples and flow-based tumor percentages. The left side shows a scatter plot comparing deconvolution and flow-based purity of CD138 + samples. Boxes drawn at 75%, 80%, and 90% show that both the deconvolved estimates and clinical estimates indicate high levels of CD138 + purity. The right side shows the same data if samples that are more than 90% CD138 + are considered "pure." (D) Comparisons between estimated deconvolution and pathologist-determined tumor percentages for 247 patients with matched pretreatment WBM and CD138 + gene-expression data, as well as microscopy and flow-based tumor percentages. (E) Shown is the deconvolved eosinophil estimates and the pathology-estimated percentage for proximal samples taken from 345 patients. Note the pathology estimates show vertical bands consistent with rounding to the nearest half a percent, impeding correlation analysis. A blue X denotes the mean deconvolution percentage of those samples with a pathology estimate of less than 6%, and the dotted blue line shows the regression line for those samples. A green X denotes the mean deconvolution percentage of the remaining samples, and the dotted green line shows their regression line. Importantly, the predictions are in the correct range: 323/345 (94%) of the samples are less than 8% eosinophils as predicted by the deconvolution. (F) Deconvolved neutrophil estimates and the pathologyestimated percentage for proximal samples taken from 356 patients. (G) Deconvolved estimates of PBMC data extracted from E-MTAB-3732: 182 samples. The percentages on the left side show the expected cell types from pathologist-purified peripheral blood mononuclear cells, and those on the right show the deconvolved estimates. The ranges on the right show the mean ± 1 standard deviation. The deconvolved estimates for cell types are calculated by adding the means of all cells of a particular type, as indicated by the color code; e.g., the pink B cell estimates (13%) are the sum of 10% MM.plasma.cells, 2% PlasmaMemory cells, and 1% B.cells. memory. (H) Shows the RMSEs for different deconvolution algorithms to predict the percentage tumor in CD138 + purified aspirates and WBM. Some of these samples were replicates that were not included in the main study. As an orthogonal approach to validate the biological relevance of low-granulocyte populations in the generation of a high-risk microenvironment, we constructed patient-specific clusters of coexpressed genes (i.e., biclusters) [31] using cMonkey2 separately on WBM and CD138 + -selected samples. cMonkey2 generates gene clusters that are often highly enriched in genes of related function and/or in related biochemical or regulatory pathways [32][33][34]. We expected biclusters from this analysis to be enriched for genes driven by specific processes in cell types that are transcriptionally related. By contrasting those biclusters between the 2 data sets, we expected to identify processes that were specific to the high-risk microenvironment. (A) cMonkey2 biclusters significantly correlated with patient PFS. Shown are GO terms enriched in cMonkey2 biclusters (out of a maximum 9,258 clusters for CD138 + data and 8,540 clusters for WBM data). The y-axis shows the number of clusters built on the pretreatment WBM samples; the x-axis shows the number of clusters built on CD138 + samples. This analysis resulted in 8,540 (potentially redundant) WBM biclusters and 9,258 CD138 + biclusters, each containing an average of 30.0 and 30.3 genes, respectively, from 187.0 and 253.6 patients, respectively. After performing functional enrichment analysis using the GO terms [35,36], we found that although many common functions (e.g., cell cycle) were enriched in similar numbers of biclusters between the 2 data sets, many more biclusters generated from the WBM data (approximately 130 biclusters, or approximately 1.5% of all biclusters) were enriched for granulocyte-and myeloid-related genes (specifically, activation or degranulation of neutrophils, leukocytes, and/or myeloid cells) than were biclusters generated from the CD138 + data (approximately 15 biclusters, or approximately 0.15%). (B) shows the top GO terms enriched in WBM clusters. The greater number of biclusters enriched for granulocyte-related functions in WBM is a clear indication of the greater expression from cells in the TME in the WBM data set, whereas the small but nonzero number of biclusters from the CD138 + data set could be readily explained by the expected 0%-10% CD138contamination in CD138 + samples. In both data sets, >90% of these granulocyte-associated clusters were positively correlated with PFS (i.e., greater expression of these genes was associated with longer PFS, on average; average correlation of 0.092) compared with 48% for the remaining WBM biclusters (average correlation of -0.01; t test difference p < 10 −16 ). Taken together, this biclustering and GO functional discriminatory analysis provides a separate analysis pipeline, which narrows down the relevant cell types present in Cluster 5 to implicate granulocytes as an important correlate of patient outcomes. ATP, adenosine triphosphate; CD, cluster of differentiation; DNA, deoxyribonucleic acid; BP, biological process; ER, endoplasmic reticulum; GO, gene ontology; ncRNA, noncoding ribonucleic acid; PFS, progression-free survival; RNA, ribonucleic acid; SRP, signal-recognition particle; TME, tumor microenvironment; WBM, whole bone marrow.