Identification of biomarkers associated with the invasion of nonfunctional pituitary neuroendocrine tumors based on the immune microenvironment

Introduction The invasive behavior of nonfunctioning pituitary neuroendocrine tumors (NF-PitNEts) affects complete resection and indicates a poor prognosis. Cancer immunotherapy has been experimentally used for the treatment of many tumors, including pituitary tumors. The current study aimed to screen the key immune-related genes in NF-PitNEts with invasion. Methods We used two cohorts to explore novel biomarkers in NF-PitNEts. The immune infiltration-associated differentially expressed genes (DEGs) were obtained based on high/low immune scores, which were calculated through the ESTIMATE algorithm. The abundance of immune cells was predicted using the ImmuCellAI database. WGCNA was used to construct a coexpression network of immune cell-related genes. Random forest analysis was used to select the candidate genes associated with invasion. The expression of key genes was verified in external validation set using quantitative real-time polymerase chain reaction (qRT‒PCR). Results The immune and invasion related DEGs was obtained based on the first dataset of NF-PitNEts (n=112). The immune cell-associated modules in NF-PitNEts were calculate by WGCNA. Random forest analysis was performed on 81 common genes intersected by immune-related genes, invasion-related genes, and module genes. Then, 20 of these genes with the highest RF score were selected to construct the invasion and immune-associated classification model. We found that this model had high prediction accuracy for tumor invasion, which had the largest area under the receiver operating characteristic curve (AUC) value in the training dataset from the first dataset (n=78), the self-test dataset from the first dataset (n=34), and the independent test dataset (n=73) (AUC=0.732/0.653/0.619). Functional enrichment analysis revealed that 8 out of the 20 genes were enriched in multiple signaling pathways. Subsequently, the 8-gene (BMP6, CIB2, FABP5, HOMER2, MAML3, NIN, PRKG2 and SIDT2) classification model was constructed and showed good efficiency in the first dataset (AUC=0.671). In addition, the expression levels of these 8 genes were verified by qRT‒PCR. Conclusion We identified eight key genes associated with invasion and immunity in NF-PitNEts that may play a fundamental role in invasive progression and may provide novel potential immunotherapy targets for NF-PitNEts.


Introduction
Pituitary neuroendocrine tumors (PitNEts) account for approximately 10-20% of intracranial tumors and are the second most common neoplasms of the central nervous system (1,2). The prevalence of PitNEts ranges from 76-116 cases per 100,000 population, and the incidence is between 3.9 and 7.4 cases per 100,000 per year (3). These tumors are classified into functional and nonfunctioning pituitary tumor subtypes according to endocrine status (4). Nonfunctional pituitary neuroendocrine tumors (NF-PitNEts) account for 36%-54% of PitNEts and are usually detected based on signs and symptoms (headache, visual disturbance, and/or hypopituitarism) related to the effects of tumor mass because of the lack of excessive hormone secretion (2,(4)(5)(6). In this context, surgery is the treatment of choice because it can rapidly achieve decompression and symptomatic improvement (7)(8)(9). Because most macro-NF-PitNEts have the potential to invade the surrounding structures, such as the cavernous sinus or the sphenoid sinus, complete resection is often challenging and is achieved in up to 60-73% of patients (10,11). Moreover, invasive tumors have an increased recurrence rate due to tumor residues, which require additional surgery or radiation therapy and thus pose a further risk of complications (12)(13)(14). As a result, it is necessary to explore the pathogenesis of invasive NF-PitNEts to optimize the treatment of this tumor.
The tumor immune microenvironment (TIME) plays a crucial role in tumor development, progression, and immunotherapy (15,16). The TIME is composed of immune cells (lymphocytes and macrophages), immune-related pathways and cytokines secreted by tumor cells or immune cells (17). Pituitary tumor cells have been shown to recruit a variety of tumor-infiltrating immune cells, such as macrophages, T lymphocytes, B lymphocytes, FOXP3+ cells, neutrophils, and NK cells, into the tumor microenvironment (18)(19)(20). Moreover, the TIME has many effector functions and may promote the proliferation, migration and invasion of pituitary tumors (18,21,22). Therefore, it is essential to comprehensively analyze immunological genes affecting the abundance of immune cells in the invasive NF-PitNEts microenvironment.
In the current study, differentially expressed genes (DEGs) were identified at the tumor invasive and immune levels. Weighted correlation network analysis (WGCNA) was used to screen immune cell-related genes. The key invasive and immunological genes were further investigated by constructing a classification model and enrichment analysis. Our study screened out critical invasive-immune associated genes, which could provide new ideas for exploring immunological studies and some potential treatment strategies for NF-PitNEts patients.

Human tissue samples and clinical data
In this study, we used two cohorts to explore novel biomarkers in pituitary tumors. The first dataset included 112 patients, and another independent test dataset contained 73 patients. Tumor specimens were obtained from patients with NF-PitNEts (n=112) who underwent transsphenoidal surgical resection at Beijing Tiantan Hospital between June 2018 and July 2019. The diagnosis of NF-PitNEts is defined as the absence of clinical and biochemical evidence of overproduction of adenohypophysis hormone. The mean age of these 112 patients was 52 years (range, 21-75), and there were 63 males and 49 females. The demographics and clinicopathological features of the patients are summarized in Table 1. Tumor cavernous sinus (CS) invasion was defined as Knosp grade 3 and 4 or intraoperative evidence (23). The expression profiles and matching clinical information of independent test datasets (n=73) were described previously (24). In addition, 16 NF-PitNEts specimens (8 invasive and 8 noninvasive) were collected from the same hospital as an independent validation cohort (Table S3), and their expression levels were verified by quantitative real-time polymerase chain reaction (qRT-PCR). This study recruitment process and protocol were approved by the Medical Ethics Committee of Beijing Tiantan Hospital, and informed consent was obtained from all individual participants.
Total RNA extraction and RNA sequencing A total of 1-3 mg RNA per sample was extracted and purified from the collected specimens of NF-PitNEts. According to the instructions provided, sequencing libraries were constructed using the NEBNext ® Ultra ™ RNA Library Prep Kit for Illumina ® (#E7530L, NEB, USA). After the library was successfully generated (effective concentration >10 nM), the index-coded samples were clustered on the cBot cluster generation system using HiSeq PE Cluster Kit v4-cBot-HS (Illumina). The library was then sequenced on an Illumina platform, and 150-bp pairedend reads were generated. Raw data were filtered with FAST-QC, and the clean reads were then mapped to the human genome hg19 sequence (GRCh37) using HISAT2 (25). HTseq was used to generate gene counts, and the RPKM method was used to determine gene expression (26).

Differential expression analysis
The Estimation of Stromal and Immune cells in Malignant Tumor tissues using Expression data (ESTIMATE) algorithm (27) was used to obtain the immune levels of 112 pituitary tumor patients. Based on the median immune score, patients were divided into high and low groups. These 112 patients were also divided into invasion and noninvasion groups based on their clinical invasion information.
The "limma" R package was used to obtain the DEGs from the high vs. low immune score groups and invasion vs. noninvasion groups. Genes with an adjusted P value < 0.05 and |log2-fold change| ≥ 0.585 were filtered as DEGs.

Weighted gene coexpression network analysis
First, the ImmuCellAI database was applied to estimate the abundance of 24 immune cell types in the 112 patients with NF-PitNEts (28). The abundance of these immune cells in the invasion and noninvasion groups was then used for WGCNA.
Second, the "WGCNA" R package was applied to build a coexpression network of immune cell-related genes. Next, several gene modules were detected based on their similar expression patterns. Finally, the abundance of immune cells was associated with these gene expression modules, and genes in the modules were selected for further analysis.

Random forest analysis
The "randomForest" R package was used to select the candidate genes associated with invasion. We first divided the 112 patients into a training dataset (n=78) and a self-test (n=34) dataset based on the "sample" function of the R program. Then, random forest analysis (RF) was performed through the "randomForest" function in the training dataset, the error rate curve was drawn and changes in the error rate of different numbers of genes selected were observed. Finally, the "ggplot" R package was used to show the MeanDecreaseAccuracy and the best RF model of these genes.
Genes with the top 20 MeanDecreaseAccuracy were used to build the SVM model through the "e1071" R package, and the "pROC" package was used to perform the classification efficiency of the model in the training dataset, self-test dataset and independent test dataset.

Consensus clustering
The R package "ConsensusClusterPlus" was applied to explore the classification efficiency of 8 crucial genes. Subsequently, the "pROC" package was used to show the classification efficiency of these genes.

Enrichment analysis
Functional enrichment analysis was carried out by the "clusterProfiler" R package, and only functions with a p value<0.05 were selected.

qRT-PCR assay
The QIAGEN RNeasy Kit and High Capacity cDNA Reverse Transcription Kit were used to extract total RNA from the verified samples and conduct reverse transcription reactions. qRT-PCR was performed in a volume of 20 µl with Power SYBR ™ Green PCR Master Mix on a QuantStudio 3 and 5 System (Applied Biosystems). GAPDH was used as a housekeeping control. The sequences of the primers are shown in Table 2.

Statistical analysis
All statistical analyses were performed by R software (version 3.6.3). The bar plot between two different groups was drawn by the "ggplot" package, and a T test was used to compare the differences. A P value <0.05 was considered significant.

Results
Identification of immune infiltrationassociated and invasion-associated DEGs in pituitary tumor patients ESTIMATE analysis was performed in the 112 NF-PitNEts expression data, and the samples were divided into high-and low-immune score groups. Next, 3152 DEGs were obtained from the high-vs. low-immune score groups, including 662 upregulated and 2490 downregulated genes ( Figure 1A and Table S1). The enrichment analysis indicated that these DEGs were involved in Th1 and Th2 cell differentiation, the relaxin signaling pathway, and the FoxO signaling pathway ( Figure 1B). Then, we obtained 525 DEGs (352 upregulated and 173 downregulated) associated with invasion ( Figure 1C and Table S2). These genes are related to the FoxO signaling pathway and Th17 cell differentiation ( Figure 1D). This finding reveals that invasion-associated DEGs are involved in the immune-related functions of pituitary tumors.
Subsequently, we predicted the immune cell abundance of these patients. Cytotoxic T cells (Tc), Th2 cells, natural killer T cells (NKT), dendritic cells (DC), B cells, monocytes and neutrophils differed between the invasion and noninvasion groups ( Figure 1E).

Immune cell-associated modules in pituitary tumors
We used the "WGCNA" package to calculate the immune cellassociated modules in pituitary tumors. The soft power of the coexpression network was calculated through the "pickSoftThreshold" function and was set to 20 (Figure 2A). At this power, the R 2 of the scale-free topology model under the soft threshold was 0.92, which indicates that the network conformed to the scale-free feature (Figure 2A). Then, the network was constructed, and seven coexpression modules (green, turquoise, blue, yellow, black, brown, and red) were built ( Figure 2B). Subsequently, the correlation between these seven modules and immune cell abundance was calculated. The green module was significantly correlated with neutrophils; the turquoise module was significantly correlated with monocytes; the blue module was significantly correlated with cytotoxic cells; and the yellow, black, brown, and red modules were significantly associated with NKT cells ( Figure 2C). These results suggest that these genes are immune-related in the tumor microenvironment of pituitary tumors.

Construction of the invasion and immuneassociated classification model (IICM)
To select the candidate crucial genes in pituitary tumors, we first set the intersection of immune-related DEGs, invasion-related DEGs and module genes and found 81 common genes ( Figure 3A). Next, we used the 81 common genes to perform RF analysis. We started the RF analysis by generating 1000 decision trees, which showed a lower error rate ( Figure 3B). Then, the random forest results indicated that, when the number of candidate genes was 20, the classification efficiency error rate of the model was the lowest. We showed the top 50 accuracies of these genes and selected the top 20 for further analysis (Figures 3C, D).
Subsequently, SVM analysis was used to construct the IICM ( Figure 4A). We found that this model could distinguish the invasive and noninvasive patients in the training dataset (AUC=0.732, Figure 4A). It also exhibited good efficacy in the self-test dataset (AUC=0.653, Figure 4B) and independent test dataset (AUC=0.619, Figure 4C). The results showed that this model could be used to predict the invasion state of pituitary tumors.

Eight genes in IICM exhibit better efficacy in pituitary tumors
We then analyzed the function of these 20 key genes. The results suggested that 8 of the 20 genes were enriched in multiple signaling pathways, such as the KRAS signaling pathway and PPAR signaling pathway ( Figure 5A). This result suggests that the function of this model is mainly driven by these 8 genes (BMP6, CIB2, FABP5, HOMER2, MAML3, NIN, PRKG2 and SIDT2).
All 8 genes were differentially expressed in the invasion and noninvasion groups ( Figure 5B), and most were differentially expressed in the high-and low-immune score groups ( Figure 6). Subsequently, the consensus clustering analysis suggested that the 8 genes could well divide the patients into two groups (Figures 7A, B). We then constructed a classifier for these 8 genes, and the results showed that the 8-gene classification model showed good classification efficiency in pituitary tumors (AUC=0.671, Figure 7C). Then, qRT-PCR was performed to validate the expression level of these genes between invasive and noninvasive NF-PitNEts ( Figure 7D). Consistent with the sequencing data, the mRNA expression of BMP6, CIB2, HOMER2, MAML3, NIN, PRKG2, and SIDT2 was significantly upregulated in invasive NF-PitNEts, while FABP5 was significantly downregulated. Overall, we filtered 8 genes that could predict the invasion status of pituitary tumors and could be used as predictors for further treatment of the tumors.

Discussion
Although NF-PitNEts are benign neoplasms, they often invade surrounding structures and cannot be cured using standard therapies (29). Moreover, invasion is known as an important prognostic factor for recurrence (30). Recent studies have reported that immune cells infiltrate pituitary adenomas and may play an important role in tumor invasion and progression (31)(32)(33)(34). Therefore, understanding the mechanisms involved in immunity with invasive NF-PitNEts could lead to the discovery of new therapeutic targets in the future. First, 3152 DEGs (662 upregulated and 2490 downregulated) were identified based on the high-and low-immune score groups, which were obtained using the ESTIMATE algorithm. We obtained 525 DEGs between invasive and noninvasive NF-PitNEts. Then, the abundance of immune cells was predicted using the ImmuCellAI database, and Tc, Th2, NKT, DC, B cell, monocyte and neutrophil cells were found to be different between the invasion and noninvasion groups. Huang X et al. (35) found that patients with invasive NF-PitNEts had significantly lower CD3-CD56+ natural killer (NK) cells than patients with noninvasive NF-PitNEts in peripheral blood.
Subsequently, the abundance of these immune cells was used to construct a coexpression network of immune cell-related genes by WGCNA. As a bioinformatics method, WGCNA clustering results (coexpression gene modules) have high biological significance and reliability (36,37). Next, the 81 interacting genes were identified by immune-related DEGs, invasion-related DEGs and module genes significantly associated with immune cells. To narrow down the number of invasion-and immune-associated genes, random forest analysis was performed. The 20 genes with the highest RF score were selected to constrict the SVM model, and the classification efficiency was verified in the training, self-test, and independent test datasets. This result indicated that these genes play an important role in invasive behavior. Then, functional enrichment analysis was performed on these 20 genes, and 8 were found to be enriched in multiple signaling pathways, such as the KRAS signaling pathway, Notch signaling pathway, and PPAR signaling pathway. Liu et al. (38) found that upregulation of secreted phosphoprotein 1 affects tumor cell proliferation, migration, and invasion via the KRAS/MEK pathway in head and neck cancer. Feng et al. (39) reported that the Notch signaling pathway was associated with the invasion of growth hormone adenomas. Finally, these 8 invasion-and immune-related genes (BMP6, CIB2, FABP5, HOMER2, MAML3, NIN, PRKG2 and SIDT2) were validated by consensus clustering analysis and verified by qRT-PCR between invasive and noninvasive NF-PitNEts. Bone morphogenetic protein-6 (BMP-6) belongs to the TGF-b superfamily (40). BMP-6 was deemed to be associated with several tumor metastases, such as breast, prostate, rectal, and thyroid carcinomas (41-45). In addition, BMP-6 changes the morphology of macrophages and induces the expression of the cytokine tumor necrosis factor (TNF)-a (46). Calcium and integrin-binding protein 2 (CIB2) is a small EF-hand protein that can bind Mg 2+ and Ca 2+ ions and participates in basic cellular functions (47). Zhu et al. (48) found that CIB2 is correlated with cell proliferation, migration, and   another independent NF-PitNEts cohort, a large-scale multicenter cohort is needed for stronger validation. Second, our study only verified expression levels in NF-PitNEts tissues by qRT-PCR, but the underlying functions and mechanism of these genes still need to be explored in vivo and in vitro.
In summary, we conducted a comprehensive bioinformatic analysis and screened out immune-related genes that were significantly correlated with invasion in patients with NF-PitNEts. The current study may provide a novel potential immunotherapy target for invasive NF-PitNEts. A B FIGURE 5 Twenty crucial genes were enriched in multiple pathways. (A) Functional enrichment analysis of the 20 genes. The results showed that 8 genes played an important role in these pathways. (B) The expression levels of the 8 crucial genes between patients with and without invasion. *p < 0.05, **p < 0.01, ***p < 0.001. The expression levels of the 8 crucial genes in high Immune-score samples and low Immune-score samples. The results are expressed as the means ± SD (Student's t test. **p < 0.01, ***p < 0.001, ****p < 0.0001). ns, no significance. The results are expressed as the means ± SD (Student's t test. *p < 0.05; **p < 0.01).

Data availability statement
The data provided in the study are deposited in the Gene Expression Omnibus (https://www.ncbi.nlm.nih.gov/geo/), the accession number is: GSE169498, and another set of data is included in the supplementary material, further inquiries can be directed to the corresponding author/s.

Ethics statement
The studies involving human participants were reviewed and approved by the Medical Ethics Committee of Beijing Tiantan Hospital. The patients/participants provided their written informed consent to participate in this study.

Author contributions
WX, YZ worked on the conception and designed research, and eventually approved the manuscript. QF and YL was contributed to collected and analyzed clinical data of patients. JW and JG were dedicated to data analysis, interpretation, and drafting. All authors read and approved the final manuscript.

Funding
This study was supported by the National Natural Science Foundation of China (Grant codes: 82071558, 82071559, 82203174, and 82103028) and the Capital's Funds for Health Improvement and Research (Grant no.2020-4-1077).