Dissecting the heterogeneity of the microenvironment in primary and recurrent nasopharyngeal carcinomas using single-cell RNA sequencing.

Nasopharyngeal carcinoma (NPC) has a 10-15% recurrence rate, while no long term or durable treatment options are currently available. Single-cell profiling in recurrent NPC (rNPC) may aid in designing effective anticancer therapies, including immunotherapies. For the first time, we profiled the transcriptomes of ∼60,000 cells from four primary NPC and two rNPC cases to provide deeper insights into the dynamic changes in rNPC within radiation fields. Heterogeneity of both immune cells (T, natural killer, B, and myeloid cells) and tumor cells was characterized. Recurrent samples showed increased infiltration of regulatory T cells in a highly immunosuppressive state and CD8+ T cells in a highly cytotoxic and dysfunctional state. Enrichment of M2-polarized macrophages and LAMP3+ dendritic cells conferred enhanced immune suppression to rNPC. Furthermore, malignant cells showed enhanced immune-related features, such as antigen presentation. Elevated regulatory T cell levels were associated with a worse prognosis, with certain receptor-ligand communication pairs identified in rNPC. Even with relatively limited samples, our study provides important clues to complement the exploitation of rNPC immune environment and will help advance targeted immunotherapy of rNPC.

Introduction

Nasopharyngeal carcinoma (NPC), endemic to East and Southeast Asia, is biologically different from head and neck squamous cell carcinoma (HNSCC) and is associated with Epstein–Barr virus (EBV) infection.[1] Although chemo-radiotherapy can achieve a high response rate in primary NPC (pNPC), 15–58% of the patients experience NPC recurrence. Currently, no long term or durable treatment options are available. Traditional treatments for local and systemic control are limited to endoscopic surgery, reirradiation, chemotherapy, or targeted therapy, and severe late complications of salvage strategies can lead to treatment-related death.[2] Based on their abundant immune cell-infiltrated microenvironment and high PD-L1 expression, NPCs are generally classified as “hot” tumors and should theoretically respond well to immune checkpoint inhibitors (ICIs).[3,4] Anti-PD-1 ICIs have been demonstrated to be efficacious in recurrent and/or metastatic NPC, without severe toxicities, regardless of monotherapy or combination with chemotherapy.[5-7] Radiotherapy is generally acknowledged to play a critical role in the immunomodulation of tumor cells and the tumor microenvironment (TME). The TME of recurrent tumors in radiation fields is believed to be distinct from that of treatment-naïve malignancies, but the corresponding investigation is lacking.[8] Thus, exploring the complex TME could help better understand the molecular mechanism of recurrence and explore new therapeutic options for recurrent NPC (rNPC).

Recent advances in single-cell RNA sequencing (scRNA-seq) have offered new insights into the TME of NPC at a single-cell resolution.[9-12] Various cell components present in the complex tumor ecosystem could be further classified, along with the elucidation of their transcriptional features. In this study, for the first time, we aim to provide deeper insights into the differences of immune and cancer cells between primary and recurrent NPC by scRNA-seq and discover the complex role of TME in recurrence.

Materials and Methods

Patients and single-cell isolation

This study was approved by the Ethics Committee of the Shanghai Cancer Center, Fudan University, China (No. 1711178–9). Written informed consent was obtained from all patients. Six nasopharynx biopsies were obtained by endoscopy from four treatment-naïve patients with pNPC and two patients with rNPC.

Each fresh tumor sample was cut into pieces of 1 mm3 and washed with 1× Dulbecco’s phosphate-buffered saline (DPBS; HyClone) immediately after collection. The tumor tissue was enzymatically digested in Dulbecco’s modified Eagle’s medium (HyClone) containing 5% fetal bovine serum (HyClone), 0.2% collagenase type IV (#C5138; Sigma–Aldrich), 0.05% hyaluronidase type I-S (#H3506; Sigma–Aldrich), and 0.002% DNase I (#DN25; Sigma–Aldrich). After 45 min of digestion at 37°C in a shaking water bath, the enzymatic hydrolyzate was filtered through a 40-μm pore size nylon mesh and centrifuged at 1,500 rpm for 5 min to obtain a single-cell suspension. Next, erythrocytes were lysed using a red blood cell lysis solution (#C3702; Beyotime) for 5 min. Finally, cells were washed twice with DPBS, and cell viability was assessed using 0.1% trypan blue staining.

scRNA-seq library construction and data preprocessing

The Cell Ranger software pipeline (version 3.1.0) provided by 10× Genomics was used to demultiplex cellular barcodes, map reads to the genome and transcriptome using the STAR aligner, and down-sample reads as required to generate normalized aggregate data across samples, producing a matrix of gene counts versus cells. Using the R package Seurat (version 3.2.1),[13] low-quality cells and likely doublets were removed; only cells with UMI/gene numbers within two standard deviations from the mean on either side and <20% of the counts belonging to mitochondrial genes were retained. After quality control, 62,164 single cells were included in downstream analyses. Using the NormalizeData function of Seurat, a log-transformed gene expression matrix was generated. After principal component analysis of the top variable genes, cells were clustered using a graph-based clustering approach and visualized in two dimensions using t-distributed stochastic neighbor embedding (t-SNE). By repeating dimensionality reduction and clustering, subclusters of T/natural killer (NK), B, myeloid, and epithelial cells were identified and annotated based on the average gene expression of canonical markers. Likely doublets expressing both epithelial and T/B/myeloid cell markers were further excluded from malignant cell analysis.

Using the FindMarkers function of Seurat, differentially expressed genes (DEGs) were identified. The scores of functional modules for each cell were calculated using the AddModuleScore function of Seurat. For normally distributed data, P values were evaluated using a two-sided Student’s t-test. Gene set variation analysis (GSVA package, version 1.30.0)[14] was used to estimate pathway activities of cell groups. Differences in pathway activities scored per cell were calculated using the LIMMA package (version 3.38.3). Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes pathway enrichment analyses and gene set enrichment analysis (GSEA) were performed using the clusterProfiler package version 3.6.0.[15] To establish a protein–protein interaction (PPI) network of DEGs identified in malignant cells, STRING[16] was used to retrieve interacting genes, and the Cytoscape software was used to visualize the network.[17] Significant modules and hub genes in the PPI network were identified using the MCODE plugin.[18]

The initial copy number variation (CNV) of each region was calculated using the inferCNV R package (version 1.2.1).[19] Epithelial cells were selected as malignant cells using immune cells as a reference. After denoizing, the CNV score of each cell was calculated as a quadratic sum of CNV regions. The CopyKAT R package (version 1.0.5)[20] was also used to confirm the CNV results without identifying the reference cell type.

Using the Monocle R package (version 2.14.0),[21] the potential lineage differentiation trajectory was determined. To select ordering genes (qval < 0.01) along the pseudotime trajectory, we converted the raw count into the CellDataSet object and used the differentialGeneTest function of Monocle. After dimensional reduction and cell ordering, the immune cell differentiation trajectory was inferred, and gene expression was plotted to track changes over pseudotime. RNA velocity was measured to investigate the potential cellular state dynamics of the myeloid lineage. Using the velocyto pipeline (https://github.com/velocyto-team/velocyto.py),[22] the spliced and unspliced reads were recounted, and cell fate decisions were represented by their transition probabilities from starting cells to neighboring cells. Velocity fields were projected onto the t-SNE plot obtained in Seurat.

Transcription factor (TF) activity was analyzed using single-cell regulatory network inference and clustering (SCENIC, version 1.1.2.2)[23] with default parameters based on the motifs database for RcisTarget and GRNboost. Cell–cell interactions between different cell types were calculated using CellPhoneDB (version 2.0)[24] to identify significant ligand–receptor pairs within each sample. Any two cell types in which a ligand was expressed in one type and its receptor was expressed in the other type were linked.

Bulk RNA-seq data analysis

NPC bulk RNA-seq raw data were downloaded from the public Gene Expression Omnibus database (accession number: GSE102349). This cohort consisted of 113 samples, and only 88 samples had prognostic information. After data processing, gene expression levels were normalized using transcripts per kilobase million. By randomly selecting 50 cells of each immune cell type, we generated a single-cell reference file and calculated the numbers of individual cell types from the bulk RNA-seq data using CIBERSORTx (http://CIBERSORTx.stanford.edu/).[25] Cell type enrichment was dichotomized using the median score. The log-rank test and Kaplan–Meier analysis were performed to evaluate the prognostic value of each cell type. All analyses were performed using R version 3.6.3.

Results

Single-cell expression atlas and cell typing of nasopharyngeal carcinoma

Following the overall workflow (Fig. 1A), we generated scRNA-seq profiles for nasopharyngeal tumors from the six patients with NPC, four treatment-naïve, and two recurrent samples. The clinical characteristics of the patients, hematoxylin and eosin tissue histology and EBV-encoded RNA data are shown in Table S1. After standard data processing and quality control (metrics and mitochondrial count distribution are shown in Table S2), we obtained transcriptomic profiles for 62,164 cells. To construct a global TME atlas, cell classification and marker gene identification were performed using Seurat.

Eight major cell types were identified using marker gene expression projected on a t-SNE plot. The cell types included T/NK cells (30,612 cells, 49.2%, marked with CD3E/G and KLRF1), B cells (15,982 cells, 25.7%, marked with CD79A and MS4A1), myeloid-derived cells (4,982 cells, 8.0%, marked with CD68, LYZ, and CSF3R), mast cells (714 cells, 1.1%, marked with TPSAB1), epithelial cells (8,268 cells, 13.3%, marked with EPCAM and KRT18/19), endothelial cells (506 cells, 0.8%, marked with VWF and PECAM1), and fibroblasts (1,100 cells, 1.8%, marked with COL1A1) (Fig. 1B and C). The t-SNE plotting of all cells showed that immune and stromal cells from samples of different origins were mixed. Conversely, malignant epithelial cells were relatively separated, indicating high intertumoral heterogeneity of malignant cells.

Immune cell infiltration was abundant in NPC. Most of the tumor stroma consisted of epithelial cells, while T/NK and B cells comprised most of the tumor immune microenvironment (Fig. 1D). All these cell subtypes were shared among patients and between pNPC and rNPC samples, although in different proportions. The above data showed the existence of a complex ecosystem in NPC. Therefore, we further focused on functional cell types and identified over 50 different cell subclusters.

Active-dysfunctional coupled program and higher immunosuppressives potential of T cells in rNPC

A total of 30,612 T/NK cells were detected as the most prevalent cell type in the NPC samples, accounting for 49.2% of all cells. The T/NK cells were grouped into 15 clusters, each cluster consisting of both rNPC- and pNPC-derived cells (Fig. 2A). According to the function-associated genes, the 15 clusters were classified into five CD8+ T cell subtypes, six CD4+ T cell subtypes, one regulatory CD4+ T cell (Treg) subtype, two NK cell subtypes, and one CD4/CD8 double-positive T cell subtype (Fig. 2B).

In a total of 13,696 CD8+ T cells, we could not clearly distinguish a dysfunctional CTL subtype. Among CD8-C1 to C5 clusters, cytotoxic molecules (IFNG, RPF1, NKG7, GZMA, GZMB, and GNLY) were co-expressed with inhibitory markers (PDCD1, HAVCR2, LAG3, TIGIT, and CTLA4). This result suggested that T cell dysfunction is a gradual but not binary state in NPC. CD8-C5 showed higher expression of proliferative genes, such as MKI67 and STMN1, along with cytokine and exhaustion signals, indicating that the exhaustion and proliferation processes were twisted in NPC (Fig. 2B). In the T cell function-associated gene set calculation, CD8-C4-GNLY/CD8-C1-LAG3 presented the highest cytotoxicity and exhaustion signature, while CD8-C2/C3 were comparatively naïve (Fig. 2 C). Monocle trajectory analysis of CD8+ T cells inferred a differentiation trajectory that mainly began with the CD8-C2 naïve cluster and bifurcated into either the CD8-C2-MKI67 cycling T cell cluster or the CD8-C1/C4 cytotoxic/activated cluster. During the differentiation process, beginning with naïve genes (CCR7, TCF7, KLRB1, and LEF1) and ending with proliferative genes (MKI67, TYMS and ZWINT), coactivators (TNFRSF14, ICOS, and TNFRSF9), coinhibitors (PDCD1, CTLA4, LAG3, HAVCR2, and TIGIT), and genes involved in the cytotoxic capacity (IFNG, GZMK, GZMA, GZMB) were all upregulated (Fig. 2D).

Compared to T cell subtypes in pNPC (41.6%), CD8+ T cell levels were elevated in rNPC (60.8%) patients. In trajectory comparison, early-stage naïve CD8+ T cells were mostly distributed in the pNPC samples, whereas in the rNPC samples, CD8+ T cells were primarily at terminal stages of the cytotoxic and exhausted states, which showed that the CTL response was triggered in rNPC (Fig. 2E). DEG analysis in CD8+ T cells showed elevated expression of exhaustion signatures, such as HAVCR2, CTLA4, TIGHT, and EOMES, along with increased expression of the cytotoxic gene IFNG and inhibitory gene SOCS3 in rNPC (Fig. 2E). Furthermore, with both active and dysfunctional signatures elevated in rNPC, the exhaustion/cytotoxicity ratio was higher in rNPC than in pNPC (regression coefficient: 0.302 vs. 0.215, respectively, P < .01) (Fig. 2F), which means that exhaustion acts as the main CTL program.

Next, GSVA was conducted to compare differences in the pathways. JAK/STAT, phospholipase D, which has previously been shown to play an important role in T cell activation,[26] and the NF-κB signaling pathway were upregulated in rNPC-derived CD8+ T cells, suggesting an activated immune response. The upregulated p53, apoptosis, and cell cycle pathways also showed activation of cell cycle activity in rNPC CTLs (Fig. 2G). A SCENIC analysis showed that cell cycle and exhaustion-related transcript factors E2F2, EOMES, IRF1/2, JUN, and STAT1/2/3 were significantly upregulated, while NFKB1, JUND, and FOXP1 were significantly downregulated in rNPC-derived CD8+ T cells. The alterations in the TFs further confirmed the coupled exhausted and activated state of CTLs in rNPC (Fig. 2G).

Among 14,908 CD4+ T cells, resident memory (CD4-C1-CD69), follicular helper (CD4-C2-CXCL13+), naïve (CD4-C3-CCR7+), conventional (CD4-C4/C5), and cycling (CD4-C6-MKI67) T cells, as well as Tregs (FOXP3+), were identified. The percentage of Treg in rNPC-derived CD4 + T cells (52.9%) was higher than that in pNPC-derived CD4 + T cells (31.3%) but not statistically significant (P= .267). These CD4+ T cell subtypes each had specific pathway activities and TFs (Fig. S1A and B).

In addition to the relatively higher proportion of Tregs, which indicated an immunosuppressive environment, the total IL2R score, inhibitory score, and costimulatory scores were all increased in rNPC-derived CD4+ T cells, suggesting that the recurrent samples had a stronger immune-suppressing potential and were in a more activated state than primary samples (Fig. 2I). GSEA showed that the T helper (Th) 17 immune response and the interferon (IFN)-γ pathway were upregulated in the rNPC samples (Fig. 2J). Even though we could not clearly distinguish a Th1 or Th17 cluster, pathway enrichment showed a positive shift toward Th1/Th17 activity in rNPC-derived CD4+ T cells. Nonetheless, no prior study has been conducted concerning the Th17 role in NPC, which warrants further investigation.

Increased cytotoxicity and higher suppression of NK cells in rNPC

In total, 2,008 NK cells were detected in the NPC samples, accounting for 6.5% of all cells. NK cells were further grouped into five clusters and classified into two main types, CD56brightCD16− and CD56dimCD16+ NK cells (Fig. 3A). CD56dimCD16+ NK cells were characterized by high FCGR3A (encoding CD16) expression and low NCAM1 (encoding CD56) expression, and acted as a cytotoxic NK subtype, with overexpression of effector genes (GNLY, GZMA, GZMB, and PRF1) (Fig. 3B).

The composition of NK cells was similar between pNPC and rNPC samples, with CD56dimCD16+ NK cells comprising 65.3% and 70.3% of pNPC and rNPC samples, respectively. However, rNPC showed elevated cytotoxicity and chemotaxis of the CD56dimCD16+ NK cell-related pathway activities. The DEG analysis showed that the chemokine CXCR4, KIR inhibitory receptors KLRD1 and KIR2DL1, and coinhibitory receptors LAG3, TIGIT, HAVCR2 were upregulated in rNPC, indicating the involvement of NK cell immune evasion mechanisms (Fig. 3 C). CXCL12/CXCR4 chemokine signaling has previously been shown to be essential for NK cell activation and development.[27]

A SCENIC analysis showed that STAT1/2, EOMES, KDM5A, and ETS1/2 were upregulated, while MXD4 and XBP1 were downregulated (Fig. 3D). KMD5A promotes NK cell activation by regulating IFN-γ production, inhibits SOCS1 expression, and promotes STAT4 activation.[28] EOMES and ETS1 control key checkpoints for NK cell maturation and terminal differentiation.[29,30] Changes in the TF expression further suggested the complicated active/suppressive NK cell function in rNPC.

Elevated immunoglobulin production and reduced antigen presentation by B cells in rNPC

Recent data have shown that B cells and plasma cells located in tumors or tumor-draining lymph nodes can play important roles in shaping antitumor immune responses.[31] A total of 15,982 B cells were detected, accounting for 25.7% of all cells. As a relatively homogeneous population in the TME, following markers reported in other scRNA studies,[32,33] B cells could be divided into four main types, namely, plasma B cells (XBP1+/MZB1+), MALT B cells (JCHAIN+/IGHA+/IGHM+), naïve B cells (MS4A1+/TCL1A+/IGHD+), and memory B cells (MS4A1+/TCL1A+/CD27+) (Fig. 3E). The B cell types were enriched in various immune regulation-related pathways and different transcription regulators (Fig. S2A and B).

Compared with the primary samples, the recurrent group had more plasma cells (34.1% vs. 53.2%, respectively) and MALT B cells (8.9% vs. 13.2%, respectively) and fewer naïve B cells (14.5% vs. 7.8%, respectively) and memory B cells (39.5% vs. 25.7%, respectively), which indicated an increased B cell response in rNPC (Fig. 3F). B cells, especially the naïve and memory subtype, can serve as antigen-presenting cells (APCs).[34,35] DEG analysis showed upregulation of the IGHG, IGLL, and IGLC immunoglobulin-encoding genes in plasma cells, suggesting elevated antibody production in rNPC. In rNPC-derived naïve B cells, HLA-DRB1 and HLA-DQA1, involved in antigen presentation were downregulated, as well as interferon-induced IFITM1, IFIT1, IFI6, and IFI44L (Fig. 3G). Direct comparison of the memory B cells of revealed strong activation of TNFA signaling via NF-κB, p53, hypoxia and apoptosis pathways among recurrent samples, indicating B cell survival, death, and differentiation; whereas antigen processing and presentation, IFN-γ signaling, IFN-α signaling and inflammatory response were activated in B cells in pNPC. All these results suggested that B cells in primary samples exhibited an inflammation-dominant gene expression pattern, while B cells in rNPC showed a stronger secretory-like phenotype.

Enrichment of rNPC in M2-polarized macrophages and LAMP3+ dendritic cells (DCs)

Next, we performed unsupervised clustering of 4,982 myeloid cells, of which 2,516 (50.5%) cells were from the four patients with pNPC, and 2,466 (49.5%) cells were from the two patients with rNPC. The average proportion of myeloid cells was higher in rNPC (11.7%) than in pNPC (6.1%), but the difference was not significant (P = .40). A total of nine clusters were identified within the myeloid lineage, including four macrophage clusters (Macro1, Macro2, Macro3, and Macro4), four DC clusters (pDCs, cDC1, cDC2, and mature cDCs), and one neutrophil cluster, based on functional markers (Fig. 4A and B). Each subcluster contained cells derived from all samples, but the myeloid cell compositions in pNPC and rNPC were different (Fig. 4B).

Macrophages were identified by the high expression of genes encoding C1Q family proteins and CD68. The percentage of Macro1 (APOE+) in the rNPC samples (54.1%) was higher than that in pNPC (38.7%), but the difference was not statistically significant. Macrophages are usually classified into the canonical proinflammatory, classically activated M1, and anti-inflammatory, alternatively activated M2 classes. We could not clearly distinguish M1 and M2 macrophages using known marker genes, such as FCGR3A (M1) and CD163 (M2), as they were expressed in both of these cells, consistently with a previous report on M1/M2 coupled activation in NPC[12] (Fig. S3A). Notably, calculation of M1 and M2 polarization scores using related gene sets[36] showed that the M1 scores tended to be similar among macrophage clusters, while the M2 score was upgraded in Macro1. With an elevated proportion of Macro1, the total M1 score in rNPC was downgraded, while the M2 score was upgraded (Fig. 4C). DEG analysis showed that the common markers of M2 polarization MMP9, CD274, PD-L2, and IL1RN were overexpressed in rNPC-derived macrophages.[37,38] (Fig. 4D). Further pathway enrichment analysis showed increased regulation of B cell, Th2-cell, and NK cell differentiation and IFN receptor activity in rNPC, indicating that the ability to shape the immune environment was enhanced (Fig. 4E). Thus, induction of a shift from M2 tumor-associated macrophages (TAMs) toward tumor-suppressive M1-type macrophages is a promising target for rNPC therapy.

As professional APCs, DCs are crucial components that induce and maintain antitumor immunity. Three distinct subsets of conventional DCs (cDCs) were identified, namely, mature cDCs (LAMP3+/CCR7+), a classical cDC1 subset (CLEC9A+), and a cDC2 subset (CD1C+/FCER1A+). Plasmacytoid DCs (pDCs; CLEC4C+) showed high expression of LILRA4, GZMB, and IL3RA (Fig. 4B). GSVA of different DC clusters showed that antigen processing and presentation, EBV infection, and cell adhesions were mostly enriched in cluster cDC1, including cells that recognize viruses and intracellular antigens. Meanwhile, apoptosis, NF-κB, p53, and cytokine–cytokine receptor interaction signaling pathways were significantly upregulated in mature DCs (Fig. 4 F). The DC subtypes showed variation between the patient groups, with mature DCs being more abundant in rNPC than in pNPC (70.1% vs. 15.0%, respectively, P < .001), while pNPC was enriched in classical cDC1/2 and pDCs (Fig. 4B and F).

Consistent with the elevated proportion of mature DCs, DEG analysis indicated that the top upregulated genes in rNPC DCs were maturation-related (LAMP3 and MARCKSL1), migration-related (CCR7, FSCN1, and SLCO5A1), and encoding chemokine ligands (CCL17, CCL19, and CCL22) that recruit immune cells, mostly Tregs (Fig. 4E).. Furthermore, immune regulatory signatures, consisting of CD274 (PD-L1), PDCD1LG2 (PD-L2), CD200, EBI3, IDO1, IL4I1, SOCS1, SOCS2, and SOCS3, were upregulated, while antigen presentation was reduced in rNPC DCs (Fig. 4G). Comparison of the hallmark pathways showed increased IL-18 production, T cell tolerance induction, and lymphocyte migration in rNPC-derived DCs. These observations suggested that rNPC DCs were more regulatory and tolerogenic cells, which restrained the activation of T cells (Fig. 4H). The TFs that were upregulated in rNPC-derived DCs were consistent with those of LAMP3+ mature DCs, including CREM, NFKB2, RELA, JUN, and KDM5A, which correlated with the NF-κB pathway activity and transcription of immunosuppressive molecules (Figs. 4H and S3B).

Using RNA velocity, which predicts the future state of individual cells through spliced and unspliced mRNAs, we obtained more insight into the potential fate of myeloid cells. On the t-SNE plot, macrophages exhibited a transition toward the M2 state, with the arrow pointing to Macro1 cluster. Similarly, RNA velocities of mature DCs were significantly higher than other DC groups, supporting their dynamic transition state, and the arrow indicated a positive shift of DCs toward a mature state (Fig. 4I).

Hallmark pathway changes in malignant cells of rNPC

As NPC originates from the nasopharyngeal epithelium, we extracted all epithelial cell clusters for downstream analysis and further excluded likely doublets expressing both epithelial and T/B/myeloid cell markers. In total, 6,795 epithelial cells from the six samples were grouped into 15 clusters after the extraction of doublets. The patient-specific expression patterns indicated high intertumoral heterogeneity (Figs. 5A and S4A). To distinguish between malignant and nonmalignant cells, we inferred large-scale CNVs based on scRNA-seq data. Using immune cells as a reference, epithelial cells showed altered CNVs. Next, we evaluated the relative expression density of genes on each chromosome, i.e., chromosome deletions or amplifications (Fig. 5B). Consistent with previously reported whole-genome sequencing results,[39] the inferred CNVs indicated deletions in chromosomes 1p, 3p, 11q, 14q, and 16q and copy number gains in chromosomes 3q, 5q, 12p, and 12q. Therefore, most epithelial cells were likely cancer cells because of copy number changes. The CNV levels varied among the patients and cell clusters (Fig. S4B and C). The inferCNV results were consistent with copyKAT (Fig. S4D).

The heatmap (Fig. 5C) shows the top 10 genes that were preferentially expressed in each malignant cell cluster. DEGs were enriched in pathways that varied across clusters, showing significant phenotypic diversity. Some clusters indicated common malignant programs shared among patients. Of note, cluster 14, originating from four patients, showed high expression of genes related to epithelial cell formation and cilia movement (CAPS, TPPP3, PIFO, and C9orf24), previously identified in lung cancer scRNA from ciliated airway epithelial cells.[40] GO analysis of marker genes of cluster 14 showed enrichment of cilium-related pathways (Fig. 5D). Histologically, the human nasopharynx is lined with a ciliated pseudostratified columnar epithelium, but cilia are not present in NPC cells. Our scRNA analysis showed intertumoral heterogeneity with some ciliated cells retained in the NPC microenvironment with relatively low CNVs (Fig. S4A). Similarly, cluster 11, which originated from two patients, showed high expression of ECM1, a glycoprotein secreted into the extracellular stroma and found to promote metastatic progression and invasion in multiple cancer types.[41] GO enrichment of marker genes of cluster 11 showed pathways related to cadherin binding and focal adhesion, which are involved in cell interactions with the extracellular matrix and are critical for the epithelial-to-mesenchymal transition (EMT).[42] A subset of highly proliferative cells was identified in clusters 5 and 7, which originated from all patients. Cell cycle-related genes (TOP2A, MKI67, and CCNB1) were highly expressed in clusters 5 and 7, with relatively high CNV levels. GO enrichment of clusters 5 and 7 indicated that these cells underwent active mitotic division and could be potential therapeutic targets.

Our analyses demonstrated that rNPC was associated with elevated proportions of Tregs, more exhausted CTL features, M2-polarized macrophages, and the LAMP3+ DC subset, which led us to hypothesize that the substantial differences between the pNPC and rNPC tumor ecosystems could be due to differences in malignant cells. We performed DEG and pathway enrichment analyses between malignant primary and recurrent tumor cells to elucidate these differences. TNFA signaling via NF-κB, p53, and apoptosis pathways was enriched in both pNPC and rNPC EBV-related malignancies. We also observed an enrichment of genes involved in immune response pathways (e.g., antigen processing and presentation, IFN-γ signaling, and IFN-α signaling) in rNPC. In contrast, genes upregulated in pNPC primarily belonged to malignant programs (e.g., EMT and MYC target V1) (Fig. 5E).

A PPI network was constructed using the STRING database to reveal the interactions between the DEGs upregulated in rNPC. Using the MCODE in Cytoscape, three modules that might play important roles in the characteristics of rNPC were detected (Fig. 5F). Module 1 correlated with antigen presentation, consisting of both MHC-I (HLA-A/B/C/E/F) and MHC-II (HLA-DQ/DP/DR) molecules. Upregulation of these genes in rNPC indicated decreased tumor immune escape and increased T cell recognition/activation, which further confirmed increased immunogenicity of rNPC cancer cells. Module 2 consisted of KRT genes encoding keratins, which are the typical intermediate filament proteins of the epithelium. KRT6/16/17 are constitutive keratins with a relatively high proliferative state and are induced upon stress, injury, or inflammation. KRT5 is commonly used to diagnose undifferentiated NPCs.[43] This module showed an altered epithelial differentiation state in patients with the complicated epithelial cell development process. Module 3 consisted of enzymes responsible for ATP synthesis (ATP5L/ATP5F1/ATP5H/ATP5G3) and the mitochondrial respiratory chain (COX5B/COX7B/COX4I1 and UQCRB/UQCRQ), consistent with increased oxidative phosphorylation (OXPHOS) and hypoxia in rNPC. With high metastatic and tumorigenic potential, cancer stem cells are more reliant on OXPHOS than the bulk and putatively nonstem components.[44] Thus, OXPHOS inhibitors could be used to target rNPC and alleviate therapeutically adverse tumor hypoxia. Apart from the hub genes listed, the heatmap of DEGs showed that some immune-related genes, such as chemokines (CXCL1, CXCL8, CXCL14, and CCL20) and inhibitors (IDO1), were highly expressed in rNPC, indicating a complex immune response of these cancer cells (Fig. S4E).

Correlation between cell composition and survival and cell communication differences

Using our scRNA-seq data as an annotated signature matrix, we estimated the immune cell abundances in 113 patients with NPC using RNA-seq data via CIBERSORTx (Fig. 6A). The patients were clustered into four groups based on immune cell composition, with cluster 4 enriched in a naïve and inhibitory TME and cluster 2 enriched in NK cells and CTLs, thus showing different immune microenvironment subtypes. Subsequently, we evaluated the hazard ratio for each immune cell type. We found that a higher abundance of Tregs correlated with worse progression-free survival (PFS), whereas other immune cell subtypes did not significantly correlate with survival (Fig. 6B).

To characterize intercellular interactions among NPC immune components, we inferred putative cell–cell interactions based on ligand–receptor signaling inferred from our high-resolution scRNA-seq data using CellPhoneDB. Intensive communications were identified across all immune cell types, with macrophages and DCs showing the greatest number of interactions (Fig. 6C). The total receptor–ligand pair numbers were similar between pNPC and rNPC samples (Fig. S5).

To further investigate how Tregs interact with other components and influence the survival outcome, we visualized the major cellular communications of Tregs between pNPC and rNPC (Fig. 6D). Some common inhibitory, costimulatory, or chemokine communications with Tregs were identified in both pNPC and rNPC samples, while some specific communications were identified in rNPC but not in pNPC.

As Tregs are chemoattracted by chemokine gradients, increased chemokine interactions could explain the enrichment of Tregs in the rNPC samples. In the recurrent group, CCL22/CCR4, expressed by cDC2 and mature DCs, and CCL17/CCR4, expressed by mature DCs, could recruit Tregs into tumor tissue, consistently with the increased chemotactic ability of rNPC-derived DCs. CD8+ T cells and NK cells in rNPC highly expressed CCL5 and CCL4, which showed ligand–receptor binding to CCR5 and CCR4 on Tregs, respectively, suggesting a potential interaction and chemotaxis between Tregs and CTLs/NK cells. For costimulatory modules, CD8+ T cells and mature DCs in rNPC were uniquely predicted to interact with Tregs via TNFSF4–TNFRSF4. A subset of TNFRSF4+ Tregs was reported based on NPC scRNA-seq data to have high activation and immunosuppressive potential, explaining Treg activation in rNPC.[9] Inhibitory signals from Tregs to immune cells were commonly identified in both pNPC and rNPC, while suppressive signals from immune cells to Tregs were different. In pNPC, PDCD1LG2/FAM3C/CD274–PDCD1 interactions were commonly identified. In the absence of PDCD1 expression on Tregs, rNPC showed increased suppressive signals from DCs and macrophages to Tregs via LGALS9–HAVCR2. As TIM-3 (encoded by HAVCR2)-positive Tregs have been shown to express higher levels of IL-10 compared with those in TIM-3-negative Tregs and a higher capacity for inhibiting IFN-γ and TNF-α secretion by effector T cells, a recent study in HNSCC showed that treatment with anti-TIM-3 concurrently with anti-PD-L1 and radiotherapy led to a significant tumor growth delay, enhanced T cell cytotoxicity, decreased numbers of Tregs, and improved survival.[45] Taken together, the different Treg communications identified in rNPC can be potential targets for immunotherapy.

The C-type lectin family members NKG2A and NKG2C bind to HLA-E and are major inhibitory receptors expressed by NK cells in patients with both pNPC and rNPC. Other NK inhibitory receptors, namely KIR2DL1, KIR2DL3, and KIR3DL1, were significantly upregulated in the rNPC samples. As many KIR-related communication pairs were identified, we propose KIR antagonists as promising in rNPC treatment.

Discussion

Although immunotherapy provides some benefit in the treatment of rNPC,[5,46,47] our knowledge of rNPC is still based on the molecular and immune features of the primary tumor. Some studies have been carried out to compare the genomic alterations between pNPC and rNPC, but differences in tumor immune ecosystems between pNPC and rNPC remained to be discovered.[48-50] In this study, for the first time, we compared the tumor ecosystem of pNPC and rNPC using several sample representatives at single-cell resolution, and focused on the dynamic changes in immune cell subtype compositions and intercellular interactions. Our study showed an altered immune ecosystem in rNPC samples that was characterized by increased fractions or features of classic immunosuppressive cell such as Tregs, M2-polarized macrophages, and mature LAMP3+ DCs. Exhausted CTL features and elevated immune respond in tumor cells in rNPC samples further indicated an immune activated and dysfunctional state.

In a previous scRNA-seq analysis of breast cancer, elevated expression of markers of cytotoxic activity (PRF1 and GZMB) and exhaustion (HAVCR2 and LAG3) of CD8+ T cells, the Th1 signature of CD4+ T cells, and MHC-I/II of cancer cells in pretreatment biopsies, positively correlated with T cell expansion and potential clinical benefits during anti-PD1 treatment.[51] A subset of PD1+ CD8 + T cells expressing T cell activation markers also predicted a good respond to ICIs in NSCLC.[52] These signatures were consistent with those that we defined in rNPC and explained the good response rate to ICIs in recurrent/metastatic NPC. Consistent with previous scRNA-seq analyses of HCC, NSCLC, CRC, and melanoma, there is a transitional process in which effector T cells gradually lose cytotoxicity and become exhausted.[53-56] Here, we found that different CD8+ T cell clusters presented with diverse co-expression of cytotoxic activity and exhaustion markers. Furthermore, an overall analysis of rNPC samples showed that it presented with an elevated active-dysfunctional coupled state. In NPC clinical trials, the copy number gain of the cytotoxic effectors, GZMB and GZMH, was associated with increased survival during anti-PD1 treatment.[57] Thus, in rNPC samples, with elevated cytotoxic and exhausted features, there were sufficient evidence to conjecture rNPC respond better in ICI than primary tumor.

In a study of 95 paired rNPC and pNPC samples, recurrent tumors displayed higher levels of FOXP3+ lymphocytes.[58] Our data also showed an elevated proportion of Tregs in rNPC samples, with the bulk RNA analysis indicating that Tregs are a risk factor for a poor PFS. A subset of TNFRSF4+ Tregs was reported in NPC to have high activation and immunosuppressive potential.[9] Although we could not clearly distinguish a TNFRSF4+ Treg subset, the elevated inhibitory and costimulatory features in rNPC samples showed high immunosuppressive Tregs. Through a cell communication analysis, CTL/NK cells and mature DCs were predicted to recruit Tregs in rNPC. The CCL17/CCL22/CCL5 ligands were reported previously in gastric cancer and pancreatic ductal adenocarcinoma to shape the immunosuppressive microenvironment by recruiting Tregs, which could be further targeted for chemokine blockade therapy in rNPC treatment.[59,60]

Identified in multiple cancers using scRNA-seq, LAMP3+ DCs are considered as a subset of regulatory or tolerogenic DCs.[61-64] LAMP3+ DCs were found to shape the immunosuppressive ecosystem in rNPC samples in our study. In HCC, the LAMP3+ DC level was positively correlated with the infiltration of exhausted CD8+ T cells and Tregs.[62] In mouse models, LAMP3+ DCs can stimulate the clonal expansion of Th1-like cells.[65] These findings were consistent with the elevated exhausted CTL features, and Treg and Th1 activity in rNPC, which showed promising treatment effects by targeting LAMP3+ DCs through blocking CCL19 chemotaxis.[62] Similarly, M2-polarized macrophage signatures in rNPC resemble the findings of the C1Q+ TAM subset, which showed polarized phagocytosis, and antigen processing and presentation.[62,65-67] Cell communication and TCGA enrichment in CRC suggested that C1Q+TAMs could recruite and regulate CD8+ exhausted T cells, Tregs, and Th1-like cells,[65] which were consistent with features we found in recurrent samples. Transforming this macrophage subtype and turning M2 polarized TAM into M1 can be of potential value in future treatment of rNPC.

The features we identified in rNPC samples here were contrary with sc-RNA study of relapse HCC, in which they found rHCC presented with innate-like CD8+T cells with low cytotoxicity and decreased Tregs.[68] That is probably due to the radiation history in NPC. Radiation is known to increase MHC-I expression on the surface of tumor cells. Radiotherapy can kill cancer cells while simultaneously inducing tumor-associated neoantigen presentation, triggering the release of inflammatory cytokines, and increasing tumor-infiltrating immune cells, thereby transforming “cold” tumors into “hot” ones.[69] Thus, we also found elevated immune respond in tumor cells in rNPC samples. The differences also indicated the complex ecosystem between primary and recurrent tumors across different malignancy type. Furthermore, the globally elevated expression of the immune checkpoint genes, CTLA4, HAVCR2, and TIGIT in T cells, NK cells, and DCs from recurrent tumors suggested that targeting these markers might be more beneficial in recurrent tumor, but not in pNPC.

Because of small recurrent sample size, results driven in our study need to be further validated in a large patient cohort. Further studies of TCR/BCR and spatial information in rNPC will also help deepen our understanding of single-cell sequencing data. In conclusion, our study provides important clues to complement the exploitation of rNPC immune environment and can be a valuable resource for developing more effective therapeutic targets and biomarkers for immunotherapy in patients with recurrent NPC.

Click here for additional data file.