Uncovering the mode of action of engineered T cells in patient cancer organoids

Extending the success of cellular immunotherapies against blood cancers to the realm of solid tumors will require improved in vitro models that reveal therapeutic modes of action at the molecular level. Here we describe a system, called BEHAV3D, developed to study the dynamic interactions of immune cells and patient cancer organoids by means of imaging and transcriptomics. We apply BEHAV3D to live-track >150,000 engineered T cells cultured with patient-derived, solid-tumor organoids, identifying a ‘super engager’ behavioral cluster comprising T cells with potent serial killing capacity. Among other T cell concepts we also study cancer metabolome-sensing engineered T cells (TEGs) and detect behavior-specific gene signatures that include a group of 27 genes with no previously described T cell function that are expressed by super engager killer TEGs. We further show that type I interferon can prime resistant organoids for TEG-mediated killing. BEHAV3D is a promising tool for the characterization of behavioral-phenotypic heterogeneity of cellular immunotherapies and may support the optimization of personalized solid-tumor-targeting cell therapies.


Uncovering the mode of action of engineered T cells in patient cancer organoids
Dekkers, Alieva et al.

Supplementary Discussion
By applying BEHAV3D, we show that, through unbiased clustering of dynamic imaging features, differences in behavior between different engineered T cell products can be uncovered, as well as changes in behavior induced by the type of PDO.Using this pipeline, we report on the broad targeting potential of TEGs for BC, poorly permissive to current immunotherapies 61 , thereby, providing evidence in favor of clinical potential against solid tumors, albeit with various responsiveness between individual donors.In addition, through behavior-guided transcriptomics we have generated, to our knowledge, the first molecular map underlying the behavioral landscape of immune cells targeted to solid tumors.By exploiting these results, we were able to design an optimal sequence of IFN-I and TEG combination therapy to boost TEG BC PDO targeting.
Different from recent studies that have mapped the activation trajectories of murine immune cells during viral infection 62 , or human immune cells in normal physiology or cancer 63 , we here reconstructed activation trajectories for engineered T cells and uniquely exploited dynamic imaging data revealing their single-cell behavior.This allowed us to dissect gene programs induced by environmental stimuli, versus induction by short or prolonged tumor engagement, and thereby identify the gene signature of TEGs that (serially) killed tumor cells.
A preferred behavioral feature of T cell therapies 31-33 that we also observed in WT1 T cells and ROR1 CAR T cells.This (serial) killing signature includes genes not previously linked to T cell function, thereby providing opportunities to potentially engineer next generation T cells with potent serial killing capability.Furthermore, multiple genes in this signature are associated with morphological plasticity.Such plasticity may underlie the remarkable cellular extensions of serial-killing TEGs observed in our 3D imaging data.Using these protrusions, TEGs intercalated between tumor cells while sequentially killing multiple tumor cells in the PDO, suggesting that morphological plasticity may be an important attribute in the targeting of solid tumors.In support of this, we observed similar morphological plasticity also in WT1 T cells and ROR1 CAR T cells targeting solid tumor PDOs, which rely on a different mode-ofrecognition (tumor-specific antigen), compared to metabolome-sensing TEGs.
By linking T cell behavior to population phenotypes, we were able to show that profound tumor targeting, including serial killing, was a predominant feature of CD8 + TEGs, making this subset the most attractive effector population for treatment.Yet, CD4 + T cells, for which we detected a specific behavioral signature characterized by high movement and short organoid engagement, are thought to display an indispensable role in supporting the cytotoxic function and persistence of CD8 + engineered T cells 64,65 .Thus, the preferred product for clinical treatment, will depend on an optimal combination between CD4 + and CD8 + subsets.
Moreover, we show that tumor-targeting by CD8 + TEGs can be further enhanced by sorting NCAM1 + cells, thereby providing proof-of-concept that the efficacy of an engineered T cell population can be improved through cell selection.Similar results could potentially be achieved through functional skewing or combination therapy, as we observed higher proportions of NCAM1 + cells in TEGs expanded in the presence of IL-15, and recombinant IL-15 is already used in a human clinical trial as a combination partner to natural killer cell immunotherapy 66 .
Type 1 IFNs have been described to be beneficial for the control of tumor growth, including in breast cancer, either by exerting direct antitumor effects 67 , or by improving the response to therapies, such as chemotherapy and checkpoint inhibition 68,69 .Yet, opposite roles in inducing treatment resistance have been described as well 42,70,71 .By using defined BC immune-organoid co-cultures, we have shown that an IFN-I signature intrinsic to tumor cells associates with TEG sensitivity, and that IFN-β primes tumor cells for more efficient targeting, rather than directly affecting TEG killing behavior.Thus, our data support the clinical use of IFN-I in combination with TEGs and possibly other cellular immunotherapies.
Adding to patient-specific drug responses observed in PDOs biobanks [14][15][16][17][18]72 , we show that not only killing efficacy, but also the underlying behavioral and molecular mechanisms of cellular immunotherapy differ between different PDO cultures. Weeven detected differences in killing dynamics between individual organoids belonging to the same PDO culture that appeared to arise from intrinsic biological differences between individual organoids.This demonstrates that our platform captures inter-and intra-patient heterogeneity, a major obstacle for treating solid tumors 73 .It is intriguing that gene signatures induced in TEGs upon organoid engagement were partly dictated by the type of PDO.In addition, the extent of IFN-β pretreatment outcome on tumor targeting differed between PDOs, with the highly resistant BC culture 100T remaining unresponsive, whereas 34T displayed the highest (4-fold) increase in targeting.Together, these findings warrant caution regarding generalizing the outcome of immuno-oncology studies that use a single tumor model, and further supports the value of human organoid technology for development of personalized therapies.
Altogether, BEHAV3D combines organoid, imaging and sequencing technologies to offer a comprehensive platform that integrates multiple single-cell readouts, including tumor death dynamics, single-cell behavior and underlying transcriptomic profiling (Supplementary Video 1).BEHAV3D may thus contribute to efforts aimed at enhancing the efficacy of solid tumor-targeting by cellular therapies.

Cell lines
Daudi (CCL-213) 6 , HL60 (CCL-240) 8 and Phoenix-Ampho (CRL-3213) cell lines were obtained from ATCC.Phil Greenberg (Fred Hutchinson Cancer Research Center, Seattle, USA) kindly provided LCL-TM.Daudi, HL60 and LCL-TM cells were cultured in RPMI media supplemented with 10% fetal calf serum (FCS) and 1% pen/strep (all from Thermo Fisher).Phoenix-Ampho cells were cultured in DMEM medium (Thermo Fisher) supplemented with 10% FCS and 1% pen/strep.All cells were passaged for a maximum of 2 months, after which new seed stocks were thawed for experimental use.Stocks were reauthenticated by short tandem repeat profiling/karyotyping analysis provided by Eurofins Genomics in 2019 (Daudi), 2021 (LCL-TM and Phoenix-Ampho) and 2022 (HL-60), respectively.Furthermore, all cell lines were routinely verified by growth rate, morphology, and/or flow cytometry and tested negative for mycoplasma using MycoAlert Mycoplasma Kit.
Peripheral blood mononuclear cells (PBMCs) were obtained from Sanquin Blood bank (Amsterdam, The Netherlands) and isolated using Ficoll gradient centrifugation methods from buffy coats.
Transduced T cells were then stimulated using a REP 8 as described above for TEGs.Transgenic TCR expression and purity of CD8 + populations was routinely assessed by flow cytometry.

ROR1 CAR T cells
The ROR1-specific chimeric antigen receptor (CAR) T cell was previously described 80 .The CAR sequence was derived from patent WO2016187216 and cloned into the pCCL-cPPT-hPGK-GFP-bPRE4-SIN lentiviral transfer vector in place of the GFP sequence, which was derived from the pRRL-cPPT-hPGK-GFP-bPRE4-SIN plasmid 81 .Lentiviral particles were produced using standard calcium phosphate transfection (Clontech) of HEK-293T cells (ATCC CRL-3216) as described previously 82 .CD8 + T cells were separated from cord blood by ficoll isopaque density gradient centrifugation (GE Healthcare Bio-Sciences AB) and magnetic bead separation (Miltenyi Biotech).Cells were expanded and transduced according to a previously established protocol 82 .Transduced cells were subsequently FACS sorted for CAR expression (anti-IgG1-FITC) on a BD FACSAria II, after which the pure CAR expressing T cell population was stimulated using a REP 8 as described above for TEGs.

T cell serial killing capacity analysis
For accurate long-term (up to 24 h) T cell tracking, TEGs were plated at an E:T ratio of 1:25.
Tracks were manually corrected where required.Tracks were divided into shorter subtracks of 160 minutes.Using the random forest classifier described above (Extended Data Fig. 3d,e), each subtrack was assigned to a behavioral signature.The following statistics were calculated for each type of behavioral signature (9 clusters): for continuous variables (square displacement; speed, T cell death) the mean, median and standard deviation of the upper quantile were calculated, and for discrete variables (organoid contact and interaction with T cells) the mean, cumulative mean, maximum and cumulative maximum were calculated.
Principal component (PC) analysis was used to reduce the dimensionality.The top 5 PCs were used to classify the change in behavioral signature over time (Extended Data Fig. 4b,c).
Equivalent to the approach that was used for the full tracks in Fig. 2b, we computed a crossdistance matrix based on the multivariate time-series data using the dynamic time warping algorithm and performed k-means clustering in UMAP space.The change in behavioral signatures was represented in a time-series color plot where each row represents one cell track and the color codes for behavioral signature (Extended Data Fig. 4c).The relative proportion of CD4 + and CD8 + TEGs in each cluster was calculated and plotted next to each long-term classification (Extended Data Fig. 4c).
To quantify the number of cells killed by a TEG we divided the killed volume by the average volume of a single 13T cell (2182 μm 3 ).

PDO bulk RNA sequencing
For bulk RNA sequencing characterization, RNA of PDOs grown in 'Type 1' culture medium was isolated according to the manufacturer's protocol using the RNeasy Mini Kit (QIAGEN).
Quality and quantity of the RNA samples and the libraries were measured with Agilent's Bioanalyzer2100 and Invitrogen™ Qubit™ 3.0 Fluorometer.Quality control was done using FastQC, alignment has been done using STAR (https://github.com/alexdobin/STAR/releases/tag/STAR_2.4.2a) and reads have been mapped to the GRCh37 version of the human reference genome.Quality control on the bams was done using Picard.Read counts were generated with Htseq-count after which normalization was done using DESeq.RPKMs have been calculated with edgeR.For the library preparation the TruSeq Stranded mRNA Library Prep kit from Illumina was used.Sequencing was performed on the nextseq500 sequencer (also Illumina) with single-end 75bp reads.PDO cultures were ranked by responsiveness to TEGs (Fig. 1d) and differentially expressed genes between the 6 most TEG-sensitive and 6 least TEG-sensitive cultures were analyzed.Genes exhibiting a more than 4-fold expression change with an adjusted p-value <0.05 after multiple hypothesis testing correction were used as input gene set enrichment analysis.

SORTseq library preparation and sequencing
All sorted plates were processed according to the CEL-Seq2 protocol with the total transcriptome amplification via poly-A RNA-capture, library preparation, and sequencing into Illumina sequencing libraries as previously described 83 .Paired-end sequencing (read1: 30 bp; read2: 120 bp) was used to sequence the prepared libraries using an Illumina NextSeq sequencer.
integration 85 .Briefly, all three datasets were normalized using SCTtransform 86 followed by selection of 5000 features for downstream integration.Shared nearest neighbor graph-based clustering was done using Seurat package's FindNeighbors and FindClusters functions with a resolution of 0.8.For cell type identification marker genes for each cluster were calculated using the FindAllMarkers function and examined to profile marker genes that correspond to known cell types.Additional support for identifying cell subpopulations similitudes was achieved by analyzing the differentially expressed genes with a cell-type annotation tool 87 .

Differential gene expression analysis of TEGs co-cultured with distinct PDO cultures
For comparison of TEGs targeting 10T or 13T PDOs (Fig. 6a-c

Gene set enrichment analysis
The functional enrichment analysis in this study for pathway and biological processes annotations for gene sets of interest was conducted using ToppFun on the ToppGene Suite 88 (Fig. 1h, Extended Data Fig. 6c, Extended Data Fig. 8b).An enrichment score was assigned based on gene enrichment ratio and log p value.For redundant annotations, the annotation with the highest gene enrichment ratio was selected.
), SORTseq dataset was used including TEGs from distinct Experimental engagement states: Non-engaged Enriched and super engager.No target control TEGs were used as a control group.SORTseq data were mapped and quantified and visualized with UMAP as described above.Differential gene expression analysis was performed with the FindMarkers function from Seurat v3.Common and specific gene sets were filtered and visualized by Venn diagram with the VennDiagram package.