Multi-omics reveal immune microenvironment alterations in multiple myeloma and its precursor stages

Abstract

Tumor immune microenvironmental alterations occur early in multiple myeloma (MM) development. In this study, we aim to systematically characterize the tumor immune microenvironment (TME) and the tumor-immune interactions from precursor stages, i.e., monoclonal gammopathy of undetermined significance (MGUS) and smoldering MM (SMM), to newly diagnosed MM, comparing these to healthy donors. Using CIBERSORT, mass cytometry (CyTOF), and single-cell RNA sequencing (scRNA-Seq), we examined innate and adaptive immune changes across these stages. We found a decrease in granulocytes in the TME predicts MM outcomes. HLA-DR is reduced in CD16⁺ monocytes and plasmacytoid dendritic cells, while myeloid dendritic cells show decreased expression of stress and immune-response genes. NK cells and CD8⁺ T cells shift from a GZMK⁺ to a GZMB⁺ cytotoxic phenotype in the TME, with increased inhibitory markers TIM3 and TIGIT. In paired samples, the proportion and gene expression pattern in patient-specific GZMB⁺CD8⁺ T cells remain largely unchanged despite MM progression. Our findings provide a comprehensive immune landscape of MM and its precursors, offering insights into therapeutic strategies. Enhancing neutrophil and NK cell cytotoxicity, tumor antigen presentation, and CD8⁺ T cell versatility in precursor stages may prevent MM progression.

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.

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 −log₁₀ 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.

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 −log₁₀ 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.

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.

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.

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.

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).