Identification of a drug-response gene in multiple myeloma through longitudinal single-cell transcriptome sequencing.

Despite recent therapeutic advances for multiple myeloma (MM), relapse is very common. Here, we conducted longitudinal single-cell transcriptome sequencing (scRNA-seq) of MM cells from a patient with relapsed MM, treated with multiple anti-myeloma drugs. We observed five subclusters of MM cells, which appeared and/or disappeared in response to the therapeutic pressure, and identified cluster 3 which emerged during lenalidomide treatment and disappeared after proteasome inhibitor (PI) treatment. Among the differentially expressed genes in cluster 3, we found a candidate drug-response gene; pellino E3 ubiquitin-protein ligase family member 2 (PELI2), which is responsible for PI-induced cell death in in vitro assay. Kaplan-Meier survival analysis of database revealed that higher expression of PELI2 is associated with a better prognosis. Our integrated strategy combining longitudinal scRNA-seq analysis, in vitro functional assay, and database analysis would facilitate the understanding of clonal dynamics of MM in response to anti-myeloma drugs and identification of drug-response genes.

Introduction

Multiple myeloma (MM) is a life-threatening, chronic, and incurable hematological malignancy, which is characterized by the expansion of neoplastic plasma cells in the bone marrow (BM). In the last two decades, novel classes of anti-myeloma drugs have been developed; proteasome inhibitors (PI), immunomodulatory drugs (IMiDs), and targeted monoclonal antibodies, which have dramatically improved the quality and length of life of patients with MM (Anderson, 2016). PI inhibits proteasome-mediated protein degradation, thus causing cell death (Adams, 2004; Hideshima et al., 2001), and plays a critical role in the backbone therapy in MM. MM cells might be highly sensitive to PI, as they produce and secrete a large amount of proteins such as immunoglobulins (Obeng et al., 2006). However, despite the recent therapeutic advances, therapy for relapsed and refractory MM has largely been empiric so far, partly because of no useful predictive biomarkers for therapeutic efficacy. Indeed, long-term exposure to drugs with no therapeutic effects often causes relapses in patients with MM. It is, therefore, desirable to identify prognostic biomarkers and even drug-response genes which confer drug sensitivity or resistance. However, given that MM cells originate from not only plasma blasts but also pre-B cells (Barwick et al., 2019), there might be marked heterogeneity of MM clones, thus hampering the identification of such biomarkers and drug-response genes. To exclude the effect of the inter-patient heterogeneity, it would be helpful to analyze a single case with longitudinal sampling at multiple time points.

MM is composed of heterogeneous cell populations, which are induced by linear and branching phylogenies. Indeed, the heterogeneity of MM cells undergoes genetic diversification and selection, thus resulting in marked changes in the population of clones and clonal evolution (Morgan et al., 2012). It is, therefore, difficult to characterize minor clones in a bulk transcriptome sequencing analysis. As genomic abnormalities might accumulate during the disease progression, it is critical to elucidate cellular heterogeneity at multiple time points during the course of treatment, so that we can understand the cellular mechanisms responsible for the drug response.

Single-cell transcriptome sequencing (scRNA-seq) is the state-of-the-art technology capable of dissecting cellular heterogeneity of transcriptome at a single-cell resolution, thereby providing new insight into the subpopulation structures of tumors (Ledergor et al., 2018). Recently, Cohen et al. reported through scRNA-seq analysis of BM cells obtained from 20 patients with MM that MM is composed of genetically diverse clones with marked clonal dynamics and transcriptional changes before and after the treatment (Cohen et al., 2021). We have speculated that scRNA-seq analysis coupled with longitudinal sampling and detailed genetic profiling would facilitate the understanding of drug-induced changes in the clonal dynamics and identification of drug-response genes.

Here, we conducted a longitudinal scRNA-seq analysis of MM cells from a single patient with relapsed and refractory MM, who had been treated with multiple anti-myeloma drugs. We found a small cell cluster with distinct transcriptomic features, which emerged during lenalidomide treatment and disappeared after the addition of a PI, ixazomib. Among the differentially expressed genes (DE-Gs), we identified pellino E3 ubiquitin-protein ligase family member 2 (PELI2) as a PI-sensitive gene, which confers the PI sensitivity in in vitro cultured cell assays. Kaplan-Meier survival analysis of datasets from the publicly accessible the database revealed that high expression of PELI2 is associated with good prognosis. Thus, our integrated strategy combining longitudinal scRNA-seq analysis of a single patient with in vitro cultured cell studies and database analysis would be useful for the identification of biomarkers and drug-response genes in hematological malignancy with marked heterogeneity.

Results

Case description

To elucidate the drug-induced changes in MM cells over a long period of time, we recruited a patient with relapsed and refractory MM, who had been treated with multiple anti-myeloma drugs; cyclophosphamide, dexamethasone, PI, IMiDs, and anti-CD38 antibody for 6 years. We collected BM samples every time of relapse over a long-term clinical course, so that we might understand the subpopulation dynamics in response to anti-myeloma drugs.

A 67-year-old man, who had undergone a routine medical examination, presented with anemia and increased levels of total protein in March 2014. He was admitted to Kyushu University Hospital, diagnosed with IgA-kappa MM and the MM cells were positive for IGH/FGFR3 by fluorescence in situ hybridization. The baseline disease characteristics of the patient are documented (Table 1). Before treatment, the first BM sample (termed P1) was collected. The clinical course and time points of BM aspiration are depicted (Figure 1).

Table 1

Summary of peripheral blood laboratory data at baseline

Variables	Value	Reference range
Peripheral blood
Red blood cell count (× 106/mm3)	3.48	4.00–5.39
Hemoglobin (g/dL)	10.8	14.0–18.0
Hematocrit (%)	32.6	36.0–52.0
White blood cell count (/mm3)	2,640	3,500–9,000
Differential count (%)
Neutrophils	44.1	40.0–70.0
Lymphocytes	41.4	18.0–53.0
Monocytes	7.3	2.0–12.0
Eosinophils	2.8	1.0–4.0
Basophils	1.0	0.0–1.0
Platelet count (× 103/mm3)	184	140–440
Reticulocytes (%)	0.6	0.5–1.5
Urea nitrogen (mg/dL)	18	8–22
Creatinine (mg/dL)	1.0	0.6–1.1
Estimated GFR (ml/min/1.73 m2)	58	>90
Total protein (g/dL)	9.8	6.7–8.3
Albumin (g/dL)	3.1	4.0–5.0
Calcium (mg/dL)	9.0	8.7–10.3
Lactate dehydrogenase (U/liter)	104	119–229
β2-microglobulin (μg/mL)	2.7	0.8–1.7
IgG (mg/dL)	321	872–1815
IgA (mg/dL)	5366	95–405
IgM (mg/dL)	9	31–190
Free kappa light chain (mg/L)	12.1	3.3–19.4
Free lambda light chain (mg/L)	0.2	5.7–26.3
Free kappa: lambda ratio	8.64	0.26–1.65
C-reactive protein (mg/L)	0.2	<0.1
Bone marrow
Cellularity	11.4% clonal plasma cells
Karyotype	46, XY
FISHa analysis (%)	IGH/FGFR3 (30.0)

Fluorescence in situ hybridization.

Figure 1

Bone marrow aspiration sampling and clinical course of the patient

(Top) Bone marrow aspiration samplings (P1-P5). P2 was not used for analysis, because the cell count was very low. (Middle) Colored boxes indicate the periods of treatments with PI (red), IMiDs (blue), mAb (yellow), CPA (green), DEX (gray). (Bottom) The IgA levels at diagnosis and for the following 6 years. The dashed gray line indicates the normal upper level of IgA. Auto-PBSCT, Autologous peripheral blood stem cell transplant; PI, proteasome inhibitors; Bor, bortezomib; Ixa, ixazomib; IMiDs, Immunomodulatory drugs; Thal, thalidomide; Len, lenalidomide; mAb, monoclonal antibodies; Dara, daratumumab; CPA, cyclophosphamide; DEX, dexamethasone.

The patient received the induction chemotherapy with bortezomib, cyclophosphamide, and dexamethasone for four cycles, followed by autologous stem cell transplantation with bortezomib and high-dose melphalan. He achieved a very good partial response as defined by the International Myeloma Working Group response criteria, which was consolidated with bortezomib, thalidomide, and dexamethasone for three cycles. Ten months after consolidation chemotherapy, the patient was diagnosed with the first relapse (P2) with elevated IgA. He was treated with lenalidomide plus dexamethasone and thereafter had stable disease. To achieve better response, his therapeutic regimen was switched to ixazomib, lenalidomide, and dexamethasone in October 2017 (P3), so that he obtained a partial response. As he had the disease progression 15 months after the start of ixazomib (P4), his therapy was switched to daratumumab, bortezomib, and dexamethasone, with a very good partial response. After 16 months, his IgA level was elevated and he was diagnosed as relapse (P5). He was treated with elotuzumab, pomalidomide, and dexamethasone, so that he obtained a partial response. However, 3 months later, he developed pneumonia and dizziness as adverse events of elotuzumab and pomalidomide, respectively. We discontinued them and he died 3 months later.

Isolation and validation of multiple myeloma cells

Among the BM samples obtained at P1-5, we did not analyze P2, because of the smaller number of cells. To elucidate the clonal dynamics during the treatment with multiple anti-myeloma drugs, we performed longitudinal scRNA-seq analysis of BM samples collected at the time of diagnosis (P1) and relapse (P3, P4, and P5). The MM cells from each sample were isolated using fluorescence-activated cell sorting (FACS) (Figure 2A and Table S1). After de-multiplexing the BM samples at diagnosis and relapse based on their specific barcode, we obtained 454, 1008, 948, and 984 cells from P1, P3, P4, and P5, respectively. Low-quality cells with a number of expressed genes less than 1000 or those with mitochondrial genes of more than 10% were excluded. After quality control and data normalization among four samples, 390 (85.9%), 984 (97.6%), 852 (89.8%), and 756 (76.8%) cells were subjected to further downstream analysis (Figure S1).

Figure 2

Characterization of the clusters with features of MM cells

(A) Flow cytometry plots showing sorting strategy (CD38++CD45−CD19−) for MM cells after doublet exclusion. Plots were generated using FlowJo software.

(B) Uni-form manifold approximation and projection (UMAP) plot depicting 2982 single bone marrow MM cells obtained at diagnosis and each relapse. Each cluster is represented by a specific color and number.

(C) UMAP plot of a selected set of key plasma cell genes and main MM driver genes over the metacell model. Normalized single-cell expression (unique molecular identifier (UMI) count, log scale) of eight representative genes.

The color and value represented the average expression level.

See also Tables S1, S2 and Figures S1–S3.

Based on the above criteria, we selected MM cells for dimensionality reduction and clustering using the Seurat R package for single-cell analysis (Stuart et al., 2019). We mitigated the effect of cell cycle heterogeneity by using cell cycle scoring in Seurat (Figure S2). Unsupervised clustering of 2982 MM cells sorted from the four samples created a detailed map comprising five transcriptionally homogeneous subpopulations (Figure 2B). We performed cell type annotation with the SingleR R package and confirmed that more than 99% of cells examined have the characteristics of plasma cells (B cells) (Figure S3 and Table S2). We validated whether all clusters express several critical gene groups for MM, including myeloma markers and tumorigenic genes. As expected, myeloma-specific surface markers, such as CD38, SDC1 (CD138), and NCAM1 (CD56) were similarly expressed in all clusters, with the downregulation of PTPRC (CD45) (Figure 2C). Transcriptional factors specific to MM cells, such as IRF4, PRDM1 (Blimp-1), and XBP1 were also detected in MM cells (Borson et al., 2002; Mimura et al., 2012; Shaffer et al., 2008). In this study, FGFR3 was positive for all clusters, which is consistent with the data from fluorescence in situ hybridization for IGH/FGFR3 (Table 1). Collectively, these observations confirmed that cells in all clusters were MM cells.

Longitudinal single-cell transcriptome sequencing analysis

We next examined how the clonal structure evolves during the clinical course of the disease. We observed five UMAP subclusters (clusters 0-4) of MM cells, which appeared and/or disappeared in response to the therapeutic pressure (Figures 3A and S4). Longitudinal scRNA-seq analysis revealed the presence of cluster 0 across P1-P5 samplings, with significant increases in clone sizes of clusters 1 and 2 from P1 to P3. We initially used the CytoTRACE R package to evaluate the transcriptional diversity of each cell cluster, which is a hallmark of developmental potential. As a large proportion of cells with the highest CytoTRACE score, which are the most immature cells, was detected in cluster 0 (Figure 3B, arrow), we defined this area in cluster 0 as the start point of the pseudotime analysis. The evolutionary trajectory derived from the pseudotime analysis using Monocle 3 R package depicted a finely resolved cascade of clonal evolution from cluster 0 to cluster 1 or 2 (Figure 3C). In this study, to understand the influence of treatment on the acquisition of drug resistance, we selected DE-Gs in cluster 0 when comparing P1 and P5, among which the pathway related to cellular responses to stress was highlighted at P5 (Figure S5). We also found that ARHGEF2 and JUNB, which are reported to be candidates as PI-resistant genes (Mosca et al., 2013; Fan et al., 2017), are among the up-regulated genes at P5 (Table S3). It is noteworthy that cluster 3, which emerged at P3, occurred at the time of relapse after 22-month treatment of lenalidomide. The clone size of cluster 3 declined significantly from P3 (7.6%) to P4 (0.4%) after the addition of ixazomib (Figure 3D). In this study, cells expressing CD38 were markedly reduced from P4 to P5 after the treatment with anti-CD38 monoclonal antibody, daratumumab (Figure S6), supporting the validity of our scRNA-seq analysis.

Figure 3

Longitudinal scRNA-seq analysis of MM cells

(A) UMAP plot of MM cell subclusters at P1, P3, P4, and P5 disease stages. The colors indicate different subclusters within every time point. Dotted circle in red indicates cluster 3 at P3.

(B) CytoTRACE analysis depicting the differentiation state of MM cells. The color indicates the state of differentiation; from more differentiated (gray) to less differentiated (red) states. A large proportion of the immature cells is detected in cluster 0 (arrow).

(C) Monocle3 pseudotime analysis showing diverging transcriptional trajectories of MM cells. The color code indicates the cellular moment in the calculated pseudotime. Dotted circle in red indicates the start point of the pseudotime analysis.

(D) Bar plot depicting the ratio of cluster 3. The numbers indicate the fraction of cluster 3.

See also Table S3 and Figures S4–S6.

Identification of differentially expressed genes in cluster 3

To identify genes associated with the ixazomib response, we selected the DE-Gs between cluster 3 and others (Absolute average log2 fold change >0.5 and Adjusted p value <0.05) (124 up-regulated and 38 down-regulated) (Table S4), among which enrichment analysis highlighted dozens of Reactome pathways with the up- or down-regulated genes in cluster 3 (Figure 4A). The most enriched pathway in cluster 3 was cytokine signaling in the immune system. We compared the transcriptomes among the clusters and found that several genes are specifically expressed in cluster 3, being almost absent in other clusters (Figure 4B). Using more stringent criteria (Average log2 fold change >1.5, expression in other clusters ≦ 15%, and Adjusted p value <0.01), we selected the DE-Gs only in cluster 3 and not in other clusters; those up-regulated significantly in cluster 3 were PELI2, ELF3, and EREG (Figure 4C).

Figure 4

Molecular pathways and gene expression pattern in cluster 3

(A) Bar charts showing the statistically significant enriched Reactome pathways in the list of differentially expressed genes (DE-Gs) in cluster 3.

(B) Heatmap showing the clustering analysis of 2982 MM cells sorted from four points, featuring normalized single-cell expression levels of the most variable genes. The blue frame shows the DE-Gs in cluster 3.

(C) Comparison of the expression levels of DE-Gs among subclusters. The rows show the expression of PELI2 (left), ELF3 (middle), and EREG (right).

See also Table S4 and Figure S7.

Importantly, cells expressing PELI2, ELF3, or EREG not only in cluster 3 but also in all the clusters disappeared in response to ixazomib (Figure S7).

Pellino E3 ubiquitin-protein ligase family member 2 as a potential mediator of proteasome inhibitor response in multiple myeloma cells

To investigate whether PELI2, ELF3, and EREG are causally related to the PI response, we established MM cell lines (KMS-34) expressing PELI2, ELF3, or EREG. We detected a substantial amount of mRNAs for PELI2, ELF3, and EREG in KMS-34 cells expressing PELI2, ELF3, or EREG, although no significant amounts of mRNAs were present in mock-transfected KMS-34 cells. In this study, the proliferation of KMS-34 cells expressing PELI2, ELF3, or EREG was roughly comparable to that of mock-transfected cells (Figures 5A and S8). Given that the clone size of cluster 3 was markedly reduced after the addition of ixazomib (Figure 3), we focused on the second line PI, ixazomib from the three PI used.

Figure 5

Cytotoxic assay for KMS-34 cells expressing PELI2, ELF3, or EREG

(A) Bar plot showing PELI2 (left), ELF3 (middle), and EREG (right) mRNA levels.

(B–D) Dose-response curves for KMS-34 cells expressing PELI2, ELF3, or EREG in the presence of Ixa (B), Bor (C), and Len (D).

Results are shown as the mean ± SEM of three independent experiments. p values (extra sum-of-squares F test) of the statistical comparison between IC50 values are indicated. Bor, bortezomib; Ixa, ixazomib; Len, lenalidomide.

See also Figure S8.

In the cytotoxic assay, cells expressing PELI2 were more sensitive to ixazomib (p < 0.0001, extra sum-of-squares F test) with IC50 of 68.03 nM, relative to mock-transfected cells with IC50 of 104.10 nM. On the other hand, forced expression of ELF3 or EREG did not alter the sensitivity to ixazomib (Figure 5B). We also found that those expressing PELI2 were significantly sensitive to another PI, bortezomib (p < 0.0001, extra sum-of-squares F test) with IC50 of 6.46 nM, relative to mock-transfected cells (IC50 of 9.72 nM), confirming the role of PELI2 in the cytotoxic effect of PI (Figure 5C). In this study, forced expression of PELI2 did not alter the sensitivity to lenalidomide (Figure 5D), suggesting the specific effect of PELI2 on the PI response.

Pellino E3 ubiquitin-protein ligase family member 2 can predict the clinical response to proteasome inhibitor treatment

To explore the clinical implication of PELI2, we examined the expression level of PELI2 in patients with MM who were registered in publicly accessible databases (Goldman et al., 2019). The Multiple Myeloma Research Foundation (MMRF)-CoMMpass database (n = 2119) using the UCSC Xena web-based tool revealed that PELI2 was expressed to a varying degree in patients with MM (Figure S9). We next used the data of gene expression and clinical outcome which were analyzed with previously published phase 2 and phase 3 clinical trials of the first line PI, bortezomib, in patients with relapsed MM (Gene Expression Omnibus Accession Number GSE9782) (Mulligan et al., 2007). We only selected patients treated with bortezomib monotherapy from this database and examined whether the expression level of PELI2 is associated with the clinical outcome. The cohort was divided into two groups, based on the median PELI2 mRNA level (Figure 6). In this setting, high expression of PELI2 was associated with the better prognosis than low expression among patients with bortezomib monotherapy (log rank test: p = 0.014; multivariate hazard ratio = 0.672, 95% confidence intervals = 0.459-0.984, p = 0.041). On the other hand, high expression of PELI2 was not associated with a better prognosis than low expression among those treated with dexamethasone monotherapy (log rank test: p = 0.416) (Figure S10).

Figure 6

Kaplan-Meier survival analysis of patients treated with bortezomib monotherapy

Kaplan-Meier plot for patients with high and low PELI2 expression. Patients with MM treated with bortezomib monotherapy (n = 188) are divided into two groups according to the median expression level. The log rank test is used for statistical significance.

See also Figures S9 and S10.

Discussion

Although a variety of medical treatments for MM have been developed, therapy has largely been empiric, based on the mechanisms of action of the drugs. Given the marked heterogeneity of MM (Szalat et al., 2016), it is of great importance to understand how MM cells would clonally evolve during the disease progression and/or treatment course and to predict which clonal cells would be drug-sensitive or -resistant, when anti-myeloma drugs are to be selected. Here, we conducted scRNA-seq analysis of BM cells from a single patient with relapsed and refractory MM, coupled with longitudinal sampling and transcriptomic profiling at the single-cell resolution. The longitudinal sampling strategies herein allowed us to eliminate the effect of inter-patient tumor heterogeneity, so that we characterized the temporal transcriptomic evolution of MM cells, and even understand the drug-induced clonal dynamics.

In this study, longitudinal scRNA-seq analysis revealed dynamic cellular heterogeneity; five transcriptionally distinct cell clusters in a single patient with MM, which is consistent with the notion that MM is composed of heterogeneous cancer cell populations (Jang et al., 2019). Based on the evolutionary trajectory obtained from the pseudotime analysis, we found clonal evolution from cluster 0 to either cluster 1 or cluster 2, during which their clone sizes increased under the therapeutic pressure. These observations suggest that cluster 0 represents drug-resistant precursor MM cells even before treatment, which serves as an origin of resistant clones during treatment with multiple anti-myeloma drugs. On the other hand, treatment with an anti-CD38 monoclonal antibody, daratumumab depleted MM cells expressing CD38 with the expansion of CD38-deficient MM cells, which are thought to be resistant to the anti-CD38 monoclonal antibody, at a single-cell resolution. As the expression level of CD38 is shown to be restored to baseline after daratumumab treatment (Nijhof et al., 2016), it is interesting to know whether the expression of CD38 directly affects the response to daratumumab. Overall, this study helps understanding the underlying mechanism of emerging drug resistance in our patients.

It is noteworthy that cluster 3 occurs at relapse after lenalidomide treatment, and the clone size of cluster 3 declined significantly after the addition of ixazomib. In this study, we found that PELI2 as a DE-Gs in cluster 3 and that forced expression of PELI2 in an MM cell line significantly enhances PI-induced cell death in in vitro cytotoxic assay, suggesting that PELI2 serves as a mediator of PI sensitivity in MM cells. PELI2 is an E3 ubiquitin-protein ligase that catalyzes the cytokine signaling pathway (Jensen and Whitehead, 2003), and is thought to be involved in the activation of NLRP3 inflammasome, thereby inducing caspase-1-mediated cell death or pyroptosis (Humphries et al., 2018). Given that NLRP3 is thought to be degraded by proteasome, PI could contribute to the accumulation of NLRP3 inflammasome (Han et al., 2015). As a forced expression of PELI2 does not affect the proliferation of MM cells, it is likely that PELI2 is capable of inducing cell death only in the presence of PI. Further studies are required to elucidate the underlying mechanism of how increased expression of PELI2 confers the PI-sensitivity in MM cells.

In this study, we found that PELI2 is expressed not only in an MM cell cluster in a single patient but also in MM cells derived from publicly accessible databases, thereby suggesting that PELI2 plays a critical role in the pathogenesis of MM. Importantly, Kaplan-Meier survival analysis of datasets from the GSE9782 database revealed that higher expression of PELI2 is associated with better prognosis in patients with MM, who had been treated with PI. These observations, taken together, suggest that PELI2 has the potential to be used as a predictive biomarker of the PI response in patients with MM.

As a majority of patients with MM relapse during the long-term follow-up, it is important to predict when and how the drug-resistant clones would emerge during the clinical course of the disease and to select which drugs are to be used in real-world clinical setting. Longitudinal scRNA-seq analysis of MM cells provides the clue to identify the drug-resistant clones, which are not reduced even during remission and thereafter are increased. On the other hand, Tirier et al. have recently reported that MM cells alter the phenotype of non-malignant cells in the BM microenvironment, which might contribute to the drug resistance of MM (Tirier et al., 2021). Further studies are required to elucidate the temporal profiles of MM cells, non-MM cells, or even the BM interstitial cells during the treatment course of the disease.

In conclusion, we performed temporal profiling of clinical and transcriptomic features and cellular heterogeneity in a single patient with relapsed and refractory MM through longitudinal scRNA-seq analysis and identified a drug-response gene; PELI2, which confers the PI sensitivity. This study also suggests that PELI2 serves as a predictive biomarker of the PI response in patients with MM. Our integrated strategy combining longitudinal scRNA-seq analysis with in vitro cultured cell assay and database analysis would facilitate the identification of drug-response genes in hematological malignancies.

Limitations of the study

Major limitation of this study was the small sample size. One may argue that PELI2 could not be generalized as a prognostic marker and sensitive gene of MM as we analyze only one patient. However, we thought PELI2 could be generalized, as our integrative approach using a publicly available database revealed that PELI2 was expressed in about half of MM. In addition, previous reports showed that the strong inter-patient heterogeneity of MM might hamper the identification of drug-response genes through the analysis of a large number of patients. Therefore, our method could be a useful tool to identify the drug-response genes of MM. Further investigations on MM using these methods in more patients could lead to their application in the research context. The other limitation is the lack of paired baseline/relapse samples in the published database. Further studies are required to evaluate how PELI2 expression affects the clinical outcome of patients with MM treated with PI.

STAR★Methods

Key resources table

Resource availability

Lead contact

Further information and requests for resources should be directed to the lead contact, Yoshihiro Ogawa (ogawa.yoshihiro.828@m.kyushu-u.ac.jp).

Materials availability

This study did not generate new unique reagents.

Experimental model and subject details

BM cell sampling

The patient was a 67-year-old man with MM at diagnosis. BM cells were aspirated at the time of initial onset and recurrence. Immediately after the aspiration, BM cells were diluted 1:2 in RPMI-1640 (Nacalai Tesque), and mononuclear cells were isolated using density centrifugation media (Lymphoprep) (Axis-Shield, Oslo, Norway). After centrifugation, the mononuclear cells were carefully aspirated and washed with RPMI-1640. The mononuclear cells were stored at −80°C using Cellbanker 1plus (Takara Bio, Shiga, Japan) until use. Written informed consent was obtained from the patient as required by the Declaration of Helsinki. The study was approved by the clinical research ethical committee of Kyushu University (No. 29–494).

Cultured cells

The human MM cell lines, KMS-34, were purchased from the Japanese Collection of Research Bioresources (JCRB, Osaka, Japan). KMS-34 cells were cultured in RPMI-1640 supplemented with 10% FBS and 1% penicillin-streptomycin (Nacalai Tesque, Kyoto, Japan). Cell lines were maintained in a humidified incubator at 37°C with a 5% CO2 atmosphere.

The KMS-34 cells expressing PELI2, ELF3, EREG, and mock-transfected cells were established as follows; cDNAs encoding human PELI2, ELF3, and EREG were separately subcloned into the pMP71CW-IRES-EGFP vectors, a kind gift from Dr. T. Chinen (Kyushu University, Japan). These vectors and empty vector were transfected into KMS-34 cells by lipofection using vesicular stomatitis virus G protein pseudotyped vector particles, respectively.

Method details

Single-cell sorting

The mononuclear cell sample was washed in ice-cold fluorescence-activated cell sorting (FACS) buffer (2% Fetal Bovine Serum (FBS) in HBSS, 4 mM EDTA, pH 8.0). Samples were incubated for 10 min with FcR Blocking Reagent (Miltenyi Biotec, Bergisch Gladbach, Germany), followed by staining for 20 min with antibodies against CD38 (Cytognos, multi-epitope, FITC, Salamanca, Spain), CD45 (BioLegend, Pacific Blue, Clone 2D1, San Diego, CA), CD19 (BioLegend, PE/Cy7, Clone HIB19), and CD138 (BioLegend, APC, Clone DL-101), and filtered through a 35 μm strainer before sorting. Single-cell sorting was performed using FACS Melody (BD Biosciences, San Jose, CA). After doublet exclusion, CD38++CD45−CD19− cells were index-sorted into 384-well plates containing 1.5 μL lysis solution and barcoded poly (T) reverse transcription primers for scRNA-seq. Immediately after sorting, each plate was spun down to ensure cell immersion was in the lysis solution and stored at −80°C.

Single-cell library preparation

Single-cell libraries were prepared as described (Hashimshony et al., 2016). Briefly, mRNA from cells sorted into 384-well plates was reverse-transcribed into cDNA, with which 12-well samples were pooled. The pooled samples were then amplified by T7 in vitro transcription, and the amplified RNA was fragmented and used to construct final libraries by reverse transcription and PCR amplification. The libraries were sequenced on Illumina HiSeq 1500 sequencers (Illumina, San Diego, CA) and the reads were mapped to the GRCh38 human reference genome assembly (release 84) using HISAT2 (version 2.2.1).

Bioinformatic analysis of scRNA-seq data

After generating gene-barcode matrices, we applied Seurat (version 3.2.3) for quality control and downstream analysis. To remove low-quality cells, cells with less than 1000 genes and more than 10% of mitochondrial genes were filtered out. The filtered data was log-normalized using “NormalizeData” function. The scale.factor option was set to 100,000. After the separation of filtered data by experiments, the "FindIntegrationAnchors" and "IntegrateData" functions were used to correct for batch effects and integrate these data. The k.filter option of “FindIntegrationAnchors” function was set to 50 according to the minimum number of cells in the split data. The “CellCycleScoring” function was used to classify the cell-cycle stages, and the heterogeneity of cell cycle scores was then regressed out by “ScaleData” function. The principal component analysis was performed using the "RunPCA" function. Using the data of the top 30 principal components as input, dimension reduction by Uni-form Manifold Approximation and Projection (UMAP) was performed using the "RunUMAP" function. Next, we performed the unsupervised clustering using “FindClusters” function (resolution of 0.4). We used the transcriptome dataset to annotate the cell population through SingleR (version 1.4.1) (Aran et al., 2019). The differential expression analysis was performed using the “FindAllMarkers” and “FindMarkers” function based on the Wilcoxon rank-sum test. The Bonferroni method was used to correct for multiple comparisons.

Trajectory analysis

The preprocessed data were analyzed by CytoTRACE (Cellular Trajectory Reconstruction Analysis using gene Counts and Expression) (Gulati et al., 2020) in R software to predict the relative differentiation state of a single-cell. The "iCytoTRACE" function was used to calculate the CytoTRACE score. Monocle3 (version 0.2.3.0) was used to infer the developmental trajectory and estimate the pseudotime as previously described (Trapnell et al., 2014). In brief, the Seurat object was converted into a Monocle3 object using “as.cell_data_set” function. Then “Preprocess_cds” function was run with top 30 dimensions (PCA) and the resolution parameter of “cluster_cells” function was set to 0.01. The starting point of the inferred trajectory was set to the region where cells with high CytoTRACE scores clustered.

Functional enrichment analysis

The pathway enrichment of DE-Gs was analyzed using the online tool, Metascape (Zhou et al., 2019). With the species limited to “Homo sapiens”, functional enrichment analysis was performed based on Reactome Gene Sets. All genes in the genome have been used as the enrichment background. Terms with a p value <0.01, a minimum count of 3, and an enrichment were collected and grouped into clusters based on their membership similarities.

Quantitative reverse-transcription PCR

Total RNA from transfected KMS-34 cells was extracted with RNeasy Plus Mini Kit (Qiagen, Hilden, Germany), and cDNA was synthesized with PrimeScriptTM RT Master Mix (Takara). Real-time PCR was performed using TB GreenTM Premix Ex TaqTM (Takara). We measured mRNA expression levels of PELI2, ELF3, and EREG using CFX Connect Real-Time PCR detection system (Bio-Rad, Hercules, CA). All PCR data were normalized against expression of glyceraldehyde 3-phosphate dehydrogenase (GAPDH) as described (Gicquel et al., 2015). The primer sequences used in this study are listed in Key resources table.

Cell proliferation and cytotoxic assays

For cell proliferation assay, KMS-34 cells were seeded at a density of 5 × 103 cells/well in 96-well plates. Cell viability was assessed every 24 h over a period of 3 days using CellTiter 96 AQueous One Solution Cell Proliferation Assay kit (Promega, Madison, WI), according to the manufacturer’s instruction.

For cytotoxic assay, KMS-34 cells were seeded at a density of 5 × 104 cells/well in 96-well and incubated with or without ixazomib (MedChemExpress, Monmouth Junction, NJ) or bortezomib (MedChemExpress) for 48 h at 37°C as reported (Cohen et al., 2021). The KMS-34 cells treated with lenalidomide (MedChemExpress) were seeded at a density of 1.25 × 104 cells/well and incubated for 96 h at 37°C. Ixazomib, bortezomib, and lenalidomide were added with a serial dilution. Control wells without cells and containing the same volumes of culture medium and MTS reagent were used as blank to subtract background absorbance. After incubation, 3-(4,5-dimethylthiazol-2-yl)-5-(3-carboxymethoxyphenyl-2-(4-sulfophenyl)-2H-tetrazolium (MTS) reagent (CellTiter 96 AQueous One Solution Cell Proliferation Assay; Promega) was added and the cells were incubated for 2 h. The absorbance was measured at a wavelength of 490 nm using an EnSpire Multimode Plate Reader (PerkinElmer, Waltham, MA) and expressed as the percentage of the values of untreated cells. MTS assay was used to establish dose-response inhibitory (variable slope) curve and determine the IC50 values of each drug. The IC50 value were calculated by non-liner regression analysis with GraphPad Prism software (GraphPad Software, San Diego, CA). The values are the average of three independent experiments. The extra sum-of-squares F test was used to test whether the IC50 values differ between the groups.

Evaluation of expression level of PELI2 in MM patients using the public database

The UCSC Xena is a functional genomics browser that provides visualization and integration for analyzing and viewing the public data hubs. Stacked area chart was generated by data mining in MMRF-CoMMpass database using the UCSC Xena browser.

Evaluation of PELI2 as a prognostic marker in clinical setting

Gene expression data and clinical characteristics of patients with relapsed MM were extracted from the GSE9782 dataset. We only selected 188 patients treated with bortezomib monotherapy from the database. For Kaplan-Meier plotter analysis, the patients were divided into two groups, based on the median PELI2 mRNA level. We examined differences in overall survival using log rank test and multivariate Cox regression analysis. Age, sex, race, immunoglobulin type and number of prior lines were included in the multivariate analysis. We also selected 76 patients treated with dexamethasone monotherapy and compared the Kaplan-Meier curves using log rank test.

Quantification and statistical analysis

Statistical analyses were performed using Graphpad Prism and R software. Data are expressed as the mean ± SEM of three independent experiments performed in duplicate, and analyzed using the Student’s t-test or Wilcoxon rank-sum test. The Kaplan-Meier method was used with log rank test to compare overall survival between groups. The Cox regression model was used to calculate hazard ratio. p < 0.05 was considered statistically significant.

REAGENT or RESOURCE	SOURCE	IDENTIFIER
Antibodies
CD38-multi-epitope-FITC Kit	Cytognos	Cat#CYT-38F2; RRID: AB_2828013Lot#1911689
Pacific Blue(TM) anti-human CD45 (Clone 2D1)	BioLegend	Cat#368539; RRID: AB_2716029
PE/Cyanine7 anti-human CD19 (Clone HIB19)	BioLegend	Cat#302215; RRID: AB_314245
APC anti-human CD138 (Syndecan-1) (clone DL-101)	BioLegend	Cat#352307; RRID: AB_10901175
Pacific Blue Mouse IgG1, κ Iso-type Ctrl Antibody(Clone MOPC-21)	BioLegend	Cat#400131
PE/Cyanine7 Mouse IgG1, κ Iso-type Ctrl Antibody(Clone MOPC-21)	BioLegend	Cat#400125; RRID: AB_2861433
APC Mouse IgG1, κ Iso-type Ctrl Antibody(Clone MOPC-21)	BioLegend	Cat#400119; RRID: AB_2888687
Bacterial and virus strains
NEB Stable Competent E.coli (High Efficiency)	NEB	C3040
Biological samples
Bone marrow samples	This paper	N/A
Chemicals, peptides, and recombinant proteins
Lymphoprep	Axis-Shield	Cat#1856
Cellbanker 1plus	Takara Bio	Cat#CB021
FcR Blocking Reagent	Miltenyi Biotec	Cat#130-059-901
Propidium Iodide Solution	BioLegend	Cat#421301
ERCC RNA Spike-In Mix	Invitrogen	Cat#4456740
Nonidet P40 Substitute solution	Sigma-Aldrich	Cat#98379-10ML-F
Deoxynucleotide (dNTP) Solution Mix	NEB	Cat#N0447
SuperScript Ⅱ Reverse Transcriptase	Invitrogen	Cat#18064022
RNaseOUT Recombinant RNase Inhibitor	Invitrogen	Cat# 10777019
Second Strand Buffer	Invitrogen	Cat#10812014
DNA Polymerase I	Invitrogen	Cat#18010025
E. coli DNA Ligase	Invitrogen	Cat#18052019
RNase H	Invitrogen	Cat#18021071
AMPure XP	Beckman Coulter	Cat#A63880
ExoSAP-IT PCR Product Cleanup Reagent	Applied Biosystems	Cat#78201.1.ML
RNAClean XP	Beckman Coulter	Cat#A63987
Phusion High-Fidelity PCR Master Mix with HF Buffer	NEB	Cat#M0531
RNase A Solution (10 mg/mL)	Nacalai tesque	Cat#30100-31
E-Gel SizeSelect II Agarose Gels, 2%	Invitrogen	Cat#G661012
PhiX Control v3	Illumine	Cat#15017666
Fetal Bovine Serum, qualified, Brazil	Gibco	Cat#10270106
KOD One PCR Master Mix(Dye-free 2xPCR Master Mix)	TOYOBO	Cat#KMM-201
EcoRⅠ	Takara Bio	Cat#1040A
SalⅠ	Takara Bio	Cat#1080A
PEI MAX	Polysciences	Cat#24765-1
DOTAP Liposomal Transfection Reagent	ROCHE	Cat#11202375001
PrimeScript RT Master Mix (Perfect Real Time)	Takara Bio	Cat# RR036
TB Green Premix Ex Taq II (Tli RNaseH Plus)	Takara Bio	Cat# RR820
Ixazomib	MedChemExpress	Cat#HY-10453
Bortezomib	MedChemExpress	Cat#HY-10227
Lenalidomide	MedChemExpress	Cat#HY-A0003
CellTiter 96 AQueous One Solution Cell Proliferation Assay (MTS)	Promega	Cat#G3580
Critical commercial assays
MEGAscript T7 Transcription Kit	Invitrogen	Cat#AM1334
MinElute PCR Purification Kit	QIAGEN	Cat#28006
High Sensitivity DNA kit	Agilent	Cat#5067-4626
NEBNext Library Quant Kit for Illumina	NEB	Cat#E7630
HiSeq PE Rapid Cluster Kit v2	illumina	Cat#PE-402-4002
HiSeq Rapid SBS Kit v2 (50 cycles)	illumina	Cat#FC-402-4022
MinElute Gel Extraction Kit	QIAGEN	Cat#28604
DNA Ligation Kit < Mighty Mix>	Takara Bio	Cat#6023
QIAprep Spin Miniprep Kit	QIAGEN	Cat#27106
RNeasy Plus Mini Kit	QIAGEN	Cat#74136
Deposited data
Raw and analyzed data	This paper	GEO: GSE193695
Mendeley data	This paper	https://doi.org/10.17632/wxyb4wt7j2.1
Human referece genome NCBI GRCh38	Genome Reference Consortium	https://www.ncbi.nlm.nih.gov/grc/human
Experimental models: Cell lines
Human: KMS-34	JCRB	JCRB1195; RRID: CVCL_2996Lot#06052015
Oligonucleotides
PELI2 Forward primer: 5′-GAACACTGCGCCCCCAATAA-3′	This paper	N/A
PELI2 Reverse primer: 5′-AAGCACCATTGTACCCGAGC-3′	This paper	N/A
ELF3 Forward primer: 5′-CACTGATGGCAAGCTCTTC-3′	This paper	N/A
ELF3 Reverse primer: 5′-GGAGCGCAGGAACTTGAAG-3′	This paper	N/A
EREG Forward primer: 5′-TGCTCTGCCTGGGTTTCCAT-3′	This paper	N/A
EREG Reverse primer: 5′-ACTGGACTCTCCTGGGATACA-3′	This paper	N/A
GAPDH Forward primer: 5′-CCTGGTCACCAGGGCTGC-3′	This paper	N/A
GAPDH Reverse primer: 5′-GCCAGTGGACTCCACGAC-3′	This paper	N/A
Recombinant DNA
pMP71CW-IRES-EGFP vectors	T. Chinen (Kyushu University, Japan)	N/A
pCMV-VSV-G	addgene	Cat#8454; RRID: Addgene_8454
Software and algorithms
Graphpad Prism 8	GraphPad	https://www.graphpad.com/; RRID: SCR_002798
FlowJo v10	Tree Star	https://www.flowjo.com/; RRID: SCR_008520
R Project for Statistical Computing v4.0.3	R software	https://www.r-project.org/; RRID: SCR_001905
Seurat v3.2.3 R package	Stuart et al. (2019)	https://satijalab.org/seurat/; RRID: SCR_016341
SingleR v1.4.1 R package	Aran et al. (2019)	https://github.com/dviraran/SingleR
CytoTRACE v0.3.3 R package	Gulati et al. (2020)	https://cytotrace.stanford.edu/
monocle3 v0.2.3.0 R package	Trapnell et al. (2014)	https://cole-trapnell-lab.github.io/monocle3/; RRID: SCR_018685
Metascape	Metascape	https://metascape.org/gp/index.html; RRID: SCR_016620
UCSC Xena	UCSC Xena	https://xena.ucsc.edu/; RRID: SCR_018938