Evolutionary and immune microenvironment dynamics during neoadjuvant treatment of oesophagael adenocarcinoma

Locally advanced oesophageal adenocarcinoma (EAC) remains difficult to treat because of common resistance to neoadjuvant therapy and high recurrence rates. The ecological and evolutionary dynamics responsible for treatment failure are incompletely understood. Here, we performed a comprehensive multi-omic analysis of samples collected from EAC patients in the MEMORI clinical trial, revealing major changes in gene expression profiles and immune microenvironment composition that did not appear to be driven by changes in clonal composition. Multi-region multi-timepoint whole exome (300x depth) and paired transcriptome sequencing was performed on 27 patients pre-, during and after neoadjuvant treatment. EAC showed major transcriptomic changes during treatment with upregulation of immune and stromal pathways and oncogenic pathways such as KRAS, Hedgehog and WNT. However, genetic data revealed that clonal sweeps were rare, suggesting that gene expression changes were not clonally driven. Additional longitudinal image mass cytometry was performed in a subset of 15 patients and T-cell receptor sequencing in 10 patients, revealing remodelling of the T-cell compartment during treatment and other shifts in microenvironment composition. The presence of immune escape mechanisms and a lack of clonal T-cell expansions were linked to poor clinical treatment response. This study identifies profound transcriptional changes during treatment with limited evidence that clonal replacement is the cause, suggesting phenotypic plasticity and immune dynamics as mechanisms for therapy resistance with pharmacological relevance.


Introduction
Oesophageal cancer is the sixth most common cause of cancer-related death worldwide, with a median overall survival of <1 year 1 . Incidence rates for oeshophageal adenocarcinoma (EAC) have risen sharply and it is now the predominant subtype in high-income countries 2 . Patients with locally advanced EAC are treated with neoadjuvant chemotherapy (CTx) or radiochemotherapy (RCTx) followed by surgical resection 3 . Although neoadjuvant treatment confers a survival bene t over resection alone, 50-60% of tumours are resistant to neoadjuvant therapy, leading to an overall poor outcome with a 5-year survival of 12.6% 3 . There is currently no way to determine who will bene t from neoadjuvant treatment.
Genetic sequencing studies have examined the genetic alterations in EACs during treatment in order to understand evolutionary dynamics that may lead to treatment resistance [4][5][6] . The translational goal of these studies is to translate measures of the evolutionary dynamics into treatment-predictive biomarkers.
However, there is growing evidence that genetic evolution alone does not fully explain resistance evolution 7 . Tumour plasticity, de ned as cellular phenotype changes in the absence of underlying genetic change, is associated with treatment resistance in multiple cancer types [8][9][10] . Furthermore, there is mounting evidence that the e cacy of chemotherapies also relies on activating antitumor immune responses 11 . Measures of the tumour microenvironment (TME) are an important predictor for the response to neoadjuvant treatment and for future perspectives in regard to clinical trials e.g. KEYNOTE-585, which evaluates combined checkpoint and chemotherapy in the perioperative setting in EAC patients 12 .
There is a major practical challenge in collecting longitudinal samples over space and time from solid tumours. Here we have performed multi-region, multitimepoint whole exome sequencing (WES) at high depth (mean 300x) of locally advanced EAC to identify mutations and track clonal dynamics across time and space, with matched RNA-sequencing to characterize phenotype changes, and imaging mass cytometry (IMC) and T-cell receptor (TCR) sequencing to characterize immune in ltrates and stromal cell dynamics. Our comprehensive analysis provides a holistic multi-omic view of the dynamic changes in the tumour and its microenvironment through neoadjuvant treatment.

Results
Longitudinal multi-omic pro ling in 10 Non-Responder EAC and 17 Responder EAC EAC specimens were obtained from the prospective clinical DKTK-funded MEMORI trial 13 . Patients underwent baseline (18)F-FDG PET followed by 1 cycle of neoadjuvant platin-based chemotherapy, and clinical response was assessed at day 14-21 by (18)F-FDG PET. Responders (REs) to the initial treatment continued to receive chemotherapy, whereas non-responders (NRs) were switched to intensi ed RCTx (for detailed information see methods part 1). EAC shows limited evidence for strong clonal selection during neoadjuvant treatment We assessed genetic changes (SNVs and CNAs) in responders and non-responders (Fig 1A, B and C). We observed an average mutational burden of 23.09 SNVs/Mb (range 7.5 -127.7 SNVs/Mb) in treatment naïve EAC, which is consistent with previously described mutation burden in EAC, accounting for our higher depth of coverage 14,15 .
The overall number of SNVs between pre-and post-treatment samples showed no signi cant difference (Fig 2A), but we noted that two NRs showed large increases in SNV burden after radiation ( Fig 2B) suggesting radiation-induced mutagenesis. Further, we assessed changes in mutation clonality. If treatment drove clonal selection, we might expect to see increased clonal and/or subclonal mutational burden between pre-and post-treatment. However, no signi cant difference was detected for either mutation class (Fig 2C).
To probe clonal dynamics, we constructed phylogenetic trees using SNVs for 15 patients with multi-timepoint and multi-region samples (≥ 3 samples) (Fig 2D, supplemental g 1B) 16 . The majority of phylogenetic trees were 'balanced' with similar clade lengths for all samples. We observed progressive clonal sweeps during neoadjuvant treatment in only one sample (NR9), wherein the clade from Timepoint C arose from B, and B arose from Timepoint A (Fig2D, supplemental g 1B). In all other cases, samples from timepoint C were either most similar to precursor clones of samples from timepoints A and B (e.g. NR21 and RE22), or most similar to a clone detected at timepoint A (e.g. RE10). A parsimonious explanation of these observations is that the cancer was a mosaic of clones, and inferred clonal relationships between timepoints and/or spatial samples were determined by spatial sampling rather than by a clonal sweep driven by treatment.
Patients NR21 ( Fig 2D) and NR22 (supplemental g 1B) who accrued many mutations during RCTx were outliers with long branches for timepoint C (resection post-treatment), potentially indicating the outgrowth of a highly mutated subclone. Indeed, we observed a highly fragmented genome in our copy number analysis post RCTx ( Fig 3A) and a dN/dS analyses to detect evidence of clonal selection indicated RCTx-induced random DNA damage (dN/dS<1) rather than SNVs undergoing positive selection (supplemental table 5).
EAC driver mutations as listed by IntOGen ©17, 18 and neoantigenic SNVs were annotated on each phylogenetic tree ( Fig 2D and supplemental g 1B). Known high-frequency driver events such as TP53 and CDKN2A were mostly truncal on the phylogenetic trees, indicating that these drivers persisted through treatment. TP53 was mutated by SNVs in 17 patients (Fig 2D, 3E, supplemental g 1B) and indels in 4 out of 27 patients (Fig 2D, 3E, supplemental g 1B) from which 13 out of 22 patients with ≥ 2 samples showed clonal TP53 alterations present in all samples. CDKN2A showed SNVs in 2 patients (Fig 3E, supplemental g 1B) and harboured indels in 5 patients (Fig 2D, 3E) with clonal genetic CDKN2A alterations in 5 out of 22 patients with ≥ 2 samples (Fig 2D,   3E), consistent with prior reports on EAC genomes 14,15,19,20 . In our phylogenetic trees (Fig2D, supplemental g 1B), additional non-silent mutations in known driver genes events were detected at the end of neoadjuvant treatment in 7 out of 12 patients with sequenced samples from pre-and posttreatment timepoints (timepoint A and C present). These new putative driver mutations were high frequency driver (genes that are detected as drivers in >5% of samples in Intogen cohorts) in only 1 patient (CDKN2A in NR22). Therefore, 11/12 patients with sequenced samples from pre-and posttreatment timepoints had either no new putative drivers post-treatment, or low-frequency drivers, suggesting no major changes in the clonal make-up for the large majority of patients during neoadjuvant treatment.
Distinct treatment regimens have been associated with speci c mutational footprints 21 . We assessed the speci c mutagenic consequences of therapy on our cohort. SNVs from our phylogenetic trees were classi ed into treatment naïve SNVs, SNVs acquired during chemotherapy exposure (CTx-SNVs) and SNVs acquired during RCTx exposure (RCTx-SNVs). In CTx-SNVs a signi cant increase in C>A mutations was observed compared to treatment-naïve SNVs (Wilcoxon rank-sum test P-value=0.03) (Fig 2E), which is consistent with previously described platin-induced mutation signatures 4,5,21 . RCTx-SNVs showed a signi cant increase in T>C (Wilcoxon rank-sum test P-value=0.009) (Fig 2E). Across subgroups from different treatment response and timepoints, we detected nine prevalent COSMIC signatures (Fig 2F), of which eight have been previously identi ed in oesophageal cancer 6,[22][23][24] . Signature 12, which had not been previously described in EAC, was at very low activity in treatment naïve samples, but increased during RCTx and might therefore represent a treatment induced rather than EAC speci c signature. REs showed an non-signi cant increase in S4 signature during treatment (Wilcoxon rank-sum test p-value=0.036 and p-adj=0.34 between Timepoints A and C), which is dominated by C>A changes and thus most likely oxaliplatin driven 21 (Fig 2F+G). NRs depicted a nonsigni cant increase of radiotherapy associated S5 signatures 25 after exposure to RCTx (Wilcoxon rank-sum test p-value=0.14 and p-adj=0.6 between Timepoints A and C) (Fig 2F+H).
Copy number evolution through treatment EACs genomes are dominated by a high burden of CNAs in the majority of cases 15,22,26 . Our WES data from 17 REs and 10 NRs allowed us to comprehensively track CNAs in coding regions throughout neoadjuvant treatment. We observed our cases had a high CNA burden with a median of altered genome of 71.4% (range 1.1 -100%) (Fig 3A). Copy number (CN) ampli cations were seen frequently on chromosomes 1q, 3q, 5p, 7p, 8q, 19, 20, while losses or loss of heterozygosity (LOH) predominated on chromosomes 3p, 4q,5q, 9p and 17p, consistent with previous work 4, 18,22 (Fig. 3A).
Treatment-naïve tumours from REs and NRs displayed similar proportions of ampli cations, deletions, loss and LOH regions (Pearson's Chi-squared test >0.05) (Fig 3B) to post-treatment samples. Chemotherapy did not increase the overall percentage of altered genome, whereas RCTx, known to cause broad single or double strand breaks, led to a slightly higher proportion of altered genome in samples post-RCTx compared to pre-treatment samples, though the difference was not signi cant (t-test p-value=0.09) (Fig 3B). The radiation-associated changes were dominated by exome-wide small fragment losses or "scrambled" chromosomes ( Fig 3A) in a manner previously described in radiation-induced CNA and small fragment deletions 25,27 . We analysed the size of clonal, subclonal and private CN fragments in RE and NR and observed that genome fragments that change their CN state during treatment were signi cantly smaller than clonal fragments ( Fig 3C). Furthermore, we observed signi cantly smaller fragment sizes of private CNAs in NRs than in REs (t-test p-value < 0.001), which is consistent with the radiation-induced DNA damage in NRs.
We examined if dynamics in CNAs correlated with early treatment response in EAC. REs showed signi cant changes in copy number alteration burden with a mean of 63% (range 36-87%) of exome changing its CNA status (to any CN state) after the rst cycle of chemotherapy, whereas copy number states in NRs were more stable with a mean of 33% (range 11-61%) of exome altering CNA state (Wilcoxon rank-sum test p-value 0.009; Fig 3D). We hypothesize that NRs have genomes that confer pre-existing resistance to the effect of neoadjuvant chemotherapy, whereas REs have patterns of CNAs that make the cells more sensitive to chemotherapy, and so the rst dose of chemotherapy prompts further CNA evolution.
Genetic and transcriptomic dynamics in EAC drivers re ect treatment response We focused on genetic alterations in the 108 genes suggested to be drivers of EAC evolution (supplemental table 6) 17,18 . Genetic alterations in 15 common EAC driver genes across are shown in Fig 3E (genetic alterations in all 108 IntOgen drivers are shown in supplemental g 2 and 3). In 19 out of 27 patients, we observed that TP53 LOH or loss was associated with SNVs or indels of the same gene, presumably leading to a loss of normal function of TP53. CDKN2A showed CN loss or LOH in 20 out of 27 patients, which cooccurred with SNVs or indels in 5 patients, consistent with the nature of the disease 14 . We further observed ampli cations containing MYC, KRAS, ERBB2, ERBB3, GNAS, TOP2A and PI3KCA and deletion of SMAD4. Interestingly, CNAs in driver genes were mostly present before treatment and persisted in their CNA category (ampli cation, LOH, loss) during treatment, although the exact CNS could vary ( Fig   3E). Both RE and NRs showed no signi cant change in the burden of genomic alterations of driver genes during treatment (supplemental g 4A-F). REs showed a mean of 89.3 genetically altered putative driver genes at Timepoint A, 92.4 at B and 79.4 at Timepoint C (supplemental g 4B). NRs had a mean of 77.1 genetically altered putative drivers at Timepoint A, 65.4 at B and 91.0 at Timepoint C. Moreover, we did not detect any speci c cancer gene that was enriched for CNAs, SNVs or indels in the RE versus NR group (Fig 3E, suppl. g 2 and 3). We note that work in colorectal cancer has shown that mutations that are drivers are "on average" may not be a driver in a particular cancer type 10 and so we do not claim that all of these mutations are necessarily driving EAC evolution.
We next explored gene expression changes in high frequency drivers using our matched RNA-seq data. Signi cant expression changes were observed for most driver genes: expression tended to increase for GNAS, PIK3CA, SMAD4, decrease for ARID1A, CDKN2A, ERBB2, ERBB3, TOP2A, SMARCA4, KRAS, KMT2D and was stable for TP53, MYC, APC and AXIN1 (supplemental g 4G). We tested whether these expression changes were likely driven by underlying genetic copy number alteration, nding little evidence that copy numbers drove expression. Copy number state and RNA expression were strongly positively correlated for GNAS and moderately correlated for SMARCA4 and TP53 (supplemental g 7), and poor correlations were observed for ARID1A, APC, AXIN, CDKN2A, ERBB2, KMT2D, SMAD4, TOP2A, and no correlations were observed for ERBB3, KRAS, MYC and PIK3CA (supplemental g 5A). Correlation analyses between SNVs in driver genes and their expression, was only possible for TP53 and CDKN2A, as analyses for other drivers were underpowered. We observed signi cantly higher CDKN2A expression in CDKN2A mutated samples and signi cantly lower p53 expression in TP53 mutated samples (supplemental g 5B).
In summary, we observed major changes in gene expression of purported driver genes. However, these changes were rarely explainable by underlying genetic alterations or change.

Neoadjuvant treatment alters the EAC transcriptome
We assessed gene and pathway expression changes during neoadjuvant treatment using our matched RNA-seq data from 26 samples from 10 NRs and 53 samples from 17 REs. Principle component analysis (PCA) using hallmark pathways 28 separated samples from Timepoint A and B from samples at the end of the neoadjuvant treatment (Fig 4A), whereas samples from REs and NRs were intermixed. Hallmark pathways were grouped into "classes" according to their biological function ("oncogenic", "immune", "stromal", "cellular stress") 29 . PCA loading assessment showed that samples from Timepoint A and B were enriched for oncogenic pathways, whereas samples from Timepoint C were enriched in immune and stromal pathways and selected oncogenic pathways such as KRAS, Hedgehog and WNT signalling ( Fig 4B). We identi ed genes that were on average signi cantly differentially expressed hallmark pathways between pre/early-treatment samples and samples at the end of neoadjuvant treatment, and performed hierarchical clustering of all samples using these genes. Expectedly, and consistent with PCA analysis, the dendrogram separated into two groups ( Fig 4C): Group 1, containing most samples from Timepoint A/B, was enriched for "oncogenic" acting pathways, such as MYC, cell cycle associated pathways, DNA repair and MTORC signalling (Fig 4C), Group 2, including mostly samples from Timepoint C, showed signi cant enrichment of EMT and stemness associated WNT signalling and promoters of stem cell like state such as TGF-β signalling and hypoxia. Moreover, immune pathways such IL6-, IL2-signaling and in ammatory response and oncogenic KRAS-, Hedgehogand p53-signaling were signi cantly enriched in the post-treatment group (Fig 4C).
When examining pathway enrichment in early and late response in the different treatment regimens, neither REs nor NRs showed signi cant changes in pathways between Timepoint A and B (supplemental table 7 for REs, no KEGG enrichments found for NRs). However, REs showed some evidence of early upregulation in immune activation pathways (q-value=0.24) (supplemental table 7), whereas NRs did not. Of note, samples from REs and NRs showed enrichment of same pathways after completion of neoadjuvant treatment with upregulation of WNT-, PI3K, RAS-, MAPK-, JAK-STAT and Hedgehog signalling and immune related pathways (Fig 4D), indicating that both chemotherapy and RCTx lead to fundamental but similar transcriptomic changes in EAC, potentially as acute wounding/wound healing response.
We used immune deconvolution tools to delineate cellular composition changes from RNA-Seq data 30,31,32 . Using Cibersort 30 deconvolution, we observed a gradual increase in absolute numbers of immune cells during neoadjuvant treatment in REs and NRs (Fig 4E). Speci cally, we observed a doubling of immune cell in ltrates in samples after completion of neoadjuvant treatment compared to pre-treatment samples, which might be due to lower tumour cellularity in bulk RNA sequencing data from post-treatment samples. Proportions of individual immune cell types stayed stable during treatment ( Fig 4E, supplemental g 6A), in which CD4 and CD8 cells represented the most abundant immune cell fraction (40-50% of all immune cells). Other deconvolution tools (Consensus TME , Syllogist) 31,32 con rmed the gradual increase in immune in ltrates during treatment (supplemental g 6B, C) but there was variability amongst tools in the estimated proportions of immune cell types (supplemental g 6B).
In summary, gene expression programs, relevant for cancer biology and the surrounding immune microenvironment, signi cantly changed during treatment. Marked transcriptomic changes with enrichment of developmental programs such as EMT, stemness-associated WNT signalling and promoters of stem cell like state such as TGF-β signalling and hypoxia were observed, despite infrequent clonal replacement. Our data therefore suggest that cancer cells alter their phenotype without clonal replacement and that phenotypic plasticity might underly those phenotype alterations.
Immune escape correlates with treatment response Neoantigenic mutations can induce antitumor immune response, through the presentation of altered peptides that are being recognized as 'foreign' by the immune system 33,34 . With the recent introduction of combined immuno-chemotherapy for treatment of advanced EAC 35 , we focussed on understanding evolutionary dynamics of neoantigentic mutations, immune escape mechanisms and the response of the surrounding immune microenvironment.
Neoantigens were called from WES data on 24 samples from 10 NRs and 47 samples from 17 REs. We did not observe any signi cant changes in the neoantigenic mutation burden during neoadjuvant treatment ( Fig 5A). However, NRs showed a slight increase in neoantigenic mutation burden after being switched to RCTx ( Fig 5B) consistent with the observed increase in total mutational burden ( Fig 2B).
To understand the selection mechanisms of neoantigenic mutations, we examined if treatment enhances the immune response against the tumour cells. Depletion of neoantigenic mutation burden relative to the overall mutation burden is a well-established signature of immune-editing 36,37 . Our data showed a stable ratio between immunogenic and total mutations during neoadjuvant treatment (Fig 5C), indicating no enhanced negative selection of neoantigenic mutations during treatment. We further did not observe a signi cant difference in the subclonal neoantigenic mutational burden between preand posttreatment samples ( Fig 5D) and for most patients comparable amounts of neoantigenic SNVs on the clades of their phylogenetic tree (Fig 2D), which did not indicate removal of neoantigenic subclones during treatment.
Next, we explored the relationship between CNAs and neoantigens. We hypothesised that gains of alleles carrying antigenic mutations were likely to experience negative selection, and so expected that gained alleles would be relatively depleted for neoantigens compared to not-gained alleles. To test this, we calculated the copy-number normalised proportion of neoantigens in each CNA segment. Diploid regions showed a ratio of ~ 1 ( Fig 5E). Gains (balanced and unbalanced) were signi cantly depleted for neoantigens ( Fig 5E) while regions of CN loss were weakly, but not statistically signi cantly, enriched for neoantigens ( Fig 5E).
We hypothesized that rather than clonal evolution leading to removal of neoantigen clones, transcriptional downregulation of neoantigens during treatment might underlie immune evasion and potentially mediate treatment response. Counter to our hypothesis, we detected a signi cant increase in neoantigen expression during treatment (Fig 5F). High neoantigen expression correlated with high CD8 and CD4 cell in ltration (Fig 5G+H), suggesting that tumour cells with a high neoantigen expression were indeed recognized as 'non-self' and attracted immune cells. Consequently, we proposed that either EAC cells develop further immune escape mechanisms during treatment in order to avoid negative selection acting on their expressed neoantigens, or that the immune in ltrate becomes non-functional during treatment, and is not able to effectively eliminate tumour cells with high neoantigen expression. Both would be arguments for inclusion of immunotherapy in neoadjuvant treatment protocols.
We next examined the presence of immune escape mechanisms 38 (LOH or mutations in HLA, B2M mutations, or the interaction with the immune microenvironment by programmed death-ligand 1 (PDL-1) overexpression) through treatment. Genetic immune escape mechanisms such as HLA LOH or B2M mutations were detected signi cantly more often prior to treatment start (early), whereas transcriptomic PDL1 overexpression occurred mostly during and after neoadjuvant treatment (late) (Fig 5I, J, K) possibly as a consequence of expression in stromal/immune cells. Consequently, genetic disruption of antigen presentation machinery did not appear to experience positive selection during neoadjuvant treatment, whereas expression of PDL-1, potentially as a mechanism of rapid adaptation, may be positively selected during treatment.
We examined the correlation between immune escape (either HLA LOH, B2M mutations or transcriptomic PDL1 overexpression present) and treatment response. Samples from NRs were signi cantly more often immune escaped than REs ( Fig 5L). Further, patients with a poor pathological regression (Becker grade 3) were signi cantly more immune escaped (de ned as a minimum of one immune-escaped sample) than patients with an intermediate (Becker grade 2) or good regression (Becker grade 1) (Fig 5M), suggestive of a role for the immune system in potentiating the chemotherapy treatment response. As immune escape correlated with a poor early response to chemotherapy and a poor overall response to neoadjuvant treatment, our data suggest it has an important role in treatment resistance, and a possibility that immune therapy could help to circumvent resistance to standard treatment.
To understand phenotypic changes in the tumour and immune landscape in REs and NRs through treatment, we assessed changes in the abundance of different cell populations during treatment ( Fig 6E, supplemental g 6E). Dynamics in NRs during treatment were dominated by a decrease in tumour cell populations (C1-2, C4-6), decreases in the CD4 + cell population (C7) and an increase in one CD8 + cell population (C19) (Fig 6E). REs also showed decreasing abundances of tumour cell populations (C1, C4, C6, C12, C16, C17) and increases in a population of CD8 + cell (C19-cluster) (Fig 6E). To further explore qualitative differences of the immune cell compartment we selected all leukocytes (CD45 + cells) and reclustered using immune cell markers of cell activation, exhaustion, antigen-presentation and proliferation. Our immune cell clustering identi ed 24 distinct leukocyte subpopulations ( Fig 6F+G). We observed signi cant decrease in one activated CD8 cluster (IC15) in REs and NRs and signi cant decreases in two exhausted CD8 clusters (IC7 and IC18) in NRs, whereby NRs started with higher abundance of exhausted CD8 cells (IC18) in pretreatment samples than REs (Supplemental g 6F+G).
Given their relevance to immunotherapy, we explored the dynamics of the CD8 + and CD4 + cells in more detail. NRs showed stable CD8 + counts during treatment and a signi cant decrease in gated CD4 + cell counts, whereas the abundance of CD4 + and CD8 + cells stayed stable in REs ( Fig 6H). Correlation of IMC data and quantitative RNA-deconvolution of CD8 + and CD4 + cells showed a moderate correlation (R CD8 = 0.42 and R CD4 = 0.38, supplemental g 6H). As we noted a marked heterogeneity in the CD4 and CD8 abundance between ROIs in samples where multiple ROIs were analysed (supplemental g 6I), we assume that the moderate correlation between RNA-Seq and IMC data at least partially resulted from comparing IMC data from small 1-2mm 2 large ROIs to RNA-derived immune deconvolution from entire tumour sections.
To understand the degree of activation and effector function of the adaptive immune system, we tested CD4 + and CD8 + cells for expression of molecules informative about T-cell differentiation, regulation and immune activation. We observed a strong decrease in CD4 + and CD8 + cell fractions expressing CD45RO in NRs, indicative of a loss of non-naive T cells during RCTx (Fig 6I), whereas in REs the abundance of non-naïve CD45RO + CD8 + and CD45RO + CD4 + cells was stable ( Fig 6I). NRs showed a signi cant drop in CD8 + fractions expressing Granzyme B and in CD4 + fraction expressing Granzyme B or HLA-DR (Fig 6I), which indicates a decrease in the highly activated effector T cell phenotype in NRs. In contrast REs showed a decrease in Granzyme B and HLA-DR expression in CD4 + cells, whereas CD8 + cells stayed in their activated phenotype state during treatment (Fig 6I). The fraction of exhausted CD8 + cells and CD4 + cells (expressing PD1 or TIM-3) showed signi cant drops during RCTx in NRs (Fig 6I, supplemental g 7A). When calculating the ratio between CD8 and CD4 cell counts expressing cytotoxic Granzyme B + and expressing the exhaustion marker PD1 + , we observed a higher ratio in RE than in NR for both at all timepoints during treatment (Fig. 6J), indicative of a T-cell phenotype with more cytotoxic activity in REs than in NRs, which might explain the better treatment response.
Moreover, patients with poor remission (Becker grade 3) showed a signi cant drop in the CD8 + Granzyme B + /CD8 + PD1 + ratio during treatment, indicative of progressive CD8 cell exhaustion, which was not present in patients with good remission (supplemental g 7B).
In summary, our IMC data indicates a loss of CD4 + and CD8 + effectiveness in NRs during treatment via quantitative downregulation and decrease of activated CD4 + and CD8 + T cells. Compared to NRs, REs had no decrease in T-cell quantity and showed a more activated CD8 + and CD4 + phenotype at all timepoints during treatment. Our ndings strongly suggest that both the quantitative and qualitative composition of the T-cell compartment might represent a central mediator and predictor in treatment response of EAC. When considered together with our genetic and gene expression analyses, we hypothesize that immune escape and a less activated T-cell phenotype in NRs allows tumour cells to resist treatment and could potentially be exploited with new treatment strategies.

Dynamics of the T-cell Repertoire
Antigen-speci c T cell responses are a central feature of the antitumoral response. The TCR repertoire, measurable by TCR sequencing , provides a way to assess the magnitude of the T cell immune response. We therefore analysed the intra-tumoral TCR repertoire in a subset of 54 samples from 9 patients, to evaluate TCR dynamics during neoadjuvant treatment in EAC.
We observed an average of 296 total alpha TCR counts and 668 total beta TCR counts in treatment naïve NRs and an average of 229 total alpha TCR counts and 380 total beta TCR counts in treatment naïve REs, with a slight (non-signi cant) increase in total TCRs and unique TCRs at the end of neoadjuvant treatment ( Fig 7A, supplemental g 7C). The number of detected TCRs moderately correlated with the absolute score of T-cells from our Cibersort immune cell deconvolution from RNA-Sequencing data (R alpha =0.52, R beta = 0.48) (Fig 7B).
We examined changes in the proportions of TCR clone sizes in the repertoire during treatment exposure. We did not observe any differences in the overall shape of TCRs clone size distribution between different timepoints (Fig 7C). In order to track individual TCR clones in our longitudinally collected samples, we focussed on TCR clones that had signi cantly expanded between two timepoints. We observed early (starting between Timepoint A and B) and persisting TCR expansions in both REs with good pathological regression (RE11, RE23), while NR10 and RE5 with poor response to CTx showed late (starting between Timepoint B and C) and intermittent expansions (Fig 7D, supplemental g 7D). We then considered different thresholds of expansions, which were reactive to treatment and persisting through treatment (=expansions between timepoint A and C) in patients with data from all three timepoints. The trend was observed across the different expansion thresholds, with TCRs that expanded at least 4 or at least 8 times throughout treatment ( Fig 7E). Next, we tracked T-cell clone expansions during treatment in patients with different pathological regression grades. Patients with excellent pathological treatment response (Becker grade 1), showed signi cantly higher numbers of TCRs that expanded at least 4 or at least 8 times between two time points during treatment than patients with poor regression (Fig 7F). In summary our TCR data are suggestive that expansions of individual T-cell clones during treatment play a role in successful treatment response.
Together, we observed a stronger immune cell response of the adaptive immune system with a persisting activated CD8 phenotype and TCR expansions during treatment in patients with good treatment response, which is further supported by the previously observed reduced immune escape in those patients compared to patients with poor treatment response.

Discussion
Despite advances in treatment, the prognosis for patients with EAC remains poor, with frequent development of resistance to therapy. Genetic analyses of preand posttreatment samples from EAC patients have recently started to shed light on tumour clonal dynamics to understand treatment resistance from an evolutionary perspective 4,5,6 . In our prospectively-designed multi-timepoint and multi-omic study, based on longitudinal sampling with pre, during and post neoadjuvant treatment samples from the MEMORI trial, we integrated genetic analyses with transcriptomic analyses and analyses of the tumour immune microenvironment, in order to generate a holistic understanding of treatment resistance.
Somewhat surprisingly, we did not observe clonal replacement for particular subclones during neoadjuvant treatment. Our results are in line previous ndings in a multi-region WES study of pre-and posttreatment EAC samples, which similarly reported clonal persistence with rarely new occurring putative drivers posttreatment 5 . We note that another study by Findlay et al did report loss of mutations (and clones) through treatment 4 . Potentially the discrepancy could be explained by low cellularity in some of Findlay et al's post treatment samples.
Recent evidence implicates a key role of non-mutational resistance mechanisms (where cells can plastically switch phenotype without underlying genetic change) in drug resistance 8,9,40,41,42 . Here, we observed marked changes in transcriptional programmes without clonal replacement during neoadjuvant treatment, which may be evidence of tumour cell plasticity. Although, we emphasise that lots of mutations were gained and "lost" through the time course, so our data and analyses do not de nitively rule out that such phenotypic change is a consequence of genetic changes with underlying polyclonal evolution. Nevertheless, the enrichment of epithelial-mesenchymal transition (EMT) gene expression programmes in our post-treatment samples is consistent with previously reported plasticity of these programmes 41,8,43,44 . Further, EMT is known to be controlled by multiple signalling pathways, such as transforming growth factor-β (TGF-β), wingless/integrated (WNT), Notch, and Ras-mitogen-activated protein kinase (Ras-MAPK) pathways 41,45,46 , which all showed a signi cant increase during treatment in EAC. We speculate that inhibition of the pathways that promote EMT, with the goal of maintaining cancer cells in a state susceptible to chemotherapy, could represent a novel avenue for treatment.
Neoantigens and their generated antitumor immune response play a major role in tumour evolution and high neoantigen burden is associated with good prognosis in patients with a variety of cancers [47][48][49] and with successful immune-checkpoint therapy [50][51][52][53] . Indeed, recently approved combined immunochemotherapy for treatment of advanced EAC patients showed a signi cant improvement in overall survival compared to chemotherapy alone 54,55 . Here, our longitudinal analyses revealed relatively stable neoantigen burdens through treatment, pointing to a lack of stringent immunoediting and consistent with our detection of genetic disruption to antigen presentation machinery. Further, cellular analysis showed remodelling of the T-cell compartment during treatment, with expansion of T-cell clones in responders, and increasing T-cell exhaustion and PDL1 expression in non-responders. Previous studies have shown the immunomodulatory effects of C(R)Tx with an increase in PD-L1 from 45% in pre-NAT EAC to 77% post-CRT 56 58 . Following completion of neoadjuvant treatment, restaging was performed. In the absence of progression to metastatic/unresectable disease, patients underwent surgical resection. 10 NRs and 17 REs were enrolled for longitudinal molecular analyses. Most patients provided 3 cancer samples with one pre-treatment biopsy (Timepoint A), one biopsy collected after the rst cycle of chemotherapy (Timepoint B) and the surgical specimen sample after completion of neoadjuvant treatment (either CTx or RCTx) (Timepoint C) (Fig 1B and C). Two patients (RE10 and RE11) had multi-timepoint and multi-region samples (≥ 6 samples per patient). Further clinical and histopathological information is provided in Supplemental table 1. Tumour tissue was xed and stabilized with the PAX-gene xation method (PAXgene Tissue FIX Container (Qiagen)), which showed improved preservation biomolecules compared to formalin-xed tissue 59,60 . For eight samples PAXgene xated tissue was not available and formalin-xed tissue was used (supplemental table 8). Blood for germline control was collected in PAXgene Blood DNA Tubes (Qiagen). All patients gave informed consent for collection and molecular analysis of their sample material within the MEMORI trial protocol.

Section preparation:
PAXgene-xed para n-embedded (PFPE) or formalin-xed para n-embedded (FFPE) tissue was cut into consecutive sections. The rst section (4µm) was used for hematoxylin and eosin (H&E) stain. Eight consecutive sections were cut at 10µm under RNAse free conditions onto PEN Membrane Glass Slides (ThermoFisher) for subsequent isolation of DNA and RNA. Then, four consecutive sections were cut at 4µm onto standard slides for IMC analysis and the last section of 4µm was stained with H&E. Both H&E stainings were reviewed by a board-certi ed pathologist to assure high tissue quality and tumour cellularity. If tumour cellularity was estimated <50%, tumour areas were annotated for later tumour cell enrichment via laser capture dissecting microscopy. Pathological . Samples with a total DNA yield higher than 10ng and DNA peak detection at > 1000bp were taken forward for WES library preparation.
RNA extractions from PFPE sections were performed using the PAXgeneTissue RNA Kit from Qiagen according to the manufacturer's protocol. For FFPE sections RNA was extracted using the High Pure FFPET RNA Isolation Kit (Roche Molecular Systems, Basel, Switzerland). RNA quantity was measured using the Qubit uorometer RNA High Sensitivity assay (Life Technologies, Paisley, UK). Quality control was performed using the High Sensitivity RNA TapeStation system (Agilent, Santa Clara, California, USA). RNA quality was assessed using the DV200 value. RNA samples with a DV200 > 50% were eligible for further analysis.

Preparation of whole exome sequencing libraries:
Libraries were prepared using the SureSelectXT Low Input Target Enrichment System with Dual Indexing (Agilent Technologies, Santa Clara, US) according to manufacturer's instructions. A short mechanic fragmentation step using a Covaris S2 Sonicator to generate 150-200bp long fragments was performed and 8 to 14 PCR cycles were used for library enrichment depending on the DNA input. After puri cation, libraries were quanti ed by Qubit and run on the Agilent Tapestation using HSD1000 screentapes. Samples with su cient library DNA yield and characteristic fragment size distribution (~200-500bp) were further subjected to deep (~300x depth) WES.

Preparation of RNA sequencing libraries:
Libraries were prepared using the QuantSeq 3'mRNASeq Library Prep Kit FWD (Lexogen GmbH), according to the manufacturer's protocol with minor modi cations. To increase the total RNA library yields, the pre-PCR incubation time was increased to 1h. After puri cation, libraries were quanti ed by Qubit and run on the Agilent Tapestation using D1000 screentapes. Samples with su cient library cDNA yield and characteristic fragment size distribution (~170-700bp) were further subjected to RNA-Sequencing.

WES and RNA Sequencing:
Sequencing libraries were multiplexed and sequenced on an Illumina Novaseq, typically using S2 ow cells. Read length and depth was varied as required by library composition. Target depth for WES was 250x and >2 million reads per sample for RNA-Seq. Sequencing was performed by the Genome Centre at Queen Mary University London.

Copy number analysis:
Copy number analysis and cellularity estimates were performed with Sequenza 69 . Pre-processing of bam les to seqz les was done using the sequenza-utils bam2seqz algorithm and binned using the sequenza-utils seqz_binning -w50 command. The per patient set of break points, binned depth-ratio and B-allele frequencies data from binned seqz les were then inputted into the sequenza algorithm (version 2.1.2) to determine allele speci c copy-numbers, ploidy Ψ and purity ρ estimates 69 . The threshold segment length for tting cellularity and ploidy parameters was set to 10 6 bp to lter out short segments, that can cause noise for cellularity and ploidy estimates. The initial parameter space searched was restricted to [ρ | 0.1 ≤ ρ ≤1] and [Ψ | 1 ≤ Ψ ≤ 5]. Upon manual review of the results, we identi ed several samples with unreasonable ts (cases where calls suggested extremely variable ploidy values across samples from the same patient). For these samples, we manually identi ed alternative solutions consistent with the other samples and somatic variant calls. Clonal CNAs were de ned as CNAs present in all samples of a patient, subclonal CNAs were de ned as CNAs present in > 1/3 of samples, but not all samples of a patient and private CNAs as CNAs present in ≤ 1/3 of samples of a patient. For patients with only 2 available samples subclonal/private CNAs were de ned as CNAs present in 1 sample. To calculate the fraction of exome, that changes its CNS between timepoint A and B (Fig 3D), 15 patients with available samples from timepoint A and B and no multi-region samples at timepoint A or B were included.

SNV detection:
Somatic mutations were rst called for each tumour sample separately against matched blood samples with Mutect2 (version 4.1.2.0) 70 . Variants detected were annotated using ANNOVAR 71 . Variants detected in any tumour sample (marked PASS, coverage AD 5 in both normal and tumour, at least 3 variant reads in the tumour and a maximum variant read of 0 in the normal reference genotype) were merged into a single list of "candidate mutations". For later MOBSTERanalysis 16 for multi-region/multi-timepoint samples (see below), SNV caller bam-readcount 72 was used to call nucleobases and coverage at each candidate mutation position in all samples of a respective patient. Candidate driver genes of EAC were annotated using 108 cancer driver genes for oesophageal and gastric cancer reported by IntOGen © 73 .

Identi cation of subclonal and clonal mutations per sample:
For each sample "candidate mutations" were classi ed as clonal or subclonal mutation using MOBSTER 74 . Based on the calculated variant allele frequency (VAF), i.e., the ratio of read counts mapping to the mutant allele, over the total coverage at the variant locus, which was previously normalized by the tumour purity and tumour copy number segments, MOBSTER clusters input mutations as clonal peak mutations (=clonal mutations) or tail mutations for neutral mutations (= subclonal mutations).

Construction of phylogenetic trees:
For phylogenetic tree analysis we analysed SNVs from multi-region/multi-timepoint samples of each patient using MOBSTER for multiple samples. To avoid previously described biases in clonality analysis of multi-region/multi-timepoint analysis 16 , MOBSTER excludes mutations belonging to the evolutionary neutral-tail for clonality analysis in multi-region/multi-timepoint samples. From non-tail mutations maximum-parsimony trees were reconstructed with the Parsimony Ratchet method 75 using the phangorn package for R 76 . Included non-tail mutations within a sample were considered to be mutated (state 1) and others considered to be non-mutated (state 0) in a given sample. The Ratchet was run for a minimum of 10 and a maximum of 10 3 iterations and terminated after 10 rounds without improvement. The acctran algorithm [76][77][78] was utilized to approximate ancestral character states. From these a set of shared mutations between samples and mutations that were uniquely mutated on each clade of the phylogeny were obtained.

Extraction of mutational signatures and dN/dS analyses:
The extraction of mutational signatures in tumour samples was carried out with SigPro ler 79,80 , using COSMIC signatures as reference signatures and the default settings. Un ltered signatures for each tumour sample are shown in the extended data (supplemental g 1C). In order to reduce the number of called mutation signatures to more prevalent ones, we rerun the analysis excluding signatures that had a COSMIC signature weight of <5% in any of the responsiveness groups at a given timepoint. Next, we evaluated mutational footprints of chemo-and radiochemotherapy in our cohort. COSMIC signatures were called separately for treatment-naïve mutations, mutations occurring under chemotherapy (= "CTx-SNVs") and mutations occurring under RCTx (="RCTx-SNVs "). Treatment-naive mutations included mutations at timepoint A in REs and NRs, CTx-SNVs included newly occurring mutations at timepoint B in REs and NRs and newly occurring mutations at timepoint C in REs and RCTx-SNVs harboured newly occurring mutations at timepoint C in NRs. dN/dS ratios were calculated for treatment-naïve mutations, CTx-SNVs and RCTx-SNVs using the dNdScv R package 81 .

Genetic alterations and RNA expression in EAC drivers:
SNVs, indels and CNA were assessed in 108 IntOGen © drivers 17 of oesophageal and gastric cancer. To assess the number of altered driver genes in each sample, driver genes were counted as altered if they harboured one or more SNVs, indels or CNAs. To assess the correlation between CNS and RNA expression of 15 high frequency driver genes, samples with matching WES and RNA-Sequencing data were analysed (n=66). RNA-Sequencing data was pre-processed as described in 4.1 and 4.2. Changes in expression levels of driver genes between treatment-naïve samples and post-treatment samples were calculated by Wilcoxon tests. For correlation analyses between the expression level and the CNS of driver genes, the expression level was adjusted for the estimated tumour cellularity. Correlation between expression levels and CNS were calculated for each driver genes by Pearson correlation. For expression analyses of driver genes in samples with mutated and non-mutated driver status, only CDKN2A and TP53 displayed a su ciently powered analysis. The expression level of the respective driver was adjusted for the estimated tumour cellularity. Expression levels between mutated and non-mutated samples were compared by the Wilcoxon test.

HLA haplotyping:
HLA haplotyping for each patient was carried out with the raw paired fastqs of blood-derived normal samples using Optitype 82 . The pair of A, B and C alleles with the highest score were taken as a patient's haplotype.

Calling immune escape:
We predicted four types of immune escape mechanisms: (1) somatic mutations in one of the HLA alleles or (2) in the B2M gene 83 ; (3) loss of an HLA haplotype through LOH 84 ; and (4) PD-L1 overexpression 85 .
Somatic mutations in the HLA locus were called using polysolver 83 . First, alleles were converted into a polysolver-compatible format (lower case, digits separated by underscore) and outputted into a patient-speci c winners.hla.txt le, following the standard output of polysolver. Next, the mutation detection script of polysolver (shell_call_hla_mutations_from_type) was run on matched tumour-normal pairs to call tumour-speci c alterations in HLA-aligned sequencing reads using MuTect (v1.16). In addition, Strelka2 (v2.9.9) 86 was run independently to detect short insertions and deletions in HLA-aligned reads as this software offers increased sensitivity over polysolver's default caller. Finally, both single nucleotide mutations and insertions/deletions passing quality control were annotated by polysolver's built-in annotation script, shell_annotate_hla_mutations. LOH at the HLA locus was assessed using the software LOHHLA 84 and sequenza 69 . First, the allele-speci c copy number as predicted by sequenza at the HLA-A, -B, -C loci was evaluated. Samples with a predicted minor allele copy number of 0 (e.g. 2:0, 3:0) were labelled as candidate LOH. Then, we ran LOHHLA with the previously generated winners.hla.txt as input. A type-I allele of a patient was annotated as "allelic imbalance" (AI) if the p-value testing the difference in evidence for the two alleles was lower than 0.01.
Alleles with AI were labelled as LOH if the following criteria held: (i) the predicted copy number of the lost allele was below 0.5 with con dence interval strictly below 0.7; (ii) the copy number of the kept allele was above 0.75; (iii) the number of mismatched sites between alleles was above 10. All LOHHLA-based LOH called were in candidate LOH samples, but some candidate LOH could not be con rmed most probably due to low purity and read density.
Mutations in B2M were called if the sample contained contained a nonsynonymous somatic, frameshift (FS), stop-loss or stop-gain mutation located inside the exonic regions of the exons of the B2M gene, as called by ANNOVAR 71 . PD-L1 overexpression was assessed using matched RNA-seq data from the samples. RNA-Seq raw data was pre-processed and then normalized as described in 4.1 and 4.2. The mean expression value of normalized gene counts (in transcripts per million mapped reads (TPM)) of PD-L1) in treatment-naïve EAC samples was calculated and considered as PDL-1 treatment-naïve baseline. PDL-1 overexpression was de ned as ≥ PD-L1 treatment-naïve baseline + 1 standard deviation.

Calling neoantigens:
Neoantigens were predicted from variant call tables and HLA types using the NeoPredPipe 87 , python-based pipeline combining Annovar and netMHCpan 4.0.
Brie y, all quality-controlled somatic mutations were annotated using Annovar 71 , and for all non-synonymous exonic mutations the mutated peptide sequence was predicted. We took any 9-and 10-mer spanning the mutated amino acid(s), resulting in either (i) a 19-aa window for SNVs or (ii) a peptide until the next predicted stop codon for FSs. These peptides were evaluated according to their novelty and predicted binding strength to the patient's six-allele HLA set. Peptides novel compared to the healthy human proteome with binding rank 2 and below (amongst the best 2% of binders compared to a large set of random peptides) were reported as neoantigens. All patient-speci c HLA alleles were used for neoantigen predictions, regardless of mutation or LOH status of the HLA locus. We considered a mutation neoantigenic if at least one of its downstream mutated peptides were a neoantigen with respect to any of the patient's six HLA alleles. The clonality of neoantigentic mutations was assessed based on MOBSTER 74 analyses (more details in 3.4). To calculate the enrichment/depletion of neoantigens in regions of different copy number states, we calculated for each sample, the ratio between the fraction of copy-number adjusted NeoSNVs present in a certain CN category (loss, cnLOH, high LOH, diploid, gain and ampli cation) and the fraction of genome showing the respective CN category.

Expression of neoantigens:
To assess the expression levels of neoantigenic peptides, we analyzed gene expression levels of genes harbouring neoantigenic mutations. RNA-Seq raw data was pre-processed and then normalized as described in 4.1 and 4.2. As a proxy for neoantigen expression, expression levels (CPM) of protein-coding genes harbouring neoantigenic mutations were assessed. Mean expression levels of neoantigens were calculated for each sample and compared between samples from different timepoints via Wilcoxon tests. For analyses of neoantigen expression levels in low and high CD4 and CD8 cell groups, CD4 and CD8 Tcell levels were deconvoluted from normalized RNA-sequencing data using CIBERSORT 30 (more details in methods part 4.5). Samples were assigned to high and low CD4 and CD8 cell groups based on the median CD4 and CD8 score of the study sample.
4. Bioinformatic analysis of RNA-Sequencing 4.1 RNA sequencing -alignment and ltering of RNA samples: After initial quality control with FastQC 88 and default adaptor trimming with Skewer 62 , paired-end reads were aligned to GRCh38 reference genome and version 28 of the Gencode GTF annotation using the STAR 2-pass method 89 . Read groups were added with Picard v.2.5.0 90 . Per gene read counts were produced with htseq-count that is incorporated into the STAR pipeline 89 . Raw gene counts were rst ltered for reads uniquely assigned to non-ribosomal protein-coding genes located on canonical chromosomes (chr1-22, X and Y). If samples had less than 1.5M of these 'usable' reads they were re-sequenced to improve coverage. Where possible, the same library preparation pool was sent again for sequencing. These top-ups' of the second library where necessary) contained <1.5x10 6 reads, samples were excluded for further downstream analysis.

Gene expression normalisation and ltering of expressed genes:
A master table including raw read counts for each samples of 58,243 genes were converted to a DGEList object via DESeq2 91 . Next, a list of expressed genes (n=11,900) was determined by edgeR 92 default ltering, to determine genes with su ciently large counts to be retained in a statistical analysis. Normalization

Immune deconvolution from RNA-Sequencing data
Immune cells were deconvoluted with CIBERSORT 30 a computational method for quantifying cell fractions from bulk tissue gene expression pro les based on support vector regression with prior knowledge of expression pro les from puri ed leukocyte subsets. For each tumour sample, the CPM expression of genes converted to Gene Symbols was used as input for CIBERSORT single sample immune deconvolution and was uploaded to the online available CIBERSORT tool (https://cibersort.stanford.edu/index.php). A LM22 signature le compromising 22 immune cell types and included in the deconvolution tool was used for deconvolution. Deconvoluted immune cell data from our samples was outputted as absolute and relative scores. For further analysis CIBERSORT output tables were processed in R. For each sample CIBERSORT scores of "Macrophages.M0", "Macrophages.M1" and "Macrophages.M2", were summed to a joint score for "Macrophages", scores of activated and resting cells of the same cell types were summed (e.g. "NK.cells.resting" and "NK.cells.activated" were summed to a "NK cell" score), scores for "B.cells.naive", "B.cells.memory" and "Plasma.cells" were summed to "Other B cells", scores for "T.cells.follicular.helper", "T.cells.regulatory..Tregs.", "T.cells.gamma.delta" were joint to "Other T cells" and scores for "T.cells.CD4.naive", "T.cells.CD4.memory.activated" and "T.cells.CD4.memory resting" were added to common score for "CD4 T-cells". For validation purposes further immune cell deconvolution methods were used. For each tumour sample, the CPM expression of protein-coding genes and expressed genes converted to entrez gene IDs (n= 10,463) was used as input for single sample gene set enrichment analysis using the Consensus TME gene sets 31 and the GSVA R package 93 in R. As a third immune deconvolution tool Syllogist 32 was used. For each tumour sample, the CPM expression of protein-coding genes converted to HUGO nomenclature (n= 39,921) was used as input for the provided Syllogist script 32 , which was run in R. Odds ratios of immune cells of interests were visualized in bar graphs generated in R.
IMC staining: IMC stainings were performed on 54 ROIs from PAXgene xed samples and 12 ROIs from FFPE samples (supplemental table 8). Incubation, rehydration and epitope retrieval was performed on 4μm thick tissue sections as described in previous work 95 . The sections were then stained with a mix of metal-labeled primary antibodies (Supplemental table 4) diluted in TBS with 0.5% BSA and incubated at 4°C overnight. Sections were then rinsed, stained with Iridium Cell-ID and air-dried as previously described 95 . Slides were stored at room temperature until image acquisition.
IMC image acquisition A 1-2 mm2 image per sample was acquired using a Hyperion Imaging System (Fluidigm). Tuning of the instrument was performed according to the manufacturer's instructions. Regions of interest were determined based on a pathologist's annotations on consecutive HE-stained section of each sample.
Tissue sections were laser-ablated spot-by-spot at 200 Hz resulting in a pixel 2 size/resolution of 1 mm. Preprocessing of the raw data was conducted using the CyTOF software v7.0. Image acquisition control was performed using MCD Viewer v1.0.560.6. For larger samples multiple ROIs were obtained.
IMC raw data processing: Raw data were exported as MCD les. MCD les were converted into tiff les and segmented into single cells using an unbiased, supervised analysis pipeline adapted from https://github.com/BodenmillerGroup/ImcSegmentationPipeline. In more detail, nucleated cells were segmented using Ilastik v1.3.3 96 and CellPro ler v4.1.3 97 . Using Ilastik, pixels were classi ed into nuclei, cytoplasm and background. A probability map for the three classi cations was generated and exported to create a cell mask with CellPro ler. Data folders containing tiff images of the 21 markers and the combined cell mask were loaded into histoCAT v1. 76. In histocat mean marker intensity of pixels and spatial features of segmented cells were analyzed and exported as a table containing single cell data for each marker. Subsequent processing of single-cell level data was performed to OMIQ. TCR sequencing: The α and β chains of the TCR repertoire of 27 samples from 10 patients were sequenced based on total RNA isolated from tumour samples. Brie y, the method introduces unique molecular identi ers attached to individual cDNA molecules to provide a quantitative and reproducible method of library preparation. Full details for both the experimental TCRseq library preparation and the subsequent computational analysis (V, J, and CDR3 annotation) using Decombinator are published in our previous work 98,99 TCR statistical analysis: Statistical analyses were performed in R. Wilcoxon rank tests were used to compare the abundance of detected TCRs between timepoints in NRs and RE. For correlation analysis between deconvoluted T-cells from Consensus TME and total TCRs, Consensus TME scores for 'CD4 T-cells', 'CD8 T-cells' and 'regulatory Tcells' were summed to a 'Total T-cell score' for each sample. Correlation between calculated 'Total T-cell scores' and abundance of detected