Single cell RNA sequencing of plasma cells from guided biopsies of osteolytic lesions and corresponding bone marrow identifies inter-patient heterogeneity
We implemented a translational workflow (Figure 1A) to obtain and purify viable PC from OL and corresponding BM samples from patients with newly diagnosed or relapsed MM. The bone marrow biopsies from the iliac crest (bone marrow sample, BM) were obtained and processed at the same time as OL biopsies since differences in sample processing times might cause changes in MM gene expression (14). We performed scRNA-seq on paired samples from 10 patients (7 with newly diagnosed and 3 relapsed/refractory MM) with one patient having 2 OL biopsied. In three patients with NDMM, we obtained additional samples after induction therapy. Patient characteristics are summarized in Table 1.
Clustering of 148,746 single cells from BM and OL (median 7712 cells/sample) created a map of distinct populations based on transcriptomes from individual patients (Figure 1B). There were no significant differences between both locations regarding the purity of isolated and sequenced PC, underlining the feasibility of our protocol and comparability of paired samples. In total, 94.8% of cells from OL (n = 70,036) and 95.7% of cells from BM (n = 71,580) were PC (Figure 1C). Since the cell cycle can introduce heterogeneity in gene expression in very homogeneous cell populations like in our study, we assigned cell cycle scores to PC and found no significant differences between OL and BM (Figure 1D).
Next, we identified marker genes for malignant PC clusters from each individual patient (Figure 1E). As expected, the genes for cyclin D1-3 (CCND1-3) were preferentially expressed in patients with IgH translocations (NDMM02, NDMM04 and NDMM06) as detected by FISH (Table 1). Furthermore, genes associated with myeloma bone disease (DKK1 and FRZB), cytokine signaling (IFI27 and IL6R) and EDNRB were identified as marker genes for PC clusters from individual patients. All of the aforementioned genes were previously described as characteristic genes elevated in the molecular subtypes of MM that were derived from bulk GEP (15, 16). Besides these known marker genes, we identified STMN1 to be preferentially expressed in small subsets of malignant PC (Figure 1E). No significant differences for the identified marker genes were found between BM and OL.
Following single gene analysis, we performed a gene set enrichment analysis (GSEA) using the curated MM subtype gene sets from the Molecular Signature Database (MSigDB). Consistent with the single gene analysis, gene-set analysis showed that patients could be grouped according to the molecular classification of MM. For example, patient RRMM02 in the low bone disease (LB) group showed higher expression of IL6R and EDNRB and NDMM02 in the MAF group (MF) showed higher levels of CCND2 and lower levels of DKK1 (Figure 1F). This demonstrates that scRNA-seq reveals inter-patient heterogeneity and identifies MM subtypes for individual patients.
Intra-patient heterogeneity based on scRNA-seq and inferred CNVs
After investigating inter-patient heterogeneity using scRNA-seq data, we aimed at deciphering intra-patient heterogeneity on the transcriptional and chromosomal level. We performed clustering and differential expression analysis of PC for each individual patient. Specifically, we looked for genes that were recurrently differentially expressed in PC clusters and that have been associated to play a role in the pathogenesis and prognostication of MM. A summary of the differential expression analysis for each individual patient is available at Supplemental Table 2.
In 9 out of 10 patients, we identified clusters of malignant PC characterized by the overexpression of genes encoding the microtubule-associated proteins STMN1 and TUBA1B (Figure 2A and Supplemental Figure 1). STMN1 is among the 15 genes associated with high-risk disease in the GEP score identified by the Integroupe Francophone du Myelome (4). Furthermore, higher expression levels of STMN1 and TUBA1B are associated with shorter progression-free and overall survival in the CoMMpass dataset of the Multiple Myeloma Research Foundation (Supplemental Figure 2). This underlines the biological and prognostic relevance of the identified subclusters. Since the respective clusters accounted only for a small number of malignant PC in each individual patient (Supplemental Figure 1) they would have been missed by bulk sequencing.
CNVs are another cornerstone in the biology of MM and we have shown recently that cytogenetically defined subclones have major prognostic implications (17). However, in routine practice, chromosomal aberrations are usually detected by fluorescence in situ hybridization (FISH). The detection of subclones by FISH is limited by the number of analyzed cells and applied probes. We leveraged the scRNA-seq data to identify genome-wide CNVs and subclones for 141,616 PC using inferCNV (18). Figure 2B provides two examples for CNVs from a patient with newly diagnosed disease (NDMM01) and a patient with relapsed disease (RRMM01). UPhyloPlot2 (19) was used to draw phylogenetic trees from inferCNV outputs. In 6 of 7 with NDMM and 2 of 3 patients with RRMM, we were able to detect subclones based single cell transcriptomes (Figure 2C and Supplemental Figure 3). No significant differences in CNVs were found when comparing OL and BM.
Limited evidence for spatial heterogeneity from whole-exome sequencing
Next, we aimed at characterizing spatial heterogeneity. Whole-exome sequencing (WES) has been utilized to investigate spatial genomic heterogeneity in MM based on banked, frozen PC from a retrospective study (7). To explore whether our study from scRNA-seq on freshly isolated PC could reveal another layer of complexity of spatial heterogeneity, we performed WES on all paired samples in our study. The comparison of PC with the matched germline normal identified a total of 1103 somatic mutations, including 1063 SNVs and 40 Indels (Figure 3A). Among these somatic mutations, 665 were predicted to cause an amino acid alteration, 72 were truncating and 366 were silent mutations. Supplemental table 2 gives an overview of the individual variant calls from the WES analysis. For each patient, we quantified the similarity between the BM and OL by calculating the Jaccard score, defined as the ratio between shared mutations and all mutations. Patient NDMM03 was excluded from calculations since only limited numbers of malignant PCs were captured from the BM.
In 8 out of 10 patients (NDMM02-07 and RRMM02-03), the percentage of shared mutations between the bone marrow and OL were ~ 80% or higher, suggesting the BM and OL were highly consistent (Figure 3 A and B). In the patient with 2 OL (NDMM06), the Jaccard scores among all 3 locations (BM and two OL) were above 97% with a clonal TP53 mutation present in all three locations.
Only in 2 out of 10 patients, WES revealed relevant spatial heterogeneity:
In patient NDMM01, 75% of all mutations were shared between the bone marrow and OL, and 25% of mutations were only present in the OL (Figure 3A) including a BRAF V600E mutation. In patient RRMM01 with an OL of the right clavicle with extramedullary spread, only 20% of all mutations were shared, with 24% of the mutations found only in the bone marrow, and 56% of the mutations unique to the OL (Figure 3A). Two distinct BRAF mutations were detected: V600E was only in BM, and another activating Class 2 BRAF mutation (G469R) in the OL. The latter mutation causes resistance to the BRAF inhibitor vemurafenib (20). Furthermore, we detected an additional NRAS mutation (G12D) in the OL. NRAS mutations have been associated to drive spatially divergent resistance to vemurafenib in BRAFmut MM (21). Both examples demonstrate that treatment with a BRAF inhibitor would most likely be ineffective against PC from different locations in individual patients.
Single cell RNA sequencing revealed another layer of complexity and links site-specific gene expression to the development of osteolytic lesions
After analyzing WES from malignant PC from the OL and BM, we investigated whether the observed similarities between both locations are also reflected by scRNA-seq data. Therefore, average gene expression of PC from the OL and BM were correlated to each other (Figure 3B). In agreement with findings from WES, we found a high concordance of average gene expression between both locations. However, in every patient outliers in both directions were observed (Figure 3B).
To identify genes that are differentially expressed in PC from both conditions, we performed an integrated analysis after anchoring datasets from OL and BM for each individual patient as described before (22). After applying the integration procedure, malignant PC were robustly detected in all datasets and the same PC clusters were identified in BM and OL (Figure 3C). In all patients, we were able to find marker genes that were differentially expressed in malignant PC from OL compared to BM. Overall, 1140 genes were identified that were differentially expressed between OL and BM (Figure 3D). Genes that have been associated with the development of myeloma bone disease such as DKK1, HGF and TIMP-1 (23, 24) were among the markers with higher expression levels in OL (Figure 3D). Furthermore, in agreement with the first scRNA-seq analysis in MM (10), we found LAMP5 to be upregulated in PC from the OL in patient NDMM03 (Figure 4B).
Genes that were recurrently downregulated in PC from OL were JUN and FOS (6 of 10 patients), dual specificity phosphatase 1 (DUSP1, 5 of 10 patients) and hemoglobin beta chain (HBB, 3 of 10 patients) (Figure 3D and E). Importantly, no somatic mutations were detected in the respective genes in PC from the BM and OL. Downregulation of JUN/FOS, DUSP1 and HBB has been connected to extramedullary spread of MM in the past (25). Furthermore, JUN/FOS are linked to the malignant transformation of B-cells (26) and dissemination of clonal PC in a preclinical model (27). GSEA confirmed that downregulation of genes in pathways associated with normal B-cell was common in PC from the OL (Figure 3F).
Additionally, lower expression levels of genes encoding for the non-restricted light (e.g. IGLC2/3 in NDMM01, Figure 4A) and heavy chain (e.g. IGHM in NDMM03, Figure 4B) were observed in PC from OL. Downregulation of the affected heavy chain was observed in PC from the OL in RRMM01 and RRMM03 providing a first link between evolving disease and the rare phenomenon of light chain escape that can be observed in heavily pretreated RRMM and EMD.
In patient RRMM01 with a PC tumor with extramedullary spread, we identified the gene for zinc-alpha2-glycoprotein (AZGP1) to be downregulated in the OL (Figure 4D). AZGP1 is a known tumor suppressor gene causing mesenchymal-to-epithelial transition and inhibiting invasion and metastasis (28). While its expression was virtually absent in PC from the OL, it was homogeneously expressed in the magnitude of PC from the BM. However, a small cluster of PC with lower expression levels compared to the majority of PC was identified in the BM (Figure 4D). Trajectory inference revealed that PC from the respective cluster underwent the largest transcriptional change, might have occurred latest in the developmental process and might have given rise to the OL in the right clavicle (Figure 4D).
Our results show that scRNA-seq adds another layer of complexity compared to WES to spatial heterogeneity in MM. The scRNA-seq data allow us to link site-specific gene expression to the development of myeloma bone disease and identify subclusters that might be the origin for OL.
Single cell RNA sequencing identifies minimal residual disease and transcriptional changes after therapy
In three patients (NDMM01, NDMM03 and NDMM06) we collected samples after 4 cycles of induction therapy. While in NDMM01 (Figure 4A) and NDMM06 (Figure 4B) we performed a regular bone marrow biopsy to assess residual disease, an imaging-guided biopsy of a residual focal MRI lesion in T8 was biopsied in NDMM03 (Figure 4C). NDMM01 and NDMM03 were in MRD-positive complete remission (CR) after 4 cycles RVD. NDMM06 had achieved a partial response (PR) after 4 cycles daratumumab-RVD. Correspondingly, we captured less malignant PC in NDMM01 (Figure 4A) and NDMM03 (Figure 4B), while almost the same amount of malignant PC after therapy were isolated and sequenced in NDMM06 compared to the BM and the two OL at primary diagnosis (Figure 4C).
To identify transcriptional changes after therapy, we performed a differential expression analysis after integrating datasets before and after therapy. In both patients in MRD-positive CR, higher expression levels of HLA-DRA, HLA-DRB1 and HLA-DPB1 were observed in malignant PC after 4 cycles RVD.
In patient NDMM06 in PR, ISG15 was upregulated in residual PC after 4 cycles of daratumumab-RVD. ISG15 encodes an ubiquitin-like protein and its expression has been linked to Carfilzomib-resistance in a recent preclinical study (29). Additionally, TXNIP and DDIT4 were downregulated upon therapy. Both genes are induced by dexamethasone and involved in glucocorticoid-induced apoptosis providing a link to steroid-refractoriness in the detected residual disease (30, 31).
This demonstrates that scRNA-seq identifies MRD in patients in remission and reveals transcriptional changes consistent with immunomodulatory effects of lenalidomide and drug resistance to proteasome inhibitors.