Immunophenotypic changes in the tumor and tumor microenvironment during progression to multiple myeloma

Abstract

Investigation of the cellular and molecular mechanisms of disease progression from precursor plasma cell disorders to active disease increases our understanding of multiple myeloma (MM) pathogenesis and supports the development of novel therapeutic strategies. In this analysis, single-cell RNA sequencing, surface protein profiling, and B lymphocyte antigen receptor profiling of unsorted, whole bone marrow (BM) mononuclear cell samples was used to study molecular changes in tumor cells and the tumor microenvironment (TME). A cell atlas of the BM microenvironment was generated from 123 subjects including healthy volunteers and patients with monoclonal gammopathy of unknown significance (MGUS), smoldering MM (SMM), and MM. These analyses revealed commonalities in molecular pathways, including MYC signaling, E2F targets and interferon alpha response, that were altered during disease progression. Evidence of early dysregulation of the immune system in MGUS and SMM, which increases and impacts many cell types as the disease progresses, was found. In parallel with disease progression, population shifts in CD8 + T cells, macrophages, and classical dendritic cells were observed, and the resulting differences in CD8 + T cells and macrophages were associated with poor overall survival outcomes. Potential ligand-receptor interactions that may play a role during the transition from precursor stages to MM were identified, along with potential biomarkers of disease progression, some of which may represent novel therapeutic targets. MIF, IL15, CD320, HGF and FAM3C were detected as potential regulators of the TME by plasma cells, while SERPINA1 and BAFF (TNFSF13B) were found to have the highest potential to contribute to the downstream changes observed between precursor stage and MM cells. These findings demonstrate that myeloma tumorigenesis is associated with dysregulation of molecular pathways driven by gradually occurring immunophenotypic changes in the tumor and TME. Trial registration: This project has been registered at EudraCT (European Union Drug Regulating Authorities Clinical Trials Database) with protocol number NOPRODMMY0001 and EudraCT Number 2018-004443-23 on 12 December 2018.

Fig 1. Single-cell atlas of BMNCs shows proportional changes in key immune cell populations.

(A) Study design. BMNCs were collected from 123 subjects across 4 cohorts (HVs, MGUS, SMM, and newly diagnosed MM) and were further analyzed using scRNAseq, ADT-seq and scBCRseq. (B) Expression profile of marker genes/proteins for each cell type. The ‘TotalSeqC’ extension designates the proteins. (C) UMAP representation of the scRNA-seq dataset, colored by identified cell types. (D) Proportions of major cell types in each patient cohort (excluding plasma cells). (E) Proportions of minor cell types in each patient cohort (excluding plasma cells). Only cell populations where significant proportional changes were observed are displayed. For both sections (D) and (E), each dot represents an individual subject. The proportional differences between each disease state versus HVs are compared, and the significance was calculated using the Mann–Whitney U test: *p < 0.05, **p ≤ 0.005, ***p ≤ 0.0005. ADT-seq, antibody derived tag sequencing; BMNC, blood marrow mononuclear cell; DC, dendritic cell; HSPC, hematopoietic stem and progenitor cell; HV, healthy volunteers; mDC, myeloid DC; MGUS, monoclonal gammopathy of undetermined significance; MM, multiple myeloma; NK, natural killer; NKT, natural killer T cell; pDC, plasmacytoid DC; scBCRseq, single-cell B-cell receptor sequencing; scRNAseq, single-cell RNA sequencing; SMM, smoldering multiple myeloma; Treg, regulatory T cell.

Fig 2. Transcriptional similarities and differences were observed in plasma cells across disease stages.

(A) Cell number proportions of plasma cell clonality within subjects, obtained by BCR-seq. The heatmap on the top of the stacked bar plot shows the number of cells per sample. Expanded clones: per-subject clones with ≥ 10 cells; Other clones: identified clones with no expansion (marked with an asterisk); NA: cells with no clonal information captured by BCR-seq. (B) UMAP representation of plasma cell transcriptome, colored by cohort (left), subject (center) and expanded clonality (right) information. True: per-subject clones with at least 10 cells; False: identified clones with no expansion; NA: cells with no clonal information captured by BCR-seq. (C) Within and between cohort diversity is represented by hallmark gene set scores, where the x-axis is categorized by cohort and sorted based on mean gene set scores, whereas the y-axis is sorted by mean gene set scores in MM subjects. (D) GSEA results where each disease state is contrasted with HVs. Bolded gene sets highlighted in both panels (C) and (D) were significant in MGUS, SMM, and/or MM, and show gradual divergence at both patient and cohort level. (E) UMAP representation of transcriptome of plasma cells colored by gene set scores of hallmark gene sets enriched in precursor stages. The figure colors range from minimum value (blue) to the 99^th percentile of the value distribution (red). (F) Differential expression results of the contrast diseased (MGUS, SMM, MM) versus non-diseased (HVs) plasma cells. Points colored in red represent significantly down- (negative logFC) and upregulated (positive logFC) genes. Significantly down- and upregulated genes were identified based on adjusted p value < 0.05, absolute value of logFC > 1, and logCPM > 1. BCR-seq, B-cell receptor sequencing; GSEA, gene set enrichment analysis; HV, healthy volunteers; IFN, interferon; logCPM; log counts per million; logFC, log-fold change; MGUS, monoclonal gammopathy of undetermined significance; MM, multiple myeloma; NA, not available; NES, normalized enrichment score; SMM, smoldering multiple myeloma.NES, normalized enrichment score; SMM, smoldering multiple myeloma.

Fig 3. Proportional and transcriptional changes in the immune microenvironment showed early dysregulation of inflammatory response.

GSEA results of hallmark gene sets where the expression profile of each disease state is contrasted by HVs per immune cell type. On the x-axis, the cell types are grouped by their corresponding major cell-type category. The y-axis is sorted based on average NES score in MM. Gene sets, enriched along disease progression, are highlighted at the top. DC, dendritic cell; GSEA, gene set enrichment analysis; HV, healthy volunteers; MGUS, monoclonal gammopathy of undetermined significance; MM, multiple myeloma; NES, normalized enrichment score; SMM, smoldering MM.

Fig 4. Increase in pre-dysfunctional phenotype within CD8 + ATCs observed during disease progression.

(A) UMAP representation of the transcriptome of CD8 + T cells and a separate UMAP including only CD8 + ATCs each colored by cell groups in per-cell-type analyses. (B) Gene set scores and GSEA plots of TNFα signaling highlighting enrichment of TNFα signaling in pre-dysfunctional CD8 + ATCs. The figure colors range from minimum value (blue) to the 99^th percentile of the value distribution (red). In the GSEA analysis, the gene expression levels of pre-dysfunctional cells are compared against other CD8 + ATCs. (C) Gene expression profiles of marker genes define CD8 + ATC subtypes. The figure colors range from minimum value (yellow) to the 99^th percentile of the value distribution (red). (D) Proportional differences among cohorts for CD8 + ATC subtypes. Each dot represents an individual subject. The proportional differences between each disease state versus HVs are compared, and the significance was calculated using the Mann–Whitney U test: *p < 0.05, **p ≤ 0.005, ***p ≤ 0.0005. (E) Density plots showing the transcriptomic landscape of CD8 + ATCs per cohort. The highlighted area represents the region where the pre-dysfunctional cells reside. A transcriptional shift with disease progression was observed. ATCs, activated T cells; GSEA, gene set enrichment analysis; HV, healthy volunteers; MGUS, monoclonal gammopathy of undetermined significance; MM, multiple myeloma; SMM, smoldering MM; TNF, tumor necrosis factor.

Fig 5. Macrophages polarized towards an M1 phenotype during myeloma disease progression.

(A) UMAP representation of transcriptome of the macrophage population, colored by subtype annotations. (B) Expression profiles of marker genes that define macrophage subpopulations. (C) Gene set scores of IFNα and IFNγ highlight IFN-responsive regions in macrophages. The figure colors range from minimum value (blue) to the 99^th percentile of the value distribution (red). (D) Proportional differences among cohorts for macrophage subtypes. Each dot represents an individual subject. The proportional differences between each disease state versus HVs are compared, and the significance was calculated using the Mann–Whitney U test: *p < 0.05, **p ≤ 0.005, ***p ≤ 0.0005. (E) Density plots showing the transcriptomic landscape of macrophages per cohort. The highlighted area represents the region where M1 macrophages reside. (F) The distribution of BAFF gene expression within the macrophage population. The cut-off value was defined by the local minima at 0.92 to separate BAFF-low and BAFF-high populations. (G) Density plots of the BAFF-low and -high populations, showing dense regions on transcriptomic landscape for both categories. The highlighted area represents the region where the M1 macrophages reside. HV, healthy volunteers; IFN, interferon; MGUS, monoclonal gammopathy of undetermined significance; MM, multiple myeloma; SMM, smoldering MM.

Fig 6. CD1C + DCs activation and maturity increased in parallel with disease progression.

(A) UMAP representation of transcriptome of the CD1C + DC population colored by CD1C + DC subtypes and gene set scores of various inflammation-related pathways. The figure colors range from minimum value (blue) to the 99^th percentile of the value distribution (red). (B) Expression profiles of marker genes define CD1C + DC subpopulations. The figure colors range from minimum value (yellow) to the 99^th percentile of the value distribution (red). (C) Density plots showing transcriptomic landscape of CD1C + DCs per cohort. Highlighted area represents the region where the activated CD1C + DCs reside. (D) Proportional differences among cohorts for CD1C + DC subtypes. Each dot represents an individual subject. The proportional differences between each disease state versus HVs are compared, and the significance was calculated using the Mann–Whitney U test: *p < 0.05, **p ≤ 0.005, ***p ≤ 0.0005. DC, dendritic cell; HV, healthy volunteers; IFN, interferon; MGUS, monoclonal gammopathy of undetermined significance; MM, multiple myeloma; SMM, smoldering MM; TNF, tumor necrosis factor.

Fig 7. Survival analyses revealed associations between proportions of pre-dysfunctional cells and M1 macrophages with OS.

(A) Distribution of the proportion of pre-dysfunctional CD8 + T cells within the CD8 + ATC population in the CoMMpass dataset among NDMM samples. The samples were categorized into pre-dysfunctional high and low groups using the median cut-off (median = 0.4). (B) Kaplan-Meier survival curve demonstrates the association between CD8 + ATC proportion and OS (pre-dysfunctional high median = 1353 days, pre-dysfunctional low median not met, p = 0.0025). (C) Distribution of the proportion of M1 macrophages within the macrophage population in the CoMMpass dataset among NDMM samples. The samples were further categorized into M1 macrophage high and low groups using the median cut-off (median = 0.55). (D) Kaplan-Meier survival curve demonstrate the association between the proportion of M1 macrophages and OS (M1 high median = 1574 days, M1 low median not met, p = 0.036). ATC, activated T cells; NDMM, newly diagnosed multiple myeloma; OS, overall survival.

Fig 8. Ligand-receptor interaction modelling revealed potential interplay between tumoral plasma cells and immune populations.

(A) The circos plot (left) shows the interactions when plasma cells are selected as senders and the immune cells as receivers. The dot plot (right) shows significantly differentially expressed ligands/receptors in sender/receiver cells (p < 0.05). (B) The circos plot (left) shows the interactions when the immune cells are selected as senders and plasma cells as receivers. The dot plot (right) shows significantly differentially expressed ligands/receptors in sender/receiver cells (p < 0.05). (C) The activity of the ligands from immune cells and their regulatory potential over significantly upregulated and downregulated genes in plasma cells are shown. The order of the ligands (top to bottom) represents their rank on the prioritization. The highlighted ligands either had significant differential expression in the cell types of interest or were observed in our population shift study. All the plots in this figure are obtained using the results from the differential expression analysis of MM versus premalignant stages (MGUS & SMM) contrast. DC, dendritic cell; MGUS, monoclonal gammopathy of unknown significance; MM, multiple myeloma; NK, natural killer; SMM, smoldering MM; Treg, regulatory T cell.

Fig 9. Proposed model of MM progression with dynamic interplay between malignant plasma cells and immune populations.

Each cell type is represented with a different color: MM, purple; DCs, turquoise; macrophages, yellow, T cells, blue. The MM plasma cell ligands and receptors are represented in green. Ligands and receptors from the immune populations are represented in dark blue. The straight arrows show the direction of the transformation. The arrows with descending color represent the direction of the effect on transformation. The orange cloud represents the inflammatory environment. DC, dendritic cell; MM, multiple myeloma.