Cross center single-cell RNA sequencing study of the immune microenvironment in rapid progressing multiple myeloma

Overview of the CD138− cell fraction

ScRNA-seq data from the CD138- cells fraction of 48 BM biopsies from 18 patients were combined into an integrated dataset of 102,207 single-cells generated at three different medical centers. Low-quality cells are filtered out by removing cells with <200 unique genes, <500 UMI reads, and >30% mitochondrial UMIs. The malignant plasma cells, identified based on marker genes expression (MZB1+, JCHAIN+, SDC1+) (Supplementary Fig. 1a, b), heterogeneity as well as CNV analysis, were filtered out. The filtering of low-quality and malignant plasma cells resulted in a final integrated dataset of 90,502 CD138- single-cells.

Bone marrow cells were clustered based on gene expression profile (Fig. 1a) and annotated to cell types from various lineages using canonical marker genes: erythrocytes (HBA+, HBD+), erythroblasts (BLVRB+, PRDX2+), T-cells (IL7R+, CD3D+), CD8+ T-cells (CD8A+, CCL4+, GZMK+), CD8+ effector T-cells (GNLY+, GZMB+, CD3D+), NK Cells (GNLY+, GZMB+, CD3D-), monocytes/macrophages (CD14+, CD68+), CD1c+ dendritic cells (CD1c+), Granulocyte-Macrophage Progenitors-GMP (ELANE+, MPO+), plasmacytoid dendritic cells-pDC (IRF8+, MZB1+), hematopoietic stem cells-HSC (CD34+), Pro B-Cells (IGLL1+), and B-Cells (MS4A1+, CD79A+) (Fig. 1a, b). All the defined 24 cell types were detected in samples from both NP and RP groups (Fig. 1c). Most of the cell type clusters contained samples from multiple patients with none of the clusters being dominated exclusively by a single patient or sample (Fig. 1d). On average, CD4+ T-cells were the largest cluster among all patients, followed by erythrocytes, myeloid cells, CD8+ T-cells, and B-cells (Fig. 1e). Some cell types such as CD8+ exhausted T-cells and Memory B-cells are slightly elevated in RP, though currently these differences are not significant as a proportion of all cells. Further generation of unbiased cell type signature based on supervised analysis (P < 0.01, FC > 2) identified top markers that generally correspond with what is expected for each cell type (e.g., GNLY for CD8+ effector T-cells, LTB for CD4+ naive T-cells) (Fig. 1f). We observed overexpression of some immune markers such as IGHA1 that correspond to samples from an individual patient and may represent contamination from apoptotic plasma cells. Additional more stringent mitochondrial cutoff (filtering out cells with >20% mitochondrial DNA content) analysis depicted similar clustering patterns across all cell types (Supplementary Fig. 2a) and within the T-cell, Myeloid, and B-cell compartments (Supplementary Figs. 35). We also observed a significantly higher ratio of exhausted T cells (P = 0.049) (Supplementary Fig. 3b) as well as a significant enrichment of M2 macrophages (P = 0.049) (Supplementary Fig. 4b) in the RP compared to the NP samples.

Fig. 1: Single-cell profiling of bone marrow from MM patients with rapid (RP) and no progression (NP).figure 1

Clinical samples with rapid and no progression were identified from the Multiple Myeloma Research Foundation (MMRF) CoMMpass study, a longitudinal genomic study of patients with newly diagnosed, active multiple myeloma (NCT01454297). In the study, 48 bone marrow aliquots from 18 patients diagnosed with Multiple Myeloma (MM) were processed for scRNA-seq at three medical centers, Beth Israel Deaconess Medical Center (BIDMC), Mount Sinai School of Medicine (MSSM), and Washington University (WashU). Patients are classified as either rapid progressors (RP) or non-progressors (NP) based on the rate of disease progression, <18 months or >4 years post-diagnosis, respectively. a Uniform manifold approximation and projection (UMAP) embedding of scRNA samples across all patients consisting of 90,502 high quality single-cells portioned into 24 cell types. Plasma cells were removed prior to embedding. These clusters are colored based on canonical cell types based on the expression of marker genes that include erythrocytes (HBA+, HBD+), erythroblasts (BLVRB+, PRDX2+), T-cells (IL7R+, CD3D+), CD8+ T-cells (CD8A+, CCL4+, GZMK+), CD8+ effector T-cells (GNLY+, GZMB+, CD3D+), NK Cells (GNLY+, GZMB+, CD3D-), monocytes/macrophages (CD14+, CD68+), CD1c+DC (CD1c+), Granulocyte-Macrophage Progenitors-GMP (ELANE+, MPO+), plasmacytoid dendritic cells-pDC (IRF8+, MZB1+), HSC (CD34+), Pro B-Cells (IGLL1+), and B-cells (MS4A1+, CD79A+). b Dot Plot depicting expression profile of markers genes used for annotating different cell type clusters. The over and under expression of specific markers is shown by red and cyan colors, respectively. c Clinical phenotype based split UMAP showing the distribution of cell types in the RP and NP groups. There is slightly elevated NK, CD8+ effector, Pro B-Cell, GMP, and pDC counts in the NP group, while RP samples show elevation of exhausted T-cells, naive and memory B-cells, and erythrocytes. d A stacked bar plot showing the relative patient contribution to each individual cell type cluster. The patients from RP and NP groups are shown with shades of red and blue, respectively. Each cluster depicted the varying levels of contribution from multiple samples of NP and RP groups. e Comparative analysis of cell types enriched in the RP and NP groups. Each bar plot depicts the mean and ±standard error of the mean in NP and RP groups. Each dot represents an individual patient sample. f A heatmap displaying the top markers expressed by each cell type. Columns represent individual cells, grouped by cell type, while rows display individual genes. Horizontal colored bars above the heatmap indicate the cell type, with the legend on the right listing the cell type for each colored bar. Cell type labels are also displayed above their corresponding bar for all cell types except for the three smallest populations (M2 macrophages, MDSCs, and stromal cells). Relative gene expression is shown in pseudo color, where blue represents downregulation, and red represents upregulation. Top markers generally correlate with well-established canonical markers for each cell type.

To assess the potential malignancy of our plasma cell population (Supplementary Fig. 1b), we performed copy number variation (CNV) analysis17 using either normal plasma cells from the Human Cell Atlas (HCA) Census of Immune Cells18 or the naive and memory B cell populations from the current study dataset as reference cells. The clustering and UMAP analysis depicted that plasma cells formed multiple patient clusters with no co-clustering with the mature B cell populations indicating heterogeneity and difference in transcriptomes profiles (Supplementary Fig. 6a). CNV analysis identified multiple amplifications and deletions in a patient-specific manner on several chromosomes as well as promiscuous deletions on chromosomes 1 and 6 in the plasma cells (Supplementary Fig. 6b). These deletions/amplifications in the plasma cells as compared to B cell populations likely correspond to the malignant phenotype of plasma cells. Further comparative analysis of plasma cells from this study with the normal plasma cells from HCA also depicted significant heterogeneity and CNVs pointing toward the malignant phenotype of plasma cells identified in this study (Supplementary Fig. 7).

CD138− microenvironment cells exhibited similar single-cell profiles for three different testing centers

To assess the center/sample processing technical variations in the single-cell profiles, samples were processed from the same patients at three different medical centers/universities (BIDMC, WashU, and MSSM). All the samples were processed with droplet-based single-cell barcoding techniques for scRNA-seq alone (WashU, MSSM) or CITE-Seq (BIDMC). Twenty samples from 18 patients were processed at BIDMC, seven samples from seven patients were processed at MSSM, and 21 samples from 17 patients were processed at WashU (Supplementary Table 1). Samples had similar age distributions between both RP and NP groups, however, RPs had higher risk assessments at diagnosis as defined by the IMWG risk class (Supplementary Table 1).

Prior to the batch correction, most of the cell types from the three centers depicted similar clustering patterns with subtle variations among processing centers (Supplementary Fig. 8a). Shannon’s entropy was computed per cell to assess the degree of mixing of samples from centers, and low entropy values, indicating poor mixing, can be observed in MSSM samples, and in the ‘CD4+ T-cell’ class from BIDMC (Supplementary Fig. 9b, d). The batch effect correction improved the entropy values indicating better mixing of cells from different centers (Supplementary Fig. 9c, e). The evaluation of batch effect and corrections in the other cellular compartments (monocytes/macrophages, B cells) also depicted that batch effect correction improved the mixing of cell types/subtypes from different processing centers effectively alleviating site-dependent batch effect (Supplementary Figs. 1011). All major cell types (Fig. 1a) are uniformly detected in samples processed at each center as well as co-embed in the clustering indicating similarity in transcriptome profiles (Fig. 2a). This includes even rare cell populations such as the stromal and CD4+ regulatory T-cells, which represent <1% of the total cells profiled. Samples from MSSM and WashU had a similar ratio for all cell types. The comparative analysis of cellular proportions among the centers depicted subtle variation, which might be due to differences in the number of cells captured for single-cell assay, sequencing depth, and assays performed. All three centers captured striking similar proportions of pDCs, HSCs, B cells, and stromal cells (Fig. 2b). The samples processed at BIDMC depicted variation in the proportion of T/NK and myeloid cells compared to MSSM and WashU.

Fig. 2: Comparison of scRNA profiles of samples processed at three different centers.figure 2

Bone marrow aspirates from the same set of patients were processed at three different centers, BIDMC, MSSM, and WashU, and analyzed using a uniform bioinformatics workflow for comparative analysis. The comparative analysis was performed on 20 samples processed at BIDMC, 7 processed at MSSM, and 21 processed at WashU. a Split UMAP based on sample processing centers of scRNA samples. All major cell types are captured in the single-cell profile from each center. Clusters are colored based on cell types identified in Fig. 1a. b Comparative analysis of cell type proportion across centers. Each bar represents the mean ratio for a given cell type for all samples processed at a specific center. Error bars show the standard error of the mean. Individual dots represent individual patient cell type ratios. Samples from MSSM and WashU had similar ratios across the cell types. BIDMC shows a higher proportion of CD4+ T-cells and a lower ratio of Myeloid cells. c Comparative analysis of canonical cell type-specific markers across three centers. Most of the cell type defining markers are concordantly expressed across cell types indicating strong similarity in the single-cell profiles generated across centers. BIDMC, which performed CITE-Seq, tends to have higher percent.mt relative to other centers this might be due to longer processing time for CITE-Seq due to antibody labeling. d Violin Plots comparing the expression of various cell markers among different centers. Overall, the level of expression of these markers are consistent among centers indicating no batch effect or center-based expression artifact. e A Circos plot showing the correlation between expression profiles of cell types profiled at different centers. The individual cell types across centers depict significant similarity in the expression profiles. Some cell types with lower correlations include CD4+ memory T-cells from BIDMC, and monocytes and CD1c+ DCs from WashU.

Samples processed at BIDMC showed a higher proportion of CD4+ T-cells and a lower proportion of myeloid cells relative to samples processed at other centers (Fig. 2b). Some of these differences may be explained by the additional processing time required for BIDMC samples (simultaneous measurement of surface marker expression along with gene expression via the CITE-seq approach). A high proportion of the CD4+ T-cells in the BIDMC group belong to the general “CD4+ T-cell” clusters (Supplementary Fig. 12a, b), which tend to consist of lower viability samples with higher mitochondrial reads compared to other clusters. Furthermore, the comparative analysis of cellular subtypes for myeloid and B cells depicted striking similarities among the three centers (Supplementary Fig. 12c, d).

Further correlation of cell type proportions of individual samples processed at different centers depicted no specific pattern (Supplementary Fig. 13). To further explore the impact of sample processing from different centers on gene expression, we performed a comparative gene expression profile analysis of cell types and subtypes obtained from the three centers (Supplementary Fig. 14). Without batch correction, correlation is primarily driven by the cell compartment over the processing center, with T-cells, myeloid cells, B cells, and Erythroid cells clustering strongly regardless of the processing center. Within these compartments, there is a similar pattern as observed in the uncorrected embeddings, that is that WashU and BIDMC samples tend to have higher correlations relative to MSSM samples. MDSCs, a very small myeloid population, show the largest differences between centers. This appears to be driven by the lower viability of this cluster based on mitochondrial expression, along with a low cell count.

Next, we assessed the similarity in canonical markers expression for various centers. The canonical markers for each cell type depicted similar expressions for most cell types (Fig. 2c), with the primary exception being the previously noted MDSC cluster. The consistency of key marker expression ensures that cell types can be identified reliably across all three centers, and comparative analysis can be performed among the samples generated at different centers. This is further evident in violin plots of select markers illustrating the consistency in both the average expression and the distribution of gene expression within a cell population (Fig. 2d). The primary difference in gene expression between the centers is higher percentages of mitochondria-associated transcripts from samples processed at BIDMC (Supplementary Fig. 15). This may be attributed to a decrease in cell viability associated with the longer processing times for cell surface protein labeling before proceeding with the single-cell gene expression library preparation protocols at BIDMC.

To further validate the consistency of cell type labeling across centers, we assessed the similarity of the differentially expressed markers for each of the 24 cell types across each center. To achieve this, cells from each center were subset, and the top differentially expressed markers for each cell type with respect to all other cells from the same center were identified. Cell types from each center were clustered hierarchically based on the binary distance between these differentially expressed markers (Supplementary Fig. 16). We observed that overall, the differentially expressed cell type markers from all three centers strongly correlated with the matching cell type from other centers. This includes closely related cell types such as CD8+ T-cells and CD8+ exhausted T-cells (Fig. 2e). A few cell types, e.g., CD4+ Memory T-cells from BIDMC, HSCs, and Erythroblasts from MSSM, and certain myeloid cell types, such as CD1C+ DC, monocytes, and MDSCs, showed weaker correlation between centers. Some of these weaker correlations may be driven by large contributions from individual patient samples processed at certain centers. For example, a relatively large proportion of CD4+ Memory T-cells from BIDMC are from MMRF1413, which has a higher degree of plasma cell contamination, notably IGHA1, relative to other samples (Supplementary Fig. 17). In summary, the single-cell profiles generated at three different centers with different approaches have considerable similarities in their overall gene expression profiles. These results demonstrate that center-independent analysis can be implemented for a large cohort of samples processed at different centers after implementing appropriate batch correction approach.

Rapid progressors depicted significant enrichment of T cell exhaustion and attenuation of CD8+ effector T-cells

Focused analysis was carried out on T-lymphocyte and NK cells to understand their role in the progression of MM (rapid or non-progression). The analysis included 43,039 cells; 25,381 cells from NP samples and 17,658 cells from RP samples, which were identified as eight different subtypes of T and NK cells following further clustering (Fig. 3a). These cells were classified based on RNA and protein/ADT data: CD4+ T-cells (CD4+), CD4+ naive T-cells (CD4+, TCF7+, CCR7+, CD45RO+), CD4+ memory T-cells (CD4+,IL7R+, CD45RA+), CD4+ regulatory T-cells (CD4+, FOXP3+), CD8+ T-cell (CD8A+, KLRB1+, IL7R+, GATA3+, GZMK low), CD8+ exhausted T-cells (CD8A+, GZMK+, TIGIT+), CD8+ effector T-cells (CD3D+, GNLY+, GZMB+), and NK cells (CD3D-, GNLY+, GZMB+) (Fig. 3b).

Fig. 3: Comparative analysis of T and NK cell subpopulations in multiple myeloma patients with rapid- and no- progression of the disease.figure 3

a A UMAP displaying the T-cell subclusters split based on clinical groups (i.e., NP, RP). Subclusters were manually labeled as CD4+ T-cells (naive, memory, regulatory), CD8+ T-cells (memory, exhausted, effector), or NK cells based on the expression of specific markers. Limited CITE-Seq data from BIDMC was used to confirm some cellular annotations. NK and CD8+ effector T-cells show elevated counts from NP samples, while RP samples contain higher counts of CD8+ exhausted T-cells. b Dot plot demonstrating the expression profile of key markers for each T-cell subtype from both the scRNA-seq and CITE-Seq (ADT) assays. Markers for cell types used were CD4+ T-cells (CD4+), CD4+ naive T-cells (CD4+, TCF7+, CCR7+, CD45RO+), CD4+ Memory T-cells (CD4+,IL7R+, CD45RA+), CD4+ regulatory T-cells (CD4+, FOXP3+), CD8+ T-cell (CD8A+, KLRB1+, IL7R+, GATA3+, GZMK low), CD8+ exhausted T-cells (CD8A+, GZMK+, TIGIT+), CD8+ effector T-cells (CD3D+, GNLY+, GZMB+), and NK cells (CD3D-, GNLY+, GZMB+). c The patient contribution to each cell type cluster indicates that most of the clusters consist of cells from multiple patients. The patients from RP and NP groups are shown with shades of red and blue respectively. d Comparative analysis of the T-cell types enriched in the RP and NP T-cell subsets. Each bar plot depicts mean and standard error of mean. Significant enrichment (P = 0.022) of the CD8+ exhausted T-cells was observed in the RP population. e Differential expression analysis of the three CD8+ T-cell subtypes (CD8+ memory T-cells, CD8+ exhausted T-cells, CD8+ effector T-cells). Columns represent individual cells, grouped by cell type, while rows display individual genes. CD8+ T-cells depicted upregulation of markers related to T-cell memory, such as IL7R, but has under-expression of cytotoxic markers such as GZMK, GZMB, or NKG7 as compared to other CD8 T-cell subtypes. CD8+ exhausted T-cells depicted upregulation of GZMK and multiple genes related to chemokine signaling, such as CCL3 and XCL2. CD8+ effector T-cells showed downregulation of GZMK and upregulation of GZMB and GZMH, along with cytotoxic markers such as PRF1 and GNLY. f Comparative analysis of the CD8+ subset between NP and RP groups. Differential expression analysis was performed based on the Wilcoxon rank sum test of NP and RP CD8+ T-cells. NP CD8+ T-cells showed upregulation of markers related to NK cells, such as GZMB and GZMH. RP CD8+ T-cells instead show upregulation of the CD8+ exhausted T-cell markers, specifically chemokines like CCL3L1 and XCL2. These differences are reflected in the average cell type ratios in the CD8+ subset in NP and RP samples. The ratio of these CD8+ exhausted T-cells to the GZMB+ effector cells is significantly higher in RP samples (P = 0.048). g Expression profile of markers of exhaustion in the CD8+ and NK subsets. CD8+ exhausted T-cells, predominantly found in the RP group, show the highest expression of TIGIT and EOMES in the RNA assay relative to other CD8+ T-cells. CD160 is detectable in CD8+ T-cells and CD8+ exhausted T-cells, though only in samples from the RP population. CITE-Seq confirms elevated expression of TIGIT and PD-1 (CD279) in the CD8+ exhausted T-cell cluster, and general enrichment of exhaustion markers in the RP group over the NP group across multiple CD8+ cell types.

In both NP and RP groups, various subtypes of CD4+ T-cells were the dominant T-cells (~60%), with the remaining cells consisting primarily of CD8+ T-cells (30%) along with some NK cells (10%) (Supplementary Fig. 18). There are no significant differences in the CD4:CD8 ratio across clinical groups. Comparative analysis of cellular proportion among RP and NP groups depicted a significantly higher proportion of CD8+ exhausted T-cells in the RP samples (P = 0.022), whereas NK and CD8+ effector T-cells were non-significantly enriched in the NP samples (Fig. 3c, d). The higher enrichment of CD8+ exhausted T-cells in the RP samples indicates that RP patients with rapid progression of the disease have a prevalence of exhausted T cells even at the baseline/diagnosis. Pathway analysis was performed to compare the general CD4+ and CD8+ populations between clinical groups. Pathway analysis on these CD4+ T-cells depicted significantly higher enrichment of metabolic pathways such as Hallmark glycolysis and fatty acid metabolism in the patients with rapid progression. In addition, RP CD4+ T-cells depicted enrichment for downstream targets of both interferon alpha and gamma (Supplementary Fig. 19a). The pathway analysis on CD8+ T subset showed that RP samples were enriched in xenobiotic metabolism, while NP samples were enriched in TNFα signaling (Supplementary Fig. 19b). Survival analysis was performed using the MMRF CoMMpass dataset in the Survival Genie tool19, and the metabolic pathways which are enriched in rapid progressors, such as, fatty acid metabolism, glycolysis, and xenobiotic metabolism, were found to be associated with poor overall survival (Supplementary Fig. 19c).

We further investigated the gene expression profiles of the three CD8+ T-cell subtypes, revealing a range of cells with various degrees of exhaustion or cytotoxic activity. The CD8+ memory subset had high expression of markers associated with memory T-cells, such as KLRB1 and IL7R, but had low expression of cytotoxic markers such as GZMK, GZMB, or NKG7. CD8+ exhausted T-cells, in addition to elevated expression of TIGIT and multiple genes related to chemokine signaling, such as CCL3 and XCL2, had enriched expression of GZMK relative to other T-cell clusters. CD8+ effector T-cells show a shift away from GZMK and towards GZMB and GZMH, along with enrichment of cytotoxic and NK cells markers such as NKG7, GNLY, and PRF1 (Fig. 3e). CD8+ T-cells from the NP group show enrichment of markers related to the NK cells (GZMB, GZMA), whereas the RP group CD8+ cells show enrichment of many of the CD8+ exhausted T-cell markers, such as chemokine signaling genes XCL2 and CCL3L1 (Fig. 3f). The ratio of these more cytotoxic cells to the CD8+ exhausted T-cells was significantly enriched in NP group (P = 0.048). As aforementioned, the CD8+ exhausted T-cell subset shows higher expression of TIGIT compared to other CD8+ cells (Fig. 3g), along with other exhaustion markers such as EOMES and CD160 which are primarily found in RP samples. Additional CITE-Seq data (ADT) further confirms enriched expression of TIGIT and PD-1 (PDCD1/CD279) in the CD8+ exhausted T-cells, and that expression of TIGIT and PD-1 are generally enriched in all RP CD8+ T-cell types compared to NP CD8+ T-cell types. Overall, the CD8+ T-cell compartment shows a shift towards more exhausted effector-memory T-cells in RPs, compared to a more cytotoxic effector phenotype in NPs. Lastly, differential expression between clinical groups within the GZMB + CD8+ effector T-cells subset showed significant enrichment of multiple genes related to cytotoxicity and T-cell activation in NP samples (GZMB, GNLY, TNF), potentially indicating a better effector function and T-cell activation in NPs (Supplementary Fig. 20, Supplementary Table 2).

Rapid progressors depicted a higher proportion of alternatively activated M2 macrophages along with activation of complement cascade and lipid processing pathways

To study the alterations in cells of the myeloid lineage, focused analysis was performed after subsetting the myeloid cells. The myeloid subset comprised 16,245 cells of which 10,517 were derived from NP samples and 5728 from RP samples. NP samples show a much higher count of most myeloid cell subtypes compared to RP samples, with a nearly 2:1 ratio of myeloid-lineage cells overall. Based on the correlation of gene expression profiles, myeloid cells clustered into seven different subtypes (Fig. 4a). Cells were annotated as Granulocyte-Macrophage Progenitors (GMP) (MPO+, ELANE+, MKI67+), monocytes (CD14+, S100A9+, S100A12+), M1 macrophages (CD14+, CD44+, VCAN+, ITGAX+, CD86+), M2 macrophages (CD163+, MRC1+), MDSCs (HLA-DRA low, ITGAM+, ARG1+), CD16+ monocytes (CD14-, FCGR3A+), or CD1c+ DCs (CD1c+, CLEC10A+, MHC-II High) (Fig. 4b). Most of the myeloid cell clusters consist of cells from both NP and RP groups without any patient-specific clusters (Fig. 4c). Comparative analysis of cellular proportion across patient groups reveals significant enrichment of the M2 macrophage cluster in the RPs (P = 0.046) (Fig. 4d). However, NPs depicted higher enrichment of immature cell types, such as GMPs, and CD1c+ DCs relative to RPs. To explore pathways level dysregulation (if any) among various subtypes of myeloid cells, pathways enrichment analysis was performed. Similar to the T-cell subsets, monocyte, and macrophages from RP samples show enrichment for pathways related to IFNα and IFNγ, whereas NP samples predominantly show enrichment for TNFα via NFkB-related signaling pathways (Fig. 4e, f).

Fig. 4: Comparative analysis of the “monocyte and macrophage” and “GMP” immune microenvironment cell subpopulations in multiple myeloma patients with rapid- and no- progression of the disease.figure 4

a A UMAP displaying the monocyte and macrophage subcluster split based on clinical groups (NP and RP). Subclusters were labeled as either Granulocyte-Monocyte Progenitors (GMP), monocyte, CD16 + monocytes, M1 macrophages, M2 macrophage, MDSCs, or CD1c + dendritic cells (DC) based on expression of specific markers. GMP and CD1c + DCs show elevated counts in NP samples. b Dot plot demonstrating the key markers for the monocyte and macrophage subtypes. Markers to identify cell types include GMP (MPO+, ELANE+, MKI67+), monocytes (CD14+, S100A9+, S100A12+), M1 macrophages (CD14+, CD44+), M2 macrophages (CD163+, MRC1+), MDSCs (HLA-DRA low, ITGAM+, ARG1+), CD16 + monocytes (CD14-, FCGR3A+), and CD1c+ DCs (CD1c+). c The patient contribution to each cell type cluster indicating most of the clusters consist of cells from multiple patients. The patients from the RP and NP groups are shown with shades of red and blue. Overall, the NP group had a higher proportion of monocytes and macrophages relative to the RP group. d Comparative analysis of the myeloid cell types in the RP and NP myeloid subset. Each bar plot depicts the mean proportion of a specific cell type across clinical groups, with error bars displaying standard error of the mean. Individual dots show individual patient samples. M2 macrophages were significantly enriched (P = 0.045) within the RP population. e Pathway enrichment analysis on the monocyte and macrophage clusters. The Violin plots display the ssGSEA enrichment score of significantly differentially enriched pathways/gene sets between RP and NP groups. The RP group showed significant enrichment of interferon alpha and interferon gamma signaling pathways, while the NP group showed enrichment for TNF signaling and epithelial-mesenchymal transition pathways. f A bar graph displaying the top differentially enriched genesets of the monocyte and macrophage clusters based on FDR analysis between NP and RP is also shown. g A heatmap, displaying the top differentially expressed markers genes for NP and RP M1 macrophages. Columns represent individual cells, grouped by the RP or NP clinical groups, while rows display individual genes. Relative gene expression is shown in pseudo color, where blue represents downregulation, and red represents upregulation. h Selected pathways that are significantly (P < 0.01) enriched in the markers differentially expressed in the RP and NP M1 macrophage groups. Each bar represents a pathway with significant activation and inhibition in the RP group based on Z-score calculated using the IPA analysis platform. The pathways that are significantly activated (Z-score > 2) and inhibited (Z-score < −2) in the RP group are shown with orange and blue bars, respectively. i A heatmap, displaying the top differentially expressed genes for M1 and M2 macrophages. Columns represent individual cells, grouped by the type of macrophage (i.e., M1, M2), while rows display individual genes. Relative gene expression is shown in pseudo color, where blue represents downregulation, and red represents upregulation. j Selected pathways that are significantly (P < 0.01) enriched in the markers differentially expressed in the M1 and M2 macrophages. Each bar represents a pathway with significant activation and inhibition in the M1 macrophages group based on Z-score calculated using the IPA analysis platform. The pathways that are significantly activated and inhibited in the M1 macrophages are shown with orange and blue bars, respectively.

Differential gene expression analysis on M1 Macrophages, the largest myeloid cell subtype, was performed between NP and RP groups. RP samples showed enrichment of multiple genes consistent with enriched interferon signaling, such as EGR1, IFI6, and IFI16, while M1 Macrophages from NP samples show enrichment of proliferative genes, such as G0S2 and CEBPB (Fig. 4g). Pathway analysis on these DEGs further confirms the enrichment of interferon signaling in RP samples (Fig. 4h).

Differential gene expression between M1 and M2 macrophages was also performed as M2 macrophages are enriched in patients with rapid progression of the disease. M1 macrophages showed enrichment of inflammatory markers, including multiple S100A genes, whereas, M2 macrophages were enriched in multiple genes involved in the immune complement cascade (C1QB, C1QA), and lipid processing (APOE) (Fig. 4i). M2 macrophages also highly expressed VCAM-1, which is known to interact with myeloma cells through VLA-420. M2 macrophages also showed enrichment of tumor-promoting pathways, such as the platelet-derived growth factor (PDGF) signaling21 (Fig. 4j).

Non-Progressors depicted higher enrichment of immature B cell lineages

To determine dysregulations in B cell repertoire, we performed focused analysis on the clusters enriched with B-Lymphoid lineage cells. This subset consisted of 8009 cells; 4094 derived from NP samples and 3915 derived from RP samples, which formed clusters corresponding to four identified cell types (Fig. 5a). Cells were identified as naive B-cells (IGHM+, IGHD+, MS4A1+), memory B-cells (MS4A1+, CD27+), pre-B-cells (IGHM+, MS4A1-Low, IGLL1−), or progenitor B-cells (IGLL1+, RAG1+, RAG2+) (Fig. 5b). The profile of RP group reflects more mature B cell types (naive, memory), compared to the NP group. On the other hand, NP samples depicted higher enrichment of immature B cells (pre B-cells, pro B-cells) (Fig. 5c, d).

Fig. 5: Progenitor and precursor B-Cells are enriched in the multiple myeloma patients with no progression of the disease.figure 5

a A UMAP displaying the B Cell subcluster split based on clinical groups (NP and RP). Subclusters were labeled as either pro B-cell, pre B-cell, memory B-cell, or naive B-cell based on expression of specific markers. b A dot plot displaying the key markers used to identify each B cell subtype. Naive B-cells were identified via the expression of MS4A1, SELL, and LTB. Memory B-cells were identified by expression of CD27. Pre B-cells were identified by low MS4A1 expression. Pro B-cells were identified by the expression of RAG1, RAG2, and IGL11. c The patient contribution to each cell type indicates most of the clusters consist of cells from multiple patients. The patients from the RP and NP groups are shown with shades of red and blue respectively. The majority of naive and memory B-cells are derived from samples of the RP group, while the majority of pre and progenitor B-cells are derived from samples of the NP group. d Comparative analysis of the B-cell types/subtypes in the RP and NP clinical groups. Each bar plot depicts average abundance across patients and error bars show the standard error of the mean. Individual dots represent the abundance of a cell type within an individual sample. Within the B-cell subtypes, there is no significant difference in the average patient ratio. e A heatmap, displaying the top differentially expressed marker genes for naive B-cells, memory B-cells, pre B-cells, and pro B-cells. Relative gene expression is shown in pseudo color, where blue represents downregulation, and red represents upregulation. f Pathway analysis was performed on the differentially expressed markers between mature and memory B-cell groups versus the pre and pro B-cells. Selected pathways that are significantly (P-value < 0.01) enriched in these markers are displayed in the bar chart. Each bar represents a pathway with significant activation and inhibition in naive or memory B-cells based on Z-score calculated using the IPA analysis platform. The pathways that are significantly activated and inhibited in the naive and memory B-cells are shown with orange and blue bars respectively.

Differential gene expression was performed across the cell types to investigate the differences between the more mature and immature B-cell subtypes. Many of the DEGs are consistent with the standard markers for each category, such as IGHD and IL4R in the naive B-cells22, AIM2 in memory B-cells23, IRF4 in pre B-cells24, and IGLL1, DNTT, and VPREB1 in pro B-cells25 (Fig. 5e). In addition, the mature B Cells, show higher expression of MHC-II receptors and receptors involved in BAFF signaling, such as TNFRSF13B. Pathway analysis performed on these top markers show that the more mature B cell types are significantly enriched in BAFF and IL6 signaling (P < 0.01), both of which are involved in B cell and plasma cell maturation and proliferation (Fig. 5f). The decreases observed in the immature B cell population in the rapid progression group are consistent with previous studies showing the association of decreased progenitor B cell populations with symptomatic or relapsed MM26,27,28.

Cellular communication predicted dysregulation in BAFF, CCL, and IL16 signaling network

To further investigate potential signaling differences within the immune micro-environment in rapid or non-progressors we performed cellular communication analysis29. The potential for cell communication was measured by comparing the average expression of various ligands and receptors among the defined cell types. Using the ligand-receptor expression, an information-flow score is computed for each ligand-receptor pathway between all cell types to indicate communication patterns that are most likely to occur. If two cell types have a high expression of both the ligand and the receptor, then the information flow score between these cell types for this pathway will be higher relative to the other pathways.

First, the overall signaling structure was compared between the NP and RP samples. In both NP and RP, the CD8+ T-cells are the primary receivers of cell–cell signaling, and B-cells, monocytes, and macrophages act as primary senders (Fig. 6a). The overall communication structure is similar, and there are no cell types where ligand-receptor interactions are exclusively restricted to one cell population. CD8+ effector T-cells show increased signaling received in the NP group relative to the RP group from all cell types, while CD4+ memory T-cells, CD1c+ DCs, and CD16+ monocytes show higher signaling received in RP samples from other myeloid or B-lymphoid cells (Fig. 6b). These differences seem to be primarily driven by MHC-II–CD8a interactions in NP CD8+ effector T-cells, and MHC-I–CD4 interactions in CD4+ RP cells (Supplementary Fig. 21).

Fig. 6: Cell

留言 (0)

沒有登入
gif