The analyzed scRNA-seq data were deposited at the Gene Expression Omnibus database under accession code: GSE150430 (14). All samples were collected at Sun Yat-sen University Cancer Center. Written informed consent was obtained from all participants. Ethical approval was obtained from the Institutional Review Board of Sun Yat-sen University Cancer Center. The inclusion criteria for NPC samples we used were (1):. pathologically confirmed NPC (2); treatment-naive (3); EBV positive as confirmed using in situ hybridization of EBV encoded small RNAs in tumor tissue (4); available baseline EBV DNA data. Ten EBV-related NPC samples and one chronic nasopharyngitis sample were included in the final analysis. The study samples’ sequencing parameters and patient characteristics are summarized in Table S1 and Table S2 , respectively. EBV DNA Sero- was defined as undetectable EBV DNA in peripheral blood. EBV DNA Sero+ was defined as detectable EBV DNA in plasma with a titer of 10³ copies/mL or greater because of the detection limits of the plasma EBV DNA assays used in the study.

The preprocessed gene expression matrices for study samples were converted to a Seurat object using the Seurat package (version 4.1.0) in R (18). Low-quality cells with unique molecular identifiers < 200, expressed genes > 9000 or < 200, or mitochondrial genes > 20%, were removed. In the remaining cells, the gene expression dataset was normalized using the NormalizeData function and scaled with linear regression using the ScaleData function. In data scaling, the variables used to regress were “S.Score” and “G2M.Score” calculated by the CellCycleScoring function (19). The included features were the top 6000 highly variable genes selected by the FindVariableFeatures function.

The selected HVGs were summarized by principle component analysis (PCA) to reduce dimensionality, and the first 20 principle components were further projected using the Uniform Manifold Approximation and Projection (UMAP) (20). The number of principal components used was determined by an Elbow plot. We performed cell clustering using the FindClusters function with default parameters. The cell clusters were annotated as different major cell types based on DEGs and well-recognized marker genes. For sub-clustering analysis, we performed the second-round PCA reduction and UMAP projection separately on cells within each major cell type. The number of principal components and resolution used in each major subtype were dataset-dependent. The second-round dimension reduction and clustering revealed 28 distinct subpopulations in T/NK, B, and myeloid cells.

DEGs for each cluster were identified with the FindAllMarkers function using the MAST test. The differentially over-expressed genes in the specific group compared to the other group were identified with the FindMarkers function using the MAST test. We used the clusterProfiler package (version 4.0) for GO terms and KEGG pathways enrichment analysis. We also applied the Gene Set Variation Analysis (GSVA) implemented in the GSVA package (version 1.40.1) to estimate the pathway activities. Differences in pathway activities between groups were calculated using a linear model provided with the Limma package (version 3.48.3).

We assessed differences in transcriptional profiles across major cell types in EBV DNA Sero- and Sero+ NPCs using the Bhattacharyya distance method implemented in the distdimscr package (version 0.0.0.0.9) (21, 22). The distdimscr quantifies these differences by examining the distance between cell types in a high-dimensional space. Here, we only measured distances between cell types that had 300 or more cells present across all samples from each tumor type (here, EBV DNA Sero- vs. Sero+) in PCA. We embedded all cell types in a new PCA space based on highly variable genes and retained the top 10 PCs for subsequent analysis. We then sampled 300 cells from each tumor type for 100 replicates. The Bhattacharyya distance between cells sampled regardless of EBV DNA status was used to generate a background distribution for statistical comparison. The Wilcoxon rank-sum test was conducted to compare the EBV DNA Sero- vs. Sero+ cell types against the random selection for each cell type.

The large-scale chromosome CNVs in epithelial cells were inferred from raw gene expression matrices using the inferCNV package (version 1.10.1) (19). We performed inferCNV analysis with the following parameters: “denoise”, default Hidden Markov Model (HMM) settings, “subclusters” mode with “random_trees” method, and a value of 0.1 for “cutoff”. Epithelial cells extracted from the normal sample were used as reference data. To distinguish malignant and non-malignant epithelial cells in NPC samples, K-means clustering analysis was performed separately on epithelial cells derived from each of the eight tumor samples and the non-tumor sample according to CNV scores. Cells of the tumor samples in the cluster that predominantly contained cells from the normal sample were regarded as non-malignant cells. In contrast, the cells with high CNV scores in other clusters were identified as malignant cells. Consulting with the genome annotation data (GRCh38 version), each CNV was further converted into a p- or q-arm level change according to its position and annotated as either a gain or a loss. After data conversion, subclones containing identical arm-level CNVs were collapsed, and trees were restructured to represent subclonal CNV architecture. The UPhyloplot2 plotting algorithm was used for the visualization of intra-tumor evolutionary trees (23, 24).

Consensus Non-negative Matrix factorization (cNMF, version 1.4) implemented by Python was used to infer gene expression programs from the scRNA-Seq data (25). A total of 50 gene expression programs from the eight tumor samples with at least 50 epithelial cells. We next determine the underlying biological functions of each of the 50 gene expression programs using Gene Set Enrichment Analysis (GSEA) based on their top 50 rank genes. Hierarchical clustering based on pairwise correlation coefficients calculated from each possible GEP pair was further performed to identify potential co-expression modules. We identified six potential co-expressed gene modules. To determine the gene signature of each co-expression module, we combined the top 50 genes of each meta-gene expression program. We calculated the average weighted score of each gene according to the NMF load. Finally, the top 50 genes with the highest weighted scores were defined as the gene signature of the corresponding co-expression module.

For scRNA-seq data, we calculated the module/signature enrichment scores on a single-cell level using the AddModuleScore function provided by the Seurat package. For bulk RNA-seq data, we calculated module/signature enrichment scores for each sample using the Mann-Whitney U statistic provided by the UCell package (version 1.44.2). Details on gene signature datasets we used were summarized in the supplement ( Table S3 ). These gene data sets were collated based on reviewing relevant literature reports.

We assessed the associations of the malignant signatures and survival outcomes in a bulk RNA-seq cohort (GSE102349), which consisted of 113 NPC samples (26). Among 113 NPC samples, 88 samples with available survival data were used for subsequent survival analysis. We determined the optimal cutpoints for gene expression signatures using the maximally selected rank statistics from the ‘maxstat’ R package (version 0.7.25).

We inferred gene regulatory networks using the SCENIC package (version 1.2.4) implemented in R (27). The gene regulatory networks are inferred based on co-expression modules and TFs motif enrichment analysis from scRNA-seq data. We calculated and ranked the activities of TFs, as measured by Regulon Specificity Score (RSS), using the AUCell package (version 3.12) (28).

The pseudotime trajectories of malignant cells and CD8⁺ T cells were inferred with the Monocle 2 package (version 2.22.0). Only genes with high dispersion across cells were used in the trajectory analysis. Dimension reduction was performed with the reduceDimension function using the “DDRTree” method and a value of 2 for “max_components”. We ordered and visualized cells using the plot_cell_trajectory function. Genes that separated cells into two differentiated cell states were identified by the Branched Expression Analysis Modeling (BEAM) implemented by the BEAM function. Branch-dependent genes with a q-value less than 10^-4 were classified with hierarchical clustering using the plot_genes_branched_heatmap function. Genes that changed along with the pseudotime were identified with the differentialGeneTest functiom and visualized with plot_pseudotime_heatmap function. The GO terms enrichment analysis for the genes in each cluster was performed to characterize the biological functions.

The CytoTRACE (Cellular (Cyto) Trajectory Reconstruction Analysis using gene Counts and Expression) computational method, another robust trajectory analysis tool provided in the CytoTRACE package (version 0.3.3), was performed to validate the Monocle 2 inferred linear transition (29). CytoTRACE can infer the direction of differentiation from scRNA-seq data without any prior knowledge.

We assessed the stemness of malignant cells using the RNA expression-based stemness index (mRNAsi). The mRNAsi was proposed by Malta et al. using the One Class Linear Regression (OCLR) algorithm based on pluripotent stem cell samples from the Progenitor Cell Biology Consortium dataset (https://www.synapse.org) (30). The generated mRNAsi is a standardized value between 0 and 1; the closer it is to 1, the stronger the cell’s stemness. We applied the model to the scRNA-seq dataset to calculate the mRNAsi of each malignant cell.

We performed ligand-receptor analyses using the CellChat package (version 1.1.3) to compute cell-cell communication from single-cell transcriptomics data. CellChat has built ligand-receptor interactions for human in the CellChat database. The database considers known interactions between ligands, receptors, and cofactors. CellChat quantifies the cell-cell communication probability and identifies significant cell-cell interactions. CellChat can compare cell-cell communication patterns between two scRNA-seq datasets and identify significant signaling changes across different conditions.

All statistical analyses were performed using R (version 4.1.0) or Python (version 3.9.5). The differences between groups were analyzed by chi-square test, student t test, Wilcoxon test, Pearson correlation test, and Kruskal-Wallis test, when appropriate. All tests were two-sided and P < 0.05 was considered as statistical significance.

We obtained single-cell transcriptome sequencing data from ten NPC samples and one non-tumor nasopharyngeal tissue for bioinformatics analysis ( Figure 1A , Table S1, S2 ). The ten NPC samples comprised three EBV DNA Sero- and seven EBV DNA Sero+ NPCs. After data processing and quality control ( Figure S1 ), a total of 28,423 single cells were subsequently analyzed. We conducted major cell-type annotation based on classical gene markers ( Figure 1B ). Four major cell types were identified: epithelial cells, T/NK cells, B cells, and myeloid cells ( Figures 1C, D ). The cell composition in each sample, along with patient characteristics, are provided in Figure 1E . The distribution of these cell types differed significantly between the EBV DNA Sero- and Sero+ NPC groups ( Figure 1F ). The proportions of epithelial cells and myeloid cells were higher in the EBV DNA Sero- group while those of T/NK cells and B cells were higher in the EBV DNA Sero+ group. We then evaluated differences in transcriptional profiles between major cell types in EBV DNA Sero- and Sero+ groups using the Bhattacharyya distance ( Figure 1G ). Epithelial cells had the most considerable differences. The mean fold change between groups and random samples varied from 21.8-fold (epithelial cells) to 6.2-fold (T/NK cells). The above data revealed significant shifts in cellular compositions and transcriptomics between EBV DNA Sero- and Sero+ NPCs, suggesting a distinct TME context may be linked to EBV DNA seropositivity.

We next analyzed eight NPC samples that each sample contained at least 50 epithelial cells from the individual sample, including two EBV DNA Sero- and six EBV DNA Sero+ samples. We used large-scale CNVs to distinguish the malignant cells from the non-malignant cells and probe the clonal structure of each tumor ( Figures 2A, B ). We identified 5,726 malignant cells from eight NPC samples for further analyses, The inferred CNV profiles indicated that deletions in chromosomes 3p, 9p, 11q, 14q, and 16q and amplifications in chromosomes 12q were frequent in NPC ( Figure 2C ), which were consistent with previously reported whole-exome sequencing results (31, 32). The clonality analysis revealed the complexity of typical and atypical CNV events in NPC ( Figure 2D ). Typical CNV events dominate the chromosomal landscape. Interestingly, we also observed multiple subclonal typical and atypical CNV events in some samples. This result indicated that canonical CNV events do not constantly occur in a single CNV event but can be generated in an evolving genome. Several previous studies have reported that tumors can evolve with the development of unknown atypical CNV subclones, which may contribute to tumor progression (24, 33).

We next investigated the intertumor and intratumor heterogeneity of malignant NPC cells. We first used cNMF to identify a total of 50 gene expression programs that were preferentially co-expressed by subpopulations of malignant cells within each tumor across all samples ( Figure 3A ). Correlation clustering combined with GSEA was used to characterize these meta-programs into underlying co-expressed gene modules. We identified six co-expression modules ( Figure 3B ). We then analyzed the relationship between gene signatures of the modules and survival outcomes in bulk RNA-seq data of an NPC cohort (GSE102349, n = 88). Cell cycle signature, Epi-Diff 1 signature, and immune-related signature were significantly associated with progression-free survival ( Figure 3C ). The co-expressed gene modules were differentially expressed in EBV DNA Sero- and Sero+ NPC cells ( Figure 3D ). Specifically, NPC cells from EBV DNA Sero+ samples have higher cell cycle scores indicating higher proliferative capacity. Differential gene expression analysis showed that CDKN2A, ATF3, CXCL10, and interferon (IFN) stimulated genes (IFI6, IFIT3, IFI44L, IFITM1, etc.) were significantly upregulated in EBV DNA Sero+ NPC cells ( Figure S2 ), suggesting more active EBV viral activity. We further used GSVA to reveal cancer hallmarks pathways. Notably, the results exhibited more cancer hallmarks pathways enriched in EBV DNA Sero+ samples, including IFN α/γ response, hypoxia, p53, and TNF-α signaling via NFκ-b pathways ( Figure 3E ). Nevertheless, a set of transcription factors (TFs) was upregulated in EBV DNA Sero+ tumor cells, evidenced by SCENIC analysis ( Figure 3F ). Among them, the expression of genes regulated by ATF3, JUN, IRF7, MYC, and JUNB showed highest regulon activity ( Figure 3G ).

We next explored the dynamic transitions of NPC cells by inferring the developmental trajectories using Monocle2 ( Figure 4A ). The pseudotime analysis showed that EBV-associated NPC cells experienced three potential differentiation states. The distribution of cell states in tumor cells with different EBV DNA statuses varied significantly, that the EBV DNA Sero- tumor cells predominated at the beginning of the trajectory path (state 1), whereas most EBV DNA Sero+ tumor cells were at a terminal state (state 3) ( Figures 4B, C ). In addition, we used the CytoTRACE algorithm and OCLR algorithm to quantify the developmental potential and stemness index of tumor cells, respectively ( Figures 4D, E ). EBV DNA Sero+ NPC cells showed low-differentiation potential and higher stemness index than EBV DNA Sero- NPC cells ( Figure 4F ), which indicated the high aggressiveness of EBV DNA Sero+ NPC. Next, we visualized modules of genes that shared similar lineage-dependent expression patterns using the radiation heatmap ( Figure 4G ). The result indicated that cell states 2 and 3 had more such genes related to ribonucleoprotein complex, mitochondrial translation, and translational termination pathways. Besides, neutrophil and T cell-mediated immunity pathways were upregulated in the differentiation trajectory of cell state 3. Furthermore, the SCENIC analysis revealed that STAT1/2, ATF4, JUNB, IRF1/7, and other TFs related to cell signal transduction and IFN response pathways were highly activated in cell state 2, while cell state 3 had upregulated TFs of THAP11, MYC, ATF3, and JUN that related to cellular stemness, stress response, and cellular senescence pathways ( Figure 4H ).

The above data revealed significant alterations in transcriptomics and differentiation between EBV DNA Sero- and Sero+ NPC malignant cells, suggesting EBV seropositivity is strongly associated with NPC malignant cell development.

A total of 10,039 T/NK cells were analyzed and grouped into 16 clusters. Based on the DEGs ( Figure S3A ) and function-related markers ( Figure S3B ), the 16 clusters were classified into four CD8⁺ T cell subtypes (Cluster 12/13: proliferating T cells; Cluster 1/14: effector memory T cells [T_em cells]; Cluster 7/10: granzyme K-positive [GZMK⁺] exhausted T cells [T_ex cells]; Cluster 6: IFN-stimulated genes-positive [ISG⁺] T cells), five CD4⁺ T cell subtypes (Cluster 3/5: naive T cells [T_n cells]; Cluster 4: ISG⁺ T helper cells [T_h cells]; Cluster 9: proliferating regulatory T cells [T_reg cells]; Cluster 2: resting T_reg cells; Cluster 8/11: suppressive T_reg cells), and two natural-killer T (NKT) cell subtypes (Cluster 15: cytotoxic NKT cells; Cluster 16: resting NKT cells) ( Figure 5A ). The proportion of T cell subtypes differed among patients ( Figure S3C ) and patients with different EBV seropositivity status ( Figure 5B ).

The overall pathway analysis in CD8⁺ T cells with GSVA revealed that pathways of cytotoxicity, NK cell functions, antigen processing, and Toll-like receptors were upregulated in EBV DNA Sero+ CD8⁺ T cells, which is consistent with activated immune defendant machinery to fight against EBV activation ( Figure 5C ). Remarkably, cytotoxic and IFN-γ signatures were elevated in EBV DNA Sero+ samples, while naïve, effector memory, and exhaustion signatures were elevated in EBV DNA Sero- samples ( Figure 5D ). According to the cytotoxic and exhausted signatures ( Figure 5E ), CD8⁺ T cells were further grouped into four functional subtypes (Cyt^LowExh^Low, Cyt^LowExh^High, Cyt^HighExh^Low, and Cyt^HighExh^High). EBV DNA Sero+ group had more Cyt^LowExh^High CD8⁺ T cells and fewer Cyt^HighExh^Low CD8⁺ T cells ( Figure 5F ), suggesting immunosuppressive TME associated with EBV DNA seropositivity status. In addition, certain classical inhibitory markers, such as CTLA4, LAG3, and HAVCR2 were upregulated EBV DNA Sero- CD8⁺ T cells, whereas other inhibitory markers, such as PDCD1, TIGHT, ENTPD1, and TOX were upregulated in EBV DNA Sero- CD8⁺ T cells ( Figure 5G ). These indicated different immunoinhibitory mechanisms employed by malignant cells depending on EBV DNA seropositivity status and may be helpful to guide the future development of immunotherapy strategy based on EBV seropositivity status.

In comparing cell composition, ISG⁺ CD8⁺ T and CD8⁺ T_em levels were elevated in EBV DNA Sero+ samples, whereas GZMK⁺ CD8⁺ T_ex and proliferating CD8⁺ T levels were elevated in EBV DNA Sero- samples ( Figure S3D ). We further explored the dynamic differentiation lineages and cell transitions of these CD8⁺ T cell subtypes by inferring the pseudotime trajectories ( Figure 5H ). This analysis showed a developmental trajectory that mainly began with the ISG⁺ CD8⁺ T cells, bifurcated into either the CD8⁺ T_em cells or the GZMK⁺ CD8⁺ T_ex, and ended with proliferating CD8⁺ T cells ( Figure 5I ). CD8⁺ T cell subtypes showed relatively similar distribution patterns along with the pseudotime in samples with different EBV DNA statuses, except that ISG⁺ CD8⁺ T cells exhibited a peaked distribution at the very beginning ( Figure 5J ). Overall, EBV DNA Sero+ CD8⁺ T cells presented a lower pseudotime score compared with EBV DNA Sero- CD8⁺ T cells ( Figure 5K ). These results indicated that the cytotoxic T-lymphocyte response may be triggered early in EBV DNA Sero+ NPC. We next examined the transcriptional changes related to the differentiation process and found that the CD8⁺ T cell differentiation lineages could be classified into four major phases with distinct enriched GO pathways ( Figure 5L ).

We next examined the cell composition and transcriptome heterogeneity of CD4⁺ T cells driven by EBV DNA status. Overall, the EBV DNA Sero+ group had a higher proportion of T_reg cells (53.3 vs. 35.0%) than the EBV DNA Sero- samples. The EBV DNA Sero+ group had more ISG⁺ Th cells (17.3% vs. 3.4%) and resting T_reg cells (22.4% vs. 9.8%) but fewer proliferating T_reg cells (8.7% vs. 45.1%) than that in EBV DNA Sero- group. A recent report suggested that ISG⁺ Th cells might represent a fraction of T_reg cells that responded to type I IFNs in TME during activation (34). The total IL2R, inhibitory, and costimulatory signatures were all elevated in the EBV DNA Sero- CD4⁺ T cells ( Figure 5M ). The result suggested that T_reg cells were more activated in EBV DNA Sero- NPC and had a more immunosuppressive context than the EBV DNA Sero+ NPC. GSEA showed that the hypoxia, IFN-α response, and KRAS signaling were significantly upregulated in the EBV DNA Sero+ samples ( Figure 5N ), whereas E2F targets, epithelial-mesenchymal transition, and G2M checkpoint were upregulated in the EBV DNA Sero- samples ( Figure 5O ).

In total, 8592 B cells were identified and grouped into seven subsets according to the DEGs and classic markers ( Figures S4A, B , Figure 6A ), including memory B cells (cluster 1), follicular B cells (cluster 2), IFN-induced B Cells (cluster 3), plasmablasts (cluster 4/7), FCRL4+ B cells (cluster 5), proliferating B cells (cluster 6), and germinal center B cells (cluster 8). The proportion of B cell subtypes showed variation across patients ( Figure S4C ) and between groups with different EBV DNA statuses ( Figures 6B, C ). The proportions of FCRL4+ B cells and memory B cells were higher in the EBV DNA Sero+ samples, whereas that of the proliferating B cells and germinal center B cells were higher in the EBV DNA Sero- samples, indicating a strong immunoregulatory potential in EBV DNA Sero+ NPC. GO enrichment analysis showed that the antiviral response, IFN response, antigen processing and presentation, and ribosome assembly were activated in the EBV DNA Sero+ group ( Figure 6D ). KEGG pathway enrichment analysis revealed strong activation of T-cell differentiation, antigen processing and presentation, and EBV infection among EBV DNA Sero+ samples ( Figure 6E ). SCENIC analysis revealed that TFs of JUN/JUNB, FOS/FOSB, IRF7/9, ELF1, STAT1 were upregulated in EBV DNA Sero+ samples ( Figure 6F ). These observations indicate that B cells in EBV DNA Sero+ NPC exhibit an inflammation-dominant antiviral immune response.

The re-clustering of 2,902 myeloid cells revealed eight populations ( Figure 6G ), including four subsets of macrophages and four subsets of dendritic cells (DCs) ( Figures S5A, B ). We could not distinguish M1 and M2 macrophages based on curated marker genes, however, we examined the expression of M1 and M2 polarization signatures of macrophages subsets ( Figure S5C ). Macro1 was a classic activated M1 macrophage with high activation of IL1B, IL8, LYZ, and IFI30 genes, characterized by intermediate M1 polarization score and low M2 polarization score. Macro2 was a proliferating macrophage marked with high expression of cell cycle-related genes and low M1 and M2 polarization scores. Macro3 showed high M1 and M2 polarization signals, consistent with previous reports on M1/M2 couple-activated macrophages in NPC (16, 35). Macro4 showed a high M2 polarization signal, which may be a selectively activated M2 macrophage. DC1 was plasmacytoid DCs (pDCs) characterized by increased expression of CLEC4 and IGJ genes. DC2 was a mature DC that showed high expression levels of LAMP3 and CCR7, yielding a high potential for migration, activation, and maturation ( Figure S5D ). DC3 represented the classical DCs that were characterized by overexpression of CLEC9A. DC4 showed increased expression of CD1A, CD207, and ID2, suggesting that these cells were Langerhans cells (LCs). The myeloid cell compositions differed among patients ( Figure S5E ) and between EBV DNA Sero- and EBV DNA Sero+ samples ( Figures 6H, I ). Notably, Macro3 (99.4%) and Macro4 (98.9%) were almost all from EBV DNA Sero+ samples, suggesting predominant macrophage activation and M2 polarization. In the meanwhile, the M1 and M2 polarization scores were significantly upgraded in EBV DNA Sero+ NPC ( Figure 6J ). GSVA suggested different patterns of enriched pathways in the myeloid compartment between different EBV DNA statuses. Several pathways involved in myeloid innate immunity were upregulated in EBV DNA Sero+ myeloid cells, including angiogenesis, Fc receptor signaling, cell cycle and apoptosis, cytokine signaling, and IFN signaling pathway ( Figure 6K ). In contrast, T cell activation and immune checkpoint signaling, cell migration adhesion, and antigen presentation signaling pathways were activated in EBV DNA Sero- myeloid cells. SCENIC analysis revealed that TFs of JUN/JUNB, FOS/FOSB/FOSL2, IRF1/2/4/7/8, and STAT1/2 were upregulated in EBV DNA Sero+ samples ( Figure 6L ). These observations indicated that the ability to shape the immune environment was enhanced in EBV DNA Sero+ NPC.

Previous studies have reported that IFN responses were significantly upregulated in the TME of NPC, especially in recurrent NPC (16, 35). Our DEG analysis suggested that genes related to IFN responses were upregulated in most cell types of EBV seropositive NPC ( Figure 7A ). We further calculated type I IFN response and type II IFN response scores for each major cell type using well-defined gene markers. The results showed that type I IFN response and type II IFN response scores were significantly upgraded in almost all cell types from EBV DNA Sero+ samples ( Figure 7B ). These observations suggest that the IFN responses are globally upregulated in the TME of EBV DNA Sero+ NPC, exhibiting an antiviral immune response-dominated context.

We next used the CellChat package (version 1.1.3) to explore the cellular communication network in NPC immune microenvironment. We detected complex communications between immune cells ( Figure 7C ). Memory B cells, macrophages, and DCs showed the highest intercellular interactions. Overall, the number and intensity of cell-cell interactions in EBV DNA Sero+ NPC were higher than those in EBV DNA Sero- NPC ( Figures 7D, E ). Macrophages represented the cell type with the most incoming and outgoing signalings, regardless of EBV DNA status ( Figure 7F ). Surprisingly, we did not observe incoming and outgoing signalings in memory B cells from EBV DNA Sero- samples, while memory B cells from EBV DNA Sero+ samples showed strong incoming signalings, indicating a stronger immunoregulatory potential. We further identified and ranked the essential signaling pathways ( Figures 7G, H ). Several signaling pathways that induce inflammatory responses and are involved in immune regulation were significantly upregulated in EBV DNA Sero+ NPC, such as SPP1, BTLA, TNFSF12, and IL6. SPP1, as a cytokine, plays a vital role in the type-I immune response by enhancing the production of IFN-γ and IL12 and reducing the production of IL10.