Targeting the PD-1/PD-L1 pathway potentiates immunoediting and changes the dynamics of neutral evolution in a mouse model of colorectal cancer

The cancer immunoediting hypothesis postulates a dual role of the immune system: protecting the host by eliminating tumor cells, and shaping the tumor by editing the genome. However, to what extent immunoediting is shaping the cancer genome is still a matter of debate. Moreover, the impact of cancer immunotherapy with checkpoint blockers on modulating immunoediting remains largely unexplored. Here we elucidated the impact of evolutionary and immune-related forces on editing the tumor in a mouse model of colorectal cancer (CRC). We first show that MC38 cell line is a valid model for hypermutated and microsatellite-unstable (MSI) CRC. Analyses of longitudinal samples of wild type and immunodeficient RAG1 knockout mice transplanted with MC38 cells revealed that upregulation of checkpoint molecules and infiltration of Tregs are the major tumor escape mechanisms. Strikingly, the impact of neutral evolution on sculpting the tumor outweighed immunoediting. Targeting the PD-1/PD-L1 pathway potentiated immunoediting and rendered tumors more homogeneous in the MC38 model. The immunoediting effects were less pronounced in a nonhypermutated/MSI- model CT26. Our study demonstrates that neutral evolution is the major force that sculpts the tumor, and that checkpoint blockade effectively enforces T cell dependent immunoselective pressure in a hypermutated/MSI model of CRC.


BACKGROUND
The concept of cancer immunosurveillance, i.e. that the lymphocytes can recognize and eliminate tumor cells was proposed almost 50 years ago 1 , but the definitive work supporting the existence of this process was published 30 years later by the Schreiber lab 2 . In this seminal work an elegant experiment was carried out using a mouse model lacking recombination activating gene 2 (RAG2), i.e. a gene that encodes a protein that is involved in the initiation of V(D)J recombination during B and T cell development. RAG2 deficient mice are viable but fail to produce mature B or T lymphocytes 3 . RAG2 deficient mice developed sarcomas more rapidly and with greater frequency than genetically matched wild type controls and tumors derived from those mice were more immunogenic than those from wild type mice 2 . These findings led to the development of the refined cancer immunosurveillance concept: the cancer immunoediting hypothesis 4 . The cancer immunoediting postulates a dual role of the immunity in the complex interactions of the tumor and the host: the immune system, by recognizing tumor-specific antigens, can not only protect the host by eliminating tumor cells, but can also sculpt the developing tumor by editing the cancer genome and producing variants with reduced immunogenicity.
Cancer immunoediting is more difficult to study in humans, but clinical data from patients with severe immunodeficiencies is supporting the notion that this process exists also in humans 5 . Indirect evidence for the existence of immunoediting for some cancers was provided by calculating the ratio of observed and predicted neoantigens, i.e. tumor antigens derived from mutated proteins 6 . Using similar approach, we recently provided additional data that support the existence of immunoediting in microsatellite instable (MSI) colorectal cancer (CRC) 7 . However, as we recently showed in a pan-cancer genomic analyses, the composition of the intratumoral immune infiltrates is highly heterogeneous and changing during tumor progression 8 , making it difficult to distinguish between genetic, immune, and other evasion mechanisms. Over and above these mechanistic questions on tumor progression, there is an urgent need to investigate cancer immunoediting also in the context of cancer immunotherapy. Cancer immunotherapy with checkpoint inhibitors like anti-CTLA-4 or anti-PD-1/-PD-L1 antibodies are showing remarkable clinical responses 9 . However, one of the biggest challenges is intrinsic resistance to immunotherapy and the development of resistant disease after therapy, i.e. acquired resistance to immunotherapy. As many patients with advanced cancers are now receiving immunotherapy, elucidating the role of cancer immunoediting as a potential mechanism of acquired resistance to immunotherapy 10 is of utmost importance.
Surprisingly, despite the recognition of the cancer immunoediting process and the widespread use of both, mouse models and next-generation sequencing (NGS) technologies, the impact of immunoediting on the cancer genome has not been well characterized. Cancer immunoediting was investigated in a mouse model of sarcoma using next-generation sequencing (NGS) of the tumor exome and algorithms for predicting neoantigens 11 . This sarcoma model showed that immunoediting can produce tumor cells that lack tumorspecific rejection antigens, but how this finding translates into common malignancies remained unclear. Later, two widely used tumor models, a CRC cell line MC38 and a prostate cancer cell line TRAMP-C1 were used to identify immunogenic tumor mutations by combining NGS and mass spectrometry 12 . However, since neither longitudinal samples of wild type and/or immunodeficient mice nor checkpoint blockade was applied, two major questions remained unanswered: 1) To what extent is T cell dependent immunoselection sculpting the cancer genome? 2) How is immunotherapy with checkpoint blockers modulating immunoediting? Quantitative evaluation of immunoediting during tumor progression as well as following therapeutic intervention using checkpoint blockers could not only provide novel mechanistic insights, but might also inform immunotherapeutic strategies that could potentially be translated into the clinic.
We therefore designed a study to investigate immunoediting of an epithelial cancer genome using wild type and immunodeficient mice, NGS, and analytical pipelines to process and analyze the data. We first characterized the genomic and transcriptomic landscape of the mouse colon adenocarcinoma cell line MC38 (mouse colon #38) that was induced by the subcutaneous injection of dimethylhydrazine in C57BL/6 mice 13 , and show that this cell line is valid model for hypermutated/MSI CRC. We then carried out experiments with wild type and immunodeficient RAG1 -/mice with transplanted tumors and analyzed longitudinal samples with respect to the genomic landscape and the immunophenotypes of the tumors. The results show the extent of immunoediting of the cancer genome in this model in relation to other selection processes. Finally, we performed experiments with anti-PD-L1 antibodies using the MC38 cell line and another CRC cell line which is a model for nonhypermutated/MSI-CRC (CT26 (Colon Tumor #26)) and show how targeting the PD-1/PD-L1 pathway modulates immunoediting.

Genomic, transcriptomic, and immunogenomic characterization of MC38 cell line
Functional studies on immunoediting require genetic tools and controls afforded by mouse studies. Since immunoediting has not been quantified using mouse epithelial cancers, we designed experiments with transplanted tumors using the murine MC38 cell line. The MC38 murine CRC cell line is a grade III adenocarcinoma that was chemically induced in a female C57BL/6 mouse and used since then as a transplantable mouse tumor model 14 . Several studies have shown that the cell line is immunogenic and can be used as a model for investigating anticancer immunity and immunotherapy 15,16,17,18 . To characterize the genome and transcriptome of the MC38 cell line, we performed whole-exome sequencing, SNP array analysis, and RNA-sequencing ( Figure 1A). We identified 7581 somatic mutations of which 3099 were nonsynonymous (2917 missense, 179 stop-gained, 3 stop-loss) and 240 indels ( Figure 1B). Of the 7581 SNVs, the majority (6037) were transversions, of which most (3252) were C>A/G>T. Human hypermutated CRC tumors containing POLE mutations show increased proportions of C>A/G>T and T>G/A>C transversions 19,20 . In contrast, it has been shown that the mouse CT26 cell line shows predominantly C>T/G>A SNVs 21 , similar to the primary human non-hypermutated CRC tumors 22 . Analysis of the data using previously published mutational signatures 19 revealed a mutational profile consisting of a combination of signatures including signature for DNA mismatch-repair (MMR) deficiency (Supplementary Figure S1).
We investigated whether known CRC driver mutations are also present in MC38. We found missense mutations in TP53, BRAF, PTEN, and mutations in the TGF beta pathway (SMAD2, SMAD4, ACVR2A, TGFB2, but not TGFBR2). Additionally, BRAF was also mutated, which is frequently associated with the MSI-high phenotype 23 . KRAS was not mutated and there was only one intron mutation in APC, however there was a truncating mutation in AXIN2 which is known to regulate ß-catenin in the Wnt signaling pathway. Recently discovered frequent mutations in SOX9, and ARID1A 22 were also present in the MC38 cell line. SOX9 is a transcription factor that inhibits Wnt signaling 24 and has a role in regulating cell differentiation in the intestinal stem cell niche 25 , whereas ARID1A is involved in suppressing MYC transcription 26 . Four driver mutations of the MC38 cell line correspond to known hotspots in human CRC, albeit at different positions 27 (Supplementary Table 1).
In a previous large-scale genomic study of human colorectal samples three subtypes of colorectal cancer were identified 22 : (1) microsatellite stable tumors (MSS), (2) tumors with MSI due to a DNA mismatch repair system deficiency and (3) hypermutated tumors that harbor mutations in the exonuclease (proofreading) domain of the DNA polymerases Pol δ (POLD1) and Pol ε (POLE). The MC38 data also showed mutations in the mismatch repair genes MSH3 as well as in POLD1, indicating that the MC38 cell line is a valid model to study human MSI and hypermutated CRC. Both MSI and hypermutated CRC are reported to have better prognosis, higher infiltration of CD8 + T cells and respond well to checkpoint blockade therapy 28 , likely due to the high number of neoantigens.
We then characterized copy number variants (CNV) of the MC38 cell line using exome sequencing and SNP arrays. The analysis of the copy number profiles inferred from the exome sequencing data using hidden Markov model algorithm (see Methods) and from the SNP array data were concordant and showed mostly diploid genome with some regions of amplifications and deletions ( Figure 1A). We identified amplifications in the regions that contain the genes MYC and ERBB2. Finally, we carried out transcriptomic analysis of the MC38 cell line in comparison to normal skin tissues. The transcriptomic data was used to: 1) identify pathways that were up-or downregulated in the cell line, and 2) to identify expressed tumor antigens including neoantigens (identified using exome sequencing data and prediction algorithm as previously described 29 ) and cancer-germline antigens (CGAs). The latter are tumor antigens that are considered to be tumor specific since these molecules are expressed only in germline cells and in 6 tumor cells. Pathway enrichment analysis identified pathways related to cell cycle, DNA replication, DNA repair, and metabolism of nucleotides (Supplementary Figure S2).
With respect to the tumor antigens, we identified a large number of expressed neoantigens ( Figure 1C) and expressed CGAs ( Figure 1D), which provides evidence for the immunogenicity of this model. Of the 3096 amino acid changes (missense and stop codon) in MC38, 1529 neoantigens were predicted to strongly bind to the C57BL/6 MHC class I molecules H2-K b and H2-D b with < 500 nM, and of these, 564 were in expressed genes. Additionally, several CGAs were highly expressed in MC38 including ATAD2, RQCD1, SPAG9, PBK, CTAGE5, CASC5, CEP55, which were also found to be expressed in the CT26 cell line 21 . It is noteworthy that these CGAs were also expressed in the skin samples.
Thus, the characterization of the genomic and transcriptomic landscape of the CRC MC38 cell line demonstrates its validity as a model for hypermutated and/or MSI colorectal cancer.

Upregulation of checkpoint molecules and infiltration of Tregs are the major tumor escape mechanisms in MC38 model of CRC
In our mouse model used to recapitulate the process of cancer immunoediting, MC38 cells were subcutaneously injected into wild type C57Bl/6 and immunodeficient RAG1 -/mice. The tumor growth was monitored regularly and samples were collected at predefined time points and subjected to detailed analysis using FACS, exome and RNA sequencing, and SNP array analysis ( Figure 2A). As expected, the tumor growth was significantly accelerated in RAG1 -/mice compared to the wild type mice ( Figure 2B).
FACS analysis revealed infiltration of both innate and adaptive immune cells including CD8 + T cells, NK cells, and M1 macrophages in wild type mice, that increased with time, although not significantly ( Figure  2C). RNA expression profiles revealed higher expression of chemoattractant molecules such as CXCL9 and CCL5 in wild type mice in comparison to immunodeficient mice (Supplementary Figure S3). However, despite the presence of tumor infiltrating lymphocytes and the slower growing tumors, the adaptive immune system failed to eliminate the tumors. Tumors may utilize several mechanisms of escape such as antigen loss, upregulation of inhibitory molecules, downregulation of major histocompatibility molecules (MHC), or by creating an immunosuppressive environment. The CD8/Tregs ratio, which is a surrogate marker for suppressive tumor microenvironment, was higher in the skin samples compared to the tumor samples at day 23 ( Figure 2C) suggesting that one escape mechanism in this model is the presence of immunosuppressive cells. The numbers of MDSCs and Tregs were comparable in both time points in wild type mice, whereas the M2 macrophages were significantly reduced. The tumor progression in wild type samples was associated with upregulation of immunoinhibitory genes, including PD-1, CTLA-4, TIM3, and LAG3 (Supplementary Figure S4). MC38 expressed low levels of PD-L1, whereas PD-L1 was slightly upregulated in RAG1 -/and more in wild type mice. Although we cannot exclude bethedging (i.e. only clones with higher PD-L1 expression survive) due to the lack of single-cell sequencing data, our finding is in accordance with previous studies showing that PD-L1 expression of MC38 transplanted tumors increases as a result of exposure to inflammatory cytokines such as INF-gamma, which is sufficient for tumor escape and immune evasion 30,31 . Analysis of the differentially expressed genes with respect to overrepresented pathways in wild type vs RAG1 -/tumors showed upregulation of several immune processes related to activation of adaptive immune system response such as costimulation by the CD28, PD-1 signaling, antigen processing and presentation, NK cell mediated cytotoxicity, TCR signaling and interferon gamma signaling ( Figure 2D and Supplementary Figure S5A). Downregulated pathways and GO terms included processes related to cell cycle, DNA replication and TNF signaling (Supplementary Figure S5B).
These data indicate that two tumor escape mechanisms are activated in this model: infiltration of immunosuppressive Treg cells and upregulation of inhibitory genes.

Neutral evolution outweighs T cell dependent and T cell independent immunoselection during tumor progression
Tumor progression is an evolutionary process under Darwinian selection 32 , a characteristic that has been attributed as the primary reason of therapeutic failure, but also as a feature that holds the key to more effective control. At the time of detection, a tumor has acquired novel somatic mutations of which only a small subset (drivers) has an evolutionary advantage. The immune system exerts also an evolutionary pressure through a T cell dependent immunoselection process by acting on a tumor clones that displays strong rejection antigens 11 , and to some extent by T cell independent immunoselection through M1 macrophages, IFNγ, and NK cells 33 . In addition to the ongoing evolutionary and immune-related clonal selection, recent study using a theoretical model demonstrated the occurrence of neutral evolution during tumor development 34 . According to this model, tumor heterogeneity in some cancers including CRC can be explained by neutral expansion and the accumulation of passenger mutations without selective sweeps.
To elucidate the impact of immunoselection on the progressing tumor, we used exome sequencing to identify nonsynonymous mutations and MHC class I binding algorithm to predict the corresponding neoantigens. The number of shared mutations was similar between the biological replicates, suggesting that the sampling bias is rather small. Analysis of exome sequencing data showed high number of mutations that were shared between the MC38 cell line and the two consecutive time points in both, wild type (2919) and RAG1 -/-(2942) samples ( Figure 3A). The number of newly generated mutations was about eight-to ten-fold higher than the number of potentially targeted mutations in both wild type and RAG1 -/mice ( Figure 3A). For example, in the wild type sample at day 23 there were 386 newly generated mutations compared to 50 mutations shared only with MC38 cell line.
According to the cancer immunoediting hypothesis, the immune system can sculpt the developing tumor by editing the cancer genome and thereby modifying the heterogeneity of the tumor: strong immunoediting would render tumors more homogeneous by eradicating immunogenic clones. In order to analyze the heterogeneity of the tumor during progression, exome sequencing data and SNP array data was used to estimate cancer cell fractions (CCF) of all point mutations and subsequently tumor heterogeneity. Analyses of the tumor heterogeneity did not reveal large differences during progression in both, wild type and RAG1 -/samples ( Figure 3B). Strikingly, the analyses showed that the variant allele frequencies (VAF) of the majority of the mutations did not change with time in both the wild type and in the RAG1 -/mouse. On average, 95% of the mutations in the wild type and in the RAG1 -/samples did not change their VAF (Supplementary Table 2), suggesting that neutral evolution rather than Darwinian evolution is driving the tumor growth in this model.
We then characterized the neoantigens using exome sequencing data (to derive somatic mutations), RNAsequencing data (for filtering expressed mutations) and an algorithm for predicting MHC binding (see Methods). In order to identify immunogenic mutations, we filtered the expressed neoantigens with the highest binding affinity (IC 50 <500 nM). In a previous study with MC38 cell line seven mutant peptides were identified using mass spectrometry, of which two elicited a T cell response 12 . In our analysis six out the seven peptides were predicted and four of them were detectable from the RNA-expression data ( Figure  3C). The large impact of neutral evolution was evident also in the Venn diagrams for the neoantigens ( Figure 3D). The number of newly generated neoantigens was comparable in all the samples (126 and 129 for the wild type samples at day 23 and 46) and was much higher than the potentially lost or targeted neoantigens (17 in both wild type samples). Since the key value to understand immunoediting is the ratio between expressed neoantigens and total number of mutations, we analyzed the data and show that the ratio was similar across all samples (Supplementary Figure S5A). Finally, to exclude the possibility that clones may have lost neoantigen-generating mutations, we determined the degree of loss of heterozygosity (LOH) in the samples. Although there was an increased number of events in the transplanted tumors in 8 wild type and in immunodeficient mice (Supplementary Figure S5B), there were no LOH events at the genomic positions of the neoantigens, suggesting that no neoantigens were lost due to LOH.
We then focused our analysis on the tumor samples taken at the same time point, day 23 for wild type and RAG1 -/samples and considered neoantigens found both in the MC38 cell line and in at least one of the RAG1 -/tumors ( Figure 3E). There were 530 neoantigens shared by the wild type, RAG1 -/-, and the MC38 cell line samples. About 3% of the neoantigens (17 out of 530) were detectable only in RAG1 -/tumors (Supplementary Table 3), out of which 16 were derived from mutations not detected or eliminated in the wild type tumors. Only one out of the 17 neoantigens was lost because of low expression. The small number of lost neoantigens imply that the impact of the T cell dependent immunoediting in this model is rather modest. Additionally, similar number of neoantigens (19) was detectable only in wild type tumors, suggesting that these neoantigens were edited by T cell independent mechanisms. Upregulation of genes related to NK cell mediated toxicity and IFN signaling further supports this observation (Supplementary Figure S6A). Analysis of the downregulated transcripts revealed genes related to DNA replication and cell cycle (Supplementary Figure S6B).
Heterogeneity analysis showed that all samples including the MC38 cell line were similarly heterogeneous ( Figure 3F and 3G). To infer how the clonal composition changes between samples, we used a Bayesian Dirichlet process to cluster clonal and subclonal mutations. The results showed that the clonal and subclonal clusters were on the leading diagonal of the plots indicating there was no change in the mutational profile and the clonal and subclonal composition between any two samples ( Figure 3F). There was a large percentage of subclonal mutations both in the MC38 and in the individual tumor samples ( Figure 3H and Supplementary Figures S7 and S8). This was evident also from the variant allele frequency (VAF) plots of the mutations found in diploid regions (Supplementary Figure S9).
Overall, the results suggest that neutral evolution outweighs both, T cell dependent and T cell independent immunoselection. Moreover, the temporal variation in subclonal architecture is largely determined by neutral evolution and to a small extent by Darwinian selection pressure.

Targeting the PD-1/PD-L1 pathway potentiates immunoediting and renders the tumors more homogeneous
We next investigated the impact of the strong immunological pressure induced by targeting the PD-1/PD-L1 axis on the cancer genome on the neoantigen landscape, and on the tumor heterogeneity. It has been previously shown that MC38 responds to different immunotherapies 16,18,35,36 . In order to identify neoantigens that would be potential targets of T cells activated by checkpoint blockade therapy, wild type C57Bl/6 mice were treated with antiPD-L1 antibodies or IgG2b antibodies as control. Treatment was started one day after tumor inoculation and then every three to four days. Samples from six tumors treated with anti-PD-L1 and six tumors treated with IgG2b were taken on day 14. Three samples of each group were used for exome sequencing, and three for RNA-sequencing.
Treatment with anti-PD-L1 antibodies reduced tumor growth in the treated mice by 65% compared to the controls ( Figure 4A), which is in line with previous studies showing that MC38 responds well to PD-1/PD-L1 blockade therapy 37,38 . This was reflected also from the RNA-sequencing data by a strong upregulation of IFNγ, perforin (PRF1), and granzyme A and B (GZMA and GZMB), as well as a number of immunomodulators and MHC molecules (Supplementary Figure S10). GO and pathway analysis showed upregulation of immune related processes such as PD-1 signaling, chemokine signaling, cytokinecytokine receptor interaction, and NK cell mediated cytotoxicity ( Figure 4B and Supplementary Figure  S11). Hence, blocking of the PD-1/PD-L1 pathway induces very strong adaptive and to a lesser extent innate mediated anti-tumor activity in this mouse model.
Analysis of the exome sequencing data showed a large fraction of mutations that were shared in all samples (2555) and 305 mutations that were detectable in the control sample and in the MC38 cell line, 9 but absent from the anti-PD-L1 treated samples ( Figure 4C). These mutations are potentially targeted by the immune system following blockade of the PD-1/PD-L1 pathway. A smaller number of mutations were detectable only in the anti-PD-L1 treated samples and the MC38 cell line (52). Overall, in the anti-PD-L1 treated samples the fraction of mutations resulting in expressed antigens was similar to the control sample (about 25%). The ratio of expressed neoantigens and total number of mutations was highest in the anti-PD-L1 treated sample (Supplementary Figure S12). Analysis of the peptides did not show any obvious patterns that could pinpoint rules defining immunogenicity of the mutations (Supplementary Table 4).
A major shift was observed in the fraction of expressed neoantigens from clonal origin in the anti-PD-L1 treated samples ( Figure 4D and Figure S13). The fraction of clonal neoantigens was 8.8, 26.8, and 10.8 in the MC38, anti-PD-L1 treated, and the control tumors, respectively. Tumor heterogeneity analysis revealed more homogenous tumors undergoing treatment with checkpoint blockers compared to the control tumors and the MC38 cell line ( Figure 4E). The same pattern can be observed in the 2-d density plots which show a shift of subclonal mutations in MC38 towards clonality in the anti-PD-L1 samples ( Figure 4F), suggesting clonal expansion because of a selective advantage of subclones.
Additionally, we investigated the effect of targeting the PD-1/PD-L1 axis on the cancer genome using another widely-used CRC cell line, CT26. Previous genomic characterization of this cell line showed mutation in KRAS and lack of mutations in MMR, POLD1/POLE, and BRAF genes, suggesting that this cell line is a better model for non-hypermutated/MSI-human CRC tumors 22 . From the 1172 point mutations in expressed genes in the CT26 cell line, 154 were in epitopes predicted to strongly bind to MHC molecules 22 , showing that the neoantigen burden is about 73% lower compared to MC38 cell line.
Treatment with anti-PD-L1 antibodies reduced tumor growth in the treated mice by 50% compared to the controls (Supplementary Figure S14A), indicating that this model is less sensitive to immunotherapy with PD-L1 blockers compared to the MC38 cell line. Our results are in line with a recent study showing that tumor growth inhibition following treatment with anti-PD-L1 blocking antibodies was twice as large in mice transplanted with MC38 cells compared to mice transplanted with CT26 cells 31 . Analysis of the exome sequencing data showed a large fraction of mutations that were shared in all samples (1745) and 41 mutations that were not detectable in the anti-PD-L1 treated samples and therefore likely to have been targeted by the immune system (Supplementary Figure S14B). A similar number of mutations was detectable only in the anti-PD-L1 treated samples and the CT26 cell line (33). Thus, the fraction of immunoedited mutations compared to all mutations in the CT26 model was fivefold smaller than in the MC38 model, confirming that the CT26 cell line is less immunogenic than the MC38 model.
Similar to the MC38 cell line, there was a decrease in the number of neoantigens from subclonal origin in the anti-PD-L1 treated sample, albeit less pronounced (Supplementary Figure S14C). The ratio of neoantigens and total number of mutations was higher in the transplanted samples (Supplementary Figure  S14D). There was no overlap of the peptides between the CT26 (Supplementary Table 5) and the MC38 model. The tumor heterogeneity in the transplanted samples was decreased compared to the CT26 cell line (Supplementary Figure S14E), indicating that the transplanted tumors were more homogenous. However, the tumor heterogeneity of the anti-PD-L1 and the control sample was comparable (Supplementary Figure  S14E), as would be expected by the similar numbers of targeted mutations (33 vs 41). Finally, there was no detectable shift of subclonal mutations in CT26 towards clonality in the anti-PD-L1 samples (Supplementary Figure S14F and Figure S15).
Overall, the analyses of this experimental data suggest that targeting the PD-1/PD-L1 pathway potentiates immunoediting and changes the evolutionary dynamics from neutral to non-neutral in the MC38 model of hypermutated/MSI CRC. Moreover, this immunotherapeutic intervention renders the tumors more homogeneous, which could possibly explain the development of resistance to checkpoint blockers. The immunoediting effects were less pronounced in the CT26 model, likely due to the less immunogenic nature of this model.

Immunoediting and acquired resistance to PD-1 blockade in melanoma
In order to test the relevance of our findings in human cancer, we analyzed genomic data from a recent study on acquired resistance to PD-1 blockade in melanoma 39 . In this work pretreatment and relapse samples from four patients with metastatic melanoma, which were subjected to anti PD-1 blockade therapy, were analyzed by exome sequencing. Sequencing data showed that two of the tumors developed loss-of-function mutations in JAK1 and JAK2, respectively, which resulted into lack of response to IFNγ. The third tumor had a mutation in the antigen-presenting protein ß2M which prevented the immune system to recognize the tumor, whereas the fourth tumor had no defined mutations which could be associated with the relapse 39 .
Using exome sequencing data, we analyzed the samples taken before therapy and after relapse with respect to the changes of the mutational landscape, the tumor heterogeneity and the clonal architecture. As can be seen in Figure 5A, large fraction of the mutations was detectable in baseline samples and in the relapse samples in all four cases, implicating that the bulk of the mutations were not efficiently targeted. Newly generated mutations ranged between 5% (case 1) and 33% (case 2). Mutations that were potentially immunoedited following PD-1 blockade, i.e. mutation detectable only in the baseline samples ranged between 4% (case 2) and 58% (case 3). Specifically, case 3 appeared to have strong immunoediting effects on the cancer genome.
With respect to the tumor heterogeneity, targeting the PD-1/PD-L1 pathway showed similar trend: relapsed tumors that acquired larger number of mutations became more heterogeneous (case 2 and case 4), whereas the tumor with lower number of acquired mutations became more homogeneous (case 3) ( Figure  5B). The analysis for case 1 did not reveal changes in the tumor heterogeneity likely due to the high number of mutations in both, baseline and relapse sample (1045). Thus, in this case the impact of newly generated mutations on the tumor heterogeneity is rather small. The analyses of the clonal architecture revealed that in all tumors there was a loss of clonal mutations in the relapsed samples compared to the baseline, ranging from 1% (Case 2) to 24% (Case 3) ( Figure 5C). Tumors that became more heterogeneous had increased number of subclonal mutations compared to the baseline (case 2 and case 4). In accordance with the immunoediting hypothesis, the relapsed sample showing strong immunoediting effect (case 3) had the largest decrease of both, clonal and subclonal mutations, and hence, was more homogeneous.
Overall, these results indicate that immunoediting can be associated with acquired resistance to PD-1 blockade in melanoma in some tumors with specific mutational phenotypes. Targeting the PD-1 pathways in these phenotypes seems to broaden the T cell repertoire in a way that both, clonal and subclonal mutations are targeted and subsequently render the tumor more homogeneous. Hence, a clone which is resistant to immune attack will ultimately dominate the population. However, given the small number of cases and the variability of the results, further studies will be necessary to investigate the effects of the checkpoint blockade on the tumor heterogeneity in relapsed tumors and confirm our findings.

DISCUSSION
With the development of immunotherapies with checkpoint blockers as well as other immunotherapeutic strategies including therapeutic vaccines and engineered T cells 40 , the interaction of the tumor and the immune system, and the question how the cancer genome is edited came into focus. Our understanding of the process of cancer immunoediting and its relevance for therapeutic intervention is still incomplete and requires comprehensive genomic analyses of longitudinal samples. Here we characterized for the first time the extent of immunoediting that tumors undergo during progression or as a consequence of the targeting the PD-1/PD-L1 axis. The quantification of cancer immunoediting using a mouse model of common cancer suggests several biological conclusions and has also important implications for clinical translation.
First, neutral evolution outweighs the effects of T cell dependent and T cell independent immunoselection on the cancer genome during tumor progression in the MC38 model of hypermutated/MSI CRC tumors. Neutral tumor evolution was only recently identified using a theoretical model that determines the expected distribution of subclonal mutations, and implies that a large number of new mutations are generated in ever smaller subclones, resulting in many passenger mutations that are responsible for intratumor heterogeneity, but have minimal or no impact on tumor expansion 34 . In this neutral evolution model all the mutations responsible for expansion are present in the founding cell and subsequent mutations are neutral. Analysis of the TCGA data showed that CRC and other cancers were dominated by neutral evolution whereas other cancers were not 34 . It should be noted that there are three major evolutionary forces that influence the tumor progression and can shape the clonal trajectory of the tumors: 1) drift, 2) positive selection for engraftment and 3) negative selection by the immune system. Drift can be due to random sampling of few related clones from the initial population, and would have the effect of increasing the number of the clonal composition of the transplanted tumors. The same effect would be observed in the scenario of positive selection for engraftment because of selection of clones that are able to survive transplantation. Our analyses show the majority of the new mutations in the transplanted tumors, both in wild type and immunodeficient mice, are subclonal, suggesting that the effects of random drift and selection due to transplantation are negligible. Lastly, we provide three lines of evidence that the contribution of the negative selection by the immune system is small: 1) the ratio between the expressed neoantigens and the total number of mutations was not changed between different samples; 2) although the number of LOH events was increased in the transplanted tumors (both in wild type and in immunodeficient mice), there were no LOH events at the genomic positions of the neoantigen-deriving mutations, suggesting that no neoantigens are lost due to LOH; and 3) there was no loss of expression of the shared neoantigens. Collectively, these data support the neutral evolution model in the MC38 model of hypermutated/MSI CRC tumors.
Second, targeting the PD-1/PD-L1 pathway effectively potentiates immunoediting and changes the dynamics of the system from neutral to non-neutral in the MC38 model. Currently, we can only speculate on the underlying mechanisms driving the strong immune response. It has been previously shown that immunotherapy with anti-CTLA-4 antibodies leads to a significant number of newly detected T cell responses 41 , which can be assigned to broadening of the T cell receptor (TCR) repertoire 42 . Our data support this model also in therapeutic strategy that blocks the PD-1/PD-L1 axis. The broadening of the TCR repertoire might be one of the mechanisms of action of anti-PD-1 treatment and could explain the success of immunotherapy in a number of malignancies. Since CTLA-4 and PD-1 have differing immunological effects on circulating T cells, further mouse and human studies are necessary in order to test the hypothesis that broadening of the TCR repertoire is a mechanism that potentiates immunoediting also in a therapeutic strategy that blocks the PD-1/PD-L1 axis. Notably, targeting the PD-1/PD-L1 pathway in a less immunogenic model, CT26, resulted in fivefold smaller immunoediting, and consequently less pronounced effects on the cancer genome. Intriguingly, similar genotype-immune response associations are observed in humans CRC tumors: MSI tumors respond to checkpoint blockade whereas MSI-are refractory 28 . Further studies are necessary to investigate this genotype-immunophenotype relationships and pinpoint genetic drivers of immunoediting, and ultimately provide possible explanation for the resistance to immune checkpoint blockers in MSI-patients despite the fact that tumor-infiltrating lymphocytes represent a strong independent predictor of relapse and survival 7 .
And third, targeting the PD-1/PD-L1 pathway renders the tumors more homogeneous in the MC38 model. While we did not carry out long-term experiments with different dosages and treatment schedules, one implication of this data is that the tumors might eventually become resistant to immunotherapy. We provide also additional data from a human study showing that in some cases tumors that relapse after PD-1 blockade are more homogeneous. Hence, cancer immunoediting could represents one mechanism of acquired immunotherapy resistance in specific mutational phenotypes. However, other mechanisms like epitope spreading 43 (immune responses to secondary epitopes) and immunodominance 44 (dominant epitopes can mask subdominant ones) could counterbalance or ameliorate this effect and thereby determine if the tumor is eradicated or not.
On the cautionary side, the model we have used has certain limitations since it is based on a cell line which has been edited and it does not recapitulate evolution of the tumor as it occurs naturally. However, as shown by others 12 and in this study, the MC38 model is immunogenic and responds to treatment with immune checkpoint blockers, suggesting that the MC38 evolved and acquired mutations that are detected by the immune system. Moreover, the resemblance of the MC38 and the CT26 models with the clinical observations of response of CRC patients treated with immune checkpoint blockers further supports the relevance of the chosen model. More sophisticated approaches using CRISPR/Cas9 technology for introducing mutations that drive spontaneous rejection of the tumor could provide additional insights into the complex evolutionary dynamics and the interaction with the immune system.
Our findings have important implications for basic research studies on mechanisms of resistance to checkpoint blockade and for clinical translation. Most importantly, given that neutral evolution, T cell dependent immunoediting, and T cell independent immunoediting are sculpting the tumor, it is of utmost importance to carry out comprehensive genomic and immunogenomic analyses of pre-and post-treatment samples. Since conventional cancer therapy as well as cancer immunotherapy are altering the genomic landscape, clones that are resistant to therapy might arise and outcompete other clones. Thus, it is an imperative to characterize the used mouse models and the evolutionary forces driving the tumor in order to dissect the contribution of individual components on shaping the cancer genome. Over and above, since MMR deficiency predicts response to checkpoint blockade in 12 different solid tumor types 45 , further studies using genetically engineered mouse models 46 could help to identify molecular determinants that make not only CRC tumors sensitive to immune checkpoint blockade, but possibly also other tumor types.
Finally, our results have important implications also for clinical research. Given the fact that some cancers including CRC, stomach, lung, and bladder are dominated by neutral evolution 34 , it will be important to study tumors over time to dissect out the impact of the immunological selection following checkpoint blockade. Neutral evolution theoretically generates greater tumor heterogeneity and hence, may facilitate adaptation after the initiation of immunotherapy. However, investigating evolutionary dynamics within human cancer is challenging since longitudinal observations are unfeasible and both, the genetic and immune landscape of cancer are highly dynamic and interwoven 8 . Use of new technologies such as single cell sequencing, as well as multi-region sequencing and higher sequencing depth together with improved computational methods as recently shown 47 will provide better understanding of the relationship between the clonal architecture of a tumor and the antitumor response of the immune system. In this context, advances of organoid technologies and gene editing will open new avenues of research and ultimately lead to the development of precision immune-oncology.

13
In summary, we demonstrated that neutral evolution is the major force that sculpts the tumor during progression, and that checkpoint blockade effectively enforces T cell dependent immunoselective pressure in mouse model of CRC. Our study adds another layer of complexity of the tumor evolution and the dynamic nature of clonal selection driven by immunological and non-immunological mechanisms. An improved understanding of how the immune system shapes the tumor progression will be fundamental to improving response to immunotherapies and combating resistance, and will require comprehensive genomic and immunogenomic analyses of both, mouse models and human samples.

Mouse experiments
Wild type C57BL/6N mice RAG1 -/-(B6.129S7-RAG1 tm1Mom /J) mice were purchased from Charles River. Mice were maintained under SPF conditions. All animal experiments were performed in accordance with the Austrian "Tierversuchsgesetz" (BGBI. Nr.501/1989 i.d.g.F. and BMWF-66.011/0061-II/3b/2013) and were approved by the Bundesministerium für Wissenschaft und Forschung (bm:wf). 5×10 4 MC38 colon carcinoma cells were injected subcutaneously (s.c.) into the left flank of 8-to 12week-old female wild type or RAG1 -/mice. Tumor growth was monitored three times per week by measuring tumor length and width. Tumor volume was calculated according to the following equation: ½(length × width 2 ). Each excised tumor was randomly divided in 3 pieces and used for either DNA or RNA isolation or for FACS analysis. For survival analysis, mice with tumors greater than the length limit of 15 mm were sacrificed and counted as dead.
Wild type C57BL/6N mice were injected s.c. with 5x10 5 MC38 scells and administered with 0.5mg of an anti-mouse PD-L1 (Clone10F.9G2; BE0101) or corresponding IgG2b (LTF-2; BE0090) control antibody (all from BioXCell, USA) every 3 to 4 days starting from day 1 of MC38 challenge according to 48 . Tumor growth was monitored as described above. DNA and RNA isolations were done from complete excised tumors.
Wild type BALB/c mice were injected s.c. with 5x10 5 CT cells and administered with 0.5mg of an antimouse PD-L1 (Clone10F.9G2; BE0101) or corresponding IgG2b (LTF-2; BE0090) control antibody (all from BioXCell, USA) every 3 to 4 days starting from day 1 of CT26 challenge. Tumor growth was monitored as described above. DNA and RNA isolations were performed from complete excised tumors.

Immunophenotyping
Mononuclear infiltrating cells were isolated from both subcutaneous tumors and skin tissue at the indicated time points 48 . Briefly, tumor and skin tissues from sacrificed mice were prepared by mechanical disruption followed by digestion for 45 min with collagenase D (2.5 mg/ml; Roche, 11088858001) and DNase I (1 mg/ml; Roche, 11284932001) at 37°C. For skin tissue Liberase (5mg/ml; Roche, 5401020001) was added to the above described digestion mix. Digested tissues were incubated 5 min at 37°C with EDTA (0.5 M) to prevent DC/T cell aggregates and mashed through a 100-µm filter and a 40-µm filter. Cells were washed, and resuspended in PBS+2% FCS.
Tumor and skin infiltrating immune cells were incubated with FcR Block (BD Biosciences, 553142) to prevent nonspecific antibody binding before staining with appropriate surface antibodies for 30 min at 4°C, washed with PBS+2% FCS, and used for FACS analysis. For intracellular cytokine staining, cells were stimulated with 50 ng/ml Phorbol 12,13-dibutyrate (PDBu, Sigma, P1269), 500 ng ionomycin (Sigma, I0634), and GolgiPlug (BD Biosciences, 555029) for 4-5h. After fixation with the FoxP3 staining buffer set (eBiosciences, 00-5523) for at least 30 minutes at 4°C, cells were permeabilized with the fixation/permeabilization buffer (eBiosciences, 00-5523) and incubated with FcR Block (BD Biosciences, 553142) before staining with specific cell surface or intracellular marker antibodies. Data acquisition was performed on a LSR Fortessa cell analyzer (Becton Dickinson). Data analysis was conducted using the Flowlogic software (eBioscience, version 1.6.0_35). Exome-and RNA sequencing DNA sequencing from the tumor, skin and MC38 cell samples was performed by exome capture using SureSelectXT Mouse All Exon capture probes (Agilent Technologies Österreich GmbH, Vienna, Austria) followed by sequencing with the Ion Proton TM System (Ion Torrent, Thermo Fisher Scientific). For RNA sequencing, total RNA was extracted, quality validated with the Agilent Bioanalyzer, and submitted to QuantSeq 3' mRNA-Seq library preparation, following the manufacturers instructions (Lexogen, Vienna Biocenter, Austria). Resulting libraries were sequenced with the Ion Proton TM System.

Exome-sequencing data analysis
The sequence reads were aligned to the mm10 reference genome with tmap (https://github.com/iontorrent/TS/tree/master/Analysis/TMAP) and preprocessed using GATK. Somatic point mutations were identified with Mutect 49 by comparing each tumor sample with the two skin samples and taking the intersection of the mutations. Insertions/deletions were called with Strelka 50 in the same way. The somatic mutations were annotated using the Ensembl Variant Effect Predictor tool 51 . Somatic copy number estimations were derived from exome-sequencing data using EXCAVATOR 52 by calculating log 2 ratios between the read depth of the tumor and two germline skin samples using the "pool" mode. The estimated log 2 ratios were then segmented by their novel heterogeneous shifting level model (HSLM). The CNAs identified using exome-sequencing data were concordant to those in the same samples by using Affymetrix SNP Array. Loss of heterozygosity (LOH) events were derived using VarScan2. MutationalPatterns (https://doi.org/10.1101/071761) was used to infer the contribution of published mutational signatures 19 .

SNP arrays
Genome-wide copy number profiles of two wt samples (day 23 and day 46), two RAG1 -/samples (day 13 and day 23), all six aPD-L1 and IgG2b samples, MC38 and skin germline DNA were obtained using the Affymetrix Mouse Diversity Array. The genotyping analyses were carried out at Eurofins Genomics (Ebersberg, Germany) using the Affymetrix Mouse Diversity Array. The SNP arrays were processed, quantile-normalized, and median-polished using the Aroma Affymetrix CRMAv2 algorithm 53 together with 351 publically available Mouse Diversity Genotyping Array CEL files which were downloaded from the Center for Genome Dynamics at The Jackson Laboratory (http://cgd.jax.org/datasets/diversityarray/CELfiles.shtml). Copy number alterations (CNAs) for each probe were computed as log 2 -ratios between the probe signal intensities of each sample and the reference skin sample and then those ratios were segmented using the circular binary segmentation algorithm implemented in the R package DNAcopy 54 .

Tumor heterogeneity
Normal contamination estimates were calculated using the homozygous point mutations in the cell line MC38. Considering that the purity of the cell line is 1, we checked the variant allele frequency of the homozygous mutations in MC38 in all the samples together with the estimated copy numbers of the corresponding region. The expected VAF of these mutations should be 1 in all samples assuming that there is no normal contamination and no new mutations appearing in the mouse samples at the same genomic position. As an estimate of the purity of the tumor, we took the mean of the VAF of those mutations found in a diploid region. These estimates were used to correct the mutation VAFs or copy number estimates in the rest of the analyses.
Mutations were filtered so that only mutations with at least 10 total reads and at least 5 alternative reads were considered. The CCF of each mutation was calculated using the approach of Yates et al 55 . Briefly, for each mutation the observed mutation copy number, n mut (the fraction of tumor cells carrying a given mutation multiplied by the number of chromosomal copies at that locus) was calculated as: where VAF is the variant allele frequency of the mutation, p is the tumor purity, and CN t and CN n are the tumor and the normal locus specific copy number. Since mutations that are present of multiple chromosomal copies will have a mutation copy number higher than 1, we determined the number of chromosomes that the mutations is residing on. This was done so that for all mutations in amplified regions with a copy number of CNt, the observed fraction of mutated reads is compared to the expected fraction of mutated reads resulting from a mutation present on 1,2,3,…,CNt copies, considering a binomial distribution. The cancer-cell fraction was then calculated as the mutation copy number divided by the value of C with the maximum likelihood.
Mutations were defined as clonal if the CCF was > 0.95, and subclonal otherwise. Subclonal clusters of mutations were identified using a previously described statistical modelling of the distribution of clonal and subclonal mutations by a Bayesian Dirichlet process [56][57][58] .

RNA-seq data analysis
The sequencing reads were first preprocessed through a quality control pipeline consisting of adapter removal with Cutadapt (DOI: http://dx.doi.org/10.14806/ej.17.1.200) and quality trimming with Trimmomatic 59 to remove bases with bad quality scores and reads shorter than 22 nucleotides. The quality trimmed reads were then mapped to the mm10 reference genome using a two-step alignment method; alignment with STAR 60 followed by alignment of the unmapped reads with Bowtie2. From the reads that mapped to multiple locations in the genome only the primary alignment was retained. Reads that mapped to ribosomal RNA locations in the genome were removed from further analysis using the split_bam.py script from the quality control package RSeQC 61 . Gene-specific read counts were calculated using HTSeq-count 62 . The R package DESeq2 63 was used for differential expression analysis. The p-values were adjusted for multiple testing based on the false discovery rate (FDR) using the Benjamini-Hochberg approach.

Neoantigens and cancer-germline antigens
All possible 8-11 mer mutated peptides generated from all the nonsynonymous mutations (missense and nonsense) were used as an input to netMHCpan to predict their binding affinity to the C57BL/6 MHC class I alleles H-2K b and H-2D b . Amongst the candidate antigenic peptides, only the strong binders with binding affinity < 500 nM, and peptides arising from expressed genes were retained for further analysis. A mutation was considered expressed if the normalized counts of the corresponding gene were greater than 5.
The list of cancer-germline antigens (CGA) was downloaded from the Cancer-Testis database 64 . Their expression level was estimated using the normalized counts from DESeq2.

Human data
Mutational data from the melanoma patients 39 were provided by Dr. Antoni Ribas. Heterogeneity analyses were performed as described above.

Data accessibility
The mouse data sets were deposited in the GEO (GSE93018) and the Sequence Read Archive (SRP095725).

Statistical analysis
For comparison of two sample groups two-tailed unpaired Student's t-test was performed. Analysis and visualization of Gene Ontology terms and pathways associated with differentially expressed genes was performed using ClueGO 65 . A p-value of <0.05 was considered statistically significant: *p<0.5; **p < 0.01; ***p < 0.001.