Immune escape leads to neutral-like evolutionary dynamics measured by dN/dS

To determine the impact of immune evasion in the selective landscape of tumorigenesis, we obtained a catalog of 88 genes involved in the antigen presenting machinery or previously associated to immune evasion [18] (defined as “escape genes,” Additional file 2: Table S1). We classified 9896 TCGA patients from 31 different cancer subtypes into escaped (escape +) and non-escaped (escape −) cohorts based on the presence of a non-silent point mutation in one of these genes (Fig. 1). These resulted into 2089 escape + individuals with an average tumor mutation burden per patient (TMB) of 426—over 4 times higher than the average TMB for the 7087 escape − patients (95 mutations per individual, Additional file 2: Table S2). Specifically, we observed that the average number of mutations per individual in escape + was 4.51 times higher for missense and 3.95 times higher for truncating mutations compared to escape − tumors. Other mutation types, such as essential splice sites, missense and nonsense events, were also higher in escape + compared to escape − (Additional file 1: Fig. S1). When considering at least a single mutation in one of these genes, there was a heterogeneous proportion of escape patients between tumor types, i.e., TCGT had the lowest proportion of only 1% escape + patients versus SKCM that had 51% of escape + patients (Additional file 2: Table S2, and Additional file 1: Fig. S2). When performing hierarchical clustering of escape gene frequencies (Additional file 2: Table S3), we observed tumors with a similar profile of immune evasion, i.e., lung and melanoma tumors (Additional file 1: Fig. S3).

Fig. 1

9896 tumors across 31 cancer subtypes from TCGA were classified into escape + and escape − based on mutations presented in the antigen presenting machinery. These two cohorts were then analyzed to detect genes under significant selection using dN/dS. Genes under significant selection can be used as molecular targets to improve cancer treatment strategies

We then calculated dN/dS on missense and truncating mutations for the escape + and escape − cohorts. We first calculated cohort global dN/dS, using all 19,562 protein coding genes, and driver dN/dS, using 366 known driver genes. At the pancancer level, global dN/dS (Additional file 1: Fig. S4, Escape + : 1.05, CI = [1.041:1.051], Escape − :1.07, CI = [1.060:1.072], Wilcoxon-Mann p-value = 4.45e − 6), and driver dN/dS (Fig. 2A, Escape + :1.223, CI = [1.189:1.257], Escape − :1.619, CI = [1.57:1.669], Wilcoxon-Mann p-value = 0.0067) were significantly lower and closer to neutrality in escape + compared to escape − , suggesting different evolutionary trajectories for each group. To control for possible mutation burden bias, we randomized the list of 88 “escape genes,” exclude patients classified as escape + , and calculate driver dN/dS. We found that driver dN/dS of the “random escape + ” was significantly higher than the driver dN/dS of the true escape + group (Additional file 1: Fig. S5, Random escape + dN/dS ~ 1.5 vs True escape + dN/dS ~ 1.22).

Fig. 2

Selective landscape in escape + and escape − tumors. A Overall driver dN/dS for escape + and escape − tumors across all TCGA tumor types. B Overall driver dN/dS for 31 cancer subtypes separated by escape status (red—escaped − , blue—escaped +). Volcano plot for gene-level dN/dS using missense mutations versus p-value (Log10). Venn diagram showing the number of significant driver genes (Q-value < 0.1) using missense mutations considering E all genes or F restricted to known driver genes when separating by escape group or with all samples together. G List of significant driver genes using missense mutations in at least one group. H dN/dS value for significant genes in escaped − versus escaped + groups

When looking at each tumor type, global dN/dS was significantly higher than one in 25 out of 31 cancer types in escape − patients compared to 19 out of 31 tumor types in escape + patients (Additional file 1: Fig. S6). Driver dN/dS was significantly higher than one in 29/31 escape − tumors compared to 20/28 escape + tumors. Overall, global dN/dS was similar between escape − and escape + , with a majority (21 out of 31) of escape − having higher global dN/dS compared to escape + . In ACC, GBM, and UVM, global dN/dS was significantly higher in the escape − group, whereas in KIRP, global dN/dS was significantly higher in the escape + group. Driver dN/dS of escape − was higher in several cancer types (BLCA, COAD, ESCA, GBM, LIHC, LUAD, LUSC, STAD, UCEC) compared to escape + (Fig. 2B). Moreover, when including deletions and point mutations in “escape genes” to classify escape status, driver dN/dS for escape + was not significantly different from one (Additional file 1: Fig. S7). Similarly, when including overexpression of PDL1 (Additional file 1: Fig. S8), as an orthogonal evasion mechanism, driver dN/dS was closer to one for escape + compared to escape − .

Next, we hypothesized that if global and driver dN/dS are different between escape + and escape − groups, the driver gene landscape would also be different. dN/dS analysis using missense and truncating mutations in all genes revealed 85 significant driver genes across pancancer (Additional file 2: Table S2). For missense mutations, there were 30 and 68 significant genes in escape + (Fig. 2C) and escape − (Fig. 2D) tumors, respectively. Seventeen out of 30 were escape + specific and 55 out of 68 were escape − , with 13 driver genes common to both groups. For truncating events, we found 64 and 41 driver genes for escape − and escape + tumors, respectively, with 33 out of 64 escape − , 10 out of 41 escape + , and 31 common to both (Additional file 1: Fig. S9). To determine whether stratifying patients into molecular subgroups revealed novel driver genes, we calculated dN/dS using all patients together. Twenty-nine out of 68 significant genes in escape − were missed when mixing patients (Fig. 2E). Similarly, 12 driver genes were only found in the escape + group, highlighting the impact of mixing patients with different evolutionary paths into cohort analysis for cancer driver discovery. Moreover, if we restrict this analysis to only known driver genes, we still observed 6 genes under significant selection only in the escape − group, which were missed when combined with escape + patients (Fig. 2F). Importantly, when combining all patients to predict driver genes, the majority still has a significant p-value (89 out of 94 genes) but lost significance after multiple testing correction (Fig. 2G), hence the importance of properly stratifying groups.

Among 55 escape − specific genes, the majority (41/55) were evolving neutrally in the escape + group despite having a similar number of mutations (Additional file 1: Fig. S10). For escape + specific significant genes, 15/17 have small signals of positive selection in the escape − group with values closer to neutrality (Additional file 1: Fig. S11). Interestingly, two genes were under significant negative selection in the escape + group: SLC12A3 and C19orf47. Among the common genes, there was a significant higher dN/dS in 9 out of 13 genes in the escape − compared to escape + , and all were previously known drivers (Additional file 1: Fig. S12). Driver dN/dS of escape − versus escape + was significantly correlated in the pancancer analysis (Fig. 2H) and in most cancer types (Additional file 2: Table S4). Few genes acted as strong drivers in one group while being completely neutral in the other. In the escape − group, these drivers included GTF2I, VHL, REG1B, SPANXD. In the escape + group, these included CPZ, CRTAP19-5, and CFAP58. GTF2I and REG1B are associated with negative regulation of angiogenesis and antimicrobial humoral immune response. VHL is involved in cell morphogenesis and the negative regulation of apoptotic process. Allele frequencies of driver genes in pancancer (Additional file 1: Fig. S13) and per-cancer (Additional file 1: Fig. S14) were significantly higher in driver compared to escape genes, suggesting that early clonal expansions precede acquisition of evasion mechanisms.

Mutations are evenly distributed across driver genes in immune-escaped patients

To determine the immune system’s impact on the driver genes landscape, we explored whether mutations occurred more often at specific driver sites in escape + and escape − groups. We first compared two significant driver genes, IDH1 and KRAS. While IDH1 dN/dS was significantly higher in the escape − group (dN/dS ~ 31 versus dN/dS ~ 3, Fig. 3A), KRAS dN/dS was higher but not significantly different (dN/dS of 75 versus dN/dS of 35, Fig. 3B) compared to escape + . Known hotspots for IDH1 (Fig. 3C, position R132) and KRAS (Fig. 3D, position G12) were the most frequently mutated in both groups. Somatic mutations were more abundant in the hotspot of escape − compared to escape + group, while the number of unique sites mutated remained similar (Fig. 3E, IDH1: chi-square p-val = 3.6e − 12, Fig. 3F KRAS: chi-square p-val = 0.01). We then performed the same test on 55 significant “de novo” driver genes together and found the same pattern of mutations occurring preferentially at specific sites (i.e., hotspots) in escape − patients while occurring more evenly distributed across the gene in escape + patients (Fig. 3G, Pandriver p-value = 2.49e − 53). This was the case for known driver genes such as BRAF (p-value = 3.01E − 08), TP53 (p-value = 3.41E − 11), EGFR (p-value = 0.00877), and GNA11 (p-value = 0.00689).

Fig. 3

Non-random mutational distribution in non-escaped tumors. dN/dS for A IDH1 and B KRAS in escape − (red) and escape + (blue) groups including number of nonsynonymous (upper number) and synonymous (lower number) mutations. Lolliplots for mutations in escape − and escape + tumors for C IDH1 and D KRAS. Chi-square test comparing mutation number and unique mutated sites for E IDH1 and F KRAS. G Chi-2 p-value for significant known driver genes comparing the number of mutations versus the number of sites in both groups

To test whether the greater number of escape − individuals was a confounding factor, we repeated our analysis by downsampling the number of escape − patients to match the number of escape + patients. We found that fewer unique sites were mutated in the escape − group in IDH1, but not in KRAS, with respect to all mutations (Additional file 1: Fig. S15). This result suggests that there is a trade-off, at least for some driver genes, between oncogenicity and immunogenicity of the mutations accumulated in certain mutational hotspots. We then removed mutations in the most common hotspots of IDH1 and KRAS (R132 and G12, respectively). We found that for IDH1, the difference was not significant, suggesting that it is the only non-immunogenic hotspot, while for KRAS, there was still a significant difference, possibly associated to multiple non-immunogenomic hotspots (Additional file 1: Fig. S16). To determine whether oncogenes and tumor suppressors would be equally affected by escape status, we compared the proportion of mutations and found that oncogenes have more mutations in escape − patients compared to escape + (Additional file 1: Fig. S17). Interestingly, we also looked for selective differences between escape + and escape − in specific molecular subtypes. We found differences in ER + but not in ER- breast cancer patients (Additional file 1: Fig. S18), a reverse signal in HPV + compared to HPV − head and neck cancer patients (Additional file 1: Fig. S19), and overall minimal differences in HBV or HCV positive and negative patients of liver carcinoma (Additional file 1: Fig. S20).

Mutational signatures associated to immune evasion in cancer

To determine the mutational signatures associated to immune evasion, we ran deconstructSigs [22] on both groups. We found a different profile in the trinucleotide substitutions between escaped + (Fig. 4A) and escaped − (Fig. 4B) tumors, especially in sites associated to C > A substitutions in a TCT context and C > T substitutions in a TCA or TCC context.

Fig. 4

Mutational signatures associated to immune evasion. Proportion of mutation substitution in 96 trinucleotide contexts for A escaped + and B escaped − tumor cohorts. Dominant signatures per patient in C escaped + and D escaped − tumor cohorts. E Distribution of frequency of the top 60 signatures in escaped − tumors and their distribution in escaped + tumors

We tested whether immune escape allows a broader repertoire of mutational signatures to occur by comparing the most dominant signatures in escaped and non-escaped individuals. Although SBS1 signature was the most dominant in both groups, a greater proportion of escape + individuals (Fig. 4C) exhibited an alternative dominant signature compared to escape − (Fig. 4D). The top signatures for the escape + group were SBS4 (smoking), SBS13 (APOBEC), SBS7a (UV exposure), and SBS2 (APOBEC). In contrast, escape − signatures followed a flatter distribution with SBS39 (unknown) and SBS5 (clock-like unknown) being the second and third more frequent signatures. We next compared the top differential signatures between escaped + and escaped − . Interestingly, we found that immune evasion was associated to APOBEC, tobacco, and UV light signatures, while escape − tumors harbored signatures associated to mismatch repair deficiency and to various chemical exposures (Fig. 4E). We also investigated differences in mutational signatures between escape groups per cancer type (Additional file 1: Fig. S21) and found that 16 out of 31 tumor types have at least one mutational signature significantly different between groups. Escape + lung adenocarcinomas had a significant higher frequency of SBS4 (tobacco associated), compared to escape − tumors, suggesting that tobacco smoke leads to neoantigen accumulation accelerating the acquisition of immune evasion mechanisms. When controlling for mutation burden, we found that 46 out 49 mutational signatures tested display a significant difference between groups (Additional file 2: Table S5).

Immune inflammation leads to better prognosis in the absence of immune escape

To determine whether classifying tumors into immune categories can reveal a difference on the clinical prognosis between immune escaped and non-escaped tumors, we classified patients into 6 categories previously defined by Thorsson et al. [23] (Additional file 2: Table S6). We found that the only pancancer category where there was a significant overall survival difference between escape − and escape + cohort was C3 (Fig. 5, p-value = 0.0001), which was characterized by an inflammatory signature. The other immune categories displayed no significant survival advantage (Additional file 1: Fig. S22), suggesting that the inflammatory response is a major factor responsible for long-term neoantigen-mediated immune surveillance. Immune-escaped tumors from the inflammatory cluster had a lower median overall survival time compared to non-immune escaped tumors possibly associated to the absence of immune control. An observation which was recently demonstrated in long-term survivors of pancreatic cancer [24], which have stronger immunoediting signals compared to short-term survivors, characterized by weak immunoediting and high intratumoral heterogeneity.

Fig. 5

Survival and proportion comparison for escape + versus escape − patients classified as Inflammatory group (referred as C3) from Thorsson et al.

In our analysis of individual cancer types, thyroid carcinoma emerged as a prominent cancer within the inflammatory group (Fig. 5). Additionally, we found that escape + patients with bladder cancer (BLCA) and cervical squamous cell carcinoma (CESC) had improved survival rates. Conversely, escape − patients with mesothelioma (MESO), thyroid carcinoma (THCA), and thymoma (THYM) had a higher survival (Additional file 1: Fig. S23). These findings point to significant survival differences between escape + and escape − groups on a global scale; however, these patterns do not hold consistently for all cancer types. This inconsistency might indicate that our current understanding of immune escape mechanisms is not entirely comprehensive, or it could suggest that immune evasion impact depends on the tissue type.

In summary, our escape classification, in conjunction with immune infiltration measures, seems to illustrate varying stages of co-evolution between somatic and immune cells, influencing cancer prognosis differently across various organs and tissues.