10 NPC tumor tissue samples were obtained from the patients who were diagnosed and treated at the Fujian Cancer Hospital between May 9, 2013, and August 2, 2016. All 10 patients met the following eligibility criteria: newly diagnosed NPC, received standardized radiotherapy and chemotherapy, ≥18 years old, adequate hematological, renal and hepatic functions, and no other malignant diseases. All the patients provided written informed consent. The study was approved by the Ethics Committee of Fujian Cancer Hospital and Fujian Medical University Cancer Hospital (approval number SQ2019-035-01). The tissue samples were stored in liquid nitrogen for subsequent RNA extraction. During the 5 years of follow-up, 5 samples were from patients with disease progression after radiotherapy and chemotherapy, containing 3 liver metastases, 1 bone metastasis, and 1 lung metastasis. While the other 5 samples were evaluated as having a complete or partial response after the treatment.

As an external validation cohort, RNA-seq data of NPC from the GEO database (https://www.ncbi.nlm.nih.gov/geo/, GSE102349) were selected to verify the reliability and applicability of the data of this study (19, 20). The single-cell dataset GSE150430 was designed to validate the accuracy of the model at the individual cell level and to probe the communication of cells and ligand receptors in the immune microenvironment of NPC. Also, the tumor tissue biopsies of 74 NPC patients treated in our hospital in 2021 and 2022 were used for the immunohistochemistry to validate CD8 T cell infiltration and immune checkpoints expression, including 11 cases in the metastatic group and 63 cases in the non-metastatic group ( Supplementary Table 1 ). The 8th edition of the American Joint Committee on Cancer (AJCC) Staging Manual was used to restage all the patients.

NPC biopsies were fixed with 10% formalin overnight and processed into 5-μm-thick paraffin sections. The slides were then analyzed by immunohistochemistry with anti-human CD8 (Cat # ab237709; Abcam), anti-human PD1 (Cat # ab52587; Abcam), and anti-human PD-L1 (Cat # ab213524; Abcam) followed by HRP secondary antibody (Cat #ab205718; Abcam) and DAB staining. Images were obtained using a microscope (BX43; Olympus, Japan). Histochemistry score (H-score) was used to evaluate the expression. H-score = (percentage of cells of weak intensity × 1) + (percentage of cells of moderate intensity × 2) + (percentage of cells of strong intensity × 3).

The R package “ggplot2” was employed to visualize DEGs from sequencing data of NPC samples in Fujian Cancer Hospital (21). The cut-off values met the following two conditions: fold-change of >2 and the p-value of<0.05. Gene ontology (GO) (22, 23) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses (24) were applied to further explore the pathways of DEGs enrichment. A false-discovery rate of<0.05 was set as the cut-off value. The immune gene data was downloaded through the ImmPort data portal (www.immport.org/immport-open/public/home/home), and 2,498 immune-related genes were obtained. Then the intersection of the DEGs and the immune-related genes was selected as differentially expressed immune-related genes. Progression-free survival (PFS) was subjected to minimum absolute shrinkage and selection operator (LASSO) Cox regression with 10-fold cross-validation to screen for DEGs with prognostic value on the basis of the univariate Cox analysis. The R package “glmnet” was employed to determine the gene signatures containing the biomarkers most helpful for prognosis (25). The prognosis risk score was established by linearly combining the following formula:

where exp denotes the gene expression value, while coef refers to the coefficient of a gene in LASSO analysis.

To assess the predictive power of our prognostic risk model, receiver operating characteristic (ROC) for 1- and 3-year survival were performed in the validation cohort GSE102349 using the R package “timeROC”. Next, the samples were divided into high-risk and low-risk groups according to the best cut-off value of the risk score from the R package “survival” for survival analysis. The survival curves were compared using the Kaplan-Meier method and the log-rank test. The univariate and multivariate Cox regression models were applied to determine whether the risk score was an independent prognostic factor.

To evaluate the infiltration of immune cells from several aspects, we adopted multiple immune scoring approaches, like TIMER and ssGSEA algorithms (26, 27). The immune scores and tumor purity were estimated by the R package “ESTIMATE” (28). From an earlier study, we retrieved a group of six inhibitory immune checkpoints that displayed immune therapeutic efficacy (29). Gene sets that displayed T cell-inflamed gene expression profile (GEP) and tertiary lymphatic structure (TLS) were acquired (30, 31). Furthermore, we assessed the enrichment of 10 oncogenic pathways using the ssGSEA method (32). The score of activation minus the score of repression represented the final score of each pathway. We used a validated set of 31 genes related to cell cycle progression (CCP) to estimate the rate of cell proliferation (33). The cluster score was calculated as the average expression level of CCP-related pathways by subtracting the mean level.

To assess the predictive efficacy of the model for immunotherapy efficacy, we collected several immunotherapy cohorts from the GEO database and the TIGER website (http://tiger.canceromics.org/#/), including nasopharyngeal carcinoma-GSE102349, melanoma-GSE91061, melanoma-PRJEB23709, NSCLC-GSE126044. We visually compared the proportion of patients with and without response to immunotherapy in high- and low-risk groups.

This study performed quality control, downscaling, and clustering of scRNA-seq data as well using Seurat (v.4.0.4) (34). To ensure data quality, genes detected in less than 3 cells and cells with less than 250 genes detected were excluded, and the percentage of mitochondria was limited to less than 35% (35). Data were processed by the logNormalize method for normalization. The nonlinear dimensionality reduction method Uniform Manifold Approximation and Projection for Dimension Reduction (UMAP) was utilized for unsupervised classification and unbiased visualization of cell populations on two-dimensional maps (36). TISCH (http://tisch.comp-genomics.org/) provides detailed cell type annotations at the single-cell level (35). After that, the “FindAllMarkers” function was configured to identify marker genes in each cluster using a filter value of absolute log2 fold change (FC) ≥ 0.3 and a minimum cell cluster fraction of 0.25.

A risk score of each single cell sample from GSE150430 was calculated by single sample Gene Set Enrichment Analysis (ssGSEA) method and was completed using the “GSVA” and “GSEABase” packages in R. We used single-cell data as a reference, apply a newly developed deconvolution algorithms (CIBERSORTx) to the bulk transcriptome data to quantitatively estimate cell-type proportions for each tumor in GEO database (37).

CellChat v1.1.3 software inferred cell-cell communication based on ligand-receptor interactions (38). Cell groups with less than 10 cells were filtered out of cell-cell communication. Pairwise tests were performed on communication probability values to assess their statistical significance.

Statistical analysis was done using R software (V.3.6.1) and SPSS software (ver. 25.0). Wilcoxon rank sum test and chi-square test were conducted for continuous and categorical variables, respectively. For all analysis, two-by-two pairs indicate statistically significant differences. *, **, *** and **** indicate, respectively <0.05, <0.01, <0.001, and <0.0001.

The schematic diagram presents the workflow of our study ( Figure S1 ). The RNA-seq profiles were generated for the NPC samples of 10 patients treated at the Fujian Cancer Hospital, 5 of whom were assigned to the non-metastasis group, while the other 5 were in the metastasis group owing to distant metastasis after chemoradiotherapy. The baseline characteristics of patients in the metastatic and non-metastatic groups could be seen in Table 1 (n=10). In general, PCA indicated distinct transcriptional profiles between the metastatic group and the non-metastatic group ( Figure S2A ). Then, 1328 DEGs were conspicuously illustrated in the volcano plot, with 520 upregulated genes and 808 downregulated genes ( Figure 1A ). The KEGG and GO analyses are the universally applicable statistical methods of enrichment analysis. The DEGs were enriched in the immune-related pathways of the bubble chart containing signal transduction, adaptive immune system, innate immune system, and hemostasis ( Figure 1B ). Simultaneously, they were also centralized in the cell periphery, plasma membrane, and immune system processes ( Figure 1C ). The expression levels of the top 154 immune-related genes selected from the DEGs can be significantly distinguished between the two groups in the heat-map ( Figure 1D ). Overall, immune gene expression and immune signaling pathway were down-regulated in the metastasis group, indicating a potential “immune-cold” tumor phenotype in the metastasis group. For the validation, Therefore, we performed immunohistochemistry staining of CD8 T cell, PD1, and PD-L1. Our immunohistochemistry analysis showed that PD1 and PD-L1 expressions were down-regulated, and the infiltration of CD8 T cells was decreased in the metastasis NPC group (n=11) compared to the non-metastasis group (n=63, Figures 1E, F ).

The LASSO logistic regression model was applied to establish the prognostic immune biomarkers, which involved 9 hub genes (A2M, APLNR, CD8B, RAC3, PRDX2, ULBP1, TMSB15B, KIR3DL2, and SEMA4F; Figure 2A ). The standard for high and low risk scores was evaluated based on cut points associated with the median risk score. Cut-off value of 1.31 for the risk model was identified, which served to divide the patients into high-risk group (with levels of risk score ≥ 1.31) and a low-risk group (with levels of risk score< 1.31). The risk scores were significantly distinguished between the clinical stages I–III and stage IV in GSE102349 ( Figure S2B ), which indicated that the clinical stage of the tumor could be one of the critical factors in assessing the effect of the treatment. The risk scores were also apparently different between the metastasis and non-metastasis groups in our hospital cohort ( Figure S2C ). Patients in the high-risk group had worse tumor metastatic presentation, which is indicative of a worse prognosis ( Figure 2B ). This finding was further validated in a cohort of patients from Fujian Cancer Hospital ( Figure 2C ).

It was found that APLNR, KIR3DL2, CD8B, and A2M were upregulated in the low-risk group, while PRDX2, ULBP1, TMSB15B, SEMA4F and RAC3 were upregulated in the high-risk group (Figure 2D). The assumption could be proposed that the former 4 genes were protective biomarkers, while the latter 5 genes were risk biomarkers. The area under the ROC curve (AUC) was 0.79 at 1-year, and 0.81 at 3-years, respectively, indicating a high predictive value ( Figure 2E ). Combining the results of univariate ( Figure 2F ) and multivariate ( Figure 2G ) Cox analysis, it appeared that risk scores could be an independent prognostic factor compared to other clinical traits.

In the mRNA sequencing data of NPC from Fujian Cancer Hospital, the expressions of the nine immune-related hub genes were apparently different in the metastasis and non-metastasis groups ( Figure 3A ). Of the 9 genes, CD8B, APLNR, A2M and KIR3DL2 were upregulated in the non-metastasis group where the patients would have lower risk and gain better outcomes. In contrast, SEMA4F, PRDX2, RAC3, ULBP1 and TMSB15B were upregulated in the metastasis group where the patients would have higher risk and suffer worse outcomes ( Figure 3A ). To verify the predictive validity of nine hub genes for prognostic outcome, survival analysis demonstrated promising prognostic differentiation ( Figures 3B–J ).

The pathway in which the gene is enriched tends to indicate that the gene plays a role in that physiological process. Using the GSEA method, the high-risk group was mainly distributed into the E2F, G2M checkpoint, and MYC targets, which were closely related to interactions on angiogenesis, extracellular matrix remodeling, and tumor cell-endothelial cell interactions ( Figure 4A ). Correspondingly, the low-risk group was mainly distributed in the INF-γ, INF-α, and inflammatory responses, which were closely related to antitumor effect in anti-tumor immune response ( Figure 4B ).

Moreover, patient samples from high- and low-risk groups showed significant differences in scores across the ten carcinogenic pathways ( Figure 4C ). Patients in the high-risk group had higher oncogenic pathogenic activity, predicting that a higher risk of cancer progression was involved. And the CCP scores corroborated this finding, with patients in the high-risk group having high CCP scores, which suggested that the tumors had stronger proliferative activity ( Figure 4D ).

Here, we estimated how the immune microenvironment differed between patients in high- and low-risk groups in terms of immune scores and levels of immune cell infiltration. The patients from the low-risk group had higher immune scores but lower tumor purity ( Figures 5A, B ). Additionally, the compositions of the 29 immune-cell types were significantly different in the high- and low-risk groups ( Figure 5C ). In the low-risk group, almost all levels of immune cell infiltration were higher than in the high one, including B cells, CD8 T cells, dendritic cells (DC), macrophages ( Figure 5D ). Moreover, there were significant statistical differences in the immune checkpoint inhibitors (CTLA-4, HAVCR2, SIGLEC15, TIGIT, PD1 and LAG3) between the high- and low-risk groups ( Figure 5E ).

We were the first to evaluate GEP and TLS score, and showed that there were higher levels of immune cell receptors in low-risk patients ( Figures 6A, B ). Subsequently, the same results were observed in numerous immune-related indicators ( Figure 6C ). These results suggested that tumors stimulate more immune cell activation and strong ligand-receptor activation in patients in the low-risk group, laying the biological foundation for a positive response to this immunotherapy. As Figure 6D–G showed, patients in the low-risk group had a higher immune response in a cohort of patients with whether nasopharyngeal carcinoma or melanoma, or non-small cell lung cancer. It was evident that the patients of the high-risk group had less chance of benefiting from immunotherapy, which represented a worse prognosis when compared with the patients of the low-risk group ( Figure 6H ).

A cluster of 29 distinct cell types in GSE150430 cohort was defined by two-dimensional spatial visualization of UMAP analysis ( Figure 7A ). Cell lineages were distributed to each cluster by gene expression with reference to the human primary cell atlas data in TISCH. As a result, cells were annotated ( Figure 7B ). We targeted the most significantly differentially expressed genes in each cluster to better understand the species of cell fascicles ( Figure S3A ). In the identified cell subsets, the GSVA and ssGSEA algorithm was employed to calculate the performance of the nine hub genes at the single-cell level. Significantly higher risk scores were observed in malignant cells than in B cells and CD8 T cells ( Figures 7C , S3B ). The same conclusion can be drawn in the cellular localization map ( Figures S3C, D ). Moreover, the percentage of B cells and CD8 T cells in the low-risk samples was notablely higher than that of the high ones; however, the percentage of malignant cells in the high-risk samples was significantly higher than that of the low ones ( Figures 7D, E ). This was consistent with previous findings indicating that the high-risk scores predicted worse biological behavior.

Next, we carried out functional exploration. The major pathways enriched for differential genes between high- and low-risk groups were those related to intercellular adhesion and immune cell activation, suggesting that the response to distant metastasis and immune resistance differed between high and low-risk groups ( Figure 7F ). Also, active pathways were observed to vary in the high- and low-risk groups, like LT, TGFb, SEMA3, KIT, FGF and CD70 pathways being active in the high group while CALCR, CD40, and SPP1 pathways being vibrant in the low group ( Figures 7G , S3E ). In Figures 7H, I, the distinction of CD70 and SPP1 signaling pathways in high- and low-risk groups was more intuitive. Finally, the intracellular expression of nine hub genes is exhibited ( Figure 7J ). It can be seen that the expression of PRDX2, TMSB15B, ULBP1, and RAC3 was specifically increased in malignant cells, and the high expression of these genes coincides with a worse survival prognosis ( Figures 3E–H ).