GSE32960, GSE70970, and GSE102349 datasets containing the expression profiles of NPC samples were accessed from the Gene Expression Omnibus (GEO) database (24), where GSE32960 and GSE70970 include miRNA expression data and GSE102349 includes mRNA expression data. We screened the samples of GEO datasets according to the following criteria: 1) remove samples without survival time and survival status; 2) convert probes into gene symbols; 3) remove a probe matching to multiple genes; 4) select the averaged expression of a gene with multiple probes; 5) for miRNA data, only human samples were remained. After preprocessing, 312 NPC and 18 normal samples were remained in the GSE32960 dataset; 253 NPC and 10 normal samples were remained in the GSE70970 dataset; and 88 NPC samples were remained in the GSE102349 dataset.

In the GSE32960 dataset, differentially expressed miRNAs (DEmiRNAs) were screened by Limma R package (25) from NPC and normal samples under conditions of P < 0.05 and |log2(fold change, FC)| > 1.2. Gene modules were identified by weighted correlation network analysis (WGCNA) (26). Firstly, samples were clustered and the co-expression network was constructed. A scale-free network was ensured under scale-free R2 = 0.85, and soft threshold (power) = 3 was determined. The co-expression network was then converted to the adjacent matrix and was further converted to the topological overlap matrix (TOM). Subsequently, we clustered genes by the dynamic cutting method and average-linkage hierarchical clustering based on TOM. Eigengenes were calculated for each gene, and gene modules were clustered and merged under parameters of deepSplit = 2, and minModuleSize = 60, height = 0.25. The module–trait relationships were assessed by Pearson correlation analysis. By overlapping the miRNAs in the NPC-associated gene modules and DEmiRNAs, NPC-associated miRNAs were determined.

Random grouping of samples from the GSE32960 dataset into training group and test group at a ratio of 3:2 was performed. Two-group differences were assessed using Student’s t test. Univariate Cox regression analysis screened NPC-associated miRNAs from the training group. MiRNAs with P < 0.01 were selected as prognostic miRNAs. Multivariate analysis (stepAIC) was introduced to measure the coefficients of prognostic miRNAs. Then, the miRNA-based risk model was defined as: risk score = Σ(coef i*expression i), where coef indicates the coefficients of miRNAs and i represents miRNAs.

Each sample obtained a risk score calculated by the risk model. The median risk score was employed in dividing samples into low-risk and high-risk groups. Receiver operating characteristic (ROC) curve analysis was used to predict the efficiency of the risk model in predicting overall survival through timeROC R package (27). The prognosis difference of two risk groups was studied by Kaplan–Meier survival analysis. Univariate and multivariate Cox regression models were used to analyze the hazard ratio of risk type. A nomogram was used to optimize the clinical use of the risk model with the rms package.

First of all, the potential target genes of miRNAs were predicted by different online tools and software including microT (28), miRanda (29), mircode (30), miRDB (31), miRmap (32), miRtarbase (33), PicTar (34), PITA (35), TargetMiner (36), TargetScan (37), RNA22 (38), and starbase (39). The target genes predicted by at least six tools were remained as key potential target genes. The WebGestaltR package (40) was utilized to annotate the enriched Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways and Gene Ontology (GO) terms. Then, we used the key target genes to establish the mRNA-based prognostic model. Prognostic target genes in 70% samples of the GSE102349 dataset were screened by univariate Cox regression analysis. Random sampling with 70% samples in the GSE102349 dataset was performed for 1,000 times corresponding with univariate analysis. The top five frequent target genes from the results of 1,000 times of univariate analysis were selected as the final target genes for constructing the mRNA-based prognostic model (defined by the same formula with the miRNA risk model).

There were 31 SECRETORY_PATHWAYS screened from the KEGG PATHWAY Database (https://www.genome.jp/kegg/pathway.html). SECRETORY_PATHWAYS scores in the GSE102349 dataset were calculated by ssGSEA (41). There were 26 SECRETORY-related genes obtained from KEGG (https://www.genome.jp/dbget-bin/www_bfind_sub?mode=bfind&max_hit=1000&locale=en&serv=kegg&dbkey=genes&keywords=Secretory&page=1). Pearson analysis was conducted to determine the correlation between SECRETORY_PATHWAYS scores/SECRETORY-related genes and RiskScore.

We used ssGSEA (41) to estimate the immune cell proportion in two risk groups. The gene sets of 22 immune-related cells, innate immune response, and adaptive immune response were obtained from Charoentong et al. (42). The ESTIMATE algorithm (43) was employed to measure immune cell infiltration and stromal cell infiltration. The ESTIMATE score represents the combined score of immune score and stromal score. MCP-counter (44) was used to evaluate 10 immune-related cells including monocytic lineage, CD3+ T cells, NK cells, CD8+ T cells, endothelial cells, cytotoxic lymphocytes, B lymphocytes, neutrophils, fibroblasts, and myeloid dendritic cells based on the expression matrix. The immune checkpoint genes were downloaded from a previous study (42).

Biological pathways were accessed from “h.all.v7.4.symbols.gmt” downloaded from the Molecular Signatures Database (MSigDB) (45). There were 13 tumor-related genes obtained from a previous research (46), which are classic cancer pathways, involved in the development and progression of cancer, including DNA damage repair (DDR), epithelial–mesenchymal transition (EMT), cell cycle, mismatch repair, CD8 T effector, FGFR3, nucleotide excision repair, base excision repair, DNA replication, homologous recombination, and WNT target. ssGSEA was performed to determine the enrichment score of biological pathways. The relation of risk score with pathways was inspected with Pearson correlation analysis.

We conducted SubMap analysis (47) to compare the expression profiles between GSE102349 and IMvigor210. IMvigor210 contains the expression profiles of patients with metastatic urothelial carcinoma treated by PD-L1 inhibitors (48). The higher similarity of samples in GSE102349 with complete response (CR) and partial response (PR) groups in IMvigor210 suggested that the samples were more sensitive to anti-PD-L1 treatment. The TIDE (http://tide.dfci.harvard.edu/) algorithm (49) was employed to predict the escape and response to immune checkpoint inhibitors, according to the score of T-cell dysfunction and exclusion, the enrichment of immunosuppressive cells. The sensitivity of two risk groups to chemotherapeutic drugs was estimated using the pRRophetic R package (50).

The performance of the prognostic model was further examined in immunotherapy datasets including IMvigor210, GSE135222, and GSE78220. GSE135222 contains non-small cell lung cancer patients treated with immune checkpoint inhibitors (51). GSE78220 includes patients suffering from metastatic melanoma treated by anti-programmed cell death protein 1 (anti-PD-1) therapy.

After resuscitating NPC cells and normal cells, they were placed in 1,640 and 5A cell medium containing 10% fetal bovine serum, respectively, at 37°C, 5% CO₂ concentration, and constant temperature to around 80%–90%, then passed and spread on a plate.

Total RNA was extracted from cells by RT-qPCR, and cDNA was synthesized by reverse transcription. After reverse transcription, samples were added according to the experimental instructions. The reaction conditions of RT-qPCR were 95°C for 30 s. 95°C, 5 s; 60°C, 30 s, 40 cycles. The mRNA expressions of miR-142-3p and miR-93 in the experimental group and control group were analyzed.

The statistical methods used in this study were performed in R software (v4.2.0). The Sangerbox platform (52) was used to provide an assistant in data analysis. The log-rank test was applied in survival analysis. Difference between two groups was examined by the Wilcoxon test. The Kruskal–Walls test was used to test the difference among over two groups. P < 0.05 was determined as statistically significant.

To begin with, we assessed the expression difference of miRNAs between normal and tumor samples in the GSE32960 dataset. The results showed that 332 miRNAs were differentially expressed between normal and tumor groups, including 168 upregulated and 164 downregulated miRNAs in tumor groups ( Figure 1A ; Table S1 ). Then, we applied WGCNA to cluster samples and dig out key gene modules based on DEmiRNAs. To ensure the co-expression network to be a scale-free network, the Pearson coefficient was selected as 0.85 and the power of the soft threshold was determined as 3 to construct an adjacency matrix ( Figures 1B–D ). Next, a topology overlap matrix (TOM) was generated based on the adjacency matrix and the genes were divided into different modules using the Dynamic Tree Cut algorithm ( Figure 1E ). Four modules were finally determined after merging the adjacent modules. Then, we analyzed the correlation of four modules with different sample groups. As a result, brown and blue modules were found to be evidently associated with sample groups and they manifested opposite correlations with two groups ( Figure 1F ). Specifically, blue and brown modules were negatively correlated with the tumor group (R = -0.35 and -0.53, respectively). The Venn plot was constructed to describe the intersection between DEmiRNAs and miRNAs of blue and brown modules ( Figure 1G ). There were 99 DEmiRNAs including 42 upregulated and 55 downregulated found in the blue module. There were 169 DEmiRNAs including 79 upregulated and 90 downregulated found in the brown module. The above total 268 DEmiRNAs were determined as potential miRNAs associated with NPC.

We randomly divided 312 NPC samples of the GSE32960 dataset at a ratio of 3:2 into the training and test groups ( Table S2 ). We used the training group to screen prognostic miRNAs from 268 DEmiRNAs according to the univariate Cox regression model. The analysis identified two miRNAs (hsa-miR-142-3p and hsa-miR-93) that were significantly associated with the overall survival (P < 0.01). Based on the multivariate coefficients of two miRNAs by stepAIC analysis, we established the risk model defined as follows ( Figure S1A ).

Risk score = -0.835*(hsa-miR-142-3p) + 0.85*(hsa-miR-93)

Moreover, the RT-qPCR analysis results showed that miR-142-3p expression was downregulated and miR-93 was upregulated in five NPC cell lines in comparison with normal NP69 cells ( Figures S1B, C ).

To validate the performance of the risk model, we calculated the risk score for each tumor sample in the training and test groups. The AUC for 1-, 3-, and 5-year survival derived from the ROC curve analysis was 0.56, 0.70, and 0.70 in the training group and 0.96, 0.71, and 0.71 in the test group, respectively ( Figures 2A, B ). Furthermore, based on the median risk score, tumor samples were classified into low-risk and high-risk groups. The two risk groups had distinct overall survival (OS) in the training and test groups, as shown by Kaplan–Meier survival analysis (P = 0.00035 and P = 0.013, respectively, Figures 2A, B ).

In the total GSE32960 dataset, the risk model still showed good performance in predicting overall survival ( Figure 2C ). In addition, we evaluated the effectiveness of the risk model in different survival times including recurrence-free survival (RFS), disease-free survival (DFS), and metastasis-free survival (MFS). High-risk and low-risk groups exhibited differential prognosis of DFS and MFS (P < 0.0001) but showed no significant difference in RFS classification (P = 0.063, Figure 2C ). In another independent dataset (GSE70970), the risk model also showed an effect classification for both OS and DFS (P = 0.017 and P = 0.022, respectively) ( Figure 2D ). Moreover, we verified the effectiveness of the risk model in the groups of different clinical characteristics. In addition to the female group and N0 group, the risk model was sufficient to distinguish samples into different risk levels in the groups of age >45, age ≤45, male, T1–T2, T3–T4, N1–N3, I–II, and III–IV ( Figure 3 ).

Using univariate and multivariate Cox regression analyses, we evaluated the relation of clinical features and risk type with prognosis. The result displayed that only gender and risk type were independent risk factors in both univariate and multivariate analyses ( Figures 4A, B ). Gender had a hazard ratio (HR) of 2.2 in both univariate and multivariate analyses (P = 0.019 and P = 0.02, respectively). Risk type had HR of 3.7 (P = 1.1e-6) and 3.6 (P = 2.1e-6) in univariate and multivariate analyses, respectively. Therefore, we included gender and risk score to construct the nomogram ( Figure 4C ). Compared with gender, risk score contributed the more total points in the nomogram. Calibration curve analysis suggested that the predicted 1-, 3-, and 5-year survival was highly accordant with the observed survival ( Figure 4D ). In addition, decision curve analysis was implemented to evaluate the benefit that patients may obtain from gender, risk score, and nomogram. As a result, nomogram performed a more favorable performance than gender and risk score ( Figure 4E ). Consequently, the nomogram based on gender and risk score was more efficient to predict the prognosis of NPC patients.

In the above sections, we identified two key miRNAs (hsa-miR-142-3p and hsa-miR-93) that had a close relation with NPC prognosis. To understand the potential molecular mechanism of the miRNAs, we applied 12 tools (starbase, PITA, TargetMiner, miRanda, microT, miRmap, miRtarbase, mircode, TargetScan, RNA22, miRDB, PicTar) to predict the potential targets of two miRNAs. The results output 39 target genes of hsa-miR-142-3p and 86 target genes of hsa-miR-93, which were visualized in ceRNA networks ( Figure 5A ). KEGG analysis demonstrated the involvement of these target genes in mitophagy and autophagy ( Figure 5B ). The top 10 enriched terms of cellular component and molecular function, as well as biological process, were visualized using GO function analysis ( Figures 5C–E ). For example, biological process terms of positive regulation of amyloid−beta metabolic process, amyloid−beta formation, and negative regulation of phosphatase activity were enriched ( Figure 5C ). Cellular component terms of clathrin-coated pit and trans-Golgi network membrane were enriched ( Figure 5D ). Molecular function terms of clathrin adaptor activity and clathrin heavy chain binding were enriched ( Figure 5E ).

We predicted a total of 125 potential target genes of hsa-miR-142-3p and hsa-miR-93. Then, we used the univariate Cox regression model to analyze the relation between target genes and overall survival. Random sampling for 1,000 times from the samples of the GSE102349 dataset was performed. The target genes closely related to overall survival were remained (P < 0.05) and were ranked by the occurring frequency from 1,000-times analysis. The top five frequent target genes were selected as key target genes for establishing the mRNA prognostic model ( Figure 6A ). The coefficients of five target genes were calculated by multivariate analysis. Finally, the prognostic model was defined as the following: mRNA risk score = 1.333* E2F1 - 1.766*KCNJ8 + 1.075*SUCO - 1.030* HECTD1 - 0.340*KIF23.

The risk score for each sample in the GSE102349 dataset was determined with the mRNA prognostic model. ROC curve analysis showed that the model was efficient in predicting 1-, 3-, and 5-year survival, with AUCs of 0.87, 0.81, and 0.79, respectively ( Figure 6B ). The median risk score value was used in classifying NPC samples into two groups, high-risk and low-risk groups. Kaplan–Meier survival curves of two risk groups showed that they had an apparently different prognosis (P = 0.0018, Figure 6B ). Hence, the five target genes could be validated by the above results to be closely involved in the prognosis.

The tumor microenvironment plays a central role in antitumor response and immunotherapeutic response. We used several tools to assess the immune cell component, as well as immune and stromal infiltration, and analyzed immune checkpoint genes for their expression levels. Two risk groups showed a significantly different immune microenvironment. We estimated an ssGSEA enrichment score of 28 immune-related cells and found that 25 of them had a differential enrichment score, such as myeloid-derived suppressor cells (MDSCs), activated CD8 T cells, regulatory T cells, natural killer cells, activated B cells, and macrophages ( Figure 7A ). Most immune cells showed a higher abundance in the group with a low risk. In addition, ssGSEA also revealed higher enrichment of both adaptive and innate immune response scores in the low-risk group (P < 0.0001, Figure 7B ). The ESTIMATE algorithm was used to evaluate stromal and immune infiltration of two groups, and not surprisingly, the low-risk group showed both higher immune score and stromal score than the high-risk group ( Figure 7C , P < 0.0001). Moreover, MCP-counter was employed to dig out similar results with ssGSEA. A total of 10 types of immune-related cells all showed a higher enrichment score in the low-risk group ( Figure 7D ). Immune checkpoint expression levels also had an extreme difference in two risk groups that most of immune checkpoints, such as CTLA4, TIGIT, LAG3, and PCDC1, were more highly expressed in the low-risk group than in the high-risk group ( Figure 7E ). Immune checkpoints’ differential expression may contribute to different antitumor immune responses. The distinct immune microenvironment suggested that the five prognostic genes may play critical roles in immune modulation.

To further understand the molecular mechanisms contributing different prognoses of two risk groups, the enrichment score of pathways from the “h.all.v7.4.symbols.gmt” file was calculated using ssGSEA. There were 30 pathways differentially activated in the two risk groups ( Figure 8A ). Immune-related pathways, for instance, interferon alpha response, interferon gamma response, inflammation response, IL2-STAT5 signaling, IL6-JAK-STAT3 signaling, and complement, were significantly more activated in the low-risk group, which was consistent with the result of immune analysis. Cell cycle-related pathways such as E2F targets, MYC target V2, MYC target V1, and G2M checkpoint were less enriched in the low-risk group. In addition, some oncogenic pathways were more activated in the high-risk group, such as Wnt signaling and Hedgehog signaling. Notably, apoptosis was more enriched in the low-risk group, which was associated with good prognosis. Risk score also manifested significant correlations with the above pathways ( Figure S2 ).

Gene sets of 13 tumor-related pathways were obtained from a previous study, and their enrichment scores were calculated using ssGSEA. Pearson correlation analysis uncovered a negative association of risk score with DNA repair-related pathways including DDR (R = -0.41), base excision repair (R = -0.25), nucleotide excision repair (R = -0.33), homologous recombination (R = -0.33), and mismatch repair (R = -0.33) ( Figure 8B ).

Next, we evaluated the response of two risk groups to chemotherapy and immunotherapy. The similarity of expression profiles between GSE102349 and IMvigor210 (treated by PD-L1 inhibitors) datasets was shown by SubMap analysis. Higher similarity between two datasets indicates higher sensitivity to PD-L1 inhibitors. The results presented that the low-risk group had a higher similarity with CR and PR groups (P < 0.05, Figure 8C ), implying that the low-risk group could obtain more benefit from immunotherapy. Furthermore, the TIDE algorithm was implemented to predict immune escape to immune checkpoint inhibitors. Two risk groups showed no significant difference of TIDE score ( Figure 8D ). However, the low-risk group exhibited more severe T-cell dysfunction than the high-risk group in which enrichment of MDSCs and M2 tumor-associated macrophages (TAM) was higher, contributing to its higher T-cell exclusion. In the response to chemotherapeutic drugs, we screened a total of 103 drugs including 46 drugs such as YM155 and vinorelbine sensitive to high-risk groups and 57 drugs such as sunitinib and temsirolimus sensitive to the low-risk group ( Figure 8E ).

Abnormal endocrine metabolism is one of the complications of tumor treatment. As shown in Figure 9A , endocrine-related genes, such as MON1B, SCAMP2, and FAM20A, were negatively associated with RiskScore; SCAMP3 was positively correlated with RiskScore. Pearson analysis between endocrine pathways and RiskScore showed that most endocrine pathways were negatively related to RiskScore ( Figure 9B ).

To further examine the exhibition of the five-gene prognostic model, we introduced three independent datasets containing patients administrated by immunotherapy (IMvigor210, GSE135222, and GSE78220). Using the same calculation method, we measured the risk score of each patient in three datasets. In the IMvigor210 dataset, the high-risk group displayed apparently worse prognosis than the low-risk group (P < 0.0001, Figure 10A ). The prediction efficiency was evaluated by ROC, with AUCs of 0.61, 0.66, and 0.64 at 0.5-, 1-, and 1.5-year survival, respectively ( Figure 10A ). In addition, we analyzed the proportion of different response groups in two risk groups. The low-risk group had a higher proportion of CR and PR groups (11% and 20%, respectively), compared with the high-risk group (6% and 9%, respectively). The CR and PR groups also showed a lower risk score than PD and SD groups. In the GSE135222 and GSE78220 datasets, we observed correspondent results ( Figures 10B, C ). The GSE135222 dataset showed AUCs of 0.82 and 0.85 at 0.5- and 1-year survival, respectively. The low-risk group exhibited a markedly higher percentage of CR/PR patients compared with the high-risk group (50% versus 8% in low- versus high-risk groups). The GSE78220 dataset showed AUCs of 0.74 and 0.71 at 1- and 2-year survival, respectively. Expectedly, the CR and PR groups were more accumulated in the low-risk group (21% and 43%) compared with the high-risk group (8% and 31%). Moreover, the CR/PR group showed a lower risk score than the PD/SD group in both GSE135222 and GSE78220 datasets, but the difference was not significant in the GSE78220 dataset. The above observation suggested that patients with a low risk score were more sensitive to immunotherapy and could attain longer survival.