A Pilot Genome-Scale Profiling of DNA Methylation in Sporadic Pituitary Macroadenomas: Association with Tumor Invasion and Histopathological Subtype

Pituitary adenomas (PAs) are neoplasms that may cause a variety of neurological and endocrine effects. Although known causal contributors include heredity, hormonal influence and somatic mutations, the pathophysiologic mechanisms driving tumorigenesis and invasion of sporadic PAs remain unknown. We hypothesized that alterations in DNA methylation are associated with PA invasion and histopathology subtype, and that genome-scale methylation analysis may complement current classification methods for sporadic PAs. Twenty-four surgically-resected sporadic PAs with varying histopathological subtypes were assigned dichotomized Knosp invasion scores and examined using genome-wide DNA methylation profiling and RNA sequencing. PA samples clustered into subgroups according to functional status. Compared with hormonally-active PAs, nonfunctional PAs exhibited global DNA hypermethylation (mean beta-value 0.47 versus 0.42, P = 0.005); the most significant site of differential DNA methylation was within the promoter region of the potassium voltage-gated channel KCNAB2 (FDR = 5.11×10−10). Pathway analysis of promoter-associated CpGs showed that nonfunctional PAs are potentially associated with the ion-channel activity signal pathway. DNA hypermethylation tended to be negatively correlated with gene expression. DNA methylation analysis may be used to identify candidate genes involved in PA function and may potentially complement current standard immunostaining classification in sporadic PAs. DNA hypermethylation of KCNAB2 and downstream ion-channel activity signal pathways may contribute to the endocrine-inactive status of nonfunctional PAs.


Introduction
Pituitary adenomas (PAs) are typically benign monoclonal neoplasms with an overall prevalence of 16.7% (14.4% in autopsy studies and 22.5% in imaging studies) in the general population. The majority of PAs, however, are small and nonfunctional tumors, and only 0.16-0.20% of them are macroadenomas $ 10 mm in diameter [1,2]. Although histologically benign, many PAs may cause significant morbidity due to their anatomical location, often resulting in tumor mass effect and neurological symptoms in addition to causing hormonal over-secretion or hypopituitarism. Invasion into surrounding anatomical structures (e.g. cavernous sinus invasion) remains a major barrier to achieving long-term tumor and disease control, especially in cases of functional PAs resulting in malignant endocrinopathies such as Cushing's disease or acromegaly [3]. With regard to the phenotype of invasion, PAs are often classified according to the Knosp grading system, in which one of five grades is assigned based on the degree of cavernous sinus invasion and relationship to the internal carotid artery [4].
Previous studies have shown that genetic mutations play an important role in the tumorigenesis of familial PAs, particularly those with heritable mutations in the multiple endocrine neoplasia I (MEN1) and aryl hydrocarbon receptor interacting protein (AIP) genes [5]. In sporadic PAs, however, it has been suggested that alterations in epigenetic regulation, and particularly DNA methylation, may play a particularly prominent role in PA tumorigenesis and invasion, likely via loss or reduced expression of tumor suppressor genes (TSGs) [6]. Many known TSGs have been shown to harbor C-phosphate-G (CpG) island hypermethylation in PAs, which is associated with and frequently results in silencing of TSGs. For example, DNA methylation within promoter and exon 1 regions of the cyclin-dependent kinase inhibitor (p16/ CDKN2A) gene was found at a high frequency in NFAs [7,8], and DNA methylation of the CpG island like cell-cycle regulatory genes of growth arrest and DNA-damage-inducible, gamma (GADD45G) [9] and apoptosis gene of rhomboid domain containing 3 (RHBDD3) have been associated with PA evolution [10]. Other studies have shown that tumor specific epigenetic silencing of cadherin 13, H-cadherin (CDH13) and cadherin 1, type1, E-cadherin (CDH1), alone or in combination, were involved in PA development and invasion [11], and CpG hypermethylation-mediated glutathione S-transferase pi gene (GSTP1) inactivation was a common finding in PAs potentially contributing to their invasive behavior [12]. Few prior studies, however, have utilized genome-scale approaches according to functional PA subtypes and phenotypical invasion status to identify candidate genes involved in these processes [13,14].
Although several targeted genes have been associated with PA tumorigenesis and progression, primary whole genome-scale epigenetic alterations remain largely unknown, especially as they pertain to invasion and histopathological PA subtype classification. In the current study, we utilized genome-scale DNA methylation technology to identify DNA methylation alterations between invasive and noninvasive PAs subtypes, and across varying functional classes of PAs.

Patients and tissue specimens
This study was conducted according to the Helsinki human subject doctrine and was approved by the Institutional Review Boards of the Keck School of Medicine of the University of Southern California (Los Angeles, USA). Written informed consent was signed and obtained from 23 participants and one guardian on behalf of a minor enrolled for tissue specimen collection and subsequent analysis. Twenty-four patients with surgically-resected PAs from the University of Southern California (USC) Keck Hospital and Los Angeles County + USC Medical Center were included in this retrospective study.
Of the 24 PAs, six (25%) were functional and 18 (75%) were nonfunctional. More specifically, these included 17 (70.8%) nonfunctional adenomas, 5 (20.8%) somatotroph adenomas, 1 (4.2%) corticotroph adenoma, and 1 (4.2%) silent corticotroph adenoma. All PAs analyzed were macroadenomas $ 10 mm in diameter, and were diagnosed based on laboratory evaluation and Magnetic Resonance Imaging (MRI), with an average tumor diameter of 23.55 TR619.36 CC mm (Table 1). Knosp invasion scoring on MRI was performed by a staff neuro-radiologist who was blinded to the genomics analysis. For the purposes of this study, PAs with Knosp scores of 0-1 were classified as noninvasive, and those with Knosp scores of 2-4 were classified as invasive. Equal numbers of invasive and noninvasive PAs (12 each) were selected.
Following transsphenoidal tumor resection, PA specimens were fresh-frozen in liquid nitrogen and stored at 280uC until used for DNA and RNA extraction. Genomic DNA purification was processed using DNeasy Blood & Tissue kit (Qiagen) according to the standard protocols of the manufacturer. Total RNA was purified using RNeasy Plus Universal Kits (Qiagen). Tissue samples were first cut into small pieces with no more than 50 mg on dry ice then moved promptly into 1.5 ml precooled tubes. Liquid nitrogen was then added into the tube, and tissue was grinded with a grinding rod before adding 0.9 ml of QIAzol Lysis Reagent. The tubes were vortexed for 60 seconds and left at room temperature for 5 minutes before following the standard RNeasy Plus protocol for RNA purification.

Genome-scale DNA methylation array
One microgram of bisulphite converted DNA (using the Zymo EZ DNA Methylation kit, Zymo Research, Irvine, CA) per sample was used for genome-scale DNA methylation profiling, and was performed in the USC Epigenome Center using the Illumina Infinium HumanMethylation450 (HM450) Beadchip platform, which interrogates 482,421 DNA methylation sites and covers 99% of Reference Sequence of National Center for Biotechnology Information (NCBI-RefSeq) genes, with an average of 17 CpG sites per gene region distributed across the promoter, TSS1500 (1500 bp within transcription start site), TSS200 (200 bp within transcription start site), 59UTR, 1stExon (the first exon), gene body, 39UTR. It covers 96% of CpG islands, with additional coverage in island shores and the regions flanking them. The DNA methylation data was submitted to Gene Expression Omnibus (GEO) with access number of GSE54415. Beta values for all HM450 probes were calculated after background correction and normalization as previously described [15]. Probes (11,648) located on chromosomes X and Y were excluded. Probes targeting non-CpG sites (3,091) and probes targeting 65 known single nucleotide polymorphisms (SNPs) were also excluded. An additional set of 86,560 probes with common SNPs at the CpG site with a minor allele frequency (MAF) greater than 1% (identified using NCBI dbSNP builds 128), or with common SNPs within 10 bp from CpG site, or within 15 bp from the CpG site lying entirely within a repeat region (from RepeatMaster and Tandem Repeat Finder databases) were excluded. Finally, 591 probes with detection P-values .0.01 were removed. Following these filtering steps, a total of 383,718 probes were utilized for analyses ( Figure A in File S1). We determined the distribution of each probe by CpG density (island, shore, shelf, etc) and by location (promoter, gene body, UTR, etc.) ( Figure A and B in File S2). A secondary analysis of probes with beta values ,0.1 across 77.8% (14/18 samples) of nonfunctional adenomas were found to be less variable and therefore could conceal the potential variance. These probes were excluded in order to investigate the underlying difference between invasive and non-invasive nonfunctional adenomas. Remaining subsets of 884 probes were selected for the in-depth analysis ( Figure B in File S1). Finally, probes with Dbeta .6 0.3 were used in differentiating various PA classes based on phenotype.
Hierarchical Clustering (HCL) was performed using the Microarray Software Suite of MultiExperiment Viewer (MeV) V4.8 [16] using probes with the top 1%, 2% and 5% of standard deviation (SD) across the set of 383,718 ''global'' probes in 24 PA samples respectively. Two-way clustering was performed using Pearson's correlation and average linkage for both the gene tree and sample tree ( Figure 1 with top 2% SD). Differential DNA methylation analyses were performed between: 1) Invasive versus noninvasive PAs (12 samples each), 2) FA subjects (n = 6) versus NFA subjects (n = 18), 3) Invasive NFAs (n = 10) versus noninva-  sive NFAs (n = 8), and 4) PAs and NFAs stratified according to the five different Knosp invasion grades (0-4). Two tailed t-test analyses were performed between FAs versus NFAs and noninvasion versus invasion PA. Kruskal-Wallis test was subsequently used to analyze DNA methylation differences among 18 NFA subjects with five different Knosp grades ( Table 1). All of the analysis were executed in globe probes with the proportion of false significant genes did not exceed 0.01.

MethyLight DNA methylation analysis
Validation of findings from DNA methylation profiling in KCNAB2, which spans the CpG site interrogated by the Illumina HM450 array at probe cg18192083, was performed using MethyLight real-time PCR assays as previously described [17]. One microgram of genomic DNA from each of the 24 PA samples and the CpG methyltransferase (M.SssI) treated reference DNA sample was converted with bisulfite using the Zymo EZ DNA Methylation kit (Zymo Research, Irvine, CA), as specified by the manufacturer. A control reaction targeting ALU repetitive elements [18] was used to normalize the input bisulfite-DNA amounts. The real-time PCR reactions were obtained from Life Technologies, Carlsbad, CA. The KCNAB2 MethyLight primer and probe sequences are: forward: 59-GGT TTT TTT ATT TGG GTT ACG CG-39; reverse: 59-CGC AAA ACT AAA AAA CCT AAC GC-39; probe: 59-6FAM-CCG TAA AAT ATC GAA ACG TAA CC-MGBNFQ-39, in which 6FAM is the fluorophore and MGBNFQ refers to a minor-groove binding non-fluorescent quencher. The ALU control reaction primer and probe sequences were previously described [18]. The data were reported as a Percent of Methylated Reference (PMR), calculated as: PMR = (((KCNAB2/ALU) sample )/(KCNAB2/ALU) M.SssI )) 6 100.
Using HM450 manifest (downloaded from gene expression omnibus, GPL13534) 86,772 promoter-associated CpG probes were selected for the pathway analysis using gene set enrichment analysis (GSEA) software [19]. A total of 1,454 gene sets of Human Gene Ontology were exported from the molecular signatures database (MSigDB) version 3.1, and after filtering, 469 of the gene sets were either larger than 500 or smaller than 25 gene set size allowed. Additional statistical analysis was performed using R software (http://www.r-project.org/) and GraphPad Prism Software (GraphPad Software Inc.).

Whole-genome gene expression sequencing
One microgram of total RNA from each of the 24 samples was assessed by Agilent Bioanalyzer to examine quality, and subsequently used for library preparation. Library quantization was performed by Kapa Biosystems qPCR assay. RNA sequencing was performed in the USC Epigenome Center using Illumina Hi-seq 2000 platform to produce 100 bp paired-end reads, by indexing four samples per lane. The Fastq files could be downloaded from Sequence Read Archive (SRA) with accession number SRP035646. Raw sequence data were first examined by FastQC as a quality control measure. We next used Tophat and Cuffdiff to examine genes with differential expression, based on the fold change of fragments of per kilo base of transcript sequence per millions base pairs (FPKM).

Genome-wide DNA methylation in invasive and noninvasive PAs
In this pilot study, twenty-four surgically-resected sporadic PAs with varying histopathological subtypes and Knosp invasion scores were characterized to determine if genome-wide and gene-specific methylation levels can serve as biomarkers of PA subtype and invasion grading. We calculated the mean DNA methylation beta value for the entire filtered probe set (383,718 probes) as a measure for global DNA methylation in noninvasive and invasive PAs (Figure 2A), across different tumor grades by Knosp categorization ( Figure 2B), and in FAs and NFAs ( Figure 2C). Based on these calculations, no significant differences in global DNA methylation levels were observed between the invasive and noninvasive groups ( Figure 2A). Furthermore, Kruskal-Wallis test analysis among PAs with five different Knosp grades using the global filtered probe set also did not show significant DNA methylation differences as a function of tumor grade ( Figure 2B), even though we relaxed the threshold to FDR ,0.05. Within the set of 18 NFA subjects, similar (nonsignificant) findings between invasive and noninvasive PAs were also seen ( Figure 2D). These data suggest that global DNA methylation profiles may not independently predict clinically significant differences in the invasive phenotype of PAs.
We further refined the probe set for unsupervised analyses to remove non-variably methylated probes. After filtering, 884 probes were selected, of which 355 were enhancer-associated, 20 were promoter-associated, and 4 were both promoter-and enhancer-associated based on the Illumina HM450 manifest for probes ( Figure B in File S1). We identified 34 CpGs (localized to 17 genes) that were independently associated with enhancers and were hypomethylated in invasive NFAs compared to noninvasive NFAs (FDR ,0.01) (Table S1). These differentially-methylated genes included the fms-related tyrosine kinase 1 (FLT1) and slit homolog 3 (SLIT3), which are two important genes known to be associated with cell motility and invasion. Although we attempted to validate methylation status in these two genes using Methy-Light, we were not able to because the probes lie in CpG-depleted regions of the genes. No DNA methylation differences in promoter regions were noted between invasive and noninvasive NFAs.

Nonfunctional PAs are globally hypermethylated compared to somatotroph PAs
The group of 24 PA samples independently clustered into two subgroups of nonfunctional and functional (all somatotroph) PAs following HCL with top 2% of SD ( Figure 1). All somatotroph adenomas were exclusively clustered from NFAs (and one additional functional corticotroph adenoma) (Figure 1), and further clustering of the top 1% (File S3) and 5% (File S4) of SD validated this finding. The mean global beta values from FA (six subjects) and NFA (18 subjects) tumors showed that NFAs were significantly DNA hypermethylated compared to FAs ( Figure 2C). However, NFAs and FAs showed similar DNA methylation trends from gene promoters through gene bodies and UTR regions ( Figure 2E and F). These findings suggest that DNA methylation analysis may provide insight into hormonal activity in functional versus nonfunctional PAs, and possibly complement current PA histopathological classification systems.
In order to interrogate which subsets of genes were involved in distinguishing specialization of nonfunctional PAs, we filtered the global panel of 383,718 HM450 probes to include those targeted in RefSeq genes for a more thorough investigation. Using a twotailed t-test, we identified 3,027 CpG sites with significant DNA methylation differences between NFAs and FAs (FDR ,0.01). In order to display the methylation status of the significant CpG sites, we graphed DNA methylation levels of these probes stratified by location relative to gene region ( Figure 3A) and as a function of CpG density ( Figure 3B). We further restricted the set of 3,072 probes to 360 CpG sites associated with promoter regions, and found substantial DNA hypomethylation in FA cases compared to NFA cases (simple t-test, P,0.01) ( Figure 3C). Furthermore, CpG sites were stratified by density ( Figure 3D) as well. The most significant gene with regard to differential DNA methylation was KCNAB2, a member of the potassium voltage-gated channel, shaker-related subfamily (FDR = 5.11610 210 ), which was hypermethylated in the NFA group. This result was subsequently validated using MethyLight real-time PCR assays ( Figure 3E). Furthermore, we found that among the top 30 most significantly differentiated genes, 4 (NME9, PSEN2, MDN1 and PCSK6) were hypomethylated in NFAs compared to FAs (Table 2). Interestingly, we also corroborated the recent finding discovered using Illumina Infinium Methylation 27K Arrays that the transcription factor AP-2 epsilon (TFAP2E) and 13 other genes (Table S2) were hypermethylated in NFAs [13] (Table 2). Epigenetic alteration of the TFAP2E gene has also been commonly reported in selected human cancers, including prostate cancer and colorectal carcinoma [20,21]. In addition, Wnt signaling plays an important role in adenoma development and tumorigenesis, including PAs [22]. Whether hypermethylation of TFAP2E could down-regulate Wnt signaling, however, was still unknown. In order to further understand the detailed DNA methylation status of KCNAB2, we performed hierarchical clustering of beta values from KCNAB2targeted probes on the HM450 array using the top 50% highest standard deviation. In nonfunctional PAs, KCNAB2 was generally hypermethylated across the promoter region (TSS1500 and TSS200), 59UTR, the first exon, and gene body. This was substantially different from the coverage profile observed in functional PAs (File S5). Since KCNAB2 is expressed abundantly in the nervous system and T lymphocytes, and appears to play an important role in K + channel activation [23] and subsequent hormone secretion [24], we further extended our investigation of this pathway.
Methylation of ion channel activity signal pathway genes may be associated with PA functional status Because gene promoter DNA hypermethylation is known to play an important role in tumorigenesis [25] and contributes to loss of gene function and silencing [26], we selected promoterassociated probes according to the Illumina HM450 manifest for further pathway analyses. We found that four gene sets were enriched (GSEA reports the gene sets with FDR ,0.25) in NFAs, and no such gene sets were detected in FAs (Table S3, File S6). The most significant gene set was related to ion-channel activity signaling. Twenty core enrichment genes are shown in Table S4 and the pathway heat map is shown in File S7. The genes enriched in the pathway encoded K + , Cl 2 and Ca 2+ ion-channels, which are known to play important roles in hormone secretion [27,28]. In order to determine whether there were any pivotal genes driving these pathways, we performed a leading edge analysis of the four significant gene sets. Two core genes, KCNMB4 (potassium large conductance calcium-activated channel, subfamily M, beta member 4) and CACNA1C (calcium channel, voltage-  CpG islands (CGI), CpG island shores (0-2 kb from island edge, N and S indicate the upstream and downstream of the island, respectively), CpG island shelves (0-2 kb from shore edge, N and S indicate the upstream and downstream of the island, respectively) and non-CpG island probes. C and D, Relative distribution and count of the promoter associated CpGs (total 360), which were singled out from the global significant CpG sites (total 3027). E, KCNAB2 (cg18192083) MethyLight assay in NFA and FA, which was consistent with the genome-wide DNA methylation analysis. doi:10.1371/journal.pone.0096178.g003 Table 2. Top 30 significant genes with differential DNA methylation between NFAs and FAs. dependent, L type, alpha 1C subunit) overlapped in the significant gene sets, and the ion-channel activity signal gene set highly overlapped with the substrate-specific channel activity signal gene set ( Figure A and B in File S8). These data suggest that ionchannel pathways may play an important role in the process of hormonal secretion in functional PAs. As such, identified hypermethylation pathway genes, and in particular KCNAB2, may provide new insight to the epigenetic differentiation of FAs and NFAs.

Association of DNA methylation with gene expression in sporadic PAs
Gene expression can be modulated in multiple ways, from transcription to post-translational modification. Epigenetic modification may affect gene expression levels at multiple steps, including DNA methylation, histone modification, chromatin remodeling and transcription factor interaction. To examine whether DNA methylation was associated with gene expression in PAs, we performed differential gene expression analysis among the 24 subjects using RNA-Seq.
Comparison between NFAs and FAs showed that gene expression was not found to be significantly different (FDR, 0.01) between the two groups. However, when relaxing the threshold to claim statistical significance, 27 genes (FDR,0.05) displayed a potentially different expression profile (Table 3). Based on the globally significant methylated CpG sites, we selected the relative hypermethylation genes in NFAs and investigated their gene expression patterns. We found that the trend of gene expression was negatively associated with DNA methylation ( Figure 4A). When matching those 27 genes with significant genes identified by DNA methylation analysis, three genes, including ribosomal protein S6 kinase, 90 kDa, polypeptide 2 (RPS6KA2), retinol dehydrogenase 10 (RKH10), and odontogenic ameloblast associated protein (ODAM) were identified with altered methylation CpG sites in the gene body ( Figure 4B). We found that genes RDH10 and ODAM, which were hypomethylated in the CpG island and enhancer region, showed increased expression in FAs compared to NFAs. However, the RPS6KA2 gene, which was hypermethylated in a CpG site in the gene body, also showed increased expression in FAs compared to NFAs. Finally, we found that ODAM was expressed in 16.7% (1/6) of FAs but showed no expression in NFAs (0/18) ( Figure 4C). These results indicate that the function of DNA methylation seems to vary with context, and the relationship between DNA methylation and transcription is more complex than previously expected [29].

Discussion
Genome-scale DNA methylation screening of PAs was performed to investigate whether DNA methylation was associated with PA invasion and histopathological subtype classification. To our knowledge, this is the first genome-scale DNA methylation analysis of sporadic PAs, demonstrating that epigenetic modification of key gene substrates may in part account for functional differentiation of PAs, and that DNA methylation analysis of key candidate genes may potentially be used to complement PA current histopathology subtype classification systems.
Tumor invasion has been observed to occur in up to 85% of surgically-resected PAs, based on microscopic analysis of dural samples [30]. Consequently, tumor invasion is perhaps the greatest barrier to achieving adequate tumor control in PAs, as complete surgical resection of noninvasive PAs is typically achieved in 70-95% of noninvasive PA cases, compared with only 20-40% of invasive PAs [31]. In our study, no significant global differences in CpG methylation were identified between invasive and noninvasive PAs. Although CpG island hypermethylation of a series of well-characterized cell cycle regulation genes, including retinoblastoma 1 (RB1) [32], CDKN2A [33,34], and GADD45G [35] have been linked to gene repression and negative regulated cell growth in PA, these genes were not demonstrated to contribute to the invasive PA phenotype in our study. Importantly, no significant differences in CDKN2A methylation were observed between noninvasive and invasive NFAs in our study, a finding which is consistent with prior observations [7], and may also suggest that CDKN2A inactivation is alternatively related to PA subtype and size [36]. The current findings suggest that genome-scale DNA methylation assessment may be utilized to identify candidate genes involved in PA invasion that would otherwise possibly be concealed in targeted gene studies. In our study, a secondary analysis showed that the angiogenic gene SLIT3 [37] and oncogene FLT1 [38] may be two candidate biomarkers for invasive PAs, if validated in future studies with a larger PA sample size. Meanwhile, although a previous study showed that promoter hypermethylation of CDH13 and CDH1 was detected in PAs but not in normal pituitary tissue, and promoter hypermethylation of CDH13 was observed more frequently in invasive PAs [11], there was no significant global methylation difference between invasive and noninvasive PAs in our study. We did, however, find gene body hypermethylation in CDH7, a type 2 classical cadherin from the cadherin superfamily, when comparing noninvasive to invasive NFAs (FDR = 0.01). Further validation in a larger PA sample population is nevertheless warranted to identify and confirm genomic and epigenomic pathways involved in tumor invasion. Meanwhile, other epigenetic mechanisms (i.e. histone acetylation) should be explored in an effort to explain the character of PA invasion and prioritize genes associated with this clinically significant phenotype.
The anterior pituitary gland (adenohypophysis) secretes six known hormones that regulate homeostasis, including adrenocorticotropic hormone (ACTH), growth hormone (GH), prolactin, thyroid-stimulating hormone (TSH), follicle-stimulating hormone (FSH) and luteinizing hormone (LH) [39]. Clinically nonfunctioning PAs are not able to synthesize and/or secrete functional hormones, although in a majority of cases they nevertheless demonstrate positive immunostaining for LH, FSH, and/or the alpha-subunit [40]. PAs are monoclonal in origin and typically benign tumors, suggesting they arise from expansion of single precursor cells that possess a unique proliferative advantage [41]. These monoclonal adenomas therefore secrete specific hormones reflective of their differentiated cell of origin. In the current study, hierarchical clustering analysis readily separated somatotroph (GH) adenomas from nonfunctional adenomas, suggesting a potential role of DNA methylation in the differentiation and/or functional regulation of these tumors. Furthermore, an unexpected and interesting finding in the current study was that the functional corticotroph adenoma causing Cushing's disease and the silent corticotroph adenoma clustered to the same hierarchical group (group 2).
The exact mechanisms by which nonfunctioning PAs retain immunostaining but fail to produce clinically-active hormones remains to be determined. Our genome-scale screening of variability in DNA methylation between functional and nonfunctional PAs suggests that nonfunctional PAs are globally hypermethylated compared to functional ones. This finding is particularly interesting in that the most differentially methylated gene, KCNAB2, encodes a potassium ion-channel that has been previously implicated in endocrine function pertaining to insulin secretion [42]. Pituitary cells resemble neurons and muscle fibers in that they also fire action potentials (APs) [43], which are mediated via expression of numerous voltage-gated sodium, potassium, calcium and chloride channels. These APs are accompanied by a rise in intracellular calcium and spontaneous electrical activity that drives intracellular calcium concentrations above the threshold for stimulus-secretion and stimulus-transcription coupling. KCNAB2 is a subunit of the shaker-related voltage dependent potassium channel, and upon binding to K + channel alpha subunits, contributes to regulation of channel excitability. In addition, KCNAB2 has been demonstrated as an functional aldoketoreductase (AKR) [44,45]. Among the primarily hypermethylated genes in NFAs, we also identified four significant genes PSEN2 [46], MDN1 [47], PCSK6 and NME9 [48], which were relatively hypormethylated in NFAs compared to FAs (Table 2). Importantly, PCSK6 has been reported to be a member of the mammalian subtilisin-like proprotein convertase family, which participates in maturation of precursor protein, and is expressed at high levels in the anterior pituitary gland [49]. This link potentially explains the significant DNA methylation differences in PCSK6 methylation between NFA and FA, and suggests that PCSK6 may therefore be an important biomarker for distinguishing FA from NFA. Collectively, these data suggest that genome-scale DNA methylation analysis may provide a practical, complementary molecular correlate to standard PA classification schemes. Further analysis in a larger sample size is required to support the validity of these findings.
Additionally, our study demonstrates that the ion-channel activity signal pathway shows a relative hypermethylated profile in nonfunctional PA subjects. In particular, the pivotal gene KCNMB4 encodes an important subunit of potassium channels and is expressed throughout the brain (especially in the thalamus and the brainstem) [50]. In addition, transcription of the calcium channel gene CACNA1C has been reported to undergo epigenetic regulation via DNA methylation [51]. Taken together, these findings suggest that modulation of ion activity pathways may contribute to defective hormone secretion in nonfunctional PAs, and possibly serve as novel candidate biomarkers in PAs.
Finally, in order to investigate the association of DNA methylation and gene expression, total RNA sequencing and differential expression analysis was performed. Although the expression of KCNAB2 and other genes in ion-channel activity signal pathway was not significantly different between FAs and NFAs, our global correlation analysis showed that DNA hypermethylation was negatively associated with gene expression. The elevated expression of ODAM in FAs provides a potentially novel insight to the pathogenesis of epithelial neoplasms, since ODAM is a developmental antigen with an essential role in tooth maturation and the pathogenesis of various epithelial neoplasms [52]. Although we explored differences in RNA expression of KCNAB2 between NFAs and FAs, the result was not significant. A larger sample size would likely be needed to validate this association. Because gene expression is modulated not only by epigenome modification, low levels of gene expression cannot be attributed solely to hypermethylation. Additionally, the intratumoral multifocality and heterogeneity of DNA methylation of the PAs may account for this abberant association of DNA methylation and gene expression [53]. There are several limitations to the current study, including a relatively small sample size and nonuniform representation of functional to nonfunctional PAs and various histopathological PA subtypes. Furthermore, MethyLight validation of only selected candidate genes was performed, mandating further validation in a larger sample size and using another technique. The mechanisms by which methylation alterations of selected candidate genes result in observed invasion or hormonal phenotypes were also not studied, and remain beyond the scope of the current study. Finally, the current study is limited in that no true control samples (normal pituitary gland) were included in the analysis. Nevertheless, the current study suggests that DNA methylation analysis can be used to provide valuable insight into the phenotype of histopathological subtype in PAs, and if validated may complement current pathological classification systems. The current study is the first genome-wide DNA methylation analysis in PAs, and can be used to appropriately design and power future studies with a larger sample size in order to validate many of the preliminary findings from our study.
In conclusion, genome-scale DNA methylation profiling and RNA sequencing of PAs identified DNA methylation variations in candidate genes associated with functional subtype (both globally and in selected genes) and possibly invasion. Hierarchical clustering analysis showed PA clustering according to functional status and immunohistochemical subtype, suggesting that DNA methylation analysis may possibly provide a clinically useful and complementary molecular correlate to standard PA classification. Differential methylation of cell motility related genes (such as FLT1 and SLIT3) require further validation prior to being considered candidate biomarkers for PA invasion. DNA hypermethylation of KCNAB2 and enrichment of DNA methylation in ion-channel activity signal pathways may be associated with the endocrine-inactive status of nonfunctional PAs.

Supporting Information
File S1 Analysis strategies. Figure A, Filtering strategy for Illumina HM450 Probes. SNP10: SNPs within 10 bp from CpG site; SNP15: SNPs within 15 bp from the CpG site lying entirely within a repeat region; ChrX: chromosome X. Figure B, Additional filters were used in the secondary analysis of the invasive phenotype of NFAs. A total of 884 probes were enrolled in further analysis, of which 355 were enhancer (Enh)-associated, 20 were promoter (Pro)-associated, and 4 were both (Bot) promoter-and enhancer-associated, and we identified 34 significant (Sig) CpGs that were independently associated with enhancers and were hypomethylated in invasive NFAs compared to noninvasive NFAs. (TIF) File S2 The distribution of probes on the Illumina HM450 DNA methylation platform. Figure A, The distribution of 383,718 HM450 probes stratified by CpG density (in and out of CpG islands). Figure