• We defined a subpopulation of malignant plasma cells associated with poor outcomes as EMICs.

  • LILRB4 mediated the migration of NCI-H929 cells and induced cell cycle arrest in the G1 phase in vitro.

Myelomatous effusion (ME) is a rare manifestation of extramedullary multiple myeloma (MM) with limited therapeutic options and poor outcomes. The molecular mechanisms underlying ME are incompletely understood. We profiled transcriptomes of bone marrow, peripheral blood (PB), and pleural effusion/ascites from 3 patients with ME using single-cell RNA sequencing analysis. We found that ME contained a higher percentage of cytotoxic T cells, whereas PB contained a higher proportion of naive T cells. Malignant cells varied within and between sites and patients in their expression of signatures. We identified a gene module highly expressed in intramedullary and extramedullary plasma cell clusters and defined cell clusters expressing this gene set as extramedullary-initiating cells (EMICs). This gene set was associated with increased cellular proliferation, involved in p53 signaling, and related to poor prognosis in MM. The transcriptional regulators E2F1, YY1, and SMAD1 were activated in EMICs. Leukocyte immunoglobulin–like receptor subfamily B4 (LILRB4) was upregulated in extramedullary EMICs. We confirmed that LILRB4 promoted MM cell migration in vitro. This study provided insight into the evolutionary mechanisms of ME and defined EMICs and LILRB4 associated with extramedullary development.

Multiple myeloma (MM) is characterized by the clonal proliferation of malignant plasma cells (PCs) and accounts for 10% of malignant hematopoietic diseases.1 The development of protease inhibitors, immunomodulatory drugs, and monoclonal antibodies has increased survival rates.2 However, MM is incurable, and disease progression results in drug resistance and poor outcomes.3 When a subclonal of PCs become independent of marrow microenvironment and egress out of bone marrow (BM), it is termed as extramedullary multiple myeloma (EMM). The overall incidence of EMM is 6% at the initial diagnosis and 6% to 20% during the follow-up.4 EMM can spread to skeletal tissue and extraosseous sites, including lymph nodes, liver, skin, central nervous system, pleura, and blood.5 Extraosseous MM is associated with worse outcomes.6,7 

Myelomatous effusion (ME) is manifested with the presence of neoplastic PCs in bodily fluids, including pleural effusion and ascites. ME occurs in <1% of MM cases.8,9 Although a few cases of remission with long-term survival have been reported,10,11 the median survival of patients with ME is only 4 months.12,13 The pathogenesis of ME has not been fully elucidated.14 The increased proliferation of PCs and the activation of antigen presentation, proteasomes, oxidative phosphorylation, glycolysis, and immune evasion in these cells are involved in the formation of EMM.15 Pleural and peritoneal lesions and hematogenous metastasis may be related to ME.16 However, the molecular mechanisms underlying ME need to be studied further.

This study included single-cell RNA sequencing (scRNA-seq) of malignant PCs as well as immune cells in 3 patients with MM who relapsed with ME. BM, peripheral blood (PB), and pleural effusion/ascites were collected. The transcriptome analysis of PCs in these tissues allowed for the identification of cell subclones involved in the migration to extramedullary sites. We determined the molecular mechanism of ME formation and the impact of target genes on myeloma invasiveness using in vitro experiments. Our findings underscore the importance of analyzing single-cell transcriptome profiles to identify PC clones able to migrate and survive in extramedullary sites.

Patients and sample preparation

Primary BM, PB, and pleural effusion/ascites samples were collected from 3 patients with MM with ME at the Department of Hematology of Jiangsu Province Hospital from January 2019 to December 2021. The main properties and quality statistics of the experiment are summarized in supplemental Materials. All procedures in this study were approved by the ethics committees of Jiangsu Province Hospital and in accordance with the Declaration of Helsinki. All the enrolled patients signed an informed consent.

Library preparation and scRNA-seq

Cells (1 × 105 cells per mL) were harvested and suspended in equal volumes of phosphate-buffered saline. Cells were loaded onto a microfluidic chip, and scRNA-seq libraries were constructed using a GEXSCOPE Single-Cell RNA Library Kit17 (Singleron Biotechnologies, Nanjing, China). Individual libraries were diluted to 4 nM and pooled, and the pools were sequenced using an Illumina NovaSeq 6000 platform to produce 150 bp paired-end reads.

Analysis of raw reads

Raw reads from scRNA-seq were processed to generate gene expression matrixes using CeleScope version 1.1.7 (https://github.com/singleron-RD/CeleScope). Low-quality reads were removed using CeleScope, and poly-A tails and adapter sequences were removed using Cutadapt2 version 1.17. Cell barcodes and unique molecular identifiers (UMIs) were extracted. High-quality reads were mapped to the reference genome GRCH38 (Ensembl version 92) using STAR version 2.6.1a.18 UMI counts and gene counts were grouped using featureCounts19 version 2.0.1 to generate expression matrixes.

Quality control, dimensionality reduction, and clustering

Cells with a range from 200 to 5000 genes, less than 30 000 UMI counts, and less than 50% mitochondrial content were retained; 37 731 cells passed quality control. Six samples from GEXSCOPE (MM01-ME, MM01-PB, MM01-BM, MM03-BM, MM03-ME, and MM03-PB) and 2 samples from 10X Genomics (MM02-BM and MM02-ME) were integrated using Harmony version 1.0. Dimensionality reduction, unsupervised clustering, and data were visualized using Seurat20 version 3.1.2. Gene expression was normalized using the Normalize Data function of the Seurat package, in which the number of UMIs of each gene was divided by the total UMIs of each cell, multiplied by 10 000, and log-transformed (ln [UMI-per-10 000 + 1]). Highly variable genes were identified using the Find Variable Genes function and analyzed using principal component analysis. Clustering with 20 principal components at a resolution of −0.8 was performed via graph-based clustering. Cells were visualized in a 2-dimensional space using a uniform manifold approximation and projection (UMAP) algorithm and the RunUMAP function of the Seurat package.

Cell type annotation

The cell types in each cluster were identified using canonical markers from literature data and the reference SynEcoSys database (Singleron Biotechnologies, Nanjing, China). Cell doublets were estimated based on the expression pattern of canonical cell markers. Clusters enriched with multiple cell type–specific markers were excluded from the analysis. Subtype analysis of major cell types is summarized in supplemental Materials.

Analysis of DEGs

Differentially expressed genes (DEGs) were determined using the FindAllMarkers function of the Seurat package (Wilcoxon rank-sum test with Bonferroni correction for multiple comparisons). For computing DEGs, genes expressed in at least 10% of PCs in either of the 2 samples being compared were selected for analysis. Cut-off criteria included the absolute value of difference in gene expression on a natural log scale of at least 0.25 and adjusted P value < .05. Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes enrichment analyses were performed using clusterProfiler621 version 3.16.1.

Survival analysis

Survival analysis was performed using the survminer package in R. Gene expression and survival data from the The Cancer Genome Atlas and MMRF CoMMpass data sets were downloaded from the Genomic Data Commons platform. The MMRF CoMMpass data set contains 859 RNA sequences from CD138+ PCs in BM samples from patients with MM. Data on overall survival or progression-free survival were available for 787 patients. The samples were divided into 2 groups (high and low gene expression) based on median values. Survival was assessed using the Kaplan-Meier method and compared using the 2-sided log-rank test.

SCENIC

Single-cell regulatory network analysis was performed using pySCENIC22 version 0.11.0 in Python, following the protocol described in the single-cell regulatory network inference and clustering (SCENIC) workflow. Coexpression gene regulatory networks were generated using the “pyscenic grn” command and the GRNBoost2 method. AUCell analysis was performed using the “pyscenic aucell” function with the following parameters: rank_threshold = 5000, auc_threshold = 0.05, and nes_threshold = 3.

The generation of LILRB4 knockout (KO) NCI-H929 cells and LILRB4-overexpressing (OE) NCI-H929 cells and procedures of cell culture, cell proliferation and viability analysis, western blot, real-time quantitative polymerase chain reaction, cell cycle analysis, transendothelial migration assay, tube formation in human umbilical vein endothelial cells (HUVECs), and flow cytometry can be found in the supplemental Materials.

Single-cell RNA-seq mapping of MM cells in pleural effusion and ascites

We evaluated 3 patients with MM who relapsed with ME. The disease course is summarized in Figure 1A, and patient characteristics are shown in Table 1. Malignant PCs in MEs were detected using flow cytometry. BM, PB, and ME were collected from each patient and analyzed using scRNA-seq. After quality control, the single-cell transcriptomes of 37 731 cells were selected (supplemental Figure 1A). In addition, 17 599 malignant and 20 132 nonmalignant cells were distinguished using 3 complementary approaches. After the integration, single-cell transcriptomes were visualized in a two-dimensional map using the UMAP method for unsupervised dimensionality reduction (Figure 1B). The BM of all patients contained PCs, immune cells, myeloid cells, and erythroblasts (Figure 1C), whereas the pleural and peritoneal effusion contained abundant PCs and few immune components. The heatmap of the top 10 DEGs for each cell type is shown in supplemental Figure 1B.

Figure 1.

Single-cell atlas of the immune microenvironment in MM. (A) Disease course in 3 patients with MM with the extramedullary disease. (B) Analysis of cell clusters using UMAP, showing the annotation and color codes for cell types in extramedullary sites; (left) 37 731 cells color-coded based on cell type; (right) 37 731 cells color-coded by sample. (C) Proportion of cell types in the BM and PB of 3 patients with MM. (D) UMAP plot of immune cell subsets in the myeloma microenvironment. Red: regulatory T cells; green: γδ T cells; orange: CD8+ effector T cells; and blue: naive T cells. (E) T-cell subtypes in the BM, ME, and PB of 3 patients with MM. (F) DEGs in 4 T-cell subtypes. ASCT, autologous stem cell transplantation; NMPs, neuromesodermal progenitors; PMCs, pleural mesothelial cells.

Figure 1.

Single-cell atlas of the immune microenvironment in MM. (A) Disease course in 3 patients with MM with the extramedullary disease. (B) Analysis of cell clusters using UMAP, showing the annotation and color codes for cell types in extramedullary sites; (left) 37 731 cells color-coded based on cell type; (right) 37 731 cells color-coded by sample. (C) Proportion of cell types in the BM and PB of 3 patients with MM. (D) UMAP plot of immune cell subsets in the myeloma microenvironment. Red: regulatory T cells; green: γδ T cells; orange: CD8+ effector T cells; and blue: naive T cells. (E) T-cell subtypes in the BM, ME, and PB of 3 patients with MM. (F) DEGs in 4 T-cell subtypes. ASCT, autologous stem cell transplantation; NMPs, neuromesodermal progenitors; PMCs, pleural mesothelial cells.

Close modal

Immune microenvironment atlas accompanying myeloma migration

Then, we investigated the immune microenvironment within or outside the BM. The predominant immune cell types were T cells (Figure 1D), indicating humoral immunodeficiency in patients with myeloma. The predominant T-cell cluster (with ∼3546 cells) contained the following 4 cell subsets: naive T cells, CD8+ effector T cells, γδ T cells, and regulatory T cells. MEs contained a higher percentage of cytotoxic T cells, whereas PB contained a higher proportion of naive T cells (Figure 1D; supplemental Figure 2A). Of note, there was a population of Treg cells in the effusion group, indicating that reduced immune response might contribute to the formation of ME (Figure 1E).

We analyzed different components of T cells in PB and ME. We found that T cells in the ME highly expressed proliferation protein (eg, MKI67), proteins involved in the regulation of the microtubule filament system (eg, STMN1), and mobility protein (eg, HMGB2), indicating that T cells in the ME were proliferating (supplemental Figure 2B). CD8+ effector T cells differed in the expression of cell cycle proteins (eg, TOP2A) and were associated with cell proliferation (Figure 1F; supplemental Figure 2C). We observed higher cytotoxic scores of the T cells in the BM than those in extramedullary lesions, whereas there were no statistically significant differences in exhaustion scores between that within and outside the BM (supplemental Figure 2D). We observed the same phenomenon in CD8+ effector T cells (supplemental Figure 2E).

Transcriptional landscape of malignant PCs

Next, we determined the scRNA-seq profile of the BM from 3 patients with MM for matching extramedullary involvement from these patients (ME and PB). scRNA-seq was not performed in 1 PB sample because the number of cells was small. A total of 17 599 malignant cells clustered based on the tumor of origin (Figure 2A). We removed immunoglobulin genes from the clustering analysis as the high expression of clones limited the assessment of tumor heterogeneity. Different clusters of malignant PCs were identified in the BM, ME, and PB (Figure 2B). Some clusters mainly derived from single samples, whereas others included samples from both the extramedullary site and BM. Every cluster contained unique signature genes, suggesting the intertumoral heterogeneity of PCs. DEGs in each cluster are shown as a heatmap (Figure 2C). There was considerable variability between the patients and their sites (BM vs ME). Clusters 1, 4, and 6 contained PCs originated from ME samples, expressing cellular stress response protein (SGK1; cluster 1), MT1X (cluster 4), and GNB2L1 (cluster 6) gene that played a role in myeloma cell survival, migration, and contact with the extracellular matrix. The myeloma cells in BM samples were clusters 2, 5, and 7 expressing the regulator of G-protein signaling family RGS1 (cluster 2), cyclin family gene CCND2 (cluster 5), and the inflammatory responses gene ANXA1 (cluster 7) that regulated the cell cycle and promoted the growth of myeloma cells in the BM. Although PC clusters were site- and patient-specific, myeloma cells in the BM and extramedullary sites shared the same subclusters, indicating the presence of a transition state of extramedullary migration. MM01 presented PC clusters 3, 8, 10, and 11, whereas MM03 presented PC clusters 9 and 11, which existed both inside and outside the BM; these clusters contained genes encoding cell cycle regulation proteins (eg, CCND2), proliferation marker proteins (eg, MKI67), chemokine receptors (eg, CXCR4), extracellular matrix proteins (eg, MMP9), and high-mobility group proteins (eg, HMGB2), respectively. Clusters 3, 8, and 10 contained genes encoding major histocompatibility complex (MHC) class II molecules, such as HLA-DRA and HLA-DRB1 (supplemental Table 1).

Figure 2.

Single-cell atlas of malignant plasma cells (MPCs) from patients with MM. (A) Analysis of 17 599 MPCs color-coded using UMAP based on the cluster. (B) Proportion of MPCs in the BM and PB of 3 patients with MM. (C) Heatmaps of DEGs for each cluster in MPCs. (D) DEGs in MPCs from BM and extramedullary sites (ME + PB).

Figure 2.

Single-cell atlas of malignant plasma cells (MPCs) from patients with MM. (A) Analysis of 17 599 MPCs color-coded using UMAP based on the cluster. (B) Proportion of MPCs in the BM and PB of 3 patients with MM. (C) Heatmaps of DEGs for each cluster in MPCs. (D) DEGs in MPCs from BM and extramedullary sites (ME + PB).

Close modal

Next, we performed differential gene expression analysis of PCs from the BM and extramedullary site. The expression of the metallothionein family (eg, MT1X and MT2A) and SGK1 increased in extramedullary PCs, whereas gene encoding S100 family of proteins (eg, S100A8 and S100A9) and hemoglobin (eg, HBA1, HBA2, and HBB) were upregulated in PCs of the BM (Figure 2D).

Signaling pathways enriched in MM cells

Signaling pathways enriched in malignant PCs from the BM and extramedullary sites were analyzed. Extramedullary disease upregulated pathways associated with cellular response to oxidative stress, protein kinase regulator activity, and reactive oxygen species metabolic process, and downregulated pathways associated with protein folding and endoplasmic reticulum (ER) response (supplemental Table 2).

The malignant PCs of MM01 were heterogeneous, which made up of clusters 3, 5, 7, 8, 10, and 11. Kyoto Encyclopedia of Genes and Genomes analysis showed that the actin cytoskeleton pathway and leukocyte transendothelial migration were upregulated in cluster 11 (supplemental Figure 3A). GO enrichment analysis indicated that the focal adhesion and cell-substrate conjunction were downregulated in the cellular component of cluster 10 (supplemental Figure 3B). These 2 clusters (10 and 11) were mainly present in the BM and PB but rarely in the ME. Cluster 8 not only existed in the BM but also occupied the majority portion in the PB and ME. Proteasome and antigen presentation signal pathway were upregulated in this cluster (supplemental Figure 3C).

MM02 presented PC clusters 3, 4, 6, and 7. Clusters 3 (13 cells in the ME of MM02) and 7 were observed in both the BM and ME. GO analysis indicated that cell cycle control and oxidative phosphorylation pathways were enriched in cluster 3, and ER and unfolded protein response (UPR) pathways were enriched in cluster 7 (supplemental Figure 4A-B). Clusters 4 and 6 were found in the ME and were associated with the upregulation of reactive oxygen species pathways (supplemental Figure 4C-D).

MM03 showed PC clusters 1, 2, 8, 9, and 11. Clusters 1 (8 cells in the BM of MM03) and 2 (2 cells in the ME of MM03) were found in the BM and ME and were associated with the upregulation of ER pathways (supplemental Figure 5A-B). In addition, cell cycle control and oxidative phosphorylation pathways were upregulated in cluster 9 (supplemental Figure 5C).

Extramedullary-initiating cells (EMICs) in MM

We performed hotspot analysis and found that the gene modules 3 and 7 were highly expressed in intramedullary and extramedullary PC clusters, suggesting that these modules were associated with cell viability and migration (Figure 3A). The gene sets in these modules were present in PCs both from extramedullary lesions and BM (Figure 3B). We observed a higher expression of modules 3 and 7 in the PCs of the BM compared to PB and ME, according to UCell Gene Set Scoring23 (supplemental Figure 6A). We speculated that the gene set was a key factor in the survival and proliferation of myeloma cells migrating into the pleural space, and we defined cell clusters expressing this gene set as EMICs. Hub genes in module 3 (TOP2A, NUSAP1, and TPX2) correlated with nuclear division, whereas hub genes in module 7 (MKI67 and HMGB1) were implicated in cell proliferation (supplemental Table 3). We further subdivided the functions of the genes and found that the functions of most of them were related to the cell cycle, P53 signal pathway, and other cellular functions (Figure 3C), suggesting that the improvement in these functions could help tumor cells adapt to the extramedullary environment. We subsequently explored the expression of modules 3 and 7 in GSE66291 and GSE39683 data sets. Single sample gene set enrichment analysis was performed on 8 newly diagnosed multiple myeloma (NDMM) cases and 8 secondary plasma cell leukemia (sPCL) cases from GSE66291 and 7 NDMM cases and 7 sPCL cases from GSE39683. We observed higher enrichment scores of modules 3 and 7 in the sPCL cases, indicating that the EMICs played an important role in the extramedullary migration of PCs (supplemental Figure 6B).

Figure 3.

Gene modules 3 and 7 play a role in myeloma cell migration and viability. (A) z scores of 10 gene modules highly expressed in PCs from the BM and extramedullary sites. (B) Scores for gene signatures 3 and 7 in MPC clusters. (C) Pathways significantly enriched in EMICs. The cut-off criterion was a false-discovery rate <0.01. (D) Overall survival of patients with different expression of gene sets 3 and 7. Statistical significance was calculated using the log-rank test.

Figure 3.

Gene modules 3 and 7 play a role in myeloma cell migration and viability. (A) z scores of 10 gene modules highly expressed in PCs from the BM and extramedullary sites. (B) Scores for gene signatures 3 and 7 in MPC clusters. (C) Pathways significantly enriched in EMICs. The cut-off criterion was a false-discovery rate <0.01. (D) Overall survival of patients with different expression of gene sets 3 and 7. Statistical significance was calculated using the log-rank test.

Close modal

There were 68 genes in modules 3 and 7. Overlapped with 50 genes of patients with relapsed/refractory multiple myeloma from the DisGeNET database, we suggested that ANP32B and PCNA might be the biomarkers of EMICs (supplemental Figure 6C).

We analyzed the relationship between gene expression and prognosis in 2 MM cohorts from the MMRF CoMMpass data set (n = 787). The high expression of genes in these modules was associated with worse overall survival in patients with MM (Figure 3D), whereas the other modules had no significant effect on prognosis (supplemental Figure 7A).

Then, we considered how EMICs relate to genetic heterogeneity and tumor extramedullary evolution. We assessed genetic differences between EMICs and other tumor cells. EMICs expressed microtubule-associated genes (TUBA1B and NUSAP1), proliferative genes (TTK and BIRC5), and cell cycle–associated genes (CCNB1 and FAM111A) (supplemental Table 4). In addition, we identified the molecular pathways in which the differential genes were expressed. The cell cycle pathway was upregulated, whereas protein processing in ER pathway was downregulated in EMICs (supplemental Figure 7B). Genetic changes in extramedullary clones and intramedullary clones indicated they are substantially heterogeneous.

Furthermore, considering that both intramedullary and extramedullary cells possess EMICs, we further examined the genetic differences between them. The expression of leukocyte immunoglobulin–like (LIL) family proteins LILRB4 and long noncoding RNAs (MALAT1 and NEAT1) were upregulated in extramedullary EMICs (supplemental Table 5). Intramedullary EMICs had the potential to migrate extramedullary after acquiring some extra genetic alterations.

Finally, we investigated the different transcription-factor expression levels between EMICs and other tumor cells. Cells that highly expressed gene signatures of modules 3 and 7 exhibited specific expression of transcription factors essential for the survival and growth of MM tumors such as E2F1, YY1, and SMAD1 (supplemental Figure 7C). These 3 genes were associated with transforming growth factor β and NF-κB signaling pathways, suggesting that the NF-κB pathway played a crucial role in myeloma migration to extramedullary sites. Given that transcription factors were important in defining cell identity, we then compared the expression levels of transcription factors of intramedullary and extramedullary EMICs and found high activities of RUNX2 and antigen presentation (MHC class II) pathway in extramedullary EMICs. Upregulated transcription factors, such as PRDM1, STAT1, and IRF family, related to the master regulator of MHC class II gene expression–MHC Class II transactivator, might reduce T-cell–mediated cytotoxicity (supplemental Figure 8A). These results suggested that EMICs were involved in tumor cell proliferation and migration and immune evasion.

LILRB4 is upregulated in ME and reduces patient survival

Gene signatures associated with myeloma cell migration were analyzed. The genes LILRB4, H3F3B, ARPC2, TMEM173, LAP3, and CD63 were upregulated, whereas FOS, SLC1A5, BCMA, DUSP1, P4HB, and JUN were downregulated at extramedullary sites (Figure 4A).

Figure 4.

Effect of LILRB4 on survival, PARP and caspase-3 cleavage expression, and cell cycle progression. (A) Venn diagrams of DEGs of PCs within or outside the BM from 3 patients. (B) Expression of LILRB4 in PCs from BM and ME. (C) Overall survival of patients with different expressions of LILRB4. (D-E) Effect of LILRB4 overexpression and KO on the protein and messenger RNA expression of LILRB4 in NCI-H929 cells. (F) Effect of LILRB4 overexpression and KO on PARP and caspase 3 cleavage in NCI-H929 cells. (G-H) Effect of LILRB4 overexpression and KO on cell cycle progression.

Figure 4.

Effect of LILRB4 on survival, PARP and caspase-3 cleavage expression, and cell cycle progression. (A) Venn diagrams of DEGs of PCs within or outside the BM from 3 patients. (B) Expression of LILRB4 in PCs from BM and ME. (C) Overall survival of patients with different expressions of LILRB4. (D-E) Effect of LILRB4 overexpression and KO on the protein and messenger RNA expression of LILRB4 in NCI-H929 cells. (F) Effect of LILRB4 overexpression and KO on PARP and caspase 3 cleavage in NCI-H929 cells. (G-H) Effect of LILRB4 overexpression and KO on cell cycle progression.

Close modal

LILRB4, an immune inhibitory receptor that contains putative immunoreceptor tyrosine–based inhibitory motifs, was more highly expressed in the ME than in the BM (Wilcoxon, P < 2.2e-16; Figure 4B). The analysis of the MMRF CoMMpass data set showed that the high expression of LILRB4 or ARPC2 was associated with poor survival, whereas H3F3B, TMEM173, and LAP3 had no significant effect on prognosis (Figure 4C; supplemental Figure 8B). Furthermore, the LILRB4 gene was significantly upregulated in extramedullary EMICs (supplemental Figure 8C).

LILRB4 promotes myeloma cell proliferation and migration in vitro

To verify the role of LILRB4 on ME development, we used the CRISPR-Cas9 strategy to knock out LILRB4 (LILRB4-KO), and overexpress LILRB4 (LILRB4-OE) via lentiviral transfection in MM cell line (NCI-H929). Western blot (Figure 4D) and real-time polymerase chain reaction (Figure 4E) were performed to verify the overexpression or KO of LILRB4. Then, we detected the level of apoptotic marker cleavage of PARP and caspase 3 in the LILRB4-OE and LILRB4-KO cell lines. The cleaved PARP and caspase 3 were upregulated in LILRB4-KO cells and downregulated in LILRB4-OE cells (Figure 4F).

The effect of LILRB4 on cell cycle progression was also assessed. LILRB4-OE induced arrest in the G1 phase, whereas the percentage of cells in the G1 phase was strikingly decreased in LILRB4-KO cells (Figure 4G-H). There was no significant difference in cell proliferation between LILRB4-OE/KO and control H929 cells using the CCK-8 assay, indicating that LILRB4 did not induce cell proliferation (Figure 5A-B). The effect of LILRB4 on myeloma cell migration was evaluated using a transwell chamber (Figure 5C) and tube formation in HUVECs. LILRB4 overexpression increased cell migration (Figure 5D-E) and promoted the tubular structure of HUVECs (supplemental Figure 9A-B). Meanwhile, we observed decreased transendothelial migration (Figure 5F) and disrupted the tubular structure of HUVECs in LILRB4-KO (supplemental Figure 9C). Moreover, LILRB4 changed the messenger RNA expression of angiopoietin 1, integrin αv, and integrin β3 (Figure 5G-I).

Figure 5.

LILRB4 contributes to MM cell migration. (A-B) Analysis of the effect of LILRB4 on cell proliferation using the cell counting kit-8 assay. (C) Cell migration assay. HUVECs expressing green fluorescence protein were seeded on the transwell membrane (8 μm pore size). After 48 hours, LILRB4-OE and KO NCI-H929 cells expressing red fluorescence protein were cultured in the upper chamber for 24 hours. (D) Number of representative migrating cells counted using confocal microscopy. Scale bar, 50 μm. (E-F) The number of migrating LILRB4-OE and KO NCI-H929 cells. Data are means and standard error of the mean of 3 independent experiments. (G-I) Effect of LILRB4 KO and overexpression on the messenger RNA expression of angiopoietin 1, integrin αv, and integrin β3.

Figure 5.

LILRB4 contributes to MM cell migration. (A-B) Analysis of the effect of LILRB4 on cell proliferation using the cell counting kit-8 assay. (C) Cell migration assay. HUVECs expressing green fluorescence protein were seeded on the transwell membrane (8 μm pore size). After 48 hours, LILRB4-OE and KO NCI-H929 cells expressing red fluorescence protein were cultured in the upper chamber for 24 hours. (D) Number of representative migrating cells counted using confocal microscopy. Scale bar, 50 μm. (E-F) The number of migrating LILRB4-OE and KO NCI-H929 cells. Data are means and standard error of the mean of 3 independent experiments. (G-I) Effect of LILRB4 KO and overexpression on the messenger RNA expression of angiopoietin 1, integrin αv, and integrin β3.

Close modal

EMM is a malignancy that presents as a soft tissue mass of monoclonal PCs. Pleural and peritoneal involvement, including ME, is rare, accounting for less than 1% of MM cases.24 ME was assessed in case reports25,26 and a systematic review.12 Nonetheless, the molecular mechanisms underlying this entity are incompletely understood. The use of scRNA-seq analysis allows for determining the transcriptional profile of immune cells in the BM and extramedullary sites.

RNA-seq analysis of the immune microenvironment showed that the genes implicated in cell proliferation, migration, and cell cycle progression were highly expressed in effector T cells from MEs. The expansion of effector T-cell subtypes does not produce a cytotoxic effect on tumor.27 Genes associated with T-cell exhaustion, including PD-1, CTLA-4, 2B4, and CD160, were expressed in T cells from MM,28 especially in extramedullary sites,15 indicating that the dysregulation of T-cell immune surveillance might play a role in ME formation.

EMM has a more aggressive clinical course than myeloma.29 ScRNA-seq helps to assess the evolution of malignant PCs.30 We found that clusters 3, 8, 9, 10, and 11, which were both present at BM and extramedullary sites, were associated with PC migration and shared proliferation marker genes (eg, MKI67) and cell migration–related genes (CXCR31,32 and HMGB33,34 family), which were related to poor prognosis in MM. Consistent with a previous study,15 these clusters were associated with the proteasome, cell cycle control, antigen presentation, and oxidative phosphorylation pathways.

Cancer stem–like cells (CSCs) contribute to cancer progression and metastasis.35-37 However, CSCs in MM have not been thoroughly defined.38 ME often occurs during relapse, and disease progression and metastasis are driven by clonal selection.39 We identified a PC subtype in the BM, PB, and ME through hotspot analysis, and this subtype survived in extramedullary sites. In addition, this subtype expressed genes related to cell cycle control, chromatin remodeling, microtubule dynamics, and other cellular functions. Cell cycle dysregulation40 and higher proliferation indices41 contribute to drug resistance and poor outcomes. In addition, we determined that the expression of modules 3 and 7 was higher in the BM of sPCL than that of NDMM using 2 other databases (GSE39683 and GSE66291). On the basis of data from a recent study, there was significant tumor heterogeneity between BM and extramedullary sites in MM.3 We speculated that EMICs, similar to CSCs, were a high-risk PC clone and more difficult to remove.42 It was necessary to identify biological markers for predicting the potential of developing ME. Two potential genes were filtered out by overlapping the signature of relapsed/refractory multiple myeloma and EMICs. ANP32B, a potent survival factor for exogenous cytokine–dependent myeloma cells, upregulated the expression of the 2 major antiapoptotic proteins in myeloma cells, Mcl-1 and Bcl-2.43 PCNA was an essential factor for DNA replication and repair. Its expression on MM cells, interacting with the NKp44 receptor, induced immuno-compromise by inhibiting NK cell function.44 

The transcription factors E2F1, SMAD1, and YY1 were upregulated in EMICs. E2F1, a gene encoding DNA-binding protein, was indirectly activated by cyclin D1 and induced the growth of MM.45,SMAD1 was a determinant of drug resistance in MM cells and an unfavorable prognostic marker,46 whereas a higher nuclear YY1 expression was associated with disease progression.47 

EMICs promoted tumor cell survival in extramedullary sites, indicating that EMICs were the prerequisite for extramedullary metastasis. Thus, we speculated that based on the presence of EMICs, further genetic alterations contributed to extramedullary infiltration. A set of genes was consistently expressed in PCs in the BM and extramedullary sites. LILRB4, an inhibitory receptor of immune checkpoint pathways,48 was upregulated in ME and was known to suppress the immune response in solid tumors49 and promote the infiltration of myeloid leukemia cells.50 Furthermore, LILRB4 induced nonsmall lung cancer invasion by increasing epithelial-mesenchymal transition and the expression of vascular endothelial growth factor A.51 The KO of LILRB4 in MM cells decreased cell migration, and LILRB4 overexpression had the opposite effect. These results suggested that LILRB4 was involved in myeloma extramedullary migration. Besides this, we also detected the association between high LILRB4 expression and poor outcomes. Of interest, we found that BCMA, an immunotherapy target for MM,52 was downregulated in extramedullary sites. The efficacy of anti-BCMA chimeric antigen receptor T-cell therapy for EMM is limited.53,54 

In summary, the transcriptome analysis of PCs and the immune microenvironment at different physiological sites revealed that patients with EMM had a unique gene expression profile, including gene sets associated with cell migration and survival.

This study has limitations. Firstly, the number of ME samples was small because the incidence of this entity was low. Notwithstanding, this study helped elucidate the mechanisms underlying MM cell proliferation and invasion. Through the analysis of intratumor heterogeneity, we believe that the development of ME needed a subpopulation of proliferating stem–like cells, whereas invasion-related genes such as LILRB4 mediated the migration of PCs.

The authors thank Xiaming Kong for helping with data analysis of single-cell sequencing.

This study was supported by National Natural Science Foundation of China (81720108002), CSCO Research Foundation (Y-Roche2019/2-0090), Science Foundation Project of Ili & Jiangsu Joint Institute of Health (yl2021ms04), and Jiangsu Province Hospital (the First Affiliated Hospital with Nanjing Medical University) Clinical Capacity Enhancement Project (JSPH-MB-2021-14 and JSPH-MB-2021-9).

Contribution: X.Q., J.L., Z.S., and J.J. contributed to the conception and design; Y.L., Y.C., and F.L. contributed to the acquisition of data; X.Q., Z.S., and J.J. contributed to the analysis and interpretation of data; X.Q., Z.S., and J.J. contributed to the writing, review, and/or revision of the manuscript; X.Q. and J.L. supervised the study; and all authors read and approved the final manuscript.

Conflict-of-interest disclosure: The authors declare no competing financial interests.

Correspondence: Jianyong Li, Department of Hematology, The First Affiliated Hospital of Nanjing Medical University, Jiangsu Province Hospital, 300 Guangzhou Rd, Nanjing, 210029, China; e-mail: lijianyonglm@126.com; and Xiaoyan Qu, Department of Hematology, The First Affiliated Hospital of Nanjing Medical University, Jiangsu Province Hospital, 300 Guangzhou Rd, Nanjing, 210029, China; e-mail: quxiaoyan@jsph.org.cn.

1.
van de Donk
N
,
Pawlyn
C
,
Yong
KL
.
Multiple myeloma
.
Lancet
.
2021
;
397
(
10272
):
410
-
427
.
2.
Silberstein
J
,
Tuchman
S
,
Grant
SJ
.
What is multiple myeloma?
.
JAMA
.
2022
;
327
(
5
):
497
.
3.
Shen
YJ
,
Mishima
Y
,
Shi
J
, et al
.
Progression signature underlies clonal evolution and dissemination of multiple myeloma
.
Blood
.
2021
;
137
(
17
):
2360
-
2372
.
4.
Varettoni
M
,
Corso
A
,
Pica
G
,
Mangiacavalli
S
,
Pascutto
C
,
Lazzarino
M
.
Incidence, presenting features and outcome of extramedullary disease in multiple myeloma: a longitudinal study on 1003 consecutive patients
.
Ann Oncol
.
2010
;
21
(
2
):
325
-
330
.
5.
Jagosky
MH
,
Usmani
SZ
.
Extramedullary disease in multiple myeloma
.
Curr Hematol Malig Rep
.
2020
;
15
(
2
):
62
-
71
.
6.
Bhutani
M
,
Foureau
DM
,
Atrash
S
,
Voorhees
PM
,
Usmani
SZ
.
Extramedullary multiple myeloma
.
Leukemia
.
2020
;
34
(
1
):
1
-
20
.
7.
Sevcikova
S
,
Minarik
J
,
Stork
M
,
Jelinek
T
,
Pour
L
,
Hajek
R
.
Extramedullary disease in multiple myeloma - controversies and future directions
.
Blood Rev
.
2019
;
36
:
32
-
39
.
8.
Byun
JM
,
Kim
KH
,
Choi
IS
, et al
.
Pleural effusion in multiple myeloma: characteristics and practice patterns
.
Acta Haematol
.
2017
;
138
(
2
):
69
-
76
.
9.
Singh
D
,
Kumar
L
.
Myelomatous ascites in multiple myeloma
.
Leuk Lymphoma
.
2005
;
46
(
4
):
631
-
632
.
10.
Sekiguchi
Y
,
Shirane
S
,
Imai
H
, et al
.
Response to low-dose bortezomib in plasma cell leukemia patients with malignant pleural effusion and ascites: a case report and a review of the literature
.
Intern Med
.
2012
;
51
(
11
):
1393
-
1398
.
11.
Yanamandra
U
,
Deo
P
,
Sahu
KK
, et al
.
Clinicopathological profile of myelomatous pleural effusion: single-center real-world experience and review of literature
.
Clin Lymphoma Myeloma Leuk
.
2019
;
19
(
3
):
183
-
189.e1
.
12.
Riveiro
V
,
Ferreiro
L
,
Toubes
ME
,
Lama
A
,
Álvarez-Dobaño
JM
,
Valdés
L
.
Characteristics of patients with myelomatous pleural effusion. A systematic review
.
Rev Clin Esp
.
2018
;
218
(
2
):
89
-
97
.
13.
Wang
Z
,
Xia
G
,
Lan
L
, et al
.
Pleural effusion in multiple myeloma
.
Intern Med
.
2016
;
55
(
4
):
339
-
345
.
14.
Rasche
L
,
Chavan
SS
,
Stephens
OW
, et al
.
Spatial genomic heterogeneity in multiple myeloma revealed by multi-region sequencing
.
Nat Commun
.
2017
;
8
(
1
):
268
.
15.
Ryu
D
,
Kim
SJ
,
Hong
Y
, et al
.
Alterations in the transcriptional programs of myeloma cells and the microenvironment during extramedullary progression affect proliferation and immune evasion
.
Clin Cancer Res
.
2020
;
26
(
4
):
935
-
944
.
16.
Rodríguez
JN
,
Pereira
A
,
Martínez
JC
,
Conde
J
,
Pujol
E
.
Pleural effusion in multiple myeloma
.
Chest
.
1994
;
105
(
2
):
622
-
624
.
17.
Dura
B
,
Choi
JY
,
Zhang
K
, et al
.
scFTD-seq: freeze-thaw lysis based, portable approach toward highly distributed single-cell 3' mRNA profiling
.
Nucleic Acids Res
.
2019
;
47
(
3
):
e16
.
18.
Dobin
A
,
Davis
CA
,
Schlesinger
F
, et al
.
STAR: ultrafast universal RNA-seq aligner
.
Bioinformatics
.
2013
;
29
(
1
):
15
-
21
.
19.
Liao
Y
,
Smyth
GK
,
Shi
W
.
featureCounts: an efficient general purpose program for assigning sequence reads to genomic features
.
Bioinformatics
.
2014
;
30
(
7
):
923
-
930
.
20.
Satija
R
,
Farrell
JA
,
Gennert
D
,
Schier
AF
,
Regev
A
.
Spatial reconstruction of single-cell gene expression data
.
Nat Biotechnol
.
2015
;
33
(
5
):
495
-
502
.
21.
Yu
G
,
Wang
LG
,
Han
Y
,
He
QY
.
clusterProfiler: an R package for comparing biological themes among gene clusters
.
OMICS
.
2012
;
16
(
5
):
284
-
287
.
22.
Van de Sande
B
,
Flerin
C
,
Davie
K
, et al
.
A scalable SCENIC workflow for single-cell gene regulatory network analysis
.
Nat Protoc
.
2020
;
15
(
7
):
2247
-
2276
.
23.
Andreatta
M
,
Carmona
SJ
.
UCell: robust and scalable single-cell gene signature scoring
.
Comput Struct Biotechnol J
.
2021
;
19
:
3796
-
3798
.
24.
Harbhajanka
A
,
Brickman
A
,
Park
JW
,
Reddy
VB
,
Bitterman
P
,
Gattuso
P
.
Cytomorphology, clinicopathologic, and cytogenetics correlation of myelomatous effusion of serous cavities: a retrospective review
.
Diagn Cytopathol
.
2016
;
44
(
9
):
742
-
747
.
25.
Jha
LK
,
Dingli
D
.
An unusual cause of ascites
.
Blood
.
2015
;
126
(
13
):
1629
.
26.
Martinez-Iribarren
A
,
Fernandez-Prendes
C
.
A case of myelomatous pleural effusion: an unusual onset of multiple myeloma
.
Blood
.
2020
;
135
(
2
):
153
.
27.
Zavidij
O
,
Haradhvala
NJ
,
Mouhieddine
TH
, et al
.
Single-cell RNA sequencing reveals compromised immune microenvironment in precursor stages of multiple myeloma
.
Nat Cancer
.
2020
;
1
(
5
):
493
-
506
.
28.
Zelle-Rieser
C
,
Thangavadivel
S
,
Biedermann
R
, et al
.
T cells in multiple myeloma display features of exhaustion and senescence at the tumor site
.
J Hematol Oncol
.
2016
;
9
(
1
):
116
.
29.
Beksac
M
,
Seval
GC
,
Kanellias
N
, et al
.
A real world multicenter retrospective study on extramedullary disease from Balkan Myeloma Study Group and Barcelona University: analysis of parameters that improve outcome
.
Haematologica
.
2020
;
105
(
1
):
201
-
208
.
30.
Liu
R
,
Gao
Q
,
Foltz
SM
, et al
.
Co-evolution of tumor and immune cells during progression of multiple myeloma
.
Nat Commun
.
2021
;
12
(
1
):
2559
.
31.
Ito
S
,
Sato
T
,
Maeta
T
.
Role and therapeutic targeting of SDF-1α/CXCR4 axis in multiple myeloma
.
Cancers
.
2021
;
13
(
8
):
1793
.
32.
Ullah
TR
.
The role of CXCR4 in multiple myeloma: cells' journey from bone marrow to beyond
.
J Bone Oncol
.
2019
;
17
:
100253
.
33.
Guo
X
,
He
D
,
Zhang
E
, et al
.
HMGB1 knockdown increases MM cell vulnerability by regulating autophagy and DNA damage repair
.
J Exp Clin Cancer Res
.
2018
;
37
(
1
):
205
.
34.
Gao
D
,
Lv
AE
,
Li
HP
,
Han
DH
,
Zhang
YP
.
LncRNA MALAT-1 elevates HMGB1 to promote autophagy resulting in inhibition of tumor cell apoptosis in multiple myeloma
.
J Cell Biochem
.
2017
;
118
(
10
):
3341
-
3348
.
35.
Lawson
DA
,
Kessenbrock
K
,
Davis
RT
,
Pervolarakis
N
,
Werb
Z
.
Tumour heterogeneity and metastasis at single-cell resolution
.
Nat Cell Biol
.
2018
;
20
(
12
):
1349
-
1360
.
36.
Dagogo-Jack
I
,
Shaw
AT
.
Tumour heterogeneity and resistance to cancer therapies
.
Nat Rev Clin Oncol
.
2018
;
15
(
2
):
81
-
94
.
37.
Huang
T
,
Song
X
,
Xu
D
, et al
.
Stem cell programs in cancer initiation, progression, and therapy resistance
.
Theranostics
.
2020
;
10
(
19
):
8721
-
8743
.
38.
Yaccoby
S
.
Two states of myeloma stem cells
.
Clin Lymphoma Myeloma Leuk
.
2018
;
18
(
1
):
38
-
43
.
39.
Abe
M
,
Harada
T
,
Matsumoto
T
.
Concise review: defining and targeting myeloma stem cell-like cells
.
Stem Cells
.
2014
;
32
(
5
):
1067
-
1073
.
40.
Salomon-Perzyński
A
,
Jamroziak
K
,
Głodkowska-Mrówka
E
.
Clonal evolution of multiple myeloma-clinical and diagnostic implications
.
Diagnostics
.
2021
;
11
(
9
):
1534
.
41.
Bergsagel
PL
,
Kuehl
WM
,
Zhan
F
,
Sawyer
J
,
Barlogie
B
,
Shaughnessy
J
.
Cyclin D dysregulation: an early and unifying pathogenic event in multiple myeloma
.
Blood
.
2005
;
106
(
1
):
296
-
303
.
42.
Pawlyn
C
,
Morgan
GJ
.
Evolutionary biology of high-risk multiple myeloma
.
Nat Rev Cancer
.
2017
;
17
(
9
):
543
-
556
.
43.
Moreaux
J
,
Legouffe
E
,
Jourdan
E
, et al
.
BAFF and APRIL protect myeloma cells from apoptosis induced by interleukin 6 deprivation and dexamethasone
.
Blood
.
2004
;
103
(
8
):
3148
-
3157
.
44.
Iraqi
M
,
Edri
A
,
Greenshpan
Y
, et al
.
Blocking the PCNA/NKp44 checkpoint to stimulate NK cell responses to multiple myeloma
.
Int J Mol Sci
.
2022
;
23
(
9
):
4717
.
45.
Liu
JL
,
Zeng
GZ
,
Liu
XL
, et al
.
Small compound bigelovin exerts inhibitory effects and triggers proteolysis of E2F1 in multiple myeloma cells
.
Cancer Sci
.
2013
;
104
(
12
):
1697
-
1704
.
46.
Wu
J
,
Zhang
M
,
Faruq
O
, et al
.
SMAD1 as a biomarker and potential therapeutic target in drug-resistant multiple myeloma
.
Biomark Res
.
2021
;
9
(
1
):
48
.
47.
Che
F
,
Ye
X
,
Wang
Y
, et al
.
METTL3 facilitates multiple myeloma tumorigenesis by enhancing YY1 stability and pri-microRNA-27 maturation in m(6)A-dependent manner
.
Cell Biol Toxicol
.
2022
.
48.
Liu
J
,
Wu
Q
,
Shi
J
, et al
.
LILRB4, from the immune system to the disease target
.
Am J Transl Res
.
2020
;
12
(
7
):
3149
-
3166
.
49.
Sharma
N
,
Atolagbe
OT
,
Ge
Z
,
Allison
JP
.
LILRB4 suppresses immunity in solid tumors and is a potential target for immunotherapy
.
J Exp Med
.
2021
;
218
(
7
):
e20201811
.
50.
Deng
M
,
Gui
X
,
Kim
J
, et al
.
LILRB4 signalling in leukaemia cells mediates T cell suppression and tumour infiltration
.
Nature
.
2018
;
562
(
7728
):
605
-
609
.
51.
Li
Z
,
Deng
M
,
Huang
F
, et al
.
LILRB4 ITIMs mediate the T cell suppression and infiltration of acute myeloid leukemia cells
.
Cell Mol Immunol
.
2020
;
17
(
3
):
272
-
282
.
52.
Yu
B
,
Jiang
T
,
Liu
D
.
BCMA-targeted immunotherapy for multiple myeloma
.
J Hematol Oncol
.
2020
;
13
(
1
):
125
.
53.
Deng
H
,
Liu
M
,
Yuan
T
, et al
.
Efficacy of humanized anti-BCMA CAR T cell therapy in relapsed/refractory multiple myeloma patients with and without extramedullary disease
.
Front Immunol
.
2021
;
12
:
720571
.
54.
Li
W
,
Liu
M
,
Yuan
T
,
Yan
L
,
Cui
R
,
Deng
Q
.
Efficacy and follow-up of humanized anti-BCMA CAR-T cell therapy in relapsed/refractory multiple myeloma patients with extramedullary-extraosseous, extramedullary-bone related, and without extramedullary disease
.
Hematol Oncol
.
2022
;
40
(
2
):
223
-
232
.

Author notes

Z.S. and J.J. contributed equally to this work.

Presented in abstract form at the 64th annual meeting of the American Society of Hematology, New Orleans, LA, 10-13 December 2022.

The single-cell RNA sequencing and transcriptomic data in this paper have been deposited in the NCBI Gene Expression Omnibus database (accession number GSE230510).

The datasets used and/or analyzed during the current study are available on reasonable request from the corresponding authors, Jianyong Li (lijianyonglm@126.com) and Xiaoyan Qu (quxiaoyan@jsph.org.cn).

The full-text version of this article contains a data supplement.

Supplemental data