The microarray data of 116 NPC samples and 17 adjacent healthy nasopharynx samples used in the present study were downloaded from Gene Expression Omnibus (GEO) with accession code GSE12452 (Sengupta et al., 2006), GSE13597 (Bose et al., 2009), GSE64634 (Bo et al., 2015) and GSE103611 (Tang et al., 2018). The expression data were measured by different platforms and protocols and have strong batch effects that would affect downstream analysis. We therefore used sva package (Leek et al., 2012) to remove the batch effect, and then merged the four datasets by removing the genes that were not present in any of the four studies. Afterward, a quantile normalization was performed to avoid the influence of outlier data points. Principal component analysis (PCA) was performed to visualize at low dimensions whether the batch effect was removed.

For clinical information, we used the grades of lymph node spread (N) and diagnosis of overt metastasis (M) to split the tumors into three progression stages: Stage 1 includes the patients with N0M0 or N1M0; Stage 2 includes the patients with N2M0, N3M0, N0M1 or N1M1; and Stage 3 includes N2M1 or N3M1. The overall survival data from the dataset GSE64634 (Bo et al., 2015) were used to analyze survival difference among the three patient groups using survival package (Therneau, 2020), and was visualized using Kaplan-Meier plot (K-M plot).

Given the stage variable x = 1, 2, 3, four denoting normal (Stage 0) and progression Stage 1, 2 and 3, respectively, and expression of each gene y, we estimated the trend of expression in different stages (x,y) in each sample using least squares fit ∑ [ y − ( β 1 x + β 0 ) ] 2 . The null hypothesis is β 1 = 0 denoting no trend, and we performed a permutation test by shuffling the group variable 10,000 times at random to generate a Gaussian distribution centered at zero. An empirical p value was calculated to denote how likely the fitted β 1 is sampled from the distribution generated by the random shuffling. We selected the genes with p < 0.05 as the set associated with NPC progression. And the gene was sorted by the corresponding β 1 value and the ranking list was input into GSEA (Subramanian et al., 2005) under “pre-rank” mode to query enriched cancer hallmarks curated by the Molecular Signatures Database (MSigDB) (Liberzon et al., 2015).

We downloaded the raw data of single-cell RNA sequencing of NPC samples at different progression stages from one previous study (Chen et al., 2020) under accession code CNP0000428 from the China National GeneBank Sequence Archive. The raw data were aligned and processed into count tables using CellRanger (Zheng et al., 2017) developed by 10x Genomics Inc. Further downstream analysis was performed using Seurat R package (Satija et al., 2015; Butler et al., 2018).

We used ConsensusTME (Jiménez-Sánchez et al., 2019) to computationally deconvolve the expression profiles of the 133 samples into 18 immunological components and one prognostic metric termed immunoscore. The parameters used in the ConsensusTME are “cancer = HNSC” and “statMethod = ssgsea”.

The formalin-fixed paraffin-embedded (FFPE) tissue slides were sectioned at 4 μm thick from the tumor specimens resected at Department of Otolaryngology, People’s Hospital of Haining, and were mounted on coated glass slides for immunohistochemistry (IHC) staining. The slides were stained by Ultra60-1,600 automatic IHC stainer (ZSGB-BIO, Beijing, China). Anti-CD20 antibodies (L26, dilution 1:2000, #ZM-0039, ZSGB-BIO, Beijing, China) were used for staining. Digital pathological images were scanned using digital whole slide scanning (HS6, 190,,011-6, SOPTOP, Ningbo, China), and were analyzed using the software Aperio-ImageScope (Leica, Version: 12.4.3). We obtained the digital images of the stained tissue sections at ×10 magnification.

We used Oncomine (Rhodes et al., 2007) to test whether inflammation-associated signature genes are significantly differentially expressed in multiple cancer transcriptomic datasets. Given one input gene symbol, Oncomine performed differential expression analyses using 715 datasets from the datasets of over 20 cancer types, and returned the number of analyses within which the input gene is significantly expressed between cancer and normal samples. In this study, the significant expression was defined as p value <0.001, log fold change >3 and p-value ascending ranking list of top 5% out of all the genes.

To identify genes associated with NPC progression, we collected the microarray data of 133 NPC at different pathological grades and normal samples from four previous studies: GSE12452 (Sengupta et al., 2006), GSE13597 (Bose et al., 2009), GSE64634 (Bo et al., 2015) and GSE103611 (Tang et al., 2018). Using the measures of lymph node spread (N) and overt metastasis (M) in TNM system, we split the 133 samples into four progression stages (Figure 1A). To rationalize our strategy in progression staging, we used the strategy to split a cohort of 112 NPC patients with survival data available from GSE64634 (Bo et al., 2015) using the same criteria, which has shown to have different overall survivals (p = 0.21 between Stage 1 and 2, p = 0.0098 between Stage 2 and 3, and p = 0.001 between Stage 1 and 3, log-rank test, Figure 1B). Unexpectedly, we found that splitting the NPC patients using solely tumor size did not yield significant difference in overall survival (p = 0.38, log-rank test, Supplementary Figure 1A), and some patients with T4-tumors exhibit better survival than those with T1-tumors.

To enhance statistical power, we collected a large dataset of NPC progression by integrating the microarray expression data from the four studies and removing the batch effect generated by different sequencing protocols. Using Principal Component Analysis (PCA) to project all the 133 samples into a 2D dimensions, we showed a well-mixed distribution among different batches of expression data (Supplementary Figure 1B). In addition, we found that the expression profiles of NPC and normal nasopharynx samples are remarkably different, while the expression of NPCs at different stages are not well-separated (Supplementary Figure 1C). To further investigate the expression changes between different stages, we developed a multiscale analytic framework consisting of four modules (Supplementary Figure 1D): 1) trend analysis for comparison of different patients at cohort scale that captures key progression modulators with increasing average expression at sequential order of progression stages, 2) single-cell analysis for comparison of different molecules at single-cell scale which identifies crucial cell-cell communications via chemokine-receptor interactions, 3) bulk-deconvolution for comparison of different immune cell compositions at patient scale that consolidate significant changes of immune cell fractions observed at single-cell scale during NPC progression, and 4) pan-cancer analysis for comparison of different cohorts at population scale that investigate common aberrant gene expressions in multiple types of cancers.

To identify the expression patterns associated with NPC progression, we developed an empirical statistical method to test whether the expression of one gene exhibits a linear trend at progressive stage order (see Methods). Using the trend statistic beta and its empirical significance, we prioritized 12,746 genes based on their increasing/decreasing trends during NPC stage progression, and identified 2,698 genes with significant trends (p < 0.05). To demonstrate the rationale of the trend analysis, we showed the expression of MKI67, a prominent histological biomarker of cellular proliferation, in different stages and observed a remarkably increasing trend from Stage 0, a healthy control, to Stage 3, the final stage with distant metastasis (Figure 2A). Visualizing the trend statistics of all the genes in Figure 2B, we discovered that the genes with top increasing trends are CXCL11, MMP1, CXCL10 and CXCL9, while the genes with the top decreasing trends are MSMB, LTF, SCGB1A1 and GSTA2. Notably, CXCL10, a chemokine responsible for chemoattraction, has been demonstrated in an independent large-scale study as a significant prognostic biomarker of distant metastasis in locoregional advanced NPC patients (Tang et al., 2018).

To uncover biological insights of the prioritized genes by the trend analysis, we performed the Gene Set Enrichment Analysis (GSEA) using the trend statistic as a pre-rank metric. Remarkably, the GSEA revealed 24 significantly enriched cancer hallmarks (FDR <0.05, Figure 2C), among which the top functional hallmarks are mainly from four categories of biological processes: cellular proliferation (E2F targets, MYC targets v1 and v2, G2M checkpoint, Mitotic spindle, dark green bars in Figures 2C, D), immunity (interferon alpha, gamma responses, IL6-JAK-STAT3 signaling, and inflammatory response, yellow bars in Figures 2C, D), development (epithelial mesenchymal transition and angiogenesis, red bars in Figures 2C,D), and signaling (MTORC1 signaling and TNFA signaling via NFkB, light green bars in Figures 2C,D). Strikingly, epithelial mesenchymal transition (EMT) as a predominant phenotype of cancer metastasis, were found to be significantly associated with NPC progression (FDR <0.0001) with MMP1 as the leading gene in the EMT gene set. In addition, another notable enriched hallmark is interferon gamma response (FDR <0.0001) with CXCL9/10/11 as the leading genes in the set. These findings suggest that beside accelerated cell proliferation, NPC progression involved the activation of EMT program in which malignant nasopharyngeal epithelial cells undergo a morphological transition into a mesenchymal-like phenotype with strong invasive and migratory potentials, and the activation of immune response to recruit lymphocytes by upregulating CXCL9/10/11.

Evidenced by the trend analysis (Figure 2B), the functional enrichment (Figures 2C, D) and previous prognostic study (Tang et al., 2018), we inferred that the chemokine CXCL10 might play a crucial role in tumor microenvironment of NPC progression. However, its biological role and dynamics in NPC progression remain largely unknown. To first investigate the cellular origin of CXCL10 in NPC progression, we collected single-cell RNA sequencing (scRNA-seq) data from a recent study with NPC patients at different stages (Chen et al., 2020). From the scRNA-seq dataset, we selected and integrated four samples from normal nasopharynx and NPC tumors at Stage 1, 2 and 3, respectively. The integrated expression profiles of 20,841 cells were projected into a 2D space (Figure 3A) using UMAP by Seurat package, and were further categorized into 14 cell types based on regular cell markers (Figure 3B). Similar to the microenvironment of other cancer types (Puram et al., 2017; Zhang et al., 2020), the NPC microenvironment consists primarily of T cells, B cells and myeloid cells (Figure 3A).

Highlighting the cells with highly expressed CXCL10 and its binding receptor CXCR3, we uncovered that CXCL10 is predominantly expressed by macrophages and its receptor CXCR3 is mainly expressed by T cells, indicating macrophages attract T cells through CXCL10-CXCR3 communication (Figure 3C). Notably, the cell proportion of CXCL10 + macrophages exhibit a significantly increase (z-test, Supplementary Figure 2A) from Stage 0 (normal) toward Stage 3 (metastasis) during NPC progression (Figure 3D), which is consistent with the trend identified by the trend analysis (Figure 2B). However, the CXCL10 + macrophages do not keep recruiting CXCR3+ effector T cells: the proportion of CXCR3+ effector T cells dramatically declines from 5.14% at Stage 2 to 0.98% at Stage 3 (z-test, Figure 3D and Supplementary Figure 2A). In contrast, the cell proportion of CXCR3+ regulatory T cells (Treg), which play a suppressive role in immune response, shows a significantly ascending trend from 0.84% at Stage 2 toward 1.97% at Stage 3 (z-test, Figure 3D and Supplementary Figure 2A). These results imply that macrophages recruit CXCR3+ effector T cells through CXCL10 expression at early stage of NPC progression, but at late stage, the NPC microenvironment might undergo a reprogramming from an active into a suppressive state, evidenced by the decline of CXCR3+ effector T cells and the increase of CXCR3+ regulatory T cells. In addition, we observed significant depletion of B cells during NPC progression (z-test, Figure 3E and Supplementary Figure 2B). Especially, memory B cells dramatically drop from 25.11% at Stage 0–0.72% at Stage 3, which suggests that like EBV-associated B-cell lymphomas, the B cells in NPC are vulnerable under EBV infection.

To further validate the dynamic changes of NPC microenvironment during progression discovered by the scRNA-seq analysis, we computationally estimated 22 immune cell proportions in the bulk microarray data of 133 NPC and normal samples (Figure 4A) using ConsensusTME (Jiménez-Sánchez et al., 2019). One of the apparent changes among the 22 immune cells is B cell depletion during the progression (Figure 4B), which is consistent with the observation of scRNA-seq analysis (Figure 3E). Estimating the cytotoxic T lymphocyte (CTL) level of each sample using the average expression of five marker genes (CD8A, CD8B, GZMA, GZMB and PRF1) (Rooney et al., 2015), we showed that the CTL levels boost up at the early stage but slightly drop at the late stage during NPC progression (Figure 4C). Similarly, the proportion of anti-tumor M1 macrophages display a trend of “increase first, decrease later” (Figure 4D). Consistent with the scRNA-seq analysis (Figure 3D), we also observed an increasing trend of Treg cell proportion in the deconvolution analysis (Figure 4E). All the evidence of the analyses pointed to a unique trend of tumor microenvironment during NPC progression: from active immune surveillance at the early stage to immunosuppressive’ at the late stage.

Enlightened by the consistent trends of B-cell depletion observed in the scRNA-seq and bulk deconvolution analyses, we conducted hematoxylin and eosin (HE) staining (Figure 5A) and anti-CD20 immunohistochemistry staining (Figure 5B) of three NPC samples at Stage 1, 2 and 3, and one adjacent normal sample collected by the Department of Otolaryngology, People’s Hospital of Haining, Zhejiang, China. Consistently, we observed a decreasing trend of the number of CD20-positive cells from Stage 1 to Stage 3 NPC samples (Figure 5C), which consolidates the computational results of B-cell depletion during NPC progression.

Commonly, the chemokines CXCL9, CXCL10 and CXCL11 are thought to be induced by interferon (IFN)-γ against inflammation (Metzemaekers et al., 2017). Using Oncomine (Rhodes et al., 2007), we investigated what types of cancers are associated with the upregulation of the three chemokines observed in NPC progression. Consistent with the observation in NPC, we found that the three chemokines are upregulated in 11 transcriptomic studies of head and neck cancer (Supplementary Figure 3A). Interestingly, the cancers in which the three chemokines are upregulated are significantly enriched with viral or bacterial infections (p < 0.05, z-test, Supplementary Figure 3A). Within the 10 cancers with the chemokine upregulation, six cancers (lymphoma, head and neck cancer, liver cancer, colorectal cancer, cervical cancer and bladder cancer) are associated with viral or bacterial infection. In contrast, only one cancer (gastric cancer) out of the other 10 cancers without chemokine upregulation is associated with bacterial infection. This suggests that the upregulation of CXCL9, CXCL10 and CXCL11 is a common event in the cancers with viral or bacterial association.