Identifying critical protein-coding genes and long non-coding RNAs in non-functioning pituitary adenoma recurrence.

Non-functioning pituitary adenoma (NFPA) is a very common type of intracranial tumor. Monitoring and predicting the postoperative recurrence of NFPAs is difficult, as these adenomas do not present with serum hormone hypersecretion. Long non-coding RNAs (lncRNAs) and protein-coding genes (PCGs) play critical roles in the development and progression of numerous tumors. However, the complex network of RNA interactions related to the mechanisms underlying the postoperative recurrence of NFPA is still unclear. In the present study, 73 patients with NFPA were investigated using high-throughput sequencing and follow-up investigations. In total, 6 of these patients with recurrence within 1 year after surgery were selected as the fast recurrence group, and 6 patients with recurrence 5 years after surgery were selected as the slow recurrence group. By performing differential expression analysis of the fast recurrence and slow recurrence groups, a set of differentially expressed PCGs and lncRNAs were obtained (t-test, P<0.05). Next, protein-protein interaction coregulatory networks and lncRNA-mRNA coexpression networks were identified. In addition, the hub lncRNA-mRNA modules related to NFPA recurrence were further screened and transcriptome expression markers for NFPA regression were identified (log-rank test, P<0.05). Finally, the ability of the hub and module genes to predict recurrence and progression-free survival in patients with NFPA was evaluated. To confirm the credibility of the bioinformatic analyses, nucleolar protein 6 and LL21NC02-21A1.1 were randomly selected from among the genes with prognostic significance for validation by reverse transcription-quantitative PCR in another set of NFPA samples (n=9). These results may be helpful for evaluating the slow and rapid recurrence of NFPA after surgery and exploring the mechanisms underlying NFPA recurrence. Future effective biomarkers and therapeutic targets may also be revealed. Copyright: © Guo et al.

Introduction

Pituitary adenomas are benign neuroendocrine tumors (1,2) that originate from adenohypophyseal cells, and they account for 10–20% of intracranial neoplasms, in epidemiological data from the United States between 2005 and 2009 (3–5). Pituitary adenomas can be divided into functional and non-functional adenomas according to clinical and biochemical characteristics, like GH-secreting adenoma characterized acromegaly caused by growth hormone abnormal rise (6). Non-functional pituitary adenomas (NFPAs) are the most common and account for 43% of pituitary adenomas, in epidemiological data from Iceland between 1955 and 2012 (7,8). NFPA is often characterized by a lack of symptoms associated with excessive hormone production, like acromegaly and Cushing's disease (6). Due to the mass effect on surrounding structures, NFPA may cause headaches, visual defects, and/or hypopituitarism (7,9). Surgical resection is the primary treatment for NFPA, although patients are often left with tumor residue, as the tumor can invade the cavernous sinus or area surrounding the internal carotid artery (10,11). In total, 12–58% of patients with NFPA with macroadenoma may experience regrowth within 5 years (12–15). Radiotherapy is often recommended for patients with tumor residue, although its long-term complications, such as visual defects and hypopituitarism, are still of concern (16,17). Therefore, surgery is still the best option for patients with tumor recurrence. Serum hormone monitoring is an approach used to detect functional pituitary adenoma (18); however, early intervention is difficult to achieve due to the absence of an effective evaluation approach for NFPA. Therefore, research on the molecular mechanisms underlying tumor recurrence and effective prognosis prediction methods is important.

Studies have shown that protein-coding genes (PCGs) are involved in the activation of pathways or key proteins and play vital roles in the biological processes of pituitary adenomas. For example, Uraki et al (19) showed that reducing the expression of MSH6 and MSH2 can directly promote the growth of pituitary tumors through the ATR-Chk1 pathway. Long et al (20) suggested that collagen α VI chain interacts with P4HA3 to inhibit pituitary adenoma cell proliferation and invasion by inhibiting the PI3K-Akt pathway. The low expression of TGF-β RII may be related to the development and invasion of NFPAs (21), and Zhu et al (22) confirmed that the expression of TGF-β1 and Wnt inhibitory factor 1 (WIF1) in recurrent tumors is higher compared with that in primary tumors, suggesting that these PCGs may be related to cell proliferation and recurrence. Compared with non-invasive NFPAs, the expression levels of WIF1 and secreted frizzled-related protein 4 are reduced in invasive NFPAs; thus, WIF1 may be a potential biomarker for the aggressiveness of NFPAs (23).

Long non-coding RNAs (lncRNAs) are a type of RNA molecule with a transcript >200 nucleotides in length, and they play an important role in regulating gene expression through epigenetic or posttranscriptional mechanisms; however, they do not encode proteins (24–26). The differential expression and dysregulation of lncRNAs is considered to be involved in carcinogenesis and cancer progression, recurrence and metastasis (24). However, the role of lncRNAs in NFPA recurrence and the regulation of cellular processes remains unknown. A study have shown that LINC00858 plays a tumor-promoting role in colon cancer by upregulating hepatocyte nuclear factor 4-α and downregulating WNK2 (27). Xu et al (28) showed that the overexpression of lncRNA PAXIP1-AS1 can upregulate KIF14, thereby enhancing human umbilical vein endothelial cell migration, invasion and angiogenesis in gliomas. Moreover, several studies found that identifying novel lncRNA-mRNA networks using microarray analyses could contribute to exploring the potential molecular mechanisms and prognosis of tumors (29–31). The aforementioned studies indicate that the dysregulation of lncRNAs and lncRNA-mRNA interactions may affect the prognosis of NFPAs.

The present study aims to screen out the critical PCGs and lncRNAs, which play an essential role in NFPAs recurrence. We obtained differentially expressed PCGs and lncRNAs by performing differential expression analyses of recurrence within 1 year after surgery (fast recurrence group) and after 5 years (slow recurrence group). Protein-protein interaction (PPI) networks and coregulatory networks between lncRNAs and mRNAs were also identified. The hub lncRNA-mRNA modules related to NFPA recurrence were further screen and the enrichment of the differentially expressed genes (DEGs) was assessed in different pathways by Gene Set Enrichment Analysis (GSEA). In addition, the ability of the hub and module genes [nucleolar protein 6 (NOL6), cyclin dependent kinase 15 (CDK15), Moloney leukemia virus 10 (MOV10), SAMM50 sorting and assembly machinery component (SAMM50), collagen type XXIV α 1 chain (COL24A1), epoxide hydrolase 1 (EPHX1) and decapping mRNA 1A (DCP1A)] to predict recurrence and progression-free survival (PFS) time was evaluated in patients with NFPA. These results may help us explore the mechanisms underlying NFPA recurrence, and may also provide future effective biomarkers and therapeutic targets.

Materials and methods

Patients and samples

Patients diagnosed for NFPA (n=73) who underwent surgical resection at Beijing Tiantian Hospital (Beijing, China) from October 2007 to July 2014 were included in the study. The study inclusion criteria were: i) Patients older than 18 years; ii) MRI/CT showed a sellar region lesion; iii) pathological diagnosis was pituitary adenoma with no hormonal excess and iv) sufficient pre- and postoperative clinical and radiologic data. The exclusion criteria were: i) Patients with functioning adenoma, including tumor which secreted ACTH, prolactin, growth hormone and/or TSH and lead to the corresponding clinical syndromes of hormone excess and ii) history of pituitary surgery or radiotherapy. The patients included 34 men and 39 women, with a median age of 52 years (age range, 25–73 years). Out-patient clinic follow-up was conducted, the minimum follow-up time was 4 months, and the median follow-up time was 60 months (range, 4–98 months). The clinicopathological characteristics of all patients are shown in Table I. All tumor samples were immediately placed into a sample tube, frozen in liquid nitrogen and stored. Among them, 6 cases of recurrence within 1 year were randomly selected as the fast recurrence group and 6 cases of recurrence after 5 years were selected as the slow recurrence group. The postoperative recurrence of NFPA refers to the increase in the maximum tumor diameter by >2 mm from the day of surgery to the end of follow-up as measured from any direction on magnetic resonance imaging. The histological subtype of tumors was defined according to the World Health Organization 2017 classification of endocrine tumors (32). The Medical Ethics Committee of Beijing Tiantan Hospital (Beijing, China) approved the study.

Table I.

Summary of non-functioning pituitary adenoma clinical characteristics.

Characteristic	Value, n
Sex
Female	39
Male	34
Age, years
≤52	41
>52	32
Invasion
Yes	47
No	26
Histological type
GA	41
SA	29
NC	3
Tumor size classification
Macroadenoma	53
Giant adenoma	20
Headache
Yes	35
No	38
Vision and visual field disorders
Yes	53
No	20
Recurrence
Yes	27
No	46

GA, gonadotroph adenoma; SA, silent adenoma; NCA, null cell adenoma.

Total RNA extraction and RNA microarray

The Phenol-free mirVana™ miRNA Isolation kit (cat. no. AM1561; Ambion; Thermo Fisher Scientific, Inc.) was used to extract and purify total RNA to generate fluorescently labeled cRNA targets (4×180 K), following the manufacturer's protocol. The labeled cRNA targets were then hybridized to a glass slide. After hybridization, the slides were scanned using an Agilent microarray scanner (Agilent Technologies, Inc.). After extracting data using Feature Extraction software 10.7 (Agilent Technologies), the Quantum algorithm was used to normalize the raw data using the limma software package of the R program (http://www.R-project.org). The analysis was performed by version 3.6.1. (http://www.rstudio.com/).

Identification of differentially expressed lncRNAs and mRNAs

A differential gene expression analysis was performed within 1 year after the initial postoperative NFPA (n=6) and 5 years later (n=6), and a significance analysis of microarrays (SAM) was performed to identify the differentially expressed PCGs and lncRNAs (DEGLs) between the two groups (33). The Biobase, multtest and siggenes packages were downloaded from Bioconductor (http://www.bioconductor.org/). Subsequently, the available data were analyzed using R (www.r-project.org), and DEGLs with fold-changes of >2 and <-2, and P<0.05, were selected for further research.

Construction of a PPI network and lncRNA-mRNA coexpression network

Cytoscape software (version 3.2.3) was used to construct, visualize and analyze the PPI network (34). The latest version of the validated human PPI dataset was downloaded from both the Human Protein Reference Database (HPRD) (www.hprd.org/; release 9) and BioGRID (www.thebiogrid.org/; release 3.4.140) (35,36). These two datasets contain 18,595 unique proteins and 174,552 interactions, and were used as parent PPIs in the present study; their reliability has been effectively verified, and they have been used extensively in disease research involving human PPI networks. The non-redundant interactions in Homo sapiens from these two data sets were manually integrated (37).

First, a PPI subnetwork was generated by mapping all the DEGs and extracting them from the PPI network. To improve reliability, network reconstruction was limited to the first interacting protein neighbors of these DEGs. Second, the DEG-adjacent protein axis was detected, and a DEG-central PPI network was constructed. Third, after mapping all DEGs to the PPI network to detect internal interactions between the DEGs, Cytoscape was used to select nodes with all edges to create a subnetwork. The single-node and self-interactions of proteins in these subnets were deleted. Pearson's correlation test was used to calculate the coexpression relationships between lncRNAs and PCGs, and the coexpression relationships with a P<0.05 and a Pearson coefficient absolute value of >0.9 were selected. Finally, an lncRNA-mRNA network related to NFPA composed of differential genes was obtained. The PPI-lncRNA network was visualized using Cytoscape, and Molecular Complex Detection (MCODE) (38) was used to identify important modules in the PPI network. The key modules and hub genes were further analyzed and visualized using the MCODE plugin in Cytoscape. The screening criteria for module genes were as follows: Degree cut-off, 2; node score cut-off, 0.2; k score, 2; and maximum depth, 100.

Functional enrichment analysis

The ClueGO (39) plugin of Cytoscape was used to perform Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses of the DEGs and the biological functions of the lncRNAs. Functional annotations with P<0.05 were considered significant. In addition, GSEA was used to identify the relevant pathways of the selected genes. GSEA was performed using GSEA (4.1.0) software (40). The gene set used in the study was downloaded from the Molecular Signatures Database (http://software.broadinstitute.org/gsea/msigdb/index.jsp, MSigDB version 4.0, released June 7, 2013). The Molecular Signatures Database contains various types of gene sets. The online pathway database includes 1,320 canonical pathways derived from pathway databases such as BioCarta, KEGG, Pathway Interaction Database and Reactome (40).

Validation and efficacy evaluation of the hub genes by survival analysis

Among the hub genes, genes of interest that have not been studied with regard to NFPA were further validated in the two groups. The PFS analysis of the hub genes and module genes was performed using Kaplan-Meier curves in the R program. P<0.05 was considered to indicate a statistically significant difference.

Validation of gene expression using reverse transcription-quantitative (RT-q)PCR

RT-qPCR was performed using another set of NFPA samples to verify the credibility of the bioinformatics analyses. According to the inclusion and exclusion criteria, five cases of NFPAs with recurrence within 1 year and four cases of NFPAs with recurrence after 5 years were randomly selected from the patients who underwent surgical resection in Beijing Tiantan Hospital between August 2009 and November 2014 as the validation set. The total RNA of validated samples was extracted and purified as aforementioned. Reverse transcription into cDNA was performed using a High Capacity cDNA Reverse Transcription kit (cat. no. 0049472; Thermo Fisher Scientific, Inc.) following the manufacturer's instructions. Next, a Power SYBR™ Green PCR Master Mix (cat. no. 4367659; Thermo Fisher Scientific, Inc.) was used for qPCR with a total reaction volume of 20 µl. Amplification was performed as follows: 95°C for 10 min, followed by 40 cycles at 95°C for 15 sec and 60°C for 60 sec, and a final extension at 72°C for 5 min. GAPDH was used as an internal control gene. All primers were synthesized by Sangon Biotech Co., Ltd. The level of mRNAs was determined using QuantStudio 3 and 5 systems (Applied Biosystems; Thermo Fisher Scientific, Inc.). For relative quantitation, expression levels were calculated expression levels were calculated using the 2−ΔΔCq method (41). The sequences of the primers are as follows: GAPDH forward, 5′-ACCCACTCCTCCACCTTTGA-3′ and reverse, 5′-CCACCCTGTTGCTGTAGCCA-3′; NOL6 forward, 5′-ATTCGGGAAGCTGTGGTCTG-3′ and reverse, 5′-ATGTCAGCATGGAGTGCCAA-3′; and LL21NC02-21A1.1 forward, 5′-CTGCCCGATCTCACCTCTTC-3′ and reverse, 5′-TCAGGGAAGGACTCCAGGTT-3′.

Statistical analysis

Data are presented as the mean ± standard deviation, unless otherwise shown. All experiments were performed in triplicate. The differentially expressed PCGs and lncRNAs were identified using an unpaired t-test. The difference in mRNA and lncRNA expression level between the two groups was assessed by the Mann-Whitney U test with GraphPad Prism 7 software (GraphPad Software). Kaplan-Meier analysis of PFS was conducted using the log-rank or the Renyi test (when there was survival curve crossover between the observed groups) using the R packages ‘survival’ and ‘survMisc’ in R (3.5.1). P<0.05 was considered to indicate a statistically significant difference.

Results

Identification of DEGLs between the fast recurrence and slow recurrence groups

Through microarray sequencing of 73 NFPA samples, 18,827 PCGs and 19,740 lncRNAs with expression values >0 were identified. The differences in PCGs and lncRNAs between 6 cases of NFPA recurrence within 1 year after surgery and 6 cases of NFPA recurrence 5 years after surgery were analyzed using the SAM test. By selecting the threshold |fold-change| >2 or adjusted P<0.05, a total of 299 differentially expressed PCGs (228 upregulated and 71 downregulated PCGs) and 214 differentially expressed lncRNAs (120 upregulated and 94 downregulated lncRNAs) were identified (Fig. 1A and B). The expression heat map further validated the results, and Fig. 1C and D shows the differential PCGs and lncRNAs with different expression trends for recurrence.

Figure 1.

Identification of differentially expressed PCGs and lncRNAs related to recurrence in >5 years or <1 year. Volcano plot of (A) mRNAs and (B) lncRNAs. Red dots indicate genes without significantly different expression. Blue dots indicate significantly regulated genes. Expression of differential (C) PCGs and (D) lncRNAs in the slow recurrence group (recurrence after >5 years) or rapid recurrence group (recurrence in ≤1 year). lnc, long non-coding.

Pathway enrichment of the DEGs by GSEA classifies the fast recurrence and slow recurrence groups

GSEA demonstrated that 30 different pathways related to 299 differentially expressed PCGs were downregulated or upregulated according to the recurrence rate. Several enriched terms are presented in Fig. 2 and Table SI, such as ‘regulation of cell death’, ‘regulation of cell adhesion’, ‘positive regulation of biosynthetic process’, ‘positive regulation of gene expression’, ‘small molecule metabolic process’ and ‘response to extracellular stimulus’. The results indicated that changes in these pathways lead to the recurrence and progression of NFPAs.

Figure 2.

Pathway enrichment of the DEGs. (A-F) Gene Set Enrichment Analysis of DEGs between the rapid recurrence and slow recurrence groups. The enriched terms including (A) positive regulation of biosynthetic process, (B) regulation of cell death, (C) positive regulation of gene expression, (D) small molecule metabolic process, (E) response to extracellular stimulus and (F) regulation of cell adhesion. DEGs, differentially expressed genes; GO, Gene Ontology.

Dysregulated lncRNA-mRNA interaction network establishment and module analysis

Based on Pearson's correlation test, a differential coexpression network of lncRNAs and mRNAs was constructed, selecting genes with P<0.05 and a Pearson's coefficient absolute value of >0.9 (lncRNA/mRNA quantity, 78/104; Fig. 3A). This network was transferred to the differential PCG PPI parent network. Subsequently, the lncRNA-mRNA interaction network was obtained by combining these two networks (Table SII). The lncRNA-mRNA network for the DEGLs contained a total of 4,490 nodes and 6,933 interactions. Fig. 3B shows that the degrees of the genes followed a power-law distribution, further illustrating that the network was similar to most biological networks, and the network is scale-free. The average path length of the network was also calculated, and it indicated that the characteristic path length of the network was much longer compared with the path length of the random network (1,000 times longer; P<0.001; Fig. 3C), which implied that the network had reduced global efficiency.

Figure 3.

Topological features of the non-functioning pituitary adenoma progression-related lncRNA-mRNA network. (A) Relation of dysregulated lncRNAs coexpressed with dysregulated mRNAs visualized by a Circos plot. (B) Degree distributions of the network; all degrees followed a power-law distribution. (C) Average path length distributions of the real network and 1,000 random networks. lnc, long non-coding.

Numerous studies have shown that PCGs and lncRNAs usually function by participating in functional modules (42–44). Through cluster analysis of the PPI network using the MCODE plugin of Cytoscape, eight important modules according to the degree of importance were obtained. Module 1 contained 46 nodes and 54 edges (Fig. 4A), module 2 contained 71 nodes and 79 edges (Fig. 4B), module 3 contained 59 nodes and 61 edges (Fig. 4C), module 4 contained 59 nodes and 65 edges (Fig. S1A), module 5 contained 55 nodes and 55 edges (Fig. S1B), module 6 contained 44 nodes and 44 edges (Fig. S2A), module 7 contained 83 nodes and 88 edges (Fig. S2B), and module 8 contained 214 nodes and 218 edges (Fig. S2C).

Figure 4.

Module analysis of the non-functioning pituitary adenoma-related lncRNA-mRNA network using Molecular Complex Detection. (A-C) lncRNA-mRNA interactions in three modules and pathway enrichment of protein-coding genes in each corresponding module. (A) Module 1, (B) Module 2 and (C) Module 3. Square nodes represent lncRNAs, and elliptical nodes represent mRNAs. Red indicates upregulation and green indicates downregulation. lnc, long non-coding.

ClueGO was used to perform an enrichment analysis of the genes in these modules. As shown in Figs. 4, S1 and S2, the GO analysis indicated that the genes in modules 1–8 were mainly concentrated in the categories ‘T cell migration’, ‘T cell chemotaxis’ and ‘T cell activation’, ‘regulation of cell-cell adhesion’ (Fig. 4A-C) and ‘regulation of cytokine production involved in the immune response’ (Fig. S1A). In addition, the KEGG analysis showed that enrichment of these module genes mainly occurred in the categories ‘cell cycle’, ‘adherens junction’, ‘TNF signaling pathway’ (Fig. 4A-C), ‘VEGF signaling pathway’ and ‘TGF-β signaling pathway’ (Fig. S1A and B). The expression of 21 hub genes and 10 module lncRNAs from the aforementioned eight modules [CCR1, CCL3, CCL4, PACSIN1, CGNL1, TRIM69, STAB2, ATP8A1, CD48, NOL6, ZFP36, KLF4, CD247, REEP6, SQRDL, KCNJ6, ANXA2, SPRY2, KCNS3, ITM2C, THBS2, CTD-2515H24.2, LL21NC02-21A1.1, LOC200772, RP11-402C9.1, LINC01203, RP11-479G22.8, RP11-615I2.1, RP1-249I4.2, RP11-116N8.2, and RP11-288L9.4) showed significant differences in the different recurrence time groups [|Log2 (fold-change)| ≥1, P<0.05; Figs. 5 and S3]. These results indicated that lncRNAs may regulate the downstream pathways in NFPA through gene modules and thus play an important role in tumor recurrence.

Figure 5.

Analysis of 16 hub genes and module gene expression levels in NFPA recurrence at ≥5 years (red boxes) or ≤1 year (blue boxes), including (A) CCR1, (B) CCL3, (C) CCL4, (D) PACSIN1, (E) CGNL1, (F) TRIM69, (G) STAB2, (H) ATP8A1, (I) CD48, (J) NOL6, (K) ZFP36, (L) KLF4, (M) CD247, (N) CTD-2515H24.2, (O) LL21NC02-21A1.1 and (P) LOC200772. FC, fold-change; NFPA, non-functioning pituitary adenoma.

Evaluation of the hub and module genes for predicting the recurrence and PFS of patients with NFPA

Next, the predictive ability of the module genes for the recurrence process was evaluated. Kaplan-Meier analysis of the central or module genes, NOL6, CDK15, MOV10, SAMM50, COL24A1, EPHX1 and DCP1A, showed that the patients could be divided into two groups with different risk in the recurrence, according to the median value of each gene expression as the cut-off value. Compared with patients at low risk, the PFS time of patients at high risk was significantly shorter (P<0.05; Fig. 6A-G).

Figure 6.

A total of seven hub genes and module genes were evaluated for predicting the progression-free survival of patients with non-functioning pituitary adenoma, including (A) NOL6, (B) MOV10, (C) CDK15, (D) DCP1A, (E) SAMM50, (F) COL24A1 and (G) EPHX1. Log-rank or Renyi test P<0.05 was considered statistically significant. (H) Relative expression of NOL6 and LL21NC02-21A1.1 by reverse transcription-quantitative PCR analysis, which results basically matched those of the microarray analyses. *P<0.05.

RT-qPCR was conducted to confirm the reliability of the expression profiles generated using the microarray and DEG analyses. Among the prognostic hub and module genes aforementioned, NOL6 and LL21NC02-21A1.1 were randomly selected for verification (P<0.05; Fig. 6H). As expected, the RT-qPCR results basically matched those of the microarray analyses. These results indicated that the bioinformatics analysis of the microarray data reliably identified critical candidate genes involved in NFPA recurrence.

Discussion

NFPAs are pituitary adenomas without clinical evidence of hormonal hypersecretion, and they have a prevalence of 7 to 41.34 cases/100,000 and an annual incidence of 0.65 to 2.34 cases/100,000, in epidemiological data from global between 1955 and 2014 (7,45,46). Transsphenoidal surgery is the recommended first-line treatment (47). However, compared with that in functioning pituitary adenoma, monitoring of the tumor recurrence of NFPA through specific serum hormone alterations is difficult. Therefore, the present study aimed to develop a new predictive signature that could be used as a prognostic prediction model to identify early recurrence. The main purpose of the study was to divide patients into high-risk or low-risk groups so that the most effective and timely treatment can be performed for NFPA.

Numerous studies have focused on the factors of tumor recurrence of NFPA to improve the prognosis of patients postoperatively. Age is recognized as an important independent factor influencing the prognosis of NFPA, and a younger age indicates a greater chance of tumor recurrence (12,48). Ki-67 is another commonly used pathological prognostic evaluation index (49), although a single indicator used in prognostic assessment has certain limitations in accurately evaluating the prognosis of each patient. A previous study tried to establish a statistical model that combined clinical features (age and tumor volume) and molecular markers (p16, WIF1 and TGF-β) to evaluate the recurrence probability of patients with NFPA postoperatively (50). Moreover, compared with a previous study (51), the current study added a temporal component to the prognostic assessment and independently assessed the prognosis of patients at different time points.

In recent years, lncRNAs have been reported in various tumors and serve as promising new molecular markers for tumor biological behavior, tumor diagnosis and prognostic evaluation (52,53). For example, lncRNA H19 is decreased in pituitary adenomas, and its overexpression could markedly inhibit the growth of pituitary tumor cells and be used as a drug resistance marker (54). Xing et al (55) identified differentially expressed mRNAs and lncRNAs in NFPA and normal pituitary tissue samples, and constructed an mRNA-lncRNA coexpression network. However, the research failed to illustrate the regulatory mechanisms of the key genes or lncRNAs and their influence on patient prognosis. The current study focused on identifying molecular markers of NFPA recurrence.

First, the DEGLs based on NFPA recurrence at <1 year and >5 years were obtained. According to GSEA, these DEGs were enriched in the ‘regulation of cell death’ and ‘cell adhesion’. The present results are consistent with those of previous studies, which have shown that intercellular adhesion and adhesion molecules play a crucial role in tumor recurrence and proliferation (56,57).

Second, a total of eight modules were identified via cluster analysis using the PPI network based on the DEGLs. GO and KEGG enrichment analyses illustrated that these module genes were mainly involved in different GO functions and pathways. For module 1, the related GO functions were ‘T cell migration’ and ‘chemotaxis’, which implied that the process of recurrence may be associated with the immune-related tumor microenvironment. Similar to the current study, Marques et al (58) found that a low CD8:CD4 ratio is associated with a higher proliferative index (Ki-67) in pituitary adenoma. In addition, the present KEGG analysis of other modules found that these genes were involved in the ‘cell cycle’, ‘TNF signaling pathway’, ‘VEGF signaling pathway’ and ‘TGF-β signaling pathway’. These pathways might participate and regulate the proliferation and recurrence processes that occur in NFPAs.

Third, the current analysis obtained hub genes and module lncRNAs with significant differential expression (CCR1, CCL3, CCL4, PACSIN1, CGNL1, TRIM69, STAB2, ATP8A1, CD48, NOL6, ZFP36, KLF4, CD247, REEP6, SQRDL, KCNJ6, ANXA2, SPRY2, KCNS3, ITM2C, THBS2, CTD-2515H24.2, LL21NC02-21A1.1, LOC200772, RP11-402C9.1, LINC01203, RP11-479G22.8, RP11-615I2.1, RP1-249I4.2, RP11-116N8.2 and RP11-288L9.4). As an example, ANXA2 is a pleiotropic calcium-dependent phospholipid-binding protein that is abnormally expressed in a variety of cancer types (59), including prostate cancer (60) and liver cancer (61). Liu et al (62) performed a meta-analysis and indicated that ANXA2-overexpression might be related to poor outcomes in patients with malignant tumors, which is consistent with the present findings. In addition, the current study reported that lncRNAs could be used as a prognostic signature. However, the functions and regulatory mechanisms of lncRNAs in NFPA have not yet been reported.

Finally, the present study also assessed the predictive ability of the module genes (such as NOL6, CDK15, MOV10, SAMM50, COL24A1, EPHX1 and DCP1A) for the recurrence process. In addition, the expression levels of NOL6 and LL21NC02-21A1.1 were validated, which were randomly selected from among the hub and module genes, using RT-qPCR. The results confirmed the accuracy of the bioinformatics analyses.

NOL6 encodes a nucleolar RNA-associated protein that is associated with the early stage of ribosome biosynthesis (63). Dong et al (64) found that NOL6 is highly expressed in human prostate cancer and that knockdown of NOL6 inhibits the proliferation and mitosis, and increases the apoptosis of human prostatic carcinoma cells (PC-3). In the current study, NOL6 was found to be upregulated in NFPAs that became recurrent within 1 year compared with those recurring after >5 years, suggesting that NOL6 could be a critical gene in prognostic development and a potential target for NFPA treatment. MOV10 belongs to the RNA helicase superfamily of proteins and could regulate mRNA stability and translation (65). Nakano et al (66) demonstrated that the mRNA and protein levels of MOV10 in cancer cells, such as human leukemia and human cervical carcinoma cells, were higher compared with those in normal cells. In addition, MOV10 has been revealed to promote the angiogenesis of glioma by binding circ-DICER1 (51). These studies indicated that MOV10 could be critical in tumorigenesis. DCP1A is a protein-coding gene for mRNA-decapping enzyme 1a, and several studies have revealed that DCP1A is upregulated in tumor tissues, such as malignant melanoma, colorectal carcinoma and gastric cancer (67–69). In addition, Tang et al (67) and Wu et al (68) reported that the high expression of DCP1A in colorectal carcinoma is correlated with poor prognosis, which is consistent with the current results, thus indicating that the other present PCGs and lncRNAs could also be prognostic indicators for NFPA.

A few limitations of the current study need to be acknowledged. First, the molecular mechanisms underlying the action of these PCGs and lncRNAs in NFPA are still unclear, and further studies might provide important information to understand their functional roles. Second, sequencing data for NFPA are limited; thus, it was not possible to verify the results in an independent validation set. Third, the limited number of samples used for RT-qPCR testing make it necessary to perform larger scale experiments in the future. Finally, the application of the present signature in clinical practice should be tested prospectively. Despite these limitations, the current study verified a certain association between PCG and lncRNA signatures and regression, which is a potentially powerful prognostic marker of NFPA.

In conclusion, to the best of our knowledge, the present study is first to integrate PCGs and lncRNAs to predict tumor recurrence in patients with NFPA. The current study may provide novel insights into prognostic evaluation and help patients benefit from early intervention.