A Role for Circular Non-Coding RNAs in the Pathogenesis of Sporadic Parathyroid Adenomas and the Impact of Gender-Specific Epigenetic Regulation

Epigenetic changes, including altered small non-coding RNAs, appear to be implicated in the pathogenesis of sporadic parathyroid adenomas (PAs). In this study, we investigated the circular RNAs (circRNAs) expression profile in sporadic PAs. Sixteen tissue samples of sporadic PAs, and four samples of normal parathyroid tissue (NPT) were investigated. Sample preparation and microarray hybridization were performed based on the Arraystar’s standard protocols, and circRNAs sequences were predicted by bioinformatics tools. We identified 35 circRNAs that were differentially expressed in sporadic PAs compared to NPT; 22 were upregulated, and 13 were downregulated, according to the pre-defined thresholds of fold-change > 2.0 and p < 0.05. In the subgroup analysis of PAs from male patients (n = 7) compared to PAs from female patients (n = 9), we also find a different expression profile. In particular, 19 circRNAs were significantly upregulated, and four circRNAs were significantly downregulated in male patients, compared to female counterparts. We show here for the first time a differential circRNA expression pattern in sporadic PAs compared to NPT, and a different expression profile in PA samples from male compared to female patients, suggesting an epigenetic role in the PA pathogenesis, and also an effect of gender in the epigenetic regulation of PAs.


Introduction
Primary hyperparathyroidism (PHP) is considered to be one of the most common endocrinopathies with a prevalence of 0.3% in the general population [1] showing a predilection for postmenopausal women [2]. PHP is characterized by hypercalcemia, hypercalciouria, hypophosphatemia, and high or Tissue samples were not suitable for one of the following reasons: no tumor cells present in tumor samples (n = 7). Tissue samples were withdrawn for one of the following reasons: patients with atypical features (n = 3), analysis unsuccessful (n = 4).
We performed an immunohistochemical evaluation of the adenomatous polyposis coli (APC) gene, Ki67, cyclin D1, and parafibromin expression, in order to definitely establish a diagnosis of PA, and to exclude cases of atypical PA and/or parathyroid carcinoma [9][10][11]. The study was approved by the Scientific Committee of the affiliated AHEPA University Hospital.

Histopathological and Immunohistochemical Analysis
Tissue samples from all patients were collected, fixed in formaldehyde 10% and embedded in paraffin blocks (Formalin Fixed, Paraffin Embeded, FFPE). Sections (3 μm) from all the blocks were

Histopathological and Immunohistochemical Analysis
Tissue samples from all patients were collected, fixed in formaldehyde 10% and embedded in paraffin blocks (Formalin Fixed, Paraffin Embeded, FFPE). Sections (3 µm) from all the blocks were taken and stained with Hematoxylin-Eosin (H&E) for morphological evaluation and verification of adenoma or normal PTH, accordingly. In addition, the H&E slides were used for marking the epithelial rich areas for microdissection and tumor microarray (TMA) construction. TMAs constructed after the necessary slides for the below described microdissection were obtained. From each patient's block, three cores were obtained, each 1 mm in diameter. The cores were taken from areas with high ratio of epithelial to stromal component. From the TMA blocks, one slide was taken and stained with H&E, in order to confirm that there was adequate tissue. Another four sections (3 µm each) from each TMA block was taken and stained for Ki67 (MIB1 clone, DAKO, Glostrup, Denmark), CyclinD1 (EP12 clone, DAKO), APC (C-20 clone, Santa Cruz Biotechnology Inc., Dallax, TX, USA), and parafibromin (2H1 clone, Santa Cruz Biotechnology Inc.). All immunohistochemical stains were performed by using Bond Max and Bond III autostainers (Leica Microsystems, Wezlar, Germany). Parafibromin was evaluated as "positive" for any percentage of nuclear positivity, and "negative" for no nuclear positivity. For APC, cytoplasmic positivity was evaluated as a percentage of positive cells. For CyclinD1 and Ki67, nuclear positivity was evaluated as a percentage of positive cells. Each stain was evaluated in three cores, and the average value was taken into account for statistical analysis.

RNA Extraction from Paraffin-Embedded Tissue
Manual microdissection was performed from each FFPE block to enrich for epithelial tissue before total RNA extraction. In brief, one H&E slide and up to five unstained slides (10 µm thickness) were generated from each FFPE block. The target lesion was marked on the H&E slide for each case, which was then used to guide the removal of non-target tissues (e.g., normal parathyroid tissue) from unstained slides. After an overnight deparaffinization, the tissue area of interest was scraped off the slide into a microcentrifuge tube to facilitate total RNA extraction. RNA extraction was performed by miRNeasy FFPE kit (Qiagen GmbH., Hilden, Germany), according to manufacturer's instructions. The kit is optimized to isolate RNA molecules that are longer than 18 nucleotides from FFPE tissue samples, and to reverse as much formaldehyde modification as possible without further RNA degradation.
The concentrations of the RNA samples were determined by OD 260 by using a NanoDrop ND-1000 instrument. The integrity of RNA was assessed by electrophoresis on a denaturing agarose gel.

Sample Labeling and Microarray Hybridization
The sample preparation and microarray hybridization were performed based on the Arraystar's standard protocols (Arraystar Inc., Rockville, MD, USA). All samples were used for microarray analysis.
Briefly, total RNA from each sample was digested with Rnase R (RNase R; Epicenter, Madison, WI, USA) to remove linear RNAs, and to enrich circular RNAs. Then, the enriched circular RNAs were amplified and transcribed into fluorescent complement RNA utilizing a random priming method (Arraystar Super RNA Labeling Kit; Arraystar, Rockville, MD, USA). The fluorescently-labeled cRNAs were quantified with a NanoDrop ND-1000 spectrophotometer (NanoDrop) and hybridized onto the Arraystar Human circRNA Array V2 (8x15K, Arraystar) ( Figure S1). After having washed the slides, the arrays were scanned by the Agilent Scanner G2505C. Agilent Feature Extraction software (version 11.0.1.1) was used to analyze acquired array images. CircRNA sequences were predicted by bioinformatic methods [12].

Quantitative Real Time Reverse Transcription Polymerase Chain Reaction (qRT-PCR)
Following RNA extraction from parathyroid adenomas and normal parathyroid tissue cDNA was synthesized with a reverse transcriptase according to the kit's instructions (miScript II RT Kit, Qiagen). In brief, we used 1 ug RNA, polyadenylated by poly(A) polymerase and converted it into cDNA by reverse transcriptase with oligo-dT priming using the miScript HiFlex Buffer. The 18S and 28S ribosomal RNAs were used as internal standards. Cycling was performed under standardized conditions with 2× QuantiTect®SYBR Green PCR Master Mix on the QIAGEN Rotor-Gene Q (Corbett Rotor-Gene 6000, Hilden, Germany) real-time PCR cycler. Relative circRNA expression was calculated with the 2 −∆∆Ct method.

Prediction of circRNA-miRNA Interactions
Circular RNAs have been shown to affect miRNA-mediated regulation of gene expression through miRNA sponging [13]. Therefore, the target prediction software, based on Circbase 2.0, TargetScan and Circular RNA Interactome (NIH), was applied to predict putative circRNA/miRNA interactions for the differentially expressed circRNAs in PAs compared to NPT, identified from the microarray hybridization. The Circ-Interactome is a computational tool that enables the prediction and mapping of binding sites for RNA binding proteins (RBP) and miRNAs on reported circRNAs. This tool searches mammalian circRNAs' names and the genomic position and the best-matching transcripts of the circRNA. Then, it retrieves genomic and mature circRNA sequences and searches RBPs binding to a circRNA, and to sequences upstream/downstream of the circRNA. Finally, it identifies RBPs binding to the circRNA junctions and miRNAs targeting a circRNA [14].
Prediction of circRNA-miRNA-Target Gene Associations Several online bioinformatics platforms such as TargetScan, CircBase, MirGator, MirBase, CircNet, and NIH Circ interactome were applied in conjunction in order to predict circRNA-miRNA-target gene associations.

Statistical Analysis
Quantile normalization and subsequent data processing were performed, using the R software limma package. Differentially expressed circRNAs with statistical significance between two groups were identified through Volcano Plot filtering. Differentially expressed circRNAs between two samples were identified through Fold Change filtering. Hierarchical Clustering was performed to show the distinguishable circRNAs expression pattern among samples. When comparing the two groups (PAs versus control), for profile differences, we computed the "fold change" between the groups for each circRNA. circRNAs having more than 2-fold changes and p-values < 0.05 were selected as being significantly differentially expressed.

Results
Mean age of the included patients was 60 ± 8 years, including nine postmenopausal females and seven males. The control group included samples of NPT from two males and two females mean age 65.7 ± 3.7 years. The demographic and clinical characteristics of the included patients are depicted in Table 1. As expected, patients with PHP had significantly higher levels of serum calcium and PTH levels compared to controls (sCa: 2.7 ± 0.2 mmol/L vs. 2.28 ± 0.09 mmol/L, p < 0.001, and sPTH: 12.68 ± 3.5 pmol/L vs. 4.42 ± 0.57 pmol/L, p < 0.001, respectively). None of the patients had significant level of hypophosphatemia, although serum phosphate was in the low normal range. After parathyroidectomy, there was a complete biochemical remission of primary hyperparathyroidism in all patients with serum PTH, calcium, and PO 4 reaching normal values within one to two months post-surgery.

Immunohistochemical Findings
Representative images of H&E-stained parathyroid tumor specimens are provided in Figure 2. All tested samples, both PAs and NPT were positive for parafibromin ( Figure 2, Table S1). Loss of APC expression was found in one case of PA, while all other cases were positive. Ki67 expression was low in all cases of PAs, reaching up to 1%, and negative in all NPT samples. Expression of cyclinD1 was significantly higher in PAs compared to NPT, as previously described [15,16]. Immunohistochemical findings confirmed the diagnosis of PAs.

Identification of Differential circRNA Profiles
The Arraystar Human circRNA Microarray was used to sample labeling and microarray hybridization of all the finally selected parathyroid-tissue samples ( Figure 3).

Identification of Differential circRNA Profiles
The Arraystar Human circRNA Microarray was used to sample labeling and microarray hybridization of all the finally selected parathyroid-tissue samples ( Figure 3).
We identified 35 circRNAs that were differentially expressed in sporadic PAs compared to normal parathyroid tissue; 22 were upregulated, and 13 were downregulated (Table 2), according to the pre-defined thresholds of fold-change > 2.0 and p < 0.05. Among upregulated circRNAs; hsa_cirRNA_051778, hsa_cirRNA_402533, hsa_cirRNA_406174, and hsa_cirRNA_0008267 showed the highest fold-change (>5-fold), whereas hsa_circRNA_032603 and hsa_circRNA_058097 showed the highest fold-change (>3-fold) among the downregulated circRNAs. Figure 3. Scatter plots and volcano plots to identify differentially-expressed circular RNAs (circRNAs). Scatter plots used to identify differentially-expressed circRNAs in (A) PAs vs NPT and (B) PAs from F patients vs PAs from M patients. The y-axis represents the mean normalized circRNA signal values for each comparator group (log2-scaled). The green fold-change lines represent 2.0× fold-changes, so the circRNAs lying above and below these green lines displayed greater than a 2.0fold upregulation or downregulation. Volcano plots used to identify differentially-expressed Scatter plots used to identify differentially-expressed circRNAs in (A) PAs vs. NPT and (B) PAs from F patients vs. PAs from M patients. The y-axis represents the mean normalized circRNA signal values for each comparator group (log2-scaled). The green fold-change lines represent 2.0× fold-changes, so the circRNAs lying above and below these green lines displayed greater than a 2.0-fold upregulation or downregulation. Volcano plots used to identify differentially-expressed circRNAs in (C) PAs vs. NPT and (D) PAs from F patients vs. PAs from M patients. The x-axis represents fold-change values (log2-scaled), while the y-axis represents p-values (−log10-scaled). The green vertical lines correspond to 2.0× upregulation and downregulation, respectively, while the green horizontal line corresponds to a p-value of 0.05. On this basis, the red rectangles represent the differentially-expressed circRNAs of statistical significance. NSV, normalized signal value; PAs, parathyroid adenomas; NPT, normal parathyroid tissue; F, female; M, male.

qRT-PCR Validation of Differential circRNAs
The four significantly upregulated circRNAs and the two significantly downregulated circRNAs with the highest fold change were selected for qRT-PCR validation in the same samples that were used for the microarray analysis. The qRT-PCR results revealed that all six circRNAs but one (namely hsa_circ_RNA_402533) were significantly differentiated between PAs vs. normal parathyroid tissue, with fold-changes mirroring those derived from the microarray data ( Figure 4).

qRT-PCR Validation of Differential circRNAs
The four significantly upregulated circRNAs and the two significantly downregulated circRNAs with the highest fold change were selected for qRT-PCR validation in the same samples that were used for the microarray analysis. The qRT-PCR results revealed that all six circRNAs but one (namely hsa_circ_RNA_402533) were significantly differentiated between PAs vs normal parathyroid tissue, with fold-changes mirroring those derived from the microarray data ( Figure 4).
We also searched for the total miRNA binding site (miRBS) that has been reported so far for the three most upregulated and the two most down-regulated circRNAs. We identified 46 miRBS for hsa_circRNA_008627; 41 for hsa_circRNA_058097; 18 for hsa_circRNA_406174; seven miRBS for hsa_circRNA_032603; one for hsa_circRNA_402533.

circRNA-Genes Interactions
To further explore the underlying molecular mechanism(s) of the differentially expressed circRNAs in sporadic PAs, we used the bioinformatics platforms CircBase, and NIH Circ-Interactome to predict circRNA-genes interactions. The predicted genes associated with the differentially expressed circ-RNAs are depicted in Table 2. Two circRNAs, hsa_circRNA_404337 and
We also searched for the total miRNA binding site (miRBS) that has been reported so far for the three most upregulated and the two most down-regulated circRNAs. We identified 46 miRBS for hsa_circRNA_008627; 41 for hsa_circRNA_058097; 18 for hsa_circRNA_406174; seven miRBS for hsa_circRNA_032603; one for hsa_circRNA_402533. The experimentally assayed circRNAs were cross-checked for the genomic position and the best transcript in circBase (http://circbase.org). ‡ ‡: The circRNA -corresponding hsa-miRNA targets with the best context+score percentile (NIH Circular RNA interactome). All identified circularRNAs are exonic arising from the exons of the linear transcript except for: **: Antisense, where their gene locus overlaps with the linear RNA, but is transcribed from the opposite strand. ***: Sense overlapping, transcribed from same gene locus as the linear transcript, but not classified into "exonic" and "intronic". has: Homo sapiens; Circ: circular; chr: chromosome.

circRNA-Genes Interactions
To further explore the underlying molecular mechanism(s) of the differentially expressed circRNAs in sporadic PAs, we used the bioinformatics platforms CircBase, and NIH Circ-Interactome to predict circRNA-genes interactions. The predicted genes associated with the differentially expressed circ-RNAs are depicted in Table 2. Two circRNAs, hsa_circRNA_404337 and hsa_circRNA_051799, interact with genes that have been shown to associate with molecular pathways that are involved in the pathogenesis of sporadic PAs, such as CDKN1B and CDK1, respectively [17,18].
Most of the genes identified are protein-coding, and associated with several cellular processes, such as organization and re-arrangement of the cytoskeleton, cell proliferation, cell metabolism, and transcription (Table S2).

Parathyroid Adenomas from Male vs. Female Patients
In order to test whether the gender of the patient has an impact in the expression profile of circRNAs in PAs, since the vast majority of PHP due to PAs appear in postmenopausal women, we performed a sub-group analysis of PAs from male patients compared to PAs from female patients. We found a different expression profile between male and female patients. In particular 19 circRNAs were significantly upregulated and four circRNAs were significantly downregulated in male patients compared to female counterparts (Table 3). It is worth noting that the identified differentially expressed circRNAs in the subgroup analysis between male and female patients were different from the ones that we have identified in our initial analysis between PA samples and controls, and there was no overlap in circRNAs expression profiles between the two analyses.  All identified circularRNAs are exonic, arising from the exons of the linear transcript, except for *: intronic, arising from an intron of the linear transcript. **: Antisense, where their gene locus overlap with the linear RNA, but transcribed from the opposite strand. ***: Sense overlapping, transcribed from same gene locus as the linear transcript, but not classified into "exonic" and "intronic". hsa: Homo sapiens, Circ: circular, chr: chromosome.
The predicted miRNAs for the significantly upregulated and downregulated circRNAs are also listed in Table 3, with hsa_miR_1184 being the most frequently targeted miRNA by 4 up-regulated circRNAs (hsa_circRNA_008636, hsa_circRNA_000763, hsa_circRNA_092388, hsa_cirRNA_044837) and one down-regulated circRNA (hsa_circRNA_100332). The predicted genes associated with the differentially expressed circRNAs in this analysis are also depicted in Table 3. The majority of these genes are protein coding (Supplemental Table S3) and involves the intracellular ubiquitin proteasome system, nucleic acid and RNA binding, protein heterodimerization and homodimerization activity, vesicle trafficking and alternative splicing of pre-mRNAs. None of the implicated genes, however, is located in the X chromosome (Table 3), and thus we cannot address whether incomplete inactivation of the X chromosome is a possible mechanism, contributing to the altered expression of some circRNAs in males vs. females.
In order to further evaluate our results in the differential expression of sporadic PAs in male vs. female patients, we also conducted a subgroup analysis between samples from males and females in the control group, to test whether the altered expression of circRNAs profile is disease-independent. We did not find any significant difference in this analysis, suggesting that the differences in circRNA expression between male and female patients may be attributed to PA, but the numbers were small.

Discussion
In this study, we investigated the circRNA expression profile in sporadic PAs, and we tested the effect of gender.
Thirty-five circRNAs, 22 up-regulated, and 13 down-regulated were identified by microarray analysis to be differentially expressed in sporadic PAs compared to NPT. Applying a set of bioinformatic tools, we constructed a network of circRNA-miRNA and genes interactions, identifying two of the differentially expressed circRNAs, hsa_circRNA_404337 and hsa_circRNA_051799, that interact with genes known to associate with molecular pathways involved in the pathogenesis of sporadic PAs, such as CDKN1B and CDK1 [4,19], respectively. In the sub-analysis within the cohort of sporadic PAs, we also report a differential circRNA expression profile in sporadic PAs from male patients compared to PAs from female counterparts with 19 circRNA being significantly upregulated and four being significantly downregulated. Primary hyperparathyroidism is most frequently present in postmenopausal women, and although the involvement of estrogen has been suggested, and parathyroid tumors are found to express ERβ1 and ERβ2 [20], the exact underlying pathophysiology remains obscure. The differences in circRNA expression profile between male and female patients, but not between male and female controls, point to a potential gender effect on the PA tumorigenesis through epigenetic regulation.
In the past decades an ever-growing list of diverse non-coding RNAs with functional capacity have been identified in eukaryotic cells [21]. CircRNAs are formed by alternative splicing of pre-mRNA, in a process known as 'backsplicing', where an upstream splice acceptor is joined to a downstream splice donor [22]. They are expressed in a variety of eukaryotic organisms, demonstrate conservation across mammals, they are abundant, expressed in a regulated manner independent of their cognate linear isoform and in a cell type-specific manner.
In cancer, circRNAs appear to have a role in cancer development and progression, by interacting with microRNAs. CircRNAs negatively regulate the inhibitory effects of miRNAs upon their target linear mRNAs by directly binding to miRNAs via MREs, thereby essentially functioning as 'miRNA sponges' within the cellular RNA network. Another potential function of circRNAs could be to compete with the splicing of an mRNA [38], while some circRNAs have also been implicated in transcriptional or post-transcriptional gene regulation of their host genes [39] (Figure 5). Compared to mRNAs, circRNAs are more stable, due to their natural resistance to exonucleases, with half-lives greater than 24-48 h, which is significantly longer than the half-life of an average mRNA (8-9 h in human cells) [40,41]; this is why they are considered as being much more promising in understanding the pathophysiology of certain diseases and the development of cancer.
The emerging role of the epigenetic component in tumorigenesis has just now started to be identified in parathyroid neoplasia [6,19].
The mature form of miRNA_24-1 was studied in the parathyroid hyperplasia of MEN1. It has been shown that miRNA_24-1 is uniquely expressed in MEN1 parathyroid adenoma tissues that conserve the wild-type allele, while it is downregulated in the MEN1 parathyroid adenoma tissues with the 11q13 LOH, as well as in sporadic parathyroid adenomas [44].
We also search for miRNA response elements (MREs) on the 35 differentially expressed circRNAs, which were then used to identify putative miRNAs that were based on their seed-sequence complementarity [13]. Hsa_miR_1248 was the most frequently encountered miRNA, and although little is known about the potential targets and pathways regulated by this miRNA it has been linked to the pathogenesis of breast [45] and non-small cell lung carcinoma [46], autism spectrum disorders [47], asthma [48], and in the aging process [49]. Genome-wide microarray and in vitro studies [48,49], have shown that miRNA-1248 regulates the expression of genes that play a role in interleukin-related inflammatory pathway, as well as, DNA repair pathways.
None of the miRNAs that were previously reported in parathyroid carcinomas, however, was identified in our analysis through interaction with the reported dysregulated circRNAs.
Apart from the dysregulated expression of non-coding RNAs, the hypermethylation of APC, the gene responsible for familial adenomatous polyposis, occurs frequently in parathyroid adenomas although aberrant APC expression has not been demonstrated [50,51]. In addition, genes that are frequently subject to epigenetic inactivation in human cancers, such as Ras association domain-containing protein 1 (RASSF1A) [52] and hypermethylated in cancer 1 protein (HIC1) [53], are also frequently hypermethylated in parathyroid adenomas, but it is unclear which of these aberrantly methylated genes may be important to the pathogenesis of parathyroid tumors.
Our study has several limitations; (i) the sample size was relatively small, (ii) our analysis was limited to parathyroid tissue, and thus results from other biosamples such as serum or plasma before and after parathyroidectomy that could strengthen our results are missing and (iii) we were not able to perform DNA analysis and investigation of somatic mutations in our samples, mainly due to inadequate material in the total of PA samples used.
Despite these limitations, however, this is the first study investigating the circRNA expression profile in sporadic PAs. We identified 35 differentially expressed circRNAs in PA samples compared to normal parathyroid tissue, suggesting a role of circRNAs in the pathogenesis of these benign tumors. In addition, we reported significant differences in the circRNA expression profile between male and female patients underlying the fact that epigenetics may be involved in the gender specific differences reported in PAs. Furthermore, our results from the bioinformatic analysis reveal a circRNA-miRNA network that could play a key role in the triggering of the tumorigenetic process in the parathyroid cells. As circRNAs regulate the epigenome in multiple ways ( Figure 5), future functional studies that can pave the way between 'circnome' and PA pathogenesis would significantly increase our knowledge for the molecular pathways in involved.
CircRNA network appears to be cell and tissue specific [37], stable, and with a significant regulatory role in pathogenesis of diseases, providing new strategic approaches to tumorigenesis, and emerge as promising diagnostic and prognostic biomarkers, and as therapeutic targets. functional studies that can pave the way between 'circnome' and PA pathogenesis would significantly increase our knowledge for the molecular pathways in involved. CircRNA network appears to be cell and tissue specific [37], stable, and with a significant regulatory role in pathogenesis of diseases, providing new strategic approaches to tumorigenesis, and emerge as promising diagnostic and prognostic biomarkers, and as therapeutic targets.