Genomic and transcriptomic profiling reveal molecular characteristics of parathyroid carcinoma

Clinical and biochemical characteristics of parathyroid carcinoma

In total, 50 thyroid tissues were collected from three groups, 12 parathyroid carcinomas, 28 parathyroid adenomas, and 10 normal parathyroid tissues, for genomic and transcriptomic profiling (Fig. 1). The detailed protocols and quality control procedures are described in the Materials and Methods section. For genomic profiling, WES of carcinoma tissues with matching blood samples was conducted. For transcriptomic profiling, RNA sequencing was conducted for all three groups, and the resulting data were used for gene and gene set-level analyses.

We first analyzed the clinical and biochemical characteristics of parathyroid tumors (Table 1). While the normal group individuals were the youngest (mean age = 38.4), the carcinoma group individuals were younger (mean age = 42.2, p = 0.0009, Mann‒Whitney test) than the adenoma group individuals (mean age = 61.3) (Supplementary Fig. 1), representing a prevalent genetic risk factor (i.e., loss of heterozygosity at the CDC73 locus31). We found a female predominance in all three groups (67–89%), without statistically significant group-specific differences. Patients with parathyroid adenoma and carcinoma exhibited higher preoperative parathyroid hormone (PTH) and serum calcium levels than normal individuals, both of which were highest in carcinoma patients. Likewise, phosphate levels were distinctive among all three groups, while they were the lowest in carcinoma patients. Clinical genetic tests (targeted sequencing of blood) identified germline mutations in CDC73 in six out of 12 carcinoma patients, showing a compatible frequency with those of previous reports (41–61.8%)10,32. Among the 12 patients with parathyroid carcinoma, three had distant metastasis, and two had local recurrence during follow-up. These findings are consistent with the known clinical characteristics and prognosis of parathyroid carcinoma1.

Table 1 Clinical characteristics of the study subjects.

Genomic profiles of parathyroid carcinoma

To investigate the profile of genomic variations in parathyroid carcinoma, we identified somatic mutations in 10 carcinoma samples (out of 12), wherein matching blood samples were available (Fig. 2a, see Methods). The number of nonsynonymous single nucleotide variants (SNVs) and insertion‒deletions (indels) ranged from 18–848, with a median of 59, which corresponded to 1.18 mutations per megabase. Except for one sample with an exceptionally high mutation count (P5), all carcinomas had relatively lower mutation counts than other cancers33,34.

Fig. 2: Genomic profiles of parathyroid carcinoma.
figure 2

a Somatic variants found in parathyroid carcinoma samples. The somatic status of each variant was determined by confirming the absence of the corresponding variant in matched normal data, and only mutated genes found in three or more samples or genes that have been reported in other PTC studies were plotted. If there was a truncating mutation in the CDC73 coding region regardless of whether it was somatic or germline, the sample was classified as CDC73Mut; otherwise, the sample was classified as CDC73WT. b Genomic positions of CDC73 mutations in the whole carcinoma cohort. All CDC73 mutations were truncating mutations, and somatic mutations were found upstream of the gene. c Allele-specific copy number status of the genomic region including CDC73. Somatic mutated copies tended to have copy number gain. d The number of variants per sample found in the CDC73Mut and CDC73WT groups. Significance between the two groups was confirmed by the Mann‒Whitney test (p = 0.0095). e Mutational signatures of parathyroid carcinoma. SBS13 was revealed to be the major signature of the CDC73Mut group, whereas SBS6 was found in the CDC73WT group. f APOBEC enrichment score of the CDC73Mut and CDC73WT groups. The score was significantly high in CDC73Mut (p = 0.0381, Mann‒Whitney test).

First, we examined the mutation patterns of CDC73. Six patients (P5, P6, P11, P22, P32, and P75) who had either germline or somatic truncating mutations (four nonsense SNVs and seven frameshift deletions) were considered CDC73-mutant (CDC73Mut). We found that the genomic loci of germline and somatic mutations were clearly separated (Fig. 2b); germline mutations were located upstream of the Ras-like domain, which is crucial for interaction with PAF1 (polymerase associated factor 1) and chromatin35, whereas somatic mutations were clustered at exon 1–2, suggesting a complete loss of transcript. Four of the six patients (P5, P6, P11, and P22) showed apparent two-hit mutations that caused biallelic inactivation of CDC73, known as Knudson’s two-hit hypothesis (i.e., one germline predisposition and one acquired somatic mutation). One patient (P32) had a somatic frameshift mutation only (CDC73 p.R52Ifs*9), but it was accompanied by a complete loss of the wild-type allele due to loss of heterozygosity (LOH), resulting in biallelic inactivation (Fig. 2c). One patient (P75), who only had a germline truncating mutation (CDC73 p.E130Gfs*11), showed borderline to early pathological classification, indicating a chance for a subclonal somatic mutation. Overall, the mutation patterns in all six patients showed truncating CDC73-directed consequent biallelic inactivation, opposing the existence of secondary hits in genes other than CDC73.

Next, we assessed the mutation profiles of four patients (P65, P66, P68, and P77) with wild-type CDC73 (CDC73WT). We found that the number of variants was lower in CDC73WT (median = 29.0) than in CDC73Mut (median = 152.5; p = 0.0095, Mann‒Whitney test) patients (Fig. 2d), indicating relatively higher genomic integrity36,37. Despite the limited sample numbers, we observed a few recurrent mutations within CDC73WT patients. Somatic missense mutations in TP53 were observed in two CDC73WT patients (TP53 p.C3F in P68 and TP53 p.H61R in P77), both of which have been previously reported in other cancers and are expected to be deleterious. Although not statistically significant (p = 0.1333, Fisher’s exact test), this may imply the confinement of TP53 mutations in CDC73-independent parathyroid carcinoma. Similarly, of the three somatic mutations in CCDC40 (coiled-coil domain containing 40), two truncation frameshift mutations (CCDC40 p.K970Nfs*51 and CCDC40 p.G987Rfs*96) were found in CDC73WT patients (P65 and P77). The CCDC40 mutation is best known as the major cause of primary ciliary dyskinesia, but its association with cancer has not yet been reported. Other mutations in KMT2D, KRAS (in-frame deletion), and FRAT2 (FRAT regulator of WNT signaling pathway 2) are potentially associated with the oncogenesis of parathyroid carcinoma; however, the evidence remains limited.

Mutational signature analysis revealed two distinct signatures for each carcinoma group: SBS13 (single base substitution signature 13) for CDC73Mut and SBS6 for CDC73WT (Fig. 2e, see Methods). SBS1 was additionally found in the CDC73Mut group, which was later confirmed to be the exclusive signature of P5 and not present in every CDC73Mut sample (Supplementary Fig. 2). SBS13, the major signature found in the CDC73Mut group, is known for its association with activated APOBEC cytidine deaminase38,39 and has been proposed as a marker for immunotherapy and targeted therapy40,41,42. Additional analysis of APOBEC enrichment also confirmed significantly high APOBEC relevance in the CDC73Mut group (Fig. 2f), which could be another clue to the different mechanisms of tumor progression in CDC73Mut and CDC73WT carcinomas.

Transcriptomic analysis of parathyroid carcinoma

Using RNA sequencing data of 49 tissues (11 carcinomas, 7 CDC73Mut and 4 CDC73WT; 28 adenomas; and 10 normal parathyroid tissues), DEGs (see Methods for criteria) were analyzed for two conditions: carcinoma vs. normal and adenoma vs. normal (Fig. 3a). We found that the overall gene expression profiles were highly conserved in adenomas, showing a strong correlation with those of normal parathyroid glands (Pearson’s r = 0.982) (Supplementary Fig. 4). In contrast, gene expression in carcinomas deviated substantially (Pearson’s r = 0.943). Therefore, the number of DEGs was larger in carcinoma tissues (n = 1,136) than in adenoma tissues (n = 33), 26 of which were differentially expressed in both. Consequently, 1,110 carcinoma- and seven adenoma-specific DEGs were identified (Fig. 3a red and green dots, respectively; see Supplementary Table 2 for a full list of DEGs).

Fig. 3: Transcriptomic analysis of parathyroid carcinoma and adenoma compared to normal tissue.
figure 3

a Two-axis DEG analysis of parathyroid carcinoma and adenoma compared to normal parathyroid tissue. padj < 0.01 and |log2-fold change | > 1 were used as cutoffs for this analysis. The X-axis represents the DEGs of adenoma to normal, and the Y-axis represents the DEGs of carcinoma to normal. Red dots represent carcinoma-specific DEGs, green dots represent adenoma-specific DEGs, and blue dots represent common DEGs of adenoma and carcinoma. bd GSEA results of all carcinomas, CDC73Mut carcinomas and CDC73WT carcinomas compared to the expression of normal parathyroid. Enrichment analysis was performed on the H collection (hallmark gene sets) of MSigDB. Each term is shown without ‘HALLMARK’.

Pathway-level analysis revealed the enrichment of many cancer hallmark pathways in carcinoma (Fig. 3b, Supplementary Table 3, Supplementary Fig. 5). In particular, pathways involved in E2F targets, G2M checkpoint, glycolysis, Myc targets, and epithelial-mesenchymal transition (EMT) were upregulated compared to normal samples (adjusted p-value < 0.01 and FDR < 0.25, see Methods). Mild upregulation of G2M and Myc target pathways was also observed in adenoma, but KRAS signaling and TNF-alpha signaling were downregulated, in contrast to carcinoma (p < 0.003). In the GO enrichment results of adenoma-DEGs, including even subtle DEGs of 1.5-fold change to normal, upregulation of genes belonging to cell growth or neuronal development was observed (Supplementary Fig. 6).

Further assessment revealed differences in pathway activation between CDC73Mut and CDC73WT (Supplementary Table 5, Supplementary Fig. 5). The upregulation of E2F targets was observed more significantly in CDC73Mut, which might suggest an activated DNA damage response against the higher mutation burden43 (Fig. 3c). While activation of E2F and Myc targets, mTORC1 and Hedgehog signaling were commonly observed, CDC73WT parathyroid carcinoma showed stronger activation of EMT and oxidative phosphorylation (Fig. 3d). EMT activation in tumor progression is well known to be significantly involved in tumorigenesis and angiogenesis. Moreover, the combination of upregulated oxidative phosphorylation can be strong proof of metastasis44,45 and, more specifically, the hybrid E/M phenotype46,47,48. Indeed, one of the patients in the CDC73WT cohort showed multiple metastases after sample preparation. Although no signs of metastasis were found in other CDC73WT patients, the possibility of metastasis cannot be ruled out considering this result.

Allelic imbalance and allele-specific expression of CDC73

As shown earlier, two-hit mutations in CDC73 result in two separable alleles: one allele that harbors germline variants (referred to as CDC73Germ) and the other that acquires somatic mutations (CDC73Som). Combined analysis of genomic and transcriptomic profiles enabled the inspection of DNA- and RNA-level imbalances, including allele-specific copy number alterations (CNAs) and expression biases between the two alleles. We detected frequent (70%, 7/10) allele-specific CNAs at 1q31.2, which included the genetic lesion of CDC73 (Supplementary Fig. 8). Notably, all four samples with CDC73 two-hit mutations (P5, P6, P11, and P22) showed elevated intratumor variant allele frequency (VAFtumor) of somatic variants, while the VAFtumor of germline variants decreased in tumor DNA, suggesting copy number gains in CDC73Som and/or losses in CDC73Germ (Fig. 4a, Supplementary Fig. 9, see Methods). Likewise, as described previously (Fig. 2c), the patient with CDC73 LOH (P32) retained CDC73Som only. Therefore, all five patients with CDC73 somatic mutations showed a relative gain in CDC73Som.

Fig. 4: Allelic imbalance and allele-specific expression of CDC73.
figure 4

a Allele frequencies of CDC73 mutations traced through normal DNA, tumor DNA, and tumor RNA and corrected to the values expected to be entirely derived from tumor cells. In the two-hit mutant group, the somatically mutated copy tended to have a higher copy number than the copy with the germline mutation. Furthermore, the expression of somatically mutated copies is commonly upregulated, and as in the case of P75, the expression of the first-hit copy is suppressed even though there is no second hit. b Hypothesis for the onset and progression of parathyroid carcinoma in terms of CDC73 mutation status. The size of the arrows indicates the transcription rate, and the light gray arrow indicates nonfunctional transcription.

Further analysis revealed additional bias at the transcriptional level. Allele-specific RNA-seq analysis (Fig. 4a, Supplementary Fig. 10, see Methods) showed 2.3- to 7.9-fold higher gene expression in CDC73Som than expected (from the DNA-level allele frequency) in the four patients with CDC73 two-hit mutations. These results indicate that the allelic gain in CDC73Som is not only retained but further intensified at the transcription level owing to allele-specific expression. Likewise, we found that the gene expression of CDC73Germ in P75, which harbored the CDC73Germ mutation only, was substantially lower than expected, also supporting the higher allele-specific expression in CDC73Som in all six CDC73Mut patients.

Based on the results, we suggest a plausible model that explains the duplex preference (genomic and transcriptomic) of CDC73Som over CDC73Germ (Fig. 4b). Born with one inactivated allele (CDC73Germ), the other intact allele (CDC73Wt) solely takes the designated role of the gene, such as cellular homeostasis and tumor suppression. This leads to a more active use of CDC73Wt before tumorigenesis, which can be achieved by either deteriorating CDC73Germ (e.g., copy loss or transcriptional suppression) or positively selecting CDC73Wt (e.g., copy gain or transcriptional activation). At the time of tumorigenesis, the second hit (somatic truncation mutations) converts CDC73Wt to CDC73Som, maintaining the genomic and transcriptomic preference over CDC73Germ, with further selective advantages acquired during tumor progression. Overall, our model explains the allelic imbalance and allele-specific expression in CDC73 based on the functional compensation for the haploinsufficiency of CDC73 caused by germline truncating mutations, as reported several times in other studies49,50.

Molecular classification of parathyroid carcinoma and adenoma

Based on the transcription profiles identified using RNA-seq, we constructed a molecular classification model for parathyroid carcinoma from adenoma and normal parathyroid glands. A total of 597 genes with strong carcinoma specificity were selected from the DEG sets by applying the most stringent filtration method (Supplementary Fig. 11g) and were used for hierarchical clustering of 49 samples (11 carcinomas, 28 adenomas, and 10 normal parathyroid tissues). In the initial clustering without any training or optimization, we found that all carcinoma samples were clearly separated from non-carcinoma samples (Fig. 5), indicating intrinsic differences in the molecular characteristics that are present in the gene set.

Fig. 5: Molecular classification of parathyroid carcinoma and adenoma.
figure 5

Hierarchical clustering result using 597 carcinoma-specific DEGs (red dots in Supplementary Fig. 11g). All carcinomas were well-divided and clustered from adenomas and normal tissues, suggesting that the feature genes are sufficient to classify carcinomas.

The clinical utility of the molecular classification was shown in two patients: P27 and P67. Both patients were initially diagnosed with atypical parathyroid neoplasm with uncertain malignant potential but were separately clustered; P27 was clustered with carcinoma, and P67 was clustered with adenoma (Fig. 5 red arrows). Further prospective follow-up identified clinical recurrence in P27, whereas no signs of pathological progression, including typical capsular or vascular invasion, were observed in P67. Further cohort-level studies are required to validate the utility of molecular classification in cases with uncertain malignant potential, which take place in 0.5–5% of parathyroid tumors1,4.

We also found that the non-carcinoma group was divided into two subgroups. Therefore, we checked whether they are biologically distinctive from each other. However, the adenomas in Group 1 and Group 2 did not show significant differences either in their transcriptome or clinicopathology (Supplementary Fig. 13, Supplementary Table 6).

WT1 as a potential marker for CDC73-mutant parathyroid carcinoma

Using whole transcriptome sequencing data, we searched for a possible single-gene marker for parathyroid carcinoma. Among the genes with carcinoma-specific expression (Fig. 4a red dots), we focused on Wilms tumor 1 (WT1) based on its functional relatedness with CDC73; WT1 is known to directly repress CDC73 and induce MYC and BCL-2 to promote cell proliferation and tumorigenesis51. In addition, WT1 has been considered a single molecular biomarker for multiple cancers due to its consistent upregulation in tumors52,53,54.

We tested the feasibility of using WT1 as a single-gene biomarker for CDC73-mutant parathyroid carcinoma. We found that the overexpression of WT1 in carcinoma was specific to CDC73Mut patients but not present in CDC73WT patients (Fig. 6a). In addition, immunohistochemical (IHC) staining of WT1 in 38 parathyroid tissues (28 adenomas, 4 CDC73WT carcinomas, and 6 CDC73Mut carcinomas) confirmed the presence of CDC73Mut–specific WT1 in parathyroid cancer tissues (Fig. 6b). That is, neither adenomas (Fig. 6c) nor CDC73WT carcinomas (Fig. 6d) were stained with the WT1 antibody, but five CDC73Mut carcinomas (out of 6, Fig. 6e) showed positive staining (Supplementary Fig. 14 for more IHC staining images). Since specific splicing alternatives of WT1 are known to be associated with certain diseases, we further checked the transcript type of WT1 with DEXSeq55, and we confirmed that the overexpressed WT1 in the CDC73Mut group was a canonical transcript (ESNT00000452863.10, Ensembl 108, Supplementary Fig. 15). We anticipate that these results will provide a basis for the future development of a clinical test for a faster, cheaper, and more accurate diagnosis of parathyroid cancer and its mutation status.

Fig. 6: Differential expression of WT1 in CDC73Mut carcinomas.
figure 6

a Groupwise comparison of WT1 expression. In CDC73Mut carcinomas, significant overexpression of WT1 was observed (p < 0.0001 in Kruskal‒Wallis test). b IHC-stained sample counts. WT1 positivity was found only in the CDC73Mut sample group, and none of the adenomas and normal samples were found to be WT1 positive. Fisher’s exact test was performed with the number of WT1-negative versus WT1-positive samples in each group. The number at the top of the bar indicates a positive ratio in each group. c WT1-negative adenoma sample (P2). d WT1-negative CDC73WT carcinoma sample (P65). e WT1-positive CDC73Mut carcinoma sample (P6).

Read more here: Source link