Single-cell analysis unveils distinct transcriptional alterations and cellular origins in IgD multiple myeloma

Multiple myeloma (MM) is the second most prevalent blood cancer, distinguished by the secretion of monoclonal immunoglobulin, leading to organ dysfunction, anemia, renal impairment, and bone lesions.¹ Despite remarkable advancements in therapeutic strategies, identifying and managing patients with poor outcomes under current regimens remains imperative. Immunoglobulin D (IgD) subtype, accounting for less than 2% of cases, exhibits inherently aggressive clinical traits: young age at presentation, higher incidence of extramedullary disease, osteolytic lesions, renal failure, etc.² Despite advancements in response evaluation criteria, IgD myeloma continues to pose a significant survival challenge, with a median overall survival (OS) of approximately 3 years, notably shorter than that of other subtypes.^3,4

To elucidate genetic alterations affecting the outcomes of IgD MM patients, we assembled 37 samples, comprising 13 IgD and 24 non-IgD MM samples. Bone marrow samples were obtained and used for sequencing. Samples were acquired with patients’ written consent in accordance with the Declaration of Helsinki with approval from the Ethics Committee of Shanghai Changzheng Hospital (2016SL019A). Eighteen samples were sorted by CD138 microbeads, and the other 19 samples were not sorted. We conducted both single-cell RNA sequencing and single-cell B-cell receptor (BCR) sequencing analyses, and found significant differences in tumor cells between IgD MM and non-IgD MM, as well as notable microenvironmental changes.

We processed and analyzed the gene expression matrix with the CellRanger pipeline and the Seurat package. After low-quality cells had been excluded, cell types were identified based on the high expression of marker genes. Copy number variant (CNV) alterations between IgD and non-IgD myeloma were identified using the InferCNV package. Consensus non-negative matrix factorization (cNMF) was used to unravel the intratumoral expression profiles and meta-profiles of MM. Single-cell V(D)J data were processed with the Dandelion pipeline. Normal cells from unsorted samples were integrated with the Harmony package to remove batch effects and explore the difference between cells from various samples. Cell–cell interaction analysis was conducted using the LIANA package.

First, we scrutinized the gene expression profiles and BCR features of tumor cells from both IgD MM and non-IgD MM. Each tumor predominantly harbored a main clone consistent with its clinical subtype (Online Supplementary Figure S1). However, beyond the main clone, numerous subclones were identified, potentially diverging from the clinical subtype identification, implying intra-tumor heterogeneity within patients. Given that MM arises from malignant plasma cells with a clonal B-cell origin, the emergence of subclones is generally uncommon. Through single-cell VDJ analysis, we confirmed the clonal nature of tumor cells, observing a predominantly monoclonal pattern in all 18 sorted samples, as expected. Remarkably, we identified instances in which cells exhibited identical VDJ regions but different immunoglobulin subtypes, suggesting the potential for antibody production after tumor formation (Figure 1A). We also analyzed the preferences of constant and variable region transcription genes in BCR across different groups, along with their physicochemical characteristics (Figure 1B). Clustering of main clonal sequences translated into corresponding amino acid sequences from different samples revealed that in the heavy chain, the CDR3 length was predominantly concentrated in four lengths: 15, 17, 18, and 20 amino acids for VDJ motifs. Additionally, motif enrichment analysis indicated exclusion of IgG subtypes from group 1, while two out of three IgD subtypes and two out of six IgA subtypes were enriched in group 1 motifs. Group 4 motifs were exclusively associated with IgG subtypes, while four out of nine IgG subtypes were enriched in group 5. Regarding the light chain, the CDR3 length was predominantly concentrated in three lengths: 12, 13, and 15 for VJ motifs. One IgD sample did not match any motif, while the other two samples were both enriched in group 3 VJ motifs. Furthermore, four out of six IgA subtypes were enriched in group 2 motifs. These unveiled that the BCR of IgD MM appeared to possess more similar physicochemical properties compared to other BCR types. However, further investigation with an expanded sample size is warranted to explore this phenomenon comprehensively.

In addition, we identified higher rates of alterations in 1q+ in IgD myeloma through the inference of CNV from single-cell RNA-sequencing expression data (Figure 2). Among MM cases with 1q+, gain(1q) refers to those with only one extra copy of 1q, resulting in three total copies, while amp(1q) denotes patients exhibiting amplification of 1q, characterized by the presence of two or more additional copies, totaling four or more copies.⁵ A recent systematic review of randomized controlled trials of MM focusing on 1q+ found that patients with 1q+ compared to those without have poorer outcomes.⁶ Notably, 1q21 amplification was present in 76.9% (10/13) of IgD patients compared to 62.5% (15/24) of non-IgD patients (P=0.47), potentially contributing to the poorer prognosis associated with this subtype.

Subsequently, we conducted a comparative analysis of the transcriptional features of MM through differentially expressed genes, cNMF-derived gene expression profiling, and enrichment analysis. Consistent with previous findings, we observed widespread involvement of endoplasmic reticulum stress and energy metabolism pathways in MM. Differentially expressed genes in IgD MM and non-IgD MM were mainly enriched in protein processing and energy metabolism pathways. Aligning with our gene expression profiling analysis, we identified a total of nine meta-programs associated with the disease, with the most significant one being meta-program 2, consistently observed in all analyzed samples, covering gene signatures related to protein processing functions (Figure 3A). However, upon examining the expression profiles of each sample in this meta-program, we found significant differences between IgD MM and non-IgD MM (Figure 3B). Notably, ribosomal proteins (RPL19, RPL28, RPLP0, RPL18A, RPS4X, RPS5, RPL7A, RPLP1) within this signature may be associated with bortezomib response⁷ (Figure 3C). This may explain why our prior studies found that IgD myeloma patients who received lenalidomide (an immunomodulatory drug) showed a better trend in median overall survival compared to patients who received bortezomib (a proteasome inhibitor),² highlighting the importance of appropriate treatment selection for this IgD subgroup. Comparing differentially expressed genes, pathways related to endoplasmic reticulum stress, autophagy, and unfolded protein response were again enriched between malignant cells in IgD and non-IgD MM (Figure 3D). Among the differentially expressed genes, the BCL-2 protein family serves as a core regulator in integrating stress signaling networks, regulating cell death, calcium homeostasis, the unfolded protein response and autophagy.⁸ Notably, we found higher expression of BCL2. Venetoclax, a BCL2-specific inhibitor, has shown promising results with higher response rates and longer progression-free survival for t(11;14) patients and those with high expression of BCL2.⁹ Thus, we speculate that IgD myeloma may be specifically sensitive to treatment regimens involving an immunomodulatory drug plus venetoclax. This area needs further clinical research.

The “seed and soil” hypothesis posits that the bone marrow microenvironment (“soil”) may create a supportive survival niche for tumor cells (“seeds”). This microenvironment consists of a complex array of cells. In our investigation, we delved into the characteristics of the bone marrow micro-environment in IgD MM. We noted a distinct composition of B cells in the bone marrow of IgD MM patients (Online Supplementary Figure S2A, B), with a higher proportion of immature B cells (P=0.014). The ratios of B-cell populations in non-IgD samples mostly aligned with previous findings, with 3% pro-B cells, 8% pre-B cells, 7% immature B cells, and 83% mature B cells.¹⁰ However, in IgD samples, the ratios were 12% pro-B cells, 40% pre-B cells, 21% immature B cells, and 27% mature B cells, indicating a higher proportion of immature B cells in IgD MM. Immature B cells in IgD MM displayed altered expression of genes characteristic of immature B cells such as BLK, BANK1, RAG2 (Online Supplementary Figure S2C). BANK1 expression is integral to pathways crucial for B-cell selection and function, such as BCR signaling and antigen processing and presentation via major histocompatibility complex (MHC) class II. The transient expression of pre-BCR and BCR signaling marks a critical checkpoint in B-cell development. Signals from the pre-BCR provide swift feedback regarding the functionality of the recombined Igμ gene, ensuring further maturation only for pre-B cells expressing a signaling-competent receptor.¹¹ These findings suggest that IgD MM status may diminish BCR signaling, thereby curtailing survival and positive selection.

Figure 1.Analysis of B-cell receptor profiles across 18 multiple myeloma samples. (A) Comparison of different classes of myeloma malignant cells from patient S_IgG_RR_02, highlighting VDJ and VJ sequences with differing amino acids marked in red. Clones sharing the same VDJ and VJ but with different heavy chain C domains are boxed, offering insights into clonal evolution and diversity within individual patients. (B) Motif analysis of VDJ and VJ sequences, with multiple myeloma (MM) subtypes, new diagnoses, and relapses annotated on the right with metadata, providing insights into the sequence motifs associated with different MM subtypes and disease stages. ND: newly diagnosed; RR: relapsed/refractory.

Figure 2.Copy number variants of all samples. Copy number variants of all samples inferred by inferCNV. On the left, a summary of key information for all multiple myeloma samples, including cytogenetic abnormalities detected by fluorescence in situ hybridization (FISH). On the right, copy number variants of all samples are depicted. An asterisk denotes cytogenetic abnormalities detected by FISH. AMP: amplification; Del: deletion; scRNA-seq: single-cell RNA sequencing; Amp Freq: amplification frequency; ND: newly diagnosed; RR: relapsed/refractory; Del Freq: deletion frequency; BCR_seq: B-cell receptor sequencing.

Figure 3.Transcriptomic characteristics of IgD and non-IgD malignant plasma cells and experimental validation. (A) Pairwise correlation clustering of 115 gene expression profiles derived from consensus non-negative matrix factorization analysis resulted in the formation of nine consensus modules. The biological significance of these modules was predicted through gene ontology (GO) analysis. (B) Comparison of the module score of M2 between IgD samples and non-IgD samples. Statistical significance is denoted by asterisks: *P<0.05; **P<0.01; ***P<0.001, ****P<0.0001. (C) Expression profiles of M2 marker genes in IgD and non-IgD samples. (D) GO enrichment analysis results for differentially expressed genes (DEG) in malignant plasma cells (IgD vs. non-IgD). (E) Differential gene expression analysis using the log-fold change expression versus the difference in the percentage of cells expressing the gene comparing IgD and non-IgD multiple myeloma in newly diagnosed patients (Δ Percentage Difference). The top ten upregulated and downregulated genes, as well as some genes related to protein processing, are labeled.

We also conducted an extensive analysis of interactions among different cell types, with a particular focus on those between B cells and tumor cells. Several factors show significant differences between IgD MM and non-IgD MM in all unsorted samples (Online Supplementary Figure S3A) as well as in newly diagnosed unsorted samples (Online Supplementary Figure S3B). Factor 4 (P<0.05) exhibited differences between the IgD and non-IgD groups for all unsorted samples and was primarily associated with myeloid cells, including pathways related to APP, GRN, CD23, MHC-I, and MHC-II (Online Supplementary Figure S3C). CD23 plays an essential role in B-cell responses and distinguishing various immature and mature B-cell subsets.¹² Remarkably, during the enrichment analysis of factor 4-related pathways, we discovered that the expression of the TNFSF13B-related pathway was downregulated in IgD samples. TNFSF13B, also known as B cell-activating factor (BAFF), enhances B-cell survival and serves as a key regulator of the peripheral B-cell population. BAFFR is closely associated with BCR signaling. As we mentioned earlier, the distinct composition of B cells in the bone marrow of IgD MM patients may be related to diminished BCR signaling and antigen processing and presentation via the MHC class II pathway, both of which are integral to pathways crucial for B-cell selection and function. Similarly, for newly diagnosed unsorted samples, factor 4 (P<0.05) presented differences between the IgD and non-IgD groups (Online Supplementary Figure S3B). It was primarily linked to B cells and included pathways related to MIF and HMGB1, which played an essential role in B-cell responses and distinguishing various immature and mature B-cell subsets^12-18 (Online Supplementary Figure S3D). The results further suggest that changes in B-cell subtype proportions in IgD MM may be related to the impact of the bone marrow microenvironment. Also, the interaction between the tumor microenvironment and tumors varies by subtypes, contributing differently to tumor development.

In summary, we have characterized both common and distinct transcriptomic and clonal features of IgD MM compared to other MM subtypes at the single-cell resolution. We found IgD MM displayed a higher frequency of 1q21 amplification. Additionally, IgD MM revealed disparities in the endoplasmic reticulum stress pathway and higher BCL2 expression independent of t(11;14), suggesting its potential sensitivity to venetoclax. IgD MM also exhibited notable microenvironmental changes, characterized by an increased presence of immature B cells, which may suggest a potential link to tumorigenesis. Furthermore, by analyzing motif patterns in VDJ/VJ, we found that the BCR of IgD MM tends to have a more similar clonotype with fewer subclones compared to non-IgD MM. While validation in larger cohorts and well-controlled randomized clinical trials is warranted, these findings illuminate the heterogeneity of MM and offer new perspectives for research and the clinical treatment of IgD MM.

Footnotes

• Received July 25, 2025
• Accepted October 21, 2025

Correspondence

*TL, XZ and XC contributed equally as first authors

J. DU - juan_du@live.com

J. LI - ljing@smmu.edu.cn

Disclosures

No conflicts of interest to disclose.

Contributions

Funding

This work was supported by the National Natural Science Foundation of China (82370206, 82170198), 2024 Shanghai Talents Plan, and the National Key R&D Program of China (2022YFA1305700).