Definition of the neoepitope Ccdc85cMUT. The Ccdc85c gene encodes a gap junction protein expressed mostly in the brain, colon, lung, kidney and testes in adult mice. The protein has no known oncogenic (driver) function. A non-synonymous (leucine to phenylalanine, Chromosome 12-108221754) somatic SNV in Ccdc85c was detected in the BALB/cJ Meth A fibrosarcoma (Fig. 1a). The mutation is heterozygous and the un-mutated as well as the mutated reads are detected in the transcripts. BALB/cJ bone marrow derived dendritic cells (BMDCs) pulsed with an 18-mer peptide with the mutant amino acid near the center (DPSSTYIRPFETKVKLLD) or un-pulsed BMDCs, were used to immunize BALB/cJ mice as described in Methods. All mice (n = 225) were challenged with the Meth A cells, and tumor rejection was monitored. A tumor rejection score (TRS) (with a maximum score of 5 indicating near 100% tumor rejection) was used to quantitate the extent of tumor rejection as described in Methods. The 18-mer peptide elicited a perfect 5.0 TRS score (Fig. 1b). Various truncated versions of this peptide as indicated in Fig. 1b were similarly tested for tumor rejection. The most and the least effective peptides along with their TRS scores are shown in Fig. 1c. Since the 10 amino acid peptide YIRPFETKVK was the shortest peptide active in tumor rejection, we consider this the precise epitope.
Evidence of presentation of YIRPFETKVK was sought by analyzing the peptides eluted from MHC I molecules purified from the Meth A cells by mass spectrometry (MS), as described in Methods. No Ccdc85c-derived peptides were detected, as expected from the low abundance of expression of this protein. In order to identify the precise peptide derived from the mutant Ccdc85c that could be cross-presented by the DCs, BMDCs were pulsed in vitro with the 18-mer peptide as previously described 11. The BMDCs were extensively washed and MHC I molecules eluted. Targeted-MS analysis of the eluted peptides in the presence of spiked-in heavy labeled synthetic peptides showed the presence of two Ccdc85c-derived peptides TYIRPFETKVK and YIRPFETKVK (Fig. 1d). These two peptides detected by cross-presentation of the 18-mer peptide were identical to the two truncated versions of the 18-mer peptide that were observed to be the most effective in tumor rejection (Fig. 1c). These peptides had very low or undetectable predicted as well as measured affinities for Kd, Dd and Ld as shown for Kd in Fig. 1e. With such low affinities, these neoepitopes would normally be considered non-binders.
The 18-mer sequence was queried for the presence of predicted Kd, Dd or Ld-binding peptides. No Dd or Ld-binding peptides were predicted; three peptides were predicted to bind Kd albeit with poor affinity (IC50 values between 692 and 864 nM) (Fig. 1e). Ironically, none of these three peptides were detected by MS among peptides eluted from MHC I of BMDCs pulsed with the 18-mer long peptide.
Immunization of BALB/cJ mice with the un-mutated Ccdc85c 18-mer peptide failed to elicit protection from tumor growth (Fig. 2a). The anti-tumor activity of Ccdc85cMUT was abrogated by depleting the mice of CD8 cells by treating the mice with the anti-CD8 antibody but not by a control antibody during the priming phase as previously described 11.
In order to determine if any peptides within Ccdc85cMUT could be presented by MHC II molecules, we analyzed the interaction of H2-Ad and H2-Ed with TYIRPFETKVK, YIRPFETKVK and IRPFETKVK as well as their wild type counterparts, using a cell-surface density assay. In this assay, β-chains of H2-Ab1d or H2-Eb1d are expressed in fusion with the peptide of interest and the amount of cell-surface MHC II, as a measure of intrinsic stability of MHC II, is quantified in engineered conditions 18. No significant difference was observed in binding of Ccdc85cMUT and Ccdc85cWT to H2-IA or H2-IE (Supplementary Table 1 and Supplementary Fig. 1).
Analysis of immune response against Ccdc85cMUT. Tumor microenvironment (TME) of Ccdc85cMUT-immunized mice was examined using single cell RNA sequencing (scRNA seq). BMDCs- immunized mice were used as controls. As a peptide control, mice immunized with a neoepitope Alms1MUT were used. During comparison of the sequences of Meth A exomes with the normal BALB/cJ exomes, we identified Alms1MUT (LYLDSKSDTTV) which was also identified among the peptides eluted from MHC I molecules from BMDCs pulsed with the 18-mer Alms1MUT peptide. Although this neoepitope has a high affinity for a mouse MHC I Kd (IC50 62.25 nM), and it can be processed and presented, immunization of mice with the 18-mer Alms1MUT failed to elicit tumor rejection (Supplementary Fig. 2). Hence, this peptide was chosen as a control peptide. Mice were immunized with Ccdc85cMUT, Alms1MUT (peptide control) or BMDCs (control) as described in Methods, and were challenged with Meth A. RNA from CD45 + cells (estimated 15,925 cells for the four libraries after QC, with an average coverage of 32,663 reads per cell and median 1,339 genes per cell) was sequenced. Data from the four libraries were pooled and clustered based on the gene expression pattern of each as described in Methods. Annotation of the clusters was informed by both differentially expressed genes (DE Genes) and per cluster highly expressed genes identified by the TF-IDF analysis as described in Methods. Two major and distinct clusters were identified namely, myeloid and lymphoid. Upregulated genes used to identify the myeloid cluster included, but were not limited to: Itgam (CD11b), Adgre1 (F4/80), Arg1 (Arginase 1), Nos2 (Nitric oxide synthase 2), Ms4a4c (Membrane-spanning 4-domains, subfamily A, member 4C), C1qa (Complement component 1) were highly expressed in the myeloid cluster. Upregulated genes used to identify the lymphoid cluster included, but were not limited to: (Cd3 (CD3), Ptprcap (CD45-AP), Nkg7 (Protein NKG7), Cd28 (CD28), Gzma (Granzyme a), Prf1 (Perforin1) etc.) (Fig. 2b, left panel). The proportion of lymphoid versus myeloid compartments in the Ccdc85cMUT library was very different compared to the Alms1MUT or control groups. The Ccdc85cMUT library was mostly composed of lymphoid cells (62.68% lymphoid and 37.32% myeloid), while the control groups and mice immunized with Alms1MUT were mostly composed of the myeloid compartment (~ 35% lymphoid and ~ 64% myeloid) (Fig. 2b, right panel). In order to study different cell types, lymphoid and myeloid clusters were re-clustered into 6 and 9 sub-clusters, respectively (Fig. 2b, bottom panel). The six identified lymphoid clusters were: CD4 T cells (CD4(1)), NKC(1), naive/early activated CD4 T cells (CD4(2), defined by a high expression of Sell, Il7r, Tcf7 and Ccr7 genes and low expression of Il2ra and lack of expression of effector and cytotoxicity genes), NKC(2) (less cytotoxic and active than NKC(1)), CD8 T cells (CD8) and proliferating CD4/CD8 T cells (Pr. CD4/CD8, defined by higher expression of Stmn1 and Mki67 genes and cell cycle gene expression analysis, described in Methods). The selected genes used as markers to annotate each lymphocyte cluster are listed in the summary heat map (Fig. 2c, right panel).
The majority of NK Cells (~ 80%) in the Ccdc85cMUT library were from the NKC(1) cluster which was the more cytotoxic and active cluster (defined by higher expression of Cd44, Tnfa and Il7r), while, the fraction of active NK cells in other libraries was about 55%.
To pinpoint differences in T cells of the four libraries, clusters 1 and 5 (activated CD4 and CD8 T cells) were computationally pooled and the expression of cytotoxicity and other effector function genes were compared between libraries. Proliferating CD4/CD8 T were excluded from further analysis because the gene expression levels in these cells could be influenced by the cell cycle effect prominent in this cluster. Interestingly, Ccdc85cMUT library had the most contribution to the aforementioned pooled cells (31% Ccdc85cMUT, 25% Alms1MUT, 21% Ctlr1, 20% Ctrl2). Also, the normalized average gene expression (described in Methods) of cytotoxicity (Gzmb, Prf and Nkg7) and other effector function (Ifng) genes were significantly higher in T cells derived from the Ccdc85cMUT library compared to the control or Alms1MUT libraries. Similarly, T cells of Ccdc85cMUT library had a significantly higher expression of genes involved in TCR engagement (Nr4a1 and Irf4). A transcription factor involved in transcription of cytotoxicity genes, Eomes, had a significantly higher expression in T cells of Ccdc85cMUT library (Fig. 2d).
In the myeloid compartment, nine distinct clusters were identified. These are: macrophage1 (Mφ1), Mφ2 (defined by a moderate expression of Arg1 and lower expression of Cd302, Ccl5, Ccl8 and C1qa), monocyte1 (Mo1), Mφ3 (defined by a lower expression of Ccl8 and a higher expression of Ly6c, Cxcl9, Il1b, H2-Ab1, H2-DMb2, Mmp14 and Cd38), Mφ4 (defined by a higher expression of Nos2, Mrc1, Itgam, Pf4, C1qa, C1qb and C1qc), DC1, Mo2 (defined by a higher expression of Itgax, Tlr7, Ace, and Adgre4), neutrophil (Ne) and DC2 (defined by a higher expression of Ccr7, Ccl5, Samsn1, Pcgf5, Gyg, Net1 and Rabgap1l) (Fig. 2c, left panel). Contribution of library Ccdc85cMUT to the most of myeloid clusters were minimal (ranging from 1.6% in Mφ3 to 13.4% in Mφ2). The only exception is for Mφ1 which Ccdc85cMUT library forms ~ 20% of this cluster (Supplementary Fig. 3).
Mutation-reversion analysis of CD8 T cell immunogenicity of Ccdc85cMUT. The studies described above examined the immunogenicity of Ccdc85cMUT when administered as a vaccine. We aimed now to analyze the role of Ccdc85cMUT in the immunogenicity of the Meth A tumor itself, and in vivo. The broader objective was to test if a poorly MHC I-binding neoepitope residing within a tumor influences the CD8 immunogenicity of the tumor. CRISPR-guided gene editing was used to generate two variants of the Meth A. For purpose of this experiment, we refer to the original Meth A tumor with the endogenous Ccdc85cMUT as MUT1. Using CRISPR, the point mutation in Ccdc85cMUT was reversed back to its WT counterpart as described in Methods and Supplementary Fig. 4; this line with the WT sequence of Ccdc85c is referred to as a Revertant (REV). The point mutation in Ccdc85cMUT was then re-introduced into the Revertant to generate the line MUT2 which recapitulates the original mutation as in Ccdc85cMUT. Thus, MUT1 and MUT2 are identical tumors (except that MUT1 is heterozygous for the mutation while MUT2 is homozygous for it) and are different from REV only with respect to the mutation in Ccdc85cMUT. Groups of mice were challenged with the three tumors individually (MUT1, REV or MUT2) and their TILs analyzed by scRNA sEq. RNA from CD45 + cells (estimated 11,961 cells for the three libraries before QC and 10,265 after QC, and with an average coverage of 61,371 reads per cell and median 2,137genes per cell) was sequenced. Data from three libraries were pooled and expression of the top average TF-IDF scoring genes was compared between the three libraries (Fig. 3a). CD8 T cells of all three libraries showed CD8 activation markers including but not limited to: Cd69 and Cd44 as activation markers, Lamp1 (CD107a) as a measure of degranulation and Tbx21 (Tbet), a transcription factor involved in transcription of cytotoxicity-associated genes (Fig. 3b, upper panel). We then compared the gene expression patterns of the total immune cell population in TILs of MUT1, REV and MUT2. The gene expression patterns in TILs of MUT1 and MUT2 showed a higher similarity to each other than to the REV: a simple hierarchical clustering (using Euclidean distance, and complete linkage) of the MUT1, REV and MUT2 libraries represented by the normalized average expression vector of top informative genes (selected by highest average TF-IDF score, see Methods) showed MUT1 and MUT2 are closer to each other (distance 1.097) than to the REV (distance 1.820) (Fig. 3a upper panel). In Fig. 3a bottom panel, where genes with variability among the three libraries are juxtaposed, it is clear that the difference between REV and MUT1/MUT2 is more pronounced than the difference between MUT1 and MUT2.
To identify differences in CD8 T cells of the three libraries, the RNA sequencing data of the combined libraries were clustered based on the gene expression pattern of each cell type as described in Methods. Annotation of the clusters was informed by both differentially expressed genes (DE Genes) and per cluster highly expressed genes identified by the TF-IDF analysis. T cells were re-clustered into 7 clusters by unsupervised clustering as described in Methods and enriched CD8 T cells populations were further analyzed (Supplementary Fig. 5a-b). Some cytotoxicity and effector function genes (Tnf (TNF), Gzma (Granzyme a) and Ifng (Interferon gamma)) had similar expression pattern in CD8 T cells of all three libraries; however, the normalized average gene expression of other cytotoxicity genes (Fasl (Fas ligand), Gzmb (Granzyme b) and Prf1(Perforin1)) as well as other effector function genes (Pcdc1 (PD1) and Tbx21(Tbet)) were significantly higher in CD8 T cells derived from MUT1 and MUT2 libraries compared to the REV library (P value < 0.001). Similarly, CD8 T cells of MUT1 and MUT2 libraries had a significantly higher expression than the REV library, of early response genes which are involved in TCR engagement (Nr4a1, Nr4a2 and Nr4a3, P value < 0.001) (Fig. 3b, bottom panel).
T cell receptors (TCRs) in the TILs of MUT1, REV and MUT2 tumors. T cell receptors (TCRs) in the TILs of the three libraries were characterized using Grouping of Lymphocyte Interactions by Paratope Hotspots (GLIPH) analysis that groups together the TCRs into specificity groups based on the global and local similarities of the CDR3 regions of the TCRs 19. Based on the GLIPH algorithm, 40 to 42.9% of all distinct clonotypes contributed to forming a network/similarity-based specificity groups in each of the libraries, while the rest were standalone clonotypes (with no similarity to other clonotypes). The similar percentage of network-based specificity groups (40- 42.9%) in all three libraries was expected because of the existence of other mutations (except Ccdc85cMUT) in all the three libraries. To further analyze the networks, we performed Louvain graph-based clustering of the specificity networks and calculated the modularity scores of the identified communities for each of the libraries (score of zero means the communities are the same and score of one refers to a perfect separation between communities). The modularity of a graph with respect to its division into communities measures how well separated (diverse) the different nodes (clonotypes) forming the communities are from each other (see Methods). In the TILs, the TCR clonotypes that form communities/specificity groups are almost identical in frequency (42.9% for REV, and 42% and 40% for MUT1 and MUT2). However, the average modularity score of the communities/specificity groups including the most frequent (expanded) clonotypes is 0.53 for the REV, 0.73 for MUT1 and 0.77 for MUT2, indicating lower diversity of TCR clonotype in the TILs of REV than those of MUT1 and MUT2.
Using GLIPH analysis, top ten clonally expanded CD8 T cells were computationally pooled and further analyzed for gene expression patterns of their cytotoxic and effector functions. The normalized average gene expression of cytotoxicity-associated genes (Fasl or Fas ligand, Gzmb or Granzyme b, Prf1 or Perforin1, Nkg7 or Protein NKG7) as well as other effector function genes (Tbx21 or Tbet, Pcdc1 or PD1 and Ifng or Interferon gamma) were significantly higher in the clonally expanded CD8 T cells derived from MUT1 and MUT2 libraries compared to the REV library (Fasl P value < 0.001, Gzmb P value < 0.001, Prf1 P value < 0.001 Nkg7 P value < 0.001, Tbx21 P value < 0.001, Pcdc1 P value < 0.001 and Ifng P value < 0.001). Similarly, the top 10 clonally expanded CD8 T cells of MUT1 and MUT2 libraries had a significantly higher expression of early response genes which are involved in TCR engagement (Nr4a1 or NUR/77 P value < 0.001, Nr4a2 or NUR-related factor 1 P value < 0.001 and Nr4a3 or Orphan nuclear receptor TEC P value < 0.001) (Fig. 3c). Interestingly, Ifng and Nkg7 which had a similar expression pattern in the pooled CD8 T cells of all three libraries (Supplementary Fig. 5c), had significantly higher expression in the top 10 clonally expanded CD8 T cells derived from MUT1 and MUT2 libraries compared to the REV library.
Molecular modeling of Ccdc85c MUT . To gain insight into how the leucine to phenylalanine mutation leads to immunogenic epitopes, we modeled the structures of the 11-mer TYIRPFETKVK neoepitope and the 10-mer YIRPFETKVK neoepitope bound to Kd. We modeled each corresponding WT peptide as well, to assess possible changes resulting from the mutation and thus infer how the neoepitopes might differ from self. We used the same stochastic, template-based modeling procedure previously applied to murine neoepitopes 11. For the TYIRPFETKVK 11-mer, the phenylalanine at position 6 is predicted to extend up from peptide near the MHC α2 helix, increasing the amount of exposed hydrophobic surface 5% over the wild type peptide and potentially allowing the aromatic phenylalanine to interact with T cell receptors (Fig. 4a). Other than the side chain replacement, no conformational changes are predicted to occur in the peptide. For the YIRPFETKVK 10-mer, the new phenylalanine at position 5 is predicted to pack between the peptide and the a2 helix, in this case reducing exposed hydrophobic surface area (Fig. 4b). Subtle structural changes are predicted for the exposed side chains at positions 6 and 9, which could be suggestive of changes not captured by static structural modeling, such as changes in peptide flexibility that lead to altered TCR recognition 10.