The Association of MUC16 Mutation with Tumor Mutation Burden and Its Prognostic Implications in Cutaneous Melanoma

MUC16 (OMIM 606154, previously known as CA125) is frequently mutated in many cancers. Although the functional role of MUC16 is extensively studied in ovarian cancer, its role as a clinical biomarker remains controversial. In melanoma, MUC16 is one of the top mutated genes according to The Cancer Genome Atlas (TCGA; Supplementary Fig. S1). However, this gene was disregarded in most genomic analyses because it is less likely to be a pathogenic or cancer driver gene. On the basis of the TCGA gastric cancer data and an independent validation cohort, a recent analysis showed MUC16 mutations are significantly associated with patient overall survival (1). The researchers also explored the potential function of MUC16 in influencing the immune system in gastric cancer.

The functional role of MUC16 in cancer is yet to be elucidated due to several challenges such as its complicated structure and lack of specific antibodies (2). However, one hypothesis is that MUC16 mutations might function through neoantigen presentation and immunoediting (3). Another plausible hypothesis is that MUC16, being large and heavily glycosylated (4), prevents the establishment of a robust immunologic synapse between T cells and MHC class I presentation of tumor antigen on the cell surface. The two main questions that are more relevant to melanoma in the clinical setting are (i) how strong are the associations between the mutation status of MUC16 and tumor mutation burden (TMB) and patient prognosis in melanoma—reflecting its potential as a surrogate biomarker for use in clinical practice; and (ii) whether the mutation or neoantigen of MUC16 is a predictive biomarker for patients with melanoma treated with immunotherapies. To investigate these questions, we analyzed MUC16 mutations and their associations with TMB and prognosis in patients from three independent melanoma cohorts. Furthermore, we investigated how MUC16 mutation is linked with neoantigen presentation and therapy response in patients with melanoma treated with anti–CTLA-4 and anti–PD-1 therapies.

Genomic and clinical data were acquired from the TCGA, Yale, and Moffitt melanoma cohorts (Their clinical characteristics are shown in Supplementary Table S1). Somatic mutation data from whole-exome sequencing (WES) and matched gene expression data from RNA sequencing (RNA-seq) from the TCGA melanoma cohort (SKCM, n = 467 patients) were downloaded from NCI GDC Data Portal (https://portal.gdc.cancer.gov). For the Yale melanoma cohort, the somatic mutation and clinical data were obtained from a previous study (5), which contained data from 132 patients. The Moffitt cohort contains 135 patients with stage 3 and 4 melanoma only, WES and RNA-seq data from the Moffitt cohort were included in this study and their deidentified clinical covariates of research interest were obtained from the honest broker for the Moffitt Cutaneous Oncology Program. Tumor-specific mutations were identified with Strelka (6) and MuTect (7), and were annotated to determine genic context using ANNOVAR and summarized using VarSifter (8). To ensure the same sensitivity in calling wild-type and mutants, mutation count in each gene in Moffitt data was adjusted by benchmarking with TCGA data on matched genes. As part of a larger study on computational drug repurposing (MCC# 19147), the study was conducted in accordance with recognized ethical guidelines (e.g., Declaration of Helsinki, CIOMS, Belmont Report, U.S. Common Rule) and was approved by Advarra Institutional Review Board. A waiver of consent was granted by Advarra.

Tumor mutation burden is defined as log₂ transformed total number of nonsynonymous mutations per megabase. We used 38 Mb as the estimate of the exome size in TCGA (9). Mutation data from germline samples in TCGA were downloaded from the Genome Data Commons legacy archive (https://portal.gdc.cancer.gov/legacy-archive/). Variant calling from sequence data of germline samples in the Moffitt cohort was based on GATK (10). Somatic mutational signature deconvolution (11) was performed using the method implemented in the R package deconstructSigs (12), which was able to determine the compositions of predefined mutational signatures within a single sample. We used 30 mutational signatures described in the COSMIC database (http://cancer.sanger.ac.uk/cosmic/signatures). Mutational signature analysis confirmed that UV signatures were dominant in melanomas in all studies. RNA-sequencing data from the same tumor and survival outcomes are also available.

HLA types and total neoantigen load for each sample from patients treated with immunotherapy were obtained from previously published studies (13–16). Putative MUC16-specific neoantigens were predicted for each patient by defining all novel amino acids—all nonsynonymous somatic mutations in MUC16 identified from the WES analysis were translated into 9-mer or 10-mer peptides flanking the mutant amino acid. The binding affinity of mutant neoepitopes to the patient-specific HLA-matched class I types was then predicted using NetMHCpan 4.0 (17). The thresholds for strong and weak binder were based on the default setting (0.5% rank for a strong binder and 2% for a weak binder).

We used the nonparametric Wilcoxon rank sum test to compare the TMB between MUC16 wild-type and mutant groups, as well as the differences in neoantigen burden. A two-sided P < 0.05 was considered statistically significant. Multivariable logistic regression was used to assess the association between MUC16 mutations and TMB while adjusting for confounding factors including sex, age, stage, and mutational signatures. In the logistic regression, TMB measures were dichotomized by median in the logistic regression. ORs and 95% confidence intervals (CI) of included covariates were reported. Kaplan–Meier curves and the log-rank test were used to assess the survival differences between patient groups with differing mutation status. Multivariable Cox regression analyses were used to assess the association of MUC16 mutations and prognosis by controlling for potential confounding factors. The association between patient survival and germline variants in the MUC16 region was tested using both single-variable Cox regression and multivariable Cox regression analysis with backward selection. For this study, we also replicated a recently published trivariate regression analysis by Lee and Ruppin (18) in predicting objective response rate (ORR) to anti-PD-1/PD-L1 therapy across 21 cancer types. We used the same dataset and followed the same statistical procedure provided in the original publication but with the total mutation burden in the model replaced by the MUC16 single-gene mutation burden. The other two variables included in the regression model were the CD8⁺ T-cell abundance (eCD8T) and fraction of high PD-1 mRNA expression (fPD1).

Of the 467 patients in the TCGA SKCM cohort, 312 (66.81%) patients carried somatic mutations in MUC16 (Supplementary Fig. S1), where BRAF mutations occurring at a frequency of 50%. TCGA SKCM samples with MUC16 mutations had a statistically significantly higher TMB than the samples of wild-type (P < 0.0001, Wilcoxon rank sum test; Fig. 1A). The difference in the median of the two groups is 2.44 at log₂ scale, which is larger than the difference reported in patients with gastric cancer (1). A significantly higher mutation rate was also observed in the independent Moffitt melanoma samples with MUC16 mutations (P < 0.0001, Wilcoxon test; Fig. 1B), as well as in a different melanoma dataset from Yale (P < 0.0001, Wilcoxon test; Fig. 1C). A consistent finding across three studies was that the median TMB (log₂ scale) in the MUC16-mutant group was around 9 and the median TMB in the MUC16 wild-type group was around 7. Of the tumors with lower mutation burdens (log₂TMB<5), most samples are MUC16 wild-type: 84.2% (101 of 120) in TCGA, 96.2% (25 of 26) in Moffitt data, and 70.0% (21 of 30) in the Yale data. On the basis of the distribution of somatic mutations in the MUC16 region in all three cohorts, no hotspot region was identified (Supplementary Fig. S2). As expected, patients with mutations in both of the main regions of MUC16 (N-terminal domain and tandem repeats) showed highest TMB in the entire cohorts. (Supplementary Fig. S3).

To adjust for possible confounding factors, we included covariates including age, sex, stage, POLE status, and common mutational signatures in melanoma (Supplementary Figs. S4–S6) in the logistic regression model. In mutational analysis, as expected, the signatures attributed to UV exposure (Signature 7) and alkylating agents (Signature 11) dominated the mutational signatures in the melanoma mutational landscape. Associations between MUC16 mutations and TMB remained statistically significant in the TCGA cohort [OR, 9.28 (95% confidence interval (CI), 5.18–17.39); P < 0.001; Fig. 2A], as well as in the Moffitt melanoma cohort [OR, 31.95 (95% CI, 8.71–163.90); P < 0.001; Fig. 2C] and Yale cohort [OR, 8.09 (95% CI, 3.12–23.79); P < 0.01; Supplementary Fig. S7A]. Beyond melanoma, the strong association between MUC16 and TMB was observed in multiple cancer types in the TCGA dataset (Supplementary Fig. S8).

In the Kaplan–Meier survival analysis of the TCGA cutaneous melanoma cohort (Fig. 1D), a significant difference was found in overall survival between the MUC16-mutant and MUC16-WT patient groups (log-rank test, P = 0.0013). MUC16-mutant patients have a median overall survival of 103 months (95% CI, 74.7–133 months) compared with 49.5 months (95% CI, 46.9–72 months) in the MUC16-WT patient group. Interestingly, as shown in Fig. 3A, the survival difference remained significant (log-rank test, P = 0.0096) in patient subgroup of lower than median mutation burden calculated on the basis of genes included in Illumina's TruSight Tumor 170 (TST170), a panel widely used in clinical setting. As shown in the Cox regression (Fig. 2B), the association between MUC16 and overall survival remained statistically significant after adjusting for potential confounding factors [HR, 0.62; (95% CI, 0.45–0.85); P < 0.01].

In the Kaplan–Meier survival analysis of the Moffitt melanoma cohort (Fig. 1E), MUC16 mutations were found associated with a better overall survival with a borderline statistical significance (log-rank test, P = 0.081). After controlling for covariates including age, sex, stage, POLE mutations, and mutational signatures (7 and 11) in the Cox regression, the association between MUC16 and overall survival was found statistically significant [HR, 0.49; (95% CI, 0.28–0.87); P = 0.014; Fig. 2D]. Mutational signature 1 was not included in the Cox regression because very few patients had a nonzero score in this signature. Because of the insufficient reporting of survival outcomes, no association was found between MUC16 and the prognosis defined in the Yale data (Supplementary Fig. S7; ref. 5). The association between MUC16 mutation and overall survival was observed in four TCGA cancer types (Supplementary Fig. S9).

In addition to the mutation status, we observed that the variant allele frequency (VAF) value (also known as mutant allele burden) of MUC16 was also associated with prognosis in both Moffitt and TCGA cohorts (Fig. 3B–C). Both data revealed a previously unreported but consistent phenomena: higher nonzero VAF levels of MUC16 were associated with worse survival in patients with cutaneous melanoma. We also investigated the mutation status and VAF in other top mutated genes (including TTN and BRAF) in melanoma. Surprisingly, MUC16 is the only gene among the top 10 most commonly mutated genes in melanoma for which a clear prognostic significance was observed in both mutation status and VAF (Fig. 3D–E and Supplementary Fig. S10), indicating that MUC16 is a selective clinical biomarker.

To investigate the predictive value of MUC16 mutations in immunotherapeutic studies, we next examined two sets of patients with melanoma treated with anti–CTLA-4 (13, 14) and two sets of patients with melanoma treated with anti–PD-1 therapy (15, 16), in which neoantigen load and patient HLA class data were available. Predicted neoantigens in the MUC16 protein-coding region were identified in 94 of the 174 patients in the anti–CTLA-4 cohorts, as well as in 61 of the 103 patients in the anti–PD-1 cohorts. The numbers of MUC16-specific neoepitopes were found significantly associated with total neoantigen burden across four studies (Wilcoxon rank sum test, P < 0.001; Fig. 4A). Associations between MUC16 mutations and TMB were also statistically significant in all four cohorts (Supplementary Fig. S11). Of the patients with MUC16-specific neoantigens, the median proportion of neoantigens attributed to MUC16 in these four studies are: 0.20% (Snyder and colleagues), 0.21% (Van Allen and colleagues), 1.77% (Hugo and colleagues), and 0.22% (Riaz and colleagues). Because only patients in Snyder and colleagues reported association of mutation and neoantigen load with therapy response, we further examined the association between MUC16-specific neoepitopes and therapy outcomes. As shown in Fig. 4B, in the combined patient sets, the long-term benefit patient group has a higher proportion of patients with MUC16-specific neoantigens (18 of 27, 66.67%) than the minimal or no benefit patient group (15 of 37, 40.54%). The discrepancy in such proportion was found larger in the discovery set defined in the original publication. However, the enrichment of patients with neoantigens was not found in patients responsive to anti–PD-1 (Supplementary Fig. S12). In terms of mutation status, we observed that MUC16 was associated with overall survival in combined anti–CTLA-4–treated patients with a borderline significance but not in combined anti–PD-1–treated patients (Supplementary Fig. S13). This result is expected as TMB itself has not be established as a successfully biomarker for anti–PD-1 response in melanoma.

The result from the replicated analysis of the trivariate regression model (18) for predicting the ORR to anti–PD-1/PD-L1 therapy across 21 cancer types is presented in Fig. 4C. The combined predictive model of MUC16 mutation burden, eCD8T and fPD1 levels achieved a comparably high prediction accuracy (R = 0.87, P < 4.1 × 10⁻⁷) compared with the original analysis where total mutation burden was used. Notably, the response rate of SKCM was found closer to the predicted value from the MUC16-based trivariate model [comparing with Fig. 2B in Lee and Ruppin (18)], providing additional evidence that MUC16 mutation is a promising candidate when designing a minimalist biomarker panel for melanoma.

The functional roles of germline SNPs in the MUC16 region remain underinvestigated. Currently, only 9 variants in MUC16 are annotated in ClinVar and most of them are classified as benign. Through multivariable Cox regression and controlling for clinical covariates, we found 34 germline SNPs in the MUC16 region that were associated with overall survival (Supplementary Table S2) in the Moffitt melanoma cohort. Three common SNPs (rs1609459, rs10402812, and rs7245949) that are in high linkage disequilibrium (LD) were associated with overall survival in the Moffitt data (Cox regression, P = 0.01) and were found associated with survival with a borderline significance in TCGA cohort (P = 0.07) while other SNPs showed a similar trend of effects (Supplementary Figs. S14 and S15). Genetic ancestries were not considered in this explorative analysis because of the limited sample size and the fact that the majority of patients with melanoma were Caucasians.

Because of its large genome size and as yet understudied gene function, MUC16 is an often-overlooked gene in cancer genomic analysis. In this report, we analyzed mutational data from cutaneous melanoma samples in multiple independent datasets. Results showed that the somatic MUC16 mutation was a reliable surrogate biomarker for global TMB across all datasets. Because the discrepancy of mutation rate between mutation groups was larger in melanoma compared with other cancer types, MUC16 can be a more effective biomarker in clinical use for melanoma. One appealing benefit in using MUC16 mutation status as a stable barometer of TMB in melanoma is that MUC16 is lowly expressed (Supplementary Fig. S16) and thus it is less likely to be subject to transcription-coupled repair. MUC16 also has important prognostic implications as MUC16-mutant and lower VAF were found associated with improved survival (but TTN-mutant was not).

To investigate whether MUC16 mutations and MUC16-specific neoantigens can be used as predictive biomarkers for immunotherapy outcomes, we analyzed patient mutation data from four studies with immunotherapies. One noteworthy finding is that the MUC16-specific neoepitopes were significantly associated with total neoantigen burden in tumors in all studies. However, the association of MUC16 mutations with therapy outcome was only observed in one study with anti–CTLA-4–treated patients, where MUC16-mutant patients were found enriched in the group of long-term benefit. In the pooled sample analysis, it was found that the MUC16 mutation status was associated with better survival outcomes in anti–CTLA-4–treated cohorts, but not in anti–PD-1 cohorts. The lack of evidence for prognostic associations in these studies may be due to limited patient sample size and varying definitions of response to immunotherapy.

We also sought to explore germline genetic variants in MUC16, which may alter the gene's function and activity in both cancer and immune cells, and the prognosis of melanoma. Although preliminary, we observed wild-type genotypes in multiple SNPs that appeared to be associated with worse prognosis. A very interesting SNP is rs7245949, which is in strong LD with two other SNPs; having homozygous wild-type genotype was clearly associated with lower expression in markers of T-cell responses (Supplementary Fig. S17) and accordingly worse survival. This missense SNP changes the codon from T (threonine) at position 2891 to I (isoleucine), which is located within the heavily O-glycosylated N-terminal region (4). This makes the hypothesis that MUC16 acts through modulating the synapse between T cells and MHC presentation more plausible. It would encourage further research on how these missense SNPs, discovered mainly in the N-terminus of MUC16, impact upon cancer immunogenicity and prognosis. As a way to further explore whether a gene is functional, we calculated the VAF of all genes in the TCGA melanoma cohort. Interestingly, it was found that the VAF of MUC16 was increased to 50% (ranked 4th in all genes) when we only look at tumors with high purity (Supplementary Fig. S18), which provides additional suggestive evidence that MUC16 might be functional and clonal despite its relatively low expression in melanoma.

Overall, the findings presented support the hypothesis that MUC16 mutations were associated with TMB and prognosis. In addition, we observed that MUC16-specific neoantigens were associated with total neoantigen load in tumors. These results support the use of MUC16 mutation as a potential surrogate biomarker for global TMB and neoantigen load analysis. Although it was hypothesized that MUC16 mutations might act through neoantigen presentation and immunoediting (3), we did not observe a direct relation between its mutation and therapy response based on existing data collected from anti–PD-1 treatments. However, we were able to replicate a predictive model for predicting the response rate to anti–PD-1/PD-L1 therapy across tumor types based on the MUC16 mutation burden, further supporting the proposal that MUC16 provides a potential gauge biomarker for global tumor neoantigens. The findings from this study encourage further studies in MUC16 and other nondriver hotspot mutations, which may contain special secondary structures such as hairpins and loops (19). Mutations in those genes regions are not necessarily associated with oncogenic driver; rather, these are mere passenger mutations, but should be better measurements of overall mutational burden than other regions of the genome.

We argue that the intertumor heterogeneity will have a negligible impact on the mutational analysis based on TCGA, Yale, and Moffitt cohorts, especially for somatic results on the TMB association. However, study heterogeneity and other unknown sources of patient heterogeneity should be carefully investigated in future large-sample analysis for the immunotherapy response prediction.

Collectively, our analysis revealed that MUC16 mutation is a reliable and independent biomarker for TMB in melanoma. Importantly, we identified that MUC16 is associated with survival outcomes across different datasets. It is also, to the best of our knowledge, the first analysis of its kind showing that MUC16 is the only gene among top mutated genes that are associated with survival outcomes in cutaneous melanoma. Although TMB is not a successful biomarker for predicting anti–PD-1/PD-L1 response in cutaneous melanoma, we demonstrated that MUC16 could still be useful as a predictor of neoantigen load and anti–CTLA-4 responses. When combined with CD8⁺ T-cell abundance and PD-1 expression, MUC16 can be used to generate a minimalist predictive panel for predicting the ORR to anti–PD-1/PD-L1 therapy across cancer types. Therefore, we conclude that MUC16 is a valid surrogate biomarker for TMB as well as a promising prognostic biomarker in melanoma. Its role in predicting response to immunotherapy should be further studied when more data have been collected.

C. Duan is a research investigator at Bristol-Myers Squibb. J.K. Teer has ownership interest (including patents) in and is a consultant/advisory board member for Interpares Biomedicine. V.K. Sondak is a consultant for Bristol-Myers Squibb, Novartis, Polynoma, Regeneron, and Merck. J.R. Conejo-Garcia is a consultant for Compass Therapeutics, Anixa Biosciences, and Leidos; reports receiving commercial research grants from Compass Therapeutics and Anixa Biosciences; and has ownership interest (including patents) in Compass Therapeutics and Anixa Bisociences. No potential conflicts of interest were disclosed by the other authors.

Conception and design: X. Wang, Y.A. Chen

Development of methodology: X. Wang, X. Yu, J.R. Conejo-Garcia

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): X. Wang, X. Yu, M. Krauthammer, W. Hugo, K.S.M. Smalley

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): X. Wang, X. Yu, C. Duan, J.K. Teer, Z.J. Thompson, K.Y. Tsai, Y.A. Chen, J.R. Conejo-Garcia

Writing, review, and/or revision of the manuscript: X. Wang, X. Yu, C. Duan, P.A. Kanetsky, J.K. Teer, D. Kalos, K.Y. Tsai, K.S.M. Smalley, V.K. Sondak, Y.A. Chen, J.R. Conejo-Garcia

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): X. Wang, V.K. Sondak

Study supervision: X. Wang, J.R. Conejo-Garcia

This study was supported in part by Moffitt Skin Cancer SPORE grant P50 CA168536 (to V.K. Sondak); Moffitt Skin Cancer SPORE Career Enhancement Program (to X. Wang); Institutional Research Grant 14-189-19 from the American Cancer Society (to X. Wang); and NIH grants R01CA157664, R01CA124515, R01CA178687, and R01CA211913 (to J.R. Conejo-Garcia). Support for shared resources was provided by Cancer Center Support Grant (CCSG) CA076292 to H. Lee Moffitt Cancer Center.

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.