The study population was selected from 551 consecutive patients with allogeneic HCT between January 2017 and May 2020 at the Department of Hematology and Stem Cell Transplantation of the West-German Cancer Center, University Hospital Essen. Patients were screened for the following inclusion criteria: administration of in vivo T cell manipulation with 1) ATG or 2) PTCy as GVHD prophylaxis for HCT from 10/10 matched-unrelated (MUD) donors (CONSORT diagram, Supplementary Figure 1A ). All patients received peripheral blood stem cells (PBSC) as graft. Given that graft failure biases the analysis of donor-derived immune reconstitution, 9 patients with graft failure were excluded prior analysis. Patients transplanted with haploidentical donors treated with PTCy were included as comparators. A total of 384 patients were eligible for downstream analysis.

GVHD prophylaxis consisted of baseline calcineurin inhibitor-based immunosuppression combined with in vivo T cell manipulation using either ATG or PTCy. ATG (Grafalon®, Neovii, Rapperswill, CH) (n=304) was applied at a dose of 10mg/kg or 20mg/kg bodyweight on three consecutive days between day -4 and day -2 before HCT based on standardized protocols, followed by ciclosporin and methotrexate starting at day -1. PTCy (n=80) was administered on day +3 and +4 (50 mg/kg body weight per day) post-HCT followed by tacrolimus and mycophenolate-mofetil (MMF) starting on day +5. Out of these 80 patients receiving PTCy as GVHD prophylaxis 35 patients (44%) were transplanted with MUD donors and 45 patients (56%) with haploidentical donors. Early supportive and follow-up care followed the same internal protocols and was considered identical for all patients. Patients were followed-up until the last documented clinical assessment or death by any cause. Surviving patients were censored at maximum follow-up of 12 months. The clinical assessment standards are detailed in the Supplementary Methods .

Baseline data concerning patient-, donor-, allogeneic hematopoietic cell transplantation (HCT) characteristics and HCT-outcome were documented prospectively in electronic forms. Laboratory results and clinical characteristics of patients after HCT were retrospectively analyzed. Clinical assessment was obtained daily for inpatients and at each outpatient visit, starting with weekly intervals. Acute GVHD (aGVHD) was defined as GVHD organ involvement of skin, gut and/or liver until 100 days post-HCT. aGVHD was clinically assessed and classified according to the consensus aGVHD (24) grading. Diagnosis of chronic GVHD (cGVHD) starting from day +100 was based on characteristic symptoms and clinical signs according to the published NIH criteria (25). Overall survival was defined as the time from HCT to the end of the 12-months follow-up period or up to death by any cause. Cumulative incidence of relapse incidence (CIR) was calculated as the time from the day of transplantation to the day of documented relapse to original disease or persistence of malignancy. For patients without relapse or persisting malignancy, non-relapse mortality (NRM) was determined as the time from day of HCT to death.

Immune reconstitution after MUD HCT was studied in peripheral blood samples from patients around months +1, +3, +6, +9, and +12 after HCT. A total of 1297 samples were analyzed by flow cytometry at the BMT Laboratory, University Hospital Essen. For flow cytometry analysis, freshly-derived patient peripheral blood samples were prepared by isolating mononuclear cells (PBMC) using an automated red blood cell lysing system (TQ-Prep, Beckman Coulter, Brea, CA), washing with fluorescence-activated cell sorting (FACS) buffer and subsequently staining with surface markers ( Supplementary Table 1 ). No samples were cryopreserved. All samples were run on the same NAVIOS flow cytometer (Beckman Coulter, Brea, CA) using the same antibodies and FACS compensation parameters using the manufacturer’s software. Adequate subset representation was ensured by analysis of a minimum of 15000 lymphocytes in each run. For each flow cytometry sample two complementary antibody panels were used. The first panel characterized immune cell subsets as follows: T Cells, CD3⁺; T helper cells, CD3⁺/CD4⁺; cytotoxic T cells, CD3⁺/CD8⁺; regulatory T cells, CD3⁺/CD4⁺/CD25⁺/CD127⁺low; conventional T cells, CD3⁺/CD4⁺/CD25^-/CD127⁺high, not including the CD3⁺/CD4⁺/CD25^-/CD127^-low fraction; naïve helper T cells, CD3⁺/CD4⁺/CD45RA⁺; memory helper T cells, CD3⁺/CD4⁺/CD45RO⁺. Given that the panel did not cover CCR7 or CD62L, we described the CD8 ⁺ population including both naïve cytotoxic T cells and effector memory T cells re-expressing CD45RA (TEMRA) as CD45RA⁺ cytotoxic T cells, CD3⁺/CD8⁺/CD45RA⁺. Memory cytotoxic T cells were characterized by CD3⁺/CD8⁺/CD45RO⁺ and B cells by CD19⁺. These subsets were gated on the CD45⁺ lymphocyte gate, excepting the regulatory- and conventional T cells, which were selected from the CD3⁺/CD4⁺ subset in the CD45+ gate. In the second panel the following immune cell subsets were gated on the CD45⁺ lymphocyte gate: Activated T cells, CD3⁺/HLA-DR⁺; NKG2D⁺-NK cells, CD16⁺/CD56⁺/CD314⁺. T cell receptor α/β, TCRα/β and T cell receptor γ/δ, TCRγ/δ were gated on the CD3⁺ gate. T cell subset counts were determined as fraction of the absolute lymphocyte count on the sampling date. For each individual point in time ( Supplementary Figure 1B ), median counts of immune subsets were compared between MUD-ATG, MUD-PTCy using Mann-Whitney U test (GraphPad Prism 9.0.0, GraphPad Software, LLC, San Diego California). Samples from patients with haplo-HCT served as comparative samples to discriminate PTCy specific effects from those specific to the donor setting.

Detailed information about the methods to analyze individual patient’ longitudinal immune reconstitution is provided in the Supplementary Methods section. In short, we defined two distinct multi-dimensional immune cell clustering models integrating two different groups of T cell subsets: 1) “GVHD-associated” T cells: CD3⁺/CD4⁺/CD25⁺/CD127^low Treg, CD3⁺/HLA-DR⁺ activated T cells, TCRα/β⁺ and TCRγ/δ⁺ T cells based on immunologic evidence of a mechanistic impact in GVHD (26, 27) and 2) “broad spectrum” T cells: CD3⁺/CD4⁺ helper T cells, CD3⁺/CD4⁺/CD45RA⁺ naïve helper T cells, CD3⁺/CD8⁺ cytotoxic T cells, and CD3⁺/CD8⁺/CD45RO⁺ memory cytotoxic T cells which represent relevant subpopulations of the CD4⁺ and CD8⁺ T cell compartment (28, 29). In order to manage the complexity of the defined models both were limited to four T cell subsets. To apply these models, the study population (n=384) was filtered for patients with 1) at least three flow cytometry measurements within +12 months post-HCT and 2) measurement of the first flow cytometry ≤ d+45 post-HCT. Data filtering resulted in a patient subgroup of 180 patients eligible for clustering analysis (n=147 MUD-ATG, n=15 MUD-PTCy, and n=18 haplo-PTCy; for technical reasons, the “broad spectrum” T cell model included 4 more patients in the MUD-ATG cohort (n=151)). Each patients’ T cell subsets time-series underwent linear interpolation between datapoints ( Supplementary equation 1 ), calculated from adjacent values to have continuous data. Individual longitudinal immune reconstitution was then studied within each patient cohort by partitional clustering of time-series data (30, 31). Here, partitional clustering was performed using dynamic time warping (DTW) as distance measure (30) with 36 different function-specific configurations tested ( Supplementary Table 2 ). The performance of clustering configurations was evaluated by the silhouette coefficient (Sil), indicating a separation of clusters between -1 and +1 with the optimum at +1 (32). Model robustness was internally validated via a 10-fold resampling approach (see Supplementary method section) examining the variance of the Sil for each configuration. Data interpolation, DTW and time-series clustering were performed using R (33) packages R stats ( 33) and dtwclust ( 34) (R version 3.6.3, R core Team, https://www.r-project.org/). Patient clusters identified using this approach were evaluated for clinical outcomes as explained in the “Clinical statistical analysis” section.

Patient baseline characteristics were analyzed with Chi-square test and one-way ANOVA where appropriate (GraphPad Prism 9.0.0). The primary clinical outcome of this study was the incidence of grades II-IV aGVHD. Secondary endpoints were the overall incidence of 100-day aGVHD, 1-year relapse and NRM, 1-year cGVHD, as well as 1-year overall survival. The cumulative incidence of all-grade aGVHD and aGVHD II-IV was analyzed in a competing risk analysis considering death before d+100 as competing event and compared by Gray’s test. Complementary competing risk analysis was performed for cGVHD, which considered death within 12 months after HCT as competing event. Furthermore, the time-to-onset of all-grade aGVHD and aGVHD II-IV in the studied subgroups was calculated with the Kaplan-Meier method, obtaining event probabilities of time-to-event intervals. 1-year OS was analyzed via Kaplan-Meier analysis (35); subgroups were compared using the log-rank test; survival hazards were calculated by a Cox proportional hazards model (36). NRM and relapse were considered as competing events to each other and analyzed by competing risk analysis. P-values <0.05 were considered statistically significant. Clinical outcome analyses were done using the R (33) packages survival ( 37), survminer (38) and cmprsk (39) (R version 3.6.3, R core Team, https://www.r-project.org/).

Study protocol approval was obtained by the institutional review board of the University Duisburg-Essen (Protocols N° 17-7675-BO and N° 18-8299-BO). All patients have given written informed consent to collection, electronic storage, and scientific analysis of anonymized HCT-specific patient data in accordance with German legislation and the revised Helsinki Declaration. We confirm that no patient can be identified through use of anonymized patient data.

The HCT cohorts included in this study were balanced for age, sex, disease, graft source, conditioning and CMV recipient/donor serostatus ( Supplementary Table 3 ). The sex mismatch proportion was lowest for MUD-PTCY patients (18%). Median study follow-up was 12 months.

At the cohort level, the 100-day cumulative incidence of grades II-IV and III-IV acute GVHD (aGVHD) did not significantly differ between MUD-ATG or MUD-PTCy (II-IV: MUD-ATG 68.3% and MUD-PTCy 55.3%, p=0.224; III-IV: MUD-ATG 19.8% and MUD-PTCy 17.2%, p=0.848; Supplementary Table 4 , Figures 1A, C ). All-grade aGVHD was numerically lower in the MUD-PTCy cohort (p=0.07, Figure 1B; Supplementary Table 4 ) and occurred relatively later (Median time 21 vs 17 days, p=0.058, Supplementary Figure 2B ). Fine and Gray competing risk regression corroborated these results with lower all-grade aGVHD subdistribution hazards for MUD-PTCy (SHR 0.69, 95%CI 0.49-0.97, p=0.032, Supplementary Table 5 ), which did not retain significance for grades II-IV aGVHD (SHR 0.74, 95%CI 0.45-1.21, p=0.220, Supplementary Table 6 ). Both, all-grade chronic GVHD (cGVHD) and moderate-severe cGVHD at 12 months after HCT were comparable between both MUD cohorts (p=0.207 and p=0.452, Figures 1E, D ). However, the incidence of severe cGVHD was numerically lower in the MUD-ATG cohort ( Figure 1F ). Overall survival (OS), non-relapse mortality (NRM) and the cumulative incidence of relapse (CIR) until 12 months did not differ between cohorts ( Figures 1G-I , Supplementary Tables 4 , 5 ). Clinical outcomes with haplo-PTCy are detailed in Supplementary Tables 5 – 7 .

Given the established role of T cells as aGVHD initiators, we investigated potential differences in immune reconstitution to detect alternative fractions involved in immune modulation by ATG or PTCy as T cell depletion and manipulation strategies. Indeed, despite the observed similar clinical efficacy ( Supplementary Table 4 ) of both protocols, comparative analysis of cellular immune reconstitution revealed significant differences in T cell subsets. Throughout the first 12 months after HCT, ATG patients had significantly lower absolute counts within the helper T cell compartment ( Figures 2A–E ) compared to patients receiving PTCy. This pattern was also observed for TCR α/β T cells up to month 6 ( Figure 2G ). Interestingly, the absolute Treg counts were also significantly higher for the first 6 months after HCT in patients receiving PTCy as T cell manipulation ( Figure 2D ) compared to the MUD-ATG cohort. The Treg/Tcon ratio did not differ between the cohorts. In PTCy patients irrespective of the distinct donor type, early immune reconstitution up to month 6 was comparable confirming a PTCy specific benefit to the helper T cell compartment. Despite this early comparability of PTCy cohorts, median absolute counts of several helper T cell subsets stagnated in MUD-PTCy patients between months 6 and 12 ( Figure 2 ) indicating a donor-type specific impact on helper T cell expansion after month 6. Contrary to the overall T cell cytopenia of the MUD-ATG cohort, its median TCR γ/δ T cell counts ( Figure 2H ) were significantly higher compared to MUD-patients receiving PTCy. The analysis of γδ T cell reconstitution further revealed broad confidence intervals for all analyzed subgroups, not originating from limited patient numbers as this was also observed for the large MUD-ATG cohort. Interestingly, the CD8⁺ subsets ( Figures 2I, J ) were not significantly affected by in vivo T cell depletion with ATG and did not account for differences between the MUD-ATG and MUD-PTCy cohorts. Of note, early NKT cell counts were also significantly higher in the MUD-ATG cohort (months 1 and 3, Supplementary Figure 3F ). The reconstitution of further subsets e.g. cytotoxic- and activated T cells subsets ( Supplementary Figures 3B-E ), as well as NK- and B cells ( Supplementary Figures 3G–H ) was comparable between MUD-PTCy and MUD-ATG patients. Despite their early increase, CD3⁺ T cell levels in MUD-PTCy patients declined between months 6 and 12 leading to CD3⁺ numbers comparable to MUD-ATG patients at the end of the observation period. This finding was consistent throughout the majority of T cell subsets, with little exceptions, equalizing the above-described early differences in immune reconstitution of MUD patients between both T cell depleting regimens in the long-term.

Following the observation of broad confidence intervals in the pooled immune reconstitution data (e.g., in TCRγ/δ) and clinical outcomes, we hypothesized a relevant interindividual heterogeneity within each study cohort. Therefore, we analyzed cellular recovery with a data science approach able to dissect such heterogeneity within cohorts and to identify patients with similar reconstitution patterns, and possibly also homogenous clinical outcomes. Here, we developed two multi-dimensional parameter models integrating longitudinally measured reconstitution data of different T cell subsets for each patient. Because of the included cell types, we coined these models the “GVHD-associated”- and “broad spectrum” T cell model. Within each patient cohort (MUD-ATG, MUD-PTCy, and haplo-PTCy), both multi-dimensional models revealed distinct patterns by time-series clustering and dissected intra-cohort heterogeneity. Based on actual cell counts and reconstitution shapes, the models produced clusters, which, when correlated to clinical outcomes, revealed differences. The methodological workflow is detailed in the Supplementary methods and illustrated in Figure 3A . OS, NRM and CIR of patients included in this multi-dimensional analysis (n=180) were representative of patients with survival beyond d+100 in the overall cohort ( Figures 3B–D ), making a selection bias unlikely. Time-series clustering of T cell subsets from the “GVHD-associated” model dissected the MUD-ATG cohort into distinct patient subgroups (labelled patient clusters, Figures 3E, F ). We determined the optimal model configuration via both a good and robust silhouette coefficient (configuration 1_1: Sil=0.524) as well as a balanced patient distribution (The MUD-ATG subgroup (n=147) was split into cluster 1: n=94 and cluster 2: n=43; Figures 4A–C ). The cell subsets that contributed most to this clustering were activated- and αβ T cells, because they revealed greater differences in shape and higher absolute counts in both clusters over time as compared to Tregs and γδ T cells. Both of these T cell subsets were comparatively illustrated at smaller scales between the cohorts in Supplementary Figure 6 . The corresponding cluster centroids ( Figure 3F ) confirmed distinct reconstitution shapes of each cluster. Patients in cluster 2, which had higher absolute counts of activated- and αβ T cells compared to cluster 1, had significantly lower NRM (p=0.032, Figure 3H ) and relapse (p=0.01, Figure 3I ) resulting in higher 1-year OS (98% vs. 79%, p=0.0023, Figure 3G ). Next, we leveraged this cluster information comparing clinical outcomes between ATG clusters and the PTCy cohorts. Interestingly, patients from both the MUD-ATG cluster 1 and the MUD-PTCy cohort had similarly significantly decreased OS compared to MUD-ATG cluster 2 and haplo-PTCy (p=0.0053, Figure 5A ). NRM and relapse (p=0.077 and p=0.057, respectively, Figures 5B, C ) were numerically increased in MUD-ATG cluster 1. Despite a numerically lower incidence of aGVHD II-IV in MUD-PTCy patients ( Figure 5D ), their cGVHD incidence was comparable to MUD-ATG patients from cluster 1 ( Figure 5E ), which were both quantitatively higher compared to MUD-ATG cluster 2 (p=0.061, Figure 5F ). We also tested the method of time-series clustering on the smaller PTCy cohorts using the “GVHD-associated” T cell model ( Figures 5G–J ). Again, the clustering successfully dissected intra-cohort heterogeneity in immune reconstitution ( Figures 5G–J ), although silhouette coefficients for optimal configurations were lower compared to the larger MUD-ATG cohort ( Figures 4D–I ). Similar to the results in the MUD-ATG cohort, the clustering of PTCy patients successfully distinguished two clusters in each setting characterized by relevant differences in absolute counts and reconstitution shape of activated- and αβ T cells ( Figures 5G–J ). A second clustering model integrating different “broad spectrum” T cell subsets, reflecting the pattern of CD4⁺ and CD8⁺ T cell reconstitution, also identified two separate clusters in the MUD-ATG cohort. Its optimal cluster configuration had a comparable silhouette coefficient (configuration 2_1: Sil=0.536) to the “GVHD-associated” T cell clustering model, and both good robustness and an appropriate patient distribution (cluster 1: n=105 and cluster 2: n=46; Supplementary Figures 4A–C ). Here, the revealed clusters were most influenced by cytotoxic and memory cytotoxic T cells, which were higher in cluster 2 ( Figures 6A, B ). Strikingly, although both models integrated biologically different T cell subsets their degree of similarity was 92.5% as revealed by cluster model comparison ( Figure 6C ). Patient re-allocation between the two models was minimal (n=11, 7.5%, Figure 6D ) and baseline characteristics were similar between the identified clusters e.g. for patient age or underlying disease ( Supplementary Tables 8 , 9 ). Consequently, clinical analysis of these MUD-ATG clusters in the “broad-spectrum” T cell model revealed analog results to the “GVHD-associated” T cell model ( Figures 6E–G ), with the exception of a marginal significance for NRM. As reported above, the clustering of PTCy patients in the “broad spectrum” T cell model, also showed overall lower and more unstable silhouette coefficients ( Supplementary Figures 4D–I ). Again, cytotoxic T cell compartments appear to be the leading contributors to PTCy clusters ( Supplementary Figures 5A–D ). Given, the previously described impact of CMV serostatus on the T cell repertoire (40) we next examined the association between CMV serostatus and the time-series clustering. While the identified clusters in the configuration 2_1 had significantly distinct proportions of patients with R+ CMV serostatus, the differences in T cell kinetics between those with CMV R-/D-, CMV R-/D+ or CMV R+/D- serostatus were minimal ( Figures 7A–D ). However, the CMV R+/D+ subpopulation revealed important differences in T cell kinetics. We consequently tested if the clustering was further able to dissect T cell reconstitution within CMV R+/D+ patients and successfully applied the same cluster configuration 2_1 identifying relevant differences in T Cell kinetics ( Figures 7E, F ) as well as in OS (p=0.051, Figure 7G ) within the CMV R+/D+ subgroup. Next, we tested if a combination of the most relevant T cell subsets from both the “GVHD-associated” and the “broad spectrum” T cell model would improve the clustering results. Indeed, the combination of CD3⁺/CD8⁺/CD45RO⁺, αβ T cells, HLA-DR⁺ T cells and CD3⁺ T cells ( Supplementary Figure 7 ) associated with significant differences in clinical outcome and with slightly better separated OS curves (p=0.0014, Supplementary Figure 7C ). Finally, we also tested a non-linear pre-processing approach instead of linear interpolation. The resulting curves had a similar, yet smoothened shape ( Supplementary Figures 8A, F ) and clinical outcome association was comparable ( Supplementary Figures 8C–E ). Taken together, both time-series clustering models integrating distinct T cell subsets were able to dissect intra-cohort heterogeneity in post-HCT immune reconstitution and to identify relevant patient subsets with distinct outcome.