The Distinctive Mutational Spectra of Polyomavirus-Negative Merkel Cell Carcinoma

Merkel cell carcinoma (MCC), or primary cutaneous neuroendocrine carcinoma, is a rare malignancy with high rates of recurrence, metastasis, and mortality. The incidence of MCC has nearly tripled in the past 20 years, and is more prevalent in the immunosuppressed and elderly. Five-year overall survival from time of diagnosis is 30% to 64% (1, 2). Previous studies to elucidate the molecular pathogenesis of MCC found that a subset of cases display TP53-inactivating mutations (14%–28%) and/or PIK3CA-activating mutations (4%–17%; ref. 1). The discovery of Merkel cell polyomavirus (MCPyV) viral DNA via digital subtraction transcriptome analysis in a majority of MCCs represented a major breakthrough (3). MCPyV may likely contribute to tumorigenesis via large T antigen (LTAg) inhibition of the tumor-suppressor RB1, and enhanced oncoprotein gene stability and mTOR activation by small T antigen (sTAg; refs. 1, 2). In MCC, MCPyV displays genomic integration and characteristic truncating mutations of LTAg, which render the virus replication deficient but preserve the RB-binding site (1, 3). In contrast, oncogenic activation events in MCPyV-negative MCC have been underexplored. No targeted therapies are currently available for MCC, although survivin, PI3K, and BCL2 inhibitors may hold promise (1, 4).

Next-generation sequencing (NGS) is a powerful, unbiased approach for identifying novel genetic aberrations in cancer, including point mutations, copy-number gains/losses, gene fusions, and viral sequences (5).Integrative sequencing incorporates data from whole-exome sequencing and whole-transcriptome sequencing to generate a comprehensive landscape of underlying genetic aberrations and outlier gene-expression changes in tumors (5). Recent exome sequencing studies on small cohorts of formalin-fixed, paraffin-embedded (FFPE) MCC samples identified recurrent RB1 mutations in MCPyV-negative tumors, as well as PDE4DIP mutations (6, 7). However, NGS studies of MCC have, thus, far been limited, and detailed somatic mutation and expression analyses of MCC by integrative sequencing have not been reported.

The objective of the MI-ONCOSEQ precision oncology study is to carry out integrative sequencing of tumors from patients with rare or refractory disease toward the goal of identifying novel therapeutic strategies (5). Here, we performed integrative sequencing of biospecimens obtained from two patients with MCC enrolled in the MI-ONCOSEQ study as well as a validation cohort of 14 MCC tumor samples.

Patient samples were procured and profiled under Institutional Review Board (IRB)-approved studies. For MI-ONCOSEQ samples, patients were enrolled and consented through a University of Michigan Hospital System (UMHS) IRB-approved protocol for integrative tumor sequencing, MI-ONCOSEQ (IRB#HUM00046018). Specimen collection through MI-ONCOSEQ has been previously described (5), with tumor purity confirmed on frozen sections by the study pathologists (L.P. Kunju or S.A. Tomlins). Additional samples were procured from the UMHS Cutaneous Surgery and Oncology Program Tumor Bank, as previously described (8), with tumor purity of >70% confirmed on frozen section by the study dermatopathologists (D.R. Fullen or P.W. Harms). For tumor bank specimens, matched normal DNA was extracted from FFPE lymph node tissue using the QIAmp FFPE DNA Extraction Kit. All cases showed classic MCC immunophenotype, including expression of cytokeratin-20 and neuroendocrine marker(s). Sample details, including age, gender, and disease stage are summarized in Supplementary Table S1.

In this study, we generated 59 NGS libraries to characterize a total of 16 MCC patients from two different cohorts, namely MiOncoseq index cases and the validation cohort (Supplementary Tables S1 and S2). Exome libraries of tumor and matched normal genomic DNA and tumor RNA were generated using the Illumina TruSeq DNA Sample Prep Kit. RNA-seq libraries were prepared according to Illumina's TruSeq RNA protocol, using 2 μg total RNA, as previously described (9). Following standard quality control measures, libraries were analyzed by RNAseq (polyA-transcriptome), exome capture-transcriptome, and/or exome capture-genome sequencing (Supplementary Table S1). Paired-end libraries were sequenced with the Illumina HiSeq 2000. Reads that passed the Illumina BaseCall software chastity filter were used for subsequent analysis. Sequence and alignment quality were assessed by FastQC and the Picard package, respectively. Whole-exome library quality parameters are described in Supplementary Table S2.

Paired-end whole Exome fastq sequence files generated by Illumina (2000/2500) were aligned using Novoalign multithreaded (version 2.08.02, Novocraft) to GRCh37 genome build. Post-processing of bam files generated by Novoalign were carried out using SAMtools (version 0.1.19; ref. 10) and Picard (version 1.93). Mutational analysis was carried out on matched normal–tumor pairs using VarScan2 algorithm (Version 2.3.2; ref. 11). The vcf files for somatic mutation were created with SNV positions having base quality phred score of at least q20, ×10 coverage in normal, less than 5% allelic fraction in 1000 Genomes, and at least 10 variant reads with two reads in each strand in the tumor library. SNVs with >5 recurrences reported in the Catalog of Somatic Mutations in Cancer (COSMIC) were also nominated to probe for known recurrent hotspot mutations. Indel analysis was carried out using Pindel (version 0.2.5). Candidate indels were further filtered by the homopolymer/repeat regions, recurrent sequencing artifacts in our compendia (N = 800), and high recurrence in 1000 Genomes, followed by manual curation. Nominated indels and SNV vcf's were then annotated using ANNOVAR package. Copy-number aberrations were quantified and reported for each gene as the copy-number ratio between each tumor and its matched normal sample, with minimum copy-number ratio cutoffs of 1.25 for gains or 0.75 for losses, as previously described (12). Copy-number analysis was performed only on index cases due to high duplication rates in some libraries of the validation cohort. Tumor content was analyzed as previously described (9).

Gene fusion and gene-expression analyses were carried out using Tuxedo suite tools [Tophat 2.0.4; ref. 13 and Bowtie (version 0.12.8)], using the “—fusion-search” option to find candidate fusions. Nominated fusions were manually inspected for annotation and ligation artifacts. Junction reads supporting fusion candidates were realigned using the BLAT alignment tool to confirm fusion breakpoints. Gene-expression analysis was performed on the accepted_hits.bam generated by Tophat as an input for the Cufflinks (version 2.0.2; ref. 14), which performs assembly of transcripts and estimates abundance in the transcriptome library. In addition, 56,369 transcripts from the Ensembl resource (Ensembl 66) were used as an annotation reference to quantify expression of individual transcripts and isoforms, as previously described (9). Unmapped sequences were used for downstream viral screening analysis using Bowtie and BLAT to align reads to all known viral genomes. Sequencing data for MIONCOSEQ cases will be deposited in dbGaP (accession: phs000673.v1.p1).

Statistical comparisons of mutation rates and C > T fraction were performed using the Student t test with Welch correction on GraphPad Prism 6 software.

NOTCH1-4, PRUNE2, GRIN2A, and HRAS somatic mutations were validated by Sanger sequencing at the University of Michigan Sequencing Core following PCR amplification. Because of limitations in frozen tissue availability, DNA for Sanger validation was extracted from archival FFPE tumor material for all cases except MO_1160, METMCC_862, and MCC345. Chromatograms were visualized using Sequence Scanner 2 software.

MCPyV copy number was quantitated by LT2 (LTAg) and SET9 (sTAg) qPCR of MCC tumor DNA along with the MKL-2 cell line as the standard as previously described (Supplementary Table S1; refs. 4, 15, 16). One sample with extremely low viral copy number (0.0006 copies/cell) was classified as MCPyV-negative, consistent with the precedent established by the previous MCC sequencing study (6). In cases with adequate tissue, MCPyV LTAg was detected by IHC with CM2B4 as previously described (15). One sample (METMCC_614) was positive for MCPyV by qPCR but negative by IHC. PCR-Sanger sequencing of LTAg was performed for this sample as previously described (15), which detected a tumor-specific truncating mutation; therefore, this tumor was classified as MCPyV positive.

Two patients with MCC were enrolled in the MI-ONCOSEQ precision medicine program. The first index patient (MO_1109) was a 60-year-old male who presented with MCC on the scalp that was treated with resection and adjuvant radiotherapy. After initial remission, the patient developed extensive metastatic disease that progressed despite chemotherapy. The patient was enrolled in MI-ONCOSEQ and underwent a biopsy of a chest wall metastasis (Fig. 1A). High-quality tumor RNA and DNA was subjected to NGS. There was no evidence of transcripts related to MCPyV or other oncogenic viruses; absence of MCPyV was confirmed by quantitative PCR, and therefore this tumor was considered MCPyV-negative (Supplementary Tables S1 and S3). In our copy-number variation (CNV) analysis, we detected 30 aberrations, including single copy loss of RB1 and a single copy gain of chromosome 1p32.2-36.3 (region that includes MYCL1; Fig. 1B; Supplementary Table S4; ref. 17). The tumor harbored 1,084 nonsynonymous somatic exonic mutations, including a splice site mutation affecting the remaining copy of RB1 (Supplementary Tables S5–S6). A predominance of cytosine to thymidine (C > T) transitions arising at dipyrimidine sites, or CC > TT tandem substitutions, is characteristic of the UV mutational signature (18). We next investigated the presence of such patterns in the SNV data and noted a high occurrence of C > T transitions (Fig. 1C; Supplementary Table S5). Examination of bases immediately 5′ and 3′ to mutated residues (18) revealed an enrichment of C > T transitions at dipyrimidine sites in MCC tumors, a pattern similar to the UV mutational signature in melanoma (Supplementary Fig. S1; ref. 18). The majority of tandem substitutions were CC > TT transitions, also consistent with UV-induced mutations (Fig. 1D; Supplementary Table S5).

The second index patient (MO_1160) was a 66-year-old male with a high-grade neuroendocrine tumor presenting in the parotid gland, favored to represent metastatic MCC of unknown primary. Metastases to bone, liver, and lymph nodes were present at diagnosis. The patient progressed despite radiotherapy and chemotherapy, and was subsequently enrolled in MIONCOSEQ. CT-guided needle biopsies were obtained from a liver metastasis (Fig. 1E). No evidence of oncogenic viruses was detected by NGS or MCPyV qPCR, leading us to conclude that the tumor was MCPyV negative. By CNV analysis, we detected 24 aberrations, including single-copy losses of RB1, TP53, and PTEN genes (Fig. 1F; Supplementary Table S4). Somatic mutation analysis revealed two PIK3CA mutations, including an activation hotspot E545K, and a P146S hotspot mutation in TP53 (Supplementary Table S5). Similar to MO_1109, MO_1160 also demonstrated a high mutation burden (with a total of 1,441 nonsynonymous somatic mutation calls; Supplementary Table S5) dominated by C > T transitions characteristic of a UV-induced mutational signature (Fig. 1G; Supplementary Fig. S1; Supplementary Table S5). CC > TT changes dominated tandem substitutions, again consistent with UV-signature mutations from a cutaneous primary tumor (Fig. 1H). On the basis of the activating PIK3CA mutation in the patient's tumor, enrollment in a PI3K inhibitor trial was recommended during the Precision Medicine Tumor Board discussion.

We next expanded our analysis to a validation cohort of seven additional MCPyV-negative and seven MCPyV-positive MCC tumors, previously classified by qPCR for MCPyV status. Transcriptome libraries were queried for transcripts corresponding to known pathogens. MCPyV transcripts were identified in all tumors where MCPyV was detected by PCR (range of RNAseq viral read counts: 2,335–26,350; Supplementary Table S3). In addition, lower MCPyV reads (range, 39–46) were detected in two tumors that tested negative for MCPyV by PCR and IHC (Supplementary Table S3). Tumor-specific truncating LTAg mutations were not detected in mapped viral reads from these two tumors. Although this finding is of uncertain significance, the much lower number of viral reads in these two samples suggests background low-level viral loads previously reported in non-MCC carcinoma, possibly representing background wild-type viral infection (19). On the basis of the negative qPCR and IHC, these tumors were categorized as MCPyV negative. Other than MCPyV, no oncogenic pathogens were identified.

Across all MCC cases, we detected 5,351 total nonsynonymous mutations, of which 356 mutations were previously reported in the COSMIC database (Supplementary Table S6). MCPyV-negative tumors demonstrated markedly higher mutation rate per megabase than in MCPyV-positive tumors (mean 10.09 ± 2.32 vs. 0.40 ± 0.09, P < 005; Fig. 2A) and significantly higher C > T fraction (mean 0.86 ± 0.01 vs. 0.41 ± 0.03, P < 0.0001; Fig. 2B; Supplementary Table S5). In MCPyV-negative tumors, the majority (mean 92.7%) of tandem substitutions were CC > TT substitutions, with only 1 CC > TT tandem substitution (12.5%) detected across all MCPyV-positive tumors (Fig. 2C; Supplementary Table S5). As in the index cases, mutation signature analysis revealed enrichment of C > T transitions at dipyrimidine sites in MCPyV-negative MCC samples (Fig. 2D; Supplementary Fig. S1; ref. 18). There was no clear mutational signature in MCPyV-positive MCC samples (Supplementary Fig. S2). Mutation rate and C > T fraction were not significantly associated with patient age or primary tumor site. These findings support a molecular dichotomy in MCC with regard to viral status and mutational burden.

To identify mutations with likely functional significance, somatic mutation calls were ranked according to known activating oncogenic mutations and recurrent inactivating tumor-suppressor mutations (Fig. 3). We found MCC tumors to be heterogeneous with regard to candidate oncogenic drivers; however, MCPyV-negative tumors harbored highly recurrent tumor-suppressor mutations.

PIK3CA mutations were present in three tumors, two of which were activating (E545K and K111E; refs. 1, 20). Interestingly, in two cases, we also identified mutations novel to MCC in KNSTRN, a kinetochore gene recently reported to undergo oncogenic mutation in 19% of cutaneous squamous cell carcinomas (SCC; ref. 21). Similar to SCC, KNSTRN mutations in MCC resided in the N-terminus, including an S24F hotspot mutation. One MCPyV-positive tumor harbored two separate HRAS mutations (G12C and G13D) in distinct subclones (Supplementary Fig. S3 and Supplementary Tables S6–S7). PIK3CA, KNSTRN, and HRAS mutations were mutually exclusive in our cohort, suggesting a potential driver role in these samples. Four tumors harbored mutations in the oncogene PREX2, a negative regulator of PTEN. Activating mutations in PREX2 are reported in 14% of melanomas (22). One tumor harbored RAC1 (P29S)-activating mutation, a small G-protein mutated in melanoma (23). Overall, candidate oncogenic-activating mutations were identified in 6 of 8 (75%) and 2 of 7 (29%) of MCPyV-negative tumors and MCPyV-positive tumors, respectively.

MCPyV-negative tumors harbored several highly recurrent mutations in tumor-suppressor genes, including TP53 (7/8, 87.5%), RB1 (5/8, 62.5%), NOTCH1 (4/8, 50%), and PRUNE2 (5/8, 62.5%; Fig. 3; Supplementary Fig. S3, Supplementary Tables S6–S7). PRUNE2 and NOTCH1 mutations have not been described in MCC. PRUNE2, a proapoptotic factor (24), is mutated in 20% of melanomas (The Cancer Genome Atlas). We identified two inactivating mutations in the glutamate receptor subunit GRIN2A (2/8, 25%): A splice site mutation and the (R902K) loss-of-function mutation described in melanoma (Supplementary Tables S6–S7; ref. 25). Furthermore, we discovered a novel missense mutation in the DNA damage response gene BRCA2 (D2819V) affecting the DNA-binding domain (Supplementary Table S6).

Interestingly, mutations affecting one or more NOTCH genes were identified in 6 of 8 (75%) MCPyV-negative MCCs (Fig. 3; Supplementary Tables S6–S7). NOTCH mutations in MCC were mainly located in EGF or ankyrin repeat regions, consistent with inactivating mutations (Fig. 4). We independently validated the somatic mutation calls by PCR-Sanger sequencing and confirmed the presence of NOTCH mutations in tumors, but not in the matched normal samples for 16 of 17 mutations (Supplementary Table S7). The majority of mutations were C > T transitions at dipyrimidine sites. Depending on cancer type, NOTCH signaling may play either an oncogenic or tumor-suppressor role (26, 27). In hematologic malignancies, NOTCH mutations or fusions that disrupt the C-terminal PEST domain result in increased NOTCH stability and aberrant signaling that promotes tumorigenesis. However, NOTCH signaling plays a tumor-suppressor role in SCC and small-cell lung carcinoma (26, 27). The clustering of NOTCH mutations in EGF and ankyrin repeat domains in our MCC cohort is consistent with loss-of-function events, suggesting that NOTCH signaling plays a tumor-suppressor role in MCC similar to other neuroendocrine malignancies.

Insertion/deletion analysis (indels) identified 32 indels (Supplementary Table S8). One tumor harbored an activating frameshift deletion at exon 6 in PPM1D, a negative regulator of p53 (28). Another tumor harbored a frameshift deletion of the DNA damage response gene ATM that is predicted to be inactivating due to loss of the kinase domain.

Fusion analysis revealed a total of 15 calls but no predicted oncogenic driver events were noted (Supplementary Table S9). However, we did identify a highly expressed fusion transcript in MCC_456 between MLH1 and SPATA4 (Supplementary Fig. S4A). This fusion results in loss of the C-terminal PMS2–EXO1 interaction domain of MLH1, a DNA mismatch repair gene, and has the potential to generate an inactive or dominant-negative form of the protein (Supplementary Fig. S4B). We independently confirmed this fusion by Sanger sequencing of the fusion PCR product from the index sample (Supplementary Fig. S4C). In addition, SPATA4 expression was restricted to the index sample, indicating 3′ partner expression driven by the underlying causal fusion event (Supplementary Fig. S4D–S4E). The functional significance of the fusion gene remains to be characterized.

Our studies found that MCCs segregate into distinct molecular classes; MCPyV-negative MCCs have high mutational burdens characterized by UV-signature events, supporting UV-induced damage as an etiology. In contrast, MCPyV-positive MCCs harbor relatively few mutations (25-fold lower than MCPyV-negative tumors) and do not display a definitive UV-signature, supporting an oncogenic role for MCPyV T antigens as primary drivers for these tumors. Our findings are analogous to those in head and neck SCCs, where tumors lacking human papillomavirus display higher mutation burdens (29). However, the observation that a minority of MCPyV-positive tumors harbor activating mutations in oncogenes such as HRAS suggests that genomic mutations may cooperate with the tumorigenic activity of MCPyV in some cases. MCPyV-negative MCCs harbored highly recurrent inactivation of candidate tumor suppressors, including TP53, RB1, PRUNE2, and NOTCH1-4 genes. Our findings suggest that the NOTCH pathway plays a tumor-suppressor role in MCC. In contrast with tumor suppressors, MCCs are heterogeneous with regard to oncogenic drivers. The most frequently perturbed pathway was PI3K, with mutation of PIK3CA and/or PREX2 identified in 5 of 15 (33%) tumors. KNSTRN, RAC1, and HRAS may represent novel oncogenic drivers in a subset of tumors. Altogether, our findings suggest that MCC pathogenesis can be molecularly divided into MCPyV-mediated and UV-mediated etiologies.

No potential conflicts of interest were disclosed.

Conception and design: P.W. Harms, P. Vats, R. Mehra, A.M. Chinnaiyan

Development of methodology: P.W. Harms, P. Vats, Y.-M. Wu, X. Cao, F. Su, R. Mehra

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): P.W. Harms, M.E. Verhaegen, D.R. Robinson, Y.-M. Wu, N. Palanisamy, J. Siddiqui, X. Cao, F. Su, H. Xiao, D.R. Fullen, C.K. Bichakjian, T.M. Johnson

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): P.W. Harms, P. Vats, D.R. Robinson, Y.-M. Wu, S.M. Dhanasekaran, L.P. Kunju, S.A. Tomlins, A.A. Dlugosz

Writing, review, and/or revision of the manuscript: P.W. Harms, P. Vats, M.E. Verhaegen, Y.-M. Wu, S.M. Dhanasekaran, J. Siddiqui, L.P. Kunju, S.A. Tomlins, D.R. Fullen, C.K. Bichakjian, T.M. Johnson, A.A. Dlugosz, A.M. Chinnaiyan

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): P.W. Harms, Y.-M. Wu, J. Siddiqui, R. Wang, H. Xiao

Study supervision: P.W. Harms, A.M. Chinnaiyan

Other (along with F. Su and R. Wang, prepared DNA and RNA, optimized the sequencing sample prep protocol and generated sequencing libraries, QC, and sequenced all the NGS samples, also helped to coordinate and manage the project): X. Cao

Other (consent, collect, annotate from charts, and distribute clinically annotated human samples from the skin cancer Biobank; direct and construct Biobank; review of article): T.M. Johnson

The authors acknowledge Jyoti Athanikar (writing and proofreading); Karen Giles (article submission); Dr. Michael Imperiale (assistance with response to reviewers); Jincheng Pan (assistance with initial fusion analysis); Angela Fullen, Lisa Johnson, and Mandy Davis (tissue microtomy); Nisha Meireles (tumor bank), Shannon Carskadon (PCR and IHC); and patients participating in the MI-ONCOSEQ clinical sequencing program and UMHS Cutaneous Surgery and Oncology Program. N. Palanisamy is a recipient of the Development Award from Melanoma Research Alliance. P.W. Harms is supported by the Dermatopathology Research Career Development Award of the Dermatology Foundation. S.A. Tomlins and A.M. Chinnaiyan are supported by the A. Alfred Taubman Medical Research Institute.

This project was supported in part by the Anatomic Pathology Project Funding Committee of the University of Michigan, Department of Pathology. The MI-ONCOSEQ study is funded by the NIH as a part of the National Human Genome Research Institute's Clinical Sequencing Exploratory Research program.