Genomic profiling reveals mutational landscape in parathyroid carcinomas

The parathyroids are small endocrine glands in the neck that regulate the circulating calcium level by producing and secreting parathyroid hormone (PTH), which acts to maintain homeostasis by regulating bone mineral turnover, renal calcium reabsorption, and dietary calcium absorption from the gut. PTH production is inhibited through a feedback loop when calcium binds to calcium-sensing receptors on the parathyroid cell membrane. Hyperparathyroidism is a state of overactive parathyroid function with excessive circulating PTH. It may cause hypercalcemia and symptoms thereof: bone pain, osteoporosis, fractures, kidney stones or other renal damage. Primary hyperparathyroidism is most commonly caused by benign parathyroid adenomas or hyperplasia, which are treated by surgery. Parathyroid carcinoma (PC), however, is an extremely rare but aggressive and life-threatening form of primary hyperparathyroidism, accounting for less than 1% of cases. Most cases of primary hyperparathyroidism are sporadic; only about 5% are associated with hereditary syndromes such as multiple endocrine neoplasia type 1 (MEN1), autosomal dominant familial isolated hyperparathyroidism (FIHP), and hyperparathyroidism-jaw tumor syndrome (HPT-JT) (1, 2).

The rarity of sporadic PC has created significant challenges to delineate its molecular etiology. Due to high incidence of biallelic inactivating mutations in tumor suppressor gene CDC73 (also known as HRPT2) in familial PCs associated with HPT-JT, somatic mutations in CDC73 were investigated in sporadic cases of PC. Mutational analysis (3–5) revealed that CDC73 somatic mutations occur in a majority of sporadic PCs. At the molecular level, parafibromin encoded by CDC73 functions as a tumor suppressor, potentially by downregulating cyclin D1 and cMyc expression (6, 7), and likely plays a key role in control of cell proliferation, apoptosis, and chromosome stability. However, owing in part to the lack of disease-relevant experimental models, the specific functional links connecting CDC73 mutations to parathyroid cancer development remain largely unknown.

The presence of somatic mutations in tumor suppressor genes of widespread importance in human oncology such as TP53, RB1, and BRCA2 has also been examined and none appear to play a major role in PC (8). However, the likely contribution of additional genes to parathyroid cancer was suggested by recurrent copy number alterations, originally detected by comparative genomic hybridization (CGH) (9–14). For example, recurrent regions of copy gain such as chromosome 1q and 16 and regions of copy loss such as chromosome 1p, 3, 13q, and 14 may contain novel driver oncogenes and tumor suppressors, respectively. CCND1 encoding cyclin D1 was first identified as a human oncogene in parathyroid adenomas through PTH-CCND1 gene rearrangement (15), in which fusion of the PTH 5′ regulatory region to the CCND1 coding region led to high-level overexpression of CCND1 in parathyroid tissues. CCND1 overexpression has been observed in over 90% of PCs (16), and while no analogous gene rearrangements have been described, a high prevalence of CCND1 gene amplification was identified in one report (17). Finally, although not directly connected to clonally selected driver mutations, epigenetic regulation has been reported to play a role in activating the Wnt pathway in PCs, specifically via APC tumor suppressor gene inactivation by DNA hypermethylation in several PCs, accompanied by accumulation of active, nonphosphorylated β-catenin (18).

Recently, comprehensive genomic sequencing has been carried out in attempts to identify novel oncogenes and tumor suppressors in PCs. Whole-genome sequencing (WGS) of a single sporadic PC and a recurrent tumor from the same patient revealed somatic mutations in well-characterized cancer genes such as MTOR, KMT2D (previously known as MLL2), CDKN2C, and PIK3CA (19). Notably, the activating PIK3CA mutation was present only in the primary tumor but not in the recurrence, suggesting that the PI3K pathway might be more important for tumor initiation than for clonal selection in recurrent disease. Additionally, whole-exome sequencing (WES) data were reported on 7 patient-matched PC and normal control DNA samples, and 1 case of 2 primary tumors and matched normal control from the same patient (20). High prevalence of CDC73 mutation was confirmed, and recurrent somatic (and 1 germline) alterations identified PRUNE2 as a candidate driver gene in 18% of the PCs analyzed by WES and by Sanger sequencing of a validation set. Analysis of mutation spectrum has revealed the presence of the APOBEC mutational signature (21) in PCs.

In the current study, we generated WES data from 10 sporadic PCs with patient-matched normal controls. Data from these 10 patients were pooled with the previously published 8 patients (20) to increase the statistical power of the analyses in seeking recurrent lesions. Somatic mutations, germline variants, and somatic copy number alterations were identified. A panel of selected genes was used to sequence an additional 7 sporadic PCs. We show that 42% (i.e., 10 of 24) of patients harbored mutations of CDC73, with 4 patients harboring both germline and somatic inactivating variants. Genes encoding the PI3K/AKT/mTOR pathway are altered in 21% of the 24 cases via activating mutations in PIK3CA and MTOR. CCND1 amplification was observed in 29% of the 17 tumors analyzed by WES. We identified the first sporadic PC cases with somatic mutations in genes in the Wnt canonical pathway. Genes involved in chromatin remodeling are also frequently mutated in our study. In addition, we have identified a new recurrent somatic mutation in PC, p.I482M in the putative kinase ADCK1, strongly suggesting that it plays an oncogenic role. Furthermore, we discovered several genes recurrently mutated including AKAP9 in 3, ZEB1 in 3, and FAT3 in 2 of the 17 cases. The new data generated from this study led to the discovery of a significant fraction of the tumors harboring somatic mutations in the PI3K/AKT/mTOR pathway, the first PCs with somatic mutation in a Wnt pathway component, and somatic alterations in recently identified tumor suppressors involved in chromatin remodeling. Pooling the new data with previously published 8 cases (20) allowed us to uncover the recurrent mutation in ADCK1 as well as recurrently mutated genes such as FAT3 as newly recognized candidate cancer genes in PC. To our knowledge, this genomic sequencing analysis represents the largest such study of sporadic PCs, and is an important step toward a full molecular characterization of this extremely rare malignancy, which in turn will inform improved and individualized treatments in the future.

Whole-exome analysis

A diagram illustrating the study design is shown in Figure 1. Briefly, we selected 10 cases of sporadic PC fulfilling stringent selection criteria — clinically sporadic presentations and demonstrating local invasion of tumor into surrounding tissues and/or distant metastasis (22–24). As high-throughput WES provides an unbiased view of the driver landscape, we generated WES data from 10 tumor samples (PC01T–PC10T) with matched normal controls (PC01N–PC10N) from these 10 patients in the current study. The new data from these cases were pooled with the raw sequencing data from 9 tumors (PC11T–PC19T; PC17T and PC18T are from the same patient) with matched normal controls (PC11N-PC17N, PC19N; PC17N is the control for both PC17T and PC18T) from 8 previously published cases (20) for joint bioinformatics analyses, thereby comprising our discovery cohort (Figure 1A). After detailed analysis, samples PC11N–PC17N and PC11T–PC18T were included for further analyses as they passed our quality control (QC) standards.

Figure 1

Flowchart of the study design.

From our discovery cohort (35 WES datasets across 17 patients), we obtained approximately 9.7 billion raw reads (mean = 3.4 × 10⁸/2.2 × 10⁸ and range = 1.8 × 10⁸ to 7.3 × 10⁸/6.5 × 10⁷ to 3.4 × 10⁸ for tumor/normal, respectively). After alignment to the hg19 reference genome, we observed a mean usable sequencing depth of 84.5× (73.6× and 94.8× for tumor and normal, respectively). Approximately 86% of the exome was sequenced to at least 20× depth, enabling high confidence variant calling (Supplemental Table 1; supplemental material available online with this article; https://doi.org/10.1172/jci.insight.92061DS1).

After variant calling, an average of 122.6 somatic (tumor-specific) mutations including single nucleotide variation (SNV) and small insertion/deletion (indel) per tumor (range 20–464; total 2,207) were observed in the discovery cohort. We noted an average mutational burden of approximately 3.6 mutations per 10⁶ bp (average mutation rate = 122.6 mutations per tumor, genome coverage = 34.0 Mbp, Figure 2A). Considering the predicted functional impact (SnpEff effect impact), variants were classified into synonymous or protein-altering (Supplemental Table 2). Filtering resulted in 1,676 protein-altering somatic variants with 1,538 SNVs and 138 small indels (Supplemental Table 3). We observed 8.7% (146 of 1676), 81.7% (1369 of 1676), 2.5% (41 of 1676), 1.7% (27 of 1676), and 5.5% (93 of 1676) nonsense, missense, splice altering, in-frame indel, and frameshift mutations, respectively.

Figure 2

Somatic mutation landscape of parathyroid carcinoma — summary of results for whole-exome and targeted gene panel sequencing assays (n = 25). (A) Frequency of synonymous and protein-altering mutations per Mbp in the discovery cohort. (B) Somatic mutational signatures and number of somatic single nucleotide variation (SNV) mutations (in parentheses) in the discovery cohort. (C) Recurrently mutated genes in the study (red = copy number gain, orange = missense mutation, and gray = truncating mutation). Genes shown are covered by the ParThy targeted panel and are recurrently mutated in the combined discovery and expansion cohorts, ordered by the frequency of mutation.

Mutational landscape in PC

Recurrently mutated genes. The landscape of somatic mutations and copy number alterations is shown in Figures 2 and 3, respectively. Assessing protein-altering somatic mutations at the gene level, there are 92 genes that are mutated in more than one case in the discovery cohort (Supplemental Table 4). We found the following cancer-associated genes most frequently mutated across our WES discovery cohort: CDC73 (8 samples), AKAP9 (3 samples), and ZEB1 (3 samples). Mutation of CDC73 (encoding parafibromin) is the only established oncogenic driver in PC (3–5). In our discovery cohort, we observe 15 distinct variants in CDC73 (9 somatic and 6 germline mutations across 8 tumor samples and 6 normal controls; Supplemental Table 5 and Figure 2C). As shown in Supplemental Figure 1, all mutations are inactivating (nonsense or frameshift) with 1 exception (missense p.L95P mutation), consistent with parafibromin’s established role as a tumor suppressor. As observed in previous studies (20), CDC73 somatic variants in our cohort have a high variant allele fraction (mean mutant allele fraction = 0.35), reinforcing the early role of CDC73 inactivation in PC. Five cases had likely biallelic inactivation of CDC73 detectable, with inactivation of 1 copy in the germline in 4 patients. Further screening of the expansion cohort led to discovery of CDC73 mutations in 3 additional sporadic tumor samples (for details, see ParThy targeted panel section). Two missense and 1 nonsense somatic mutations were identified in the AKAP9 gene. AKAP9 encodes a member of the A-kinase anchor proteins (AKAPs) regulating cellular localization and function of protein kinase A. Although AKAP9 is frequently mutated in oral squamous cell carcinoma (25) as well as gastric and colorectal cancers with high microsatellite instability (26), it is unclear if those mutations are oncogenic. We also found somatic mutations in ZEB1 in 3 cases. ZEB1 is a transcriptional repressor that mediates tumor invasion and metastasis by promoting epithelial-mesenchymal transition (EMT) (27). However, activating somatic mutations of ZEB1 have yet to be reported.

Figure 3

Somatic copy number alterations (SCNAs) in the discovery cohort (n = 17). (A) Log2 of the sequencing depth ratio of tumor versus normal sample is plotted across the samples, with red depicting gain and blue depicting loss. The color intensity is proportional to the ratio. (B) Recurrent SCNA events as predicted by GISTIC2.

We discovered a recurrent missense p.I482M mutation in AarF domain–containing kinase 1 (ADCK1), in 2 tumors in our discovery cohort (PC09T and PC18T; Supplemental Figure 2). This is the only recurrent somatic mutation (at either the DNA or amino acid level) identified in these WES data. ADCK1 encodes a relatively understudied putative kinase with no published biochemical data, although there is evidence of protein expression in parathyroid tissue from The Human Protein Atlas project (http://www.proteinatlas.org/ENSG00000063761-ADCK1/tissue/parathyroid+gland). Among other genes mutated in 2 samples, FAT3 is of particular interest as FAT3 belongs to a family of FAT atypical cadherins including 2 known tumor suppressors, FAT1 and FAT4. We detected mutations of FAT3 in 2 tumors (PC04T and PC14T; Supplemental Figure 3). These 2 mutations are nonsense mutations, with the premature stop codon occurring near the N-terminal amino acid sequences, clearly indicating they are inactivating. In addition, there are several known cancer genes that are mutated in 2 cases, and they are described below.

Somatic mutation in known cancer genes. To determine known cancer genes that are mutated in the discovery cohort, we intersected the list of mutated genes with Cancer Gene Census (version 77, downloaded 05/26/2016) (28) (Supplemental Table 6). Notably, we detected a functionally established activating mutation in PIK3CA (p.K111E) (29) in tumor PC10T (Figure 4), and a known activating mutation in MTOR (p.L1460P) (30) in PC02T (Supplemental Figure 4). In another case, PC15T, we discovered a missense mutation (p.Q2524L) near the kinase domain of MTOR. Although it has yet to be biochemically validated as an activating mutation, this mutation has been previously identified as a somatic alteration in clear cell renal carcinomas (31) and endometrial carcinomas (32). Furthermore, this mutation is predicted to have high functional impact (dbNSFP radial support vector machine [SVM] score = 0.98, 97th percentile; see Methods). Recurrence of this somatic mutation in other cancers and high functional impact score suggest it is likely an activating mutation. Mutations in canonical Wnt pathway genes were found in 2 cases. In PC02T, a missense mutation (p.E1284K) in APC was identified. Located in a mutation cluster in codons 1248–1522 of the APC gene, this mutation was previously discovered in a small-cell lung cancer (33). In addition, we detected an inactivating mutation in RNF43 (p.G659fs) in case PC05T. RNF43 encodes a transmembrane E3 ubiquitin-protein ligase that negatively regulates the Wnt signaling pathway by selectively ubiquitinating frizzled receptors (34). Frequent RNF43 inactivating mutations have been discovered in several cancer types characterized by Wnt pathway dysregulation such as colorectal, endometrial, and gastric cancers (35).

Figure 4

Mutation distribution for the PIK3CA gene in the study cohort (n = 25) (top) and all cancer studies in cBioPortal (bottom) along the body of the gene. The colored rectangles represent the known functional domains of the translated protein (Pfam database; http://pfam.xfam.org/).

We observed frequent mutations in genes that directly regulate chromatin structure, including truncating mutations in a recently discovered tumor suppressor, ARID2, implicated in hepatocellular carcinomas and melanomas (36–38); a missense mutation (p.L1549V in PC09T; dbNSFP radial SVM score = –0.83, 54th percentile) in a histone demethylase KDM5C (also known as JARID1C) that is mutated in 9% of clear cell renal cell carcinomas and its inactivation triggers genomic instability (39). In addition to ARID2 and KDM5C, there are several other chromatin remodeling genes that are not included in the Cancer Gene Census but are mutated in our discovery cohort: a frameshift indel (p.K512fs) in another ARID family gene ARID4A in PC04T; missense mutations or in-frame indel in 3 other histone demethylases (KDM4C, p.R919K in PC18T; KDM4E, p.R100H in PC04T; JMJD1C, p.K593_S601del in PC04T); a missense mutation (p.S1128L in PC02T) in the SET domain histone methyltransferase SETD1B that is frequently mutated in gastric and colorectal cancers (40). Cumulatively, 29% (5 of 17) of the discovery cohort in this study harbor mutations in a chromatin remodeling pathway.

Functional studies and recent genetic data suggested that NOTCH1 is an oncogene in hematological malignancies and a tumor suppressor in squamous cell carcinomas (41, 42). We found NOTCH1 mutations in 2 tumors: p.T194P in PC10T, and p.Q439* in PC11T. As a protein-truncating mutation, p.Q439* is very likely inactivating. The p.T194P mutation has been previously identified in head and neck squamous cell carcinomas (43) where NOTCH1 has been postulated as a tumor suppressor (44, 45), suggesting it is also likely inactivating. Therefore, our findings suggest that NOTCH1 is a potential tumor suppressor in 2 PCs in this study.

Although not included in the Cancer Gene Census list, the prune homolog 2 (PRUNE2) gene has recently been implicated in PC as a potential tumor suppressor (20). In the discovery cohort we observed the previously identified somatic stop-gain mutations (p.E474* and p.E537*) in 1 sample (PC16T), but found no additional somatic mutations in PRUNE2. Additionally, germline missense variants were found in PC12N (p.E2570A), PC14N (p.V452M), and PC16N (p.L378V) normal control samples.

Mutation signatures. The somatic SNV distributions are shown in Figure 2B. Typically, G:C>A:T and A:T>G:C somatic transitions are abundant in cancers (21). However, in PC we observe G:C>C:G (n = 13) and G:C>T:A (n = 3) transversions at elevated rates (ranked higher than expected), suggesting the presence of mutator phenotypes causing uncommon mutation types. This finding corroborates an earlier finding of an apolipoprotein B mRNA editing enzyme, catalytic polypeptide-like (APOBEC) mutational signature (21) in a subset of PC specimens (20).

Somatic copy number alterations. The spectrum of somatic copy number alterations in the cohort is summarized in Figure 3, with details shown in Supplemental Table 7. We observed recurrent loss of chromosomes 1p (8 samples), 3 (8 samples), and 13 (10 samples) and recurrent gain of chromosomes 1q (6 samples), 16 (9 samples), and 20 (9 samples). These findings corroborate observations from earlier studies (9–14). GISTIC2 (46) was used to identify focal recurrent copy number events. We observed significant (GISTIC q value < 0.25) gain of 6p11.2, 7q22.1, 10q26.2 and 11q23.2 and loss of 1p36.32, 3q29, 6q25.3 and 7p21.3 in our cohort (Figure 3B). We detected gain of genomic regions encompassing CCND1 in 5 of the 17 (29%) cases in the discovery cohort (Supplemental Table 7). Notably, 4 of the 5 cases are mutually exclusive from the 8 cases harboring CDC73 somatic mutations (Figure 2; P value = 0.11, one-sided Fisher’s exact test).

Biallelic gene inactivation. Next we investigated genes that underwent likely biallelic inactivation, i.e., either one copy of the gene was mutated in the germline and the other in the tumor or both copies were altered in the tumor. The findings are listed in Supplemental Table 8. Briefly, biallelic inactivation of a gene was observed in 15 of 17 cases from our discovery cohort. We observed dual protein-altering mutations across multiple samples for CDC73 (PC01, PC02, PC03, PC09, PC17), AKAP9 (PC02, PC07), and FAT3 (PC04, PC09).

ParThy targeted panel sequencing

Combining our findings from the WES discovery cohort with the literature, we developed a targeted cancer gene panel, called the ParThy panel, covering all coding exons (UniProt canonical isoform) of 16 genes implicated in PC (Figure 1B). The genes to be included in the panel were chosen based upon recurrence of mutations in our discovery cohort, biological significance, and/or literature review (Supplemental Table 9). In addition to being a complementary technology to Illumina, Ion AmpliSeq was chosen as the selected sequencing platform for the panel as it can yield sequencing libraries from low-input DNA derived from formalin-fixed, paraffin-embedded (FFPE) samples, has a relatively short turnaround time, and provides ultradeep coverage to detect low allelic fraction variants based on our prior experience (47). As shown in Figure 1B, we screened 10 samples from the WES cohort and an additional 12 tumor samples (no matched normal samples; expansion cohort) using the panel. From the tumor-only samples, 7 of 12 samples were completed successfully, i.e., yielding data passing our QC metrics. Sequencing depth, uniformity, and other QC metrics are displayed in Supplemental Figure 5 and listed in Supplemental Table 10. Importantly, all somatic mutations detected from the WES assay were confirmed by the panel results.

Furthermore, we observed additional mutations in known and putative PC genes from the panel results (Supplemental Table 11). Truncating or frameshift mutations in CDC73 were observed in samples PC20T, PC28T, and PC31T in the expansion cohort, furthering its role as a driver in PC. PIK3CA variants were detected in 3 cases in the expansion cohort (Figure 4), with 2 of these being known activating mutations, p.H1047R (48) and p.E545A (49). The 2 activating PIK3CA mutations are mutually exclusive from the 3 inactivating CDC73 mutations in the expansion cohort. In the combined discovery and expansion cohort of 24 cases, the 11 tumors harboring CDC73 mutations and the 3 tumors harboring the activating PIK3CA mutations are mutually exclusive, although the limited case numbers prevent this trend from reaching statistical significance (Figure 2; P value = 0.14, one-sided Fisher’s exact test). Additional mutated genes in the expansion cohort include FAT3, SETD1B, and PRUNE2 in 4, 2, and 1 cases, respectively.

PC is one of the rarest human cancers, generally accounting for less than 1% of all cases of primary hyperparathyroidism worldwide, although significant regional variation may exist (2). This rarity has created a tremendous challenge to studying the genetic basis of disease etiology by sequencing patient tumor samples. While The Cancer Genome Atlas (TCGA) and The International Cancer Genome Consortium (ICGC) have generated whole–genome/exome sequencing data of approximately 25,000 tumor samples covering more than 50 cancer types (50, 51), including 10 rare cancers (http://cancergenome.nih.gov/cancersselected/RareTumorCharacterizationProjects), they do not include PC, and thus to the best of our knowledge, genomic sequencing has been reported on only 9 sporadic and familial PC patients (19, 20). In this study, we generated exome sequencing data from paired tumor-normal samples of 10 sporadic PCs. In order to increase the statistical power, data analysis was performed by combining the 10 WES cases with 7 of the previously reported cases (20), representing the largest genomic sequencing study of sporadic PC to date. Furthermore, in PC samples from 7 additional patients, we also sequenced 16 PC-related genes on a high-throughput sequencing panel informed by our WES data. Our results have revealed newly recognized driver mutations and potential genetic mechanisms in PC development, in addition to confirming CDC73 as a key tumor suppressor. We applied 2 complementary approaches to elucidate the somatic mutations identified in our study. A knowledge-driven approach was utilized to determine if well-known oncogenic pathways are genetically altered in our study cohort, and a data-driven approach was used to identify recurrently mutated genes for discovery of novel oncogenic drivers.

The PI3K/AKT/mTOR pathway plays a critical role in regulating cell cycle progression and is also a key regulator of survival during cellular stress. Dysregulation of the PI3K/AKT/mTOR pathway profoundly disturbs the balance between cell proliferation and apoptosis, and ultimately leads to a competitive growth advantage, metastatic competence, angiogenesis, and therapy resistance in cancers (52). Frequent genetic alterations such as somatic mutations and copy gain/loss of key components in the PI3K/AKT/mTOR pathway have been reported in many solid tumors as well as hematological malignancies. In a meta-analysis of cancer genes across 21 tumor types, PIK3CA (encoding the catalytic subunit of PI3 kinase) and PTEN (encoding an inhibitor of the pathway) are among the top 3 most significantly mutated genes (53). In the current study, we detected functionally established activating PIK3CA mutations in 3 carcinomas in the combined discovery and expansion cohort. Mutual exclusivity of the 3 activating PIK3CA mutations from CDC73 somatic mutations in our study further suggests they are independently and crucially oncogenic. Notably, in the previous WGS study of a single sporadic PC, the primary tumor sample harbored an activating PIK3CA mutation and has wild-type CDC73 (19), consistent with our findings of mutual exclusivity between PIK3CA and CDC73 mutations. Furthermore, MTOR gene mutations were present in 2 other carcinomas, mutually exclusive from the activating PIK3CA mutations. Therefore, the PI3K/AKT/mTOR pathway is very likely activated in 21% (5 of 24) of the PCs in this study, representing a major oncogenic pathway and potentially a therapeutic target for sporadic PCs. Moreover, PIK3CA/MTOR mutations could potentially provide value in the diagnosis of PCs, if their still-unknown prevalence in benign parathyroid neoplasms proves to be sufficiently rare or nil.

The Wnt pathway has been recognized as a potentially key oncogenic pathway in PCs. In a study of 5 sporadic PCs, tumor tissues of all 5 cases exhibited increased nuclear accumulation of active, nonphosphorylated β-catenin in comparison with adjacent normal control samples (18). It was suggested that this activation of the Wnt pathway was due to epigenetic regulation, as the APC gene was hypermethylated in these PC tumor samples. Although parafibromin/CDC73 might play a role in regulating Wnt signaling, as dephosphorylated parafibromin stabilizes β-catenin in a gastric cancer cell line (54), somatic mutations or copy number variations (CNVs) of canonical Wnt pathway components such as activating mutations in β-catenin or loss of APC have not previously been identified in sporadic PCs. In this study, we detected a likely inactivating APC somatic mutation (p.E1284K) in tumor PC02T and an inactivating RNF43 somatic mutation (p.G659fs) in tumor PC05T. RNF43 is mutated in more than 18% of colorectal adenocarcinomas and endometrial carcinomas, most prevalently in microsatellite-unstable (MSI-H) colorectal tumors (80%) and MSI-H endometrial tumors (51%) (35). Among these mutations, 73%–75% were truncating events and the p.G659fs mutation observed in our study accounted for 42% to 48% of the RNF43 mutations in colon and endometrial cancers. Moreover, truncating mutations of RNF43 are mutually exclusive from inactivating APC mutations in colorectal adenocarcinomas, strongly suggesting that the inactivating RNF43 mutation is a bona fide oncogenic event (35). To the best of our knowledge, the APC and RNF43 mutations in this study represent the first sporadic PCs with somatic mutations putatively mediating Wnt pathway activation, which also complements previously reported epigenetic mechanisms activating Wnt pathway in PCs (18).

Amplification of CCND1 was detected in 5 of the 17 (29%) cases in the discovery cohort. Mutual exclusivity of 4 of the 5 cases from the 8 cases harboring CDC73 somatic mutations strongly suggest that gene amplification in these cases is an alternative genetic mechanism to CDC73 inactivation to upregulate CCND1 expression. Previously it has been reported that CCND1 amplification was detected in 5 of 7 PCs based on a multiplex ligation-dependent probe amplification (MLPA) assay (17). Our study provides important confirmation that CCND1 is indeed frequently amplified in PCs. It should be noted that CCND1 overexpression is frequent in benign parathyroid adenomas. Thus, in the carcinomas, other mutations (e.g., PIK3CA, etc.) that occur concomitantly with CCND1 amplification probably synergize in yielding the malignant phenotype.

In addition to delineating somatic alterations of well-established oncogenic pathways, we also focused on several highly recurrently mutated genes for discovery of novel cancer drivers in PCs. We have identified a new recurrent somatic mutation in PCs — p.I482M in ADCK1. ADCK1 is a putative kinase with unknown molecular function. Intriguingly, a kinome-wide functional screening has indicated that ADCK2, a kinase closely related to ADCK1, enhances the stability of HIF-1α, a transcription factor that plays a central role in tumor progression by regulating genes involved in cancer cell survival, proliferation, and metastasis (55). We detected mutations of FAT3 in 6 cases. While a functional role of FAT3 in cancers has not been described, recurrent somatic mutations of the closely related gene FAT1 have been identified in glioblastomas, head and neck squamous cell carcinomas, and colorectal carcinomas (56). Moreover, it was shown that inactivation of FAT1 via mutations promotes Wnt signaling and tumor development, and is associated with poor survival (56). FAT4 was frequently mutated in gastric adenocarcinomas and breast cancers, and it exhibited tumor-suppressor activity in several functional assays (57, 58). Given that both FAT3 somatic mutations in the discovery cohort are truncating, it is quite conceivable that FAT3 is also a tumor suppressor, possibly regulating Wnt signaling. In addition, AKAP9, a regulator of protein kinase A, and ZEB1 promoting tumor metastasis are each mutated in 3 cases in the discovery cohort. The consequence of these mutations with respect to protein functions warrants further investigations to assess their potential roles as oncogenic drivers in PC development.

While the new genomic data generated from 17 patients in our study were predominantly responsible for the major discoveries in this report, it is worth pointing out that our approach of analyzing them together with primary data files from the previous report (20) clearly enhanced the statistical power to detect the recurrent nature of a specific mutation (p.I482M in ADCK1) as well as recurrently mutated genes (e.g., FAT3), which would not have been apparent from either dataset in isolation. That said, 4 of the 5 activating mutations in PIK3CA or MTOR are from the new WES and ParThy panel sequencing data generated in this study, allowing us to reveal PI3K/AKT/mTOR as a likely major oncogenic pathway in sporadic PCs. We also identified the first 2 PC cases with somatic mutation of a Wnt pathway component from the new WES data. Moreover, we were able to identify from the new data possibly inactivating mutations in several recently characterized tumor suppressors involved in chromatin remodeling such as ARID2, SETD1B, and KDM5C.

Our discovery of recurrent, known activating mutations in the PI3K/AKT/mTOR pathway and frequent CCND1 amplifications in PCs immediately points to possible therapeutic avenues. However, the rarity of PC limits the possibility of a randomized, controlled clinical trial to test the efficacy of PI3K/AKT/mTOR inhibitors or CDK4/6 inhibitors specifically in this disease. To overcome this obstacle, 2 strategies may be considered. First, PC patients could be considered for enrollment into basket trials (59), where they may be assigned to a targeted treatment arm based on their genomic alterations (e.g., activating PIK3CA mutations) rather than their tumor type. Second, PC patients could undergo genomic testing of their tumor using a multiple-gene panel that includes at a minimum PIK3CA, MTOR, and CCND1, and possibly other genes in pathways considered actionable; the content of such a panel would be expected to evolve with the outcome of relevant trials and with the development of new therapeutic agents. As stated in a joint recommendation by the Association for Molecular Pathology, American Society of Clinical Oncology, and College of American Pathologists (60), such genomic testing may provide an interpretive comment on detected genetic alterations that puts the alteration in clinicopathologic context to inform management decisions. Given the rarity of this tumor type, accumulation of individual case reports on outcomes of genomic testing and therapeutics would greatly facilitate validating the utility of genetic alterations as drug response markers to ultimately guide routine clinical practice.

We recognize several limitations in our study. As noted above, due to the rarity of diagnosed sporadic PCs, we only analyzed exome sequences of 17 cases. The average mutation load is 3.4 mutations per Mb. At this background mutation rate, it would require 160 tumor samples to detect genes mutated at 10% frequency, and approximately 470 samples to detect genes mutated at 5% frequency to achieve 90% power for 90% of genes (53). Thus, because of the limited sample size available for this rare tumor, we may well have missed some mutated genes and pathways of significance in small subsets of PC. Also, the ParThy targeted panel we designed is necessarily selective and of limited scope, and thus is unlikely to reflect a complete spectrum of key oncogenic pathways in PC tumorigenesis. In addition, ParThy panel sequencing in the 7 expansion cases was performed in tumor samples only. Although we removed common genetic variants (>2% allele frequency in ESP6500, 1000Genomes, and ExAc; see Methods), it is possible some of the mutations we identified in the expansion cohort (Supplemental Table 11) are rare or private germline variants rather than somatic mutations. Another limitation of our study is that due to restricted sample availability, we were unable to generate genomic data in other platforms such as gene expression profiling. The genetic alterations identified through exome sequencing would be even better interpreted if integrative data analysis could have been performed with genomic data based on multiple platforms, as demonstrated in TCGA. Gene expression data would allow us to corroborate and extend CNV analysis results such as CCND1 amplification. Nevertheless, the current study represents the largest genomic sequencing effort and analysis for an extremely rare malignancy. Our results have indeed revealed important insights into the genetic landscape of PCs. The PI3K/AKT/mTOR pathway has been identified as a potentially major oncogenic pathway. The highly specific recurrent mutation in ADCK1 is highlighted as a key putative oncogene, and several recurrently mutated genes such as FAT3, AKAP9, and ZEB1 are also revealed as excellent candidates that warrant further functional assessment. Furthermore, to facilitate future studies of PCs and to ultimately develop novel therapeutic strategies, we advocate an expanded consortium-based collaboration to quickly assemble a large sample repository in order to achieve improved statistical power and analysis on multiple platforms.

Patients and tumor specimens. Primary, recurrent, or metastatic PC specimens (tumors PC01T–PC10T) were obtained from patients who had been treated surgically for primary hyperparathyroidism. Patient characteristics and tumor selection for cases PC11T–PC19T are described previously (20). Immediately after surgical resection, tumor samples were either frozen in liquid nitrogen for storage at –80°C until use (n = 12), or were FFPE, sectioned, and mounted on slides, with 2 or 3 sections per tumor, and then dissected from the slides for DNA extraction (n = 10). Patients were diagnosed with PC according to stringent clinicopathological criteria (22–24), namely evidence of either local invasion into surrounding tissues (including vascular invasion beyond the parathyroid tumor capsule) and/or distant metastasis. Paired peripheral blood or other nontumor tissue was available from 10 patients to serve as a source of matched germline DNA.

The following clinical information pertains to the 10 patients in the discovery cohort from which the WES data were generated in this study, and the 12 patients in the expansion cohort. At initial parathyroidectomy, patients ranged in age from 20 to 70 years, mean 47.7 (data available for n = 19). At the time of the surgery that yielded the study specimens, patients ranged in age from 26 to 75 years, mean 52.3 (data available for n = 18). Seven patients were women and 15 were men. The patients’ preoperative serum calcium levels ranged from 12 to 20.6 mg/dl (n = 14). The patients’ preoperative serum levels of PTH were 2.5- to 31-fold above the upper limit of normal for the assay used (n = 10). Eight of the specimens represented primary tumors, 6 represented local recurrences, and 8 represented metastases (3 to lung, 2 to lymph nodes, 3 to bone). Symptoms were noted in available records for 16 patients: 9 had bone pain, 2 had bone fractures, 3 had nephrolithiasis, and 2 had renal failure. At least 17 patients underwent multiple operations for recurrent cancer.

All patients were considered to have sporadic disease upon original presentation, except for one patient who described a vague family history of isolated hyperparathyroidism and had no clinical evidence of HPT/JT. One patient had previous external beam radiotherapy to a recurrent PC lesion other than the one excised and analyzed in this study. One other patient had adjuvant radiation therapy to the cervical region that, 5 years later, yielded the study sample. No patient had been treated with chemotherapy.

Genomic DNA (gDNA) was isolated from tumor and nontumor tissues using either proteinase K digestion followed by phenol-chloroform extraction and ethanol precipitation for frozen tissues and blood, or with the Qiagen Supplemental protocol for purification of gDNA from FFPE tissue using the QIAamp DNA FFPE Tissue Kit (Qiagen, 56404) and Deparaffinization Solution (Qiagen, 19093).

WES. gDNA samples were assessed for quantity by Qubit fluorometry (Life Technologies) and for quality by the 2100 Bioanalyzer system or 2200 TapeStation (Agilent). Initial shearing of 0.5–1 μg gDNA to a mean of 200- to 300-bp fragments was performed using the Covaris E210 focused acoustic energy system. Whole-genome libraries were prepared using the NEBNext DNA Library Prep kit according to the standard manufacturer’s protocol (New England Biolabs). Illumina-compatible paired-end adapters were used and the adapter-ligated DNA fragments were amplified by ligation-mediated PCR (KAPA Biosystems) using a reverse PCR primer containing a 6-nt barcode that allowed for multiple samples to be pooled and sequenced in the same run. The library was enriched for exonic sequences with either the SeqCap EZ Human Exome Library v3.0 capture system (Roche NimbleGen) or the SureSelect Human All Exon V5 system (Agilent) using the manufacturer’s protocols. The libraries were then sequenced with a 100-bp paired-end protocol on the Illumina HiSeq 2500 instrument according to the standard manufacturer’s protocol (Illumina). For the SureSelect-based preparations, the adaptors and amplification primers that came with the SureSelect kits were used. Approximately 8–14 pmol of the whole-exome libraries were clustered and run on the High Output HiSeq 2500 flowcell and sequenced for 100-bp paired-end reads according to the standard manufacturer’s protocol (Illumina).

Identification of genetic alterations. All exome sequencing data were processed through an internal QC pipeline and variant calling and filtering were carried out as previously described (47). All FASTQ files (both normal and tumor; both internal and external data) were combined into a cohort and run through an in-house pipeline (61) to yield BAM and VCF files with germline and somatic variant calls (SNVs and sufficiently small indels). Briefly, this in-house pipeline implements Genome Analysis Toolkit (GATK) (62) version 3.2 best practices for alignment (to hg19 human genome assembly), base quality recalibration, variant calling (using HaplotypeCaller), and variant quality score recalibration (using VQSR set to 99.5% sensitivity) (63). For somatic variant calling, MuTect (64) (version 1.1.6-10b1ba92, HC+PON mode with default settings, using COSMIC [ref. 32] version 65, dbSNP [ref. 65] version 137, and using variant calls from patient-matched normal control as the “panel of normals” setting) and Varscan2 (66) (used only for indels; version 2.3.5, with flags --tumor-purity 0.7 and --min-var-freq 0.07) were used. All variant calls were loaded into a custom MySQL database using in-house scripts and annotated with RVS (67). All germline calls marked “PASS” by VQSR were retained. Somatic calls whose population allele frequency in either ExAC, ESP6500, or 1000Genomes (as described previously in ref. 67) exceeding 2% were discarded on the presumption that they are any combination of: contamination, a variant present but missed in a normal sample, a low-level artifact, could not be a cancer driver because too common in general population. Variant calls were manually reviewed in IGV (68) and the UCSC Genome Browser (69) to inspect supporting alignment quality and alignability of the genomic region in the hg19 human genome assembly. Uncertain calls were manually rejected at this step. In order to identify germline variants for potential disease-predisposing mutations, we filtered variant calls from GATK HaplotypeCaller using the following criteria: high or moderate impact as predicted by SnpEff version 4.0b (described previously in ref. 67), in a canonical isoform according to UniProt v2014_11 (70) or ENSEMBL v78 (71), less than 2% population frequency in ExAc, ESP6500, and 1000Genomes.

We used dbNSFP v2.5 (72) SVM-based ensemble prediction score (RadialSVM), which incorporates 10 scores and the maximum frequency observed in the 1000Genomes populations, to determine variant impact for SNVs. A larger value means the SNV is more likely to have an impact and the scores range from –2 to 3. RadialSVM scores were ranked among all RadialSVM scores in dbNSFP where the rankscore equals the ratio of the rank of the score over the total number of RadialSVM scores in dbNSFP and ranges from 0 to 1.

Somatic copy number alterations were identified using the saasCNV v0.3.1 pipeline implemented in R (73). Briefly, from WES data, the pipeline begins with input heterozygous SNV calls from the normal control sample, and gathers their coverage and allelic fraction in the normal and tumor samples. Joint circular binary segmentation is performed on the 2 signals: the log-ratio of coverages (intensities) and the log-ratio of mirrored allelic fractions (mBAFs) in tumor versus normal. Identified segments are then classified according to loss, copy-neutral loss of heterozygosity, gain, normal (no change), or undecided categories. For WES data, as the GATK variant calling workflow combines normals and tumors as a cohort, heterozygous variants in the normal sample(s) are implicitly force-called by HaplotypeCaller in the tumor sample, yielding their depth-by-allele in both normal and tumor in the output VCF. Variant calls that pass VQSR, have a 0/1 genotype in normal, have a mapping quality greater than 30, and have non-zero coverage in tumor, were input into saasCNV.

Targeted panel sequencing. Qubit fluorometry (Life Technologies) was performed on tumor and normal gDNA samples using 1 μl of each sample to quantify the concentration and mass. An Agilent Bioanalyzer DNA 12000 chip was utilized to assess the DNA integrity in the 100-bp to 17,000-bp range prior to library preparation to eliminate any low quality input material. The DNA samples that pass qualification proceed into library construction using the manufacturer’s (Thermo Fisher Scientific) Ion Torrent Ampliseq library preparation protocol. Once libraries were completed and equalized, the Ion Chef (Thermo Fisher Scientific) was used to concentrate and load the 318 PGM chip by first placing all consumables and cartridges onto the deck of the Chef per the Ion PGM Chef Kit User Guide. While allowing the reagents to reach room temperature for at least 40 minutes, the Torrent Chef and run conditions were set up for the ParThy cancer panel. Once consumables reached room temperature, the appropriate 318 chips were placed onto each of the Chef’s centrifuge buckets and the PGM configuration was run according to the manufacturer’s suggested protocol until completion. Variants from targeted panel sequencing were called using the TorrentSuite variantCaller module v4.6.0.7 using permissive settings (Somatic – Low-Stringency preset; all COSMIC variants provided as hotspot file for force calling), exported as VCFs, and loaded into a custom MySQL database in the same way as above for the Illumina workflow. All variants were manually reviewed in IGV. If both targeted panel and WES sequencing were available, concordance analysis was done to ensure the 2 assays agreed in their somatic variant calls on regions covered by both assays by design. QC statistics, multiplexing details, variant caller version, and chip type are given in Supplemental Table 10.

Data availability. Somatic variant calls are available publicly via COSMIC (COSP42872).

Statistics. To determine mutual exclusivity of genetic alternations in gene pairs, 1-sided Fisher’s exact test was performed. A P value less than 0.05 was considered significant.

Study approval. For sequencing study of cases PC11T–PC19T human studies review and approval (MREC) was as stated in Yu et al. (20). Sequencing study of samples PC01T–PC10T and PC20T–PC31T was exempted from IRB approval as all samples were anonymous and completely de-identified.