Abstract Background Pituitary adenoma (PA) is a prevalent intracranial tumor. Metabolites differ between pituitary tumor and healthy tissues or among different tumor subtypes. However, the transcriptional changes in metabolic enzymes, which are usually seemed as targets for metabolic therapy, remain unidentified. Methods Using microarray data for 160 samples from the Gene Expression Omnibus database, across the four most common tumor subtypes, we present the integrated identification of differentially expressed genes (DEGs) between tumors and controls. Results Subtype-specific DEGs revealed 1081 prolactin tumor-specific DEGs, 437 nonfunctioning tumor-specific DEGs, and 217 common DEGs among the four subtypes. Functional enrichment showed that a lot of biological functions related to metabolism had changed. Twenty-one prolactin and twenty-three nonfunctioning tumor-specific metabolic-related DEGs are mainly involved in fatty acid and nucleotide metabolism, redox reaction, and gluconeogenesis. Eighteen metabolic-related DEGs enriched in the metabolism of xenobiotics by the cytochrome P450 pathway, sulfur metabolism, retinoid metabolism, and glucose homeostasis were abnormal in all subtypes of PA. Conclusion Based on a comprehensive bioinformatics analysis of the available PA-related transcriptomics data, we identified specific DEGs related to metabolism, and some of them might be new attractive therapeutic targets. Especially, PDK4 and PCK1 might be new attractive biomarkers and therapeutic targets. Keywords: pituitary adenoma, differentially expressed genes, metabolic dysregulation, microarray data, metabolism pathway Introduction Pituitary adenomas (PAs) are the second most common primary intracranial tumor with an overall prevalence of about 16.7% in the general population.[32]^1^,[33]^2 According to hormone secretion, PAs are grouped into clinical nonfunctioning (NF) adenomas and functional adenomas which include prolactin (PRL)-secreting adenomas, growth hormone (GH)-secreting adenomas, and adrenocorticotropic hormone (ACTH)-secreting adenomas.[34]^3 Although PAs are benign tumors, they significantly disturb hormone secretion and affect the quality of life and life span. Of the four most common PA subtypes, those associated with Cushing’s disease (ACTH-secreting), acromegaly (GH-secreting), and NF adenomas are generally treated with transsphenoidal surgery as first-line therapy. However, some PAs can invade the cavernous sinus or surrounding bone in the sellar region, and so it is sometimes difficult to achieve total surgical resection.[35]^3 Moreover, currently available pharmacotherapies are not satisfactory. Therefore, there is an urgent need for a new treatment strategy. Metabolic abnormality is considered as one of the hallmarks of tumor cells.[36]^4 Although been a benign tumor, PAs have also shown to undergo metabolic remodeling. Several studies have confirmed differences in metabolites between the pituitary tumor and healthy tissues or among different tumor subtypes by using mass spectrometry or nuclear magnetic resonance spectrometry.[37]^5^,[38]^6 Based on the metabolomics profiles of PAs, several differential metabolite-related enzymes were regarded as the specific therapeutic targets clinically. We previously reported that lactate dehydrogenase A (LDHA), a key glycolysis enzyme, was overexpressed in invasive PAs and was a promising therapeutic target.[39]^7 Along with the continuous development of transcriptomics techniques, a large amount of datasets were released in public databases, such as Gene Expression Omnibus database (GEO),[40]^8 ArrayExpress,[41]^9 and so on. These valuable resources bring an opportunity for exploring the molecular pathogenesis underlying the complex disease. For instance, the systematical investigation for various metabolic reprogramming of tumors was performed using the published transcriptomics data.[42]^10^,[43]^11 To date, there is almost no such transcriptomics research identifying metabolic-related genes and the mechanisms underlying PA development. In this study, several datasets on PAs from GEO databases were collected and systematically integrated into a large gene expression dataset. Subsequently, differentially expressed genes (DEGs) were screened, annotated, and analyzed statistically. As a result, screened DEGs were enriched in dozens of pathways associated with metabolism, and next, many metabolism-related genes were identified. These results showed that dysregulated metabolisms may be crucial roles in PA pathogenesis. This is the first transcriptome analysis of metabolism-related pituitary tumors, including four hormone subtypes and the normal pituitary. The strategy proposed, along with the metabolism-related DEGs identified in this study, provides a theoretical framework for understanding metabolic changes and developing metabolism-based therapeutics for PAs. Materials and Methods Identification of Eligible Gene Expression Datasets Published studies on the expression profile of human pituitary adenomas were screened by searching “pituitary”, “pituitary adenoma”, “pituitary tumor”, “prolactinoma”, “growth hormone”, and “Cushing disease” in the GEO database ([44]http://www.ncbi.nlm.nih.gov/geo/). Datasets containing the gene expression data of normal human pituitary or the four main subtypes of pituitary tumors were collected and for which raw data were not available were excluded. Data Preprocessing and Normalization For raw probe-level microarray datasets, those from Affymetrix platform were normalized using the robust multi-array average method,[45]^12 while those from Agilent and Illumina were quantile-normalized and then log[2]-transformed.[46]^13^,[47]^14 It is worth noting that the negative value was replaced with the smallest positive value in the gene expression matrix before log[2] transformation. By using array probe-gene annotation files from GEO, the microarray-specific probe identifiers (IDs) were parsed into gene Entrez IDs, and then the mean expression value was calculated as the expression value in repeated genes according to Entrez IDs. Data Integration and Batch Effect Removal To combine these datasets from diverse sources, we extracted the shared genes of all microarray platforms used in this study. Based on the shared genes, the multiple normalized datasets were merged into a large dataset by gene Entrez IDs. Then, we used a combat algorithm in sva package to remove the batch effect so that the comparability between these datasets or samples was enhanced.[48]^15 Identification of Differentially Expressed Genes DEGs between tumors and healthy controls for each subtype of PAs were identified using limma package.[49]^13 Specifically, linear model and empirical Bayes (eBayes) statistics were utilized to evaluate the differential expression for all genes, and finally, those with fold change (FC) >1.4 and P value <0.1 were considered as the DEGs. Subtype-Specific and Common Tumor Markers To pick up subtype-specific and common tumor markers, the intersection and union of subtype DEGs were computed and visualized using the Venn Diagram package.[50]^16 The genes appearing in only one subtype were considered as subtype-specific candidate markers, while those appearing in the four subtypes were considered as common candidate tumor markers. Gene Ontology and Pathway Enrichment Analysis To investigate the functional annotation of specific and common markers in each subtype, the enrichment analysis for gene ontology (GO) were performed using the DAVID web server (version 6.8) with a 0.05 cutoff criterion.[51]^17 Enriched pathway analysis was carried out using clusterProfiler package with a 0.05 cutoff criterion,[52]^18 based on Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway database. To determine the changed tendency of pathways in tumors, the Z-score was calculated in each term using the following formula:[53]^19^,[54]^20 graphic file with name M1.gif The N[up] N[down] separately represents the number of up- and down-regulated genes between tumors and normal controls, and the count is the number of DEGs belonging to this term. Specific Marker in the PPI Network Protein–protein interaction (PPI) networks were constructed using the STRINGdb package.[55]^21 The interaction pair in PPI networks with a combined score >0.7 was identified. The networks were visualized and analyzed in Cytoscape software (version 3.6).[56]^22 The Molecular Complex Detection (MCODE) plugin Cytoscape was applied to screen notable modules from the PPI network,[57]^23 and the following parameters were set in the advanced options: degree cutoff = 2, node score cutoff = 0.2, and K-core = 2. Significant modules were identified with scores and nodes > 2. The whole process steps are diagramed in [58]Figure 1. Figure 1. [59]Figure 1 [60]Open in a new tab The workflow of microarray data integrating and subsequently analyzing in this study. Nine datasets were analyzed in the study after normalization and batch effect removal, the DEGs were extracted by pairwise comparison (GH-NP, ACTH-NP, PRL-NP, NF-NP), and the visualization of clustered results showing with heatmap, respectively. The specific markers and common markers obtained from the DEGs overlap by the Venn diagram. Each specific and common marker was functionally characterized using gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) database. Then, protein–protein interaction network and network modules were constructed. Results Data Integration and Identification of DEGs Nine gene expression profiles were collected from the GEO database. A total of 160 samples included 52 ACTH-secreting PAs, 38 GH-secreting PAs, 40 NF PAs, 9 PRL-secreting PAs, and 21 normal pituitary (NP) samples ([61]Table 1). Table 1. Microarray Datasets for PA Collected for Analysis in This Study Subtypes of PA PA Samples Platform Reference of the Studied Datasets [62]GSE72490 ACTH 12 [63]GPL5175 Affymetrix [[64]47] [65]GSE93825 ACTH 40 [66]GPL18281 Illumina [[67]48] [68]GSE46311 GH 16 [69]GPL6244 Affymetrix [[70]49] [71]GSE51618 NF 10 [72]GPL6480 Agilent [[73]50] [74]GSE36314 PRL vs NP 4 vs 3 [75]GPL8300 Affymetrix [[76]51] [77]GSE26966 NF vs NP 14 vs 9 [78]GPL570 Affymetrix [[79]52] [80]GSE119063 PRL vs NP 5 vs 4 [81]GPL13607 Agilent [[82]53] [83]GSE62866 GH vs NF 9 vs 9 [84]GPL6480 Agilent [[85]54] [86]GSE63357 GH vs NF vs NP 13 vs 7 vs 5 [87]GPL570 Affymetrix [[88]55] [89]Open in a new tab After datasets were combined, 6818 shared genes were obtained. The large dataset was then performed global normalization and batch effect correction ([90]Supplemental Figure 1A), and the expression values of PAs and NPs were compared to identify the DEGs related to the subtype. To extract more DEGs from each subtypes, we set relax cut‐off criterion as p < 0.1 and |logFC| >1.4, provided by Xinxia Peng’ research,[91]^24 leading to the identification of 512 DEGs in GH-secreting PAs, 392 DEGs in ACTH-secreting PAs, 1581 DEGs in PRL-secreting PAs, and 1149 DEGs in NF PAs ([92]Supplementary Materials S1-4). Cluster analysis showed that tumor and normal samples were well separated ([93]Supplemental Figure 1B–E). Identification of Markers and KEGG Pathway Analysis of DEGs Twenty-two GH-specific DEGs, 1081 PRL-specific DEGs, 437 NF-specific DEGs, and 217 common DEGs were identified and visualized using the Venn diagram ([94]Figure 2A), except for ACTH-specific DEGs. Detailed information of specific and common DEGs are displayed in [95]Supplementary Materials S1-5. The top five up- and down-regulated specific DEGs are listed in [96]Table 2. Figure 2. [97]Figure 2 [98]Open in a new tab Overlap of DEGs and KEGG biomarker enrichment. (A) Traditional Venn diagram showing overlap between DEGs across four tumor types. (B–D) Top 20 pathways in PRL-specific, NF-specific, and common markers. * Full name: AGE-RAGE signaling pathway in diabetic complications. Table 2. Top Five Up-Regulated and Down-Regulated Specific DEGs Subtypes Up-Regulated logFC p Value Down-Regulated logFC p value GH TRAF4 1.736627 0.004794 ZKSCAN4 1.492069 0.000137 COL9A1 1.43918 0.008345 POP1 1.538518 0.001843 GLI3 1.539545 0.009605 IL12RB1 1.755727 0.001845 FOXD1 1.648576 0.028342 COPS7A 1.630502 0.003096 MXRA8 1.411838 0.036798 MFN1 1.402943 0.003139 PRL INHBB 2.231295 4.83E-07 DMWD 1.941849 0.000137 SPARCL1 3.132262 6.58E-07 PTPRN 2.405829 0.000160 TGFBR3 2.640177 1.22E-06 PPP1R11 1.576218 0.000304 HSPA2 2.851662 1.25E-06 GRIK2 1.643703 0.000387 VAMP8 2.507604 2.20E-06 UTF1 3.123438 0.000528 NF TIAM1 1.896465 2.02E-07 MAD2L1 1.499621 7.02E-10 TEAD3 2.207000 1.06E-06 MCM4 2.159558 1.71E-08 PIK3C3 1.647165 1.33E-06 PPIH 1.919924 4.98E-08 HPCAL1 2.498038 2.30E-06 ATP2B2 1.701433 7.86E-08 TCIRG1 1.733527 3.19E-06 AKAP6 1.533483 1.12E-07 [99]Open in a new tab Based on KEGG pathway mapping, PRL-specific DEGs were enriched mostly in osteoclast differentiation, toxoplasmosis, inflammatory bowel disease, and apoptosis ([100]Figure 2B). NF-specific DEGs were primarily enriched in glycosphingolipid biosynthesis-globo and isoglobo series, cAMP signaling pathway, proteasome, and serotonergic synapse ([101]Figure 2C). Common DEGs were enriched mostly in the PI3K-Akt signaling pathway, extracellular matrix–receptor interaction, TGF-β signaling pathway, and cytokine–cytokine receptor interaction ([102]Figure 2D). Metabolism KEGG Pathway and Gene Ontology Analysis To identify metabolic-related DEGs, subtype-specific DEGs and common DEGs enriching in metabolism GO terms and pathways were further analyzed. Based on GO analysis, PRL-specific DEGs were enriched in 28 metabolism-related GO terms ([103]Figure 3A), except for cellular response to retinoic acid (GO: 0071300, Z-score = −0.30), and all the other terms are up-regulated. However, NF-specific DEGs were enriched in 7 metabolism-related GO terms, and all terms are down-regulated ([104]Figure 3B). Detailed information about metabolism-related GO terms was provided in [105]Supplementary Materials S6. Figure 3. [106]Figure 3 [107]Open in a new tab GOCircle of metabolism-related GO terms in PRL-specific DEGs (A) and NF-specific DEGs (B). GOCircle was performed by using GOplot package.[108]^19 There are six metabolism pathways enriched in PRL subtype, including pentose phosphate pathway (Z-score = 0.35), fructose and mannose metabolism (Z-score = 3.16), fatty acid degradation (Z-score = 1.17), tryptophan metabolism (Z-score = 2.64), beta-alanine metabolism (Z-score = 2.12), and carbon metabolism (Z-score = 0.81) ([109]Figure 4A). Twelve metabolism pathways were enriched in NF subtype, including steroid biosynthesis (Z-score = −1.73), arginine biosynthesis (Z-score = −0.57), purine metabolism (Z-score = −1.66), tyrosine metabolism (Z-score = −1.00), cysteine and methionine metabolism (Z-score = −1.34), glycosaminoglycan biosynthesis-heparan sulfate/heparin (Z-score = −0.57), linoleic acid metabolism (Z-score = 0), arachidonic acid metabolism (Z-score = −0.18), glycosphingolipid biosynthesis-globo and isoglobo series (Z-score = −2.00), alpha-linolenic acid metabolism (Z-score = 0.57), pantothenate and CoA biosynthesis (Z-score = −0.57), and sulfur metabolism (Z-score = 0) ([110]Figure 4B). Figure 4. [111]Figure 4 [112]Open in a new tab Bubble plot of metabolism-related pathways in PRL-specific (A) and NF-specific (B) DEGs. Common DEGs showed enrichment in only two KEGG metabolism pathways: metabolism of xenobiotics by cytochrome P450 and sulfur metabolism. GO terms showed enrichment in glucose homeostasis and retinoid metabolic process ([113]Table 3). Table 3. Metabolism-Related KEGG Pathway and GO Terms in Common DEGs ID Description Genes Count p Value KEGG: 00980 Metabolism of xenobiotics by cytochrome P450 EPHX1, HSD11B1, GSTM5, ADH1A, ADH1C 5 5.92E-03 KEGG: 00920 Sulfur metabolism TST, PAPSS2 2 9.59E-03 GO:0001523 Retinoid metabolic process RBP4, RARRES2, RLBP1, AKR1C3, GPC3 5 6.99E-03 GO:0042593 Glucose homeostasis RBP4, PDK4, PCK1, NGFR, PAX6, POMC 6 8.60E-03 [114]Open in a new tab PPI Construction and Module Analysis To explore the role of each gene in the metabolic pathway and the mutual regulatory mechanism, PPI networks for specific DEGs in the PRL-enriched, NF-enriched and common-enriched in metabolism pathways were constructed, respectively. PPI module analysis showed that genes involved in metabolic processes and pathways were highly correlated ([115]Supplementary Materials S7). Based on PPI network Module analysis. The PRL-specific PPI network containing 21 nodes and 93 edges could be divided into 3 modules ([116]Figure 5A), modules 1, 2, 3 with scores of 12, 3, 3, respectively. In module 1 of the PRL-specific network, ACAA1 (acetyl-CoA acyltransferase 1), ACOX1 (acyl-CoA oxidase 1), ALDH2 (aldehyde dehydrogenase 2), CAT (catalase), ECHS1 (enoyl CoA hydratase, short chain, 1), and TALDO1 (transaldolase 1) were mainly responsible for fatty acid metabolism and oxidation–reduction process. Genes in module 2 were mainly responsible for the pentose phosphate pathway and nucleotide metabolism. Genes in module 3 were primarily responsible for glutamate decarboxylation and tryptophan metabolism. Figure 5. [117]Figure 5 [118]Open in a new tab Protein–protein interaction network of PRL and NF metabolic-specific markers. The dark purple nodes represent a higher tumor degree, and white nodes represent a lower tumor degree. Red circles represent up-regulated genes, and green circles represent down-regulated genes. The NF-specific PPI network containing 23 nodes and 33 edges was grouped into 3 modules ([119]Figure 5B), modules 1, 2, 3 with scores of 3.667, 3.5, 3, respectively. In module 1 of the NF-specific network, LDHC (lactate dehydrogenase C), LTA4H (leukotriene A4 hydrolase), MDH1 (malate dehydrogenase 1), GOT1 (glutamic-oxaloacetic transaminase 1), and BCAT1 (branched chain amino-acid transaminase 1) were most enriched in gluconeogenesis and transamination pathway; MAOA (monoamine oxidase A), TH (tyrosine hydroxylase), and COMT (catechol-O-methyltransferase) were the most enriched in dopamine metabolism; NPR2 (natriuretic peptide receptor 2), ENTPD1 (ectonucleoside triphosphate diphosphohydrolase 1), PDE4B (phosphodiesterase 4B), PDE6C (phosphodiesterase 6C), and PDE6D (phosphodiesterase 6D) were the most enriched in purine metabolism. These interactions between genes that play roles in different metabolic processes reveal the complex metabolic disturbances of the NF PAs. It is noticeable that most genes in the PRL subtype were up-regulated (n=19), while only three genes were down-regulated, containing FBP2 (fructose-1,6-bisphosphatase 2), MAOB (monoamine oxidase B), and PRPS1L1 (phosphoribosyl pyrophosphate synthetase 1-like 1). In contrast, most genes in the NF subtype were down-regulated (n=17), and six genes, including BCAT1, TH, NPR2, NOS3 (nitric oxide synthase 3), ALOX15 (arachidonate 15-lipoxygenase), and PLA2G2A (phospholipase A2, group IIA), were up-regulated. The PPI network on common DEGs contained 85 nodes and 179 edges. Among the identified metabolic target candidates for PA therapy, PDK4 (pyruvate dehydrogenase kinase 4) and PCK1 (phosphoenolpyruvate carboxykinase 1) were the top two metabolic genes with the highest connectivity to other genes ([120]Figure 6). Moreover, PDK4 connected to STAT5b (signal transducer and activator of transcription 5B), one of the hub genes in the network. Figure 6. [121]Figure 6 [122]Open in a new tab Protein–protein interaction network was constructed based on 217 common DEGs. The diameter of each protein was defined by its network degree, and blue color wi identified to be “metabolic-related genes” in this study. Green cycle means down-regulated. Discussion This is the first study to systematically analyze metabolism-related DEGs specific to the main PA subtypes. We analyzed gene expression in 9 microarray datasets based on the NP and four major subtypes of human PAs. A total of 139 tumor samples and 21 normal samples were included. Except for the ACTH subtype, 22 GH-specific DEGs, 1081 PRL-specific DEGs, 437 NF-specific DEGs, and 217 common DEGs were identified. We utilized bioinformatics methods to explore the biological function of these DEGs, including GO and KEGG pathway enrichment, PPI network construction, and metabolic-related DEG annotation. KEGG pathway enrichment of common DEGs revealed many previously identified therapy targets. Signaling pathways such as PI3K-Akt, TGF-β, JAK-STAT, MAPK, and Hippo are well known to be involved in the pathophysiology and treatment of pituitary adenomas.[123]^25^–[124]^30 The PI3K-Akt and Hippo signal pathways were also revealed by ubiquitin proteomics analysis of human pituitary adenoma and normal pituitary tissues.[125]^31 In addition to these signal transduction pathways, the neuroactive ligand–receptor interaction pathway was shown to be related to PA tumorigenesis in studies using microarray data.[126]^32^,[127]^33 The extracellular matrix (ECM) delivers signals through membrane receptors and is involved in PA cell proliferation and hormone production.[128]^34^,[129]^35 These results indicate the reliability of the present data; however, studies of these related potential therapeutic targets are mostly in the laboratory stage. Although recent efforts have renewed awareness of cancer as a metabolic disorder, and the causal relationship between an oncogene and metabolic reprograming is still controversial.[130]^36 However, there is a general consensus regarding the idea that the change of metabolic enzyme expression ultimately alters the metabolic flow. Recently, transcriptomics data, in conjunction with the current biochemical understanding, have been exploited to construct genome-scale metabolic workflows.[131]^11^,[132]^37 In the functional enrichment of common DEGs, we found that the metabolism of xenobiotics by the cytochrome P450 pathway ([133]Figure 2D) is correlated with two Phase III clinical trial drugs, levoketoconazole ([134]NCT03277690) and osilodrostat ([135]NCT02697734). Another enriched KEGG metabolism pathway was sulfur metabolism, which was also enriched by PAPSS2 (3ʹ-phosphoadenosine 5ʹ-phosphosulfate synthase 2) and TAT (thiosulfate sulfur transferase) in the common DEGs. Sulfation is a common modification of exogenous (drugs and xenobiotics) and endogenous (carbohydrates, lipids, and proteins) compounds.[136]^38 Although few studies have examined the two genes in PA, they are suitable as therapeutic targets in future studies. Metabolism-related GO terms of common DEGs mentioned retinoid metabolism and glucose homeostasis. Our previous study about LDHA indicated that glucose was the most critical energy source of PA.[137]^7 Retinoid acid and isotretinoin had been proved to be effective and safe alternative therapeutic options for patients with Cushing disease.[138]^39^,[139]^40 Altogether, these results indicate the possibility of metabolism-related therapy for PA and also the reliability of the present data. Additionally, in the PPI network of common DEGs, we identified PDK4 and PCK1 as two key molecules involved in PA glucose metabolism and other related signaling pathways in PA development. Inhibition of PDK4 reportedly increases colon cancer cellular apoptosis and reduces cellular migration and invasion.[140]^41 Proteins connected to PDK4 in the PPI network such as PON3, CEBPB, PCK1, and STAT5B have also been reported to be involved in cancer development.[141]^42^–[142]^45 Moreover, STAT5B was the downstream molecule of PRL, GH1, and GH2, and its mutation had been proved to be involved in growth hormone insensitivity.[143]^46 Overexpression of STAT3, another gene of STAT family, had been proved to induce growth hormone hypersecretion in pituitary somatotroph adenomas.[144]^27 However, the roles of STAT5B, PDK4, and PCK1 in PA are still unknown. Aside from these observations, the current status is that till date, no drugs have been approved for treating non-functional PA; however, the present study provides some potential metabolic therapeutic targets. In conclusion, based on a comprehensive analysis of the available PA-related transcriptomics data, we identified specific DEGs in different hormone-secreting PA subtypes and found that several biological functions related to metabolism had changed. PRL and NF tumor-specific DEGs were mainly involved in fatty acid and nucleotide metabolism, redox reaction, and gluconeogenesis. Metabolism of xenobiotics by the cytochrome P450 pathway, sulfur metabolism, retinoid metabolism, and glucose homeostasis were abnormal in all subtypes of PA. PDK4 and PCK1 might be considered as targets for the development of anticancer strategies and therapies. However, the raw transcriptomics data came from different sources, the clinical or ethnicity variables might affect the results of the particular research, and more solid scientific work was needed to verify these findings. Acknowledgments