Comparative Gene Expression Profiles in Parathyroid Adenoma and Normal Parathyroid Tissue

Parathyroid adenoma is the main cause of primary hyperparathyroidism, which is characterized by enlarged parathyroid glands and excessive parathyroid hormone secretion. Here, we performed transcriptome analysis, comparing parathyroid adenomas with normal parathyroid gland tissue. RNA extracted from ten parathyroid adenoma and five normal parathyroid samples was sequenced, and differentially expressed genes (DEGs) were identified using strict cut-off criteria. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were performed using DEGs as the input, and protein-protein interaction (PPI) networks were constructed using Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) and visualized in Cytoscape. Among DEGs identified in parathyroid adenomas (n = 247; 45 up-regulated, 202 down-regulated), the top five GO terms for up-regulated genes were nucleoplasm, nucleus, transcription DNA-template, regulation of mRNA processing, and nucleic acid binding, while those for down-regulated genes were extracellular exosome, membrane endoplasmic reticulum (ER), membrane, ER, and melanosome. KEGG enrichment analysis revealed significant enrichment of five pathways: protein processing in ER, protein export, RNA transport, glycosylphosphatidylinositol-anchor biosynthesis, and pyrimidine metabolism. Further, PPI network analysis identified a densely connected sub-module, comprising eight hub molecules: SPCS2, RPL23, RPL26, RPN1, SEC11C, SEC11A, RPS25, and SEC61G. These findings may be helpful in further analysis of the mechanisms underlying parathyroid adenoma development.


Introduction
Primary hyperparathyroidism (PHPT) is a common endocrine disease characterized by inappropriate excessive secretion of parathyroid hormone (PTH) and consequent hypercalcemia [1,2]. The prevalence of PHPT is estimated to be 21.1-65.6 per 100,000 person-years [3][4][5][6]. PHPT is more common in women than men, particularly those ≥ 45 years old; the ratio of incidence in women to men is almost 2:1 [3], and PHPT incidence increases rapidly with age. Long-standing PHPT causes serious complications, including urolithiasis, fracture, deterioration of renal function, neuropsychiatric complications, and gastrointestinal symptoms, including constipation and peptic ulcer [7]. Since the introduction of routine serum calcium measurement, health check-ups can identify individuals with completely asymptomatic PHPT.
The most common cause of PHPT is sporadic parathyroid adenoma, which accounts for 85% of PHPT. Parathyroid adenoma is a benign monoclonal tumor, usually involving a single parathyroid gland. Besides sporadic parathyroid adenomas, various hereditary types of hyperparathyroidism are caused by mutations of specific genes, including: multiple endocrine neoplasia type 1 and familial isolated hyperparathyroidism, caused by alterations in the MEN1 gene; familial hypercalcemia with hypercalciuria, attributable to CASR gene changes; and hyperparathyroidism-jaw tumor syndrome and familial isolated hyperparathyroidism, caused by HRP2 gene mutations [8][9][10].
To date, MEN1 and CCND1 (respectively, a tumor suppressor and a proto-oncogene) are the most solidly established molecular drivers of sporadic parathyroid adenoma. Loss of heterozygosity at the MEN1 locus on chromosome 11q13 is the most common genetic aberration in parathyroid adenoma, occurring in 26.2-50.0% of cases [11][12][13][14]. In addition to MEN1 and CCND1, other genes that participate in the development of parathyroid adenoma include those encoding cyclin-dependent kinase inhibitors, along with CTNNB1, EZH2, ZFX, GCM2, and CASR [15]. Moreover, recent whole-exome sequencing studies of parathyroid adenoma samples revealed mutations in POT1 and RAP1B [11,16]. Additional genes important for parathyroid tumorigenesis are likely to be identified by next-generation sequence analysis; however, data from transcriptome analysis, to compare gene expression in parathyroid adenomas with that in healthy parathyroid glands, are lacking.
In this study, we performed transcriptome analysis to compare parathyroid adenomas and normal parathyroid glands, with the aim of identifying differentially expressed genes (DEGs). We subsequently used the Database for Annotation, Visualization and Integrated Discovery (DAVID), the Kyoto Encyclopedia of Genes and Genomes (KEGG), and Gene Ontology (GO) databases to systematically extract meaningful biological information from the DEGs identified. We also constructed a protein-protein interaction (PPI) network to ascertain functional relationships between proteins in parathyroid adenoma.

Sample Acquisition
All patients with PHPT had a single parathyroid adenoma, and PHPT was diagnosed based on serum calcium and intact PTH levels, in the absence of other possible causes of hypercalcemia. Ultrasound and Sestamibi scans were performed routinely in all patients to localize the tumors and surgical pathology was used to determine whether or not they were parathyroid adenomas. All the PHPT patients were female and the mean age was 53 ± 4.8. The mean levels of serum calcium and intact PTH were 12.7 ± 1.0 mg/dL and 460.4 ± 492.9 mmol/L, respectively. Serum magnesium level was not routinely checked. The mean size of a parathyroid adenoma was 2.3 ± 1.1 cm. None of the patients had evidence of familial disease, a history of neck irradiation, renal insufficiency, or cardiovascular disorders. Fresh frozen parathyroid adenoma tissue was obtained during parathyroidectomy at Seoul National University Hospital, between July 2014 and September 2016.
Normal parathyroid glands were obtained from thyroid carcinoma patients who underwent thyroidectomy at Seoul National University-Seoul Metropolitan Government (SNU-SMG) Boramae Medical Center, between July 2016 and January 2017. All the patients who underwent thyroidectomy were female and the mean age was 50 ± 7.4. The mean levels of serum calcium and intact PTH were 8.7 ± 0.5 mg/dL and 34.7 ± 17.0 mmol/L, respectively. The normal parathyroid tissues were collected after intraoperative frozen biopsy for microscopic tissue confirmation, due to macroscopic ambiguity, as previously described in the literature [17]. Tissue samples were snap frozen in liquid nitrogen and stored at −80 • C. This study was approved by the institutional review board of SNU-SMG Boramae Medical Center (L-2018-311).

RNA Sequencing
Total RNA samples were extracted using an RNeasy mini kit (Qiagen, Hilden, Germany), according to the manufacturer's recommendations. RNA quality was assessed by analysis of rRNA band integrity, using an Agilent RNA 6000 Nano kit (Agilent Technologies, Santa Clara, CA, USA). RNA libraries were constructed using a TruSeq RNA Access Library Prep kit (Illumina, San Diego, CA, USA), according to the manufacturer's protocol. The size and quality of libraries were evaluated by electrophoresis using an Agilent High Sensitivity DNA kit (Agilent Technologies); fragments were 350-450 bp. Subsequently, libraries were sequenced using an Illumina HiSeq2500 sequencer (Illumina). Sequencing data were deposited in NCBI's Sequence Read Archive (SRA) and are accessible through BioProject accession number PRJNA516535.

Identification of DEGs
Quality control of FASTQ data was performed using FastQC [18]. Paired-end reads (101 bp) were aligned to the human genome (GRCh37/hg19 assembly) using STAR (v. 2.4.1d, https://github. com/alexdobin/STAR) mapper [19]. Gene expression levels were then calculated as fragments per kilobase per million reads (FPKMs) using cufflinks [20], with LNCipedia and refFlat databases for region annotations [21,22]. All FPKM values were assigned a pseudocount of 1, then transformed to the log2 scale. Gene data were filtered out if median log transformed expression was <3 in both tumor and control samples. DEGs were then defined as genes with fold-change ≥2 and false discovery rate-adjusted p values < 0.001.

Functional and KEGG Pathway Enrichment Analysis
Functional analyses were performed using DAVID (http://david.abcc.ncifcrf.gov/), a functional annotation tool that allows investigators to unravel the biological meaning behind an input gene list [23]. Based on the extracted DEGs, GO associations were analyzed in three categories: biological process (BP), cellular component (CC), and molecular function (MF) [24]. The KEGG pathway database was used to identify biological pathways enriched for the identified DEGs. Statistical significance was evaluated using the Fisher exact test [25], and p < 0.05 was regarded as significant.

Protein-Protein Interaction (PPI) Network Construction and Module Analysis
To ascertain functional interactions between proteins in parathyroid cells and extract functional modules, PPI networks were constructed using the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) database, based on identified DEGs (medium confidence score < 0.5) [26]. For further pruning and identification of a core module from the constructed PPI network, the MCODE algorithm, a recursive vertex weighting scoring scheme, that measures local network density, assigning higher scores to nodes whose immediate neighbors are more interconnected, was used [27]. Default pruning thresholds (degree cut-off = 2, node score cut-off = 0.2, k-core = 2, and max depth = 100) were implemented to identify a core module. A detailed scoring formula had been described previously [27]. Finally, the constructed PPI network and identified core modules were visualized using Cytoscape, a tool for network data integration and visualization [28].

Characteristics of Study Subjects and Mutation Screening
Ten parathyroid adenomas and five normal parathyroid tissue samples were analyzed. The characteristics of the study subjects are presented in Table 1. Somatic mutations in genes associated with hyperparathyroidism, including MEN1, CASR, AP2S1, GNA11, CDC73, CDKN1A, CDKN1B, CDKN2C, RET, PTH, CCND1, AIP, CTNNB1, EZH2, ZFX, CDC73, and FGF23, were not detected in either parathyroid adenoma or normal parathyroid tissues.

Genes Differentially Expressed between PHPT and Normal Parathyroid Tissue
The mRNA expression values for genes evaluated in the ten parathyroid adenoma and five normal parathyroid tissue samples are shown in Supplementary Table S1. A total of 247 DEGs (45 up-regulated and 202 down-regulated) were identified in parathyroid adenomas (Supplementary Table S2). The top DEGs in parathyroid adenoma, compared with normal parathyroid gland, are presented in Table 2.
In addition, we investigated the expression levels of selected genes that have been implicated in parathyroid function ( Table 3). Levels of Klotho (KL) and PTH mRNA were significantly lower in parathyroid adenoma than in normal parathyroid tissue, while those of CASR, FGFR1, FGFR2, and VDR did not differ significantly.

Gene Ontology (GO) Functional and KEGG Pathway Enrichment Analysis
Using the 45 genes up-regulated in parathyroid adenoma as input for analysis using DAVID, 16 GO terms that satisfied the cut-off criteria (p < 0.05) were identified (Table 4); six, two, and eight GO terms were from the BP, CC, and MF ontology categories, respectively. Nucleoplasm, nucleus, transcription DNA-template, regulation of mRNA processing, and nucleic acid binding were the most significantly enriched GO terms among up-regulated genes. Using the 202 genes down-regulated in parathyroid adenoma as input, 78 GO terms were extracted (Supplementary Table S3), of which 32, 37, and 9 were from the BP, CC, and MF ontology categories, respectively. The top five most significantly enriched of these terms in each GO category are presented in Table 5. Extracellular exosome, membrane endoplasmic reticulum (ER), membrane, ER, and melanosome were the most significantly enriched GO terms among down-regulated genes.

Protein-Protein Interaction (PPI) Network Construction and Module Analysis
The PPI network of genes up-regulated in parathyroid adenoma, consisting of 13 nodes and 8 edges (average node degree = 1.231; average local clustering coefficient = 0.00) is illustrated in Supplementary Figure S1. The PPI network constructed using down-regulated genes, which comprised 120 nodes and 230 edges (average node degree = 3.833; average local clustering coefficient = 0.229) was shown in Supplementary Figure S2. The most significant PPI sub-module was selected from sub-modules sorted by interaction score using clustering analysis ( Figure 1); it included eight hub molecules (SPCS2, RPL23, RPL26, RPN1, SEC11C, SEC11A, RPS25, SEC61G), which suggests these hub molecules might play crucial roles in development of parathyroid adenoma. Of them, RPL23 had the highest score, indicating that it was the most strongly interconnected with its neighbor proteins.

Protein-Protein Interaction (PPI) Network Construction and Module Analysis
The PPI network of genes up-regulated in parathyroid adenoma, consisting of 13 nodes and 8 edges (average node degree = 1.231; average local clustering coefficient = 0.00) is illustrated in Supplemetal Figure 1. The PPI network constructed using down-regulated genes, which comprised 120 nodes and 230 edges (average node degree = 3.833; average local clustering coefficient = 0.229) was shown in Supplemental Figure 2. The most significant PPI sub-module was selected from submodules sorted by interaction score using clustering analysis ( Figure 1); it included eight hub molecules (SPCS2, RPL23, RPL26, RPN1, SEC11C, SEC11A, RPS25, SEC61G), which suggests these hub molecules might play crucial roles in development of parathyroid adenoma. Of them, RPL23 had the highest score, indicating that it was the most strongly interconnected with its neighbor proteins.

Discussion
To better understand the genetics involved in parathyroid adenoma development, we characterized the global expression profiles in parathyroid adenomas and normal parathyroid tissues by transcriptome analysis. A total of 247 DEGs were identified by comparison of ten sporadic parathyroid adenoma and five normal parathyroid tissue samples. Among those genes up-regulated in parathyroid adenomas, some, such as MED12, KMT5A, BMP2K, and ATAD2, are implicated in cell proliferation and transcriptional regulation, supporting the notion that they potentially contribute to tumorigenesis.
The MED12 protein encoded by MED12 is a component of the CDK8 subcomplex, which is essential for CDK8 kinase activation [29]. In addition, recent studies revealed frequent mutations in MED12, exon 2, in both benign and malignant tumors [30,31]. Further, the protein encoded by KMT5A is a protein-lysine N-methyltransferase that can monomethylate Lys-20 of histone H4 (H4K20 methyltransferase) to repress transcription of various genes. This enzyme is an important regulator of the cell cycle [32], and KMT5A expression is elevated in different types of cancer tissues and cancer

Discussion
To better understand the genetics involved in parathyroid adenoma development, we characterized the global expression profiles in parathyroid adenomas and normal parathyroid tissues by transcriptome analysis. A total of 247 DEGs were identified by comparison of ten sporadic parathyroid adenoma and five normal parathyroid tissue samples. Among those genes up-regulated in parathyroid adenomas, some, such as MED12, KMT5A, BMP2K, and ATAD2, are implicated in cell proliferation and transcriptional regulation, supporting the notion that they potentially contribute to tumorigenesis.
The MED12 protein encoded by MED12 is a component of the CDK8 subcomplex, which is essential for CDK8 kinase activation [29]. In addition, recent studies revealed frequent mutations in MED12, exon 2, in both benign and malignant tumors [30,31]. Further, the protein encoded by KMT5A is a protein-lysine N-methyltransferase that can monomethylate Lys-20 of histone H4 (H4K20 methyltransferase) to repress transcription of various genes. This enzyme is an important regulator of the cell cycle [32], and KMT5A expression is elevated in different types of cancer tissues and cancer cell lines, including those of bladder cancer, non-small cell and small cell lung carcinoma, chronic myelogenous leukemia, hepatocellular carcinoma, papillary thyroid cancer, and pancreatic cancer [33,34]. A recent study demonstrated that cancer cell growth is significantly suppressed by a reduction or loss of KMT5A-mediated methylation of proliferation cell nuclear antigen (PCNA), a widely recognized cell proliferation marker of tumor progression including that of parathyroid adenoma [35], suggesting that KMT5A-dependent PCNA methylation might promote the development of parathyroid adenoma. Therefore, KMT5A should be further validated as a possible biomarker for parathyroid adenoma progression or development.
The active form of vitamin D, 1,25-dihydroxyvitamin D (1,25(OH) 2 D 3 ), is a potent regulator of parathyroid proliferation. A recent study revealed that 1,25(OH) 2 D 3 treatment induces miR-1228, a BMP2K targeting factor, in a dose-dependent manner in human osteoblasts [36]. Though there is still no evidence for the parathyroid gland, vitamin D is presumed to interact with the BMP2 pathway during parathyroid proliferation or growth. In addition, a large body of evidence shows that 1,25(OH) 2 D 3 alters the cell cycle, ultimately affecting pRB proteins and the E2F family of transcription factors [37]. ATAD2, an E2F target gene, binds to the MYC oncogene and stimulates its transcriptional activity [38]. We speculate that ATAD2 could link E2F and MYC pathways and contribute to parathyroid tumor development.
The KL and PTH genes, which play significant roles in calcium regulation, were down-regulated in parathyroid adenoma compared with normal parathyroid tissue samples. KL encodes the Klotho protein, which is a type I transmembrane protein, and FGF23 co-receptor. Our results are consistent with previously published data, demonstrating decreased expression of Klotho mRNA and protein in parathyroid adenomas [39,40]. KL is strongly expressed in tissues that require abundant calcium transport, including distal convoluted tubule cells in kidneys, choroid plexus, and parathyroid glands [41]. As deletion of KL functions in parathyroid gland hyperplasia, this gene may also contribute to the development of parathyroid adenoma [42]. By contrast, FGF23 a Klotho ligand, was not identified as a DEG in this study.
Interestingly, the PTH gene exhibited reduced expression, despite the higher serum PTH levels in patients with parathyroid adenomas. This suggests that parathyroid adenoma cells are less efficient at producing PTH relative to normal parathyroid cells; hence the higher levels of serum PTH in individuals with parathyroid adenoma may simply reflect that these tumors contain more cells that produce PTH than normal parathyroid tissue. This is consistent with previous studies, which reported that parathyroid adenomas have weaker PTH mRNA intensity relative to normal parathyroid tissue, and that disease severity is dependent on tumor weight [43,44].
Proteins function through interaction with other proteins; thus investigation of PPI networks is essential for understanding biological processes. Hub nodes, identified by computational measurement of their degree of relationship to neighbor nodes, are considered to have pivotal roles in networks. Hub molecules identified in this study were associated with signal peptidase complex subunit (SPCS2), ribosomal proteins (RPL23, RPL26, RPN1, RPS25), and the ER membrane (SEC11C, SEC11A, SEC61G).
PTH is secreted by exocytosis via secretory vesicles. Most secretory proteins are proteolytically processed by signal peptidases, to remove their signal peptides, and proper secretory protein cleavage is critical for the efficient function of many protein-secreting cells [45,46]. PreproPTH is synthesized as a 115-amino acid precursor and converted to pro-PTH upon cleavage of its signal peptide within the ER. The results of our PPI analysis suggest the importance of the proteolytic process of signal peptide cleavage in the ER in the development of parathyroid adenoma. Moreover, GO analysis also suggested that down-regulated genes in parathyroid adenoma are implicated in signal peptide processing, based on BP ontology terms. Notably, PTH was associated with KL and GJA1 (connexin 43) proteins in PPI analysis.
This study has limitations. First, the analysis included a relatively small number of samples; hence there is potential for the results to have been greatly influenced by each sample or individual patient characteristics. To minimize the effects of sex and age, we selected both parathyroid adenoma and normal parathyroid samples from middle-aged female patients. PHPT is observed predominantly in women, although the reason for this is still unclear. It would be interesting to check the influence of gender using a linear regression-based statistical model in a further study with a larger sample size [47]. Moreover, to rule out the effects of critical genes in the development of parathyroid adenoma, we tested for somatic mutations in all genes reported as associated with hyperparathyroidism. The second limitation of this study was the lack of validation testing of DEGs using an independent test set. Such validation analysis was not possible, as there is no publicly available gene database for parathyroid adenoma and normal parathyroid, and because of budget limitations. We attempted to mitigate this limitation by applying strict cut-off criteria for identification of DEGs.

Conclusions
We identified DEGs in parathyroid adenoma compared with normal parathyroid tissue and GO terms, which may contribute to the development or progression of parathyroid adenoma. Further, PPI network analysis highlighted a densely connected PPI sub-module, consisting of eight hub molecules. Additional validation of these findings, in a larger number of samples and using functional experiments, is required to confirm the results of this study.

Supplementary Materials:
The following are available online at http://www.mdpi.com/2077-0383/8/3/ 297/s1; Table S1: mRNA expression values for genes evaluated in the ten parathyroid adenomas and five normal parathyroid tissue sample, Table S2: List of differentially expressed genes (45 up-regulated and 202 down-regulated), and Table S3: GO terms extracted using down-regulated DEGs, Figure S1: PPI network of genes up-regulated in parathyroid adenoma consisting of 13 nodes and 8 edges (average node degree = 1.231; average local clustering coefficient = 0.00), showing potential functional interactions between up-regulated genes, Figure S2: PPI network constructed using down-regulated genes, which comprised 120 nodes and 230 edges (average node degree = 3.833; average local clustering coefficient = 0.229). PTH is connected with KL and GJA1 proteins (depicted in red squares), showing PTH protein functionally interacts.

Conflicts of Interest:
The authors declare no conflict of interest.