Alterations in the Transcriptional Programs of Myeloma Cells and the Microenvironment during Extramedullary Progression Affect Proliferation and Immune Evasion

Multiple myeloma, a plasma cell malignancy in the bone marrow, begins as a monoclonal gammopathy of undetermined significance (MGUS) and eventually progresses to extramedullary disease (1). Extramedullary myeloma commonly affects the pleura, lymph nodes, soft tissues, and liver and is an aggressive disease associated with treatment resistance and shorter patient survival times (2). Early extramedullary or ascites involvement at diagnosis is rare and particularly detrimental.

Cytogenetic and genetic analyses have been performed to identify the molecular characteristics associated with aggressive extramedullary myeloma progression. Cytogenetic abnormalities involving t(4;14) or t(14;16) translocations are early genetic events with poor prognosis, whereas the 17p deletion, which is the most significant genetic prognostic factor in multiple myeloma, occurs during disease progression (3). In addition, MYC translocations and mutations in TP53, CCND1, ATM, ATR, and ZFHX4 have a negative prognostic impact (3). Extramedullary myeloma has a higher incidence of the t(4;14) translocation, 17p deletion, and MYC overexpression and often manifests with different cytogenetics than matched bone marrow myeloma (4), indicating genomic instability. Gene-expression profiling–based risk stratification using many proliferation genes classifies patients with extramedullary myeloma as a high-risk group (5).

The cytogenetic, genetic, and gene-expression profiles of extramedullary myeloma are enriched for poor prognostic alterations, suggesting an intrinsically aggressive nature. As extramedullary myeloma progression requires bone marrow–independent growth and survival, it must gain an alternative source of growth/survival signals and immune suppression signals that would normally be provided by the bone marrow niche (6). The inextricable relationship between myeloma and its microenvironment requires the characterization of both components to understand the progression of myeloma to non-bone marrow compartments.

In this study, we used single-cell RNA sequencing (RNA-seq) to carry out the genomic characterization of diverse cell populations without experimental separation (7). We combined two single-cell RNA (scRNA)-sequencing approaches, high-throughput RNA-seq for the simultaneous analysis of myeloma, and the immune microenvironment and full-length RNA-seq of purified myeloma cells for molecular pathway analyses. We show that extramedullary myeloma cells exhibit autonomous and multiple growth factor–receptor expression that can drive massive cell proliferation in a bone marrow–independent manner. Second, extramedullary myeloma cells possess heterogeneous immune evasion mechanisms against NK and cytotoxic T cells accompanied by the functional impairment of immune cells. By analyzing the myeloma transcriptome aligned with tumor progression, we identified transcriptional programs controlling the progression of aggressive extramedullary myeloma.

This study was approved by the Institutional Review Board (IRB) of Samsung Medical Center (IRB approval no. SMC2013-09-009-012) and carried out in accordance with the principles of the Declaration of Helsinki. Informed written consent was obtained from each subject. The study subjects were 15 Korean patients diagnosed with multiple myeloma at the Samsung Medical Center (Seoul, Korea). Four bone marrow aspirates, one ascites, and two pleural effusion samples were subjected to red blood cell lysis (Qiagen) and cryopreserved in CELLBANK1 (Zenoaq) before 3′ scRNA sequencing. Nine bone marrow aspirate samples were subjected to a Ficoll-Paque PLUS (GE Healthcare) gradient and magnetic separation with anti-CD138 antibody microbeads (Miltenyi Biotech). Ascites (2 patients) and pleural effusion (1 patient) samples were collected and subjected to the same cell isolation process. CD138-positive myeloma-enriched fractions were immediately used for full-length scRNA sequencing, except for one sample (MM26 ascites) that was cryopreserved. Genomic DNA and RNA were purified from CD138-positive and CD138-negative fractions using an ALLPrep kit (Qiagen). Corresponding blood DNA was isolated using a QIAamp DNA Blood Kit (Qiagen).

Genomic DNA from the bone marrow and corresponding blood samples was sheared by Covaris S220 (Covaris) and libraries were constructed using a SureSelect XT Human All Exon v5 and a SureSelect XT Reagent Kit, HSQ (Agilent Technologies) according to the manufacturer's protocol. After multiplexing, the libraries were sequenced on a HiSeq 2500 sequencing platform (Illumina) using the 100 bp paired-end mode of the TruSeq Rapid PE Cluster Kit and TruSeq Rapid SBS Kit (Illumina).

Sequencing reads were aligned to the UCSC hg19 reference genome (http://genome.ucsc.edu) using Burrows–Wheeler Aligner version 0.6.2 with the default options. PCR duplications and data cleanup were marked using Picard-tools-1.8 (http://picard.sourceforge.net/) and GATK-2.2.9 (https://www.broadinstitute.org/gatk/). Single-nucleotide variants of paired samples were identified by MuTect (https://github.com/broadinstitute/mutect) and annotated using ANNOVAR.

Whole-transcriptome sequencing libraries were generated for the RNAs from CD138-negative cell fractions using a TruSeq RNA Sample Preparation v2 Kit (Illumina). Sequencing was carried out using the 100 bp paired-end mode of the TruSeq Rapid PE Cluster Kit and TruSeq Rapid SBS Kit (Illumina).

For full-length scRNA sequencing, CD138-positive myeloma cells were loaded onto a 5 to 10 μm integrated fluidic circuit mRNA sequencing chip for microscopic examination, lysis, reverse transcription, and cDNA amplification in the C1 Single-Cell Auto Prep System (Fluidigm), according to the manufacturer's instructions. RNA spike-ins 1, 4, and 7 from ArrayControl RNA spikes (Thermo Fisher Scientific) were added to the lysis solution according to the manufacturer's instructions. For MM02EM, 96 External RNA Controls Consortium spike-in sequences (Thermo Fisher Scientific) were used instead. Sequencing libraries were constructed using amplified cDNAs with the Nextera XT DNA Sample Prep Kit (Illumina) and sequenced using the 100 bp paired-end mode of the TruSeq Rapid PE Cluster kit and TruSeq Rapid SBS Kit (Illumina). Reads from bulk and C1 RNA-seq (single cells and whole transcriptome) were mapped on to the GRCh37.75 human reference genome using STAR (https://github.com/alexdobin/STAR/releases) version 2.4.0 and transcripts were quantified as transcript per million (TPM) using RSEM version 1.2.18 (http://deweylab.biostat.wisc.edu/rsem/).

For 3′ scRNA sequencing, each single-cell suspension was loaded into a GemCode system version 2 (10X Genomics) targeting 5,000 cells. Sequencing libraries were generated and sequenced on a HiSeq 2500 (Illumina) using the 100 bp paired-end mode. Sequencing reads were mapped on to the GRCh38 human reference genome using the Cell Ranger toolkit (version 2.1.0). scRNA sequencing data for normal bone marrow were downloaded from the HCA Data Portal (https://preview.data.humancellatlas.org/). Before downstream analyses, we applied three quality measures of mitochondrial gene percentages (<20%), unique molecular identifiers (UMI, over 1,000), and gene counts (over 20) calculated from the raw gene-cell–barcode matrix.

For the 3′ scRNA sequencing data, UMI counts were normalized to TPM-like values and used in the scale of log₂ (TPM+1), as described previously (8). Variable genes were selected using the R package Seurat R toolkit (ref. 9; https://satijalab.org/seurat/) and used to compute principal components (PC). Significant PCs were used for cell clustering and t-distributed stochastic neighbor embedding (tSNE) visualization. Before malignant plasma cell (MPC) clustering, the RegressOut function was applied to remove the sample-specific batch effect. Before immune cell clustering, immunoglobulin genes were excluded from the analysis. Differentially expressed genes (DEG) were selected for each cell cluster using the Seurat package with default parameters and Student t tests. Immune cell identity annotations were defined for each cluster by reference-based single-cell annotations in the R package SingleR (https://github.com/dviraran/SingleR; ref. 10). Because SingleR calculates the correlation to the reference transcriptome for each cell, clusters were annotated using the dominant cell type within the cluster.

For the full-length scRNA sequencing data, reference sequences were modified by merging with external spike-in sequences (Thermo Fisher Scientific). Gene quantification data were filtered by eliminating zero-expressed genes and applying the following criteria: (i) in genes, the number of zeros should be less than 449 (mean zero count) across all samples, and (ii) in cells, the number of genes expressed in each sample should be over 1,000. To visualize the myeloma transcriptomic landscape, tSNE and k-means clustering was performed using 11,664 genes and curated gene sets in the molecular signature database (v6.0 MSigDB). To analyze pathway activation, we performed gene set variation analysis (GSVA; ref. 11) and gene set enrichment analysis (GSEA; ref. 12) using hallmark, curated (chemical and genetic perturbations, CGP), and gene ontology gene sets in the MSigDB. For trajectory analysis, monocle 2 (13) was utilized with the default parameters.

Peripheral blood was obtained from 5 healthy donors (average age 60.6) and 10 newly diagnosed myeloma patients (average age 64.8) within the IRB protocol SMC2013-09-009-012. Informed written consent was obtained from each subject. After red blood cell lysis (Qiagen), the whole-cell pellet was cryopreserved in CELLBANK1, rapidly thawed, and stained with PerCP-anti-CD3 (clone UCHT1), PE/Cy7-anti-CD56 (clone 5.1H11), and Alexa647-anti-GZMB (clone GB11) mAbs (BioLegend) according to the manufacturer's protocol. Flow cytometry data were obtained by FACSverse and analyzed by FlowJo software (BD Biosciences).

To understand the cellular dynamics associated with myeloma progression, we first applied massively parallel single-cell RNA-seq (7) to myeloma-stricken bone marrow or extramedullary sites (Fig. 1 and Table 1). The aggregated transcriptome analysis of seven myeloma samples (four bone marrow, one ascites, and two pleural effusion) revealed patient-specific clusters of MPCs and diverse immune cell clusters shared by patients (Fig. 2A and B). The immune cell clusters were mainly comprised of mature T-cell types, NK cells, and monocytes/macrophages, as well as small numbers of B-cell–lineage cells, erythroid cells, and granulocytes (Fig. 2B; Supplementary Fig. S1A; Supplementary Table S1). Myelomatous pleural effusion and ascites samples contained a large number of MPCs and relatively few immune cells likely due to the microenvironment and highly advanced disease state (Fig. 2C, left).

The major immune cell types in the myelomatous bone marrow were similar to those in the healthy bone marrow (7, 14) but contained fewer B-cell populations (B cells and progenitors), indicating deregulated B-cell development (Fig. 2C, right; Supplementary Fig. S1B; ref. 15). Within the T and NK cell types, more CD4⁺ naïve T-cell populations were recovered from the healthy bone marrow, whereas the myeloma bone marrow was populated with more cytotoxic CD8⁺ T cells and NK cells (Fig. 2C, right). To determine whether the expansion of cytotoxic immune cell populations could be systemically detected, we assessed GZMB (granzyme B) expression on the peripheral blood cells (Supplementary Fig. S2). The number of GZMB⁺ T cells or NK cells in the patients with myeloma was highly variable, with CD8⁺ T cells and NK cells from both healthy and myelomatous bone marrow also expressing variable levels of cytotoxic mediator genes (Fig. 2D, top). Moreover, regulatory and exhaustion gene expression (16) was slightly elevated in the patients with myeloma (Fig. 2D, middle), tipping the balance between exhaustion and cytotoxicity toward an exhausted state (Fig. 2D, bottom). These data suggest the activation and persistence of cytotoxic T and NK cells in the patients with myeloma, likely undergoing exhaustion.

The in silico separation of MPCs and subsequent clustering revealed that myeloma subpopulations (Fig. 2E) demonstrated many gene-expression traits implicated in myeloma progression (Supplementary Table S2). The majority of bone marrow myeloma cells shared transcriptome features forming multipatient clusters (Fig. 2E; clusters 0, 7, 12). These bone marrow myeloma–enriched clusters expressed AP-1 (Fig. 2F; cluster 0; JUN, FOS; ref. 17) or TNFRSF17 (clusters 7, 12; ref. 18), genes that promote myeloma growth and survival in the bone marrow microenvironment. In contrast, primary refractory (MM173BM) or extramedullary (MM173EM, MM26EM, and MM36EM) myeloma cells formed patient-specific clusters with heightened intratumoral heterogeneity (Fig. 2E, right and bottom). The clustering results partially reflected the expression patterns of myeloma driver genes (Supplementary Fig. S3A). MM36EM ascites myeloma contained subpopulations with proliferative and poor prognostic (Fig. 2F; cluster 6; CKS1B; ref. 19) gene expression, as well as inflammatory and angiogenic (cluster 1; CCL3 and TMSB4X; refs. 20, 21), immune suppressive (cluster 5; TNFSF10; ref. 22), or potentially dexamethasone resistant (cluster 11; NEAT1; ref. 23) marker genes (Fig. 2F; Supplementary Table S2). MM173BM and MM173EM myelomas shared MYBL2^hi subpopulations (Fig. 2F; clusters 2, 8, 10) with proliferative (cluster 10; UBE2S; ref. 19) or metabolic (cluster 8; PKM; ref. 24) gene expression, while the MM173 metabolic subpopulation exclusively emerged during extramedullary progression (Supplementary Fig. S3B). MM26EM effusion myeloma contained WFDC2⁺ (cluster 3) subpopulations recently identified by single-cell RNA-seq (25) and additional ADM⁺ (cluster 4) hypoxia-induced subpopulations (26). We suspect that the number of myeloma subpopulations and their transcriptional states likely differ for individual patients and their disease status. Nonetheless, intratumoral heterogeneity suggests that myeloma survival and progression is a collaborative process between heterogeneous myeloma subpopulations that drive cancer cell proliferation, migration, immune evasion, and/or drug resistance.

To understand the molecular pathways activated during the progression of aggressive myeloma, we focused on the single-cell transcriptome analysis of CD138-purified myeloma cells (Fig. 1 and Table 1). We obtained bone marrow biopsies from nine patients with stable or primary refractory disease and extramedullary samples of pleural effusion or ascites from four patients. Bone marrow and extramedullary samples were serially collected from two patients. Single-cell transcriptome data were generated using the C1 Single-Cell AutoPrep System, which allowed full-length RNA-seq with a higher coverage of lowly expressed genes than 10X system (Supplementary Fig. S4A). The mutational landscape of our patient cohort was determined by whole-exome sequencing, which identified KRAS and NRAS mutations as the most frequent genetic aberrations (Table 1). Many patients carried oncogenic translocations such as t(4;14) IgH–WHSC1, t(11;14) IgH–CCND1, or t(14;20) IgH–MAFB fusions. Thus, most patients had well-characterized genetic aberrations frequently found in multiple myeloma (27).

Following data quality filtration, we acquired full-length RNA-seq data from 492 cells, with an average of 4.7 million sequencing reads and 4,924 genes in each myeloma cell. Externally spiked-in RNA sequences (Supplementary Fig. S4B) show a constant level of detection across all cells and patients, ensuring high-quality single-cell sequencing with a low batch effect. Individual myeloma cells expressed clonotypic immunoglobulin sequences within a single patient (Supplementary Fig. S4C left; ref. 28) and plasma cell–specific transcription factors (PRDM1, IRF4, and XBP1) with a lack of mature B cells (PAX5 and BACH2; Supplementary Fig. S4C, right; ref. 29). Transcripts for plasma cell-surface molecules, such as SDC1 (CD138), CD38, and SLAMF7, were consistently expressed in most myeloma cells and are candidate therapeutic targets for multiple myeloma (30).

Using the full-length single-cell transcriptome data, we performed clustering analysis on all the bone marrow and extramedullary myeloma cells from 11 patients (Fig. 3A). We eliminated immunoglobulin genes from the clustering analysis as their high levels of clonotypic expression hinder the identification of other molecular pathways causing intertumoral or intratumoral heterogeneity. tSNE projection visualized extramedullary myeloma cells with distinctive patient-specific clusters, whereas most of the bone marrow cells formed common multipatient clusters (Fig. 3A). The expression of the myeloma driver genes partially contributed to the clustering results (Fig. 3B). Correlation analysis at the gene expression (Fig. 3C, left) or pathway level (right) substantiated the overall separation of the extramedullary myeloma cells from the bone marrow myeloma cells. At the pathway level, the primary refractory bone marrow myeloma cells (MM17 and MM34) resembled those of the extramedullary myeloma cells. To extract features of aggressive extramedullary myeloma progression, we made pair-wise comparisons between the bone marrow and extramedullary myeloma cells (Fig. 3D and Supplementary Table S3). We identified transcriptional changes in the upregulation of proliferation/cell-cycle progression, glycolysis, oxidative phosphorylation, proteasome, and antigen presentation genes upon aggressive extramedullary myeloma progression, while bone marrow myeloma cells expressed higher levels of “TNFα-induced NF-κB pathway” genes. All pathways have also been identified by other single-cell analysis methods (Supplementary Fig. S4D; ref. 31). Intriguingly, we found gradual unidirectional changes in most pathway signatures of bone marrow cells from stable disease toward primary refractory and extramedullary progression (Fig. 3E). The NF-κB pathway was the exception, being upregulated in primary refractory bone marrow cells but downregulated in most extramedullary cells. Complex regulation of the NF-kB pathway was observed from correlation analysis between the pathways (Fig. 3F), suggesting that the NF-κB pathway plays a differential role in myeloma progression within the bone marrow or extramedullary microenvironment.

To assess the importance of five gene-expression signatures (proliferation, glycolysis, oxidative phosphorylation, antigen processing and presentation, and proteasome) selected from the single-cell transcriptome data, we determined their impacts on patient survival using two cohorts with extensive clinical annotation: GSE24080 (32) and GSE9782 (33). In the newly diagnosed multiple myeloma GSE24080 cohort, the proliferation signature had poor prognostic value for overall survival (Fig. 4A; Supplementary Table S4A), while in the relapsed multiple myeloma GSE9782 cohort, antigen presentation, proteasome, glycolysis, and oxidative phosphorylation signatures predicted poor prognosis (Fig. 4B; Supplementary Table S4B). The proliferation signature had less of an impact on overall survival in relapsed myeloma, suggesting that increased proliferation is involved in primary relapse, whereas upregulated antigen presentation, proteasome, glycolysis, and oxidative phosphorylation have a more profound effect on the later aggressive progression and/or treatment resistance of relapsed myeloma.

To elucidate the mechanisms driving massive proliferation in aggressive extramedullary myeloma cells, we examined the expression levels of myeloma growth and survival factors. Of the 148 genes known to affect myeloma cell growth (6, 34), we detected 32 genes expressed by the myeloma cells (Fig. 5A). Most of the growth factors, including CCL3 and IL6, were more prominently expressed in the extramedullary myeloma cells (Fig. 5B). Cognate receptor expression for these growth/survival factors may stimulate tumor growth in an autocrine or paracrine manner independently of the bone marrow microenvironment. Indeed, many extramedullary myeloma cells simultaneously expressed CCL3 and CCR1 or IL6 and IL6R, suggesting a potential auto-/paracrine loop in the growth factor signaling axis (Fig. 5C). The expression of the CCL3 and IL6 genes was also assessed in the 10X data set (Supplementary Fig. S4E).

VLA4 (α4β1), formed by the ITGA4 and ITGB1 gene products, is the major integrin in the bone marrow and mediates myeloma cell homing and survival by interacting with the bone marrow microenvironment (35). The primary ligands of VLA4 are VCAM1 and fibronectin (FN1 gene). We found that ITGA4 and ITGB1 are often coexpressed by extramedullary myeloma cells (Fig. 5C right), whereas FN1 expression was detected in the CD138-negative nonmyeloma cell fractions of myelomatous ascites (Supplementary Fig. S4F), suggesting that myeloma growth signals are delivered to VLA4 by the extramedullary microenvironment (36). Altogether, these multiple auto- and paracrine growth signals could alleviate the dependency of myeloma growth on the bone marrow microenvironment (3, 34).

The bone marrow niche provides myeloma survival signals and an immune-suppressive microenvironment (6). Hematopoietic cells, including myeloid-derived suppressor cells and regulatory T cells, are known to be strongly involved in creating an immune-suppressive microenvironment (37). In Fig. 2, we demonstrated a significant number of cytotoxic T lymphocytes and NK cells in the myelomatous bone marrow and ascites, despite the elevated expression of exhaustion molecules. In myeloma, various mechanisms of immune evasion were detected at the gene-expression level. First, MPCs demonstrated consistent classic (HLA-A, B, and C) and nonclassic (HLA-E and CD1D) MHC class I molecule expression (Fig. 5D), which could provide a mechanism of protection from NK-mediated cytotoxicity (38, 39). Sustained classic and nonclassic MHC I expression in myeloma contrasts with the frequent loss observed in other types of solid cancers (40). The expression of leukocyte immunoglobulin-like (LIL) family proteins LILRB1 and LILRB4 may also deliver inhibitory signals to NK cells (41, 42). Immunomodulatory gene expression by myeloma cells was also assessed in the 10X data set (Supplementary Fig. S4E), demonstrating that myeloma progression elicits multiple and heterogeneous NK-evasion pathways for tumor cells.

Overall, patients with myeloma demonstrated a shift in T-cell composition from a naïve toward a cytotoxic/exhausted phenotype, suggesting T-cell activation followed by a gradual waning (Fig. 2D). In myeloma samples, we found several gene-expression characteristics associated with evasion from cytotoxic T lymphocytes (CTL; Fig. 5D). First, MHC class I gene expression was variably regulated in the paired extramedullary samples, suggesting that myeloma cells flexibly balance MHC class I levels to avoid NK- or T-cell–mediated cytotoxicity. Second, some extramedullary myeloma cells expressed the TRAIL gene (TNFSF10), which may engage TRAIL receptors on T cells and trigger their aberrant activation (22). Finally transcription factor hypoxia-inducible factor-1-alpha (HIF1A) gene expression was upregulated in extramedullary myeloma; the exposure of cancer cells to hypoxia is known to induce the expression of programmed cell death ligand-1 (PD-L1), which increases tumor cell resistance to CTL-mediated lysis (43). Overall, these data suggest that heterogeneous mechanisms of immune evasion accompany myeloma progression.

After primary bone marrow biopsy, most patients received a combination treatment including the proteasome inhibitor bortezomib (Fig. 1). All primary refractory bone marrow myeloma (MM17, MM34, and MM173) and extramedullary samples were collected during the refractory (relapsed) phase. Multiple mechanisms can confer bortezomib resistance, among which proteasome activation occurs most frequently, involving the overexpression of proteasome components or mutations in the PSMB5 gene (44). GSEA detected the overexpression of proteasome components in myeloma subpopulations (Fig. 3D). Some also expressed high levels of proliferation-associated genes, which is the most validated myeloma signature for poor prognosis (45). Therefore, myeloma cell subsets with both proliferation and proteasome signature expression may be highly aggressive and drug resistant (Fig. 6A and B). Single-cell analysis revealed that heterogeneity in gene expression reflects the state of myeloma subpopulations and potentially drug resistance and disease progression (46).

To examine the relationship between tumor heterogeneity and evolution in extramedullary disease, we performed a trajectory analysis for paired bone marrow–extramedullary samples, MM02(EM), MM34(EM), and MM173(EM; Fig. 6C). Both bone marrow and ascites/effusion samples from MM34 and MM173 patients were obtained during the highly aggressive refractory disease state, with an extremely short interval between the two biopsies (Fig. 1). In comparison, the MM02 bone marrow biopsy was obtained during an early disease phase and the effusion sample was collected during the relapsed refractory phase two years later. Trajectory analysis (13) of each patient revealed 3 to 5 states and visualized the evolution of the myeloma transcriptome (Fig. 6C). Pseudotime analysis indicated a branched evolution from the bone marrow to the extramedullary myeloma. Intriguingly, one of the extramedullary branches was more remotely related to bone marrow myeloma and was enriched in proliferation or proteasome gene expression compared with the others (Fig. 6C, bottom). These results suggest that intratumoral heterogeneity is linked with the evolution of more aggressive myeloma subpopulations and disease progression.

By combining massively parallel 3′ scRNA sequencing (10X platform) and low-throughput full-length RNA-seq (C1 platform), we characterized the molecular events accompanying aggressive extramedullary myeloma progression in the context of the tumor-immune microenvironment. Massively parallel RNA-seq illustrated the overall tumor landscape encompassing the myeloma and immune microenvironment, whereas in-depth analyses of purified myeloma cells allowed the molecular characterization of individual myeloma cells. The two approaches also compensated limitations in each method, the small number of cells and cell size limitations in the C1 or extremely sparse gene expressions in the 10X platform.

Visualization of the global transcriptome landscape demonstrated the presence of diverse myeloma subpopulations surrounded by cytotoxic CD8⁺ T cells and NK cells. In previous studies, T-cell and NK-cell defects were demonstrated in numbers, phenotypes, and function for advanced myeloma patients (47). Importantly, therapeutic strategies reversing NK or T-cell dysfunction or enhancing their activity have gained much attention and achieved clinical benefits in some cases. By comprehensive transcriptome profiling at a single-cell resolution, we demonstrated differences in T cells and NK cells, with a shift from naïve toward an activated/exhausted state. We also identified dynamic regulation of MHC class I gene expression in myeloma cells and suggested additional mechanisms of NK and T-cell inhibition, including LILRB1/4-mediated NK inhibition, TRAIL (TNFSF10)-mediated immune cell regulation, and HIF1A-mediated resistance from cytotoxic T-cell killing. These potential mechanisms varied among patients, suggesting differential immune cell status and the susceptibility of cytotoxic T cells and NK cells to reactivation strategies.

The most dramatic changes accompanying aggressive myeloma progression were the autocrine and pleiotropic gene-expression phenotype for myeloma growth factor receptors. In particular, the vast majority of extramedullary myeloma cells and refractory bone marrow myeloma cells expressed IL6 receptor transcripts, unlike the sparse expression among primary myeloma cells in the bone marrow (48). Many extramedullary myeloma cells simultaneously expressed IL6, suggesting autocrine activation and growth independency. In the bone marrow, IL6 is mainly provided from the myeloma niche by a paracrine source (34). IL6 production and downstream IL6R signaling can be amplified via other molecular interactions; for instance, the VLA4 integrin (ITGA4 and ITGB1 gene expression) in myeloma cells and fibronectin (FN1 gene expression) in the CD138-negative compartment of extramedullary myeloma can augment IL6-mediated growth (49). Overall, these data demonstrate multiple mechanisms of IL6 signaling reinforcement involving myeloma cells and the hematologic microenvironment during extramedullary progression. Despite the prominent role of IL6 in myeloma progression, a clinical trial adding the anti-IL6 monoclonal antibody siltuximab to the bortezomib–melphalan–prednisone regimen did not improve the complete response rate or long-term outcome of transplant-ineligible patients with newly diagnosed multiple myeloma (50). Our study suggests this may be due to heterogeneity of the IL6 receptor expression among myeloma cells in those patients. For aggressive extramedullary myeloma, therapeutic strategies targeting multiple growth factors would be required to control pleiotropic growth stimulation.

In summary, we were able to demonstrate distinct gene-expression characteristics in refractory (relapsed) myeloma from pleural effusion or ascites consisting of pleiotropic auto- and paracrine-driven proliferation accompanied by the weakening of and evasion from the activated immune microenvironment. Analysis of intratumoral heterogeneity further demonstrated that preexisting myeloma cells had the potential for tumor progression and drug resistance. Our study has limitations of the small number of patients and zero-inflated gene-expression patterns in the single-cell RNA-seq data, in addition to the confounding factors for bone marrow versus extramedullary myeloma groups including the nature of microenvironment and the drug treatment for the extramedullary myeloma. Despite these limitations, our comparative profiling of myeloma and the immune microenvironment revealed aspects and limitations of therapeutic interventions that could be used to design molecular therapeutic strategies for myeloma patients with aggressive extramedullary involvement.

The scRNA sequencing and transcriptomic data in this paper have been submitted to the NCBI Gene-Expression Omnibus (GEO; http://www.ncbi.nlm.nih.gov/geo/) under accession numbers GSE106218 and GSE110499.

No potential conflicts of interest were disclosed.

Conception and design: H.-O. Lee, K. Kim, W.-Y. Park

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): S.J. Kim, H.-J. Kim, H.-O. Lee, K. Kim

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): D. Ryu, S.J. Kim, Y. Hong, N. Kim, H.-J. Kim, H.-O. Lee, K. Kim, W.-Y. Park

Writing, review, and/or revision of the manuscript: D. Ryu, S.J. Kim, H.-O. Lee, K. Kim, W.-Y. Park

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): A. Jo, K. Kim, W.-Y. Park

Study supervision: H.-O. Lee, K. Kim, W.-Y. Park

This work was supported by grants from the Ministry of Education (NRF-2017R1D1A1B03032194), the Ministry of Science and ICT (NRF-2016R1A5A1011974, NRF-2017M3A9A7050803, NRF-2017M3C9A6044636, and 2018M3C9A6017315), and the Ministry of Food and Drug Safety (16173MFDS004), Republic of Korea.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.