From October 2007 to July 2014, patients who were diagnosed with NFPA and underwent surgical resection at Beijing Tiantan Hospital were included in this study (n = 66). The mean age of these 66 patients was 51.5 years (range, 25–73), there were 34 males and 32 females, and the median follow-up was 76.5 months (range, 5–122). The clinical and pathological characteristics of all the patients are shown in Supplementary Table S1. Cavernous sinus (CS) invasion was defined by the Knosp grading scale (grade 3 or 4) on preoperative enhanced magnetic resonance imaging (MRI) (Knosp et al., 1993). Postoperative tumor recurrence was defined as recurrence identified from any direction on enhanced MRI from the day of surgery to the end of the follow up; the maximum tumor diameter needed to increase by > 2 mm. According to tumor size, NFPA were divided into microadenoma (<10 mm in diameter), macroadenoma (≥10 mm) and giant adenoma (≥40 mm). The local Ethics Committee approved this study, and informed consent was obtained from each subject.

According to the instructions provided, total RNA was extracted and purified from collected samples using the phenol-free mirVana™ miRNA Isolation Kit (Cat # AM1561; Ambion; Thermo Fisher Scientific, Inc.). A Thermo Scientific^TM NanoDrop 2000 was used to quantify and assess purity of the extracted RNA.

RNA samples were used to generate fluorescence-labeled cRNA targets for the SBC human ceRNA array V1.0 (4 × 180 K) and were subsequently hybridized with slides and scanned in an Agilent Microarray Scanner (Agilent Technologies, Santa Clara, CA, United States) to obtain the data. The raw data was extracted using feature extraction software 10.7 (Agilent Technologies, Inc.). Then, the quantile algorithm provided by the “limma” package (http://bioconductor.org/packages/limma/) of the R program was used to normalize the data.

The “sample” function of R progrom (www.r-project.org/) was used to randomly divided 66 NPFA patients into a training set (n = 33) and a testing set (n = 33). In the training group, univariable Cox proportional hazards regression analysis was performed to determine the association between recurrence-free survival (RFS) and lncRNA expression in each patient. We used a machine learning approach, random survival forests-variable hunting (RSFVH) algorithm, to narrow the scope of the gene set through an iteration procedure, discarding the bottom quarter of lncRNAs (the least important lncRNAs) at each step. In total, nine lncRNAs were selected (Mogensen et al., 2012; Li et al., 2014; Ishwaran and Lu, 2019).

The selected lncRNAs was used to construct a risk prediction score model as follows (Ritchie et al., 2015; Guo et al., 2016).

To verify the existence of the lncRNA signature, twelve samples were randomly selected from the entire group for RT-PCR and agarose gel electrophoresis. LncRNA reverse transcription was performed using a High Capacity cDNA Reverse Transcription Kit (0049472, Thermo Fisher). Next, PCR was performed using I-5TM High-Fidelity Master Mix (I5HM, 200MCLAB). PCR was conducted as follows: 2 min of initial denaturation at 98°C, 32 cycles of 10 s at 98°C, 58°C for 10 s and 72˚C for 10s, and final extension step for 5 min at 72˚C. GAPDH was used as an internal control gene. The PCR products were run on 2% agarose gel and visualized using a UV transilluminator. The primer sequences are presented in Supplementary Table S2.

The survival distribution of different groups was evaluated and compared using Kaplan-Meier survival analyses and two-sided log-rank tests. The chi-square test was used to analyzed the associations with clinical signatures. p < 0.05 was considered to indicate statistical significance. All analyses were performed using R program 3.6.1. The packages were downloaded from Bioconductor, including the survival, ROC, and randomForestSRC packages.

To investigate the potential function of the lncRNAs in the signature, Pearson correlation tests were used to identified protein-coding genes (PCGs) coexpressed with the prognostic lncRNAs. The genes with a p < 0.05 and an absolute value of the Pearson coefficient > 0.6 were selected for Gene Ontology (GO) (Ashburner et al., 2000; The Gene Ontology, 2017) and Kyoto Encyclopedia of Genes and Genomes (KEGG) (Kanehisa and Goto, 2000; Kanehisa et al., 2016; Kanehisa et al., 2017) enrichment analyses. The GO and KEGG analyses were performed with the clusterProfiler package (Yu et al., 2012) of the R program.

A total of 19,741 lncRNAs were extracted from the 66 NFPA expression profiles. The flow chart of this study is shown in Figure 1. The patient information of all patients is summarized in Table 1.

Initially, in the training set, univariate Cox proportional hazards regression analysis was used to obtain RFS-related lncRNAs. The 1,214-lncRNA set was identified using recurrence as the dependent variable, and the signature was significantly associated with patient recurrence (Supplementary Table S3, p value < 0.05, Figure 2A).

Secondly, to further reduce the number of prognostic lncRNAs, the random forest supervised cassification (RFSC) algorithm was employed to analyze the 1,214-lncRNA set, and the nine lncRNAs most related to recurrence were obtained according to the permutation important score calculated with the RFSC algorithm (Figure 2B; Supplementary Figure S1).

Thirdly, based on the nine types of lncRNA, we constructed a risk-score model of 2⁹–1 (511) types of lncRNA set combinations, which contained different lncRNA numbers from 1 to 9. To screen for a better prediction signature, we conducted a time-dependent ROC analysis that used recurrence status as a lable and signature risk scores as a variable in the training group and compared the sensitivities and specificities (Supplementary Table S4).

According to the AUC values of all 511 signatures (Supplementary Table S4), we identified the lncRNA combination composed of LOC101927765, RP11-23N2.4, and RP4-533D7.4 as the most promising, as it had strong ability to predict recurrence and the smallest node and the largest AUC value of 0.87 (Figure 2C; Table 2). RT-PCR was used to confirm the reliability of microarray sequencing. Consisting with the microarray data, the three lncRNAs were detected in 12 tumor tissues (Figure 2D), which revealed that the lncRNA are stable and can be used as prognostic maker.

The risk score of the signature was calculated as follows: risk score = (3.41 × expression value of LOC101927765) + (1.90 ×  expression value of RP11-23N2.4) + (−3.43 ×  expression value of RP4-533D7.4).

Each patient obtains a risk score according to the risk score model. Then, the patients from the training group were divided into a high-risk group (n = 16) and a low-risk group (n = 17) based on the cutoff point, which was the median risk score. Kaplan-Meier survival analysis was performed to determined the difference in RFS between the two risk groups. The median RFS time was significantly shorter in the high-risk group (4.44 years) than in the low-risk group (6.74 years) (p < 0.001; log-rank test, Figure 3A). Moreover, the recurrence rate of the high-risk group was higher than that of the low-risk group (>60% vs. < 1%). In a similar manner, patients from the test group were also divided into two risk groups. The results of Kaplan-Meier analyses for the high-risk (n = 16) and low-risk (n = 17) groups in the test dataset were plotted and are shown in Figure 3B (median RFS time: 5.51 vs. 6.82 years; log-rank test, p = 0.016), and the RFS rates were approximately 52.25 and 87.40%, respectively. In addition, patients in the entire group were similary divided into high-risk (n = 32) and low-risk (n = 34) groups, and Kaplan-Meier analysis further confirmed the ability of the lncRNA signature to predict recurrence (median PFS time: 4.97 vs. 6.79 years; log-rank test, p < 0.001, Figure 3C).

Figures 4A–C intuitively shows the risk score, survival status and expression pattern of lncRNAs in the training, testing, and independent datasets. For patients with low risk scores in the three datasets, RP4−533D7.4 was highly expressed, while LOC101927765 and RP11-23N2.4 was expressed at low levels; the opposite patterns for each lncRNA were seen in patients with high risk scores.

After proving the recurrence prediction ability of the lncRNA signature, we explored the correlation between the signature and clinical characteristics in the entire dataset (n = 66) to understand the clinical significance of the lncRNA signature.

Table 3 shows that there was an association between the lncRNA signature and age in the entire group (chi-suqre test, p = 0.03, Table 3). In addition, we further assessed whether the prognostic value of the three-lncRNA signature was independent of other clinical factors. Univariate and multivariate Cox regression analyses of factors including age, sex, tumor size classification, CS invasion, and the signature were performed. In the entire dataset, age (HR = 0.33, 95% CI = 0.12–0.93, p = 0.04) and the signature risk score (HR = 1.50, 95% CI = 1.24–1.82, p < 0.001) were significantly associated with the RFS of patients (Table 4). Moreover, the three-signature score was also an indenpent prognostic factor associated with RFS in the training (HR = 2.06, 95% CI = 1.36–3.12, p < 0.001) and test set (HR = 6.96, 95% CI = 1.21–40.16, p = 0.03). Hence, the results indicate that the three-lncRNA signature is an independent prognostic factor for NFPA RFS.

It has been reported that age is associated with a risk of tumor recurrence (Losa et al., 2008). ROC analysis was performed to determine the predictive power of the lncRNA signature and age. The results showed that in the training/testing/entire group, the AUC values of the lncRNA signature were lager those of age (AUC = 0.87/0.726/0.798 vs. AUC = 0.683/0.676/0.679, Figures 5A–B), indicating that the signature had high accuracy and important clinical significance. In addition, time-dependent ROC analysis was performed on the three datasets to further understand the signature prediction capabilities for 3-, 4- and 5 year RFS. The signature AUC values in the training/test/entire group at 3, 4, and 5 years, as shown in Figures 5D–F, indicated a strong predictive power of the signature for RFS (AUC = 0.767/0.818/0.833, 0.651/0.723/0.713, and 0.688/0.774/0.769, respectively).

The PCGs correlated with the lncRNAs in our prognostic signature were obtained by Pearson correlation analysis in all 66 patients, and their potential biological function were explored. The expression of 1,056 PCGs was highly correlated with that of at least one of the LncRNAs (Pearson correlation coefficient > 0.60, p < 0.05, Supplementary Table S5). Next, we performed GO and KEGG analyses and found these the genes were enriched in 99 different terms (Supplementary Table S6), such as mRNA processing, RNA splicing and oxidative phosphorylation (Figure 6).