The shared frameshift mutation landscape of microsatellite-unstable cancers suggests immunoediting during tumor evolution

This is a self-archived version of an original article. This version may differ from the original in pagination and typographic details. Author(s): Title: Year: Version: Copyright: Rights: Rights url: Please cite the original version: CC BY 4.0 https://creativecommons.org/licenses/by/4.0/ The shared frameshift mutation landscape of microsatellite-unstable cancers suggests immunoediting during tumor evolution © The Author(s) 2020 Published version Ballhausen, Alexej; Przybilla, Moritz Jakob; Jendrusch, Michael; Haupt, Saskia; Pfaffendorf, Elisabeth; Seidler, Florian; Witt, Johannes; Hernandez, Sanchez Alejandro; Urban, Katharina; Draxlbauer, Markus; Krausert, Sonja; Ahadova, Aysel; Kalteis, Martin Simon; Pfuderer, Pauline L.; Heid, Daniel; Stichel, Damian; Gebert, Johannes; Bonsack, Maria; Schott, Sarah; Bläker, Hendrik; Seppälä, Toni; Mecklin, Jukka-Pekka; Ten, Broeke Sanne; Nielsen, Maartje; Heuveline, Vincent; Krzykalla, Julia; Benner, Axel; Riemer, Angelika Beate; von Knebel Doeberitz, Magnus; Kloor, Matthias Ballhausen, A., Przybilla, M. J., Jendrusch, M., Haupt, S., Pfaffendorf, E., Seidler, F., Witt, J., Hernandez, S. A., Urban, K., Draxlbauer, M., Krausert, S., Ahadova, A., Kalteis, M. S., Pfuderer, P. L., Heid, D., Stichel, D., Gebert, J., Bonsack, M., Schott, S., . . . Kloor, M. (2020). The shared frameshift mutation landscape of microsatellite-unstable cancers suggests immunoediting during tumor evolution. Nature Communications, 11, Article 4740. https://doi.org/10.1038/s41467-020-18514-5 2020 ARTICLE The shared frameshift mutation landscape of microsatellite-unstable cancers suggests immunoediting during tumor evolution Alexej Ballhausen 1,2,3,16, Moritz Jakob Przybilla 1,2,3,16, Michael Jendrusch1,2,3,16, Saskia Haupt 4, Elisabeth Pfaffendorf1,2,3, Florian Seidler 1,2,3, Johannes Witt1,2,3, Alejandro Hernandez Sanchez 1,2,3, Katharina Urban1,2,3, Markus Draxlbauer1,2,3, Sonja Krausert1,2,3, Aysel Ahadova1,2,3, Martin Simon Kalteis 1,2,3, Pauline L. Pfuderer 1,2,3, Daniel Heid1,2,3, Damian Stichel1,5, Johannes Gebert1,2,3, Maria Bonsack 6,7,8, Sarah Schott 9, Hendrik Bläker10, Toni Seppälä 11, Jukka-Pekka Mecklin12,13, Sanne Ten Broeke14, Maartje Nielsen14, Vincent Heuveline 4, Julia Krzykalla15, Axel Benner15, Angelika Beate Riemer 6,7, Magnus von Knebel Doeberitz 1,2,3 & Matthias Kloor 1,2,3✉ The immune system can recognize and attack cancer cells, especially those with a high load of mutation-induced neoantigens. Such neoantigens are abundant in DNA mismatch repair (MMR)-deficient, microsatellite-unstable (MSI) cancers. MMR deficiency leads to insertion/ deletion (indel) mutations at coding microsatellites (cMS) and to neoantigen-inducing translational frameshifts. Here, we develop a tool to quantify frameshift mutations in MSI colorectal and endometrial cancer. Our results show that frameshift mutation frequency is negatively correlated to the predicted immunogenicity of the resulting peptides, suggesting counterselection of cell clones with highly immunogenic frameshift peptides. This correlation is absent in tumors with Beta-2-microglobulin mutations, and HLA-A*02:01 status is related to cMS mutation patterns. Importantly, certain outlier mutations are common in MSI cancers despite being related to frameshift peptides with functionally confirmed immunogenicity, suggesting a possible driver role during MSI tumor evolution. Neoantigens resulting from shared mutations represent promising vaccine candidates for prevention of MSI cancers. https://doi.org/10.1038/s41467-020-18514-5 OPEN 1Department of Applied Tumor Biology, Institute of Pathology, University of Heidelberg, Heidelberg, Germany. 2Collaboration Unit Applied Tumor Biology, German Cancer Research Center (DKFZ), Heidelberg, Germany. 3Molecular Medicine Partnership Unit (MMPU), Heidelberg University Hospital and EMBL, Heidelberg, Germany. 4Engineering Mathematics and Computing Lab (EMCL), Interdisciplinary Center for Scientific Computing (IWR), Heidelberg University, Heidelberg, Germany. 5Clinical Cooperation Unit Neuropathology, German Cancer Research Center (DKFZ), Heidelberg, Germany. 6Immunotherapy and Immunoprevention, German Cancer Research Center (DKFZ), Heidelberg, Germany. 7Molecular Vaccine Design, German Center for Infection Research (DZIF), partner site Heidelberg, Heidelberg, Germany. 8Faculty of Biosciences, Heidelberg University, Heidelberg, Germany. 9Department of Obstetrics and Gynecology, University Hospital Heidelberg, Heidelberg, Germany. 10 Institute of Pathology, University Hospital Leipzig, Leipzig, Germany. 11 Department of Gastrointestinal Surgery, Helsinki University Hospital and University of Helsinki, Helsinki, Finland. 12 Department of Education and Research, Central Finland Central Hospital, Jyväskylä, Finland. 13 Department of Sports and Health Sciences, University of Jyväskylä, Jyväskylä, Finland. 14 Department of Clinical Genetics, Leiden University Medical Center, Leiden, The Netherlands. 15 Division of Biostatistics, German Cancer Research Center (DKFZ), Heidelberg, Germany. 16 These authors contributed equally: Alexej Ballhausen, Moritz Jakob Przybilla, Michael Jendrusch. ✉email: [email protected] NATURE COMMUNICATIONS | (2020) 11:4740 | https://doi.org/10.1038/s41467-020-18514-5 | www.nature.com/naturecom munications 1 1234567890():,; DNA mismatch repair (MMR) deficiency is a major mechanism causing genomic instability in human cancer. MMR-deficient cancers accumulate an exceptionally high number of somatic mutations. These mutations encompass certain types of single-nucleotide alterations but mostly insertion/ deletion (indel) mutations at repetitive sequence stretches termed microsatellites (microsatellite instability (MSI))1,2. About 15% of colorectal cancers (CRCs), up to 30% of endometrial cancers (ECs) and multiple other tumors display the MSI phenotype3. MSI tumors can develop sporadically or in the context of Lynch syndrome, the most common inherited cancer predisposition syndrome. Due to this very specific process of genomic instability, the pathogenesis of MSI cancers can be precisely dissected4: indel mutations of coding microsatellites (cMSs), which almost exclusively affect coding mononucleotide repeats5–7, in genomic regions encoding tumor-suppressor genes are considered major drivers of MSI tumorigenesis. Importantly, the same indels that inactivate tumor-suppressor genes simultaneously cause translational frameshifts, thereby generating unique frameshift peptides as a major source of neoantigens4,8. For the recognition of neoantigens by the immune system, processing through the cellular antigen processing machinery and presentation by human leukocyte antigen (HLA) class I molecules on the tumor cell surface are essential prerequisites. These HLA class I molecules consist of a heavy chain and a non-covalently bound light chain (encoded by the Beta-2-microglobulin [B2M] gene), both of which are essential for functional antigen presentation. The likelihood of HLA binding for a defined peptide depends on the HLA genotype, as every individual harbors six alleles (HLA-A,HLA-B,HLA-C, two alleles each) that encode for HLA class I heavy chains9. The specific mutational steps required for malignant transformation during the evolution of MSI tumors are thus also responsible for their pronounced immunogenicity. MSI tumors are commonly associated with dense lymphocyte infiltration and pronounced local responses of the adaptive immune system at the tumor site10–14. Immune recognition of MSI tumor cells is not only responsible for a comparatively favorable clinical course but also reflected by the high sensitivity of advanced-stage MSI cancers toward immune checkpoint blockade (ICB)15,16. However, some patients do not respond to ICB treatment. We and others previously showed that specific T cell responses against a few MMR deficiency-induced frameshift neoantigens occur prior to and after ICB15,17,18. However, the landscape of frameshift peptides and potential epitopes in MSI cancer has not been described systematically. One important reason for this gap of knowledge is the fact that short-read next-generation sequencing (NGS) approaches have a limited sensitivity for the detection of indel mutations at homopolymer sequences such as frameshift peptide-inducing cMS19–21. Here we map the frameshift peptide landscape of the two most common MMR-deficient cancer types, CRC and EC, using a tool for the quantification of cMS mutation patterns (ReFrame, REgression-based FRAMEshift quantification algorithm) combined with NetMHCpan 4.0, a state-of-the-art in silico epitope prediction tool9. We reveal and functionally validate a set of previously unknown shared frameshift neoantigens in MSI cancers. Our results provide evidence for continuous immunoediting during MSI tumor evolution and underline the potential of neoantigen-based vaccines against MSI cancers. Results cMS mutation frequencies in MSI CRC and EC. Short-read NGS approaches are not ideally suited for mutational and frameshift peptide profiling of MSI cancers19–21, showing a high variability in mutation frequency regarding the detection of mutations in different cMS candidates like those located in the genes TGFBR2,SLC35F5,orTFAM (Supplementary Table 1). Importantly, cMS repeats of increased length, which are most susceptible to mutations and therefore encompass the most important mutational targets during MMR-deficient tumorigenesis, are missed by NGS technology that is in common use today3,5,19,22–26. To fill this gap and precisely quantify cMS mutation patterns and their resulting translational reading frames in MMR-deficient cancers, we developed an algorithm based on fragment length analysis as the current gold standard for the detection of MSI. ReFrame, our REgression-based FRAMEshift quantification algorithm, allows unbiased quantitative detection of indel mutations by solving a linear system of mathematical equations to remove stutter band artifacts, which result from polymerase slippage events during PCR amplification and subsequent nucleotide gains and losses similar to MSI-induced indels (Supplementary Fig. 1). We used ReFrame in a series of MSI CRCs (n=139; Supplementary Table 2) to screen for mutations in 41 cMS residing in 40 target genes derived from the first comprehensive cMS database (SelTarbase, Version 201307)27. In addition, we investigated mutation profiles in a cohort of MSI ECs (n=28). In agreement with previous reports23,27, our results show that the load of indels at cMS in MSI CRC and EC is high and that multiple concomitant indels at several cMS in the same tumor are very common. Although most CRCs and ECs were distinguishable based on the cMS mutation patterns, a large set of cMS mutations were shared by the majority of MSI CRCs and/or MSI ECs (Fig. 1and Supplementary Fig. 2). Moreover, we observed a significant variation of the number of mutations per tumor, ranging from 8 to 29 (median: 20) out of 41 analyzed cMS in MSI CRC and from 8 to 25 in MSI EC (median: 18). The observed variation suggests potential differences in the frameshift peptide load of MSI tumors. Potential clinical consequences, e.g., for the sensitivity toward ICB, should be assessed in future clinical studies28–30. ReFrame is not only able to quantify mutation frequency but also to distinguish indel mutation types, which is crucial for the prediction of the frame of the resulting frameshift peptides. As the translation of nucleotide into amino acid sequences is based on three base codons, every mutation in a homopolymer region can either result in a simple deletion or insertion of amino acids or in two entirely different frameshift peptide reading frames: deletions of one nucleotide (further referred to as minus 1 (m1)) or insertions of two nucleotides (plus 2 (p2)) will result in a shift to a frame here referred to as minus one (M1), while deletions of two nucleotides (minus 2 (m2)) or insertions of one nucleotide (plus 1 (p1)) will result in a shift to a frame referred to as minus two (M2) (Fig. 2, Supplementary Fig. 3, and Supplementary Table 3). The results demonstrate that m1 mutations, resulting in M1 reading frames, were the predominant mutation type (77% in MSI CRC, Fig. 2). The M1/M2 distribution varied significantly across distinct cMS, with significantly elevated numbers of M2 mutations in BANP,TAF1B, and ELAVL3, whereas in ACVR2A, HPS1,SLC35F5, and TCF7L2 there were significantly more mutations leading to an M1 frameshift than expected by chance (Bonferroni corrected binomial test, p< 0.05; Supplementary Table 4). Landscape of frameshift peptides and predicted major histocompatibility complex (MHC) ligands. Following the detection of shared indel mutations in MSI CRC and EC, we evaluated the ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-18514-5 2NATURE COMMUNICATIONS | (2020) 11:4740 | https://doi.org/10.1038/s41467-020-18514-5 | www.nature.com/naturecommunications immunogenic potential of the frameshift peptides and predicted MHC ligands resulting from antigen processing. We used NetMHCpan 4.0, a state-of-the-art MHC ligand prediction tool based on artificial neural networks, to predict neopeptides that are possibly presented as epitopes by HLA class I antigens encoded by the most important HLA supertypes9,31,32. Applying commonly accepted IC 50 thresholds, we distinguished between three classes of peptides with high (IC 50 < 50 nM), low (50 nM < IC 50 < 500 nM), and very low (500 nM < IC 50 < 5000 nM) predicted HLA-binding affinity31,33.Asafirst step, we analyzed all possible frameshift peptide sequences derived from the M1 and M2 frameshifts of the 41 cMS. We then complemented this set to cover all possible frameshift peptides (n=524) derived from 264 cMS with a length of ≥8 nucleotides published in SelTarbase (Supplementary Data 1)27. Our results indicate multiple frameshift peptides resulting from M1 or M2 frameshift mutations that are potentially recognized by the immune system. We detected a wide range of variability with regard to the number of predicted putative epitopes maximally contained within a defined frameshift peptide. The highest number of predicted putative high-affinity epitopes within a frameshift peptide was 23 (for the M1 frame of P4HB), (low affinity: 92 predicted putative epitopes in M1 SPINK5; very low affinity: 375 predicted putative epitopes in M1 P4HB). Other cMS mutation-induced frameshift peptides showed a complete lack of the predicted epitopes (Fig. 3and Supplementary Data 2). For HLA-A*02:01, the most common HLA allele in the USA European Caucasian population34, one or more high-affinity peptides were predicted for 19.8% of the frameshift peptides. HLA-A*02:01 epitopes with lower affinity were present in 39.5% (≤500 nM) and 59.7% (≤5000 nM) of candidates (Supplementary Fig. 4 and Supplementary Table 5. Plots depicting predicted HLA-binding peptides for all frameshift peptides and HLA alleles are available in the Source Data folder.). To make the potential impact of certain cMS candidates more tangible and to identify frameshift peptides with potentially highest relevance for immune recognition, we defined a general epitope likelihood score (GELS; see “Methods”section “Computation of immunological scores”). GELS accounts for MHC ligand prediction and the prevalence of the respective HLA allele in a defined population, as the latter influences the probability of a frameshift peptide to encompass an MHC ligand recognized by CRC and EC dendrogram Tumor type (CRC, EC) Tumor dissimilarity a bc Tumor identifiers (CRC) Tumor identifiers (EC) Percentage mutated [%] Coding microsatellites PCA dimension 1 PCA dimension 2 EC CRC CRC and EC PCA Candidate gene mutation rate by tumor sample 100 80 60 40 20 0 Fig. 1 Mutation frequencies of cMS in MSI CRC and EC. a The relative frequency of mutant alleles is shown for 41 cMS (rows) in CRC and EC tumor samples (columns) derived using ReFrame. cMS were sorted top–down according to their mutation frequency indicated by blue boxes of different intensity. Dark blue represents high mutation frequency, whereas pale blue represents low mutation frequency and white absence of mutations. Black boxes indicate missing data points. The respective tumor samples are shown below each column. cMS were analyzed for both CRC and EC and are depicted separately for each tumor type (left panel: CRC, right panel: EC). The complete dataset is provided in Supplementary Fig. 2. Source data for Supplementary Fig. 2 are provided as a Source Data file. bDendrogram of CRC and EC samples with respect to their mutation patterns where the color bar below indicates the tumor type (CRC, blue; EC, red). cPrincipal component analysis of EC (red) and CRC (blue) samples with respect to their mutation patterns. NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-18514-5 ARTICLE NATURE COMMUNICATIONS | (2020) 11:4740 | https://doi.org/10.1038/s41467-020-18514-5 | www.nature.com/naturecom munications 3 a c d b p2 M1 (231) M1 (231) M2 (312) M2 (312) M0 (123) Frame MSI CRC MSI EC cMS MSI colorectal cancers (CRC); n = 139 MSI endometrial cancers (EC); n = 28 Gene name BANP Type Length %mut %wt<0.5 m3 m2 m1 wt p1 %mut %wt<0.5 m3 m2 m1 wt p1 BANP n = 124 n = 120 n = 137 n = 123 n = 113 n = 127 n = 125 n = 121 n = 111 n = 107 n = 28 n = 24 n = 24 n = 27 n = 28 n = 26 n = 27 n = 25 n = 27 n = 27 CASP5 LTN1 MARCKS MYH11 SLC22A9 SLC35F5 TCF7L2 TGFBR2 TTK CASP5 LTN1 MARCKS MYH11 SLC35F5 SLC22A9 TCF7L2 TGFBR2 TTK M0 (123) p1 wt Mutation m1 m2 m3 m1 m3 m2 m4 m3 m2 m1 wt p1 p2 p3 p4 p1 M2 M0 M1 M1 M2 M0 M1 p2 p1m2 m3 T 12 0.09 0.23 0.24 0.01 0.04 0.02 0.06 0.01 0.01 0.01 0.03 0.01 0.41 0.54 0.08 0.44 0.5 0.29 0.37 0.71 0.71 0.75 0.29 0.42 0.04 0.07 0.00 0.15 0.07 0.00 0.15 0.00 0.09 0.21 0.02 A 10 0.19 0.02 0.35 0.01 0.14 0.24 A 11 0.3 0.01 0.13 0.3 0.03 0.3 A 11 8 0.3 0.05 0.36 0.01 0.19 0.02 C 0.09 0.02 0.34 0.09 0.05 A 11 10 9 10 0.24 0.01 0.08 0.32 0.01 0.31 T 0.27 0.06 0.4 0.21 0.03 A 0.15 0.01 0.4 0.02 A 90.5 0.88 0.54 0.67 0.76 0.4 0.8 0.83 0.63 0.88 0.62 0.41 0.01 0.08 0.41 0.38 A 0.07 0.04 0.3 0.65 0.47 0.58 0.53 0.58 0.58 0.56 0.55 0.59 0.42 0.14 0.85 0.62 0.98 0.76 0.67 0.86 0.78 0.68 0.61 0.68 0.02 Fig. 2 Mutational pattern distribution in cMS based on ReFrame analysis. a Scheme of frameshift mutations (m3, m2, m1, wt, p1, p2) on the left and their corresponding frames (M0, M1, M2) on the right. The numbers 1–3 indicate base triplets in the corresponding frame. Arrows mark the shift between the wt frame and the alternate frame. bDistribution of all cMS frameshift mutations (m3–p2) and their corresponding frames (M0, M1, M2) in MSI CRC quantified using ReFrame. Mutations in all MSI CRC samples were classified in corresponding reading frames (M0, black; M1, magenta; M2, green) and their overall allele ratios quantified. cThe detailed mutational patterns of ten representative cMS are depicted with their respective frequency of mutation for all possible resulting frameshift mutations (m4–p4) in MSI CRC and EC (see Supplementary Fig. 3 for complete dataset). Each row constitutes one analyzed tumor sample with its related allele ratios. For each cMS, tumors were sorted by the proportion of wild-type alleles top to bottom. The number of samples analyzed for a certain candidate is indicated below each candidate’sfigure. Color indicates the resulting reading frames as in b. Intensities represent ReFrame-calculated ratios from white (0%) to full-intensity magenta/green/black (100%) according to the resulting reading frame of the column. The annotated solid lines (first horizontal line top down) show the end of the non-mutated tumor samples (cutoff: 15%), while the dotted lines (second horizontal line top–down) mark the beginning of tumors being >50% mutated, associated with biallelic hits within the respective sample. Source data are available as a Source Data file. dCalculated mutation frequencies and mean allele ratios of most common indel mutation types (m3–p1) resulting from ReFrame analysis in ten representative cMS candidates sorted by length. %mut proportion of mutant tumors, %wt<0.5 proportion of tumors with biallelic hits. Columns m3–p1 display average frequencies of the respective indel mutations over all the samples tested. ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-18514-5 4NATURE COMMUNICATIONS | (2020) 11:4740 | https://doi.org/10.1038/s41467-020-18514-5 | www.nature.com/naturecommunications the immune system in a patient of this population9,34.We calculated GELS for all frameshift peptides using HLA allele frequencies for USA European Caucasians (calculations for additional ethnic groups are provided in Supplementary Data 3). Accounting for a potential relation between immunogenicity and mutation frequency, we noticed that the most commonly mutated cMS located in the ACVR2A gene showed a very low GELS (p mut =91%, GELS =5.1%), whereas very high GELS candidates seemed to be associated with a low mutation frequency (i.e., TMEM97,p mut =27%, GELS =91.1%; SPINK5, p mut =26%, GELS =91.1%; RUFY2,p mut =16%, GELS =90.4%; p binding =50% in USA European Caucasian population; Supplementary Data 3). Hierarchical clustering of cMS candidates on all tumor samples revealed the existence of three distinct populations of cMS (Fig. 4a), which was retained in B2M-wild-type (n=99) but not in B2M-mutant (n=33) tumors (Fig. 4b). Immunoselection during MSI carcinogenesis. In order to systematically evaluate whether these observations may result from immunoediting, i.e.., counterselection of emerging cancer cell clones that harbor highly immunogenic cMS mutations (high GELS frameshift peptides), we analyzed potential differences between the observed and expected distribution of cMS mutations, first in an HLA type-independent approach. Already here we observed a significant inverse correlation between GELS and mutation frequency with Pearson’sr=−0.45, p=0.0078 at n= 41 cMS for endometrial tumors and r=−0.42, p=0.0149 for colon tumors, with a conservative estimate of predicted HLA- binding probability of p binding =50%, indicating that a high GELS was related to lower mutation frequency (Fig. 4c). The correlation remained significant even at the lowest epitope fidelity levels of p binding =10%, with p=0.0145 for endometrial and p=0.0031 for colon cancers. The observation suggests that emerging tumor cell clones with highly immunogenic frameshift peptides are counterselected (Fig. 5), showing that immunoediting may leave its traces in frameshift peptide/cMS mutation patterns in MSI cancers35–38. Interestingly, the significant inverse correlation was only detected among B2M-wild-type tumors. B2M-mutant tumors, in which immune selection on the basis of HLA class I antigen presentation should not apply, only a trend was observed (Fig. 4d), which possibly reflects effects of immune surveillance prior to B2M mutation. To account for confounding factors potentially influencing this observation, we investigated a potential relationship between the cMS length and our observed negative correlation. cMS length is a well-known factor influencing the likelihood of indel mutations on the observed mutation frequency23,27,39 (Supplementary Fig. 5). However, our analysis demonstrated that GELS was not related to the length of the corresponding microsatellite. Moreover, we replicated the correlation analysis with length-adjusted relative mutation frequencies (relative p mut , computed by subtracting the length-specificM1 average mutation frequency from the observed M1 mutation frequency for each microsatellite), the negative correlation between GELS and relative mutation frequency in B2M-wild- type tumors was retained (Supplementary Fig. 8a, b). In order to specifically address HLA-type dependence of immunoediting, HLA-A*02:01 status was determined in MSI CRC and EC40,41. HLA-A*02:01 was present in 34 (47.9%) of 71 MSI CRC and 13 (47.1%) of 27 MSI EC samples, in line with proportions reported for German populations (prevalence of HLA-A*02:01: 46.2%, n=39,689, www.allelefrequencies.net). A significant negative correlation (Pearson’sr=−0.15, p< 0.001) was observed between cMS mutation frequency and the likelihood of HLA-A*02:01 ligands resulting from the respective cMS mutation in the HLA-A*02:01-positive but not in the HLA- A*02:01-negative group of MSI CRC (Fig. 5a). Interestingly, this significant HLA-A*02:01-restricted negative correlation was also observed in MSI EC (Fig. 5b), independently supporting the concept of HLA class I-related immunoediting of MSI tumors in the colorectum and endometrium. Again, replication of the analysis using relative mutation frequencies confirmed the negative correlations (Supplementary Fig. 8c). Despite the statistically significant negative correlation between GELS and mutation frequency, we also observed some outliers (Fig. 4c). We hypothesize that these outliers may reflect distinct effects that potentially influence the probability of a certain cell clone harboring a defined mutation to survive and thrive during tumor evolution. In addition to potential enhancement of immunogenicity, cMS mutations in tumorsuppressor genes are predicted to lead to a growth advantage, at least in cancer or pre-cancer cell clones not directly under Position with FSP [aa] Position with FSP [aa] Coding microsatellite Mutation frequency (M1) Mutation frequency (M2) # Overlapping binding peptides (M1) IC50 < 50 nM IC50 < 500 nM IC50 < 5000 nM IC50 < 50 nM IC50 < 500 nM IC50 < 5000 nM # Overlapping binding peptides (M2) Fig. 3 MHC ligand predictions for a dataset of 82 analyzed frameshift peptides. The figures display the predicted epitopes for the M1 and M2 frameshift peptides derived from all 41 cMS candidates over the length of the respective peptide. All epitopes predicted for three binding affinity thresholds (IC 50 < 50 nM, IC 50 < 500 nM, and IC 50 < 5000 nM) are shown for the M1 frameshift peptides (left, magenta) and the M2 frameshift peptides (right, green). All candidates are sorted in alphabetical order with their respective frequency of mutation according to the ReFrame results. The gray field represents background and indicates C-terminal ends of the respective frameshift peptides. NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-18514-5 ARTICLE NATURE COMMUNICATIONS | (2020) 11:4740 | https://doi.org/10.1038/s41467-020-18514-5 | www.nature.com/naturecom munications 5 attack of the immune system. Such cMS candidates with high GELS and mutation frequencies should be of great relevance for the interaction between the immune system and MMR- deficient tumor cells. To simultaneously account for mutation frequency and GELS as factors influencing the likelihood of frameshift peptides being presented to the immune system, we defined an immune relevance score (IRS), which combines GELS with the mutation frequency in tumors computed via ReFrame (see “Methods”section “Computation of immunological scores”). Gene dissimilarity PCA dimension 1 PCA dimension 2 0 10 20 30 40 50 60 70 80 GELS M1 [%] GELS = 37% GELS = 49% GELS = 63% B2M wt 12 10 8 6 4 2 0 12 10 8 6 4 2 0 PCA dimension 2 Candidate genes Cluster hierarchy by B2M status 0 10 20 30 40 50 60 70 80 GELS M1 [%] B2M mut Gene dissimilarity Candidate genes GELS = 37% GELS = 49% GELS = 63% PCA dimension 1 PCA dimension 2 0 10 20 30 40 50 60 70 80 10 0 20 30 40 50 60 70 80 100 90 GELS M1 [%] Combined cluster hierarchy Candidate genes B2M mut B2M wt IRS M1 [%] 20 50 40 30 20 10 0 10 0 pmut M1 [%] GELS M1 [%] Candidate genes by mutation frequency and GELS correlation by B2M status Pearson's r (GELS M1, pmut) ac b d Agglomerative clustering (k = 3) 0.0 –0.1 –0.2 –0.3 –0.4 –0.5 –0.6 –0.7 Fig. 4 Evidence of immune selection from cMS mutation patterns and GELS. a Hierarchical clustering of cMS candidates on all tumor samples using frameshift abundance features. The full clustering hierarchy is displayed as a dendrogram, showing three clusters (yellow, blue, pink). The same clusters are visualized using RBF kernel PCA with two principal components and colored by their GELS. The three clusters display a trend in their mean GELS, with increasing values from pink to blue to yellow. bHierarchical clustering of cMS on features split by B2M status, with wild type at the top and mutated at the bottom. The full hierarchy is displayed as a dendrogram for both feature sets. In contrast to B2M wild type, B2M-mutated features show no clustering at the computed clustering dissimilarity threshold of 5.7. The same data are again shown using RBF kernel PCA with two principal components. cMutational frequency resulting from one base pair deletions (m1) is shown on the yaxis against the GELS of the resulting M1 frameshift peptides (xaxis). For GELS, all predicted epitopes (IC 50 < 500 nM) were accounted for, with an assumed probability for a binder to be a true positive of p binding =50%. Every bubble represents one candidate. The gradient intensity of the bubbles shows IRS, with white color representing a low IRS, while dark red displays a high IRS. All candidates with an IRS of ≥10% are annotated. dCorrelation between the number of predicted epitopes in cMS mutation-induced frameshift peptides and the frequency of the respective cMS mutations in MSI colorectal cancer separated by B2M mutation status. Pearson’sris shown on the yaxis, while the different groups of tumors are shown on the xaxis. Centers indicate Pearson’sr, whiskers indicate 95% confidence intervals. A significant inverse correlation was observed showing r=−0.42, p=0.0149 at n=41 candidates for 99 MSI colorectal cancers with wild-type B2M, with a conservative estimate of predicted epitope fidelity of p binding =50%, indicating that high GELS was related to lower mutation frequency (for analysis of relative mutation frequencies, see also Supplementary Fig. 8). ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-18514-5 6NATURE COMMUNICATIONS | (2020) 11:4740 | https://doi.org/10.1038/s41467-020-18514-5 | www.nature.com/naturecommunications The M1 frameshift neoantigen derived from TGFBR2, the first described cMS driver mutation in MSI cancer and also the first ever frameshift neoantigen characterized for its immunological properties in MSI cancer in pioneering studies18,42,43, displays the highest IRS (28.57%). In addition to this well-characterized frameshift neoantigen, our study uncovered various frameshift peptide candidates with predicted importance for the immune biology of MMR-deficient cancers. The candidates LTN1, SLC22A9,SLC35F5,CASP5,TTK,TCF7L2,MYH11,MARCKS (all M1), and BANP (M2) all displayed an IRS >10% (Fig. 4c and Supplementary Data 3). The spatial distribution of predicted MHC ligands within these high-IRS frameshift peptides is visualized in Supplementary Fig. 6. To validate the immunogenicity of the predicted HLA- A*02:01-restricted epitopes, HLA-A*02:01-transgenic mice44 were vaccinated with the respective frameshift peptides. Here interferon gamma (IFNɣ) ELISpot demonstrated reactivity corresponding to HLA-A*02:01 epitope prediction for 8 out of 9 tested candidates, confirming processing and presentation of frameshift peptide-derived HLA-A*02:01 ligands for 6 out of 7 candidates (LTN1,SLC35F5,TGFBR2,SLC22A9,TTK,MYH11; Fig. 5c). For CASP5, processing and presentation of HLA- A*02:01-restricted epitopes had been proven previously45. Interestingly, candidate genes with a possible tumor-suppressor function were common among the high-IRS genes: CASP5 (apoptosis induction; IRS: 17.15%), TTK (maintenance of chromosomal stability; IRS: 12.38%), TCF7L2 (beta-catenin signaling; IRS: 11.32%), MYH11 (cell structure and proliferation; IRS: 11.11%), and BANP (migration and invasiveness; IRS: 10.73%) were all previously reported in the literature46–54. This observation may suggest that highly immunogenic frameshift peptides are tolerated preferentially if the cells gain a ab CRC EC c Pearson’s r ** ** # of predicted HLA-A*02:01-binding epitopes HLA-A*02:01 other HLA-A*02:01 other 0.000 0.00 –0.05 –0.10 –0.15 –0.20 –0.025 –0.050 –0.075 –0.100 –0.125 –0.150 –0.175 400 300 200 IFN γ spots per 1 × 106 cells 100 0 MARCKS MYH11.2 TCF7L2 LTN1 SLC35F5 CASP5 TGFBR2 SLC22A9.2 SLC22A9.1 MYH11.1 Positive control TTK Splenocytes—no peptide Splenocytes—irrelevant peptide Splenocytes—specific peptide 000447710121319 Fig. 5 HLA-A*02:01-specific analyses. Correlation between the probability of at least one putative binder truly binding to HLA-A*02:01 within an M1 frameshift peptide and the respective M1 mutation rate in B2M-wild-type HLA-A*02:01-positive and HLA-A*02:01-negative CRC (a) and EC (b). The Pearson’srfrom the correlation test is shown on the yaxis, while the different groups of tumors are shown on the xaxis. Centers indicate Pearson’sr, whiskers indicate 95% confidence intervals. A highly significant inverse correlation was observed showing r=−0.15, p=4.4 × 10−16 at n=34 HLA- A*02:01-positive CRCs with all parameters kept from Fig. 4d. For n=37 HLA-A*02:01-negative CRCs (other), significance did not persist at p=0.29. Similarly, a significant inverse correlation between HLA-A*02:01-binding scores and mutation frequency was detected for HLA-A*02:01-positive MSI EC (n=13, r=−0.15, p=0.00002) but not for HLA-A*02:01-negative MSI EC (n=14) (see also Supplementary Fig. 8c, d). cCellular immune response of M1 frameshift peptides in HLA-A*02:01-expressing A2.DR1 mice. IFNɣELISpot was performed with splenocytes from immunized A2.DR1 mice (n=3) restimulated with the respective, specific M1 frameshift peptides LTN1, MARCKS, SLC22A9, SLC35F5, MYH11, TTK, TCF7L2, CASP5 (magenta), irrelevant peptide (green), or no peptide (black) as controls. Long frameshift peptides SLC22A9 and MYH were synthesized in two parts (split frameshift peptides are marked by asterisks), with the N-terminal peptide labeled 1 and the C-terminal peptide labeled 2. Each dot represents the number of IFNɣspots per 1 × 106cells for one mouse. Mean values are indicated by horizontal lines. The M1 frameshift peptides are sorted according to the number of predicted HLA- A*02:01 ligands (IC 50 < 500 nm). Source data are available as a Source Data file. NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-18514-5 ARTICLE NATURE COMMUNICATIONS | (2020) 11:4740 | https://doi.org/10.1038/s41467-020-18514-5 | www.nature.com/naturecom munications 7 compensatory survival advantage from the mutation by switching off a tumor-suppressive pathway, supporting their role of propelling MSI tumor evolution (Fig. 6and Supplementary Fig. 9). Discussion MMR-deficient tumors, due to their well-defined mechanism of genomic instability, represent an ideal tumor type to study the evolution of solid cancer development and the role of the immune system during this process. By analyzing a broad spectrum of cMS-encompassing genes that are susceptible to mutation in MMR-deficient cells, we were able to identify recurrent mutations and frameshift peptides and to provide first evidence suggesting immunoediting during MSI cancer development. The results of our study demonstrate that, in contrast to neoantigens in many other cancer types, which typically differ between tumors or even occur as “private”mutational neoantigens, MMR-deficient cancers share a large pool of frameshift peptides: the dominance of alterations of the M1 type resulting from one-base pair deletions (m1) emphasizes that MMR- deficient cancers not only share similar sets of genes inactivated by MMR deficiency-induced mutations but also precisely the Driver non-immunogenic Driver immunogenic Bystander non-immunogenic Bystander immunogenic cMS mutation pattern Mutation frequency Low High Positive selectionNegative selection Driver effect Immunogenicity Legend Driver non-immunogenic Driver immunogenic Bystander non-immunogenic Bystander immunogenic cMS mutation pattern Driver effect HLA class I-proficient MSI cancer HLA class I-deficient MSI cancer Immunogenicity B2M mutation Initiation (MMR deficiency) HLA-dependent immune surveillance • tumor-driver effect • low immunogenicity Immune selection Outgrowth of cells Harboring cMS with No immune selection Outgrowth of cells Harboring cMS with • tumor-driver effect cMS mutation pattern a b T cells B2M-deficient MSI cancer cells HLA class I MSI cancer cell FSP epitope Fig. 6 Influence of immunoediting on cMS mutation patterns in MSI CRC. a Coding mononucleotide mutations in MSI cancer can lead to the inactivation of tumor-suppressor genes and therefore have a driver effect on tumor development. Simultaneously, the same mutations can trigger the generation of frameshift peptides with a certain immunogenicity. In the scheme, cMS are located in a two-dimensional lattice depending on the strength of the driver effect resulting from a cMS mutation (yaxis) and the immunogenicity of the respective frameshift peptide (xaxis). Color illustrates cMS mutation frequency, with dark blue representing high mutation frequency and white absence of mutations (according to heat map in Fig. 1a). The combination of positive selection (blue arrows) and negative selection (white arrows) pressure in B2M-wild-type tumors leads to highest frequencies for nonimmunogenic driver cMS mutations and lowest frequencies for immunogenic non-driver cMS mutations (left panel). In B2M-mutant tumors, negative HLA class I-dependent selection is absent; therefore, expected mutation frequencies are independent from immunogenicity of frameshift peptides. bScheme of MSI CRC evolution. B2M mutation can occur as a mechanism of immune evasion during tumor progression. B2M mutations lead to a breakdown of HLA class I-mediated antigen processing and therefore interferes with T cell-mediated elimination of MSI cell clones harboring immunogenic frameshift peptides (also see Supplementary Fig. 9). ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-18514-5 8NATURE COMMUNICATIONS | (2020) 11:4740 | https://doi.org/10.1038/s41467-020-18514-5 | www.nature.com/naturecommunications