Screening and identi�cation of key microenvironment-related genes in non-functioning pituitary adenoma

Abstract

construct a protein-protein interaction (PPI) network.Potential functions and pathways of intersection DEGs were then analyzed by GO (gene ontology) and KEGG (Kyoto Encyclopedia of genes and genomes).
Finally, the prognosis value of these genes was evaluated.

Results
The Kaplan-Meier curve indicated that the group of high immune scores were signi cantly related with poor recurrence-free survival.Then, we identi ed 497 intersection DEGs (478 over-expressed and 19 under-expressed genes) based on high-vs.low immune/stromal scores groups.Function enrichment analyses of 497 DEGs and hub genes from PPI network showed these genes mainly related with the immune/in ammatory response, T cell activation, and phosphatidylinositol 3 kinase-protein kinaseB (PI3K-Akt) signaling pathway.Among the IDEGs, 23 genes were further veri ed signi cantly expressed between invasive group and non-invasive group and 5 genes were highly associated with prognosis of NFPAs.

Conclusion
We screened out a series of critical genes associated with the tumor microenvironment in NFPAs.These genes may play a fundamental role in the development and prognosis of NFPA and may develop new therapeutic targets.

Background
Non-functioning pituitary adenomas (NFPAs) are benign pituitary neoplasms.They account for 14-54% of pituitary adenomas which are the second most common primary brain tumors [1][2][3][4][5].The delay in the diagnosis of NFPAs is related to the clinical manifestations of the lack of excessive hormone secretion.
Aside from accidental ndings, NFPAs is usually detected when the tumor is large enough to cause pressure effects on surrounding structures, with symptoms including headache, visual eld defects, decreased libido and hydrocephalus [6,7].Surgery is preferred when treatment is needed, whether or not combined with adjuvant radiotherapy [8].Although NFPAs are benign neoplasms, most NFPAs are invasive growth at the time of diagnosis, manifested as local in ltration, and have an increased risk of postoperative recurrence [9].However, the current understanding of the molecular mechanism that causes NFPAs and the mechanism of its invasion is still limited.Therefore, exploring the pathogenesis of NFPAs is essential for early detection of the disease and early prevention of the development of invasive NFPAs.
Recently, more and more studies have shown that the tumor microenvironment (TME) plays a vital role in tumor progression and treatment [10][11][12][13].The TME is the environment within which the tumor exists, consisting of blood vessels, lymphatic vessels, extracellular matrix (ECM), stromal cells, immune/in ammatory cells, RNAs, secreted proteins and small organelles [14].Among the components of TME, immune cells (lymphocytes and macrophages) and stromal cells (e.g.broblasts) are the main components of non-tumor cells and its could be impact on cancer progression and could be a new area of therapeutic target [15][16][17].For example, the poor prognosis of renal cell carcinoma is related with the in ltration of baseline CD8 T-cell [18,19].Research by Kashima, H. et al. shown that the accumulation of cancer-associated broblasts (CAFs) in esophageal squamous cell carcinoma could enhance the lymph node metastasis and promote tumor progression [20].Besides, CAFs could secret IL-6 to enhance epithelial-to-mesenchymal transition (EMT) in cancer cell resulting in chemoresistance among non-small cell lung cancer [21].According to the speci c gene expression signatures of immune and stromal cells, Yoshihara, K. et al. proposed a new algorithms called ESTIMATE (Estimation of Stromal and Immune cells in Malignant Tumor tissues using Expression data), which infers the proportion of stromal and immune cells in tumor samples by calculated immune and stromal scores [22].Some researchers have used this method in their studies to identify diagnostic and prognostic makers and immune related therapeutic targets in their studies [23][24][25].
In the current study, we used the ESTIMATE algorithm to calculate the stromal score and immune score in the tumor microenvironment of NFPAs to determine prognosis and target genes for treatment of NFPAs.
Based on the high and low immune/stromal score, we obtained a series of microenvironment related differential expression genes and the association between these genes and invasiveness and prognosis of NFPAs was also investigated.

Patients and samples
We collected 73 specimens from patients diagnosed with non-functioning pituitary adenoma and undergoing surgery at Beijing Tiantan Hospital from October 2007 to July 2014.There were 39 males and 34 females in these patients, the median age at diagnosis was 52 years (range 25 − 13 years).In addition, the median follow-up was 60 months (range 4-98 months).The invasive NFPAs was de ned by Knosp grading scale (grade 3 or 4) on preoperative MRI [26] or destruction of sellar oor, clivus or anterior skull base on preoperative CT.The clinical data of the adenomas are described in Table 1.Tumor specimens were stored in liquid nitrogen within half an hour after resection until RNA was extracted.This study was approved by the local Ethics Committee and informed consent was obtained from each subject.Total RNA extraction and RNA microarray Phenol-free mirVana ™ miRNA Isolation Kit (Cat # AM1561, Ambion, Austin, TX, USA) was used to extract and purify total RNA for the generation of uorescently labeled SBC human ceRNA array V1.0 cRNA targets (4 × 180 K).The labeled cRNA target is then hybridized with the slide.After hybridization, slides were scanned on an Agilent microarray scanner (Agilent Technologies, Santa Clara, California, USA).After extracting the data using feature extraction software 10.7 (Agilent Technologies, Santa Clara, California, USA), the Quantile algorithm was used to normalize the raw data by using the Limma software package of the R program.
Differentially Expressed Genes (degs) Based On Immune/stromal Scores By applying the ESTIMATE algorithm in R program (version 3.5.3) to calculated the immune and stromal scores of each sample [22].Differential expression genes were identi ed by using R program between high-immune/stromal score group and low-stromal/immune score group, and further study was considered with the |fold change| > 2 and the p value < 0.05.Venn diagrams was used to perform the consensus differentially expressed genes in both groups.

Construction of PPI Network
The construction, visualization and analysis of protein-protein interactions (PPI) are completed by the software of Cytoscape [27].The PPI datasets were obtained from BioGRID (http://thebiogrid.org/)(Release 3.4.140)and HPRD (http://www.hprd.org/)(Release 9), which reliability has been evaluated and widely used in several studies.The degree distribution was analyzed by Cytoscape plug-in and Network Analyzer.

Functional enrichment analysis
All of the intersection DEGs were enrichment analyzed by Gene ontology (GO) and Kyoto encyclopedia of Genes and genomes (KEGG) to predict the biological function, which was used the ClueGo[28] of Cytosccape plugin (version 3.2.3).GO terms analysis included biological process, cellular components, and molecular functions terms and KEGG pathway with adjusted p-value < 0.05 was considered signi cant.

Evaluation Of Intersection DEGs
Box-plots were used to perform the intersection gene expression levels between invasive and non-invasive groups.Besides, the statistical signi cance of the difference was calculated by the t-test.Kaplan-Meier plots and log-rank test p-values was used to clarify the correlation between each intersection genes and prognosis performance.Forest plot visualized the independent risk factors of intersection genes for recurrence-free survival by Cox regression analysis.

Results
The correlation of immune/stromal scores with clinical characters in NFPA We obtained the complete gene expression pro les by high throughput sequencing and identi ed 18827 protein-coding genes expressing in Non-functioning pituitary adenoma (NFPA) with expression value > 0 in the 73 samples.The stromal scores ranged from − 1648.8 to 1273.7, immune scores ranged from − 1429.1 to 2376.3, and ESTIMATE scores ranged from − 2866 to 3650 calculated by using the ESTIMATE algorithm.The distributions of immune and stromal scores did not vary with Knosp classi cation (Kendall p-value = 0.29/0.37),Hardy classi cation (Kendall p-value = 0.51/0.37),or pathological type (Kendall p-value = 0.12/0.13)(Fig. 1).
To investigate the relationship between prognosis and immune/stromal score, 73 NFPAs were classi ed into high-score groups and low-score groups and analyzed by Kaplan-Meier survival analysis and logrank test.The results reveled that the low immune scores signi cantly correlated with better recurrencefree survival (p = 0.01, Fig. 2a).Consistently, the low stromal score group was correlated with better RFS, although not statistically signi cant (p = 0.17, Fig. 2b).
Differentially expressed genes of immune/stromal scores in NFPA Compared on the gene expression pro les of high-and low-immune score group, we identi ed 879 immune-related DEGs (cut-off criteria as p < 0.05 and |fold change| > 2), including 818 upregulated genes and 61 downregulated genes (Fig. 2c).Similarly, based on the high-and low-stromal score group, we obtained 694 upregulated and 58 were downregulated stromal-related genes (Fig. 2d).In addition, Venn diagrams analysis revealed that there were 478 intersection upregulated genes and 19 intersection downregulated genes in both immune and stromal groups (Fig. 2e-f).These intersection DEGs were selected for subsequent analysis.

Functional Enrichment Analysis of Intersection DEGs
Based on the ClueGO gene annotation tool, we conducted GO analysis to further evaluate the biological functions of the 497 intersection DEGs (Fig. 3).Three subontologies of these DEGs were analyzed (Fig. 3a), as for biological processes (BP), these intersection DEGs mainly enriched in the T cell activation, in ammatory response, cell migration, cytokine production and so on; for cellular component (CC), these DEGs primarily clustered in the extracellular matrix, MHC protein complex and adherens junction; with regard of molecular function (MF), they are mainly associated with cytokine activity, cell adhesion molecule binding, chemokine receptor binding.Moreover, enrichment analysis of KEGG pathway suggested the they are mainly related to PI3K-Akt signaling pathway, cell adhesion molecules (CAMs), Cytokine-cytokine receptor interaction, and so on (Fig. 3b).

PPI Network Construction of Intersection DEGs
Protein-protein interaction (PPI) network was constructed based on 497 intersection genes aim to explore hub genes and develop a thorough picture of these genes at the systems level and it contains 6260 nodes and 12946 interactions (Fig. 4a).Figure 4a showed that these degrees of nodes followed a power-law distribution according to degree analysis of the network, which means the network was scale free.In addition, it revealed that the characteristic path length of the network was substantially larger than that of random networks by calculated average path length for the network (1000 times random network, P < 0.001, Fig. 4a), which implicated that the network had reduced global e ciency.Figure 4b shows the top 20 DEGs based on degree distribution in PPI network.The topological analysis suggested Interferoninducible 16 (IFI16, degree = 284) was the relatively key gene in the PPI network, which could interact with most genes as it had highest degree (Fig. 4b).Moreover, functional enrichment analysis of these top 20 DEGs also showed that they were mainly associated with the T cell activation, apoptosis signaling pathway, and angiogenesis (Fig. 4c).
The expression of stromal-immune score related IDEGs in NFPA Among the 497 intersection genes (478 upregulated genes and 19 downregulated genes), there are 23 genes which is signi cantly differentially expressed in invasive group and non-invasive group (p < 0.05).Although these genes are from the upregulated DEGs by analysis high-and low-stromal-immune score, they performed low expressed in the invasive group (|Foldchange|<1, P < 0.05).Figure 5 shows 8 genes (ABCA8,ANGPTL4 DLK1 EFCC1 PDLIM4 PREX2 TFPI2 GPC3) signi cantly differentially in invasive and non-invasive group by boxplot, these genes selected randomly from 23 genes (P < 0.05).

Survival analysis of intersection DEGs
The correlation of intersection DEGs with recurrence-free survival (RFS) was evaluated by constructed Kaplan-Meier survival curves using the 73 NFPAs.Among 497 intersection DEGs, 5 were predictive of the RFS (P < 0.05) by analyzed by log-rank test (Fig. 6a-e).Among 5 prognostic intersection DEGs, patients with upregulated ADGRG6, CD52, GPR183, and NNMT was associated with shorter RFS time than those with downregulated patients, while the upregulated TBX3 was associated with favorable outcomes.Moreover, forest plot performed the prognostic effects of these 5 DEGs by Cox regression analysis (Fig. 6f).

Discussion
Despite advances in endoscopic techniques, most NFPA patients are at increased risk of postoperative recurrence and lack of therapeutic response due to their invasive growth in surrounding normal tissues at the time of diagnosis [9,29].Therefore, it is indispensable to explore the pathogenesis of NFPAs and its invasiveness, which will help us understand the behavior of NFPAs invasion and recurrence and provide systematic treatment methods.TME (tumor microenvironment) is implicated in the initiation and development of various tumours [14,30], but research on the NFPAs is still scarce.In the present study, we devoted to screen out the TME-related genes that act on NFPAs RFS and invasiveness in our database.
First, the immune scores and stromal scores was obtained by using ESTIMATE algorithm.Although the relation between the immune/stromal score and pathological type, Knosp classi cation, and Hardy classi cation was not signi cant in our results, the low immune scores were signi cantly correlated with better recurrence-free survival, suggesting that the immune components in the TME may have deep effect in tumor progression.This is consisted with previous studies that tumor-associated macrophages (TAMs), T regulatory cells (Tregs) and mast cells are associated with poor prognosis in many type of cancer [14,[30][31][32].Although our analysis showed the relationship between stromal scores and prognosis was not signi cant, this does not mean that stromal-speci c genes have nothing to do with prognosis.
Second, we further analyzed 497 intersection DEGs (478 overexpressed and 19 underexpressed) which obtained from high vs. low immune/stromal scores groups, since the intersection genes are the most highly conserved.Among the 497 intersection genes analyzed by GO terms, we found that the most of them enrichened in the immune/in ammatory response (BP), cytokine activity (MF) and extracellular matrix (CC).KEGG pathway enrichment analysis shows that these IDEGs mainly clustered in PI3K-Akt signaling pathway, cell adhesion molecules (CAMs), Cytokine-cytokine receptor interaction.The PI3K-Akt signaling pathway is a signal transduction casacde involved in cell growth and metabolism [33].Long, R. et al. study revelas that COL6A6 could block PI3K-Akt-pathway to supress the growth and metastasis of pituitary adenoma [34].Moreover, we were able to construct a protein-protein interaction module based on these 497 intersection genes.Similarity, these hub genes are related with the B cell activation, T cell activation, angiogenesis and so on, which indicated that these genes were related to immune/in ammation and stromal response by function enrichment analysis.
Third, we analyzed 497 interaction genes and identi ed 23 genes that were signi cantly differentially expressed between invasion group and non-invasive group.Several identi ed genes in our results have previously been reported to play important roles in several types of cancer.Wang, Y. et al. study found that the expression of ABC transporter gene (ABCA8) was signi cantly down-regulated in prostate cancer as compared to noncancerous prostate tissues [35].This was consistent with our result that the expression level of ABCA8 was signi cantly decresed in invasive group vs. non-invasive group.Although several studies demonstrated that Angiopoietin-like 4 (ANGPTL4) was overexpressed in colorectal cancer, breast cancer, prostate cancer and gastric cancer [36][37][38][39], downregulated by promoter methylation in gastric cancer and mammary carcinoma [40,41].Similarly, Angiopoietin-like 4 (ANGPTL4) was downregulated in invasive group based on our data.Therefore, the gene of ANGPTL4 could play an important role in tumor microenvironmental, either as a tumor suppressor or oncogene.In addition, PDLIM4 (PDZ and LIM domain 4) was found signi cantly underexpressed in invasive group compared with non-invasive group.It also has been validated downregulation of PDLIM4 is associated with aggressive tumor features and poor prognosis in several cancer such as prostate cancer [42] and ovarian cancer [43].These genes may be potential biomarkers for early diagnosis and related to the prognosis of NFPAs.However, the function of these identi ed genes is still needed to investigate.
Finally, 5 TME-related genes were screen out which was signi cant related with the RFS.Among them, the upregulated four genes ADGRG6, CD52, GPR183 and NNMT are correlate with unfavorable outcomes of NFPAs.ADGRG6 is a member of the adhesion G protein-coupled receptor family and the depletion of ADGRG6 expression in urothelial bladder carcinomas cells compromised their abilities to recruit endothelial cell and induce tube formation [44].The upregulated CD52 is correlated with poor survival in patient with Lung adenocarcinoma by analysis [45].CD52 is a glycosylphosphatidylinositol (GPI)anchored protein expressed on the surface of normal T and B lymphocytes and also known as CAMPATH-1 antigen[46].In addition, NNMT(Nicotinamide N-methyltransferase) was found to be overexpressed in gastric carcinoma tissue compared with adjacent tissues and is related with the poor prognosis [47]

Conclusion
In summary, the tumor microenvironment-related genes might be correlated with the characteristic of tumor invasion and prognosis.Our study may provide a novel predict signature and therapy target for immune/stromal component in NFPAs. Figures

Figure 3 Function
Figure 3

Table 1
Summary of patient demographics and clinical characteristics.
Mult, Clinical non-functioning pituitary adenomas with multiple hormones positive immunostaining; GH, growth hormones.
[49]ich is a phase II metabolizing enzyme, mainly catalyzes the methylation of nicotinamide and other pyridines into pyridinium ions[48].Moreover, Wang, Y. et al. found that the overexpressed NNMT is correlated with poor survival and chemotherapy response in breast cancer patients who received chemotherapy[49].As explained above, a serious of genes we identi ed through bioinformatics methods are related to tumor microenvironment, and these genes are related to patients' prognosis.Besides, although the function and role of some of these genes have been reported in multiple types of tumor, further validation with biological experiments in our patient's samples are still needed.