The ZUMA-7 trial was conducted at 77 sites worldwide. Eligible patients were at least 18 years of age and had histologically confirmed large B-cell lymphoma, according to the World Health Organization 2016 classification criteria, that was refractory to first-line treatment or that had relapsed from complete remission no more than 12 months after the completion of first-line chemoimmunotherapy; patients intended to proceed to high-dose chemotherapy with autologous stem-cell transplantation. Refractory disease was defined as a lack of complete response to first-line therapy, and relapsed disease occurring no more than 12 months after the completion of first-line therapy.

All the patients provided written informed consent. The trial was conducted in compliance with the principles of the Declaration of Helsinki. Full details around trial design and oversight were previously reported (1, 2).

The primary end point was event-free survival (defined as the time from randomization to the earliest date of disease progression according to the Lugano classification, the commencement of new therapy for lymphoma, death from any cause, or a best response of stable disease up to and including the response on the day 150 assessment after randomization) according to blinded central review. Secondary end points included progression-free survival (defined as the time from randomization to disease progression or death from any cause), assessed centrally or by investigator and overall survival, durability of response (defined as the time between the first objective response to disease progression and the start of new lymphoma therapy), and the incidence of adverse events, as previously reported (1, 2, 4).

Endpoints utilized for the analyses reported in this manuscript are EFS, PFS, and DOR per central review from the primary EFS analysis data cut, as well as PFS per investigator assessment (as indicated in the Result section and figures) from primary overall survival analysis, which occurred 5 years after the first subject was randomized, previously published (1, 2).

Classification of patients into molecular subtypes of cell of origin (GCB vs. non-GCB), double/triple hit [high-grade B-cell lymphoma (HGBL)], MYC rearrangement, and double-expressor lymphomas (overexpression of MYC and BCL2 proteins) was based on central laboratory analysis and previously reported (1, 3).

Details regarding ethics, patient samples, and efficacy endpoints were described previously (3). In this post hoc analysis, evaluable samples collected pretreatment (axi-cel, n=134; SOC, n=122) or at progression, post axi-cel treatment (n=17), from patients in ZUMA-7 study were analyzed. Gene expression profiling was performed using the NanoString nCounter^® PanCancer IO 360™ Panel. RNA extraction and processing into the nCounter platform for the NanoString IO360 assay was conducted at NeoGenomics. Nanostring RCC and RLF files were imported on nSolver Analysis software (v.4.0). Raw data were further analyzed with nCounter Advanced Analysis (v.2.0.134), and normalized linear count output were used for all further analysis, as previously described (3). Predefined IO-360 NanoString signatures were also analyzed.

To identify genes whose transcripts were associated with clinical outcomes including duration of response (DOR), event-free survival (EFS), and progression-free survival (PFS), penalized Cox proportional hazards models were utilized (5). These multivariate models consider the expression of multiple genes simultaneously to select key predictors while preventing overfitting. The elastic net penalty was used to combine the feature selection capability of Lasso (L1 penalty) with the stability of Ridge (L2 penalty), improving the robustness of biomarker identification, as per below formula:

where γ was set to 0.9, and the optimal α was determined using fivefold cross validation based on concordance index (C-index), measuring predictive accuracy. Models were fitted individually for three clinical outcomes (DOR, EFS, PFS), aiming to identify overlapping transcripts to ensure a robust final signature. Non-zero coefficients were identified from the best-performing model to derive gene signatures for each clinical outcome. Transcripts/genes with negative coefficients (hazard ratio [HR]<1) were termed favorable transcripts (associated with improved outcome), whereas genes with positive coefficients (HR >1) were termed unfavorable transcripts (associated with worse outcome).

A gene expression signatures (GES) score was assigned to each sample using the mean scaled gene expression with the following formula:

where m is the number of genes within a specific signature, C o u n t i is the normalized linear count from NanoString or linear TPM value when derived from RNA-Seq (see below about technical replication) associated with g e n e i within the given signature for the given sample, u i   is the mean normalized count or TPM value calculated from all samples for g e n e i , and σ i is the standard deviation of normalized count or TPM values derived from all samples for g e n e i .

The patients were subsequently divided into two groups based on their scores relative to the median. Those with scores above the median were placed in the “high” group, whereas those with scores at or below the median were assigned to the “low” group. This median split method ensures an even distribution of participants across the two groups, facilitating a balanced comparison in subsequent analyses.

The predictive value of the identified signatures was replicated using a subset of ZUMA-7 samples profiled by RNA sequencing (RNA-Seq; axi-cel, n=124; SOC, n=125). Most patient samples overlapped between the NanoString gene expression and RNA-Seq analysis conducted herein (axi-cel, n=119; SOC, n=115).

Gene expression profiling was carried out from FFPE tumor samples by performing next-generation sequencing (NGS). Total RNA libraries were prepared using Illumina Stranded Total RNA Prep, Ligation with Ribo-Zero Plus kit (catalogue # 20040529). This kit utilizes the enzymatic rRNA depletion and ligation-based addition of adapters and indexes for whole transcriptome libraries, capturing coding and non-coding RNAs.

Following Illumina’s library prep protocol, 100 ng or 10 ng of total RNA, depending on the available amount, was used for rRNA depletion followed by fragmentation, cDNA synthesis, and library amplification. Resulting libraries were quantified using Qubit, and quality control was performed on TapeStation 4150 (Agilent). All libraries were diluted to equimolar concentration and pooled for sequencing on Illumina’s NovaSeq 6000 with a read length of 101-bp paired end. The BBDuk tool (6) was applied to trim the adapters and the low-quality reads. Pair-ended RNA-seq reads were mapped and quantified using Salmon v1.9.0 against the GENCODE V43 human transcriptome (https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_43/gencode.v43.transcripts.fa.gz). Gene-level counts and transcripts per million (TPM) values were generated using Tximeta v1.16.1. TPM values were used for downstream analysis. A total of 311 samples had available transcriptional data from either screening or post-infusion time point. RNA-sequencing quality control filtering was applied, removing samples with ≤10% reads mapped to coding sequences, leaving 296 samples, 277 pre-infusion samples, and 19 post-infusion samples across both arms of ZUMA-7 (axi-cel and SOC). Samples from patients who were not part of the safety subset (patients enrolled in the study but not treated) were excluded, resulting in 266 evaluable pre-infusion samples. To remove duplicated subjects, for those with multiple samples collected prior to axi-cel treatment, the sample with the latest collection date was kept. If collection dates were the same, the sample with the lower Molecular Batch ID (sample analyzed first) was retained. After duplicate removal, 249 pre-infusion samples remained available for the downstream analysis (124 from the axi-cel arm and 125 from the SOC arm of ZUMA-7).

To determine the association between the signatures identified in this study and outcome to first-line R-CHOP or -CHOP-like treatment, we utilized two publicly available datasets (7, 8). The dataset from Schmitz et al. (8) included 229 patients who received immunochemotherapy (R-CHOP or R-CHOP-like therapy), whereas 522 patients who received R-CHOP as initial therapy were selected from the Reddy et al. dataset (7). The linear-scale TPM from the online dataset were scored into signature following the same formula/method described above. The available clinical outcome readout from the literature was PFS.

For this analysis, the RegParallel (version 1.12.0) R package (https://github.com/kevinblighe/RegParallel) was utilized to examine how the expression of genes influenced the rate of progression-free survival. The function RegParallel was run to fit the Cox proportional hazards regression model to gene expression to independently test the association between survival time and each gene. The output was analyzed with TIBCO Spotfire (version 11.4.3).

All analyses described herein were exploratory and retrospective, performed without adjustment for possible confounding factors. Log−rank test was used to determine statistical significance between the high and low groups in time-to-event analyses presented as KM curves, where high means > median value and low means ≤ median value. When further subgrouping was performed (e.g., patient segregated by the GES median threshold and further into GCB and non-GCB subsets), the median GES values were calculated from the entire dataset. Where indicated, an unstratified Cox regression model was used to provide an estimate of the hazard ratio and associated P-value. Wilcoxon rank-sum test was used when comparing two groups in categorical analysis. Spearman’s rank-order correlation was used to evaluate the association between variables. For these post hoc analyses, all P values were descriptive and P < 0.05 was considered significant. No adjustments for multiplicity testing were performed.

Transcriptomic analysis, generated via Nanostring nCounter technology, aimed to identify robust gene expression signatures predictive of three clinical outcomes, duration of response (DOR), event-free survival (EFS), and progression-free survival (PFS), following second-line axi-cel therapy. Penalized Cox regression models were fitted individually for each of the three endpoints to identify gene transcripts with non-zero coefficients (see methods).

Such analysis identified two distinct gene transcript lists, encompassing transcripts favorably or unfavorably linked to the three clinical outcomes (DOR, EFS, and PFS) in patients receiving axi-cel, as represented in Figure 1A . The full list of transcripts with the non-zero coefficient for each outcome is provided in Supplementary Table 1 , split into two tabs with favorable or unfavorable transcripts. A list of six transcripts, herein referred to as 6-GES, representing the intersection of favorable transcripts for each outcome included CD19, CD45RA, CCL22, KLRK1, SIGLEC5, and SOX11. Conversely, a list of 17 transcripts, herein referred to as 17-GES, representing the intersection of unfavorable transcripts, included IL18R1, GPC4, KIR3DL2, ITGB8, PSMB5, RPS6KB1, BCL2, TNFSF4, SERPINA9, DUSP5, NBN, GLUD1, ESR1, CD45RO, ARID1A, KLRB1, and SLC16A1.

The gene expression values of the transcripts included in these two lists (6-GES or 17-GES) were then used to generate a scaled mean signature value of either 6-GES or 17-GES for each patient.

The association of the 6-GES and 17-GES with efficacy readouts was evaluated by stratifying patients into “high” and “low” groups based on the median signature scores. Kaplan–Meier curves illustrate the significant associations (descriptive P< 0.05) between the 6-GES or 17-GES and time-to-event efficacy readouts (DOR, EFS, or PFS; Figure 1B ). High expression of the favorable 6-GES correlated with improved DOR (HR: 0.29, 95% CI: 0.17–0.52), EFS (HR: 0.27, 95% CI: 0.16–0.44), and PFS (HR: 0.27, 95% CI: 0.16–0.46) in axi-cel-treated patients, whereas high levels of the unfavorable 17-GES negatively correlated with DOR (HR: 7.59, 95% CI: 3.95–14.60), EFS (HR: 6.12, 95% CI: 3.57–10.50), and PFS (HR: 7.47, 95% CI: 4.11–13.57) ( Figure 1B ; P< 0.0001 for all readouts). These associations were specific to the axi-cel arm of ZUMA-7 and not observed in patients receiving standard-of-care (SOC) treatment ( Figure 1C ), underlying the potential predictive value of these signatures in the axi-cel treatment setting. These signatures also strongly associated with PFS per investigator assessment ( Figure 2A ).

Figure 2A presents the univariate associative analysis of 6-GES or 17-GES or each individual transcripts or predefined GES from the NanoString IO-360 panel with PFS outcome (per investigator), following axi-cel treatment in ZUMA-7: as shown, 6-GES and 17-GES outperformed all individual genes and predefined GES for their association with PFS ( Figure 2A ). Notably, of the 6 or 17 transcripts identified by the multivariate analysis based on the non-zero coefficient around DOR, EFS, and PFS readouts, only CD19, CD45RA, ESR1, and DUSP5 were significantly associated with PFS in univariate analysis, consistent with the notion that 6-GES and 17-GES can carry predictive value not recapitulated by individual transcripts. Figure 2 also shows that there could be added value in combining 6-GES and 17-GES: the two GES significantly (P = 0.01), but modestly (Spearman R = −0.16) correlated (negatively) with each other ( Figure 2B ; Supplementary Figure 3 ), and subgrouping based on a combination of high and low 6-GES and 17-GES could further improve risk stratification of axi-cel patients ( Figure 2C ). Supplementary Figure 1 shows a correlation matrix of the 6-GES and 17-GES with the predefined IO360 GES (from NanoString) and the stromal and immune-suppressive index (SII), previously reported (3). This shows that 6-GES correlated with the favorable B-cell signature and represents a less immune infiltrated TME, whereas 17-GES correlated with a more complex and immune-infiltrated TME, including enrichment of the SII (3).

The predictive value of 6-GES and 17-GES was further evaluated in patients treated with axi-cel, stratified by the cell-of-origin subtype. As shown in Figure 2D , the Kaplan–Meier curves illustrate the PFS based on GES and cell-of-origin (COO) classification, defined as GCB vs. non-GCB (ABC + unclassified) (3). In the germinal center B-cell (GCB) subgroup, patients with a high 6-GES had significantly (P< 0.05) longer PFS compared to those with a low 6-GES. The 17-GES was also strongly associated with PFS in the GCB subgroup, with patients having a low 17-GES experiencing better PFS outcomes. Among non-GCB patients, the 17-GES remained associated with PFS, whereas the 6-GES did not significantly stratify patients in this subgroup, albeit there was a trend (P=0.088) with a clear separation between the KM curves. Consistently, Supplementary Figure 2 shows that the 6-GES and 17-GES do not associate with COO (GCB vs. non-GCB; Supplementary Figures 2A, B ). The 6-GES is also not associated with double-/triple-hit status (high-grade B-cell lymphoma (HGBL)), double-expressor (overexpression of MYC and BCL2 proteins), or MYC rearrangement status, whereas the 17-GES was lower in the HGBL and MYC rearrangement subgroup, compared to LBCL not otherwise specified (not applicable indicates LBCL not belonging to the other molecular subgroups; Supplementary Figures 2C, D ).

To corroborate the association between the favorable 6-GES and unfavorable 17-GES with axi-cel efficacy in ZUMA-7, we analyzed an RNA-seq dataset derived from a subset of ZUMA-7 tumor samples. A total of 21 transcripts overlapped between RNA-seq and NanoString quantification; CD45RA and CD45RO (CD45 isoforms) were excluded from the RNA-Seq analyses because RNA-seq transcript quantification was performed at the gene level to retain data robustness. Hence, in the context of RNA-Seq datasets, a 5-GES and a 16-GES were generated, consisting of the 6-GES minus the CD45RA transcript and 17-GES minus the CD45RO transcript. Spearman correlations across NanoString- and RNA-Seq-derived GESs is shown in Supplementary Figure 3 . The 5-GES and the 16-GES generated from RNA-Seq strongly associated with efficacy endpoints, including DOR, EFS, and PFS ( Supplementary Figure 4 ). Patients with a high 5-GES value exhibited improved DOR, EFS, and PFS compared to those with a low 5-GES value (P< 0.05). Conversely, patients with a high 16-GES score demonstrated inferior outcomes across all three endpoints relative to those with a low 16-GES score.

The ability to recapitulate these findings in a dataset generated through an orthogonal technology reinforces the robustness of the favorable 6-GES (or 5-GES for RNA-Seq) and unfavorable 17-GES (or 16-GES for RNA-Seq) as potential predictive biomarkers of outcome following axi-cel therapy in R/R LBCL patients in the second-line setting.

Next, we evaluated the association of the 5-GES (6-GES without CD45RA) or 16-GES (17-GES without CD45RO) with PFS outcome in two publicly available LBCL RNA-Seq datasets where patients were treated with 1^st Line R-CHOP or R-CHOP-like therapy and have corresponding PFS data; for this, we utilized the datasets from Schmitz et al. and Reddy et al. (7, 8). The dataset from Schmitz et al. (8) included 229 patients, whereas the dataset from Reddy et al. included 522 patients (7).

The 5-GES and 16-GES while being associated with outcome to axi-cel in ZUMA-7 (P< 0.05; Supplementary Figure 4 ) did not show an association with PFS to 1^st Line R-CHOP/R-CHOP like ( Supplementary Figure 5 ), indicating their predictive (specific to CAR T-cell therapy), rather than prognostic, value.

The 6-GES and 17-GES were analyzed in evaluable patients (available tumor biopsy) who progressed following axi-cel treatment (n=17). The boxplots in Figure 3 depict the distribution of 6-GES and 17-GES scores before treatment (n=256) and at disease progression following axi-cel treatment. At the time of progression, the favorable 6-GES score was numerically lower, albeit not significant by descriptive statistics (P=0.082), whereas the unfavorable 17-GES was significantly increased (P=0.031), compared to pretreatment levels. These findings corroborate that the 6-GES and 17-GES are likely representative of response and resistance mechanisms and further support previous observations (9) that the TME and the biomarkers associated with outcome can change through lines of therapy.