Single-cell transcriptomics reveals regulators underlying immune cell diversity and immune subtypes associated with prognosis in nasopharyngeal carcinoma

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.

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).

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.

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).

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⁺ T_conv, conventional CD4⁺ T cell; CD8⁺ T_dys, 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).

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.

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. CD4T_conv, conventional CD4⁺ T cell; CD8T, CD8⁺ T cell; CD8T_dys, dysfunctional CD8⁺ T cell; DC, dendritic cell; GCB, germinal center B cell; MAC, macrophage; MON, monocyte.

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.