Identification of a clinically efficacious CAR T cell subset in diffuse large B cell lymphoma by dynamic multidimensional single-cell profiling

Abstract

Chimeric antigen receptor (CAR) T cells used for the treatment of B cell malignancies can identify T cell subsets with superior clinical activity. Here, using infusion products of individuals with large B cell lymphoma, we integrated functional profiling using timelapse imaging microscopy in nanowell grids with subcellular profiling and single-cell RNA sequencing to identify a signature of multifunctional CD8⁺ T cells (CD8-fit T cells). CD8-fit T cells are capable of migration and serial killing and harbor balanced mitochondrial and lysosomal volumes. Using independent datasets, we validate that CD8-fit T cells (1) are present premanufacture and are associated with clinical responses in individuals treated with axicabtagene ciloleucel, (2) longitudinally persist in individuals after treatment with CAR T cells and (3) are tumor migrating cytolytic cells capable of intratumoral expansion in solid tumors. Our results demonstrate the power of multimodal integration of single-cell functional assessments for the discovery and application of CD8-fit T cells as a T cell subset with optimal fitness in cell therapy.

Extended Data 1.. Phenotype characteristics of infusion products measured by flow cytometry.

(A) Comparisons of CAR+ T cells recorded by flow cytometry for all sixteen patients (n). There is no significant difference in the CAR frequency between CR and PR/PD. (B) Comparisons of CD4+ T cells recorded by flow cytometry for all sixteen patients (n). There is no significant difference in the CD4 frequency between CR and PR/PD.

Extended Data 2.. T cells from CR are enriched for serial killing, increasing mitochondrial and lysosomal size, and persistent migration.

(A) Schematic of a killing event at an E:T of 1:1 in which a CAR T-cell conjugates and kills a NALM-6 cell. The plot on the right shows the killing rate comparison between T cells from CR and PR/PD within all 1E:1T nanowells. Each dot (n) represents the frequency of killer T cells for each IP. Micrograph showing an example of 1E:1T killing event through the 6-hours (hh:mm) of time-lapse imaging from a CR IP. (B) Schematic of a killing event, showing the interaction parameters. tSeek defined as the time for CAR-T cell to find and conjugate to the NALM-6 cell. tconjugation is defined as the duration of CAR-T cell in stable conjugation with NALM-6 cell. tDeath is the time interval between the start of the conjugation and the apoptosis of the NALM-6 cell. Plots show the comparison between T cells from CR and PD/PR for these parameters. Each dot (n) represents the average value for all T cells within each IP. (C) Representative violin plot shows the interaction parameters from one responder IP: tSeek (n=96 events), tconjugation (n=94 events), tDeath (n=30 events). The bar graph shows the killing frequencies of the same responder IP at an E:T of both 1:1 and 1:2. (D) Schematics and examples of serial killing, mono killing and no-killing events in nanowells with an E:T of 1:2. (E) Unsupervised hierarchical clustering based on parameters from TIMING, and confocal microscopy. Serial killing, migration, increasing mitochondrial and lysosomal volume were features associated with T cells from CR patients. (F) Cytotoxicity and motility correlations with CD4/CD8 ratio. Each point represents the average parameter for each IP (n=9 patient IPs). Cytotoxicity is defined as the frequency of 1E:1T killing from TIMING. The Pearson’s correlation coefficient was calculated for CR and PR/PD IP.

Extended Data 3.. T-cell phenotypes defined using scRNA-seq.

(A) Uniform Manifold Approximation and Projection (UMAP) for 21,469 cells from nine IPs. Bar graph showing the distribution of T cells from CR and PD/PR among 10 clusters determined using unsupervised clustering. (B) Bubble plot showing key genes associated with T-cell migration and exhaustion phenotypes. (C) Pseudotime trajectory analysis for two clusters enriched in PD (CD8–1: n=572 cells, CD8–2: n=1,031 cells) and one cluster enriched in CR (CD8–6: n= 1,253 cells). Necklace plots show CD8–2 (Central Memory dominant) differentiate into Effector/Effector Memory dominant CD8-fit (CD8–6) and CD8–1. (D) Violin plot (left) showing the ssGSEA score for AMPK activation calculated for two clusters enriched in PD (CD8–1: n=572 cells, CD8–2: n=1,031 cells) and CD8-fit (CD8–6: n= 1,253 cells) cluster enriched in CR. Violin plot (right) showing the ssGSEA score for the TCF7 signature for the three clusters. The black bar represents the median and the dotted lines denote quartiles. P values were computed using two-tailed Welch’s T-test. (E) Schematic overview of the experimental process used to identify signatures of persistent CAR T cells. Single-cell gene expression and T-cell receptor (TCR) datasets were generated by sequencing pre- (GMP: good manufacturing practice facility) and post-infusion CD19 CAR T cells from blood and bone marrow samples of pediatric patients with B-ALL. (F) Violin plots showing the transcriptome similarities between the CD8⁺ T cells from our datasets and CD8⁺ GMP effector precursors. P values were calculated using two-tailed t-test. The black bar represents the median and the dotted lines denote quartiles. P values were computed using two-tailed Welch’s T-test. (G) Gene set enrichment analysis (GSEA) of CD8⁺ GMP effector precursors gene signatures within cells from CD8-fit cluster compared with cells from all other CD8 clusters. The effector precursors gene signature is based on differentially expressed genes between the CD8⁺ effector precursors clusters and all other GMP CD8⁺ T cells from part E.

Extended Data 4.. Matrix binding genes are significantly upregulated in the CD8-fit population.

Violin plots showing the expression of matrix binding genes enriched in CD8-fit cluster (CD8–1: n=572 cells, CD8–2: n=1,031 cells, CD8-fit: n=1,253 cells). The black bar represents the median and the dotted lines denote quartiles. P values were computed using two-tailed Welch’s T-test.

Extended Data 5.. CD8-fit cells can be identified in healthy donor derived T cells and in the premanufacture PBMCs of patients treated with CAR T cells.

(A) Overview of external dataset study GSE201035. (B) Uniform Manifold Approximation and Projection (UMAP) for 6,713 cells from two donors. Bar graph showing the distribution of CD8⁺ T cells among 8 clusters determined using unsupervised clustering. (C) Violin plot showing the CD8-fit ssGSEA score comparison between the 8 healthy donor CD8+ T-cell clusters. For the violin plot, the black bar represents the median and the dotted lines denote quartiles. P values were computed using one-way ANOVA with Holm-Šídák’s multiple comparisons test. (D) Validation of the association between CD8-fit and clinical responses in pre-manufactured T cells. Single-cell gene expression datasets were generated by sequencing pre-manufactured T cells from patients with B-cell lymphoma. ssGSEA-derived migration scores between CD8⁺ T cells from CR (n=13,930) and PD (n=3,679) were computed. For the violin plot, the black bar represents the median and the dotted lines denote quartiles. P value was computed using two-tailed Welch’s T-test.

Extended Data 6.. CD4⁺T-cell phenotypes defined using scRNA-seq.

(A) Comparisons of the T-cell migration scores between all CD4⁺ T cells (n=12,527) from 9 CR and PD/PR IPs. The black bar represents the median and the dotted lines denote quartiles. P values were computed using two-tailed Welch’s T-test. (B) UMAP for CD4⁺ T cells (n=12,527). Nine clusters were identified using unsupervised clustering. (C) Heat map of two CD4⁺ T-cell clusters generated by unsupervised clustering. CD4–4 mostly cells from PR/PD while CD4–1 are enriched with CR cells. A color-coded track on top shows the cells from infusion products of CR (green) and PR/PD (red). The track below the heatmap, shows the sample origin for each cell. (D) Bubble plot showing key genes differentially expressed among CD4⁺ T clusters. P value was calculated using the Wilcoxon rank sum test with Bonferroni correction.

Extended Data 7.. The impact of AMPK inhibition on T cell antitumor function revealed by TIMING.

(A/B) The migration and polarization of 19–28z T cells treated with Compound C (CC). All data representative of three independent experiments with 19–28z T cells from three healthy human donors at an E:T of 1:1. The black bar represents the median and the dotted lines denote quartiles. The P value was computed using a two-tailed Welch’s t-test. (C) Comparisons of the killing frequency of vehicle treated (DMSO) or CC treated 19–28z CAR T cells. Each data point represents a single cell. P value was computed using two-tailed log-rank test. (D) The cumulative frequency of T cells conjugating to tumors cells over 8 hours. P value was computed using two-tailed log-rank test. (E) The duration of conjugation between individual T cells and tumor cells. Each data point represents a single cell. The black bar represents the median and the dotted lines denote quartiles. The P value was computed using a two-tailed Welch’s t-test. (F) The impact of AMPK inhibition on migrating capabilities of expanded tumor infiltrating lymphocytes (TILs) from melanoma patients. The migration speed of individual TILs as measured by TIMING (1E:0T). Each data point represents average migration for a single cell. The black bar represents the median and the dotted lines denote quartiles. P value was computed using two-tailed Welch’s t-test.

Extended Data 8.. Enrichment and functional characterization of migratory 19–28z T cells.

(A) Comparisons between the migration of migrated (migratory) and unsorted cells as measured by TIMING (1E:0T). The black bar represents the median and the dotted lines denote quartiles. P value was computed using two-tailed Welch’s t-test. (B) The frequency of conjugation of T cells to tumor cells comparing migratory and unsorted 19–28z T cells evaluated using TIMING (1E:1T). P value was computed using a two-tailed chi-square test. (C) Killing percentage of migratory and unsorted 19–28z T cells evaluated using TIMING (1E:1T). P value was computed using a two-tailed chi-square test. (D) Phenotyping of matched 19–28z and migratory 19–28z T-cell populations. This data is representative of at least four healthy donor-derived T-cell populations measured by flow cytometry. Gating strategies are shown from the side scatter and forward scatter panels on the far left.

Extended Data 9.. Migratory 19–28z T cells showed enhanced antitumor activity compared to the unsorted 19–28z T cells in suboptimal dose model.

(A) Design of mice experiments to determine the relative efficacy of the 19–28z populations. Mice were treated with either 19–28z T cells or migratory 19–28z T cells five days after the injection of ffLuc expressing EGFP+NALM-6 cells. The mice were euthanized on day 31 and five mice from each of the two T-cell treated groups was used to quantify both the tumor cells and persisting T cells within the spleen and bone marrow of mice. (B) False-colored images illustrating the photon flux from ffLuc expressing EGFP⁺NALM-6 cells treated with suboptimal doses of 19–28z T cells. (C) Time course of the longitudinal measurements of NALM-6 derived photon flux from the three separate cohorts of mice (n= 5 in each group). The dotted line marks the day (Day 5) where 19–28z T cells were introduced into the mice. The background luminescence was defined based on mice with no tumor. Error bars represent SEM and P values were computed using a two-tailed Mann-Whitney t-test.

Extended Data 10.. Quantifying the link between migration and functionality in diverse CARs.

(A-C) The polarization (reflective of the morphology of migratory cells) of individual killer and non-killer CAR T cells without and with conjugation to tumor cells. All data from an E:T of 1:1. (A) shows the data for 19–8-28z T cells tested against NALM-6 cells. (B) and (C) show the data for two different constructs of tri-specific CAR⁺ T cells tested against patient-derived tumor cells. The black bar represents the median and the dotted lines denote quartiles. All P values were computed using two-tailed t-tests and each data point represents a single effector cell.

Figure 1.. Study design for integrated single-cell multi-omic profiling of patients’ infusion products.

Overview of the experimental design for profiling the residual CAR T-cell infusion products of 16 DLBCL patients (9 CR, 7 PR/PD). Cells were used for scRNA-seq analysis, confocal microscopy, and Timelapse Imaging Microscopy In Nanowell Grids (TIMING).

Figure 2.. T cells from CR patients were enriched for migration, serial killing, and mitochondrial volume, compared to T cells from PR/PD.

(A) Schematic of a serial killing event wherein a CAR T-cell conjugates and kills two NALM-6 cells. The plot shows the comparison of serial killing by T cells from either CR (n=9 patient IPs) or PR/PD (n=7 patient IPs) within all 1E:2T nanowells. Each dot (n) represents the frequency of serial killing for a single IP. Micrograph showing a serial killing event through the 6-hours (hh:mm) of timelapse imaging. The dotted line represents the median in the violin plot and the p-value was computed using a two-tailed Mann-Whitney test. (B) Schematic of the migration for a CAR T-cell On the left plot, each dot (n) represents the average T-cell migration for each IP (CR: n=9, PR/PD: n=7 patient IPs) [1E:0T], while the plot on the right shows the comparison between all T cells by response. Micrograph showing a T-cell with high (2 μm/min) and low migratory capacity (0.2 μm/min).* (C) Schematic of a CAR T-cell migrating with/without conjugation to a NALM-6 cell (1E:1T). Plot showing the comparison between the migration of CR and PR/PD CAR T cells within all 1E:1T nanowells. Micrograph showing a CAR T-cell migrating before/during conjugation. * (D) Schematic of a CAR T-cell migrating before conjugation with NALM-6 cell (1E:1T) in killer and non-killer CAR T cells. The plot on the right shows the comparison between migration of killer and non-killer CAR T cells within all 1E:1T nanowells. Micrographs showing killer and non-killer CAR T cells. * (E) Comparisons of the size of mitochondria and lysosome between T cells derived from either CR or PR/PD IP. The confocal 3D image illustrates a representative example of the nucleus (blue), total mitochondria(green), and lysosome (red). * (F) Plots showing the correlation between average organelle size and average migration (1E:1T, without conjugation) of T cells. Each dot represents the average value for T cells from a single IP (CR: n=9, PR/PD: n=7 patient IPs). The P value and Pearson correlation were calculated for linear regression. Error bars represent SEM. *The black bar represents the median, and the dotted lines denote quartiles in violin plots. P values were computed using two-tailed Welch’s t-test.

Figure 3.. Molecular profile of CD8-fit, a subset of CAR T infusion products, associated with clinical response revealed by single-cell RNA sequencing.

(A) UMAP for CD8⁺ T cells (n=7,439 cells from 9 patient IPs). Seven clusters were identified using unsupervised clustering. (B) Heat map of three CD8⁺ T-cell clusters generated by unsupervised clustering. A color-coded track on top shows the cells from infusion products of CR (green) and PD (red), followed by tracks showing ssGSEA scores of TCR, actin cytoskeleton regulation, PGC1A, and RHO pathways, respectively. B: BIOCARTA and K: KEGG show the source for the pathways. The track below the heatmap shows the sample origin for each cell. p values were computed using a two-tailed Wilcoxon rank sum test with Bonferroni correction. (C) Bubble plot showing key genes differentially expressed among the three CD8⁺ T clusters. P values were computed using a two-tailed Wilcoxon rank sum test with Bonferroni correction. (D) Violin plots showing the ssGSEA score comparison between three CD8⁺ T-cell clusters (CD8–1: n=572 cells, CD8–2: n=1,031 cells, CD8-fit: n=1,253 cells) with the signatures of killer and multifunctional serial-killer T cells. * (E) Validation of the association between CD8-fit and clinical responses in independent datasets. ssGSEA-derived CD8-fit scores between CD8⁺ T cells from CR and PD were computed. * (F) Comparison of the migration scores between CD8⁺ T cells from CR and PD.* (G) Violin plot showing the ssGSEA score for migration calculated for two clusters enriched in PD (CD8–1 and CD8–2) and CD8-fit cluster enriched in CR.* (H) Violin plots comparing the similarity between CD8-fit and signatures of T cells that traffic to the bone marrow and spleen from reference dataset of T cells from 12 organ donors. (I) Schematic overview showing the parameters defined by STARTRAC analysis. (J-L) Molecular signatures of CD8-fit overlap with the signatures of tumor resident T cells identified by STARTRAC analysis in three separate studies (CD8–1: n=572 cells, CD8–2: n=1,031 cells, CD8-fit: n=1,253 cells) [Supplementary Table 4]. The title of each violin plot corresponds to the cluster with its functional annotation as defined by STARTRAC analysis in the studies. * *For violin plots, the black bar represents the median and the dotted lines denote quartiles. P values were computed using two-tailed Welch’s test.

Figure 4.. Marker free transwell migration assay enables the enrichment of CD8-fit cells.

(A) Modified transwell assay for the isolation of migratory and unmigratory 19–28z T cells. (B) Percentages of 19–28z T cells that migrated across transwell membrane after 4–6 hours. The individual points (n=5) are the percentages for separate transwell experiments from four separate donors (biological replicates) with the error bars shown as SEM. (C) Uniform Manifold Approximation and Projection (UMAP) for 16,158 cells from technical replicate experiments for two separate donors. Bar graph showing the distribution of CD8⁺ T cells from migratory and non-migratory among 12 clusters determined using unsupervised clustering. (D) Bubble plot showing key genes consistent with CD8-fit (Figure 2C) among four CD8⁺ T clusters. P values were computed using a two-tailed Wilcoxon rank sum test with Bonferroni correction. (E) Violin plots showing the ssGSEA score comparison between migratory and non-migratory CD8⁺ T-cell clusters with the signatures of CD8-fit, killer, T-cell migration, and PGC1a regulation. The black bar represents the median and the dotted lines denote quartiles. P values were computed using two-tailed Welch’s test (F) Basal OCR levels were measured for migratory and unmigratory 19–28z T-cell populations. P values were calculated at each timepoint (n=4 technical replicates) comparing the migratory and non-migratory subsets using two-way ANOVA with Bonferroni correction. **** represents P values less than 1e-15. (G) A confocal 3D image of a migratory 19–28z T-cell. Nuclei are shown in blue and mitochondria in green. The plot shows the number of mitochondria per cell compared between migratory and non-migratory 19–28z T cells.

Figure 5.. The phenotype and in vivo efficacy of migratory CAR T cells.

(A/B) The phenotype (A) and Granzyme B (GzB) expression (B) of the migratory and unsorted 19–28z T cells as determined by flow cytometry. (C) Representative false-colored images illustrating the photon flux from ffLuc expressing EGFP⁺NALM-6 cells. (D) Time course of the longitudinal measurements of NALM-6 derived photon flux from the three separate cohorts of mice (n= 10 mice in each group). Representative data from two repeats is shown. The background luminescence was defined based on mice with no tumor. Error bars represent SEM and p values are computed using the two-tailed Mann-Whitney test. (E) On day 31, four mice from each group were euthanized, and tissues (bone marrow and spleen) were harvested and analyzed by flow cytometry for expression of human CD3 (human T cells) and EGFP (gated within hCD19 cells). The flow data is representative from one mouse in each group.

Figure 6.. Quantifying the link between migration and functionality in diverse CARs.

(A, C and E) Schematic illustrating the CAR structure, manufacturing and expansion, and the target cells used for profiling functionality of individual CAR T cells using TIMING. (B, D and F) The migration of individual killer and non-killer CAR T cells with or without conjugation to tumor cells. All data from an E:T of 1:1. The black bar represents the median and the dotted lines denote quartiles. All P values were computed using two-tailed Mann-Whitney tests and each data point (n) represents a single effector cell.