Abstract
Nasopharyngeal carcinoma (NPC) is an aggressive malignancy with extremely skewed ethnic and geographic distributions. Increasing evidence indicates that targeting the tumor microenvironment (TME) represents a promising therapeutic approach in NPC, highlighting an urgent need to deepen the understanding of the complex NPC TME. Here, we generated single-cell transcriptome profiles for 7581 malignant cells and 40,285 immune cells from fifteen primary NPC tumors and one normal sample. We revealed malignant signatures capturing intratumoral transcriptional heterogeneity and predicting aggressiveness of malignant cells. Diverse immune cell subtypes were identified, including novel subtypes such as CLEC9A+ dendritic cells (DCs). We further revealed transcriptional regulators underlying immune cell diversity, and cell–cell interaction analyses highlighted promising immunotherapeutic targets in NPC. Moreover, we established the immune subtype-specific signatures, and demonstrated that the signatures of macrophages, plasmacytoid dendritic cells (pDCs), CLEC9A+ DCs, natural killer (NK) cells, and plasma cells were significantly associated with improved survival outcomes in NPC. Taken together, our findings represent a unique resource providing in-depth insights into the cellular heterogeneity of NPC TME and highlight potential biomarkers for anticancer treatment and risk stratification, laying a new foundation for precision therapies in NPC.
Introduction
Nasopharyngeal carcinoma (NPC) is a unique subtype of head and neck cancers with an extremely unbalanced endemic distribution; it is highly prevalent in East and Southeast Asia, especially southern China, and North and East Africa.1 Etiologically, NPC is highly correlated with Epstein–Barr virus (EBV) infection and is pathologically characterized by heavy infiltration of immune cells around and within tumor lesions,2,3 suggesting the existence of a remarkably complex tumor microenvironment (TME) in NPC.
These specific features of the TME indicate the potential benefits of immunotherapy in NPC. Indeed, the inhibition of the immune checkpoint axis programmed cell death protein 1 (PD-1)/PD-ligand 1 (PD-L1) has achieved clinical benefit in NPC patients; nevertheless, the response rate to anti-PD-1 therapies in NPC patients is only 20%–30%.4,5,6 This highlights an unmet clinical need to deepen the understanding of the TME in NPC to identify therapeutic targets and reliable biomarkers for risk stratification. However, current strategies for genomic/transcriptomic analyses of NPC are primarily based on bulk samples from these often small-sized tumors,7,8,9,10 and these approaches therefore lack the resolution and accuracy to depict the complex heterogeneity of the TME.
Recent advances in single-cell RNA sequencing (scRNA-seq) are opening a new way for evaluating transcriptional features with cellular resolution in various types of cancers.11,12,13,14 This study represents a comprehensive and unique resource revealing the cellular heterogeneity of NPC TME at single-cell resolution. We uncovered malignant signatures capturing intratumoral heterogeneity. We also identified transcriptional regulators underlying immune cell diversity and established immune subtype-specific signatures associated with prognosis, thereby providing in-depth insights into the multicellular ecosystem of NPC.
Results
Single-cell expression atlas and cell typing of the TME in NPC
We generated scRNA-seq profiles for primary tumors from 15 patients with treatment-naïve NPC; normal nasopharyngeal epithelial tissue from one patient with chronic nasopharyngitis was also collected (Fig. 1a; Supplementary information, Fig. S1 and Table S1). The fresh biopsies, which were collected during endoscopy, were rapidly digested into single-cell suspensions and analyzed using droplet-based single-cell transcriptome profiles (10× Genomics Chromium system). We also obtained whole-exome sequencing (WES) and bulk RNA-seq data from 12 of the 15 NPC primary tumors (Supplementary information, Table S1).
Fig. 1: Dissection of the tumor microenvironment in NPC with scRNA-seq.
a Workflow diagram showing the collection and processing of fresh biopsy samples from 15 primary NPC tumors and one normal tissue for scRNA-seq. b t-SNE plots of cells from the 16 samples profiled in this study, with each cell color coded to indicate the associated cell types. The panels on the right show the expression of curated gene sets in the cell types defined in the left panel. c Chromosomal landscape of inferred large-scale CNVs distinguishing malignant epithelial cells from non-malignant epithelial cells. The P03 tumor is shown with individual cells (y-axis) and chromosomal regions (x-axis). Amplifications (red) or deletions (blue) were inferred by averaging expression over 100-gene stretches on the indicated chromosomes. Inferred CNVs are concordant with the calls from WES (bottom). The inferred CNV pattern of the normal epithelial cells from N01 is also shown (top). d t-SNE plots showing the subclusters of malignant cells, myeloid cells, T cells, and B cells. e Data of the 13 malignant cell subclusters and the 23 immune cell subclusters of 47,866 cells from 16 samples (from left to right): the fraction of cells originating from each patient, the number of cells, and box plots of the number of UMIs and genes (with the box plot center, box, whiskers, and points corresponding to the median, interquartile range, 1.5× interquartile range, and outliers, respectively).
After quality filtering and doublet removal (see Materials and methods), we profiled a total of 48,584 single cells from the NPC tumor samples and normal tissue sample (Supplementary information, Table S2). Following gene expression normalization, we applied principal component analysis to evaluate variably expressed genes and subsequently used a graph-based clustering method15,16 to classify the cells into coherent transcriptional clusters. Then, we annotated the cell clusters by the average expression of curated gene sets and identified immune cells (i.e., myeloid, T/natural killer (NK), and B cells) and cancer-associated fibroblasts (CAFs) in addition to epithelial cells (including malignant and non-malignant cells) (Fig. 1b). No pronounced effects of cell dissociation on gene expression17 were noted (Supplementary information, Fig. S2a, b), and the bulk RNA-seq data correlated well with the scRNA-seq data (r = 0.71; Supplementary information, Fig. S2c).
We then distinguished malignant from non-malignant cells within the epithelial clusters in two steps. First, as NPC arises from the nasopharyngeal epithelium, we extracted all the epithelial cells based on the initial clustering and cell-type identification described. Second, we inferred large-scale chromosomal copy-number variations (CNVs) in each epithelial cell.18,19,20 These inferred CNVs, which were consistent with the WES data (Fig. 1c), were used to separate the malignant cells from the non-malignant cells with normal karyotypes. We then performed hierarchical clustering of the CNV profiles of the epithelial cells from the normal sample and each of the 15 tumor samples separately21; the non-malignant cells of each tumor sample were identified as those in the cluster that predominantly contained cells from the normal sample, while the cells with deletions and amplifications of entire chromosomes in the other cluster were identified as malignant cells. Overall, 7581 malignant cells from 11 tumor samples containing at least 50 malignant cells were identified and retained for further analyses.
Then, we performed normalization, dimensionality reduction and clustering within each major cell type (i.e., malignant, myeloid, T/NK, and B cells) to identify subclusters (see Materials and methods); CAFs were removed from analyses as fewer than 100 cells were detected. This approach revealed a complex TME, containing 13 malignant cell subclusters and 23 different immune cell subclusters (Fig. 1d, e); subsequently we identified signature genes for each of the immune subclusters, successfully developing 17 immune subtype-specific signatures (Supplementary information, Table S3). The high abundance of major immune cell subtypes in NPC was then confirmed by multiplex immunofluorescence analysis (Supplementary information, Fig. S3), especially for the T/NK cells, and the mean percentages of CD4+ T, CD8+ T, and CD57+ NK cells in all cells were 19%, 15% and 2%, respectively. Furthermore, we revealed the disparities in distributions of these cell subtypes within tumors (Supplementary information, Fig. S3c, d). Generally, the immune cells were more enriched in the tumor stroma compared to the tumor epithelial nests, in particular for the CD20+ B cells, CD57+ NK cells and IL3RA+ plasmacytoid dendritic cells (pDCs). Notably, the infiltration of CD8+ T cells in tumor epithelial nests was relatively higher than that of other immune cell types, suggesting spatial heterogeneity of different immune cell subtypes within the tumors. It should be noted that few CAFs and endothelial cells were detected in the scRNA-seq, while the mean percentages of these cell types detected by immunofluorescence staining were 4% and 3%, respectively (Supplementary information, Fig. S3b). The differences reflect the well-known disparities in single-cell dissociation efficiency that influence the recovery of individual cell types,22 especially for fibroblasts and endothelial cells, which are more deeply embedded in the extracellular matrix and basement membrane than immune cells, and hence are more difficult to dissociate.13,22
Identification of signatures reflecting the intratumoral transcriptional heterogeneity of malignant cells
Clustering of malignant cells revealed tumor-specific clusters, suggesting a high degree of intertumoral heterogeneity (Fig. 2a). Over 700 genes were preferentially expressed in the individual tumors (Fig. 2b). Differentially expressed genes were enriched within pathways that varied across the tumors, showing significant phenotypic diversity (Supplementary information, Fig. S4). We then applied dimension reduction analyses, specifically non-negative matrix factorization (NMF), and identified a total of 44 metagenes that were preferentially co-expressed by subpopulations of malignant cells across tumors (Supplementary information, Table S4). Then, hierarchical clustering was used to characterize these 44 metagenes into gene expression signatures, and high concordance was shown among five signatures (Fig. 2c; Supplementary information, Table S5).
Fig. 2: Malignant cell clusters and common malignant signatures revealed in NPC.
a t-SNE plot of 7581 malignant cells from 11 patients (indicated by colors) reveals tumor-specific clusters. b A heatmap shows genes (rows) that are differentially expressed across 11 individual primary tumors (columns). Red: high expression; blue: low expression. Selected genes are highlighted. c A heatmap depicts the pairwise correlations of 44 metagenes derived from 11 tumors. Clustering identified five coherent malignant gene expression signatures across the tumors. d Each panel (from top to bottom) shows violin plots that depict the scores for one of the five malignant signatures for malignant cells from the 11 tumors. e Changes in gene expression for the five malignant signatures in response to different EBV infection statuses (95 positive vs 14 negative, detected by in situ hydridization with the EBV-encoded small RNAs) in NPC Cohort A are shown (P values were based on the Wilcoxon rank-sum test). The box plot center corresponds to the median, with the box and whiskers corresponding to the interquartile range and 1.5× interquartile range, respectively. f A bar plot shows the direction and statistical significance (P values were based on the Spearman correlation test) of the associations between each of the malignant signatures and stromal/intratumoral TILs in NPC Cohort A. g Kaplan–Meier curves for progression-free survival in the 88 patients in NPC Cohort A stratified according to high vs low expression of the cell cycling signature. Cox regression HR and 95% CI obtained after correcting for age, sex, smoking history and disease stage are shown; the corresponding Cox regression P value is also shown. h Prognostic values of the malignant signatures in the 88 patients in NPC Cohort A. Forest plots show HRs (blue/red squares) and CIs (horizontal ranges) derived from Cox regression survival analyses for progression-free survival in multivariable analyses adjusted for age, sex, smoking history and disease stage; the corresponding Cox regression P values are also shown. Significant results are indicated with red squares.
Of these signatures, the three epithelial differentiation signatures consisted primarily of epithelial genes, such as cytokeratins (e.g., KRT6), and small proline-rich proteins (e.g., SPRR2 and SPRR3) (Supplementary information, Table S5). The fourth signature was related to cell cycling, consisting of cell cycle-related genes, and the final signature was enriched for genes associated with cell secretion (Supplementary information, Table S5). The scores for the malignant signatures varied across the malignant cells from different tumors (Fig. 2d). Intriguingly, the signatures were preferentially, but not exclusively, co-expressed by subpopulations of malignant cells (Supplementary information, Fig. S5a). When clustering using only genes from the malignant signatures, the malignant cells were separated into different clusters corresponding mainly to the malignant signatures (Supplementary information, Fig. S5b, c). These results indicated that the malignant signatures we established capture common patterns of intratumoral transcriptional heterogeneity across nasopharyngeal tumors.
Correlations of malignant signatures with tumor characteristics and survival
Then, we analyzed the correlations between the malignant signatures and clinicopathological features in two independent cohorts with bulk gene expression data: NPC Cohort A (n = 113), which was previously described by Zhang and colleagues,10 and NPC Cohort B (n = 128), which was profiled by an Affymetrix HTA 2.0 microarray at Sun Yat-sen University Cancer Center. Significantly higher Epi-differ 2 scores were observed in EBV-positive tumors than in EBV-negative tumors (Fig. 2e). We then identified genes whose expression was strongly associated with the expression of Epi-differ 2 signature (Supplementary information, Table S6).23 Intriguingly, genes related to viral entry into host cell (i.e., TMPRSS2, PVRL4, CTSB, and SLC20A2) were revealed to be highly associated with the Epi-differ 2 signature, suggesting their potential roles in promoting epithelial differentiation by regulating EBV infection in NPC. Stromal tumor-infiltrating lymphocytes (TILs) were positively associated with the cell secretion signature (Fig. 2f), suggesting that high expression of the cell secretion signature may promote lymphocyte recruitment to the TME in NPC. In contrast, the negative correlation between the intratumoral TILs and the cell cycling signature may reflect an artifact of the tumor–immune interaction (Fig. 2f), with reduced infiltration of intratumoral TILs hampering the recognition and destruction of mutant tumor cells, in which mutations of G1/S cell cycle transition genes are commonly observed in NPC.7 Correlations between the malignant signatures and smoking history and T (tumor) category are also shown (Supplementary information, Fig. S6a). Of note, the Epi-differ 3 signature was positively associated with smoking history as well as the expression of MALAT1 (r = 0.46; P < 0.001), which has been reported to mediate lung carcinogenesis induced by cigarette smoke extract.24 This observation suggests the potential role of Epi-differ 3-related genes in smoking-mediated epithelial differentiation in NPC.
Furthermore, survival analyses revealed that high cell cycling scores significantly predicted poor survival in NPC Cohort A (Fig. 2g, h). This association with survival were validated in NPC Cohort B (Supplementary information, Fig. S7a), suggesting that the cell cycling signature predicts the aggressiveness of malignant cells in NPC. Additionally, to assess the prognostic role of the cell cycling signature in a different subtype of head and neck cancers, we analyzed the expression of the corresponding signature genes in a head and neck squamous cell carcinoma (HNSCC) cohort using a dataset from The Cancer Genome Atlas (TCGA HNSCC cohort), which comprises 366 HNSCC samples and no NPC samples25 (Supplementary information, Table S7). Intriguingly, the cell cycling signature did not predict survival in HNSCC (Supplementary information, Fig. S7b), indicating that this may not be the critical expression program in malignant cells that decides the aggressiveness of HNSCC.
Novel myeloid cell subtypes and regulators underlying monocyte-macrophage differentiation in NPC
In total, 5191 myeloid cells were clustered into six separate subsets (Fig. 3a): macrophages, monocytes, pDCs, CLEC9A+ DCs (DC1), CD1C+ DCs (DC2), and CCR7+ DCs (DC3) (Fig. 3b). Notably, DC1 expressed many of the classic markers of conventional type 1 DCs (cDC1), including CLEC9A, IRF8, and IDO1 (Fig. 3c).26 cDC1 constitute the key DC subtype responsible for cross-priming antitumor CD8+ T cells and are critical in antitumor immunity.27 DC3 expressed markers, such as CCR7, LAMP3, FSCN1, and CCL19, that have been associated with DC activation. Of these markers, CCL19 is known to attract naïve T cells as well as other DCs by interacting with CCR7 (Fig. 3c).28,29 DC2 preferentially expressed a series of classic marker genes of conventional type 2 DCs (cDC2), including CD1E, FCER1A and CLEC10A (Fig. 3c). Immunofluorescence confirmed the presence of distinct populations of DC1 and DC3 cells in NPC tissue samples (Supplementary information, Fig. S8). GZMB, which encodes the pro-apoptotic enzyme granzyme B, was highly expressed in pDCs (Fig. 3c), suggesting that the cells in this population may exhibit GZMB-dependent cytotoxicity.30 Interestingly, a relatively high proportion of pDCs was observed in our scRNA-seq data (2% in all cells, and 16% in myeloid cells) compared to other cancer types (1%, 3%, and 10% in myeloid cells for lung cancer,13 breast cancer,31 and HNSCC,32 respectively). This observation was validated by immunofluorescence analysis showing a mean percentage of 3% in all cells for pDCs, with a highly notable aggregation in the tumor stroma (Supplementary information, Fig. S3d).
Fig. 3: Myeloid cell clusters in NPC.
a t-SNE plot of 5191 myeloid cells color-coded by their associated clusters. b t-SNE plot, color coding for the expression of the marker genes (gray to red) for the indicated cell subtype. c Heatmap of genes with differential expression (rows) among the myeloid cell subtypes. d Differences in pathway activities scored per cell by GSVA among the different myeloid cell subtypes. The scores of pathways are normalized. e Heatmap showing the activity of TFs in each myeloid cell subtypes. The TF activity is scored using AUCell. f Expression of BACH1 and RUNX1 in primary human monocytes and monocyte-derived macrophages. Data are presented as the means ± SEM of three independent experiments (n = 4). g Expression of NR1H3 and TFEC in primary human monocytes and monocyte-derived macrophages. Data represent the means ± SEM of two independent experiments (n = 4). h Primary human monocytes were transfected with negative control (NC), NR1H3-specific, or TFEC-specific siRNAs, followed by immunoblot analysis to determine protein expression of NR1H3 or TFEC. β-actin is the loading control. The experiments were repeated independently for three times with similar results. i, l Left, representative histograms of CD14 expression levels in siNC and siNR1H3 (n = 6 donors, n = 4 independent experiments) (i), and siTFEC (n = 6 donors, n = 4 independent experiments) (l) monocytes stimulated by M-CSF. Right, expression levels of CD14 in siNC and siNR1H3 (n = 6) (i), and siTFEC (n = 6) (l) monocytes stimulated by M-CSF, which was determined by FACS (MFI, mean fluorescent intensity). j, m Left, representative histograms of CD86 and HLA-DR expression levels in siNC and siNR1H3 (n = 6 donors, n = 4 independent experiments) (j), and siTFEC (n = 6 donors, n = 4 independent experiments) (m) monocytes stimulated by M-CSF. Right, expression levels of CD86 and HLA-DR in siNC and siNR1H3 (n = 6) (j), and siTFEC (n = 6) (m) monocytes stimulated by M-CSF, which was determined by FACS (MFI). k, n Kaplan–Meier curves showing progression-free survival in the 88 patients in NPC Cohort A stratified according to high vs low expression of NR1H3 (k) and TFEC (n). Cox regression HRs and 95% CIs obtained after correcting for age, sex, smoking history and disease stage are shown; the corresponding Cox regression P values are also shown. *P < 0.05, **P < 0.01 (paired Student’s t-test).
We further characterized the functions of different myeloid cell subtypes by comparing pathway activities. Although macrophages and monocytes exhibited some common activated pathways, pathways involved in interferon (IFN)-α and IFN-γ responses were relatively upregulated in the macrophages, whereas angiogenesis-, NF-κB-mediated TNF-α signaling-, and hypoxia-related pathways were upregulated in the monocytes (Fig. 3d). These results indicate the potential tumor-promoting features of monocytes in NPC,33 consistent with their high expression of S100A8 and S100A9 (Fig. 3c).34 We then applied Single-Cell Regulatory Network Inference And Clustering (SCENIC) analysis35 to correlate transcription factors (TFs) with gene expression differences among cell types. This analysis identified a set of TFs implicated in the biology of different myeloid cell subtypes in NPC (Fig. 3e). Interestingly, macrophages and monocytes shared similar expression patterns for many TFs (Fig. 3e), suggesting that the macrophages may be derived from monocytes recruited to the NPC TME.36 Of note, the expression of genes regulated by BACH1 and RUNX1 was specifically upregulated in monocytes, whereas expression of NR1H3 and TFEC was prominent in macrophages (Fig. 3e). Next, we stimulated primary human monocytes with M-CSF to obtain macrophages, and observed significantly decreased expression of BACH1 and RUNX1, and increased expression of NR1H3 and TFEC in the monocyte-derived macrophages (Fig. 3f, g). Though previous reports have suggested the function of BACH1 and RUNX1 in promoting monocyte development,37,38 the roles of NR1H3 and TFEC remain unknown in the monocyte-to-macrophage differentiation process. To investigate the function of NR1H3 and TFEC, we subsequently transfected primary human monocytes with negative control (NC) or NR1H3/TFEC-specific siRNAs (Fig. 3h). After being triggered by M-CSF, the expression of CD14 was significantly higher in siNR1H3 (Fig. 3i) or siTFEC (Fig. 3l) cells than in siNC cells; the induction of CD86 and HLA-DR was also impaired in siNR1H3 or siTFEC cells (Fig. 3j, m). These results support that NR1H3 and TFEC facilitate the differentiation and maturation of monocyte-derived macrophages, suggesting that these regulators may promote antitumor immunity in NPC. Interestingly, further survival analyses revealed the association between high expression levels of NR1H3 and TFEC and improved outcomes in NPC patients (Fig. 3k, n). Our data provide a resource for further investigation of regulators underlying myeloid diversity in NPC.
T/NK cell phenotypes suggest regulators modulating antitumor immune response in NPC
T/NK cells (n = 17,741) represented the most prevalent cell type in NPC. Reclustering revealed 10 clusters: CD8+ T-1, CD8+ T-2, CD8+ T-3, dysfunctional T (CD8+ Tdysfunctional) cells, conventional CD4+ T (CD4+ Tconv)-1, CD4+ Tconv-2, CD4+ Tconv-3 cells, regulatory T cells (Tregs)-1, Tregs-2, and NK cells (Fig. 4a, b). Next, we compared the expression levels of selected T cell function-associated genes among the different T cell subtypes (Fig. 4c). Among the CD8+ T-1, CD8+ T-2 and CD8+ T-3 subtypes, high expression of cytokines and effector molecules was exhibited in CD8+ T-1, followed by CD8+ T-2 and CD8+ T-3 (Fig. 4c), suggesting that CD8+ T-1 represented the main cytotoxic population in NPC. CD4+ Tconv-2 represented naïve T cells, while exhaustion markers such as PDCD1 and TIGIT were highly expressed in CD4+ Tconv-3 (Fig. 4c). Monocle trajectory analysis of CD4+ T cells inferred a differentiation trajectory that mainly began with the CD4+ Tconv-2 naïve cluster and bifurcated into either the CD4+ Tconv-1 activation cluster or the CD4+ Tconv-3 dysfunctional/terminal differentiation cluster (Fig. 4d, e). For conventional CD8+ T cells, the differentiation trajectory also exhibited a branched structure, starting with CD8+ T-3 and T-2 that bifurcated into either a cytotoxic/activation lineage (mainly CD8+ T-1) or a dysfunctional lineage (mainly CD8+ Tdysfunctional) (Fig. 4f, g). Notably, we observed high expression of immune checkpoint molecules, except CTLA4, in the CD8+ Tdysfunctional cells, and LAG-3 and HAVCR2 were the most prominently expressed molecules (Fig. 4c, h). Of note, the expression levels of these immune checkpoint molecules were generally comparable in CD8+ Tdysfunctional cells in HNSCC12 (data not shown), suggesting distinctive features of dysfunctional T cells in NPC TME compared to those in other head and neck cancers. Intriguingly, exhaustion gene expression was highly associated with the expression of cytotoxicity markers (Fig. 4c), as LAG3 and HAVCR2 were found to be correlated with T cell activity, which was measured by the mean expression of granzymes (GZMA, GZMB and GZMH) (Fig. 4h). The expression of LAG3 and HAVCR2 also exhibited significant positive associations with intratumoral TILs (Supplementary information, Fig. S9a), and immunofluorescence analysis confirmed the higher percentage of LAG-3+ and HAVCR2+ CD8+ T cells in the total CD8+ T cells in tumor epithelial nests compared with those in the tumor stroma (Supplementary information, Fig. S9b). A trend of better survival was also shown for higher expression of LAG3 and HAVCR2 (Supplementary information, Fig. S9c). These results reflected an activation-dependent exhaustion expression program13,14 especially for LAG-3 and HAVCR2 in NPC, which represent potentially appealing checkpoint molecules for therapy.
Fig. 4: T/NK cell clusters in NPC.
a t-SNE plot showing 10 clusters of 17,263 T/NK cells (indicated by colors). b t-SNE plot, color coding for the expression of the marker genes (gray to red) for the indicated cell subtypes. c Average expression of selected T cell function-associated genes of naïve markers, inhibitory receptors, cytokines and effector molecules, co-stimulatory molecules, and Treg markers in each cell cluster. d Potential developmental trajectory of CD4+ T cells (n = 5694) inferred by analysis with Monocle 2. Arrows show the increasing directions of certain CD4+ T cell properties annotated with the signatures shown in e. e Traceplots of (left) CD4+ T cell activation signature along activation component and (right) terminal differentiation signature along terminal differentiation component for the CD4+ T cells. Cells are projected along the component, with the blue line indicating the moving average of the expression of signatures (a sliding window of length equal to 5% of the total number of CD4+ T cells was used), and the shaded area displaying SEM. Signatures used are presented in Supplementary information, Table S9. f Potential developmental trajectory of CD8+ T cells (n = 6975) inferred by analysis with Monocle 2. Arrows show the increasing directions of certain CD8+ T cell properties annotated with the signatures shown in g. g Traceplots (as in e) of (left) CD8+ T cell activation signature along activation component and (right) terminal differentiation signature along terminal differentiation component for the CD8+ T cells. Signatures used are presented in Supplementary information, Table S9. h Spearman correlation between the activity of CD8+ T cells (n = 6975), as measured by average granzyme expression (GZMA, GZMB and GZMH), and the expression of CD8+ T cell-specific genes. Genes encoding known immune checkpoint molecules are highlighted in blue. CD4+ Tconv, conventional CD4+ T cell; CD8+ Tdys, dysfunctional CD8+ T cell; NK, natural killer. i Heatmap showing the activity of TFs in each T/NK cell subtype. The TF activity is scored using AUCell. j Peripheral CD8+ T cells and NK cells were transfected with negative control (NC), EOMES-specific, RUNX3-specific, or XBP1-specific siRNAs, followed by immunoblot analysis to determine protein expression of Eomes, Runx3, or XBP1. β-actin is the loading control. The experiments were repeated independently for three times with similar results. k Left, representative histograms depicting the expression of GZMB and perforin on peripheral CD8+ T cells transfected with NC and EOMES-specific (n = 6 donors, n = 5 independent experiments), RUNX3-specific (n = 6 donors, n = 6 independent experiments), and XBP1-specific siRNAs (n = 6 donors, n = 5 independent experiments). Right, percentage of GZMB+ or perforin+ cells in siNC and siEOMES (n = 6), siRUNX3 (n = 6), and siXBP1 (n = 6) CD8+ T cells. l Left, representative histograms depicting the expression of GZMB and perforin on peripheral NK cells transfected with NC and EOMES-specific (n = 6 donors, n = 5 independent experiments), RUNX3-specific (n = 6 donors, n = 5 independent experiments), and XBP1-specific siRNAs (n = 6 donors, n = 5 independent experiments). Right, percentage of GZMB+ or perforin+ cells in siNC and siEOMES (n = 6), siRUNX3 (n = 6), and siXBP1 (n = 6) NK cells. *P < 0.05, **P < 0.01 (paired Student’s t-test).
SCENIC analysis revealed potential TFs underlying the regulation across subtypes (Fig. 4i). For example, NFATC1 was identified as a candidate TF underlying the gene expression differences in CD4+ Tconv-3 (Fig. 4i). NFATC1 is reported to regulate PD-1 expression following T cell activation,39 which is consistent with the upregulated PD-1 expression observed in CD4+ Tconv-3 and implicates NFATC1 as a potential regulator of CD4+ Tconv cells in their way to become dysfunctional in NPC. Interestingly, NK cells shared upregulated expression patterns of TFs such as EOMES, RUNX3, and XBP1 with cytotoxic CD8+ T cell subtypes (Fig. 4i), indicating the potential roles of these TFs in regulating the cytotoxic activity of both CD8+ T and NK cells in NPC. Previous studies have suggested that EOMES and RUNX3 contribute to the development and homeostasis of effector and memory T cells as well as NK cell differentiation.40,41,42 Nevertheless, the role of XBP1 in regulating cytotoxic lymphocytes is less well understood,43,44,45 and the effects of these regulators in immune responses have not been addressed in NPC. To validate the SCENIC results, we isolated CD8+ T cells and NK cells from the NPC tissue samples by flow cytometry (FACS), and identified a significantly higher percentage of cytotoxic GZMB+ or perforin+ phenotypes in the EOMES+, RUNX3+, and XBP1+ CD8+ T cells and NK cells (Supplementary information, Fig. S10a, b). Furthermore, reducing the expression of EOMES, RUNX3, and XBP1 in primary human CD8+ T cells and NK cells by transfection with corresponding siRNAs also downregulated the expression of GZMB or perforin compared with siNC cells (Fig. 4j–l), suggesting their function in regulating immune cytolytic activity. These data reveal compelling TFs enhancing the cytotoxicity of both CD8+ T cells and NK cells in NPC.
Diverse B cell subtypes identified in NPC
We detected 17,353 B cells that could be divided into seven clusters (Fig. 5a, b). Of these clusters, one consisted of rarely reported FCRL4+ memory B cells, in which elevated CCR1 expression supports the homing of these cells to inflamed tissues.46 One cluster comprised germinal center (GC) B cells, whereas another comprised plasma cells; it has been reported that GC B cells with relatively high affinity could be directed to become plasma cells.47
Fig. 5: B cell clusters in NPC.
a t-SNE plot showing 10 clusters of 17,353 B cells (indicated by colors). b t-SNE plot, color coding for the expression of the marker genes (gray to red) for the indicated cell subtypes. c Potential developmental trajectory of B cells (n = 17,353) inferred by analysis with Monocle 2. Arrows show the increasing directions of certain B cell properties annotated with the signatures shown in d. d Traceplots of (left) B cell proliferation signature along proliferation component and (right) antigen secretion signature along antigen secretion component for the B cells. Cells are projected along the component, with the blue line indicating the moving average of the expression of signatures (a sliding window of length equal to 5% of the total number of B cells was used), and the shaded area displaying SEM. Signatures used are presented in Supplementary information, Table S9. e Comparison of pathway activities (calculated based on GSVA) among different B cell subtypes. The scores of pathways are normalized. f Heatmap showing the activity of TFs in each B cell subtypes. The TF activity is scored using AUCell. g t-SNE plots of B cells color coded according to the expression of BCL6 and ATF4 or the AUC of the estimated regulon activity of these transcription factors, which corresponded to the degree of expression regulation of their target genes. h Left, FACS analysis of BCL6 expression in GC and other B cells from tumor tissues. The results represent five independent experiments (n = 6 donors). Right, association between GC B cells and BCL6+ B cells (n = 6) from tumor tissues. **P < 0.01 (paired Student’s t-test). Error bars, SEM. i Left, FACS analysis of ATF4 expression in plasma cells and non-plasma B cells from tumor tissues. The results represent four independent experiments (n = 5 donors). Right, association between plasma cells and ATF4+ B cells (n = 5) from tumor tissues. *P < 0.05 (paired Student’s t-test). Error bars, SEM.
Then, we analyzed all B cells to infer the potential developmental trajectory. The results showed that plasma cells congregated at the furthest terminus of the antigen secretion components, whereas FCRL4+ B cell/B cell-3/B cell-4 clusters were enriched at the end of the proliferation components, suggesting trajectories of functional divergence among these cell types (Fig. 5c, d). Differences in pathway activities and TFs among the different B cell subtypes are shown in Fig. 5e, f, respectively. GC B cells exhibited an activation of the pyruvate metabolism and MYC pathways, while the glycolysis pathway was upregulated in plasma cells (Fig. 5e). In addition, high activities of macrophage chemotaxis pathways were also detected in plasma cells, suggesting the potential roles of these pathways in promoting macrophage recruitment to the TME in NPC (Fig. 5e). Notably, the expression of genes regulated by TFs such as BCL6 and ATF4 was specifically upregulated in GC B cells and plasma cells, respectively (Fig. 5f, g). BCL6, which has been reported to be critical for GC B cell development, reactivates the B cell transcriptional program to enhance responses to external stimuli,48 whereas ATF4 plays a role in immunoglobulin production.49 These transcriptional phenotypes were confirmed in the NPC tissue samples by FACS (Fig. 5h, i). These results identified candidate regulators underlying the gene expression differences between GC B cells and plasma cells in NPC.
Complex intercellular communication networks in the NPC TME
Next, we used CellPhoneDB50 to identify ligand–receptor pairs and molecular interactions among the major cell types (Fig. 6a). Broadcast ligands for which cognate receptors were detected demonstrated extensive communication between tumor and immune cells (Fig. 6a, b).
Fig. 6: The dense network and multiple regulatory immune responses in the TME of NPC.
a Capacity for intercellular communication between malignant cells and immune cells. Each line color indicates the ligands expressed by the cell population represented in the same color (labeled). The lines connect to the cell types that express the cognate receptors. The line thickness is proportional to the number of ligands when cognate receptors are present in the recipient cell type. The loops indicate autocrine circuits. The map quantifies potential communication but does not account for the anatomical locations or boundaries of the cell types. b Detailed view of the ligands expressed by each major cell type and the cells expressing the cognate receptors primed to receive the signal. Numbers indicate the quantity of ligand–receptor pairs for each intercellular link. c–f Overview of selected ligand–receptor interactions of tumor cells (c), dysfunctional CD8+ T cells (d), macrophages (e), and the three types of DCs (f, DC1, DC2, and DC3). P values are indicated by circle size, with the scale to the right (permutation test). The means of the average expression levels of interacting molecule 1 in cluster 1 and interacting molecule 2 in cluster 2 are indicated by color. Assays were carried out at the mRNA level but were used to extrapolate protein interactions. CD4Tconv, conventional CD4+ T cell; CD8T, CD8+ T cell; CD8Tdys, dysfunctional CD8+ T cell; DC, dendritic cell; GCB, germinal center B cell; MAC, macrophage; MON, monocyte.
Notably, tumor cells expressed relatively high levels of chemokines (e.g., CCL20, CCL19, and CXCL10), while the corresponding receptors were widely expressed in immune cells, suggesting that these chemokines play significant roles in enhancing immune cell infiltration into NPC tumor tissue (Fig. 6c). Inhibitory interactions between tumor cells and T cells were commonly observed. In addition to CD274–PDCD1, other receptor–ligand pairs that have rarely been reported in NPC (e.g., NECTIN1–CD96, NECTIN2–TIGIT, PVR–CD96 and LGALS9–HAVCR2) were also identified (Fig. 6c, d). Specifically, we evaluated the interactions of LAG-3 and its ligands between T cells and tumor cells (Supplementary information, Fig. S11a). The previously reported FGL1–LAG-3 interaction51 was not obvious in NPC, primarily due to the low expression of FGL1 in the tumor cells (Supplementary information, Fig. S11b); MHC-II was predominantly responsible for the inhibitory function of LAG-3 induced by the interaction between dysfunctional T cells and tumor cells. We also identified putative communications between tumor cells and macrophages, such as the inhibitory or stimulatory interactions produced by CD47–SIRPA and CSF1–CSF1R (Fig. 6c, e).52
In addition, we identified specific ligand–receptor complexes predicted to regulate multiple immunomodulatory pathways coordinated by diverse subtypes of immune cells (Fig. 6d–f). High chemokine expression was detected in macrophages, DC1, DC2 and DC3, especially in macrophages and DC3, indicating their important roles in recruiting diverse types of immune cells (Fig. 6f). Interestingly, inhibitory communications between DC3 and T cells were also widespread (Fig. 6f), corresponding to the recent study suggesting this DC phenotype to be the most active immune-regulators of lymphocytes.53 Of note, DC1 expressed relatively high levels of XCR1, and the ligands of XCR1 (XCL1 and XCL2) were overexpressed by NK cells (Fig. 6f), suggesting the presence of functional interactions between NK cells and DC1. NK cells mediated the recruitment of DC1 to the TME of NPC, which is consistent with previous reports in a mouse model.54
Correlations of immune subtype-specific signatures with survival
Finally, we correlated the 17 immune subtype-specific signatures (Supplementary information, Table S3) derived from the 23 immune cell subtypes with clinicopathological features of NPC patients. Five immune signatures (i.e., pDC, DC1, NK cells, and FCRL4+ and plasma B cells) were associated with EBV infection, with all showing significantly lower scores in EBV-positive tumors than in EBV-negative tumors (Fig. 7a). We further explored their correlations with plasma EBV DNA concentrations. CLEC9A+ DC1 showed a significantly negative correlation with the EBV DNA level (r = −0.38; P < 0.001), suggesting that DC1 may be affected by active EBV replication (Fig. 7b). CLEC9A is reported to induce antitumor responses55,56 and modulate cytotoxic T lymphocyte responses to virus infection in mice57,58; however, the mechanisms underlying the influence of EBV infection on the CLEC9A+ DC-related immune response in NPC remain to be elucidated. Although DC1, FCRL4+ B, GC B and plasma cells were significantly associated with a relatively high number of mutations, most immune cell subtypes generally correlated poorly with mutational load (Fig. 7c), probably due to the overall low rate of non-synonymous mutations in NPC.7,9 With respect to tumor staging, a general reduction in many immune signatures was shown in patients with relatively high T and N categories and clinical stages (Supplementary information, Fig. S6b).
Fig. 7: Correlations of immune subtype-specific signatures with clinicopathological features and survival in NPC.
a Changes in gene expression for the indicated five immune signatures (pDC, DC1, NK, FCRL4+ B and plasma cells) with significant associations with the EBV infection status (95 positive vs 14 negative, detected by in situ hydridization with the EBV-encoded small RNAs) in NPC Cohort A. The box plot center corresponds to the median, with the box and whiskers corresponding to the interquartile range and 1.5× interquartile range, respectively. P values were based on the Wilcoxon rank-sum test. b Density dot plot and Pearson’s correlation analysis (r) of the gene expression for the DC1 signature and EBV DNA level in the NPC Cohort B (n = 128). c Bar plot showing the direction and statistical significance (P values were based on the Spearman correlation test) of the association between each of the immune cell subtypes and the number of mutations in NPC Cohort A. Significant associations are shown for the immune signatures DC1 and FCRL4+ B, GC B and plasma cells, which were positively correlated with the mutational burden. d Kaplan–Meier survival curves for progression-free survival in the 88 patients in NPC Cohort A stratified according to high vs low expression of six immune signatures (macrophage, pDC, DC1, NK cell, and plasma cell). Cox regression HRs and 95% CIs obtained after correcting for age, sex, smoking history and disease stage are shown; the corresponding Cox regression P values are also shown. e Prognostic values of immune signatures in the 88 patients in NPC Cohort A. Forest plots show HRs (blue/red squares) and CIs (horizontal ranges) derived from Cox regression survival analyses for progression-free survival in multivariable analyses adjusted for age, sex, smoking history and disease stage; the corresponding Cox regression P values are also shown. Significant results are indicated by red squares.
Subsequently, we explored the prognostic roles of these immune signatures in NPC. Generally, a consistent association between high scores for different immune signatures and improved survival was observed in NPC (Fig. 7d, e; Supplementary information, Fig. S7a). Importantly, positive associations were significant between survival and the signatures of macrophages, pDCs, DC1, NK cells, and plasma cells in NPC Cohort A (Fig. 7d, e). Furthermore, these associations were also evident in NPC Cohort B (Supplementary information, Fig. S7a), validating their reliability as prognostic biomarkers in NPC. Next, we assessed the prognostic value of these five immune signatures (Supplementary information, Table S7) in the TCGA HNSCC cohort. As CLEC9A+ DC1 was not identified in HNSCC studies, the NPC DC1 signature was used instead to infer its prognostic value in HNSCC. No credible prognostic values were identified for the immune signatures in HNSCC (Supplementary information, Fig. S7b), reflecting a difference in the prognostic value of these signatures between NPC and HNSCC.
Discussion
Through revealing the cellular heterogeneity of the TME, we established the intratumoral and immune subtype-specific signatures in NPC. High scores for the cell cycling signature were associated with poor survival, suggesting that this signature may reflect a highly proliferative and aggressive status in malignant cells. Importantly, high expression of the signatures of macrophages, pDCs, DC1, NK cells, and plasma cells significantly predicted a favorable prognosis in NPC. In this study, we observed high expression of CXC chemokine ligands, such as CXCL9 and CXCL10, that recruit CXCR3+ NK and CD8+ T cells into the TME, and upregulated activity of the IFN-α and IFN-γ response pathways in macrophages, suggesting potential antitumor capacity of these macrophages in NPC.59 The role of pDCs in immune regulation has been extensively debated as they tend to be tolerogenic and are associated with worse clinical outcomes in certain types of cancer, while they also participate in the priming of immunogenic adaptive responses.60 The favorable prognostic value of pDCs in NPC suggests the potential role of these cells in inducing antitumor immune responses in this disease. Of note, the signatures did not show a pronounced prognostic value in HNSCC patients. The differences in the TME between NPC and other head and neck cancers may account for the differences in the prognostic value of the signatures. In HNSCC, for example, the approximate proportions of DCs and B cells in immune cells are ~3% and ~10%, respectively,12,32 but their proportions in NPC according to our study are ~7% and ~40%, respectively. It should be noted that CLEC9A+ DC1 have not been reported previously in HNSCC, although the specificity of this population in NPC over other HNSCC merits further exploration. As immune cells are distributed unevenly within tumors, future studies with emerging technologies, such as spatial transcriptomics, are warranted to comprehensively elucidate the spatial heterogeneity of NPC TME.
The distinctive features of the TME may provide a foundation for the design of therapies targetting NPC. For instance, analyses of myeloid cells revealed compelling TFs regulating the monocyte-to-macrophage differentiation. Similarly, SCENIC analysis in CD8+ T cells and NK cells identified regulators that can enhance their cytotoxicity in NPC. This resource would help deepen our understanding of immune cell diversity and provide directions for future exploration of the transcriptional regulators in immune cell biology. Our data also suggest potential immunotherapeutic targets that may be used in combination with PD-1/PD-L1 blockade to further improve response rates in NPC.61,62 LAG-3 and HAVCR2, which have rarely been studied in NPC to date, were highlighted as the most prominent immune checkpoint molecules in CD8+ Tdysfunctional in NPC. Moreover, the cell–cell interactions suggest that tumor cells communicate extensively with immune cells and that inflammatory and adaptive immune responses are mediated by various immunomodulatory molecules (e.g., TIGIT, CD96, CD47, CSF1R, and XCL1/2-XCR1) that have rarely been studied before in NPC.
In summary, to the best of our knowledge, this work represents a unique resource providing a comprehensive single-cell transcriptome atlas of the multicellular ecosystem of NPC TME. This study lays a new foundation for the development of precision therapies in NPC.