Front. Endocrinol.Frontiers in EndocrinologyFront. Endocrinol.1664-2392Frontiers Media S.A.10.3389/fendo.2023.1221060EndocrinologyOriginal ResearchSpecific genetic aberrations of parathyroid in Chinese patients with tertiary hyperparathyroidism using whole-exome sequencingLiLei12†ShengQixuan1†ZengHuajin1†LiWei1WangQiang1MaGuanjun1XuXinyun1QiuMing1ZhangWei1*ShanChengxiang1*1Department of Thyroid, Breast and Hernia Surgery of Changzheng Hospital Affiliated with Naval Military Medical University, Shanghai, China2Department of General Surgery, Shanghai General Hospital, Shanghai Jiao Tong University School of Medicine, Shanghai, China
Edited by: Jeff M. P. Holly, University of Bristol, United Kingdom
Reviewed by: Hao Zhang, Shanghai Sixth People’s Hospital, Shanghai Jiao Tong University, China; Pierre De Meyts, Université catholique de Louvain, Belgium
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
Background
Tertiary hyperparathyroidism (THPT) is a peculiar subtype of hyperparathyroidism that usually develops from chronic kidney disease (CKD) and persists even after kidney transplantation. Unlike its precursor, secondary hyperparathyroidism (SHPT), THPT is characterized by uncontrolled high levels of calcium in the blood, which suggests the monoclonal or oligoclonal proliferation of parathyroid cells. However, the molecular abnormalities leading to THPT have not yet been fully understood.
Methods
In this study, we analyzed DNA samples from hyperplastic parathyroid and corresponding blood cells of 11 patients with THPT using whole-exome sequencing (WES). We identified somatic single nucleotide variants (SNV) and insertions or deletions variants (INDEL) and performed driver mutation analysis, KEGG pathway, and GO functional enrichment analysis. To confirm the impact of selected driver mutated genes, we also tested their expression level in these samples using qRT-PCR.
Results
Following quality control and mutation filtering, we identified 17,401 mutations, comprising 6690 missense variants, 3078 frameshift variants, 2005 stop-gained variants, and 1630 synonymous variants. Copy number variants (CNV) analysis showed that chromosome 22 copy number deletion was frequently observed in 6 samples. Driver mutation analysis identified 179 statistically significant mutated genes, including recurrent missense mutations on TBX20, ATAD5, ZNF669, and NOX3 genes in 3 different patients. KEGG pathway analysis revealed two enriched pathways: non-homologous end-joining and cell cycle, with a sole gene, PRKDC, involved. GO analysis demonstrated significant enrichment of various cellular components and cytobiological processes associated with four genes, including GO items of positive regulation of developmental growth, protein ubiquitination, and positive regulation of the apoptotic process. Compared to blood samples, THPT samples exhibited lower expression levels of PRKDC, TBX20, ATAD5, and NOX3 genes. THPT samples with exon mutations had relatively lower expression levels of PRKDC, TBX20, and NOX3 genes compared to those without mutations, although the difference was not statistically significant.
Conclusion
This study provides a comprehensive landscape of the genetic characteristics of hyperplastic parathyroids in THPT, highlighting the involvement of multiple genes and pathways in the development and progression of this disease. The dominant mutations identified in our study depicted new insights into the pathogenesis and molecular characteristics of THPT.
tertiary hyperparathyroidismwhole-exome sequencingsomatic mutationdriver mutationcopy number variantssection-in-acceptanceMolecular and Structural EndocrinologyIntroduction
Hyperparathyroidism (HPT) is a common endocrine disorder, characterized by aberrant parathyroid gland enlargement, excessive circulating parathyroid hormone (PTH), and disturbed bone and mineral metabolism. Tertiary hyperparathyroidism (THPT), a complex and unique subtype of HPT (1, 2), arises predominantly from secondary hyperparathyroidism (SHPT) during the prolonged course of chronic kidney diseases (CKD) and continues even after kidney transplantation, while showing minimal response to therapeutic agents (3).
In contrast to SHPT, which is its precursor, THPT is characterized by uncontrolled hypercalcemia that resembles the autonomous form of HPT. This suggests monoclonal or oligoclonal proliferation of parathyroid cells in THPT, rather than a widespread hyperplastic expansion of all parathyroid cells in response to external growth stimuli. Although X-inactivation assays have already shown that monoclonal outgrowths can occur against a backdrop of supposed generalized non-clonal parathyroid hyperplasia in the context of CKD (4), the pathogenesis of THPT remains incompletely defined.
To fully understand the development of parathyroid tumors caused by clonal growth, it is important to identify the potential acquired and selected driver mutations involved in the process. In cases of nonhereditary or sporadic parathyroid adenomas (PAs), specific mutations in certain genes had been shown to be responsible for their development. Two genes, MEN1 and CCND1, an archetype tumor suppressor and a prototypical proto-oncogene, respectively, had been unequivocally established as paramount tumorigenic drivers in PAs (5, 6). MEN1 was frequently associated with somatic mutations, occurring in 12% to 35% of sporadic PAs, while rearrangements of the CCND1 locus appeared to occur in up to 8% of PAs (7). Additionally, novel ZFX mutations at residues p.R786 and p.R787 were discovered in 4.6% of a parathyroid adenoma cohort (8, 9). Moreover, confirmed somatic mutations of CDNK1B were observed in less than 1% of sporadic PAs (10, 11). Still, other genes, some involved in pathways leading to MEN1 inactivation (12), CCND1 activation (13), WNT activation (14), or the PTH dysregulation (15), had been and continued to be explored as candidate parathyroid tumor driving genes. All of these findings could provide alternative potential targets for the treatment of PAs.
Despite the vast knowledge on genetic alterations in sporadic primary hyperparathyroidism (PHPT), the underlying molecular mechanisms of THPT remained largely elusive, particularly in the Chinese population. To bridge this gap, our study employed a comprehensive approach utilizing whole-exome sequencing (WES) analysis of matched tumor and blood cell DNA samples from Chinese THPT patients. We aimed to identify potential genetic aberrations implicated in the pathogenesis of THPT and to assess the mutation frequency of established genes in the Chinese cohort.
Materials and methodsTHPT diagnostic criteria
To fulfill the diagnosis of THPT, a triad of conditions must be met according to the following established criteria (16): (i) A documented history of renal insufficiency concomitant with SHPT, either prior to or after renal transplantation; (ii) Elevated circulating levels of PTH that are biochemically abnormal, concomitant with uncontrolled hypercalcemia; (iii) Appearance of abnormal enlargement of parathyroid glands, in conjunction with pathologically confirmed nodular or adenomatous hyperplasia.
Clinical samples
The study set consisted of DNA samples isolated from pathological parathyroid hyperplasia tissue, as well as matched blood samples from each patient, following standard procedures. Agarose gel electrophoresis was employed to analyze the extent of degradation and the presence of RNA contamination. The purity of the extracted DNA was evaluated using Nanodrop, and the DNA concentration was precisely measured using Qubit. Samples over 1.5ug in contents as well as with an optical density (OD) value ranging from 1.8 to 2.0 were selected for further analysis.
The study protocol was approved by the Committee of Ethics Review Board of Changzheng Hospital affiliated to Naval Military Medical University. All patients provided informed consent for inclusion in the study according to the guidelines at our institution.
Whole-exome library construction and high-throughput sequencing
Whole-exome of the DNA samples was captured and enriched strictly deferring to the protocol of commercially available Agilent SureSelect v6 Kit on Agilent Liquid Chip Capture System. Specifically, prepped DNA fragments were hybridized with biotinylated RNA library “baits”, and then mixed with streptavidin coated magnetic beads. Beads-conjunct fractions were selectively captured by magnet and put to washing and digestion of RNA “baits” subsequently. The fractions obtained were then amplified by PCR for the following sequencing. Captured libraries were firstly testified for its quality including insert size and effective concentration by Agilent 2100 and Q-PCR respectively. Qualified libraries were sequenced on illumina Novaseq 6000 sequencing platform with reads length of 150bp. Raw image data were converted into raw sequence reads through Base Calling analysis.
Bioinformatics analysis
Sequencing data were processed for base quality distribution and GC content distribution by statistical methods. Raw reads were filtered out of reads with adapters or with low quality to get clean reads, which were mapped to the reference genome Hg 19 using BWA (http://bio-bwa.sourceforge.net/) software. Reads generated by PCR-duplication were discarded using Picard software. Analysis of the raw data from WES was performed according to GATK Best Practices recommendations (www.broadinstitute.org/gatk/). Single nucleotide variants (SNV) and insertions or deletions (INDEL) variants detection and filtration were performed using the genome analyzer Mutect2 software (17). Copy number variants (CNV) were analyzed by CNVkit (http://github.com/etal/cnvkit), and structural variants (SV) analyzed by Manta (https://github.com/Illumina/manta). The identified somatic variants were further annotated using VEP software (18). Functional impact predictions of each verified somatic variant were performed using both Polyphen2 (http://genetics.bwh.harvard.edu/pph2/) and SIFT (http://sift.jcvi.org). Driver mutations were analyzed by OncodriveCLUSTL (19). GO (Gene Ontology, http://www.geneontology.org) and KEGG (Kyoto Encyclopedia of Genes and Genomes, http://www,genome.jp/kegg/) facilitated the functional enrichment analysis of variants-related genes and were performed by KOBAS (http://kobas.cbi.pku.edu.cn/ ).
Confirmatory experiments
In the original 11 THPT samples and their corresponding normal blood samples, driver mutated genes were selected and analyzed. The mRNAs from these 11 pairs of samples were extracted for reverse transcription, and primers were designed for the target genes. The qRT-PCR was performed to detect the expression differences of six genes, including PRKDC, TBX20, ATAD5, ZNF669, NOX3, and UBR3.
Statistical analysis methods
Driver gene analysis performed by OncodriveCLUSTL generated three p-values per genomic elements (GEs), including an empirical p-value, an analytical p-value and a second analytical p-value (19). All resulting p-values were adjusted (q-values) using the Benjamini-Hochberg method and GEs below 1% false-discovery rate (FDR) are considered significant. The results shown here were based on ranking of GE score analytical p-values. The driver gene graph, CNV heatmap and Circos graphs were generated using R Studio(version 4.0.2). The result of GO and KRGG analysis were submitted to hypergeometric test and then Fisher’s exact test, and graphical presented by KOBAS. Functional enrichment with p value <0.05 is considered significant.
Clinical baseline data of patients were presented as (Mean ± SD). Statistical analyses of the result of qRT-PCR were performed and graphics were generated using GraphPad Prism (version 9.5.1). Results were presented in graphical format showing (Mean ± SEM) of three repeated experiments. Data conforming to Gaussian distribution were subjected to unpaired t-test. P value <0.05 is considered statistically significant.
ResultsClinical information of patients
Our study included 11 Chinese patients, consisting of 4 males and 7 females, who underwent parathyroidectomy for THPT between November 2020 and April 2021. All patients suffered from renal failure, with 9 patients relying on hemodialysis, one patient relying on peritoneal dialysis and one patient having undergone renal transplantation. The average age of the patients was (49.27 ± 6.83) years old, with dialysis duration ranging from 1 to 18 years. Their average preoperative serum calcium level was (2.65 ± 0.11) μmol/L (reference value 2.25-2.4 μmol/L), and average serum phosphorus level (2.01 ± 0.63) mmol/L (reference value 0.81-1.45 mmol/L). The mean preoperative PTH level in 10 uremic patients on dialysis was higher than 1000 pg/mL (reference value 15-65 pg/mL), at the level of (1577.10 ± 537.70) pg/mL, while in the post-transplantation patient, the PTH level was 380 pg/mL. The serum Creatinine level averaged at (822.09 ± 373.36) μmol/L. Nodular hyperplasia or adenomatous hyperplasia of the parathyroid gland was observed in all patients. The detailed baseline characteristics of included patients were provided in Table 1.
Clinical and molecular features of participant patients.
The sequencing accuracy was assessed based on the quality value, revealing an average Q30 of 92.90% ± 0.38%. Additionally, the GC content was found to be 48.42% ± 0.80%. A total of 73.32 (range, 56.00 to 100.50) million reads were obtained per tumor sample and 83.26(range, 62.80 to 113.60) million reads per blood sample. The tumor samples exhibited a 73.55% ± 7.35% fraction of genome with 30× coverage and an average median sequencing depth of 45.55×, while the blood samples showed an 81.27% ± 4.97% fraction of genome with 30× coverage and an average median sequencing depth of 53.82×.
Identification of tumor-specific somatic variants
According to our filter criteria, a total of 17401 mutations containing 2687 existing known variants were discovered in 11 clinical paired samples, including four most common variants types, 6690 (38.45%) missense variants, 3078 (17.69%) frameshift variants, 2005 (11.52%) stop gained variants and 1630 (9.37%) synonymous variants. A large number of mutations were evenly distributed across the chromosomes, with the exception of the sex chromosomes. The specific proportion of each type of mutation is shown in Figure 1A. The detailed distribution of Single Nucleotide Polymorphism (SNP), INDEL, and CNV in the two samples with the most and fewest mutations (4109 and 447) were visualized using Circos and displayed in Figure 1B. The distributions in the other nine samples were shown in Supplementary Figure 1. Synonymous variants, also known as silent mutations, resulting in no change to the amino acid sequence of a protein, were excluded from further analysis. CNV analysis demonstrated the copy number deletion was observed in chromosome 22 (Chr22) in 6 samples, with loci starting from 15,941,981 to 51,237,548, covering 582 genes. Besides, copy number gain was observed in Chr12 and Chr7 in 2 samples respectively, and in Chr9 in one sample (Figure 1C). The detailed CNVs detected in each sample was shown in Supplementary Figure 2.
The statistical data of mutation types and distribution on chromosomes in THPT patients. (A) The specific proportion of each type of mutation after Quality Control. Missense variants, the most common type, accounted for 38.45% of all the mutations, followed by frameshift variants (17.69%), stop gained (11.52%) and synonymous variants (9.37%). (B) Circos plots of the WES data of two patients with the most and fewest mutations. The outermost ring showed the chromosomes, with each chromosome represented as a separate segment of the circle. The second ring showed the SNP along each chromosome, with SNV types represented as colored blocks. The next ring showed Indel, with frameshift, nonframeshift and others as different colored blocks. The innermost ring showed CNV, with gain, loss and normal shown in different colors. (C) The distribution of copy number gain and loss in chromosomes of 11 samples. Except for sex chromosomes, copy number gain was observed in Chr12 in THPT tumor sample 09 and 11, and in Chr7 in THPT tumor sample 01 and 05, and in Chr9 in THPT tumor sample 02. Copy number loss was observed in Chr 22 in 6 samples, namely THPT tumor sample 02, 05,06,07 and 11.
Using genetic mutations that lead to malfunction in the coding region as a background, a number of 179 genes emerged as statistically significant during driver mutation analysis. The top 10 significant mutated genes were identified as TOPORS, TBX20, ATAD5, ZNF669, DTNA, TTN, RP1-127H14.3, SMEK2, NOX3, and UBR3 (Figure 2A). Among these mutations, TBX20, ATAD5, ZNF669, and NOX3 were found to have recurrent missense mutations in 3 different patients, respectively.
The results of driver mutation analysis and KEGG and GO analyses of the mutated gene set. (A) A number of 179 genes emerged as statistically significant during driver mutation analysis. The top 10 significant mutated genes were identified as TOPORS, TBX20, ATAD5, ZNF669, DTNA, TTN, RP1-127H14.3, SMEK2, NOX3, and UBR3. (B) Two pathways, non-homologous end-joining and cell cycle, were enriched in KEGG analysis with only one gene involved. (C) Various cellular components and cytobiological processes were enriched by GO analysis.
After conducting Mutect2 analysis and filtering the mutated genes, KEGG and GO analyses were employed to investigate the biological functions and pathways of the identified genes. However, the outcomes of KEGG analysis indicated only two pathways, non-homologous end-joining and cell cycle, which were enriched by a sole gene, PRKDC (Figure 2B). In contrast, the GO analysis demonstrated significant enrichment of various cellular components and cytobiological processes, such as positive regulation of developmental growth, protein ubiquitination, and positive regulation of apoptotic process (Figure 2C). Interestingly, the GO terms were found to be associated with four genes, namely ENTPD7, PRKDC, SPRED1, and UBR3, signifying their potential role in regulating cellular processes and functions.
The mutation data of the aforementioned genes were exhaustively compiled and summarized in Table 2. Specifically, a total of 8 mutations in PRKDC were identified, encompassing recurrent stop gained in exon 21 in 5 samples (5/8), repeated frameshift variant in exon 72 in 2 samples (2/8), and a missense variant in exon 84 in one sample (1/8). The PRKDC mutation was predicted to be one of the foremost driver mutations, ranking 15th in the list. In contrast, a repeated missense variant in exon 3 of ENTPD7 was detected in 5 samples, which was predicted to be tolerant by SIFT and benign by PolyPhen in terms of its effect on the protein sequence. Frameshift variants in 5 samples and missense variants in 2 samples were observed in SPRED1, which were not considered significant in driver gene analysis and thus were eliminated from further consideration. Regarding UBR3, the mutations comprised repeated stop gained in exon 4 in 5 samples and a different stop gained variant in exon 9 in one different sample, frameshift variants in exon 21 in two CDS positions in one sample, and splice donor variants in exon 22 and 15 in 2 samples respectively. It was deemed a significant driver gene in tumors in diver gene analysis, ranking 10th in the results.
Mutated genes with repeated variants in over 3 samples with detailed mutation information and functional analysis.
miRNA regulation of DNA damage response;Cell cycle;Retinoblastoma gene in cancer;Fas ligand pathway and stress induction of heat shock proteins;Pathways affected in adenoid cystic carcinoma;ATM signaling in development and disease;DNA IR-double strand breaks and cellular response via ATM;DNA IR-damage and cellular response via ATR;Non-homologous end joining;DNA damage response;DNA repair pathways, full network
B cell lineage commitment
T cell receptor V(D)J recombination
T cell lineage commitment
negative regulation of immunoglobulin production
ectopic germ cell programmed cell death
Frameshift
72/87
c.10087-10088gga/gAATCCTCTGAACTGga
p.G3363ESSELX
2
THPT11
Not applicable
Not applicable
nonhomologous end joining complex
regulation of smooth muscle cell proliferation
Missense variant
84/87
c.11650C>T
p.R3884W
1
THPT05
deleterious
Probably damaging
U3 snoRNA binding
maturation of 5.8S rRNA
regulation of epithelial cell proliferation
telomere capping
positive regulation of double-strand break repair via nonhomologous end joining
negative regulation of cellular senescence
response to gamma radiation
activation of innate immune response
T cell differentiation in thymus
positive regulation of erythrocyte differentiation
spleen development
small-subunit processome
intrinsic apoptotic signaling pathway in response to DNA damage
telomere maintenance
protein-DNA complex
somitogenesis
thymus development
response to activity
positive regulation of fibroblast proliferation
protein destabilization
positive regulation of type I interferon production
double-strand break repair via nonhomologous end joining
double-strand break repair
peptidyl-threonine phosphorylation
regulation of circadian rhythm
rhythmic process
regulation of hematopoietic stem cell differentiation
negative regulation of protein phosphorylation
positive regulation of translation
cellular protein modification process
double-stranded DNA binding
cellular response to insulin stimulus
nuclear chromosome, telomeric region
peptidyl-serine phosphorylation
heart development
protein kinase activity
transcription regulator complex
brain development
protein domain specific binding
cellular response to DNA damage stimulus
Based on the comprehensive analyses of above genes, PRKDC and UBR3 stood out as the significant ones that might play a role in carcinogenesis. The expression and mutation of PRKDC and UBR3 in different malignancies were primarily explored and displayed in Supplementary Figure 3. Generally, overexpression of PRKDC was observed in most cancers except for renal cancer, in which the level of PRKDC was lower than that in normal tissues (Supplementary Figure 3A). On the contrary, the differentiated expression of UBR3 were illustrated to be downregulated in a series of cancers, while merely upregulated in esophageal carcinoma (ESCA) and stomach adenocarcinoma (STAD) (Supplementary Figure 3B). Besides, the mutation frequencies of PRKDC in malignancies were dramaticly high, with over 10% of whole tested in uterine corpus endometrial carcinoma (UCEC, 16.57%), skin cutaneous melanoma (SKCM,13.24%), and colon adenocarcinoma (COAD, 10.10%) (Supplementary Figure 3C). By contrast, the percentage of cancer samples with UBR3 mutation was relative less, with mutation rate at 9.98% in UCEC being the top one among all the cancers(Supplementary Figure 3D).
The mRNA expression levels of the target genes with mutations
Based on the aforementioned analyses, a subset of six genes, namely PRKDC, TBX20, ATAD5, ZNF669, NOX3, and UBR3, were selected for further investigation of their mRNA expression levels in 11 samples. Primers designed and applied were illustrated in Supplementary Table 1. The qRT-PCR results indicated a general downregulation of PRKDC, TBX20, ATAD5, and NOX3 in hyperplastic parathyroids of THPT compared to blood samples, while the expression of ZNF669 was found to be comparable to that in blood (Figure 3A). Additionally, the expression level of UBR3 was so low that qRT-PCR could not detect it. Notably, the expression levels of PRKDC, TBX20, and NOX3 in hyperplastic parathyroids of THPT with exon mutations were relatively lower compared to those without mutations, although the difference did not reach statistical significance (Figure 3B).
The mRNA expression levels of specific genes in 11 samples tested by qRT-PCR. (A) PRKDC, TBX20, ATAD5, and NOX3 were down-expressed in hyperplastic parathyroids compared to normal blood samples, while the expression of ZNF669 was comparable between hyperplastic parathyroids and normal blood samples. (B) The expression levels of PRKDC, TBX20, and NOX3 in with exon mutations were relatively lower compared to those without mutations, while the levels of ATAD5 and ZNF669 were unchanged. The difference did not reach statistical significance. *p<0.05, **p<0.01, ****p<0.0001, ns, nonsense, p>0.05.
Analysis of known mutated genes in sporadic PAs and parathyroid carcinoma
MEN1, a tumor suppressor that are confirmed frequently mutated in sporadic parathyroid adenoma, was detected with no mutation in these THPT samples. Similarly, no mutation was detected in the proto-oncogene of PAs, CCND1, either. Loss-of-function mutations of Cell Division Cycle 73 (CDC73), the tumor suppressor gene playing a role as a major driver genetic defect in the carcinogenesis of parathyroid carcinoma, was observed with a missense mutation in exon 16, predicted to be tolerated and benign in merely one sample in this study. The gene involved in PTH synthesis or regulation pathway, Calcium-sensing Receptor (CaSR), was shown mutated in the form of a novel missense variant in one sample, while it was predicted to be tolerated by SIFT and benign by PolyPhen. Besides, known mutations in PAs revealed via WES analysis, including EZH2, EZH1, ZFX, mTOR and FAT1, were compared with the sequencing data herein. Among these, different variant sites in FAT1 gene were observed in 3 samples, with existing missense variations detected in exon 2 in two distinctive CDS positions in 2 samples. Frameshift variants in exon 11 were observed in ZFX gene in 2 samples. A missense variant in exon 7 in EZH2 was detected in one sample, predicted to be tolerated by SIFT and benign by PolyPhen, and no mutation in exon in EZH1 was detected. A missense variant in exon 4 and 2 frameshift variants in exon 11 in mTOR gene were verified in 2 samples. Further, PI3K-AKT-mTOR pathway, integrating both intercellular and extracellular signals to regulate cell growth, proliferation, metabolism and survival, was also checked. Two existing missense variants in exon 3 were observed in PIK3CA gene in 2 samples. MEP3K1, involved in the Ras-Raf-MEK-ERK pathway, was demonstrated with replicated missense variants in both exon 14 and 17 in 2 different samples respectively.
Discussion
THPT is a condition characterized by the persistent release of PTH and loss of feedback regulation to blood calcium levels, resulting in hypercalcemia and elevated PTH levels, originally described in 1970s (20, 21). This condition typically arises when treatment for secondary hyperparathyroidism (SHPT) is inadequate, leading to the formation of abnormal parathyroid structures that are unresponsive to medical interventions. Certain gene mutations may contribute to the development and progression of THPT, yet research in this area is limited. In light of this knowledge gap, we sought to conduct a comprehensive investigation of the genetic features of hyperplastic parathyroids in THPT, with the aim of enhancing our understanding of the underlying pathogenesis and molecular characteristics of this condition.
In this study, we examined abnormal parathyroid samples from patients with THPT and found that all displayed adenomatous or nodular hyperplasia as the pathological features. Using WES and CNV analysis, we compared the exome mutations of these hyperplastic parathyroids to corresponding blood samples. After quality control and mutation filtering, we identified a total of 17,401 mutations, with missense variants being the most common, followed by frameshift variants, stop-gained variants, and synonymous variants. Missense, frameshift, and stop-gained variants in genes can alter gene expression and lead to loss or alteration of protein function. Furthermore, recurrent mutations in hyperplastic parathyroids may contribute to the pathogenesis and development of THPT. Following a comprehensive analysis and verification, our study has uncovered three significant genes with variants that are the candidates for driving the pathogenesis of THPT, including PRKDC, TBX20, and NOX3.
Prior researches have endeavored to explicate the molecular characteristics of PHPT stemming from parathyroid adenoma and carcinoma. A prevalent genomic aberration observed in parathyroid adenoma is the loss of heterozygosity of the MEN1 tumor suppressor gene, which governs cell growth and cell cycle. In one study, the incidence of MEN1 mutations was found to reach up to 35% in parathyroid adenoma (22). CCDN1, an oncogene that encodes Cyclin D, has been unequivocally established as a pivotal contributor to the pathogenesis of sporadic parathyroid adenoma and carcinoma (23, 24). Conversely, the inactivation of the tumor suppressor gene CDC73 has been found to be a crucial determinant in the carcinogenesis of parathyroid carcinoma, with a prevalence of up to 70% in sporadic cases (23). In the present study, no recurrent mutations of the aforementioned three genes were detected in hyperplastic parathyroids of THPT, indicating a distinctive driving mechanism of adenomatous hyperplasia of parathyroids in TPHP as compared to PHPT. Additionally, recent research has uncovered somatic mutations in several other genes at a relatively high frequency using WES, including EZH2, EZH1, ZFX, CaSR, and FAT1 (25, 26). While the samples analyzed in this study exhibited only one or two dependent variants with no discernible patterns. The incongruity of gene mutations in hyperplastic parathyroids of THPT and sporadic parathyroid adenoma and parathyroid carcinoma reveals the fundamental disparities of pathogenetic driving genes in THPT, also accounting for the ineffectiveness of medical therapies in treating THPT, despite their definite efficacy in treating PHPT and SHPT. It is noteworthy, however, that mutations in mTOR and PIK3CA genes were identified in two TPHT samples. The overactivity of the PI3K-ATK-mTOR pathway has been confirmed to be involved in reducing apoptosis and promoting proliferation in various cancers, including parathyroid carcinoma (27). The presence of mutations in the mTOR and PIK3CA genes suggests that this renowned intracellular pathway, which regulates the cell cycle, may also contribute to the development or progression of parathyroid adenomatous hyperplasia in TPHT.
To be specific, the mutation frequency of TBX20, NOX3, and PRKDC was equal to or above 3/11. Recurrent missense variants in exon 6 of the TBX20 gene were predicted to be deleterious and damaging by SIFT and Polyphen. We verified that these mutations downregulate the expression of TBX20 mRNA primarily using qRT-PCR. TBX20 is an essential transcription factor in the process of heart development, participating in the regulation of multiple signaling pathways linked to the differentiation and proliferation of cardiac cells. Mutations in this gene have been associated with diverse cardiac defective diseases (28, 29). Among skin cutaneous melanoma cases, TBX20 was detected with missense and splice mutations at a frequency of 5.63% (25 cases), which is the highest alteration frequency among the PanCancer Atlas based on data obtained from the TCGA database. In our study, we found that the expression of TBX20 was positive in hyperplastic parathyroids, although at a lower level than that in blood cells. Samples with mutations expressed relatively lower levels of TBX20, indicating a damaging and negative effect of the mutations on the expression of this gene.
NOX3 is a gene that encodes a member of the NOX family of NADPH oxidases, which play a role in generating superoxide and reactive oxygen species (ROS) and are involved in the oxidative stress response pathway. The gene is generally expressed at low levels in various tissues, except in the inner ear where it is related to the biogenesis of otoconia. In chronic kidney disease, ROS and oxidase stress response were demonstrated to inhibit PTH receptor signaling and trafficking, contributing to the acceleration of osteoporosis (30, 31). Another study probed into the role of PTH in endothelial dysfunction and found that PTH increased ROS and mediated oxidative stress to impair endothelial angiogenic competence in vitro experiments (32). NOX3 was found to be dominantly lowly expressed in the hyperplastic parathyroid of THPT compared with blood samples, and samples with missense mutations of the gene expressed relatively lower mRNA levels than those without mutations. The role of NOX3 in generating ROS might explain the remarkably high expression level in blood samples in THPT. Mutations in NOX3 were investigated to be involved in noise-related hearing loss via a Genome-Wide Association Study (33), and the abnormalities in cerebellar structure and function resulting from neuronal precursor cell proliferation induced by ROS in animal models (34). Further investigation is required to determine whether the deleterious mutation of NOX3 in exon 13 plays a role in the hyperplasia of parathyroid and excessive synthesis and secretion of PTH. Additionally, missense mutations in the other two genes, ATAD5 and ZNF669, exerted no effects on gene expression levels, indicating that they were inactive mutations in the exons.
Further, KEGG pathway and GO analyses were performed on a filtered gene set to identify significant enriched pathways and GO terms that shed light on the molecular mechanism of THPT. The results showed that the PRKDC gene played a significant role in various biological processes of the disease. Located on chromosome 8q11.21, the PRKDC gene encodes the DNA-dependent protein kinase catalytic subunit (DNA-PKcs) and is involved in DNA double-strand breaks (DSBs) repair and the maintenance of chromosome stability. Mutations in PRKDC and defects in DNA-PKcs expression have been shown to enhance the process of cell apoptosis due to the generation of DNA DSBs, leading to mutation accumulation and tumorigenesis. Based on TCGA, PRKDC had a high mutation rate in several tumors, including UCEC, SKCM, and COAD, while its expression levels varied in different cancers, as was shown in Supplementary Figure 2 (35). In colorectal cancer, the mutations in PRKDC was deemed to be one of the primary founder mutations, which leads to increased mutation load, as well as increased tumor heterogeneity (36). In hepatocellular carcinoma, a high mutation frequency (11.5%) was observed in a study cohort, and the matched drugs of DDR-mutant patients were significantly more than those of wild-type patients (37). In this study, recurrent Loss of Function mutations were confirmed in PRKDC in five out of eleven samples, indicating a high mutation frequency in THPT. Furthermore, the expression levels of PRKDC were relatively downregulated in THPT compared to blood samples. Tumor samples with stop-gain mutations in Exon 21 of PRKDC had a lower level in this gene, although the statistical significance was nonsense due to the small sample size and interclass numerical differences. These findings suggest a negative mutation in exon 21 of PRKDC in tumor samples of THPT, which may contribute to the uncontrolled dysplasia of parathyroid by affecting DNA damage repair and cell cycle.
UBR3, ubiquitin protein ligase E3 component n-recognin 3, has been shown to regulate a variety of biological processes through ubiquitination, playing a role in several sensory pathways (38), as well as genetic stability (39). Though the mutation of UBR3 has been detected in several malignancies revealed by the data from TCGA database, few studies delve into the specific role of UBR3 in tumors. Herein, we detected repeated stop gained in exon 4 in 5 samples in UBR3 gene, while the expression level of this gene was too minor to be tested by qRT-PCR. Despite of the predicted role of being the tumor driver gene, the mutation of UBR3 might not contribute dominantly to the pathogenesis of THPT.
What’s more, CNV is a type of genomic variation that alters the copy number of DNA segments in the genome, and it can affect various biological processes, including development, drug response, and disease susceptibility. In this study, CNV analysis revealed frequent copy number deletions in Chr22 among the samples (6/11), providing insights into the susceptibility of chronic kidney dysfunction patients to deteriorating into THPT. Specifically, significant genes located in Chr22 play a role in negative regulation of cell growth or proliferation and act as tumor suppressor genes, such as NF2 (40), CHEK2 (41), and DEPDC5 (42). Deficiencies in these genes, related to CNVs, may contribute to abnormal proliferation of the parathyroid and THPT and warrant further exploration.
In summary, this study presents the genetic characteristics landscape of hyperplastic parathyroids by comparing the genetic profiles of pathological parathyroids of THPT to normal cells. Our findings reveal several distinct genetic underpinnings in THPT, providing novel insights into the pathogenesis and molecular characteristics of THPT. However, in order to fully comprehend the molecular and pathway involved, further functional studies will be necessary in the future.
Limitations
The findings of this research reveal the presence of several limitations that must be acknowledged. First and foremost, the study’s sample size was limited, which restricted the statistical significance of the results. Additionally, the exclusion of functional studies from the analysis represents a unignored demerit, as such studies in vitro and in vivo are essential for verifying the roles of the studied variants. By recognizing and addressing these demerits, future research can build upon these findings to arrive at a more comprehensive understanding of the subject.
Data availability statement
The data presented in the study are deposited in the NCBI repository, accession number PRJNA988842.
Ethics statement
The studies involving human participants were reviewed and approved by Ethics Review Committee of Naval Medical University. The patients/participants provided their written informed consent to participate in this study.
Author contributions
Research design: CS; Experiments performing and analysis: LL, QS; Data collection and analysis: QS and HZ; Initial manuscript writing: CS, LL, QW; Manuscript revision: GM, MQ, WL; Final approval of manuscript: CS, WZ. All authors contributed to the article and approved the submitted version.
Funding
The study was funded by the Pyramid Talented Personnel Project of Changzheng Hospital through an award to CS. The funding party played no role in the design of the study, the collection, analysis, and interpretation of the data or the writing of the manuscript. The study was supported by Innovative Clinical Research Project of Changzheng Hospital, 2020YLCYJ-Y05.
Acknowledgments
We appreciate the technological guidance and help of Origingene during the Whole-exome sequencing.
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fendo.2023.1221060/full#supplementary-material
ReferencesAlfieriCMattinzoliDMessaP. Tertiary and postrenal transplantation hyperparathyroidism. Endocrinol Metab Clin North Am (2021) 50(4):649–62. doi: 10.1016/j.ecl.2021.08.004MessaPAlfieriCM. Secondary and tertiary hyperparathyroidism. Front Horm Res (2019) 51:91–108. doi: 10.1159/000491041DreamSChenHLindemanB. Tertiary hyperparathyroidism: why the delay? Ann Surg (2021) 273(3):e120–2. doi: 10.1097/SLA.0000000000004069ArnoldABrownMFUreñaPGazRDSarfatiEDrüekeTB. Monoclonality of parathyroid tumors in chronic renal failure and in primary parathyroid hyperplasia. J Clin Invest (1995) 95(5):2047–53. doi: 10.1172/JCI117890CromerMKStarkerLFChoiMUdelsmanRNelson-WilliamsCLiftonRP. Identification of somatic mutations in parathyroid tumors using whole-exome sequencing. J Clin Endocrinol Metab (2012) 97(9):E1774–81. doi: 10.1210/jc.2012-1743NeweyPJNesbitMARimmerAJAttarMHeadRTChristiePT. Whole-exome sequencing studies of nonhereditary (sporadic) parathyroid adenomas. J Clin Endocrinol Metab (2012) 97(10):E1995–2005. doi: 10.1210/jc.2012-2303YiYNowakNJPacchiaALMorrisonC. Chromosome 11 genomic changes in parathyroid adenoma and hyperplasia: array CGH, FISH, and tissue microarrays. Genes Chromosomes Cancer (2008) 47(8):639–48. doi: 10.1002/gcc.20565SoongCPArnoldA. Recurrent ZFX mutations in human sporadic parathyroid adenomas. Oncoscience (2014) 1(5):360–6. doi: 10.18632/oncoscience.116ArnoldAStauntonCEKimHGGazRDKronenbergHM. Monoclonality and abnormal parathyroid hormone genes in parathyroid adenomas. N Engl J Med (1988) 318(11):658–62. doi: 10.1056/NEJM198803173181102AlvelosMIVinagreJFonsecaEBarbosaETeixeira-GomesJSobrinho-SimõesM. MEN1 intragenic deletions may represent the most prevalent somatic event in sporadic primary hyperparathyroidism. Eur J Endocrinol (2013) 168(2):119–28. doi: 10.1530/EJE-12-0327Costa-GudaJMarinoniIMolatoreSPellegataNSArnoldA. Somatic mutation and germline sequence abnormalities in CDKN1B, encoding p27Kip1, in sporadic parathyroid adenomas. J Clin Endocrinol Metab (2011) 96(4):E701–6. doi: 10.1210/jc.2010-1338Romero ArenasMAFowlerRGSan LucasFAShenJRichTAGrubbsEG. Preliminary whole-exome sequencing reveals mutations that imply common tumorigenicity pathways in multiple endocrine neoplasia type 1 patients. Surgery (2014) 156(6):1351–7. doi: 10.1016/j.surg.2014.08.073TominagaYTsuzukiTUchidaKHabaTOtsukaSIchimoriT. Expression of PRAD1/cyclin D1, retinoblastoma gene products, and Ki67 in parathyroid hyperplasia caused by chronic renal failure versus primary adenoma. Kidney Int (1999) 55(4):1375–83. doi: 10.1046/j.1523-1755.1999.00396.xBjörklundPLindbergDAkerströmGWestinG. Stabilizing mutation of CTNNB1/beta-catenin and protein accumulation analyzed in a large series of parathyroid tumors of Swedish patients. Mol Cancer (2008) 7:53. doi: 10.1186/1476-4598-7-53VarshneySBhadadaSKSaikiaUNSachdevaNBeheraAAryaAK. Simultaneous expression analysis of vitamin D receptor, calcium-sensing receptor, cyclin D1, and PTH in symptomatic primary hyperparathyroidism in Asian Indians. Eur J Endocrinol (2013) 169(1):109–16. doi: 10.1530/EJE-13-0085TangJAFriedmanJHwangMSSalapatasAMBonzelaarLB. Parathyroidectomy for tertiary hyperparathyroidism: A systematic review. Am J Otolaryngol (2017) 38(5):630–5. doi: 10.1016/j.amjoto.2017.06.009CibulskisKLawrenceMSCarterSLSivachenkoAJaffeDSougnezC. Sensitive detection of somatic point mutations in impure and heterogeneous cancer samples. Nat Biotechnol (2013) 31(3):213–9. doi: 10.1038/nbt.2514McLarenWGilLHuntSERiatHSRitchieGRThormannA. The ensembl variant effect predictor. Genome Biol (2016) 17(1):122. doi: 10.1186/s13059-016-0974-4Arnedo-PacCMularoniLMuiñosFGonzalez-PerezALopez-BigasN. OncodriveCLUSTL: a sequence-based clustering method to identify cancer drivers. Bioinformatics (2019) 35(22):4788–90. doi: 10.1093/bioinformatics/btz501LudwigGDKyleCGde BlancoM. Tertiary” hyperparathyroidism induced by osteomalacia resulting from phosphorus depletion. Am J Med (1967) 43(1):136–40. doi: 10.1016/0002-9343(67)90155-6BayVKuhlencordtFSchneiderCSeemannNZukschwerdtL. 26. New diagnostical, physiopathological and technical-surgical problems in the surgery of parathyroid glands. A contribution to tertiary hyperparathyroidism. Langenbecks Arch Chir (1967) 319:202–6. doi: 10.1007/BF02659257Costa-GudaJArnoldA. Genetic and epigenetic changes in sporadic endocrine tumors: parathyroid tumors. Mol Cell Endocrinol (2014) 386(1-2):46–54. doi: 10.1016/j.mce.2013.09.005MariniFGiustiFPalminiGPerigliGSantoroRBrandiML. Genetics and epigenetics of parathyroid carcinoma. Front Endocrinol (Lausanne) (2022) 13:834362. doi: 10.3389/fendo.2022.834362BrewerKCosta-GudaJArnoldA. Molecular genetic insights into sporadic primary hyperparathyroidism. Endocr Relat Cancer (2019) 26(2):R53–72. doi: 10.1530/ERC-18-0304HuYZhangXWangOCuiMLiXWangM. Integrated whole-exome and transcriptome sequencing of sporadic parathyroid adenoma. Front Endocrinol (Lausanne) (2021) 12:631680. doi: 10.3389/fendo.2021.631680TaoXXuTLinXXuSFanYGuoB. Genomic profiling reveals the variant landscape of sporadic parathyroid adenomas in Chinese population. J Clin Endocrinol Metab (2023) 108(7):1768–75. doi: 10.1210/clinem/dgad002PandyaCUzilovAVBellizziJLauCYMoeASStrahlM. Genomic profiling reveals mutational landscape in parathyroid carcinomas. JCI Insight (2017) 2(6):e92061. doi: 10.1172/jci.insight.92061DasSMondalADeyCChakrabortySBhowmikRKarmakarS. ER stress induces upregulation of transcription factor Tbx20 and downstream Bmp2 signaling to promote cardiomyocyte survival. J Biol Chem (2023) 299(4):103031. doi: 10.1016/j.jbc.2023.103031ChenYXiaoDZhangLCaiCLLiBYLiuY. The role of tbx20 in cardiovascular development and function. Front Cell Dev Biol (2021) 9:638542. doi: 10.3389/fcell.2021.638542LiXGarciaJLuJIrianaSKalajzicIRoweD. Roles of parathyroid hormone (PTH) receptor and reactive oxygen species in hyperlipidemia-induced PTH resistance in preosteoblasts. J Cell Biochem (2014) 115(1):179–88. doi: 10.1002/jcb.24648ArduraJAAlonsoVEsbritPFriedmanPA. Oxidation inhibits PTH receptor signaling and trafficking. Biochem Biophys Res Commun (2017) 482(4):1019–24. doi: 10.1016/j.bbrc.2016.11.150GambardellaJDe RosaMSorrientoDPreveteNFiordelisiACiccarelliM. Parathyroid hormone causes endothelial dysfunction by inducing mitochondrial ROS and specific oxidative signal transduction modifications. Oxid Med Cell Longev 2018 (2018) p:9582319. doi: 10.1155/2018/9582319LavinskyJCrowALPanCWangJAaronKAHoMK. Genome-wide association study identifies nox3 as a critical gene for susceptibility to noise-induced hearing loss. PloS Genet (2015) 11(4):e1005094. doi: 10.1371/journal.pgen.1005094MazzonettoPCArizaCBOcanhaSGde SouzaTAKoGMMenckCFM. Mutation in NADPH oxidase 3 (NOX3) impairs SHH signaling and increases cerebellar neural stem/progenitor cell proliferation. Biochim Biophys Acta Mol Basis Dis (2019) 1865(6):1502–15. doi: 10.1016/j.bbadis.2019.02.022CeramiEGaoJDogrusozUGrossBESumerSOAksoyBA. The cBio cancer genomics portal: an open platform for exploring multidimensional cancer genomics data. Cancer Discov (2012) 2(5):401–4. doi: 10.1158/2159-8290.CD-12-0095PalinkasHLPongorLBalajtiMNagyANagyKBekesiA. Primary founder mutations in the PRKDC gene increase tumor mutation load in colorectal cancer. Int J Mol Sci (2022) 23(2):633. doi: 10.3390/ijms23020633YangZLiuJXueFZhangLXueHWuY. Genomic landscape of Chinese patients with hepatocellular carcinoma using next-generation sequencing and its association with the prognosis. Ann Hepatol (2023) 28(2):100898. doi: 10.1016/j.aohep.2023.100898TasakiTSohrRXiaZHellwegRHörtnaglHVarshavskyA. Biochemical and genetic studies of UBR3, a ubiquitin ligase with a function in olfactory and other sensory systems. J Biol Chem (2007) 282(25):18510–20. doi: 10.1074/jbc.M701894200MeisenbergCTaitPSDianovaIIWrightKEdelmannMJTernetteN. Ubiquitin ligase UBR3 regulates cellular levels of the essential DNA repair protein APE1 and is required for genome stability. Nucleic Acids Res (2012) 40(2):701–11. doi: 10.1093/nar/gkr744WuJMinikesAMGaoMBianHLiYStockwellBR. Intercellular interaction dictates cancer cell ferroptosis via NF2-YAP signalling. Nature (2019) 572(7769):402–6. doi: 10.1038/s41586-019-1426-6StolarovaLKleiblovaPJanatovaMSoukupovaJZemankovaPMacurekL. CHEK2 germline variants in cancer predisposition: stalemate rather than checkmate. Cells (2020) 9(12):2675. doi: 10.3390/cells9122675PadiSKRSinghNBearssJJOliveVSongJHCardó-VilaM. Phosphorylation of DEPDC5, a component of the GATOR1 complex, releases inhibition of mTORC1 and promotes tumor growth. Proc Natl Acad Sci U.S.A. (2019) 116(41):20505–10. doi: 10.1073/pnas.1904774116