Variable Mutation Expression in Human Cancers: A “Hide-and-Seek” Mechanism Linked to Differential MHC-I Presentation Dynamics

Accurately identifying tumor-specific alterations in a patient's tumor is a preferred, if not mandatory, condition to ensure an optimized selection of gene-targeted therapies by the oncologist (1–5). Still, the confirmed presence of a pathogenic variant in a tumor sample may not always lead to a better outcome for the patient treated, as demonstrated by the high variability of efficacies observed in the latest large-scale precision oncology clinical trials (6–14). While several mechanisms of primary and secondary resistance to anticancer treatments have been proposed (15), the possibility that a genomic variant is not expressed in the tumor (either at the RNA or at the protein level) and therefore does not constitute a suitable target, is often underestimated.

The democratization of high-throughput DNA- and RNA- sequencing methods now allows a better understanding of mechanisms underlying oncogenesis, and several authors recently reported a non-negligible frequency of genomic alterations that are silenced at the RNA level (e.g., up to 20% of somatic variants are not expressed at the RNA level in genitourinary tumors, and close to 18% in lung adenocarcinomas; refs. 16–18).

The concept of variable variant expressivity is not novel in human genetics. Indeed, for many inherited diseases, including familial cancer syndromes, the same mutation is not always expressed in all individuals who carry it (concept of genetic penetrance; refs. 19, 20); moreover, when the mutation is expressed, it is not always expressed in the same way (concept of genetic expressivity; ref. 20), ultimately leading to high phenotypic unpredictability between patients. The same observations may hold true when considering the genetics of sporadic tumors.

While somatic variant silencing could explain some of the differences observed in terms of cancer diagnosis, prognosis and response to therapy, it may also be linked to differential immune responses. Immunosurveillance is the process by which cells of the immune system seek out, recognize and suppress foreign pathogens, as well as precancerous and cancerous cells in the body (21, 22). Tumor antigen–specific cytotoxic T lymphocytes (CTL) are the major effectors of the immune response against tumor cells, and CD8 (+) T lymphocytes are specifically able to recognize tumor-associated antigens and tumor neopeptides (peptides originating from mutated proteins) presented by the cell through the MHC-I moieties (23). Unfortunately, T-cell pressure can sculpt the antigenicity of tumor cells, resulting in the emergence of tumor clones that lack defined mutant antigens and therefore escape immunosurveillance; a phenomenon referred to as immunoediting and loss of immunogenicity (21, 22, 24).

In this study, our aim was to describe the prevalence of somatic mutation loss of expression (“variant silencing”) in a large collection of human samples, and the potential impact of such variant silencing on immunogenicity. We, therefore, analyzed the association between variant silencing and tumor characteristics and/or immune environment settings, before conducting an extensive MHC-I presentation evaluation (from endogenous protein processing to T-cell recognition features) of the corresponding native and mutated peptidomes. We observed that silenced variants are those generating peptides that are more likely to be positively processed or presented by the MHC-I machinery, or later recognized by the T-cell machinery. Hence, variant silencing is likely a mechanism that results in immune evasion or, alternatively stated, tumors that flourish may do so because silencing of immunogenic peptides has occurred (cancer selection by immune pressure).

Data originated from The Cancer Genome Atlas (TCGA) project (https://www.cancer.gov/about-nci/organization/ccg/research/structural-genomics/tcga). Only samples for which RNA expression was available in both tumor and normal tissue were included. Whole-exome mutation data were available for 9,104 tumor samples. Sequencing-based mRNA expression was available for 11,060 samples including 737 normal adjacent tissues. Samples corresponding to metastatic and recurrent lesions were excluded, leaving only primary tumors. Only samples with complete mutation description, mRNA primary tumor plus paired-normal tissue expression profiles, and MHC-I haplotype were kept for analysis. A total of 636 unique patients diagnosed with 21 distinct tumor types (all solid tumors) remained after the above filtering.

This study respected TCGA's Human Subjects Protection and Data Access Policies (https://cancergenome.nih.gov/abouttcga/policies/tcga-human-subjects-data-policies).

Somatic mutations (single nucleotide variants and small insertions/deletions), gene-level mRNA expression data and sample immune subtypes (as described in ref. 25) were downloaded from the University of California Santa Cruz (UCSC) Xena portal (https://xena.ucsc.edu/). Expression counts were previously normalized (by UCSC) in order to avoid potential batch effect and discrepancies between tumor types.

MHC-I haplotypes, corresponding to the description of 2 HLA-A, 2 HLA-B, and 2 HLA-C 4-digit allelic genotypes (i.e., a total of 6 alleles per sample) available for 8,912 tumor samples, were performed by (25) and downloaded from the TCGA data portal.

Silent mutations (N = 24,178), complex mutations (splice sites, in-frame or frame-shift deletions/insertions, intron and large deletion variants - N = 10,759), and variants intrinsically able to dysregulate mRNA and/or protein expression (i.e., 3′-UTR/5′-UTR, nonsense, nonstop, translation-start site and RNA altering variants (meaning alterations that change the RNA splicing process – N = 17,096) were excluded from the set of genomic variants initially described (N = 117,105). The last filter was applied to avoid potential confounding factors related to well-known silencing mechanisms such as nonsense-mediated decay (NMD) or nonstop decay (NSD). In fine, a total of 65,072 missense mutations were kept for the MHC-I processing and presentation analysis.

The mutated protein sequences produced by each genomic variant were retrieved from the Cancer Mutant Proteome Database (in silico protein translation, http://cgbc.cgu.edu.tw/cmpd/; ref. 26).

Full protein sequences and MHC-I haplotypes were then submitted to the epitope prediction tools available from the Immune Epitope Database (IEDB) and accessible by the dedicated web API (http://tools.iedb.org/main/tools-api/). Processing and MHC-I binding scores were computed using the NetMHCpan v.3.0 algorithm on 9-mer peptides. For each variant/protein, the maximum constitutive and IFNγ-induced proteasome cleavage scores (corresponding to the logarithms of the total amount of cleavage sites leading to the liberation of the C-terminus peptide), maximum TAP score (corresponding to the negative logarithm of the IC₅₀ value for the binding to the TAP transporter), minimum MHC-I binding IC₅₀ (corresponding to the concentration of peptide required to displace 50% of a probe from the respective MHC-I allele), as well as the total number of high to intermediate MHC-I binding affinity peptides (considering a MHC-I binding IC₅₀ threshold of 500 nmol/L) were evaluated.

The complete set of 9-mer peptides was then submitted to the IEDB T-cell immunogenicity predictor (http://tools.iedb.org/immunogenicity/), using the default settings (first, second, and C-terminus amino acids masked for prediction). For each variant/protein, the maximum immunogenicity score (corresponding to the probability for each peptide to elicit an immune response) was kept for evaluation.

Combined scores estimating the MHC-I presentation process in its entirety were also computed: maximum (peptide processing) score (corresponding to the association of both proteasome cleavage and TAP transport scores), maximum (processing followed by MHC-I binding) score, as well as maximum (processing, MHC-I binding followed by T-cell recognition) score.

Finally, an hydrophobicity index corresponding to the sum of all amino acid hydropathy scores (as defined in ref. 27), as well as the number of hydrophobic amino acids per peptide (namely isoleucine, valine, leucine, phenylalanine, cysteine, methionine and alanine; as defined in ref. 27) were considered for all possible 9-mer peptides. For each variant/protein, the maximum hydrophobicity index and the maximum number of hydrophobic residues per peptide were calculated.

A specific computational workflow allowing the automatization of all steps was created using the KNIME Analytics Platform (https://www.knime.org/knime-analytics-platform) and is available upon request. KNIME Analytics Platform is an open-source software allowing the design of data science workflows, using preprogrammed analytic nodes and a graphical user interface facilitating the management of data flow by a drag-and-drop system. KNIME Analytics Platform is available at https://www.knime.com/.

A genomic variant was considered silenced when the mRNA of the corresponding gene showed a −80% decreased expression in the tumor in comparison to the normal adjacent tissue. All other variants were grouped as nonsilenced variants. The distribution of silencing variants between different types of genomic variants (mutation subtypes, nucleotide substitution subtypes, and genes impacted by the variants) was studied using a χ² test. The silencing rate of each sample (defined as the percentage of silenced variants over all described variants in one tumor sample) were compared between tumor types, immune subtypes (as previously defined in ref. 25) and tumor mutation burden (“TMB high” corresponding to tumors in the 75th–100th percentile in term of total number of mutations, “TMB intermediate” corresponding to tumors in the 50th–75th percentile, and “TMB low” corresponding to tumors in the 1th–50th) categories using a one-way (single factor) ANOVA.

MHC-I presentation scores and hydrophobicity features were compared between silenced and nonsilenced variants using a bidirectional Student t test, assuming nonequal variances.

For all tests, P values ≤ 0.05 were considered significant. Statistical analyses were carried out using SAS University Edition (https://www.sas.com/en_us/software/university-edition/download-software.html).

Data originated from The Cancer Genome Atlas (TCGA) project and are available online at https://www.cancer.gov/about-nci/organization/ccg/research/structural-genomics/tcga.

The overall analytic workflow is described in Supplementary Fig. S1.

Several one-way ANOVA were conducted to estimate the effect of tumor histologic and immunologic specificities on the mutation silencing rates observed within each of the samples of interest. In this subset of 636 samples encompassing 21 distinct histology diagnoses, the average percentage of variants silenced per sample at the mRNA level was 7.84% [95% confidence interval (95% CI), 7.44–8.24)].

There was a significant effect of the tumor type on the silencing rate observed in each sample for the 21 conditions (noncorrected P = 2.5 × 10^–11; Bonferroni-corrected P value considering multiple comparisons = 5.3 × 10^–9, the number of comparisons being | $C_2^{21}\ = \ 210$ |⁠). The tumor types presenting significantly higher silencing rates (in comparison to the entire subset of tumors and considering tumor types with reasonable sample size) were kidney chromophobe carcinomas [N = 25 samples; mean (95% CI) = 12.93 (10.28–15.58) silenced variants per sample] and lung squamous cell carcinomas [N = 45; 10.13 (8.93–11.32)]. The tumor types presenting significantly lower silencing rates were thyroid [N = 57; 4.81 (3.03–6.60)] and prostate carcinomas [N = 52; 5.71 (4.43–7.00); Fig. 1].

Interestingly, the silencing rates per sample were also significantly associated with the type of immune infiltration (non-corrected p-value = 1.5 × 10^–5; Bonferroni-corrected p-value considering multiple comparisons = 3.2 × 10^–4, the number of comparisons being | $C_2^7\ = \ 21$ |⁠), as well as the overall tumor mutation burden (TMB) (a molecular feature associated with distinct immune responses) (non-corrected p-value = 6.8 × 10^–6; Bonferroni-corrected p-value considering multiple comparisons = 2.0 × 10^–5, the number of comparisons being | $C_2^3\ = \ 3)$ |⁠. Indeed, the higher silencing rates were observed in lymphocyte-depleted tumors (N = 58; 10.44 [9.13–11.75] %) and samples presenting a high TMB (≥75^th percentile in term of total number of mutations; N = 160; 9.24 [8.58–9.80] %); while the lowest silencing rates were observed in inflammatory tumors (N = 247; 6.64 [5.95–7.33] %) and samples presenting a low TMB (<50^th percentile in term of total number of mutations; N = 317; 6.90 [6.25–7.35] %) (Fig. 1).

Several χ² analyses were conducted to compare distributions of silencing rates between variant subtypes and gene localization. The entire set of variants considered (originating from 636 different samples) presented an overall silencing rate of 9.06% (N = 10,604 silenced/117,505 total variants). As previously mentioned, the average silencing rate per sample was 7.84% (Fig. 1).

Surprisingly, missense mutations, a category of variants not usually described as leading to intracellular mRNA silencing and/or decay mechanisms, presented a silencing rate significantly higher than non-missense mutations (9.30% vs. 8.75%, P = 0.0013). Nonstop mutations and RNA altering mutations presented higher silencing rates, albeit nonsignificant when compared with all other mutation types. The lack of significance of such categories may be possibly due to a low number of variants observed. In-frame deletions, frameshift deletions, and 3′-untranslated region (UTR) mutations presented significantly lower silencing rates than all other mutation types (4.63%, 7.89%, and 7.99%, respectively; all P values < 0.0100; Supplementary Figs. S2 and S3A).

The type of nucleotide substitution implicated for each missense mutation also had an impact on the variant silencing rates. Indeed, nucleotide transversions (a change from a purine to a pyrimidine nucleotide or vice versa) presented significantly higher silencing rates than nucleotide transitions (a change from a purine to another purine or from a pyrimidine to another pyrimidine; 9.84% vs. 8.77%; P = 2.5 × 10^–6; Supplementary Fig. S2). Transversions G>T and C>A were associated with the highest silencing rates (11.25% and 10.87%, respectively) when compared with all other substitution types. The detailed silencing rates for specific nucleotide changes are provided in Supplementary Fig. S3B.

We further explored the effect of variant localization on the silencing rates, notably when these variants arise in well-known protein-coding oncogenes or tumor suppressor genes. In this cohort of treatment-naïve primary tumor samples, NTRK1/2/3, FLT3, and ERBB4 variants were the most likely to be silenced (3/8, 4/5, 4/12, 3/5, and 8/19 variants, respectively); while silencing of frequent variants localized within ubiquitous oncogenes such as PIK3CA, BRAF, and the tumor suppressor TP53 was never observed (0/81, 0/45, and 0/150 variants; respectively; Fig. 2). High silencing rates were also observed in other frequently mutated genes, such as APOB, SPTA1, LRP1B, and LRP2 (Supplementary Fig. S2C).

Further MHC-I peptide processing investigation was performed for all 11,986 genomic variants for which protein sequence prediction was available (5,589 silenced and 6,397 nonsilenced variants). For each virtual protein, the entire set of possible 9-mer peptides was analyzed, and the best scores obtained for all individual processes (proteasome cleavage and transport, MHC-I binding and T-cell recognition – calculations were performed using the T-cell epitope prediction tools available from http://tools.iedb.org/main/tcell/) were retained.

Peptides generated by mutated proteins undergoing tumor-specific silencing were found to be more hydrophobic than those not silenced (hydrophobicity index = 25.4 vs. 23.6, respectively; P = 1.6 × 10^–6). These peptides were more likely to encompass amino acids considered hydrophobic by the Kyte-Doolittle scale (ref. 27; up to 7.8 vs. 7.6 hydrophobic residues per peptide; P = 5.0 × 10^–38; Table 1; Fig. 3). This physicochemical feature is important because hydrophobicity has been previously considered a hallmark of both MHC-I antigen processing and T-cell recognition (28–30).

Indeed, silenced proteins were found to generate more peptides able to bind to MHC-I than nonsilenced proteins (93 vs. 83 MHC-I–binding antigens per protein when considering an IC₅₀ lower than 500 nmol/L as a criterion for MHC-I binding; P = 1.5 × 10^–6). However, the silenced peptides were not associated with a better overall higher affinity to MHC-I moieties as compared with the nonsilenced peptides [IC₅₀ (95% CI) = 11.4 (10.2–12.5) vs. 11.6 (10.8–12.5) nmol/L; P = 0.7090], denoting a quantitative difference (rather than a qualitative difference) in MHC-I–binding capacities between silenced and nonsilenced peptides. Interestingly, the silenced peptides presented significantly higher proteasome processing (i.e., the intracytoplasmic process by which a protein is cleaved into smaller peptides) scores (P = 0.0181 in a constitutive environment and P = 0.0028 in an IFNγ-induced environment), a more effective endoplasmic reticulum transport (P = 1.5 × 10^–5), as well as better T-cell recognition abilities (P = 8.5 × 10^–6; Table 1; Fig. 3).

Assuming that a neoantigen requires optimized scores for all of the steps above (and not only for one step of the process) to be ultimately recognized by the host immune system, we also measured several combined scores. These combined scores correspond to the association of the proteasome processing score, followed by optimized reticulum transport score, followed by optimized MHC-I–binding score, and finally followed by an optimized T-cell recognition score. As expected, peptides generated by silenced proteins presented better overall combined MHC-I presentation pathway features than those originating from nonsilenced variants; this was true in both a constitutive (noninflamed) environment (P = 2.8 × 10^–7) and in an IFNγ-induced environment (P = 1.8 × 10^–8; Table 1).

Finally, we analyzed the differences in silencing rates of neo-epitopes among tumors presenting a silencing of HLA-A, HLA-B, and HLA-C moieties (evaluated at the mRNA level and considering a decrease of expression ≥ 80% comparatively to the normal tissue – N = 629 tumors evaluated) and found that the silencing rate was significantly higher in tumors presenting a lower expression of HLA-C alleles (12% vs. 8%; P = 0.037). The difference was not significant for HLA-A and HLA-B alleles, possibly due to a lower number of tumors presenting a decrease of expression for these mRNAs (N = 9 and 3, respectively; Supplementary Fig. S4).

In this study, we demonstrated that a proportion of somatic variants are not expressed at the tumoral level, in comparison to the normal adjacent tissue (more than 9% of all variants). Interestingly, missense variants, a category of variants not associated with canonical RNA decay mechanisms, presented a silencing rate significantly higher than all other mutation types (9.30% vs. 8.75%, P = 0.0013). This first observation may be of importance when considering the soaring development of anticancer treatments targeting selective tumor alterations (2–6, 31). Indeed, the recent progress in precision medicine remains limited by efficacy in only a subgroup of patients, sometimes explained by tumor biology complexity and/or heterogeneity (6–14, 31). But one may also consider the effective presence of the alteration at the protein level to be a necessary component of drug effectiveness.

Missense silent variants were significantly related to proteins presenting a better intracellular peptide processing (P = 0.0028 for the score corresponding to the IFNγ-induced inflammatory proteasome cleavage), better transport to the endoplasmic reticulum [ P < 0.0001 for the transporter associated with antigen processing (TAP) transport score], subsequent higher number of peptides able to bind with a high affinity to MHC-I moieties (P < 0.0001), and better recognition by lymphocyte T cells (P < 0.0001). In addition, these silenced proteins were processed to peptides significantly more hydrophobic than the peptides generated by nonsilenced proteins; the peptide hydrophobicity being a feature of immune recognition already described by other authors (29, 30, 32). Overall, we demonstrated that the mRNA silencing of tumor-specific genomic missense mutations is highly correlated with the ability of the mutated protein to generate hydrophobic epitopes (with hydrophobicity generally correlating with immunogenicity; refs. 28, 29).

The balance between immunodominant (i.e., presented to the immune system) and cryptic (i.e., hidden from the immune system) peptides is a subtle cell mechanism that can be disturbed, consequently resulting in undesirable autoimmunity issues (when cryptic peptides are too well presented) or, on the contrary, to infection and/or tumorigenesis events (when immunodominant peptides are no longer presented; ref. 33). MHC-I–bound immunodominant peptides are required to generate a strong immune response against cells invaded by intracellular pathogens such as viruses. But those viral peptides can be mutated, thus allowing evasion from the mandatory adaptive immune response and later active disease perpetuation, as has been demonstrated with the human immunodeficiency virus (HIV; refs. 34, 35). Similarly, peptide immunodominance disequilibrium can also have implications in cancer development. As the immunodominant epitopes are mutated and/or lost in the tumor, the immune response is no longer potent enough to battle the undesirable cells (36). Cancer immunodominant epitope losses may lead, over time, to the natural selection of tumor clones escaping the immune response and favoring negative outcome of the disease (36–38). Several immune-editing and selection pressure mechanisms have been recently described in early-stage untreated non–small cell lung cancer, such as copy-number loss of previously clonal neoantigens, loss of heterozygosity in human leukocyte antigens (HLA), or even promoter hypermethylation of genes that contain neo mutations (39). Further analyzes will be needed to fully elucidate these processes.

Here, we showed that variant silencing was more likely associated with lymphocyte-depleted tumors, possibly supporting the hypothesis of epitope loss as an immune evasion mechanism. In addition, it must be noted that immunodominance is an individual characteristic: people express specific and almost unique MHC-I haplotypes; and patients may, therefore, present different epitopes of the same tumor protein, rendering the individual prediction of immune surveillance efficacy and tumor evolution even more complex.

Differential mRNA expression may be explained by canonic tissue specificities (e.g., NTRK1/2/3 receptors are mostly expressed in brain, adrenal, thyroid tissues, and FLT3, ERBB4, ALK, and KIT are also highly expressed in brain tissues; ref. 40). However, the mechanisms underlying the silencing of mRNA and proteins carrying oncogenic mutations in a tumor cell in comparison with the normal adjacent tissue are still unknown. One may suspect the implication of tumor differentiation processes, differential translation efficiencies due to oncogene codon-usage bias (41), as well as an immune-mediated mRNA expression regulation for cells presenting specific neoepitopes (immunoediting process; ref. 42). Thorough in vitro analyses of neoepitope immunogenicity in different tumor types are needed to elucidate the relationship between tumor microenvironment and somatic mutation silencing process. Another limitation of the article is that mRNA expression data in tumor tissues did not specifically discriminate between wild-type and mutated alleles, and therefore did not acknowledge for possible allelic imbalances. However, the normal tissue harbored wild-type alleles. Finally, functional testing is needed to show that specific hydrophobic mutants are immunogenic, although it has previously been demonstrated that, in general, hydrophobicity correlates with immunogenicity (28, 29).

Finally, we are suggesting that the silencing of somatic events may enable potential resistance to anticancer therapies when those treatments aim at targeting specific mutated proteins that are no longer expressed in the cell. We demonstrated that NTRK1/2/3, FLT3, ALK, KIT, and ERBB4 mutations – all actionable by FDA-approved drugs – are often silenced in cancer cells. This mechanism is, as for now, ignored by oncologists practicing personalized medicine. In conclusion, this study suggests that, in the future, therapeutic choices should be molecularly informed by both the presence of the mutation at the genomic level and by its expression at the mRNA and/or protein level.

In this study based on a large cohort of pan-cancer tumors, we demonstrated that more than 9% of all genomic mutations are not expressed at the mRNA level in human cancer samples. Such variable expressivity may explain differences observed in cancer development, prognosis, and tumor immune settings. Indeed, silenced variants were significantly related to proteins presenting better peptide processing, MHC-I binding, and T-cell recognition (all antigenic and immunogenic features being extensively studied herein). Not surprisingly, this loss of immunogenicity was also more likely observed in lymphocyte-depleted tumors. Hence, we are suggesting that somatic variant silencing may participate in tumor resistance by clonal selection and immune evasion, (as well as primary resistance to molecularly targeted therapy).

R. Kurzrock reports grants, personal fees, non-financial support, and other support from R. Kurzrock has received research funding from Biological Dynamics, Boehringer Ingelheim, Debiopharm, Foundation Medicine, Genentech, Grifols, Guardant, Incyte, Konica Minolta, Medimmune, Merck Serono, Omniseq, Pfizer, Sequenom, Takeda, and TopAlliance; as well as consultant and/or speaker fees and/or advisory board for Actuate Therapeutics, AstraZeneca, Bicara Therapeutics, Biological Dynamics, Daiichi Sankyo, Inc., EISAI, EOM Pharmaceuticals, Iylon, Merck, NeoGenomics, Neomed, Pfizer, Prosperdtx, Roche, TD2/Volastra, Turning Point Therapeutics, X-Biotech; has an equity interest in CureMatch Inc., CureMetrix, and IDbyDNA; serves on the Board of CureMatch and CureMetrix,and is a co-founder of CureMatch. outside the submitted work; and R. Kurzrock has received research funding from Biological Dynamics, Boehringer Ingelheim, Debiopharm, Foundation Medicine, Genentech, Grifols, Guardant, Incyte, Konica Minolta, Medimmune, Merck Serono, Omniseq, Pfizer, Sequenom, Takeda, and TopAlliance; as well as consultant and/or speaker fees and/or advisory board for Actuate Therapeutics, AstraZeneca, Bicara Therapeutics, Biological Dynamics, Daiichi Sankyo, Inc., EISAI, EOM Pharmaceuticals, Iylon, Merck, NeoGenomics, Neomed, Pfizer, Prosperdtx, Roche, TD2/Volastra, Turning Point Therapeutics, X-Biotech; has an equity interest in CureMatch Inc., CureMetrix, and IDbyDNA; serves on the Board of CureMatch and CureMetrix,and is a co-founder of CureMatch. No disclosures were reported by the other authors.

A. Boichard: Conceptualization, data curation, formal analysis, validation, investigation, visualization, methodology, writing–original draft, writing–review and editing. R. Kurzrock: Conceptualization, resources, formal analysis, supervision, funding acquisition, validation, methodology, project administration, writing–review and editing.

R. Kurzrock was funded in part by the Joan and Irwin Jacobs Fund and NIH P30 CA023100.

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.