Rational Selection of Syngeneic Preclinical Tumor Models for Immunotherapeutic Drug Discovery

Recent clinical successes treating tumors with immunotherapies, including approval of the immune-checkpoint blockade antibodies ipilimumab (anti–CTLA-4) and nivolumab and pembrolizumab (anti–PD-1), demonstrate the potential to transform treatment paradigms and improve patient outcomes (1–3). These treatments represent a shift in the approach to cancer therapy as they do not target the tumor cells but instead target the immune system to circumvent inhibitory pathways that attenuate effective antitumor immune responses. Despite these successes, responses to immunotherapy usually remain restricted to a subpopulation of patients (4, 5). In order to broaden the cancer patient population benefiting from immunotherapy, a greater understanding is needed of the factors that affect response and the potential for combination of different therapies. It is clear that T-cell infiltration varies greatly between individual tumors, patients, and disease types, with some considered to harbor more immunogenic (e.g., “hot”/“inflamed”) tumors (6, 7), characterized by greater T-cell infiltration and Th1 cytokine expression, and overlapping drivers of immunosuppression. In contrast, other tumors may be characterized by a sparse T-cell infiltrate (e.g., immunologically “cold”), potentially as a consequence of reduced immunogenicity. The phenotype of the tumor immune infiltrate correlates with both patient prognosis (8) and outcome following immunotherapy (9). Intrinsic tumor characteristics such as neoantigen load can also affect response to immunotherapy (10, 11), possibly by modulating tumor immunogenicity (12).

Human xenograft tumor models, in which human tumor cell lines are implanted in mice, have played a critical role in understanding traditional cytotoxic or targeted cancer therapies. However, in the context of immunotherapy, these routinely used and well-characterized models are not suitable, because they lack an intact immune system. Several immunocompetent mouse model systems can be used to study immunotherapies, but each brings with it a series of challenges and limitations (reviewed in refs. 13 and14). For example, genetically engineered mouse models (GEMM) recapitulate the anatomical location and encompass some disease-specific mutations frequently observed in human cancer, but they require large colonies of mice, have extended latency periods and, in contrast to the clinical setting, often display limited mutational burden and minimal genetic mosaicism (15). Alternatively, models based on subcutaneous or orthotopic implantation of syngeneic tumor cell lines have short latency periods, are reproducible, and high-throughput; on the basis of this, they have been the workhorse of cancer immunology for several decades (13). These tumor models have been invaluable in providing preclinical proof of concept for candidate immunotherapeutic drugs (16, 17), as well as building an understanding of mechanism of action and evaluating potential biomarkers of response (18, 19). However, to date, the majority of studies have been performed in a small number of models (compared with xenografts, ref. 14) and, despite their widespread use, surprisingly little is known about the genotypes and phenotypes of these syngeneic murine tumor models (14). Ultimately, a better understanding of these models is required to enable appropriate model selection and to permit both data interpretation and extrapolation to the clinic (13, 14).

Here, we describe a comprehensive characterization of the genomic, transcriptomic, and immunologic composition of several murine syngeneic tumor models. Using flow cytometry, we characterized the tumor immune infiltrate, the spleen and tumor-draining lymph nodes (TDLN) in a panel of our most frequently used models. This enabled us to examine whether the genotype and gene expression profile of the tumor cells associated with their immunophenotype in vivo and determine how this related to response following immune-checkpoint blockade. This study highlights the need for extensive characterization of models used for preclinical immunotherapy research, and provides data that will, based on the proposed mechanism of action of the therapy being evaluated, support investigators in selecting appropriate models using hypothesis-driven rationales. It also forms the basis of a dataset that will increase the translational relevance of studies, by allowing parallels to be drawn between models and human disease phenotypes. Moreover, these data can be readily applied to expedite the discovery and development of novel immunotherapies by increasing the efficiency of preclinical drug development.

An overview of the experiments performed is shown in Supplementary Fig. S1. All in vivo experiments were performed in accordance with the UK Animal (Scientific Procedures) Act 1986 and the EU Directive 86/609, under a UK Home Office Project License and approved by the Babraham Institute Animal Welfare and Ethical Review Body, using guidelines outlined by Workman and colleagues (20). C57BL/6 and BALB/c mice were supplied by Charles River UK at 8 to 10 weeks of age and >18g and housed under specific pathogen-free conditions in Tecniplast Green Line IVC Sealsafe cages holding a maximum of 6 animals with irradiated aspen chip bedding, Nestlets nesting material, a cardboard tunnel and wooden chew blocks. Mice were housed on a 12/12 light/dark cycle, with ad libitum UV-treated water and RM1 rodent diet.

Cells (100 μL) in PBS were subcutaneously injected into the right flanks of mice (unless otherwise stated; details of the cell lines and cell numbers are given in Supplementary Table S1). Cells did not undergo any in vivo passaging (except for the B16F10 AP-3 cell line) and were maintained under limited passage from original stocks (typically under 5). We did not undertake additional independent validation. Tumor volume was measured using the formula (width² × length)/2, and tumors were collected when reaching an average of 150 mm³. For the Pan02 cell line, 50 μL of cells in ice-cold matrigel (Corning) were surgically implanted into the pancreas tail. For tumor growth studies, mice were injected intraperitoneally with 10 mg/kg of either anti–CTLA-4 (mouse IgG2b, clone 9D9; Biolegend) or anti–PD-L1 (rat IgG2b, clone 10F.9G2; BioXCell) or the respective isotype controls (αNIP; MedImmune).

Microarray services were provided by Almac Diagnostics. Genomic DNA was extracted from one sample of each mouse tumor cell line using the DNeasy Blood and Tissue Kit (Qiagen) and quality controlled (OD A260/230 ratio ≥ 1.8) using a 2100 Bioanalyzer (Agilent). Samples were 2-color labeled using the genomic DNA labeling kit PLUS (Agilent), hybridized on Agilent Mouse genome CGH 1 × 224K oligo microarrays using the Oligo aCGH Hybridization kit (Agilent) and scanned on an Agilent Microarray scanner. Reference mouse genomic DNA (Merck) was used as a control. All microarray data are deposited into GEO (accession number GSE85509). After normalization to the control mouse genome, probe log₂ ratios for autosomal genes were converted to copy-number variation (CNV) values (Supplementary Table S2 and Supplementary Dataset).

DNA was extracted from cell lines using the DNeasy Blood and Tissue kit (Qiagen). DNA samples were evaluated using an E-Gel (Invitrogen) and PicoGreen fluorometry to measure quality and quantity, respectively. DNA samples were then physically sheared to the desired size using a Covaris E220 Focused-ultrasonicator. Library preparation and enrichment were carried out using an Agilent SureSelectXT Mouse All Exon 49.6Mb design, followed by sequencing on an Illumina HiSeq 2000. Basecall files (*.bcl) were de-multiplexed and converted to fastq.gz format using CASAVA v1.8.2 (Illumina). CrossMap (21) was used to lift the BED files over to mm10 reference for variant calling. The reads were aligned using BWA (22) and variants were called using FreeBayes (23) and VarDict (24). Copy number was inferred from the exome data using Seq2C and CNVKit (25) for approximately 1500 immune system and oncology-related genes. Variants were strain-specifically annotated using data downloaded from ftp://ftp-mouse.sanger.ac.uk/REL-1505-SNPs_Indels/strain_specific_vcfs/ for C57BL_6NJ, BALB_cJ, and DBA_1J. Sequencing data were deposited in the European Nucleotide Archive (accession PRJEB12925).

DNA was isolated from fresh-frozen tumor tissue utilizing the Qiagen MagAttract kit and concentrations determined by Nanodrop spectrophotometry. Focused amplicon sequencing of 64 genes commonly mutated in cancer patients (based on proprietary datasets and the COSMIC database; ref. 26) was carried out utilizing a custom AmpliSeq panel (Life Technologies) to confirm their mutation status. Two pools of amplicons were generated and 20 ng of DNA was utilized for each pool. The manufacturer's protocol was followed for target amplification, adapter ligation, and purification. Libraries were sequenced on an Ion PGM instrument utilizing the 318 chip. Reads were aligned with Bowtie v2.0 (27) using the mouse reference version mm10 (B6 strain). Single nucleotide variants were called using VarScan v2.3.2 (28) with minimum average quality = 30, minimum variant allele read count = 5, and Fisher exact t test P < 0.01. Minimum variant allele frequencies of 7% were called with a minimum of 200× total depth at the locus. Annotation was conducted using Annovar (29).

Microarray services were provided by Almac Diagnostics. RNA was extracted from the tumor, spleen, and lymph node (LN) tissues from each model using RNAStat 60 (Amsbio), and quality of total RNA was determined using a 2100 Bioanalyzer (Agilent). After amplification using the WT-Ovation Pico amplification kit (NuGEN Technologies), cDNA was generated using WT-Ovation Exon module. After quality testing, cDNA was fragmented, labeled and hybridized using the FL-Ovation cDNA Biotin Module V2 (NuGEN Technologies) to Mouse Exon 1.0 ST arrays (Affymetrix). All microarray data are deposited into GEO (accession number GSE85509). Microarray data underwent Robust Multichip Algorithm preprocessing and normalization in Omicsoft ArrayStudio. Exon probes (1,192,934) remained with a probe detection P ≤ 0.05. Quality control excluded two outliers. Log₂ transcript expression values were determined as the median of detected core probes.

To enable cross-comparison of tumor samples between models, five simulated pooled samples (SPS) were generated by proportional resampling with replacement of the log₂ expression intensities for each gene (R and Omicsoft ArrayStudio). The SPS tumor control group was compared with tumor samples from each model to identify differentially expressed genes (DEG) using empirical Bayesian analysis [including vertical (within a given comparison) P value adjustment for multiple testing; Fios Genomics]. LN and spleen samples were compared with their tissue, strain, and sex-matched controls. DEGs were selected using both FDR adjusted P value ≤ 0.05 (Benjamini and Hochberg method) and fold change ≥2.0 or ≤−2.0.

DEGs underwent functional enrichment and network analyses using Ingenuity Pathway Analysis (IPA, Ingenuity Systems). Enrichment analysis results were cross-compared using the IPA comparison tool. The differentially regulated non-disease canonical pathways with the top ten highest and lowest activation Z scores per model were ranked by summed absolute Z score for all 15 models and plotted in a heatmap. Normalized log₂ gene expression intensity per tumor sample for 96 cell type–annotated genes from the Nanostring nCounter PanCancer Immune Profiling Panel gene-list underwent hierarchical cluster analysis (HCA; normalized, linkage = Ward and distance = Euclidean) to cluster samples only, with genes ordered by immune cell type annotation (Matlab).

Tumors were disaggregated using the GentleMACS Mouse Tumor Dissociation kit (Miltenyi Biotech). Spleens and TDLNs were dissociated through a 70-μm nylon cell strainer. Spleens were resuspended in red blood cell lysis buffer (Sigma) for 1 minute. All cells were then stained with a viability dye (Thermo Fisher) and blocked with antibodies to CD16/CD32 (eBioscience) before staining with fluorescence-conjugated antibodies (Supplementary Table S4) in flow cytometry staining buffer with Brilliant Stain Buffer (BD Biosciences). Intracellular staining was performed using the FoxP3/transcription factor staining buffer set (eBioscience) and cells were fixed in 3.7% formaldehyde/PBS. Counting beads (123Count eBeads; eBioscience) were added to the samples before acquisition on an LSRFortessa (BD Biosciences) and analysis using FlowJo (TreeStar). Gating strategies are shown in Supplementary Fig. S2. To enable comparisons between tumors, cell counts were normalized by dividing the cell count obtained using the counting beads by the tumor volume.

Tumors were snap-frozen in liquid nitrogen, then lysed on a TissueLyser II (Qiagen) in 1% Triton-X/PBS with phosphatase inhibitors (PhosSTOP; Roche), protease inhibitors (cOmplete; Roche), and tungsten carbide beads (Qiagen). Lysates were freeze-thawed, centrifuged and diluted to 0.5 mg/mL before analysis of 10 cytokines using the mouse Proinflammatory Panel 1 V-PLEX immunoassay (Meso Scale Discovery). Data were log transformed and plotted using Matlab.

Flow cytometry data were analyzed in GraphPad Prism using one-way ANOVA with Tukey correction for multiple comparisons. For tumor growth studies, group sizes were determined using power analyses based on the variability of the models in pilot studies. Tumor growth data were log₁₀ transformed, and the effectiveness of the therapy (Tr) with respect to the baseline treatment performance was assessed with a linear mixed-effect model (30, 31). The Yij, representing log₁₀-transformed ith tumor volume observed at jth assessment point (T), follows the linear growth model: Yij = a0i + a1i * Tj + eij, where a0i and a1i denote individual intercept and slope parameters, respectively, and eij ∼ N(0,σ) represents model error. Both intercept and slope are assumed to express random effects:

a0i = g00 + g01*Tri + u0i, a1i = g10 + g11*Tri + u1i, with u0i ∼ N(0, σ0) and u1i ∼ N(0, σ1). The parameters g00, g10, and g01, g11 represent the parameter's fixed effects; σ, and σ0, σ1 correspond to intra- and intertumor variance, respectively. Models were defined as “responsive” if the growth kinetics of the treated group compared with the control group was significantly different (P < 0.05).

During our preclinical investigation of CTLA-4 and PD-L1 as targets for immune-checkpoint blockade, we tested the antitumor activity of surrogate antibodies in several murine syngeneic tumor models (32). Significant tumor growth inhibition was seen following anti–CTLA-4 treatment in the CT26 (P ≤ 0.0001) and RENCA (P ≤ 0.0001) models (Fig. 1A), whereas PD-L1 demonstrated activity only in the CT26 model (P ≤ 0.0001; Fig. 1B). This highlighted a need to better understand the underlying immunobiology of murine syngeneic tumor models in order to identify potential drivers of response and enable rational selection of appropriate models for preclinical activity testing.

We performed genome-wide aCGH analysis on 16 in vitro murine tumor cell lines to gain an understanding of broad amplification and deletion events. CNV levels were not affected by the method of tumor cell line generation because CNV frequency of carcinogen-induced CT26, genetically induced TRAMP-C2, and spontaneous LL/2 cell lines were not markedly dissimilar (5.74%, 5.75%, and 5.04% respectively; Fig. 2). No difference in the overall burden of genomic aberrations was seen in the CT26 cell line compared with the others (5.74% of genes studied had CNV in CT26 vs. a mean of 6.05% across the other 15 cell lines). In order to assess CNVs against a paired background strain as opposed to a reference genome, we also investigated somatic gene copy number using WES (Fig. 3A–D). Gain of a single gene copy was frequent (>5% of genes studied) in the RENCA and P815 lines, but amplifications of more than one extra copy much less so. We observed few deletions in the cell lines, although the most common feature across the models was the heterozygous or homozygous deletion of the Cdkn2a tumor suppressor gene (found in 9 of 11 cell lines; Fig. 3B).

We also characterized the somatic mutational status of these genes using WES on a subset of the murine in vitro cell lines and further investigated 64 prominent cancer genes using targeted deep sequencing (Supplementary Table S3 and Supplementary Fig. S3). Our analysis revealed that several tumor cell lines did not carry mutations prevalent in the matched clinical disease. For example, Pan02 is often used as a pancreatic cancer model but lacked mutations in the Kras gene, which are observed in >90% of human pancreatic cancers (33) so it may only be a relevant model for a rare subtype of the disease. CT26 cells harbored mutations in the Apc and Kras genes, which are frequently observed in human colorectal cancer (34), but lacked the frequent mutation in Trp53. In contrast, LL/2 had many of the major mutations found in lung cancer, such as Trp53 and a Cdkn2a deletion (35). The B16 lines bore a mutation in the Braf gene, as is frequently observed in human melanoma (36); however, this mutation was not analogous to the V600E mutation that underlies sensitivity to BRAF inhibitors (37). Comparison of the mutational profile of the murine tumor cell lines with frequently mutated genes found in the cancer type from the same tissue of origin revealed that, although the cell lines carried at least one recurrent mutation seen in the corresponding cancer type, 4 out of 9 models analyzed lacked mutations in the gene most frequently mutated in the corresponding human cancer type (Supplementary Fig. S4). The mutational status of immune-relevant genes was also evaluated across the cell lines to investigate whether cell-intrinsic genetic changes were present that could affect the ability of cell lines to recruit an immune response (Fig. 3D). We found that the most prevalent aberrations affected the complement system, including 6 of the top 50 most mutated or amplified genes in the dataset. Tumor-derived complement has been previously shown to both promote tumor growth but also contribute to the immune surveillance of tumors (reviewed in ref. 38) and even impact the tumor immune infiltrate by affecting the number of infiltrating CD8⁺ T cells and myeloid cells, although these immunomodulatory effects are very context-dependent (39, 40). Several other genes with aberrations in the cell lines have also been shown to affect tumor growth such as fibronectin (Fn1; ref. 41), Cd40 (42), or Nfatc3 (43). This analysis also revealed that CT26 and MC38 (both derived from carcinogen-induced colorectal tumors) had the highest somatic nonsynonymous mutational burden, with 2,955 and 3,018 mutations, respectively (Fig. 3E). ID8, derived from a spontaneous tumor, had the lowest mutational burden (102 mutations), followed by TRAMP-C2, derived from a tumor in a genetically engineered mouse (104 mutations). Thus, mutational burden appeared to be associated with the method of cell line generation in this panel, with cell lines derived from carcinogen-induced tumors having on average a higher mutational burden (1,827 mutations) than spontaneous (1,220 mutations) or GEMM-derived tumors (104 mutations).

Given that clonal selection and evolution can take place during in vivo tumor development, we performed WES of CT26 cells grown in vitro and in vivo (Fig. 3F and G). In the tumor, we detected 2,580 out of 2,955 (87.3%) somatic mutations found in the cell line. A high correlation (r = 0.766) was observed between the mutant allele frequency in vitro versus in vivo (Fig. 3F). Disparities in these data may be due to the use of an unfractionated tumor preparation for this analysis. There were genomic alterations in 57 cancer and immune-relevant genes found in the CT26 cell line which were not detected in the tumor, but these were mostly copy number gains that are more difficult to accurately quantify in heterogeneous in vivo samples (Fig. 3G). There were only four somatic mutations detected in the tumor which were not in the cell line (Erbb3, Ikbkb, Mnx1, and Notch2). These could have arisen due to selective pressures at play in the tumor microenvironment.

To investigate the baseline characteristics of the syngeneic tumor models that could underlie their differential responsiveness to immune-checkpoint blockade, we analyzed the transcriptomes of untreated established in vivo tumors. Pathway enrichment of differentially expressed genes in tumors revealed activation of multiple immune pathways with a high degree of heterogeneity across the models. These included innate immune–related pathways, such as pattern recognition receptors, IL8, TREM1, phagocytosis, and IL6 signaling (Fig. 4A). These pathways revealed reduced activation in tumors from B16 cell lines compared with other models. However, a limited number of pathways (PPAR, RXR, LXR, and PTEN) were more highly activated in the B16 tumors compared with the other models. These data supported the hypothesis that tumors from B16 cell lines are more immunologically “silent” compared with other tumor types, potentially underlying their reduced responsiveness to immune-checkpoint blockade. Hierarchical cluster analysis (HCA) revealed that the pathway activation profile in CT26 tumors was not substantially different from those of other tumor types, suggesting that transcriptomic differences do not completely define response to immune-checkpoint inhibitors (Fig. 4A). We also specifically investigated the expression of MHC class I and II pathway genes in the tumors and found that the B16 tumors, TC1 and LL/2 had less expression of these genes compared with the other models. The TRAMP tumors had low expression of MHC class I pathway genes, whereas the MC38 tumors had low expression of MHC class II pathway genes (Fig. 4B).

To determine whether differences in the tumor immune infiltrate could be affecting responsiveness to immune-checkpoint blockade, we investigated the expression of immune cell type–specific genes. Unsupervised HCA revealed clustering by tumor model and clear changes in expression, suggesting variation in tumor immune infiltrates across the models (Fig. 5). Several findings supported the validity of this approach, such as the elevated expression of a subset of mast cell–specific genes including Kit (20-fold increase compared with the SPS tumor control group), Ctsg (22-fold increase), and Cma1 (2.5-fold increase), which occurred only in the P815 mastocytoma tumors where these genes are expected to be expressed in the tumor cells, as well as potentially in the immune infiltrate. B16 tumors again appeared immunologically “silent” with low expression of immune cell type–specific genes, apart from some NK-cell–specific genes. Elevated gene expression was observed for cytotoxic cell–specific genes Gzma (5.2-fold increase) and Klrd1 (1.8-fold increase) in the CT26 tumors, potentially linking effector function of tumor-infiltrating immune cells and responsiveness to immune-checkpoint inhibitors.

Given the diversity in the expression of immune-related genes, we further profiled the tumor immune infiltrate across both our immune-checkpoint blockade responsive and nonresponsive models by flow cytometry. Tumors were collected at an average volume of 150 mm³, as this often corresponds to the tumor volume at initiation of dosing in activity studies. The immune infiltrate was profiled by staining for 9 nonoverlapping innate and adaptive cell phenotypes and their frequency determined both as a proportion of the CD45⁺ cells and as absolute cell count per mm³ of tumor. Immune cell numbers per mm³ of tumor revealed striking variations across models (Fig. 6A and B). B16F10 AP-3 tumors contained the smallest immune infiltrate with roughly 800 CD45⁺ cells per mm³, whereas LL/2 tumors had the largest infiltrate with 5,500 CD45⁺ cells per mm³ (P < 0.0001). Moreover, the relative proportions of immune cell types also differed profoundly. For example, 4T1 tumors predominantly contained granulocytic myeloid-derived suppressor cells (gMDSCs, >40% of their immune infiltrate), whereas LL/2 and MC38 tumors contained mainly monocytic MDSCs (mMDSCs, 33% and 47%, respectively). Notably, CT26 tumors contained the highest number of NK cells (520 per mm³) and an increased proportion of granzyme B⁺ NK cells compared with other models (Fig. 6C). In addition, the CT26 and RENCA models contained the highest proportions of T cells. Within the CD8⁺ T-cell population, the proportion of reinvigorated cells [expressing markers of T-cell exhaustion (PD-1 and Eomes) but having undergone a reversal of this exhaustion leading to upregulation of Ki67 and granzyme B; ref. 44)] was highest in CT26 tumors (Fig. 6D). The absolute cell counts revealed that CT26, and to a lesser extent RENCA, tumor immune infiltrates were rich in cytotoxic immune cells whereas 4T1, B16F10 AP-3, LL/2, and MC38 are predominantly composed of cell types considered to be immunosuppressive.

In tumor lysates, the highest expression of several proinflammatory cytokines such as IL2, IFNγ, TNFα, and IL1β were found in the CT26, MC38, and RENCA models (Fig. 6E). In contrast, higher levels of Th2-associated cytokines such as IL4 and IL10 were found in the MC38, 4T1, and B16F10 AP-3 models. Notably, KC/GRO, a chemoattractant for neutrophils and MDSCs, was elevated in the MDSC-rich 4T1, MC38, and LL/2 tumors.

To expand our observations, we also assessed spleens and tumor-draining lymph nodes from these animals (tissues from non-tumor–bearing mice were included to highlight changes from baseline). An >8-fold increase in splenic gMDSCs/neutrophils was observed in 4T1 tumor-bearing mice compared with tumor-naïve mice (P = 0.0028; Fig. 7A), mirroring the predominance of these cells in the tumor (Fig. 7E). A common finding was that spleens from tumor-bearing mice had smaller percentages of T cells (30%) relative to tumor-naïve controls (42%; P < 0.0001; Fig. 7A and B). In tumor-naïve animals, a comparison of both genetic strains showed that BALB/c LNs had a larger proportion of CD4⁺ T cells compared with C57BL/6 LNs (49.7% vs. 32.3% P < 0.0001), which had more B cells (13.9% vs. 23.6% P = 0.0003; Fig. 7C and D). In TDLNs, an expansion of B cells was seen versus control tumor-naïve LNs (P = 0.0003) and this was most marked in the LL/2 and MC38 models, where it was accompanied by a corresponding decrease in the proportions of T cells and NK cells. Despite the striking contrasts in the composition of the tumor immune infiltrate (Fig. 7E and F), this was not mirrored in the TDLNs.

Murine syngeneic tumor models play a central role in the advancement of novel immunotherapies; however, there is a need to fully elucidate their distinct molecular and immunologic characteristics. In the current study, we provide a resource to rationalize the selection of syngeneic models to test specific hypotheses; increasing the value of such studies and reducing the numbers of animals used in scientific research. In addition to characterizing the tumor microenvironment, we provide a molecular annotation of these models to describe the specific driver mutations that can be leveraged to guide immunotherapy combinations with molecularly targeted therapies. Indeed, such combinations hold considerable promise to broaden the patient population benefiting from immunotherapy.

We have identified models in which the dominant immunologic phenotype of the tumor is associated with myeloid immunosuppression, for example, increased gMDSC and mMDSC infiltration in the 4T1, MC38, and LL/2 models. We have also identified models, particularly the B16F10 AP-3, that are poorly infiltrated by immune cells. In contrast, other models are both rich in tumor immune cell infiltration but also harbor more balanced populations of effector and suppressive immune cell populations, for example, CT26 and RENCA. This immunologic diversity broadly recapitulates the diversity of tumors observed clinically and therefore provides a powerful resource to evaluate different immunotherapies in the context of these distinct microenvironments. An example of this could be the selection of an immunologically “cold” tumor model such as B16F10-AP3 to evaluate the capacity of treatment to transform the tumor microenvironment leading to enhanced immune priming and immune cell infiltration into tumors, for example, with radiotherapy, agonism of pattern-recognition receptors, or through vaccine approaches. Likewise, modeling of immunologically “hot” tumors is also translationally relevant and the CT26 or RENCA models may provide suitable settings to evaluate experimental therapeutics.

Our analysis of the somatic mutation profile of the murine cell lines reveals genetic disparities when compared with the prevalent mutations in the analogous clinical disease. Examples of these include Trp53 in CT26 or Kras in Pan02. However, many of the genes that are commonly mutated in the clinical setting are also mutated across the preclinical models tested, offering opportunities to understand the impact of these mutations in a syngeneic setting. In particular these data will enable researchers to identify models with both the relevant pathway addiction for a defined molecular target, for example, MEK, BRAF (45), and overlay this with a characterized microenvironment to permit pertinent and appropriate combinations of immunotherapy and small-molecule inhibitors to be evaluated.

Our study also highlights differences between the activity of immunotherapies in a preclinical setting and clinical activity in the analogous human cancer type. An example of this is the contrast between the strong preclinical responses to anti–CTLA-4 and anti–PD-L1 treatment in the CT26 model, derived from a murine colorectal tumor, and clinical responses in colorectal cancer in which, with the exception of patients with mismatch repair–deficient disease, limited activity with immune-checkpoint blockade has been observed to date (46). Poor responses to treatment with either anti–PD-L1 or anti–CTLA-4 were observed in the B16 model, derived from a murine melanoma, which contrasts with the encouraging responses reported for melanoma in the clinic (47). However, when these responses are assessed in the context of the distinct tumor microenvironments observed across the preclinical models, our results are broadly in agreement with clinical findings, which demonstrate an association between CD8⁺ T-cell infiltration into the tumor and outcome following immune-checkpoint blockade (9, 48, 49). Indeed, our work also revealed that tumors such as CT26, in which the immune infiltrate contains greater numbers of cytotoxic immune cells such as CD8⁺ T cells and NK cells, were more responsive to immune-checkpoint blockade than those with comparatively little immune or predominantly immunosuppressive infiltrate. Thus, we conclude that the power of the models lies in recapitulating both a distinct phenotype of tumor microenvironment and mutation profile, rather than directly relating to the analogous human disease setting.

An interesting observation was that although the composition of the tumor microenvironment across the syngeneic models tested was strikingly varied, this contrasted starkly with the relative homogeneity in the immune phenotype observed in secondary lymphoid organs such as the tumor-draining lymph nodes and spleens. This homogeneity was also reflected by gene expression, with markedly fewer differentially expressed gene pathways observed in these tissues when compared with those in the tumors (Supplementary Fig. S5). Overall, these data suggest that the greatest differences between models may result from the preferential chemoattraction, retention, and or differentiation of immune cells within the tumor microenvironment rather than from skewing of the systemic immune response, and therefore the likely factors that are driving response reside within the tumor rather than within the peripheral immune compartments. Again, this finding is also in keeping with the majority of data from the clinic to-date, where the most promising predictive biomarkers have been found in tumor tissue, rather than in peripheral blood (49).

Several of our findings are concordant with previous observations, such as the predominant loss of Cdkn2a among the cell lines, which has been described to result from cell culture–selective pressures (50) and is also prevalent across multiple human cancer types (51). Previous studies have also shown that different syngeneic tumor models can respond differently to treatment (16, 49, 52). It has been hypothesized that these differences could be due to variability in immunogenicity of the models (13) but, until now, the similarities and differences between these models have not been comprehensively characterized. The most extensive study, to our knowledge, included 6 commonly used syngeneic models with analysis limited to the expression of 27 immune-related genes and staining for 4–6 immune cell types by IHC and flow cytometry (52). We also saw a high level of agreement in somatic mutation results with previous genetic characterization of the CT26 cell line, with >60% of mutated genes matching those in the study by Castle and colleagues (53, 54). Although another study described fewer somatic mutations in CT26 and 4T1, 80%–90% of those mutated genes matched our results (54) and the relative proportions of somatic mutations was conserved, with >10-fold more mutations in CT26 compared with the 4T1. Likewise, comparison of our somatic mutation results in the ATCC B16F10 line with those in Castle and colleagues (55) showed over 90% concordance in the mutated genes found in both studies including Brca2, Trp53, Jak3, Atm, Pten, and Mdm1. Disparities in somatic mutations between studies could be due to differences in cut-off levels, sensitivity and methodology, as well as cell line divergence. Our profiling of CNV, using both aCGH and WES analysis paired with the BALB/c mouse background, did not show the same large regions of triploidy and tetraploidy in CT26 seen by Castle and colleagues (53). However, for such multichromosomal regions, it is unlikely that this would be due to divergence in the CT26 cell lines, but rather due to methodological differences such as stringency of CNV calling.

In conclusion, we provide extensive characterization of a range of commonly used murine syngeneic models to rationalize model selection based on the biology of the tumor cell and the tumor microenvironment. Moreover, a greater understanding of model biology allows more robust alignment between response, mechanism of action, and the biology of human cancer subtypes that may ultimately improve the efficiency of drug discovery.

P. Brohawn is an employee of AstraZeneca. R. Stewart has ownership interest (including patents) in AstraZeneca. No potential conflicts of interest were disclosed by the other authors.

Conception and design: S.I.S. Mosely, R.C.A. Sainson, R. Leyland, J. Coates Ulrichsen, P. Brohawn, M. McCourt, J.A. Harper, M. Morrow, V. Valge-Archer, R. Stewart, R.W. Wilkinson

Development of methodology: S.I.S. Mosely, R.C.A. Sainson, R. Leyland, A. Watkins, M. McCourt, J.A. Harper, M. Morrow, R. Stewart, S.J. Dovedi

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): S.I.S. Mosely, D.M. Greenawalt, S. Mullins, L. Pacelli, D. Marcus, J. Anderton, A. Watkins, J. Coates Ulrichsen, P. Brohawn, J.A. Harper, R.W. Wilkinson

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): S.I.S. Mosely, R.C.A. Sainson, J.-O. Koopmann, D.M. Greenawalt, M.J. Ahdesmaki, R. Leyland, D. Marcus, A. Watkins, P. Brohawn, B.W. Higgs, V. Valge-Archer, S.J. Dovedi, R.W. Wilkinson

Writing, review, and/or revision of the manuscript: S.I.S. Mosely, R.C.A. Sainson, J.-O. Koopmann, D.M. Greenawalt, A. Watkins, P. Brohawn, B.W. Higgs, J.A. Harper, M. Morrow, V. Valge-Archer, R. Stewart, S.J. Dovedi, R.W. Wilkinson

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): R. Leyland, J. Anderton, S.J. Dovedi

Study supervision: R.C.A. Sainson, J. Coates Ulrichsen, M. McCourt, M. Morrow

Thanks to Athula Herath, Robert Kozarski, and Wen Wu for advice on experimental design and statistical analysis, to Tristan Lubinski for data analysis and to the biological services unit and core tissue culture facility at MedImmune for technical support. We also acknowledge the assistance provided by the Oncology Research team, in particular Olivia Harris, Marianna Papaspyridonos, and Amy Popple, for their input into the development of flow cytometric analysis of mouse tumor models.

All studies were funded by MedImmune and AstraZeneca.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.