Time-Resolved scRNA-Seq Tracks the Adaptation of a Sensitive MCL Cell Line to Ibrutinib Treatment

Abstract

Since the approval of ibrutinib for relapsed/refractory mantle cell lymphoma (MCL), the treatment of this rare mature B-cell neoplasm has taken a great leap forward. Despite promising efficacy of the Bruton tyrosine kinase inhibitor, resistance arises inevitably and the underlying mechanisms remain to be elucidated. Here, we aimed to decipher the response of a sensitive MCL cell line treated with ibrutinib using time-resolved single-cell RNA sequencing. The analysis uncovered five subpopulations and their individual responses to the treatment. The effects on the B cell receptor pathway, cell cycle, surface antigen expression, and metabolism were revealed by the computational analysis and were validated by molecular biological methods. The observed upregulation of B cell receptor signaling, crosstalk with the microenvironment, upregulation of CD52, and metabolic reprogramming towards dependence on oxidative phosphorylation favor resistance to ibrutinib treatment. Targeting these cellular responses provide new therapy options in MCL.

Keywords: drug resistance; ibrutinib; mantle cell lymphoma; scRNA-seq.

Figure 1

Overview of scRNA-seq approach and evaluation of reproducibility. (A) REC-1 cells were cultured for 48 h in total, one sample stayed untreated (Ctr), the second was treated for 6 h and the third for 48 h with ibrutinib (400 nM). Scatter plot (SSC-A vs. fluorescence intensity (propidium iodide)) of 48 h is shown representative for all treatments, cells of gate P1 (red) entered scRNA-seq (relative proportion of included cells are indicated for each sample). Droplet-based single-cell RNA sequencing was performed using the 10x Genomics platform with encapsulation of single cells in GEMs (gel bead-in-emulsions). The data was analyzed implementing Cell Ranger (10x Genomics), and R packages such as Seurat, clusterProfiler, and SCENIC; (B) uniform manifold approximation and projection (UMAP) representation visualized the five subpopulations in the aggregated (replicate 1 and 2) and cell cycle regressed data set of Ctr and the corresponding dot plot shows the expression of selected top 10 marker genes for every subpopulation; size of dots refers to the percentage of cells expressing the gene, color intensity represents the average expression level.

Figure 2

Heterogeneity of REC-1 cell line on the single-cell level. (A) Uniform manifold approximation and projection (UMAP) representation visualizing 7 clusters at resolution 0.4 in the integrated and cell cycle regressed data set of untreated cells (Ctr) and the corresponding heatmap showing selected top 10 marker genes (see Table S2 for cluster markers); (B) gene set enrichment analysis for Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways and gene ontologies (biological processes) including genes of Ctr with a log fold change > 0.1 and adjusted p-value < 0.001 (no. = number); no matches were detected for cluster 6 due to few differentially expressed genes (DEGs); (C) violin plot of CD52 expression in Ctr (clustering is shown in (A); and (D) bar plot showing the proportions of predicted cell cycle phases (G1, G2/M, or S) across clusters of Ctr by Seurat’s cell cycle scoring.

Figure 3

Evolution of subpopulations across ibrutinib treatment. (A) Uniform manifold approximation and projection (UMAP) representations of the combined data set including replicate 1 and 2 of Ctr, 6 h, and 48 h after cell cycle regression (clusters are shown for resolution 0.4 with additional subclustering of subpopulation D at resolution 0.3); (B) proportions of subpopulations across treatment; (C) heatmap of selected marker genes of Ctr, 6 h, and 48 h (see Table S3 for treatment markers); (D) distribution of predicted cell cycle phases (G1, G2/M, and S) across treatment by Seurat’s cell cycle scoring; (E) distribution of cell cycle phases acquired by flow cytometry (sensitive REC-1 compared to resistant MAVER-1 cell line, D = DMSO control, 48 h = 48 h 400 nM ibrutinib treatment; n = 3, * p ≤ 0.05, ns = not significant); (F) violin plot of LDHA expression levels of the single-cell sequencing data; (G) Western blot of LDHA expression after ibrutinib treatment (DMSO as control, 400 nM ibrutinib for 2 d, 3 d, and 4 d) in sensitive REC-1 and in resistant MAVER-1, β-actin served as loading control, relative expression (Rel. Expr.) to DMSO control was calculated after normalization to β-actin (Western blot and relative expression values are shown representative for three independent replicates); and (H) extracellular flux analysis of 3 d ibrutinib (400 nM) or DMSO (control) treated cells (sensitive REC-1 compared to resistant MAVER-1) by Agilent Seahorse XF 96 Analyzer; the ratio of oxygen consumption rate (OCR) to extracellular acidification rate (ECAR) is shown (n = 3, ** p ≤ 0.005, ns = not significant).

Figure 4

Altered gene regulatory networks during ibrutinib treatment. (A) Heatmap representing the alterations in gene regulatory networks in the eleven clusters of the combined analysis (Ctr, 6 h, and 48 h) (clustering is shown in Figure 3A), numbers in brackets indicate the amount of genes forming the gene regulatory networks of the indicated transcription factor and (B) uniform manifold approximation and projection (UMAP) representations in the first row show the regulon activity of the indicated transcription factors, UMAPs underneath display the gene expression (log normalized) of the transcription factors.