Bioinformatic analysis of gene expression profiles of pituitary gonadotroph adenomas.

The aim of the present study was to identify genes, microRNAs (miRNAs/miRs) or pathways associated with the development of pituitary gonadotroph adenomas. The array data of GSE23207, which included 16 samples of multiple endocrine neoplasia-associated rat pituitary homozygous mutations and 5 pituitary tissue samples from healthy rats, were downloaded from the Gene Expression Omnibus database. Differentially expressed genes (DEGs) were analyzed prior to functional enrichment analysis and protein-protein interaction (PPI) network construction. miRNAs associated with DEGs were predicted, and an miRNA-target regulatory network was constructed. A total of 187 upregulated and 370 downregulated DEGs were identified in the pituitary gonadotroph adenoma group compared with the healthy (control) group. Cyclin-dependent kinase (Cdk) 1 exhibited the highest degree in the PPI network. The upregulated DEGs were predominately enriched in 'neuroactive ligand-receptor interaction' pathway, and downregulated DEGs were mainly enriched in 'cell cycle'. The DEGs in module were predominately enriched in the 'cell cycle', whereas DEGs in module b and c were enriched in 'neuroactive ligand-receptor interaction'. miR-374, -153, -145 and -33 were identified as important miRNAs in the regulation of the DEGs. Cdk1, cyclin (Ccn) A2, Ccnb1, 'cell cycle' and 'neuroactive ligand-receptor interaction' pathways may serve important roles in the development of pituitary gonadotroph adenomas; Ccna2 and Ccnb1 may contribute to this development via an effect on the 'cell cycle' pathway. Furthermore, miR-374 and -145 may contribute to the development of pituitary gonadotroph adenomas via regulation of the expression of target genes.

Introduction

Pituitary adenomas are generally benign, arise from adenohypophyseal cells and are one of the most frequent types of intracranial tumor (1). It is reported that pituitary gonadotroph adenomas account for 15–22% of all pituitary adenomas and include the majority of clinically nonfunctional pituitary adenomas (2,3). Surgery remains the first line of treatment for the majority of patients, but it cannot effectively control invasive pituitary adenomas (4). Therefore, novel and effective therapeutic approaches are required. Understanding the molecular mechanisms of pituitary gonadotroph adenoma development may provide new insights for their treatment.

Heaney et al (5) proposed peroxisome proliferator-activated receptor-γ ligands as candidates for the management of nonfunctioning pituitary tumors. Chesnokova et al (6) indicated that Forkhead box gene transcription factor L2 (FOXL2) activated the clusterin promoter in gonadotroph pituitary cells. Chesnokova et al (7) reported that clusterin and FOXL2 regulated the growth of pituitary gonadotroph adenoma. Lee et al (8) suggested that somatostatin receptor 3 may be a possible target for the treatment of pituitary gonadotroph adenomas. Furthermore, the downregulation of miRNA targeting high mobility group AT-hook 1 (HMGA) 1 and 2 and E2F transcription factor 1 may contribute to pituitary tumorigenesis (9). The downregulation of miR-23b and miR-130b expression may contribute to pituitary tumorigenesis (10). In addition, targeting phosphoinositide 3-kinase/mechanistic target of rapamycin signaling may activate antitumor effects against nonfunctioning pituitary adenomas (11). Despite results demonstrated in these previous studies, knowledge about the underlying molecular mechanisms of pituitary gonadotroph adenoma development may be insufficient and further research is required.

Lee et al (12), the contributors of the GSE23207 microarray dataset, demonstrated that the multiple endocrine neoplasia (MENX) rat model could be used as an experimental tool to study the pathological mechanisms for human pituitary tumorigenesis. Zhang et al (13) used the GSE23207 dataset to identify that differentially expressed genes (DEGs) associated with ‘cell cycle’, ‘cell division’, ‘neuroactive ligand-receptor interaction’, ‘pituitary gland’, ‘adenohypophysis’ and ‘endocrine system’ may serve important roles in the pathogenesis of pituitary adenomas via DEG screening, gene ontology (GO) and pathway enrichment analysis, and protein-protein interaction (PPI) network construction. The present study used the microarray GSE23207 dataset in addition to the above techniques; however, miRNA-DEG regulatory network analysis was also performed. Important genes, miRNAs and pathways associated with the development of pituitary gonadotroph adenomas were identified, in order to aid the clarification of the molecular mechanisms of pituitary adenomas.

Materials and methods

Microarray data

The raw microarray data from GSE23207, as deposited by Lee et al (12), and GPL6247, as produced with the Affymetrix Rat Gene 1.0 ST array [transcript (gene) version; Thermo Fisher Scientific, Inc., Waltham, MA, USA] were downloaded from the Gene Expression Omnibus (GEO; http://www.ncbi.nlm.nih.gov/geo) database. From these datasets, data from a total of 21 pituitary samples, including 16 from MENX rats with homozygous mutations (p27Kip1/Cdknb1) and 5 healthy rat pituitary tissue samples were used in the present study.

Data preprocessing

The processing of the raw data, including format conversion, the supplementation of missing values, background correction and quartile normalization, was performed using the affy package in R (14).

Screening of DEGs

The DEGs in the pituitary gonadotroph adenoma group compared with the control group were screened with the limma package (15). The P-values for DEGs were calculated by a t-test in the limma package. Then, the P-values were adjusted to false discovery rate (FDR) values using the Benjamini-Hochberg procedure (16). |log2 fold change (FC) |≥1 and FDR <0.05 were set as cut-off criteria for all DEGs. The heatmap for DEGs was drawn with the pheatmap package in R (17).

Functional enrichment analysis tools

Gene ontology (GO) annotation (18), including the categories of molecular function (MF), biological process (BP) and cellular component (CC), and Kyoto Encyclopedia of Genes and Genomes (KEGG) (19) pathway enrichment analyses were performed to identify upregulated and downregulated DEGs. The functional enrichment analysis was performed using the Multifaceted Analysis Tool for Human Transcriptome (MATHT; www.biocloudservice.com). The Gene Set Function-Functional Enrichment-mRNA Enrichment module, based on Fisher's test, was used for the enrichment analysis of gene sets with a cut off of P<0.05.

PPI network

The Search Tool for the Retrieval of Interacting Genes database was used to predict and analyze the interactions of proteins encoded by DEGs (20). PPI networks were constructed using Cytoscape software (version 3.2.0) (21). The input gene sets were the identified DEGs, and the species was set to Rattus norvegicus. DEGs were constructed into a PPI network with a weight threshold of ≤0.4 (medium confidence). The top 20 genes were added one-by-one according to degree value from high to low, starting with Cdk1, in the cluster analysis until the top 20 genes clustered correctly. Clustering analysis was performed using molecular context detection Cytoscape plug-in (MCODE; version 1.4.1; http://baderlab.org/Software/MCODE).

miRNA-target genes regulatory network analysis

The miRNA that potentially regulated the DEGs were predicted using WebGestal (22,23), and a miRNA-target regulatory network was constructed using Cytoscape. Experimentally verified miRNA-DEG interactions were downloaded from Mirwalk2 (24).

Verification

The GSE26966 dataset, including 10 gonadotroph tumors and 9 normal human pituitaries from autopsy samples, were downloaded from the GEO database. Data normalization and DEG analysis was performed with the affy and limma packages, respectively. |log FC |>1 and P<0.05 were used as cut-off criteria. GO-BP and KEGG pathway enrichment analyses were performed using the Database for Annotation, Visualization and Integrated Discovery (25). FDR<0.05 was set as a cut-off value. In addition, the DEGs in GSE26966 were compared with the DEGs in GSE23207.

Results

DEGs analysis

In total, 557 DEGs, including 187 significantly upregulated and 370 significantly downregulated DEGs, were identified in the pituitary gonadotroph adenoma group compared with the control group (Fig. 1).

Figure 1.

Heatmap and hierarchical clustering analysis of differentially expressed genes between the GSE23207 (PL) and GPL6247 (NP) datasets.

Functional enrichment analysis

The results of functional enrichment analysis for the up- and downregulated DEGs are included in Fig. 2. The upregulated DEGs were predominately enriched in ‘neuroactive ligand-receptor interaction’ and the ‘GnRH signaling pathway’; downregulated DEGs were enriched in ‘cell cycle’.

Figure 2.

Functional enrichment analysis of up- and downregulated DEGs. (A) Enriched Kyoto Encyclopedia of Genes and Genomes pathways, including (i) upregulated and (ii) downregulated. (B) Encircled Gene Ontology terms, including the top 5 for (i) upregulated and (ii) downregulated DEGs. DEGs, differentially expressed genes; BP, biological process; CC, cellular component; MF, molecular function.

There were 413 nodes and 1,601 protein pairs in the PPI network (Fig. 3). The top 20 nodes with the highest degrees are included in Table I. The top 5 nodes were Cdk1 (degree, 52), topoisomerase (DNA) II α (Top2α; degree, 51), Ccna2; (degree, 50), Fraser extracellular matrix complex subunit 1 (Fras1; degree, 46) and Ccnb1 (degree, 46). The expression of the top 20 nodes from the PPI network is illustrated in Fig. 4.

Figure 3.

Protein-protein interaction network for DEGs. Circle, upregulated DEGs; diamond, downregulated DEGs. The higher the degree value, the larger the node. DEGs, differentially expressed genes.

Table I.

Top 20 nodes with the highest degrees in the protein-protein interaction network.

Node	Degree
Cdk1	52
Top2a	51
Ccna2	50
Fras1	46
Ccnb1	46
Plk1	45
Mapk12	44
Cdc20	44
Kif11	41
Racgap1	37
Bub1	37
Prc1	36
Bub1b	36
Kras	35
Kif20a	35
Kif2c	34
Vegfa	34
Dlgap5	33
Erbb4	32
Nuf2	32

Figure 4.

Heatmap and hierarchical clustering analysis for the top 20 nodes.

Network cluster module

A total of 3 cluster modules were obtained (A-C; Fig. 5). Module a contained 29 nodes and 365 edges, module b contained 9 nodes and 35 edges and module c contained 9 nodes and 16 edges. All genes in module a were upregulated, and all of the top 20 genes were included in this module. The KEGG pathways significantly enriched by these modules are included in Table II. The DEGs in module a were enriched in the ‘cell cycle’ and ‘oocyte meiosis’ pathways. The DEGs in module b were enriched in the ‘chemokine signaling pathway’, ‘calcium signaling pathway’, and ‘neuroactive ligand-receptor interaction’. DEGs in module c were enriched in ‘neuroactive ligand-receptor interaction’.

Figure 5.

Three network cluster modules. (A) Module A contained 29 nodes and 365 edges. (B) Module b contained 9 nodes and 35 edges. (C) Module C contained 9 nodes and 16 edges. Circle, upregulated DEGs; diamond, downregulated DEGs. DEGs, differentially expressed genes.

Table II.

The Kyoto Encyclopedia of Genes and Genomes pathways significantly enriched in the three identified modules.

A, enriched terms associated with module a
Pathway ID	Pathway name	Count	P-value	Genes
rno04110	Cell cycle	9	4.71×10−13	CCNB1, CCNB2, PLK1, BUB1, BUB1B, TTK, CDC20, PTTG1, CCNA2
rno04114	Oocyte meiosis	6	3.34×10−7	CCNB1, CCNB2, PLK1, BUB1, CDC20, PTTG1
rno04914	Progesterone-mediated oocyte maturation	5	6.50×10−6	CCNB1, CCNB2, PLK1, BUB1, CCNA2
rno04115	p53 signaling pathway	3	4.69×10−3	CCNB1, CCNB2, RRM2
B, enriched terms associated with module b
Pathway ID	Pathway name	Count	P-value	Genes
rno04062	Chemokine signaling pathway	4	5.26×10−4	AGTR1B, GRK5, PLCB1, CCL6
rno04020	Calcium signaling pathway	4	6.32×10−4	AGTR1B, TACR2, PLCB1, F2R
rno04080	Neuroactive ligand-receptor interaction	4	1.83×10−3	AGTR1B, TACR2, GNRHR, F2R
C, enriched terms associated with module c
Pathway ID	Pathway name	Count	P-value	Genes
Brno04080	Neuroactive ligand-receptor interaction	3	1.23×10−2	GABRA1, GABRA4, GABRB1

miRNA-target regulatory network analysis

In total, 4 miRNAs were identified (miR-374, −153, −145 and −33); miR-374 and −145 had been experimentally verified to be able to regulate DEGs (26,27). A total of 50 nodes and 60 edges were included in this regulatory network (Fig. 6). The KEGG pathway enriched by target genes in this network was ‘GnRH signaling’.

Figure 6.

The miRNA-target regulatory network. Triangle, miRNA; circle, upregulated DEGs; diamond, downregulated DEGs. miRNA, microRNA; DEGs, differentially expressed genes.

A total of 1,432 DEGs (including 651 upregulated, 781 downregulated) were identified in the pituitary gonadotroph adenoma group compared with the control group. In total, 13 common upregulated and 38 common downregulated DEGs were obtained from comparing the GSE26966 and GSE23207 data. The significantly enriched GO-BP and KEGG pathways for GSE26966 were also identified; the top 10 upregulated GO-BP terms were ‘ion transmembrane transport’, ‘positive regulation of synapse assembly’, ‘adherens junction organization’, ‘central nervous system development’, ‘membrane depolarization during cardiac muscle cell action potential’, ‘chloride transmembrane transport’, ‘heterophilic cell-cell adhesion via plasma membrane cell adhesion molecules’, ‘sodium ion transport’, ‘potassium ion transport’ and ‘negative regulation of the mitotic cell cycle’.

Discussion

A total of 557 DEGs, including 187 significantly upregulated and 370 downregulated DEGs, were identified in the pituitary gonadotroph adenoma group compared with the control group. Cdk1 (degree, 52) had the highest degree in the PPI network. Ccna2 (degree, 50) and Ccnb1 (degree, 46) were also enriched in the ‘cell cycle’. Furthermore, the upregulated DEGs were predominately enriched in ‘neuroactive ligand-receptor interaction’ and downregulated DEGs were mainly enriched in the ‘cell cycle’. DEGs in module were also mainly enriched in ‘cell cycle’. DEGs in module b and c were enriched in ‘neuroactive ligand-receptor interaction’. Furthermore, miR-374, −153, −145 and −33 were identified as miRNAs that may have regulated the DEGs.

Cdk1, Ccna2 and Ccnb1 demonstrated the highest connectivity degrees in the PPI network. A number of previous studies have suggested that Ccnb1 and Ccna2 are abnormally expressed in pituitary adenomas (10,28). Another study suggested that the downregulation of miR-410, which targets Ccnb1, contributes to pituitary gonadotroph tumor development (4). The overexpression of Ccna2 serves a key function in pituitary tumorigenesis (10). p27 belongs to the Cdk1 family, and is downregulated in invasive and recurrent adenomas as well as pituitary carcinomas (29). These studies suggest that Cdk1, Ccna2 and Ccnb1 may have important roles in the development of pituitary gonadotroph adenomas, which supports the results of the present study.

Furthermore, in this study, the upregulated DEGs were predominately enriched in ‘neuroactive ligand-receptor interaction’ and downregulated DEGs were enriched in ‘cell cycle’. DEGs in module a were also mainly enriched in ‘cell cycle’. DEGs in modules b and c were enriched in ‘neuroactive ligand-receptor interactions’. Previous research has demonstrated cell cycle dysregulation in pituitary adenomas (30–32). Quereda et al (31) indicated that the dysregulation of the cell cycle in pituitary disease had a significant effect on treatment. Cell cycle dysregulation is the main mechanism by which HMGA proteins cause the development of pituitary adenomas (33). The damage to cell cycle regulation in pituitary adenomas caused by Wee1 downregulation may cause the G2/M checkpoint to be lost, leading to the accumulation of DNA damage and the induction of tumor development (34). Thus, the dysregulation of the cell cycle may be associated with the progression of pituitary gonadotroph adenomas. In addition, Zhang et al (13) indicated that genes associated with neuroactive ligand-receptor interaction, including the γ-aminobutyric acid type a receptor α1, α4 and β1 subunits, may serve essential roles in the pathogenesis of pituitary adenomas. Although the role of neuroactive ligand-receptor interaction in pituitary adenoma development have not been fully characterized, the present study identified that the upregulated DEGs were predominately enriched in ‘neuroactive ligand-receptor interaction’, and that DEGs in module b and c were particularly enriched in ‘neuroactive ligand-receptor interaction’. Thus, it may be that a neuroactive ligand-receptor interaction is critical in the pathogenesis of pituitary gonadotroph adenomas.

Ccna2 and Ccnb1 were also enriched in ‘cell cycle’ in the present study. Therefore, Ccna2 and Ccnb1 may be associated with the development of pituitary gonadotroph adenomas via an effect on the cell cycle.

The results also identified that miR-374, −153, −145 and −33 may have been important miRNAs for regulating DEGs. miRNAs function as antisense regulators of gene expression, and contribute to the progression and development of cancer (35). One study demonstrated that miR-374 was upregulated in pituitary gonadotroph adenomas compared with the normal pituitary by reverse transcription-quantitative polymerase chain reaction (RT-qPCR) analysis (4). Palumbo et al (36) further confirmed that miR-145 was downregulated in pituitary tumors compared with normal tissues by performing RT-qPCR. Sachdeva and Mo (37) suggested that miR-145 could be regarded as a potential cancer biomarker and a novel target for the therapy of cancer. Furthermore, miR-374 and −145 have been verified to be able to regulate DEGs by previous experiments (26,27). Thus, the present study corroborates previous research by demonstrating that miR-374 and −145 may contribute to the development of pituitary gonadotroph adenomas via regulating the expression of target genes. To the best of our knowledge, no previous experimental studies concerning the association between miR-153 and −33 and pituitary adenomas have been performed; therefore, this topic should be further researched.

The upregulated and downregulated DEGs identified using the GSE26966 human dataset were not significantly enriched in the ‘neuroactive ligand-receptor interaction’ pathways or ‘cell cycle’ pathways. However, as the neuroactive ligand-receptor interaction pathway is a collection of receptors and ligands associated with intracellular and extracellular signaling pathways on the plasma membrane (38), a total of 5 of the top 10 upregulated GO-BP terms (including ‘ion transmembrane transport’, ‘positive regulation of synapse assembly’, ‘chloride transmembrane transport’, ‘sodium ion transport’ and ‘potassium ion transport’) were therefore associated with the neuroactive ligand-receptor interaction pathway.

The downregulated DEGs were predominately enriched in the GO-BP terms ‘negative regulation of cell proliferation’, ‘regulation of insulin-like growth factor receptor signaling pathway’, ‘aging’ and ‘regulation of cell growth’, among which ‘negative regulation of cell proliferation’, ‘aging’ and ‘regulation of cell growth’ were associated with ‘cell cycle’. Other enriched BP terms identified in this study include ‘cell differentiation’, ‘positive regulation of osteoblast differentiation’, ‘fat cell differentiation’, ‘positive regulation of cell proliferation’ and ‘cell cycle arrest’ were also associated with ‘cell cycle’. Therefore, the association of neuroactive ligand-receptor interactions and the cell cycle pathway were indirectly verified for the human array data of GSE26966, suggesting that the mechanisms associated with neuroactive ligand-receptor interactions and the cell cycle pathway may also serve important roles in the development of human pituitary gonadotroph adenomas.

Although the array data, GSE23207, which was used to perform the present analysis, had been analyzed by a previously published study (12), there were several differences between the present study and that of Zhang et al (13). First, as described in the introduction section, in addition to the screening of DEGs, PPI network analysis and GO and KEGG pathway enrichment analysis, miRNA-target gene regulatory network analysis was also performed. Secondly, verification was achieved using the human array data from GSE26966. Therefore, the present study provides new information to advance the understanding of the development of pituitary gonadotroph adenomas.

In conclusion, Cdk1, Ccna2, Ccnb1, ‘cell cycle’ and ‘neuroactive ligand-receptor interaction’ pathways may serve important roles in the development of pituitary gonadotroph adenomas, and Ccna2 and Ccnb1 may be involved in the development of pituitary gonadotroph adenomas via the cell cycle pathway. Furthermore, miR-374 and −145 may contribute to the development of pituitary gonadotroph adenomas via regulating the expression of target genes. These genes, pathways and miRNAs are potential targets of pituitary gonadotroph adenoma therapy. However, there were several limitations in this study, including the relatively small sample size and the lack of miRNA array data for human pituitary gonadotroph adenomas for the verification of the up- and downregulated miRNAs. Further studies concerning the association between miR-153 and −33, and pituitary gonadotroph adenomas/pituitary adenomas, are required.