Coordinated immune networks in leukemia bone marrow microenvironments distinguish response to cellular therapy

To dissect immune cell states in the BM microenvironment of patients with post-HSCT relapsed AML treated with DLI, we profiled the transcriptomes of individual marrow cells from 30 biopsy samples [five responders (Rs) and four nonresponders (NRs); Fig. 1A and data file S1]. Patients in both groups had similar baseline, transplant, and DLI characteristics; however, there were more female recipients in the R group (fig. S1, A to C, and data file S2). A total of 47,372 viable bone marrow mononuclear cells (BMMCs) were profiled by scRNA-seq, single-cell T cell receptor sequencing (scTCR-seq), and cellular indexing of transcriptomes and epitopes by sequencing (scCITE-seq) for proteomic characterization. We further included seven DLI product samples [43,919 cells from four AML responders (AML-Rs) and three AML nonresponders (AML-NRs)]; two post-HSCT samples from nonrelapsed patients with AML (7617 cells); and, as further controls, six samples (27,131 cells) from three patients with post-HSCT relapsed AML who entered remission after chemotherapy without DLI.

To expand our BMMC compendium for myeloid malignancies, we merged the AML data with our prior dataset of 42 marrow specimens from 16 patients with CML who received DLI and 2 nonrelapse controls (11) and extended the analysis from T cells to the entire BM microenvironment. After batch correction and clustering (25), we constructed an integrated atlas of 451,453 mononuclear cells (Fig. 1B). Fifty-three distinct cell clusters were annotated using both reference-based and manual annotation (Fig. 1B; fig. S1, D and E; and data files S3 and S4). Expression of lineage-defining genes confirmed major cell type annotations (fig. S1F).

Distinguishing leukemic cells transcriptionally from healthy hematopoietic stem cells (HSCs) and myeloid progenitor cells is challenging because of partial overlap in differentiation trajectories and interpatient heterogeneity (23, 26, 27). We annotated leukemia cells on the basis of a combination of gene expression (fig. S1G) (28), chromosome Y genes for sex-mismatched donor/recipient pairs (fig. S1H), and copy number variants (CNVs) (fig. S1, I and J) (29) to create a composite donor/recipient map of cells (fig. S1K). To enable analysis of shared phenotypes in leukemic clusters specific to patients with AML, we combined myeloid clusters into seven metaclusters (MCs) and then subclustered on the basis of donor (MC1-7) or recipient origin (MC1-7-leukemia) (fig. S1L).

We next asked whether CML and AML had divergent evolutionary trajectories from each other and from normal myeloid development. We applied Decipher (27), a deep generative model for reconstructing and aligning cell trajectories. Decipher revealed three continuous trajectories corresponding to shared normal myeloid differentiation in CML and AML specimens expressing VCAN, FCN1, CD14, and S100A8 (30) and divergent CML and AML leukemic evolutionary paths, with higher expression of CXCL8, RACK1, B2M, TRIM8, and IL18 for AML (31–34) and KCNH2, HIST1H4C, and HSPA8 (35) for CML (Fig. 1C and data file S4). Nonleukemia clusters were overall evenly distributed by response (P = 1.0, t test), dataset [P = 1.0, analysis of variance (ANOVA)], and donor/recipient (P = 1.0, t test; Fig. 1, D and E, and figs. S1K and S2, A to C), indicating robust correction for batch and library size. After normalizing for total cell number for AML and CML, few clusters were dominated by cells from one disease [n = 7 (13%) for AML, n = 10 (19%) for CML] (Fig. 1F). Together, construction of this cell state atlas revealed distinct disease phenotypic trajectories and identified both shared and disease-specific immune cell populations.

To identify cell subsets associated with response to DLI, we evaluated individual cluster compositions of marrow cells by disease (AML versus CML), time point (before DLI versus after DLI) and response (Rs versus NRs) (data file S4). We hypothesized that populations expanding after DLI may represent cellular mediators of a GvL response. We also assessed clusters enriched in the BM microenvironment before DLI in Rs or NRs as candidate predictors of treatment response. One population trended toward expansion in Rs but not NRs with CML: T cell C2 (P = 0.06 and P = 0.70, respectively; Fig. 2, A and B). This cluster expressed CD4, SELL, CCR7, TCF7, and IL7R (fig. S2, D to H), aligning with our previously identified T_PE cell population (11). This cluster was not expanded in AML-Rs, suggesting involvement of alternate cell types in the GvL response in AML. A tumor-infiltrating natural killer (NK) population (cluster C36) was found to be enriched in CML Rs before DLI (figs. S2, G and H, and S3A).

AML-Rs were characterized by expansion of T cell clusters C0 (P = 0.003) and C40 (P = 0.04) with a trend toward expansion in T cell C5 (P = 0.07), whereas AML-NRs were characterized by expansion of T cell clusters C21 (P = 0.01) and C26 (P = 0.001) (Fig. 2C, fig. S3B, and data file S5). None of these clusters were expanded in patients with CML. NK cluster C1 was expanded in AML-Rs (P = 0.002) but not in CML-Rs (P = 0.94) or AML-NRs (P = 0.49); however, AML-NRs had a wide range of C1 NK cell proportions across samples. T cell clusters C0 and C40 displayed higher CD8A and effector/cytotoxicity gene expression (e.g., ZNF683/Hobit, GZMB, PRF1, and B3GAT1/CD57; fig. S2, D and E) than other CD8⁺ T cell clusters. Compared with C40, the C0 cluster had even higher expression of ZNF683 and CD45RA by scCITE-seq (P = 0.0004, Kolmogorov-Smirnov test; fig. S2F), consistent with a CTL or T effector memory reexpressing CD45RA (T_EMRA) cell profile (36, 37), and was thus designated as CD8⁺ZNF683^Hi CTLs. The profile of T cell cluster C40 was most consistent with CD8⁺ T effector memory (T_EM) cells.

NK cluster C1 displayed high expression of B3GAT1, FCGR3A, and GZMB and higher scCITE-seq expression of CD57 compared with other NK clusters, suggesting a cytolytic phenotype (fig. S2, G and H) (38). The C0: CD8⁺ZNF683^Hi CTL, C1: cytolytic NK cell, and C40: CD8 T_EM cell clusters did not expand in chemotherapy-only (no DLI) controls (fig. S3C), and their median proportions were similar to that of nonrelapse controls. To account for possible confounding variables that could overrepresent expansion of C0: CD8⁺ZNF683^Hi CTLs, we used several additional analyses. Expansion of C0: CD8⁺ZNF683^Hi CTLs was observed even after (i) calculating their proportion from only lymphocytes to account for possible higher leukemia burden in NRs (P = 0.007; fig. S3D), (ii) limiting samples to the first 300 days after DLI to account for longer follow-ups in Rs (P = 0.003; fig. S3E); and (iii) down-sampling cells to account for a greater number of cells in R samples than NR samples (P = 0.002; fig. S3F). Given that our prior CML cohort (11) was treated with CD8-depleted DLI, we additionally evaluated pre- and post-DLI samples from four patients with CML treated with conventional (i.e., CD8-replete) DLI (Fig. 1A). Similar to patients with CML treated with CD8-depleted DLI (39), we did not detect post-DLI expansion of C0: CD8⁺ZNF683^Hi CTLs in Rs (fig. S3, G and H).

We asked whether any cell types were enriched in pre-DLI samples. Erythroid [C6: CML-enriched (P = 0.01) and C7: CML-enriched (P = 0.03)] and immature populations [C16: myeloerythroid progenitors (MEPs) (P = 0.03) and C37: MEPs (P = 0.02)] were elevated in AML-NRs compared with AML-Rs before DLI. No cell types were enriched in AML-Rs before DLI. Together, these analyses support the notion that multiple immune cell subtypes in the AML BM microenvironment undergo dynamic changes after DLI.

We reasoned that the dynamic patterns of response might arise from interactions and cellular cross-talk within the leukemic BM microenvironment over time. However, current computational methods for inferring cell-cell interactions from scRNA-seq data generate a static model of predicted interactions. We therefore applied DIISCO (dynamic intercellular interactions in single-cell transciptomics) (40), our newly devised Bayesian model, which uses a Gaussian process regression network (41, 42), to identify dynamic networks of cell-cell interactions within the BM microenvironment as “interactomes,” incorporating time relative to treatment as a variable (fig. S4A). In CD8-depleted CML DLI, DIISCO identified a central role for expanding C2: CD4⁺ T_PE cells in Rs after DLI, predicted to inhibit CML1 cells and HSCs (C15), through coordination with C13: mature NKs, C0: CD8⁺ZNF683^Hi CTLs, and C20: plasma B cells (Fig. 2B and fig. S4, B to D). This C2-centralized inhibitory role was not observed in CML nonresponders (CML-NRs) (fig. S4D). Predicted ligand:receptor complexes correlated with inferred interactions between C2: CD4⁺ T_PE cells and CML1 cells included the protein products encoded by TNF:LTBR and CD226:NECTIN2, which have been associated with effective T cell responses (fig. S4E) (43, 44). These predicted networks support the central role of C2: CD4⁺ T_PE cells in mediating GvL in CML.

For AML-Rs, the DIISCO analysis, in contrast, revealed a cascading response with multiple cell types expanding and interacting after DLI at different time points centering around C0: CD8⁺ZNF683^Hi CTLs, C40: CD8 T_EMs, and C5: CD8⁺ exhausted T (T_EX) cells and including C1: cytolytic NK cells (Fig. 2, C to E, and fig. S5A). A strong inhibitory interaction between C0: CD8⁺ZNF683^Hi CTLs and AML (MC1-leukemia) was observed in Rs (Fig. 2E), consistent with a productive GvL-driven contraction of MC1-leukemia cells. DIISCO identified positive interactions among C0: CD8⁺ZNF683^Hi CTLs and C40: CD8 T_EM cells along with C1: cytolytic NK cells, supporting the notion that these cell subsets formed a coordinated immune network supporting enhanced GvL activity (movie S1). To confirm that these predictions were not primarily attributable to increased immune infiltration, we retrained DIISCO on only immune clusters, resulting in recapitulation of this coordinated immune network (fig. S5B). These analyses reinforced the role of C0: CD8⁺ZNF683^Hi CTLs, C40: CD8 T_EM cells, and C1: cytolytic NK cells in the DLI response in AML.

Given that the proportion of C0: CD8⁺ZNF683^Hi CTLs varied inversely with the proportion of leukemia cells in AML-Rs over time (Fig. 2D and data file S2), we postulated that this cluster might interact with leukemia, leading to T cell activation in AML-Rs but T cell inhibition in AML-NRs. We observed inhibitory signals from C0: CD8⁺ZNF683^Hi CTLs toward AML-MC1 (e.g., GZMA:F2R, CCL4:SLC7A1) associated with CD8 T cell cytotoxicity and tumor suppression (Fig. 2F) (45–47). Intramarrow interactions in AML-Rs after DLI were driven predominantly by C0: CD8⁺ZNF683^Hi CTLs (P < 0.0001, Mann-Whitney U) (fig. S5B). Beyond lymphoid cells, DIISCO revealed inhibitory links in Rs between C0: CD8⁺ZNF683^Hi CTLs and myeloid MC2 (dendritic cells and classical monocytes; Fig. 2E, fig. S5B, and data file S4). Conversely, in AML-NRs, DIISCO revealed the absence of notable activating interactions between T cells and other immune cell types after DLI.

We further asked whether CD8⁺ZNF683^Hi CTLs could also be identified in peripheral blood. Analysis of patient R3501, whose scRNA-seq and scTCR-seq data from date-matched peripheral blood mononuclear cells (PBMCs) and BMMCs revealed T cell clonotypes in common, suggested that at least some T cell clonotypes were detectable in both compartments (fig. S5C). Moreover, the phenotype of the circulating post-DLI PBMC T cells was consistent with BMMC C0: CD8⁺ZNF683^Hi CTLs (high expression of ZNF683 and GZMH and scCITE-seq expression of CD45RA). We therefore evaluated the proportion of circulating CD8⁺ T_EMRA cells in PBMCs from 53 independent patients similarly treated with DLI for post-HSCT relapsed AML (29 Rs and 24 NRs). Multiparameter flow cytometry of PBMCs from this larger cohort had been performed in real time before and after DLI at our institution (data file S6). We detected marked expansion of the CD8⁺ CTL (T_EMRA) population after DLI in Rs (fold change: 2.56 after DLI versus before DLI, P < 0.001, Wilcoxon rank sum) but not NRs (fold change: 1.3, P = 0.47, Wilcoxon rank sum) (Fig. 2G and fig. S5D). Excluding patients on active therapy (i.e., receiving systemic antileukemia therapy within 30 days of the analyzed sample) confirmed this finding (fig. S5E). A similar increase in AML-R after DLI was also seen in a naive B cell population (fig. S4, E and F). In our scRNA-seq dataset, C3: naive B cells skewed toward expansion in four of five AML-Rs (figs. S3D and S5, F and G). This independent PBMC confirmation of our BM scRNA-seq findings in a larger independent validation cohort supported a key role of CD8⁺ CTLs in effective GvL in AML. Together, application of DIISCO confirmed key roles of distinct T cell clusters in CML and AML and postulated a model of coordinated immune interactions driving effective GvL in AML centered around C0: CD8⁺ZNF683^Hi CTLs.

Organized networks of immune cells within solid TMEs have been linked to improved outcomes after immunotherapy (5–7, 48–50). We hypothesized that the BM microenvironment of Rs might be spatially organized into lymphoid networks involving C0: CD8⁺ZNF683^Hi CTLs, potentially promoting GvL. To confirm our DIISCO predictions that CD8⁺ CTLs function as a central hub for response to DLI in AML, we used codetection by indexing (CODEX) (51) to spatially assess the organization of immune cells at single-cell resolution. This was applied to six pre-DLI (four R and two NR) and six post-DLI (three R and three NR) BM core biopsy specimens obtained concomitantly with aspirates analyzed by scRNA-seq (data file S1). Upon quantifying the proportion of T_EMRA cells (coexpressing CD3, CD57, and granzyme B) in five R and two NR samples containing >50 segmented cells, we confirmed T_EMRA cell expansion relative to all cells in Rs (Fig. 3A).

To systematically evaluate the immune networks defining response and resistance, we identified shared spatial neighborhood patterns across samples. Cells were aggregated across all R samples (four pre-DLI and three post-DLI) and separately across all NR (two pre-DLI and three post-DLI) samples and then clustered by composition of major cell types in their neighborhoods. Given that our marker panel (data file S7) was designed to resolve T cell subsets and other major immune cell types rather than leukemia cells, we focused on defining niches according to immune cell colocalization. After filtering, 98,130 cells (75,350 for Rs and 22,780 for NRs) were assigned to one of four major cell types (i.e., T and NK cells, B cells, HSC or leukemia cells, or myeloid cells), and the number of cells assigned to each major cell type within its neighborhood was quantified. Recurring colocalization patterns were designated as “niche types” using Phenograph (Fig. 3B). Together, we identified 37 distinct niche types in Rs (excluding two sample-specific niches, r0 and r1) and 19 distinct niche types in NRs (data file S8).

The R niche types (i.e., r2 to r37) typically comprised diverse immune cell populations over time (Fig. 3C), whereas the NR niche types (nr1 to nr19) lost diversity after DLI (quantified by Shannon entropy, P = 0.01, Mann-Whitney U; Fig. 3, C and D). We observed a trend toward higher cellular density (P = 0.06, Mann-Whitney U) in pre-DLI–enriched niches in Rs compared with NRs, suggesting preexisting diverse immunologic neighborhoods as a distinguishing feature of response. In addition, we observed higher abundance of T and B cells dominating pre-DLI R niche types versus NR niche types (P = 0.004 and P = 0.04, respectively; Mann-Whitney U), whereas myeloid cells exhibited higher proportions in NRs versus Rs before DLI (P = 0.002; Fig. 3, C and D). T cell abundance was further increased in Rs after DLI versus before DLI in parallel with decreased niche diversity (P = 0.0004; Fig. 3D).

Deep annotation of T cell phenotypic states across time identified that expanded T cells in Rs shifted phenotype from inhibitory [expressing T cell immunoglobulin and mucin-domain containing-3 (TIM3) and cytotoxic T lymphocyte associated protein 4 (CTLA4)] to effector [expressing lymphocyte-activation gene 3 (LAG3), granzyme B, CD57, and OX40] after DLI (Fig. 3E). Post-DLI niches in Rs had higher T cell expression of CD57 and OX40 (Mann-Whitney U test, P < 1 × 10⁻¹⁵), whereas NRs displayed enrichment of niches expressing the inhibitory marker T cell immunoreceptor with immunoglobulin and immunoreceptor tyrosine-based inhibitory motif (ITIM) domains (TIGIT) (Fig. 3F). Although CD57 is sometimes considered a marker of proliferative senescence, numerous studies have demonstrated effector function of CD8⁺CD57⁺ T cells (52–55). Down-sampling analysis to account for differences in cell numbers between Rs and NRs confirmed these marker expression differences (fig. S6, A and B). This transition to effector and cytotoxic T cell states in the BM microenvironment of Rs complemented the expansion of C0: CD8⁺ZNF683^Hi CTLs and C40: CD8 T_EM cells found in our temporal scRNA-seq analysis with DIISCO. Spatial mapping of pre-DLI R cells further underscored the colocalization of T_EX cells with other immune cells, including B cells and myeloid and leukemia cells (Fig. 3G). Niche types dominating Rs after DLI included CD8⁺ CD57⁺ CTLs colocalized with other cytotoxic and T_EM cells and myeloid cells enriched in dendritic cell marker CD11c (Fig. 3, E and H, and fig. S6C). In contrast, niche types in NRs reflected low diversity and dominance of myeloid cell types (Fig. 3I).

Our CODEX analysis was limited to a protein marker panel designed to capture T cell subsets. To extend the phenotypic analysis of spatially resolved marrow-resident immune cell subpopulations, we performed targeted spatial transcriptomic profiling (153 transcripts; Fig. 4A; fig. S7, A and B; and data file S9) paired with eight protein surface markers (56) on a subset of formalin-fixed paraffin-embedded BM samples from our CODEX discovery cohort and additional samples from our flow cytometry validation cohort (data file S6). Cells were clustered and manually annotated using differentially expressed gene (DEG) analysis (Fig. 4B and fig. S7C). T and NK cells were subclustered and finely annotated to identify CD8⁺ T_EMRA cells (Fig. 4C and fig. S7D). In this manner, we confirmed that R samples after DLI globally had higher cellular diversity than NR samples after DLI (Shannon diversity index, P = 0.04, t test; Fig. 4D). In addition, we observed a higher diversity in the local neighborhood of CD8⁺ T_EMRA cells and other immune cell types in Rs after DLI than in NRs (P < 0.05; Fig. 4E and fig. S7E), corroborating the paradigm of coordinated immune networks driving GvL. The higher diversity in CD8⁺ T_EMRA neighborhoods in Rs was independent of the choice of local neighborhood radius (fig. S7F).

To further examine local immune composition, 16 distinct cellular niche types were identified across both R and NR samples using an unsupervised neighborhood analysis approach similar to that used for our CODEX data (Fig. 4F). R samples had greater niche type diversity than NR samples after DLI (Shannon diversity, P = 0.003; Fig. 4G). Niches 2, 3, and 13 contained higher fractions of CD8⁺ T_EMRA cells and were enriched in Rs versus NRs (Fig. 4, H and I). Conversely, niches 6 and 10 were more enriched in NRs after DLI than Rs and exhibited a higher fraction of AML and myeloid cells (Fig. 4I). Niche 15 contained the highest fraction of CD8⁺ T_EMRA cells and showed an increase in fraction in Rs after DLI, but statistical significance could not be established because of the sample size (P = 0.32; Fig. 4J). The spatial analyses of samples from both the discovery and validation cohorts thus provided orthogonal confirmation of post-DLI expansion of CD8⁺ T_EMRA cells and enhanced cellular diversity of immune niches in the AML-R BM microenvironment.

Our multimodal analyses demonstrated that expansion of effector T cells is a central feature in the AML response to DLI. To better characterize the T cell phenotypes in AML-Rs and AML-NRs associated with response and resistance, we once again used Decipher (27) to map trajectories spanning all T cell clusters. We identified one key latent factor describing a trajectory toward terminally differentiated states (fig. S8A) and a second discriminating between CD4 and CD8 subsets (Fig. 5A). When projecting cells on these latent factors, we observed a trajectory recapitulating T cell differentiation on the basis of gene expression profiling within our CD8⁺ T cell clusters. T cells from AML-Rs skewed toward effector states (enriched in GZMH, CCL5, and NKG7; Fig. 5B and fig. S8, B and C), whereas those of AML-NRs skewed toward terminal differentiation with higher expression of coinhibitory molecules (TIGIT, KLRG1, and TCF7) (distribution shift in terminal differentiation trajectory, P = 2.2 × 10⁻¹¹, Kolmogorov-Smirnov test) (Fig. 5, A and B, and fig. S8D).

Comparison of our T cell clusters with those in previous TME studies confirmed the CD8⁺CTL and T_EMRA characteristics of C0: CD8⁺ZNF683^Hi CTLs (Fig. 5C). Gene expression signatures of our six CD8⁺ clusters compared with T cell MCs from a large scRNA-seq pan-cancer atlas of tumor-infiltrating CD8 T cells curated from 316 patients across 21 cancer types (57) showed that C0: CD8⁺ZNF683^Hi CTLs most closely aligned with the CD8 MC signature c06.Temra.CX3CR1 (Fig. 5C, left). A key difference between this MC and C0: CD8⁺ZNF683^Hi CTLs, however, was the high expression of ZNF683 in the latter (57). We also evaluated profiles of marrow-derived T and NK cells from healthy volunteers and from patients with AML treated with immune checkpoint blockade (ICB) (58). C0: CD8⁺ZNF683^Hi CTLs most closely aligned with CD8⁺ CTL populations from the marrow of patients with AML but not healthy marrow T cells (Fig. 5C, middle). Last, we compared our T cell populations with a recent study of Richter syndrome, which identified a population of ZNF683^Hi CD8⁺ T cells associated with response to PD-1 ICB (37). C0: CD8⁺ZNF683^Hi CTLs were most transcriptionally similar to the ZNF683-intermediate (ZNF683^Int) CD8⁺ T cell population, followed by the ZNF683^Hi CD8⁺ subset (Fig. 5C, right).

Examination of our C0: CD8⁺ZNF683^Hi CTLs demonstrated expression of some genes more commonly associated with NK cells (e.g., NKG7, B3GAT1, and FCGR3A), raising the possibility that this cluster had NK-like properties. A recent report identified an NK-like population of CD8⁺ T cells, expanded by cytomegalovirus (CMV) exposure and defined by low BCL11B expression, with effector function against leukemia cell lines (59). To determine whether our C0: CD8⁺ZNF683^Hi CTLs were compatible with this previously described NK-like population, we compared our NK and T cell clusters with that of the NK-like population but did not find concordance between expression profiles (fig. S8E), confirming C0: CD8⁺ZNF683^Hi CTLs as distinct from the previously described CMV-reactive NK-like CD8 population.

We asked whether C0: CD8⁺ZNF683^Hi CTLs were enriched in the post-HSCT setting rather than other AML treatment settings and therefore reanalyzed two scRNA-seq datasets generated from marrow and blood. We reclustered marrow-derived T cells collected from patients with post-HSCT–relapsed or HSCT-naive–relapsed myelodysplastic syndrome (MDS) or AML after treatment with the DNA methyltransferase inhibitor decitabine and ipilimumab and identified four clusters with ZNF683^Hi expression (cluster IDs 2, 3, 4, and 10) similar to C0: CD8⁺ZNF683^Hi CTLs (fig. S8, F to H) (60).These cells were more prevalent in post-HSCT patients than in non-HSCT patients (P = 0.01, Mann-Whitney U), supporting the notion that they arise in the setting of immune reconstitution after HSCT (Fig. 5, D and E). We then reanalyzed the BM T cells from the aforementioned ICB-treated AML cohort (58) and identified a subset of ZNF683^Hi CD8⁺ T cells, most consistent with C0: CD8⁺ZNF683^Hi CTLs (fig. S8, I and J). We detected a small increase in this population in Rs after ICB (P = 0.27) but not NRs (P = 0.85, paired t test; fig. S8K).

Given that C0 was present in both AML-Rs and AML-NRs, we hypothesized that distinct properties of the cells might explain their association with response. DEG analysis among C0: CD8⁺ZNF683^Hi CTLs identified ZNF683 as highly differentially expressed in AML-Rs versus AML-NRs (P < 0.0001; Fig. 6A and data file S10). AML-R C0: CD8⁺ZNF683^Hi CTLs also displayed higher expression of activation genes (NKG7, STAT1, JUNB, and IFNG). In contrast, AML-NR cells highly expressed KLRB1 (encoding CD161), an inhibitor of tumor-infiltrating lymphocyte (TIL) cytotoxicity (61), and multiple metallothionein pathway genes (MT1X, MT2A, MT1F, and MT1E) associated with TIL dysfunction (Fig. 6A) (62).

To better define potential pathways of interest in AML-R versus AML-NR C0: CD8⁺ZNF683^Hi CTLs and leukemia cells, we applied CellPhoneDB (63) and CellChat (64) to identify possible receptor-ligand (R-L) pairs between these clusters (data file S11). AML-NRs demonstrated inferred interactions between the coinhibitory marker TIGIT on CD8⁺ CTLs and NECTIN2 on AML leukemia cells (Fig. 6B and data file S11) (65–67). The TIGIT axis is a well-described immune checkpoint pathway for T and NK cell activation (with CD96 or CD226) and T cell inhibition [with Nectin-2 or poliovirus receptor (PVR)] (44, 68). Accordingly, potential NECTIN2 or PVR interactions with CD226 but not TIGIT were inferred by CellPhoneDB in AML-Rs, corroborating a more activated profile in Rs (Fig. 6B). Thus, C0: CD8⁺ZNF683^Hi CTLs from AML-Rs are predicted to be more highly activated and cytotoxic after DLI compared with cells from AML-NRs that are more inhibited and dysfunctional, consistent with their potential roles in GvL response and resistance, respectively.

To test whether CD8⁺ CTLs had functional antileukemic activity, we devised an in vitro coculture assay to evaluate leukemia-specific T cell activation of CD8⁺ CTLs after encounter with leukemia cells from the same patient (fig. S9A). CD8⁺ T_EMRA CTLs from six AML-Rs were more activated and cytotoxic, expressing higher CD137 (P = 0.024, t test, corrected with Benjamini-Hochberg) and interferon-γ (IFN-γ) (P = 0.003, t test) after coculture with primary patient-matched leukemia cells than CD8⁺ T_EMRA cells from three AML-NRs (Fig. 6, C to F). CD8⁺ CTLs cultured in the absence of leukemia cells (Fig. 6, C to F) or in the presence of patient-matched (recipient) BM fibroblasts (fig. S9, B and C) did not show elevated expression of either CD137 or IFN-γ, suggesting leukemia-specific activation in Rs but not NRs. This activity was not observed in CD8⁺ T cell subsets aside from T_EMRA cells (fig. S9, D and E), indicating leukemia-specific activation of CD8⁺ CTLs, consistent with a potential functional role of these cells in AML DLI Rs.

To determine whether the activated and cytotoxic-appearing in vitro population was consistent with the C0: CD8⁺ZNF683^Hi CTLs implicated in GvL by scRNA-seq, we performed the same coculture experiment on samples from three additional Rs. Cells were then sorted by flow cytometry into CD8⁺CD137⁺ and CD8⁺CD137⁻ populations and analyzed by scRNA-seq and scTCR-seq (fig. S9, F and G). The CD8⁺CD137⁺ population was enriched for a CD8⁺ T_EMRA/CTL gene expression signature compared with the CD8⁺CD137⁻ population (fig. S9, H and I), demonstrating similarity to C0: CD8⁺ZNF683^Hi CTLs (Fig. 6G and fig. S9J). Together, these characterizations provided independent support for a likely role of ZNF683-expressing CD8⁺ CTLs in mediating antileukemia responses.

Because expansion of the C0: CD8⁺ZNF683^Hi CTL population is central to the antileukemia response, we asked whether the increased proportion of C0: CD8⁺ZNF683^Hi CTLs in AML-Rs after DLI was due to a higher proportion in the DLI product. Analysis of cell type proportions within the DLI products of four AML-Rs and three AML-NRs revealed no differences in C0: CD8⁺ZNF683^Hi CTL percentages (P = 0.88) in other expanding BM clusters (Fig. 6H and fig. S10A) or in other major cell types (fig. S10B). Alternatively, CD8⁺ZNF683^Hi CTLs may have undergone greater clonal expansion after infusion because of increased support from a more immunologically diverse BM microenvironment. In support of this hypothesis, scTCR-seq analysis of T cell clonotypes within BM and DLI products of AML-Rs demonstrated marked clonal expansion primarily in C0: CD8⁺ZNF683^Hi CTLs (Fig. 7, A to C). Most (63%) clonally expanded cells derived from the post-DLI BM, with 78% of the post-BM expanded clones found in C0: CD8⁺ZNF683^Hi CTLs. Expanded clones expressed more ZNF683, whereas unexpanded clones expressed more TIGIT, corroborating our R versus NR DEGs (Fig. 7D). The C28: CD8 T_Unconv cluster was noted to be highly clonal. Examination of this cluster revealed that 91.7% of the cluster (609 of 664 cells) comprised mucosal-associated invariant T (MAIT) cells, explaining their clonal/invariant nature. To determine whether clonal expansion correlated with T cell activation or inhibition, we analyzed the fold change of genes expressed in C0: CD8⁺ZNF683^Hi CTLs versus C5: CD8 T_EX cells and their relationship with clonal expansion (Fig. 7E). Activation genes (e.g., ZNF683, GZMB, PRF1, and CX3CR1) were more highly expressed in C0: CD8⁺ZNF683^Hi CTLs, associating with increased clonal expansion, whereas coinhibitory markers (e.g., TIGIT, LTB, and GZMK) were more highly expressed in C5: CD8 T_EX cells, anticorrelated to clonality. These data corroborate our in vitro findings that C0: CD8⁺ CTLs become activated and clonally expanded in response to AML cells.

We next examined whether the C0: CD8⁺ZNF683^Hi CTLs originated from the recipient’s BM or the DLI product on the basis of T cell receptor (TCR) clonotype tracking. Clonal expansion and activation of C0: CD8⁺ZNF683^Hi CTLs were primarily due to expansion of T cells from the DLI product, not the pre-DLI BM (Fig. 7F and fig. S10, C to E) (11, 17, 69). This finding was corroborated by down-sampling analysis to control for the variable number of T cells profiled (Fig. 7F, fig. S10F, and data file S12). Furthermore, higher TCR diversity was observed in the post-DLI BM of Rs compared with NRs (P = 0.01, t test on Gini coefficient; fig. S10G). TCR diversity within DLI products or the pre-DLI BM did not differ between AML-Rs and AML-NRs (P = 0.59 and P = 0.6, respectively, t test; fig. S10G), suggesting that the phenotypic state of T cells and their communication with other immune cells in the BM microenvironment were the primary driver of C0: CD8⁺ZNF683^Hi CTL activation and clonal expansion after DLI, rather than differences in the T cell clonality or composition of the DLI product. Thus, our findings support the notion that DLI infusion in the less exhausted, more immunologically diverse R BM microenvironment promotes clonal expansion of highly activated effector C0: CD8⁺ZNF683^Hi CTLs with GvL cytotoxic activity (Fig. 7G). Conversely, in NRs, infusion into the exhausted, less-diverse BM microenvironment (i.e., TIGIT-high) prevents expansion and cytotoxicity of C0: CD8⁺ZNF683^Hi CTLs, leading to TIL dysfunction and therapy resistance.

The primary aim of this study was to identify a cell population(s) associated with the response to DLI and effective GvL for relapsed myeloid malignancy after HSCT. We performed scRNA-seq, scTCR-seq, and scCITE-seq on cryopreserved BM biopsy specimens from nine patients with relapsed AML after HSCT who received DLI (five Rs and four NRs) and merged these data with previously generated data from similarly treated patients with relapsed CML after HSCT. Patients with AML had similar leukemia burden at relapse and at the time of DLI and were similar in other disease and treatment characteristics. We further confirmed our major findings with spatial profiling using CODEX and spatial transcriptomic profiling, as well as by flow cytometry in a larger independent cohort. All newly sequenced patient samples were obtained from patients receiving standard of care therapy with conventional CD8-replete DLI (i.e., not in a clinical trial). All samples analyzed were obtained from patients who had provided written informed consent for use of clinical data and samples on a biobanking protocol approved by the Institutional Review Board (protocol number 01206) at the Dana-Farber Cancer Institute/Harvard Cancer Center. These studies were conducted in accordance with the Declaration of Helsinki.

The patients with AML in this cohort were treated with DLI at the Dana-Farber Cancer Institute, Boston, between 2004 and 2021. Patients had similar baseline characteristics including HLA (human leukocyte antigen) match of donor, recipient and donor sex, time from HSCT to relapse, time from HSCT to DLI, and timing of pre-DLI samples (fig. S1, A and B, and data file S2). Rs had a longer duration of post-DLI marrow sampling because of longer survival after DLI. Patients in both groups had comparable disease burden at relapse, defined as percentage of BM blasts (fig. S1C). Three of five Rs and two of four NRs received cytoreducing therapy between relapse and DLI in an attempt to achieve remission before DLI (data file S2). Two Rs and one NR had minimal disease (i.e., <5% blasts but marrow dysplasia and/or progressive cytopenias and/or recurrence of mutations by next-generation sequencing) at relapse and did not receive therapy before DLI given the minimal disease burden. All patients except one NR had <5% blasts at the time of DLI. The exception was NR3512, whose disease at relapse was considered as smoldering, and therefore DLI was pursued without prior therapy despite 16% marrow blasts.

All patients received between one and three doses of DLI (median of 1 for Rs and 1.5 for NRs), with a median dose of 1 × 10⁷ CD3⁺ cells/kg for Rs versus 1.6 × 10⁷ CD3⁺ cells/kg for NRs (data file S2). Response was defined as durable morphologic and molecular remission for at least 1 year after DLI, whereas NRs lacked posttreatment reduction in disease burden. A median of three (range two to six) serial marrow biopsy specimens collected before and after DLI were evaluated per patient (Fig. 1A and data file S1). Response was defined as durable morphologic and molecular (i.e., cytogenetic and/or next-generation sequencing detection of disease-associated mutations) remission for at least 1 year after DLI, whereas nonresponse was defined as a lack of post-DLI BM blast reduction. Two patients who never relapsed were included as nonrelapse controls. Pre– and post–chemotherapy treatment samples from three patients who had post-HSCT relapsed AML and entered durable (>1 year) remission with chemotherapy but did not subsequently receive DLI were also included as chemotherapy-only (no DLI) controls. The presence of acute and chronic graft-versus-host disease (GVHD) was graded by standard consensus criteria (82); grades 0 and 1 acute GVHD were considered clinically equivalent (83). Characteristics of patients with CML were previously described (11).

BM aspirates and peripheral blood samples were collected and banked before and after DLI. BMMCs or PBMCs were isolated via Ficoll-Hypaque density gradient centrifugation, cryopreserved with 10% dimethyl sulfoxide, and stored in vapor-phase liquid nitrogen until the time of sample processing. For DLI products, 1 ml of the cryopreserved DLI apheresis product was obtained from quality control (QC) vials stored at the time of apheresis.

CML data were obtained from previously described sequencing runs, and all cells were analyzed for the current study, whereas we previously reported only the subset of data originating from T cells (11). For preparation of AML samples, cryopreserved primary BMMCs or PBMCs (for the matched sample from R3501) were thawed on the day of sequencing at 37°C and dispensed into a warmed solution of 10% fetal bovine serum (FBS) and 10% DNase I (StemCell Technologies, catalog no. 07900) in phosphate-buffered saline (PBS). The cell suspension was centrifuged at 200g for 10 min at room temperature. Viable cells were negatively selected using the MACS Dead Cell Removal Kit (Miltenyi Biotec, catalog no. 130-090-101). Collected live cells were resuspended in warmed cell staining buffer (BioLegend). Cells were stained with TotalSeq C hashtag and CITE-seq antibodies (data file S7) according to the manufacturer’s instructions. For batch AML1, after washing three times, cells from four samples were combined into one pool (data file S3, “Sequencing Run” column) with ~20,000 cells from each sample and diluted to a final concentration of 1000 cells/μl in 0.04% UltraPure bovine serum albumin (BSA) (Thermo Fisher Scientific). Seven total pools were processed in this way. Subsequent samples [i.e., batch AML2, DLI products, and control samples (data file S3)] were run individually. Cells were then taken immediately for scRNA-seq, scTCR-seq, and scCITE-seq. Approximately 50,000 cells from pooled samples or 10,000 cells from individual samples were loaded into one lane of a 10X Genomics Chromium instrument according to the manufacturer’s instructions. scRNA-seq libraries were processed using a Chromium Single Cell 5′ Library and gel bead v2 kit (10X Genomics). Coupled scTCR-seq libraries were obtained using a Chromium Single Cell V(D)J Enrichment Kit (human T cell) (10X Genomics), and coupled scCITE-seq libraries were obtained using a Chromium Single Cell Feature Barcode Kit. QC for amplified cDNA libraries and final sequencing libraries was performed using the Bioanalyzer High Sensitivity DNA Kit (Agilent). The scRNA-seq, scTCR-seq, and scCITE-seq libraries were normalized to a 4-nM concentration and pooled in a volume ratio of 4:1. The pooled libraries were sequenced on an Illumina NovaSeq S4 platform. The sequencing parameters were as follows: read 1 of 26 base pairs (bp), read 2 of 90 bp, index 1 of 10 bp, and index 2 of 10 bp. The sequencing data were demultiplexed and processed as described below.

We imaged seven R samples (four pre-DLI and three post-DLI) from four R and five NR samples (two pre-DLI and three post-DLI) from four NRs (data file S1). Four total samples (two from Rs and two from NRs) were imaged with a confocal microscope only, whereas all other samples were imaged with confocal and wide-field microscopy as well to better capture a low signal in channel 3 (750 nm). Each sample was imaged over 18 cycles, with 33 markers across three channels, excluding 4′,6-diamidino-2-phenylindole (DAPI), which was used in channel 0 (405 nm) for every cycle. Channel 1 (561 nm) markers included CD14, CD34, granzyme B, CD44, CD38, CD21, myeloperoxidase (MPO), CD19, TIGIT, CD8, and Ki67. Channel 2 (647 nm) markers included OX40 [tumor necrosis factor receptor superfamily, member 4 (TNFRSF4)], CD11A, CD152, S100A4, CD74, CD11C, CD68, inducible T cell costimulator (ICOS), CD4, forkhead box P3 (FOXP3), CD3E, LAG3, programmed death-ligand 1 (PDL1), CD34, TIM3, and CD11B. Channel 3 (750 nm) markers included galectin-9, CD57, CD31, CD47, CD20, and vimentin. The entire area of the BM biopsy specimen was considered in the analysis.

From six AML DLI Rs (two from the scRNA-seq discovery cohort and four from the flow cytometry validation cohort) and three NRs (one from the scRNA-seq discovery cohort and two from the flow cytometry validation cohort), we obtained peripheral blood after DLI (in remission for Rs) and isolated T cells from PBMCs (human Pan T Cell Isolation Kit, Miltenyi). A total of 1 × 10⁵ T cells were cultured with 1 × 10⁵ leukemia cells from peripheral blood or BM at the time of relapse after HSCT at an effector-to-target ratio of 1:1 or with T cells alone as a negative control. Equal numbers of cells were used for both Rs and NRs. After 18 hours of in vitro coculture, cells were stained for 20 min at 4°C in a cell staining buffer (BioLegend) with surface antibodies (CD34, CD3, CD4, CD8, CD62L, and CD45RA) (data file S7) at a 1:100 dilution. T cell functional state was assessed by IFN-γ catch assay according to the manufacturer’s instructions (Miltenyi). In a separate aliquot of cells from the same well, cells were stained with the same surface antibody cocktail as before along with CD137 to assess for antigen-specific T cell activation (70, 84). Cells were then washed three times before resuspending in a fluorescence-activated cell sorting (FACS) buffer (5% v/v FBS in PBS), and data were acquired using a BD Fortessa II instrument and then analyzed using the FlowJo (Tree Star) software. Cell subsets were defined as follows: CD8⁺ T_EMRA/CTL (CD3⁺CD8⁺CD62L⁻CD45RA⁺) and CD8⁺ non-T_EMRA/CTL (CD3⁺CD8⁺CD62L⁻CD45RA⁻ and CD3⁺CD8⁺CD62L⁺CD45RA^+/−) (fig. S9A). For coculture experiments using BM fibroblasts, cryopreserved BM samples for three R patients were thawed and plated on tissue culture flasks in a medium consisting of Dulbecco’s modified Eagle’s medium (Gibco) supplemented with 20% heat-inactivated FBS, 1% penicillin/streptomycin, 1% Hepes buffer, 1% nonessential amino acids, and 1% sodium pyruvate. The medium was changed the next day to remove nonadherent cells. Adherent cells (fibroblasts) were allowed to grow until confluent, with the medium changed every 3 days. On the day of the coculture experiments, cells were lifted with 0.05% trypsin and plated into three wells of a 96-well plate at a density of 1 × 10⁵ cells per well for each patient. Post-DLI T cells were isolated from peripheral blood by negative selection (Miltenyi) and added to wells containing fibroblasts (or without fibroblasts, for the T cell–only negative control). The coculture experiment was then performed as described above. For the single-cell analysis of CD137⁺ versus CD137⁻ sorted cells, peripheral blood T cells from three Rs in remission after DLI were obtained and cocultured with patient-matched leukemia cells using the same protocol described above. After 18 hours of coculture, cells were washed in a sterile FACS buffer, stained with CD8 and CD137 antibodies along with a viability dye (Zombie Aqua) for 20 min, and then washed three times and taken for flow cytometry sorting. Live cells were collected from the CD8⁺CD137⁺ fraction and the CD8⁺CD137⁻ fraction. There were too few cells after sorting to perform scCITE-seq staining. Cells were then washed once and resuspended in 0.04% BSA and then taken immediately for scRNA-seq and scTCR-seq using the 10X Genomics instrument and protocol described above.

Statistical and graphical analyses were performed using Python (3.8.3) and scipy packages (ttest_ind, mannwhitneyu, f_oneway, and kstest) and the GraphPad Prism software (version 10.0). ANOVA and paired t test were used for analyzing normally distributed data and comparisons with few data points. Mann-Whitney U and Kolmogorov-Smirnov tests were used for nonparametric data. When available, two-sided tests were used. For Fig. 6G and fig. S2F, Bonferroni correction was used to perform multiple-testing correction for determining statistical significance, with alpha set to 0.000001 and 0.001, respectively. For all other analyses, results were considered statistically significant for P < 0.05. For DIISCO analysis, the area represents 1 SD (16th to 84th percentile). For box plots of scRNA-seq data, boxes denote lower and upper quartiles, and whiskers show the full range of data. Outliers are shown as dots and were determined with the seaborn package using a method based on the interquartile range.