Abstract

   The heterogeneous nature of tumour microenvironment (TME) underlying
   diverse treatment responses remains unclear in nasopharyngeal carcinoma
   (NPC). Here, we profile 176,447 cells from 10 NPC tumour-blood pairs,
   using single-cell transcriptome coupled with T cell receptor
   sequencing. Our analyses reveal 53 cell subtypes, including
   tumour-infiltrating CD8^+ T, regulatory T (Treg), and dendritic cells
   (DCs), as well as malignant cells with different Epstein-Barr virus
   infection status. Trajectory analyses reveal exhausted CD8^+ T and
   immune-suppressive TNFRSF4^+ Treg cells in tumours might derive from
   peripheral CX3CR1^+CD8^+ T and naïve Treg cells, respectively.
   Moreover, we identify immune-regulatory and tolerogenic LAMP3^+ DCs.
   Noteworthily, we observe intensive inter-cell interactions among
   LAMP3^+ DCs, Treg, exhausted CD8^+ T, and malignant cells, suggesting
   potential cross-talks to foster an immune-suppressive niche for the
   TME. Collectively, our study uncovers the heterogeneity and interacting
   molecules of the TME in NPC at single-cell resolution, which provide
   insights into the mechanisms underlying NPC progression and the
   development of precise therapies for NPC.

   Subject terms: Cancer microenvironment, Head and neck cancer, Tumour
   heterogeneity, Tumour immunology
     __________________________________________________________________

   Nasopharyngeal carcinoma is a diverse cancer characterised by a
   heterogeneous microenvironment. Here, the authors use single cell
   sequencing to analyse the tumour microenvironment in 10 nasopharyngeal
   carcinoma tumours and identify different cell types including
   immune-suppressive T regulatory, tolerogenic dendritic, and exhausted
   CD8 T cells.

Introduction

   Nasopharyngeal carcinoma (NPC) is a distinct type of head and neck
   cancer, which has been closely linked with the infection of
   Epstein-Barr virus (EBV)^[64]1. NPC has a remarkable ethnic and
   geographic prevalence, where high incidence rate of 15–50 cases per
   100,000 people was reported in Southern China and Southeast Asia, as
   compared to 0.4 per 100,000 in western populations^[65]2. In general,
   patients with NPC are diagnosed with advanced stages largely due to
   non-specific symptoms. Radiotherapy is the primary treatment modality
   for patients with NPC because of the radiosensitive nature of its
   tumour cells^[66]2. Survival outcomes of patients with NPC improve
   substantially, reaching a 5-year overall survival rate of
   85.6%^[67]3,[68]4, mainly benefited from the evolution of radiotherapy
   techniques and the addition of platinum-based chemotherapy in patients
   with loco-regionally advanced disease. Nevertheless, more than 10% of
   patients develop recurrent and metastatic NPC^[69]2, for whom recent
   studies showed an overall response rate of 11.7% and 25.9–34% to
   targeted therapy with the addition of an inhibitor for epidermal growth
   factor receptor^[70]5 and immunotherapies using the immune checkpoint
   blockade^[71]6,[72]7, respectively. These variations in treatment
   responses and survival outcomes indicate the heterogeneous nature of
   NPC.

   Individual genetic makeup and genomic instability foster genetic
   diversity of cancer cells that contribute to tumour heterogeneity.
   Genome sequencing analyses have revealed diverse profiles of somatic
   alterations in NPC tumours, with high mutational frequencies at CYLD,
   NFKBIA, TP53, and CDKN2A/B, as well as accumulated mutations in MHC
   class I genes and chromatin modification genes, which were associated
   with poor overall survival of the patients^[73]8,[74]9. Apart from
   heterogeneous cancer cells, tumours exhibit another dimension of
   heterogeneity, which contain diverse normal cells creating the tumour
   microenvironment (TME) for the maintenance of cancer hallmarks^[75]10.
   Heterogeneous immune cells and stromal cells have been characterised
   using transcriptional profiling at single-cell resolution in several
   cancers, revealing that certain subtypes of immune cells and gene
   signatures in TME are important for tumour progression and sustained
   treatment responses^[76]11,[77]12. Profound infiltration of lymphocytes
   has been observed in histological biopsies of NPC, amid other stromal
   cells and tumour cells of different morphology^[78]13,[79]14. Moreover,
   high density of tumour infiltrating lymphocytes (TILs) was associated
   with favourable survival outcomes of patients with NPC^[80]15,[81]16.
   However, the composition of diverse cell populations in the TME has not
   been well illustrated in NPC. Two recent studies demonstrated the
   presence of T cells with various functional states and different immune
   cells in NPC tumours, using single-cell transcriptome
   analysis^[82]17,[83]18.

   In this work, we aim to provide a comprehensive global view of tumour
   heterogeneity of NPC, by analysing the single-cell transcriptional
   profiles of 176,447 cells from 10 treatment-naïve patients with NPC.
   With the combination of T cell receptor (TCR) repertoire information
   and individual tumour-blood sample pairs, we further characterise the
   clonality and migrations of T cells. In addition, we generate a
   potential cellular interaction network of cell populations in the TME
   of NPC.

Results

Landscape view of cell composition in tumour biopsy and PBMC in patients with
NPC

   To shed light on the complexity of tumour microenvironment (TME) in
   NPC, we performed single-cell RNA sequencing in combination of TCR
   repertoire sequencing on viable cells derived from tumour biopsies and
   matched peripheral blood mononuclear cells (PBMC) for 10 patients with
   EBV-positive NPC prior to any anti-cancer treatment (Fig. [84]1a,
   Supplementary Fig. [85]1a, b, and Supplementary Table [86]1). On
   average, we obtained more than 380 million sequencing reads for each
   sample, with the median of sequencing saturation (covering the fraction
   of library complexity) at 90.75% (75.90%–94.50%; Supplementary
   Table [87]2). After strict quality control filters (low expression of
   representative genes and inferred doublets; see Methods; Supplementary
   Fig. [88]1c and Supplementary Table [89]3), a total of 176,447 cells
   were identified from the 10 patients (including 82,622 and 93,825 for
   tumours and PBMC, respectively; Supplementary Data [90]1). We obtained
   about 1500 genes and 4950 unique molecular identifiers (UMIs) on
   average for each cell, indicating sufficient coverage and
   representative of transcripts (Supplementary Fig. [91]1d and
   Supplementary Data [92]1).

Fig. 1. The landscape profiling of single cells in NPC tumours and matching
PBMC.

   [93]Fig. 1
   [94]Open in a new tab

   a An experimental scheme diagram highlights the overall study design.
   Single viable cells were collected using flow cytometry sorting (FACS)
   and subjected for cell barcoding. The cDNA libraries of 5’-mRNA
   expression and TCR were constructed independently, followed by high
   throughput sequencing and downstream analyses. b UMAP plot of 176,447
   single cells grouped into six major cell types (left panel) and the
   normalised expression of marker genes for each cell type (right panel).
   Each dot represents one single cell, coloured according to cell type
   (left panel), and the depth of colour from grey to blue represents low
   to high expression (right panel). c UMAP plot of the above single cells
   coloured according to their origins from peripheral blood or tumour
   (left panel), and the fraction of cell types originating from each
   patient (right panel). Each dot represents one single cell, coloured
   according to sample origin.

   Next, to define groups of cells with similar expression profiles, we
   performed unsupervised clustering analysis implemented in Seurat
   software^[95]19. The distribution of cell clusters for each patient
   matched well with that of other patients, suggesting that the potential
   variation of expression due to batch effect of sample processing was
   minimal (see Methods; Supplementary Fig. [96]1e). Each cluster was
   further identified as a specific cell subpopulation according to the
   expression of the most variable genes and the canonical markers,
   including CD4^+ T cells (gene markers: PTPRC, CD3D, and CD4), CD8^+ T
   cells (PTPRC, CD3D, and CD8A), myeloid cells (CD14, ITGAX for CD11C),
   malignant cells (EPCAM and KRT5), B cells (CD19 and MS4A1), and NK
   cells (FCGR3A and NCAM1; Fig. [97]1b). Besides, we detected 56
   fibroblasts and seven endothelial cells with sparse distribution among
   TME cells (Supplementary Fig. [98]1f), which might reflect their
   intrinsic nature in NPC or their low representation due to technical
   limitations. All these cell types were widespread in tumour samples,
   indicating the heterogeneous cell composition of TME in NPC
   (Fig. [99]1b), consistent with a recent single-cell transcriptome study
   of NPC^[100]17. We observed that the proportions of CD8^+ T and B cells
   were increased with 1.34 and 2.33 times on average, respectively, while
   the NK cells were decreased in the tumours compared to the PBMC,
   suggesting two distinct immune landscapes between tumour and peripheral
   blood (Fig. [101]1c). Moreover, we compared cell compositions between
   NPC and other types of cancers with single-cell data publicly
   available, including non-small-cell lung cancer (NSCLC), colorectal
   cancer (CRC), and pancreatic ductal adenocarcinoma (PDAC). We observed
   the common occurrence of infiltrating immune cells amid individual
   heterogeneity of cell composition in all types of cancers
   (Supplementary Fig. [102]1g). Notably, we observed a significantly
   higher proportion of T cells in NPC compared to any other cancers
   (Supplementary Fig. [103]1h), which is consistent with a previous
   finding that the tumour infiltrating leucocytes was the main
   characteristic of NPC stroma^[104]20.

Heterogeneity of T cells and the diversity of TCR repertoire

   Considering the abundance of T and NK cells in NPC tumour samples and
   their anti-tumour capabilities, we explored the intrinsic structure and
   potential functional subtypes of the overall T and NK cell populations.
   We grouped all 141,875 T and NK cells into 32 subgroups using
   clustering analysis, of which the majority were CD4^+ and CD8^+ T cells
   (Fig. [105]2a, Supplementary Fig. [106]2a, Supplementary Data [107]1
   and [108]2). To identify any gene with a specific expression on a cell
   type, we performed differential gene expression (DEG) analysis of T
   cell clusters. We observed that CD4^+ and CD8^+ T cell clusters in
   tumour samples had widespread overexpression of exhaustion markers
   (LAG3, TIGIT, PDCD1, HAVCR2, and CTLA4) and effector molecules (GZMB,
   GZMK, INFG, NKG7, GNLY, and IL2; Supplementary Fig. [109]2b), with
   remarkably high expression of the proliferative signatures for
   CD8_C10_MKI67 and Treg_C3_MKI67 (Supplementary Fig. [110]2c and
   Supplementary Data [111]3). The co-expression of exhaustion and
   effector genes in tumour infiltrating T cells has been also
   demonstrated recently in NPC^[112]17. Together, these observations
   suggest that T cells exhibited anti-tumour effects against antigens,
   but their effector functions were somehow suppressed in the TME of NPC.
   By contrast, we observed naïve gene signatures (high expression of
   TCF7, SELL, CCR7 and LEF1) in the resting T cells in PBMC, including
   CD4_C1_LEF1, CD8_C1_LEF1, CD8_C2_TCF7, and DN_LEF1 (CD4^−CD8^−) cell
   clusters, and partially in Treg_C1_SELL and DP_TCF7 (CD4^+CD8^+;
   Supplementary Fig. [113]2b). Besides, we identified two clusters of
   CD4^+ Th1-like cells in tumours, including Th1_like_C1_CCR7 and
   Th1_like_C2_TNF, with specific expression of naïve T cell markers and
   pro-inflammatory cytokines, respectively, as well as a common
   expression of Th1-like cell markers^[114]21 (CXCL13, BHLHE40, and
   CXCR3; Supplementary Fig. [115]2d).

Fig. 2. Expression profile and development of CD8^+ T cells.

   [116]Fig. 2
   [117]Open in a new tab

   a UMAP plot of 141,875 T and NK cells grouped into 32 cell types. Each
   dot represents a cell, coloured according to the cell types indicated
   at the right legend. b Violin plots showed the normalised expression of
   CD8^+ T cell markers (rows) in each CD8^+ T cell cluster (columns).
   Cell clusters and the expression level of a gene are indicated at the
   x- and y-axis, respectively. c Pseudotime trajectory analysis of
   selected CD8^+ T cells (CD8_C5, CD8_C7, CD8_C8, CD8_C9, CD8_C10, and
   CD8_C11; n = 10,000) with high variable genes. Each dot represents one
   single cell, coloured according to its cluster label. The inlet plot
   showed each cell with a pseudotime score from dark blue to yellow,
   indicating early and terminal states, respectively. For CD8^+ T cell
   clusters, 10,000 cells were randomly selected for the analysis. d Box
   plots showed the transition-index scores of exhausted CD8^+ T cells
   (CD8_C11_PDCD1) and other CD8^+ T cells (n = 10). Comparison was made
   using a two-sided Wilcoxon test. Cell clusters and transition-index
   scores are indicated at the x- and y-axis, respectively. Endpoints
   depict minimum and maximum values; centre lines denote median values;
   whiskers denote 1.5× the interquartile range; coloured dots denote each
   patient. e Box plots showed the expansion- (top panel) and PBMC-Tumour
   migration-index (bottom panel) scores of each CD8^+ T cell cluster
   (n = 10). Each comparison was made using either a two-sided Wilcoxon
   test (top panel) or Kruskal–Wallis test (bottom panel). Cell clusters
   are indicated at the x-axis, and the y-axis shows the expansion- and
   PBMC-Tumour migration-index scores at the top and bottom panel,
   respectively. Endpoints depict minimum and maximum values; centre lines
   denote median values; whiskers denote 1.5× the interquartile range;
   coloured dots denote each patient.

   Next, we performed T cell receptor (TCR) repertoire analysis based on
   the sequences of α and β chains of TCR, which revealed 38,720 (32.97%;
   out of 117,447) T cells with detectable TCR α-β pairs or clonotypes
   after strict quality control (Supplementary Fig. [118]2e and
   Supplementary Data [119]4). We observed no sharing of any identical TCR
   clonotype among different patients with NPC, although they had certain
   preferences of V and J fragments as well as V^–J pairs (Supplementary