A decreased granulocyte proportion predicts a poor outcome in NDMM

To better understand the characteristics of the BM microenvironment during MM progression, we applied MGSM27 signature matrix, a gene-expression–based computational technique [27], to mRNA gene expression profiles (GEP) of BM biopsy from patients with MGUS (n = 122), SMM (n = 122) and NDMM (n = 704), as well as healthy donors (normal BM (NBM), n = 67) (Supplementary Tables 1). This so-called CIBERSORT [26, 28] algorithm reconstructed 27 BM cell types and revealed the estimated cell proportions from the BM biopsy. We removed plasma cells (PC) from the 100% total cell proportions to avoid bias resulting from different tumor burdens. CIBERSORT revealed six types of innate immune cells, neutrophils, eosinophils, mast cells, dendritic cells (DCs), macrophages, and monocytes. It revealed that the proportions of neutrophils, mast cells, and monocytes account for the majority of innate immune cells and were decreased in NDMM and its precursor stages when compared with NBM (Fig. 1a). We noticed a 45% decrease of M2 macrophages in NDMM compared with NDMM, but not in the precursor stages (Fig. 1a). We performed correlation analysis of immune cell proportions with the time of progression of MGUS and SMM patients. We observed high neutrophil proportions correlated with short progression time in SMM patients (Fig. 1b). However, NDMM patients with high proportions of neutrophils had superior overall survival (OS) and event-free survival (EFS) (Fig. 1c). In support of that, the neutrophil percentage was decreased in patients with high-risk status based on the 70-gene Prognostic Risk Score (GEP-70) or at stage III of International Staging System (ISS) (Fig. 1d, e). We observed similar findings related to mast cells (Fig. 1f, g), indicating that a decreasing proportion of granulocytes predicts inferior survival in NDMM patients. In addition, MM patients with high-risk GEP-70 had a higher M2 macrophage proportion (Supplementary Fig. 1a). M2 macrophages were increased in MM patients with ISS stage III (Supplementary Fig. 1b). Collectively, these data suggest that the alterations of innate immune cells proportion start early in MGUS stage and are predictive of outcome. Notably, we observed MM patients with high or low MM plasma cell percentages displaying distinct immune signatures and that the immune signatures in MM patients with low MM plasma cell percentages resembled that of MGUS or SMM.

Fig. 1: A decreased neutrophil proportion predicts a poor outcome.

a Bar plot showing proportions of the 10 innate immune cell types of NBM (normal bone marrow, n = 67), MGUS (monoclonal gammopathy of unknown significance, n = 122), SMM (smoldering multiple myeloma, n = 122), and NDMM (newly diagnosed multiple, n = 704) from CIBERSORT analysis. Error bars represent mean ± standard error mean (SEM). p values for each cell type were calculated using Kruskal–Wallis with Dunn’s multiple comparisons test. b Correlation analysis of neutrophil proportions with the time of progression in SMM patients. Pearson r = -0.572, p < 0.0001. c Kaplan–Meier analyses of overall survival (OS, up) and event-free survival (EFS, down) in MM patients with high and low neutrophil proportions. Hazard ratios OS = 22.69, p < 0.0001; Hazard ratios EFS = 2.59, p < 0.0001. High vs. Low groups were determined by the optimal cut point of EFS survival. d Violin plot showing the proportions of neutrophils in MM patients with GEP70-Low (n = 598) and GEP70-High (n = 105) (GEP70 cut off 0.66). p values were calculated using Mann–Whitney test. e Violin plot showing the proportions of neutrophils in MM patients in R-ISS stages I (n = 88), II (n = 193), and III (n = 42) (right). p values were calculated using the One-way ANOVA with Tukey’s multiple comparisons test. f Kaplan–Meier analyses of overall survival (OS, left) and event-free survival (EFS, right) in MM patients with high and low mast cell proportions. Hazard ratios OS = 1.35, p = 0.0018; Hazard ratios EFS = 1.28, p = 0.0068. High vs. Low groups were determined by the optimal cut point of EFS survival. g Violin plot showing the proportions of mast cells in MM patients with GEP70-Low (n = 598) and GEP70-High (n = 105) (GEP70 cut off 0.66). p values were calculated using Mann–Whitney test.

Reduced HLA-DR surface expression in CD16+ monocytes and pDCs is detected early at the MGUS stage

To evaluate the immune cell ratio and function based on protein expression, we performed deep phenotyping of BM aspirates of patients with MGUS (n = 28), SMM (n = 30), NDMM (n = 33), as well as 20 NBM using CyTOF (Supplementary Tables 2, 3). Clustering based on cell lineage marker expression using the FlowSOM algorithm [29] identified 17 metaclusters (MCs, Fig. 2a, Supplementary Fig. 2a). We observed an enrichment of megakaryocytes, CD16+ monocytes, CD8+ Temra (terminally differentiated effector), CD4+ Tcm (central memory), and CD56+ MM in the TME, as well as a decrease of plasmacytoid DCs (pDCs), immature granulocytes, and CD27+ normal plasma cells (NPC) (Fig. 2b). In addition to these immune cell alterations, we also observed an enrichment of NK cells, CD4+ Tem (effector memory) and CD4+ Temra in the TME after excluding the two malignant plasma cell (MPC) clusters (Supplementary Fig. 3). CyTOF showed 65% decrease of an immature granulocytes in the MGUS stage compared to NBM, 53% decrease of pDCs in SMM and 51% decrease in NDMM compared to NBM, and 5.6-fold increase of CD16+ monocytes in MGUS, 2.0 -fold increase in SMM and 1.6-fold increase in NDMM compared with NBM (Fig. 2c). Next, we investigated the functional marker expression among innate immune cells. Complement inhibitory proteins CD55 and CD59 showed cell type specificity in monocytes and immature granulocytes, respectively, but only CD55 showed increased intensity in MM compared with its precursor stages (Fig. 2d). A decrease of HLA-DR expression was observed in pDCs and CD16+ monocytes in both MM and its precursor stages (Fig. 2d), indicating that a decrease of myeloid cell antigen presentation function begins early in the MGUS stage.

Fig. 2: Reduced HLA-DR surface expression in CD16+ monocytes and pDCs is detected early at the MGUS stage.

a viSNE plot showing 17 metaclusters (MC) in the normal bone marrow (NBM) and tumor microenvironment (TME) were identified by FlowSOM analysis. b Volcano plot showing MC changes between NBM and TME. For each cell type the log fold change in mean cell fraction between tumor and normal samples, with −log10 Benjamin–Hochberg-corrected, two-sided Wilcoxon rank-sum p values on the y-axis is shown. c Bar plot showing innate cell composition changes among different disease stages. p values for each cell type were calculated using Kruskal–Wallis with Dunn’s multiple comparisons test. d Heatmap showing median intensity of indicated channels in each innate immune cell types. p values for each cell type were calculated using Kruskal–Wallis with Dunn’s multiple comparisons test.

Decreased numbers of mDCs and their expression of stress and immune response genes are observed in MM and its precursor stages

To reduce the interference of individual heterogeneity and to further explore gene expression associated with immune cell functional alterations during disease progression, we collected sequential BM aspirate samples that consisted of two pairs of stable MGUS samples, one MGUS progressed to SMM, one MGUS progressed to MM, and 3 pairs of SMM progressed to MM, which include a total of 6 MGUS, 4 SMM, and 4 MM, along with 5 NBM controls (Supplementary Tables 4), and performed single-cell RNA sequencing (scRNA-Seq; 10x Genomics). By clustering cells based on gene expression, we identified 12 cell clusters (Fig. 3a, Supplementary Tables 5). Consistent with previous data [20], we observed a decrease in pDCs, myeloid DCs (mDCs), immature granulocytes, and CD14+ monocytes in MM and its precursor stages, as well as an increase in T cells, CD16+ monocytes, megakaryocytes both before and after excluding PCs (Fig. 3b, Supplementary Fig. 4). MM and NBM showed a distinct pattern of immune cell cluster distribution among total immune cells, with an overall decrease of the proportion of innate immune cell clusters and an increase of NK and T cells (Fig. 3c). The stable MGUS pairs resembled NBM in terms of the proportion of CD14+ monocytes, immature granulocytes and NK cells when compared to the progressed MGUS and SMM samples (Fig. 3c).

Fig. 3: Decreased numbers of mDCs and their expression of stress and immune response genes is observed in MM and its precursor stages.

a UMAP plot showing 12 cell clusters identified in the BM aspirate samples. b Volcano plot showing cell cluster changes between normal bone marrow (NBM) and tumor microenvironment (TME). For each cell type the log fold change in mean cell fraction between tumor and normal samples, with −log10 Benjamin-Hochberg-corrected, two-sided Wilcoxon rank-sum p values on the y-axis is shown. c Heatmap showing the fractions of different cell clusters in non-tumor BM cells for individual patients and healthy donors. Cell cluster fractions are z-standardized across patients. d Dot plot showing gene expression in myeloid dendritic cells for individual patients and healthy donors. e UMAP plots of mDC subtypes (left) and stacked box plots showing their proportions (right). p values for each cell subtype were calculated using Benjamin–Hochberg corrected, two-sided Wilcoxon rank sum test.

As myeloid DCs showed a significant decrease in cell number and proportion between NBM and tumor microenvironment (TME), we further explored the alterations of their gene expression among patients and healthy donors. We found the upregulated genes were highly enriched in IFN-responsive pathway (OAS2, IFITM3, ISG15, IFI44) and showed increased expression in all MM (4/4) and 2/4 SMM samples, whereas all MGUS samples (6/6) remain unchanged. The downregulated genes were highly enriched for protein folding and cellular stress response (HSPA1B, HSPA1A, HSPA6, DNAJB1), inflammation (TNFAIP3, IL1R2, FCGR2B), antigen presentation (MRC1, HLA-DQA2), and immune cell recruitment and activation (CCL3L1, CXCL8, CCL3) in TME (Fig. 3d, Supplementary Fig. 5a). These findings suggest that mDCs in the TME have a decreased ability to respond to cellular stress and to execute their immune function. Subclustering of mDCs based on their gene expression revealed 6 mDC subsets, which were annotated by their marker genes (ADAM8+, S100A8+, GINS2+, CCNA2+, VSIG4+, and XCR1+, Fig. 3e, Supplementary Fig. 5b). The proportion of VSIG4+ mDC subset showed a trend of decrease in the TME (Fig. 3e). VSIG4, a B7 family related protein, is a negative regulator of T cell activation [30]. Whether mDCs in MM and its precursor stages exhibit impaired function in regulating T cell activation needs further investigation.

NK cells are overactivated and exhausted in TME, and are characterized by high levels of CD57, GZMB, and TIM3

Next, we examined NK cell proportion and function during MM progression. CyTOF analysis showed an increase of immunomodulatory (CD56+CD16-) NK cells in NDMM compared to other stages of disease, while cytolytic (CD56dimCD16+) NK cells were increased in both MGUS and NDMM compared to NBM (Fig. 4a). Both cytotoxic markers (CD57 and CD161) and immune inhibitory checkpoints (TIM3 and TIGIT) had increased intensity among the cytolytic NK cells. CD57, also a marker of NK cell terminal differentiation, anergy and senescence [31], was increased in NDMM compared with other groups. TIM3 was increased in MGUS, SMM, and MM compared with NBM, and TIGIT was elevated in SMM and NDMM compared with NBM and MGUS (Fig. 4b).

Fig. 4: NK cells are overactivated and exhausted in the TME, characterized by high levels of CD57, GZMB and TIM3.

a Box plot showing NK cell composition changes among different disease stages. p values for each cell type were calculated using Kruskal–Wallis with Dunn’s multiple comparisons test. b Heatmap showing median intensity of indicated channels in each NK cell types. p values for each cell type were calculated using Kruskal–Wallis with Dunn’s multiple comparisons test. c Volcano plot showing differential gene expression in immature granulocytes between normal bone marrow (NBM) and tumor microenvironment (TME). d Dot plot showing gene expression of NK cells for individual patients and healthy donors. e UMAP plots of NK cell subtypes (left) and stacked box plots showing their proportions (right). p values for each cell subtype were calculated using Benjamin–Hochberg corrected, two-sided Wilcoxon rank sum test. f Density map on the UMAP plot showing gene expression levels among NK subtypes.

Next, we investigated NK cell clusters using scRNA-Seq analysis. Similar to what we found in the mDCs, gene expression of NK cells showed an upregulation of IFN response genes and a downregulation of stress response genes in the TME compared with NBM (Fig. 4c, d). Interestingly, with regards to the immune response genes, we observed a decrease of LYZ, GZMK, XCL1, CCL3, CXCR4, and an increase of CX3CR1, GZMH, GZMB, HAVCR2 (encoding TIM3) in the TME (Fig. 4c, d), indicating a switch of immune response features of NK cells during MM progression. Subclustering of NK cells revealed three populations, CX3CR1+GZMB+, XCL1+GZMK+, and a less frequent IL1R1+KIT+ subset (Fig. 4e, f). The IL1R1+KIT+ NK cell subset expressed high levels of TCF7, SELL, and IL7R, indicating those cells were naïve. The CX3CR1+GZMB+ NK cell subset was also characterized by cell surface marker FCGR3A+ (encoding CD16a), high levels of cytolytic genes (GZMH, PRF1), and the exhaustion marker HAVCR2 (encoding TIM3) (Fig. 4e, f). Notably, we observed a proportion shift from the IL1R1+KIT+ and XCL1+GZMK+ NK cell subsets towards the CX3CR1+GZMB+ NK subset in the TME compared with NBM (Fig. 4e, f). Collectively, CyTOF and scRNA-Seq data suggested that despite an increase in the proportion of cytolytic NK cells in the TME, they may be overactivated and thus express exhaustion markers.

Accumulation of PD1+CD4+ and TIGIT+CD8+ inhibitory T cells during MM progression

CIBERSORT analysis of the BM biopsy GEP data revealed a 53% increase in CD8+ T cells, a 63% decrease in naïve CD4+ T cells, and a 42% decrease in activated CD4+ memory T cells in NDMM compared to NBM (Fig. 5a). We observed that high proportions of CD8+ T cells while low proportions of γδ T cells correlated with short progression time in SMM patients (Fig. 5b). We observed low proportions of γδ T cells were also correlated with short progression time in MGUS patients (Supplementary Fig. 6). Consistently, high-risk SMM patients showed higher proportions of CD8+ T cells and lower proportions of naïve CD4+ T cells and γδ T cells (Fig. 5c). However, we did not find that T cell subpopulations were associated with OS or EFS in NDMM (Supplementary Fig. 7).

Fig. 5: Accumulation of PD1+CD4 and TIGIT+CD8 inhibitory T cells during MM progression.

a Bar plot showing the proportions of T cell types of NBM (normal bone marrow, n = 67), MGUS (monoclonal gammopathy of unknown significance, n = 122), SMM (smoldering multiple myeloma, n = 122), and NDMM (newly diagnosed multiple, n = 704) from CIBERSORT analysis. Error bars represent mean ± standard error mean (SEM). p values for each cell type were calculated using Kruskal-Wallis with Dunn’s multiple comparisons test. b Correlation analysis of CD8 T cell (left, Pearson r = −0.308, p = 0.001.) and γδT cells (right, Pearson r = 0.219, p = 0.021) proportions with the time of progression in MGUS patients. c Violin plot showing the proportions of the CD8 T cells (left), γδT cells (middle), and naïve CD4 T cells (right) among high-risk (n = 35) and low-risk (n = 77) SMM patients from CIBERSORT analysis. Error bars represent mean ± standard error mean (SEM). p values were calculated using Mann–Whitney test. d Box plot showing T cell composition changes among different disease stages by CyTOF analysis. p values for each cell type were calculated using Kruskal-Wallis with Dunn’s multiple comparisons test. e Heatmap showing median intensity of indicated channels in each T cell types. p values for each cell type were calculated using Kruskal-Wallis with Dunn’s multiple comparisons test.

To investigate the phenotype and function of T cell populations, we analyzed the 5T cell MCs identified by FlowSOM analysis of CyTOF data, including a population of CCR7+CD45RA- central memory T cells (Tcm), two CCR7-CD45RA- effector memory T cells (Tem), and two CD57+ terminal differentiation T cells (Temra). We observed an increase of CD8+ Tem, CD8+ Temra, CD4+ Tem, and CD4+ Tcm proportion at the precursor stages (Fig. 5d). Functional marker analysis identified an increase in CD25 and FOXP3 intensity among CD4+ Tcm in NDMM and its precursor stages (Fig. 5e), suggesting CD4+ Tcm exhibit immune inhibition in the TME. We evaluated immune checkpoint (PD1, TIGIT, TIM3, and LAG3) markers during different disease stages. PD1 was elevated in CD4+ Temra, while TIGIT was elevated in CD8+ Temra and CD8+ Tem. Both PD1 and TIGIT showed increased levels during disease progression (Fig. 5e). These data revealed an accumulation of inhibitory T cells in the TME.

GZMB + cytotoxic T cell subsets are enriched and show individual specificity in the TME

Consistent with CIBERSORT and CyTOF analyses, scRNA-Seq also revealed an overall increase of total T cells in the TME (Fig. 3c, d). To delve deeper into the T cell subtypes and function, we performed additional subclustering at higher resolution using genes variably expressed among T cells, resulting in 15 subclusters (Fig. 6a). The non-cytotoxic subsets consist of four populations of CCR7+TCF7+ naïve CD4 subsets, one population of CCR7+TCF7+ naïve CD8 subset, one population of CD4 helper subset, and one IFN-responding CD4 subset. The cytotoxic subsets contain four populations of GZMB+ CD8 subsets, two GZMK+ CD8 subsets, one GZMB+ CD4 subset and one GZMB+ γδ T subset (Fig. 6a, b). The non-cytotoxic subsets were characterized by high expression of immune modulatory genes (LTB, IL7R). GZMB+ cytotoxic subsets expressed a great number of cytotoxic genes, such as GZMB, GZMH, GNLY, NKG7, KLRD1, while GZMK+ cytotoxic subsets expressed GZMK (Fig. 6c). Consistent with other reports, the proportion of the CD8 GZMK+ subset showed a trend of decrease while the CD8 GZMB+ subset showed a trend of increase in the TME (Fig. 6d). Interestingly, we observed that two CD8 GZMB+ subsets, one CD4 GZMB+ subset and one γδ T GZMB+ subset was predominantly derived from individual patients and were detected at both time points (Fig. 6d), suggesting that GZMB+ cytotoxic T cells may perform anti-tumor specific immune response in the TME.

Fig. 6: GZMB+ cytotoxic T cell subsets are enriched and show individual specificity in the TME.

a UMAP plots of T cell subtypes. b Density map on the UMAP plot showing gene expression levels among T subtypes. c Dot plot showing marker gene expression of T cell subtypes. d stacked box plots showing the proportions of noncytotoxic T cell subtypes (left) and cytotoxic T cell subtypes (right). p values for each cell subtype were calculated using Benjamin–Hochberg corrected, two-sided Wilcoxon rank sum test.

We next examined the differential gene expression between TME and NBM. Similar to the findings with innate immune cells, we observed an upregulation of IFN responding genes and a downregulation of stress response genes in TME (Supplementary Fig. 8a, b). In addition, the expression levels of cytotoxic genes LYZ, GZMK were decreased while PRF1, GZMH were increased in the TME (Supplementary Fig. 8a, b). Interestingly, consistent with our previous finding that the GZMB+ cytotoxic subsets were specific in individual patients, the gene expression pattern of T cells also showed patient specificity instead of disease stage specificity (Supplementary Fig. 8b). Next, we examined the gene expression of T cell functional markers among T cell subsets. T cell activation markers CD27, CD28, and CD69 showed high expression in the non-cytotoxic subsets. Genes encoding for MHC type II proteins, which have been reported to inhibit T cell proliferation and cytokine production, showed high expression in the cytotoxic subsets. The immune checkpoint gene TIGIT was expressed in cytotoxic CD8+ T cell subtypes. CTLA4 was expressed in CD4 helper cells. LAG3 was sporadically expressed in the cytotoxic subsets (Supplementary Fig. 8c). LAG3 had increased expression in patients at late stage compared with NBM in both the CD8 GZMK+ and CD8 GZMB+ subsets (Supplementary Fig. 8d). However, we did not see a difference in the expression of TIGIT and MHC class genes. Of note, TIGIT showed increased protein levels in our CyTOF analysis (Fig. 5f), thus TIGIT may be post-transcriptionally regulated in the TME. These data indicate that T cell subsets are altered in MM and its precursor stages, with a shift from GZMK+ CD8 to GZMB+ CD8 subsets. The GZMB+ CD8 subsets exhibit patient specificity in the TME, whether those patient-specific T cell subsets are related to specific tumor clones due to antigen activation needs further investigation.

Composition alteration of MM subclones in paired samples shows patient-specific heterogeneity in progressed disease

Lastly, we studied the evolution of PC subclones during MM progression. FlowSOM analysis of CyTOF data identified two NPC MCs (CD27- and CD27+) and two MPC MCs (CD27-CD56- and CD27-CD56+). The proportion of CD27-CD56- MPCs showed a 55% increase in SMM and a 10.6-fold increase in NDMM compared to NBM. The proportion of the CD27-CD56+ MPCs increased 10.0-fold in SMM and 24.3-fold in NDMM compared with NBM (Fig. 7a). It suggests that the CD27-CD56+ MPCs drive early tumor expansion. No significant difference was observed in the expression of complement inhibitory proteins CD55 and CD59, and immune checkpoint ligand PD-L1 of MPC MCs among NDMM, SMM, MGUS, and NBM (Supplementary Fig. 9), which may be due to the remarkable heterogeneity in the expression of these markers within the groups. The gene expression pattern of PCs cluster in scRNA-Seq analysis showed that PCs in healthy donors exhibited poly-PC signatures, but PCs in diseased patients showed an increase in percentage and intensity of genes encoding for either immunoglobulin heavy chain IgG, IgA or IgD, which were consistent with their pathological examination results (Fig. 7b, Supplementary Fig. 10, Supplementary Table 4). We also observed that PCs showed increased IFN-response genes in NDMM and part of the SMM patients (Fig. 7b). Moreover, we observed that genes related to cell division (EDNRB, CCND1, CCND2), Wnt signaling (DKK1, FRZB), and cell adhesion (CNTN5, NCAM1 (encoding CD56), PCDH9, CADM1) were significantly increased in individual patients (Fig. 7b), indicating that MM and its precursor stages can be grouped into different subtypes based on their gene expression. Notably, many of those genes have been identified in CD138+ MM cell GEP analysis for the characterization of 7 MM groups [32]. Next, we subclustered plasma cells into 24 plasma cell subtypes based on their gene expression (Fig. 7c, Supplementary Fig. 11). Stable MGUS paired samples (P1, P2) showed few compositional alterations of the MM clone. Two of the five paired progression samples showed expansion of the dominant subclone (P3, P6); 3/5 Progressive paired samples showed alteration of dominant subclone (P4, P5, P7) (Fig. 7d, Supplementary Fig. 12). Whether the uncontrolled clonal expansion and composition alterations of MM cells in the paired progression samples can be explained by the lack of adaptation of cytotoxic T cells during disease progression needs biological examination.

Fig. 7: Composition alteration of MM clones in paired samples shows patient-specific heterogeneity in progressed disease.

a Box plot showing plasma cell (PC) composition changes among different disease stages by CyTOF analysis. p values for each cell type were calculated using Kruskal–Wallis with Dunn’s multiple comparisons test. b Dot plot showing gene expression of plasma cells for individual patients and healthy donors. c UMAP plots of plasms cell subtypes. d UMAP plots showing plasma cell subclusters in paired samples (left). Density map on the UMAP plot showing representative marker gene expression among plasma cell subtypes (right).