Identification of key genes in invasive clinically non-functioning pituitary adenoma by integrating analysis of DNA methylation and mRNA expression profiles

Background Tumor surrounding the internal carotid artery or invading to the cavernous sinus is an important characteristic of invasive pituitary adenoma, and a pivotal factor of tumor residue and regrowth. Without specific changes in serum hormone related to the adenohypophyseal cell of origin, clinically non-functioning pituitary adenoma is more likely to be diagnosed at invasive stages compared with functioning pituitary adenoma. The underlying mechanism of tumor invasion remains unknown. In this study, we aimed to identify key genes in tumor invasion by integrating analyses of DNA methylation and gene expression profiles. Method Genome-wide DNA methylation and mRNA microarray analysis were performed for tumor samples from 68 patients at the Beijing Tiantan Hospital. Differentially expressed genes and methylated probes were identified based on an invasive vs non-invasive grouping. Differentially methylated probes in the promoter region of targeted genes were assessed. Pearson correlation analysis was used to identify genes with a strong association between DNA methylation status and expression levels. Pyrosequencing and RT-PCR were used to validate the methylation status and expression levels of candidate genes, respectively. Results A total of 8842 differentially methylated probes, located on 4582 genes, and 661 differentially expressed genes were identified. Both promoter methylation and expression alterations were observed for 115 genes with 58 genes showing a negative correlation between DNA methylation status and expression level. Nineteen genes that exhibited notably negative correlations between DNA methylation and gene expression levels, are involved in various gene ontologies and pathways, or played an important role in different diseases, were regarded as candidate genes. We found an increased methylation with a decreased expression of PHYHD1, LTBR, C22orf42, PRR5, ANKDD1A, RAB13, CAMKV, KIFC3, WNT4 and STAT6, and a decreased methylation with an increased expression of MYBPHL. The methylation status and expression levels of these genes were validated by pyrosequencing and RT-PCR. Conclusions The DNA methylation and expression levels of PHYHD1, LTBR, MYBPHL, C22orf42, PRR5, ANKDD1A, RAB13, CAMKV, KIFC3, WNT4 and STAT6 are associated with tumor invasion, and these genes may become the potential genes for targeted therapy.


Background
Pituitary adenomas (PA) are benign neuroendocrine tumors arising from adenohypophyseal cells and account for 10-15% of all primary intracranial neoplasms [1,2]. Compared with functioning pituitary adenoma (FPA), clinically non-functioning pituitary adenoma (NFPA) often shows no specific serum hormone changes [3]. When the mass effect appears, it is possible that the tumor may have surrounded the internal carotid artery and invaded the cavernous sinus, which is the main cause of residual tumor and postoperative regrowth [4].
The overall rate of tumor invasion into the cavernous sinus is 35%, and the understanding of the mechanisms underlying the pituitary adenoma invasiveness is still far from comprehensive [5]. Nevertheless, studies have indicated that the pituitary tumorigenesis and tumor invasiveness may not be driven by gene mutations [6,7] but rather by other regulatory mechanisms. Abnormalities of DNA methylation status of cytosine phosphate guanosine (CpG) islands in gene promoter regions are recognized to play a key role in regulating transcription and are closely associated with diverse diseases [8][9][10][11]. For example, some studies have reported methylation status alterations in NFPA tumorigenesis, but failed to uncover the changes in methylation involved in tumor invasiveness [12,13]. Although these studies have identified the methylation changes in NFPA, limited work have been performed to connect these changes with whole-transcriptome changes. Furthermore, the blueprint of epigenetics cooperating with the transcriptome in NFPA has not been fully evaluated.
In this study, integrated analyses of paired wholegenome DNA methylation and gene expression microarray were performed. We found eleven genes that play an important role in the invasive behavior of NFPA. Pyrosequencing and RT-PCR were then utilized to validate the DNA methylation status and gene expression level, respectively. We aimed to identify gene expression epigenetically regulated by methylation changes in invasive NFPA, and hope to obtain a better understanding of the invasive behavior of NFPA.

Patients and samples
This study retrospectively enrolled 68 patients aged 25-74 years who were diagnosed with NFPA from October 2007 to July 2016. All patients underwent enhanced head MRI scan before and after trans-sphenoidal surgery in order to assess the maximum tumor diameter, anatomical location and tumor resection grade. Tumor invasion was defined by tumor extension beyond the lateral tangent of the intra-and supra-cavernous internal carotid artery (Grade 3 and 4) as defined by Knosp et al. [14].
Only patients who were clinically diagnosed with and pathologically confirmed as having gonadotroph adenoma and null-cell adenoma were included in the study. All tumor samples were collected from the Neurosurgery Department of the Beijing Tiantan Hospital. Tumor samples were immediately placed into a sample tube, frozen in liquid nitrogen and stored. A total of 68 pituitary adenoma samples, including 46 invasive and 22 non-invasive tumors, were used for following whole-genome DNA methylation and mRNA microarray analysis.

Whole-genome DNA methylation microarray
DNA was extracted and purified with DNeasy Blood & Tissue Kit (Qiagen, Germany). The bisulfite conversion of the DNA was performed using EZ DNA Methylation-Gold ™ Kit (Zymo, USA). The total DNA methylation status of 68 NFPA samples was assayed by Illumina Infinium MethylationEPIC 850K BeadChip. The DNA methylation microarray experiment was performed at Shanghai Biotechnology Corporation following the manufacturer's instructions of Illumina. Probes located on the sex chromosomes, bound to multiple genomic locations or associated with a known SNP were excluded [15]. Probes that failed to be detected above background were also removed from the research. The bio-conductor R package minfi was used to control quality and normalize the raw data [16]. The methylation status of different CpG sites was calculated with the average-difference β-value (Δβ), where beta-value (β) is between 0 (unmethylated) and 1 (completely methylated). Differentially methylated probes (DMPs) that showed a |∆β| of 0.1 and an adjusted p-value < 0.05 were regarded as significantly differentially methylated. Then a separate average DNA methylation β-value of the DMPs in the gene promoter and non-promoter region was employed to represent the methylation status, respectively.

Whole-genome mRNA microarray
RNA was isolated and purified with the mirVana ™ miRNA Isolation Kit (Ambion, USA), then amplified and labeled with Low Input Quick Amp WT Labeling Kit (Agilent Technologies, USA), following the manufacturer's instructions. Then, purified RNA was used to generate fluorescence-labeled cRNA targets for SBC human ceRNA array V1.0 (4 × 180 K). The labeled cRNA targets were then hybridized and scanned with an Agilent Microarray Scanner (Agilent Technologies, USA). The data were extracted with Feature Extraction software 10.7 (Agilent Technologies, USA), and the R package limma was used to normalize the raw data [17]. The microarray experiments were performed at Shanghai Biotechnology Corporation following the protocol of Agilent Technologies. Differential gene expression was analyzed with a t-test in limma of the R package [18]. The fold change method was applied to estimate the differential significance of mRNAs. In our study, invasion-related differentially expressed genes (DEGs) in 68 NFPA samples were defined as the genes with a fold change > 1.5 or < 0.67, a p-value < 0.05 and a FDR value < 0.25.

Integrated analysis
A two-step analysis process was performed to integrate promoter DNA methylation in the promoter region with gene expression. First, we located all DMPs (unfiltered) for targeted genes to identify differentially methylated genes (DMGs) between invasive and non-invasive tumors. Second, for the both differentially methylated and differentially expressed genes, we used Pearson correlation analysis to examine whether there was a strong association between the DNA methylation status and expression levels. In this study, we applied Pearson correlation coefficient (PCC) and p-values to investigate the correlation and significance between mRNA expression and DNA methylation with the function cor.test in R. A significantly negative correlation was considered if PCC < − 0.2 and p-value < 0.05.
To infer potential biological processes and pathways of invasion-related genes, the DAVID Bioinformatics Tool (v6.8) was employed to perform functional enrichment analysis using the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG). Biological processes and pathways were considered significant at p-value < 0.05.

Pyrosequencing assay
Pyrosequencing assay was used to access the level of DNA methylation of the promoter regions of target genes. Genomic DNA from NFPA samples was extracted with QIAamp DNA Mini Kit (Qiagen, Germany). For each sample, 0.5 μg DNA was used to for bisulfite conversion with the EpiTect Bisulfite kit (Qiagen, Germany). 1 μl of 10 μl of eluted bisulfite converted DNA was used for performing PCR with PyroMark PCR Kit (Qiagen, Germany), according to the manufacturer's instructions. The methylation status of each gene was assessed as the percentage of average methylation at targeted CpG sites. The PCR primer designed for pyrosequencing is shown in Additional file 1: Table S8.

Real-time quantitative reverse transcription polymerase chain reaction (RT-PCR)
Extraction and purification of total RNA were carried out as described above. RT-PCR analyses were performed with PrimerScript RT reagent Kit (Takara, China). Quantitative real-time PCR was performed using LightCycler 480 SYBR Green I Master (Roche, Switzerland) with the ABI 9700 PCR system (Applied Biosystems, USA). Reactions were run at 95 °C for 10 min followed by 40 cycles of 95 °C for 10 s and 60 °C for 30 s. Gene expression was measured by the standard curve method and normalized to the level of β-actin; the 2 −ΔΔct method were used for calculations [19]. The PCR primer sequences are shown in Additional file 1: Table S9.

Statistical analysis
Statistical analyses were performed in R (version 3.1) and IBM SPSS Statistics for Windows (v 21.0). Student's t test was applied for comparison of two groups and for pyrosequencing and RT-PCR assay, respectively. All data are presented as the mean ± standard deviation (SD). Differences were considered to be statistically significant if p < 0.05.

Whole-genome DNA methylation analysis
The flowchart of the study is shown in Fig. 1. Patients were separated into invasive and non-invasive groups based on the criteria mentioned above. The clinical characteristics of these 68 patients are shown in Table 1 and detailed information of each patient is provided in Additional file 1: Table S1. We compared the DNA methylation patterns in those two groups by whole-genome DNA methylation microarray (invasive vs non-invasive tumor). After filtering the raw data (602,698 probes remained) and performing statistical analysis, 8842 probes that showed significant changes with |∆β| > 0.1 between invasive and non-invasive groups were included for downstream analyses.
Distributions of DMPs were analyzed based on their genomic regions as well as the locations relative to CpG islands. Among all 8842 DMPs, 22% were found to be located in the TSS1500 (transcription start sites 1500), 7% in the TSS200 (transcription start sites 200), 12% in the 5′ UTR (5′ untranslated region), 4% in 1st exon, 52% in gene body and 3% in the 3′ UTR (Fig. 2a, Additional file 1: Table S10). When considering distribution of DMPs relative to the location of the CpG island, 52% were located in the open sea, 13% in the CpG S_Shore, 2% in the CpG S_Shelf, 14% in the CpG island, 16% in the CpG N_Shore and 3% in the CpG N_Shelf (Fig. 2b, Additional file 1: Table S10).
Some differences between overall DNA methylation levels for invasive and non-invasive tumors were observed when the average β-values were compared according to the gene-related regions and annotations to the CpG islands, respectively. Significant variances in DNA methylation level were found in TSS1500, TSS200, 1st exon, 5′ UTR and gene body between the two groups, though there was no significant difference in 3′ UTR region (Fig. 2c). When focusing on the diversity of DNA methylation level according to CpG island location, significant differences were observed in the open sea, S_Shore, island, N_Shore and N_Shelf but not for the S_Shelf (Fig. 2d).
Among all 8842 DMPs, there were 6062 hypermethylated and 2780 hypomethylated probes between invasive and non-invasive tumors ( Fig. 3a and Additional file 1: Table S2). The heatmap showed the methylation status of all significant probes based on β-values of these CpG sites across all 68 samples (p < 0.05, |∆β| > 0.1, Fig. 3b). We then related the DMPs to specific genes according to their location. A total of 8842 DMPs were related to 4582 genes; for 1904 genes, the DMPs were in the promoter region, and for the remaining 2678 genes, the DMPs were in the non-promoter region (Additional file 1: Table S3).

Identification of differentially expressed genes
The differential mRNA expression status was also assessed between invasive and non-invasive tumors. A total of 661 genes were differentially expressed (p < 0.05, FDR < 0.25 |log 2 fold change| > log 2 1.5) between two groups, of which, 206 genes were found to be upregulated and 455 genes were downregulated ( Fig. 3c and Additional file 1: Table S4). The heatmap and volcano plot of the DEGs are presented in Fig. 3d.

Integrated analysis of DMGs and DEGs
The DNA methylation changes in the gene promoter region have been demonstrated to be closely associated with gene expression regulation [20,21]. However, the connection between non-promoter region DNA methylation and gene expression is still in controversial [22][23][24][25]. Thus, our ensuing analyses only focused on the DMGs in promoter region with alterations in DNA methylation.
Using the differentially expressed genes based on DMG analyses, we observed 115 of the 661 genes to be both differentially methylated and expressed (Fig. 4a). Enrichment analysis showed that these genes are related to the plasma membrane, integral component of membrane and transmembrane signaling receptor activity etc. (Fig. 4b, Additional file 1: Table S6).
We further studied correlations of the methylation status of the 115 genes and their expression levels using Pearson correlation analysis and found 58 genes with a negative correlation (p < 0.05, R < − 0.2, Fig. 4c, Additional file 1: Table S5). Enrichment analysis of these 58 genes also revealed the relevance of tumor oncogenesis (Fig. 4d, Additional file 1: Table S7).

Validation of the methylation and expression levels of the candidate genes
Genes with notably negative correlations between DNA methylation and gene expression levels, are involved in various gene ontologies and pathways, or play an important role in different diseases were regarded as candidate genes. Based on the above criteria, we selected 19 genes to validate the invasive behavior of NFPA and performed pyrosequencing and RT-PCR to validate their DNA methylation and expression levels. A significant increase in DNA methylation levels in the invasive tumors compared with non-invasive tumors was observed for the following eleven genes: PHYHD1, LTBR, C22orf42, PRR5, ANKDD1A, RAB13, CAMKV, KIFC3, WNT4, STAT6, GBGT1 and GALNT14 (Fig. 5a, b, d-k, Additional file 2: Figure S1e, f ). Decreased methylation levels were observed in MYBPHL (Fig. 5c). Conversely, no significant DNA methylation changes were found for SLC12A4, PDE9A, PDLIM4, CHD5, LHX3 and BATF2 (Additional file 2: Figure S1a-d, g, h).

Discussion
Unlike FPA, such as somatotroph adenoma and lactotroph adenoma, there is no targeted medicine for clinical therapy of NFPA, and surgery has been suggested to be the most effective therapeutic approach [5]. Tumor surrounding the internal carotid artery or invading the cavernous sinus are the most frequent factors of a residual tumor, which could be an important risk factor for tumor regrowth [2]. For the above reasons, we focused on alterations in DNA methylation status and corresponding gene expression changes between invasive and noninvasive NFPA, in order to identify key genes involved in invasive behavior and potential candidates for targeted therapy. By integrating whole-genome DNA methylation and gene expression analyses, we identified eleven genes that may correlate with tumor invasion, which may facilitate an understanding of underlying mechanism of pituitary adenoma invasion.
Although germline and somatic mutations are thought to be related to the tumorigenesis of some subtypes of pituitary adenoma, the drivers of NFPA tumorigenesis remain unknown, as are the underlying mechanisms of tumor invasion [26]. PIK3CA mutations have been identified in pituitary adenomas, including lactotroph adenoma, corticotroph adenoma and NFPA. Indeed, such mutations are observed in 8.8% of invasive pituitary tumors, whereas no mutations have been detected in noninvasive tumors [27]. Other studies have found that some genes are closely related to the invasive behavior of pituitary. ESM-1 encodes endocan, which is highly expressed in pituitary adenomas and thought to be related to angiogenesis and strongly associated with tumor invasion in pituitary adenoma [28]. Additionally, Yang et al. [29] found that Caveolin-1 (Cav-1), a principal structural protein of caveolae, promotes pituitary adenoma cell migration and invasion by regulating the interaction between EGR1 and KLF5. Furthermore, ADAM12 overexpression is associated with tumor invasion in pituitary adenoma via the EGFR/ERK signaling pathway [30]. However, to some extent, these studies illustrate the mechanism of tumor invasion from a genetic standpoint, the underlying mechanism has not been fully elucidated. Previous studies have also illustrated the strong association between pituitary adenoma tumorigenesis or invasion and regulation of DNA methylation. Duong et al. [31] used infinium methylation 27 K arrays to examine DNA methylation in various subtypes of pituitary adenoma and found twelve genes showing methylation alterations, of which three genes, EML2, RHOD and HOXB1, also exhibited significantly decreased expression levels. Ling et al. [32] used genome-scale profiling of DNA methylation in FPA and NFPA and attempted to explore associations with tumor invasion and histopathological subtype. They calculated the mean DNA methylation beta value for the entire filtered probe and did not find significant differences in DNA methylation levels. However, methylation alteration in promoter region is able to affect gene expression and gene body methylation is not always consistent with promoter methylation alteration. As a result, the average beta value of an entire gene may disguise the methylation changes that genuinely affect the gene expression. Although they found that KCNAB2 may contribute to the endocrine-inactive status of NFPA, no DNA methylation changes in were detected between invasive and non-invasive NFPA. Compared with their research, our study focused on the methylation changes in the promoter region between invasive and non-invasive tumor and successfully find eleven genes by integrating analyses of DNA methylation and expression microarray. Some researchers have tried to seek the specific DNA methylation changes in clinically non-functioning pituitary adenoma. For instance, by comparing the DNA methylation profiles between normal pituitary tissue, non-invasive and invasive pituitary adenomas, Kober et al. [33] found that promoter hypermethylation and decreased expression levels of five genes in NFPA compared with normal pituitary. However, differences in the methylation profiles between invasive and noninvasive NFPA were not detected. Additionally, hypermethylation of GALNT9 was found to be significantly downregulated in invasive NFPA, possibly affecting cell adhesion [34].
Differential analysis of all 853,307 probes between invasive and non-invasive tumors enabled us to study the methylation pattern from a more comprehensive perspective. Kober et al. [33] used Illumina HM450K BeadChip to assess the DNA methylation status in NFPA (n = 34) and normal pituitary (n = 4). Meanwhile, they divided 34 NFPA into invasive (n = 18) and non-invasive (n = 16) to explore the DNA methylation alterations. Their research finds methylation and expression differences between NFPA and normal pituitary, which may explain the pathogenesis of NFPA. However, they do not find the methylation differences between invasive and non-invasive NFPA. In our study, we focused on the methylation differences between invasive (n = 46) and non-invasive (n = 20) NFPA and we found hypermethylation appeared to be a more frequent phenomenon in invasive tumors with regard to the distribution of DMPs in genomic regions and the locations of CpG islands. In our study, we used Illumina HM850K BeadChip which contains much more information than HM450K Bead-Chip, and our chip could find more potential targets. Also, Kober's research did not perform a whole-genome transcriptome analysis. After they selected target genes from methylation microarray analysis, qRT-RCR was performed to get the expression level of target genes. In our study, we simultaneously perform whole-genome methylation and whole-genome transcriptome analysis which is more comprehensive.
Changes in DNA methylation can be observed in thousands of genes, and only genes with methylation changes that influence expression were considered as candidates. Based on this criterion, we performed whole-genome DNA methylation and gene expression analyses. In this study, we selected key genes showing notably negative correlations between DNA methylation and gene expression levels, are involved in various gene ontologies and pathways, or played an important role in different diseases. We identified PHYHD1, LTBR, MYBPHL, C22orf42, PRR5, ANKDD1A, RAB13, CAMKV, KIFC3, WNT4 and STAT6 as having alterations in DNA methylation and gene expression in invasive NFPA. Most of our genes were reported in large data analysis research, but the significance of methylation changes has not yet been reported. LTBR plays a role in lymphoid development signaling and the immune response, and the hypermethylation of LTBR is observed in odontogenic keratocysts, though its expression is not affected [35]. The relationship between methylation and gene expression is complicated and it can be affected by many factors. In our study, we observed a negative correlation between the hypermethylation of LTBR and its downregulated expression. In addition, hypermethylation of ANK-DD1A in glioblastoma has been reported; this gene interacts with FIH1 and decreases HIF1-α stability to inhibit cell autophagy in a hypoxic microenvironment [36,37], but its function in pituitary adenoma or the significance of promoter methylation is still unclear. RAB13 has been predominantly studied in epithelial cells and its activity regulates tight junction assembly and stimulates cell invasion and migration. Regardless, the methylation regulation of RAB13 has not been reported [38,39]. A large-scale transcriptomic study found CAMKV to be present in the synaptic neuropil, and quantitative proteomics analysis revealed that CAMKV is downregulated at the synapse upon sensory deprivation, though changes in methylation in tumors have not been explored yet [40,41]. The overexpression of WNT4 has been observed in various subtypes of pituitary adenoma, but it displays an inverse correlation with tumor invasion [42]. In the present study, decreased expression of WNT4 was observed and our results showed that DNA methylation may play a role in its regulation. In general, alterations in methylation of most of the genes identified in our study have not been examined in pituitary adenoma or other solid tumors, and further research on their function and mechanism in tumor invasion is essential.
Promoter methylation alterations have been proved to be closely related to gene expression. In pituitary adenoma, the methylation alteration of LAMA2, GALNT9 and DAP methylation has been proven to be related to tumor invasion [34,43,44]. In our study, alteration in promoter DNA methylation of eleven genes was associated with NFPA invasion. However, recent studies have reported that the methylation changes in gene body regions may also result in the gene expression changes. Yang et al. [25] found that gene body DNA methylation increases gene expression and this regulation is dependent on the presence of the DNMT3B. Moreover, hypermethylation in the gene body region accompanied by upregulated expression was also observed in a hepatocellular carcinoma mouse model [22], and Chen et al. [23] found gene body hypermethylation to be significantly associated with silencing of the tumor-related genes in kidney cancer. Though whole-genome bisulfite sequencing analysis of multiple individuals, Lou et al. [24] revealed that gene body hypermethylation leads to reduced transcription efficiency. Regardless, the correlation between gene body methylation and gene expression is still controversial. Therefore, we selected genes that showed a negative correlation between promoter methylation and gene expression and did not further investigate differential probes located in gene body regions. This may be one limitation of our study, and our future research may focus on the regulation mechanism between gene body methylation and expression in NFPA.

Conclusions
In conclusion, this study integrates DNA methylation and gene expression and successfully identified key genes in invasive NFPA. We found alterations of DNA methylation in the promoter region and expression changes for eleven genes between invasive and non-invasive NFPA. Our study may promote the understanding of NFPA invasion and facilitate targeted therapy for patients with invasive tumors.
Additional file 1: Table S1. Clinical information of 68 patients with NFPA. Table S2. List of differentially methylated probes. Table S3. List of differentially methylated genes in the promoter region. Table S4. List of differentially expressed genes. Table S5. Pearson analyses results of 69 significant genes. Table S6. GO and KEGG analyses of 115 genes. Table S7. GO and KEGG analyses of 58 genes. Table S8. Pyrosequencing primers of 19 genes. Table S9. RT-PCR primers of 20 genes (including thr control). Table S10. Distribution of differential methylated probes.