Immune Landscape of Viral- and Carcinogen-Driven Head and Neck Cancer

Abstract

Head and neck squamous cell carcinoma (HNSCC) arises through exposure to environmental carcinogens or malignant transformation by human papillomavirus (HPV). Here, we assessed the transcriptional profiles of 131,224 single cells from peripheral and intra-tumoral immune populations from patients with HPV^- and HPV⁺ HNSCC and healthy donors. Immune cells within tumors of HPV^- and HPV⁺ HNSCC displayed a spectrum of transcriptional signatures, with helper CD4⁺ T cells and B cells being relatively divergent and CD8+ T cells and CD4+ regulatory T cells being relatively similar. Transcriptional results were contextualized through multispectral immunofluorescence analyses and evaluating putative cell-cell communication based on spatial proximity. These analyses defined a gene expression signature associated with CD4⁺ T follicular helper cells that is associated with longer progression-free survival in HNSCC patients. The datasets and analytical approaches herein provide a resource for the further study of the impact of immune cells on viral- and carcinogen-induced cancers.

Keywords: cancer immunology; head and neck cancer; immunotherapy; multispectral immunofluorescence MC; mutagen-driven cancer; single-cell RNAseq; tertiary lymphoid structures; transcriptomics; viral-induced cancer.

Figure 1.. Assessment of bulk changes in transcriptional profiles between patients and overall clustering and identification of single cells.

A total of 131,224 single cells were recovered from 63 samples. (A) Live CD45⁺ cells were sorted from PBMC and single-cell suspensions prepared from tissues samples, and were subjected to single-cell RNAseq. (B) Clustering of samples by pseudobulk expression profiles reveals strong separation between PBMC and TIL, and distinct patterns of clustering of TIL samples by the tissue of origin. (C) FItSNE visualization and DRAGON clustering of all single cells identified 26 unique clusters across all samples. (D) The same FItSNE plot as (C), but with all immune cell types identified. (E) The same FItSNE embedding as (C) and (D) with all cells from each sample type shown. Shifts in density of cells are evident between PBMC, tonsil, and TIL samples, reflective of differences in transcriptional profiles between PBMC and tissues. (F) Quantification of differences between major immune lineages in HPV⁻ and HPV⁺ TIL. Each dot represents a subsample of 500 cells from PCA space for HPV⁻ and HPV⁺ TIL or a sample of 500 cells regardless of sample type (i.e. random), and the height of the bar is the mean of the subsamples. All comparisons were statistically significant due to 100 replicates of testing, but the mean fold-change between HPV⁻ and HPV⁺ TIL and random samples varied from 3.9-fold (B cells) to 1.1-fold (CD8⁺ T cells).

Figure 2.. CD8⁺ T cells share a differentiation trajectory towards co-expression of inhibitory receptors between HPV⁻ and HPV⁺ TIL.

A total of 32,734 CD8⁺ T cells were recovered across all samples. (A) Clustering of CD8⁺ T cells by DRAGON revealed a total of 8 clusters across all samples. (B) Enrichment of clusters by sample types showed that clusters 1-4 were predominantly of TIL origin, while clusters 5-7 were composed of mixtures of PBMC, tonsil and TIL, and cluster 8 was predominantly TIL. (C) Differential gene expression analysis revealed signature genes associated with each cluster. Specifically, cluster 1 expressed genes associated with the cell cycle, cluster 2 was associated with interferon responses, and clusters 3 and 4 expressed immune checkpoints. (D) Gene set enrichment analysis to evaluate biological functions of each cluster. Interestingly, the 2 clusters of exhausted CD8⁺ T cells showing distinct biological functions. (E) Diffusion mapping embedding of CD8⁺ T cells revealed an activation trajectory beginning with naïve peripheral CD8⁺ T cells, and progressing to terminally differentiated CD8⁺ T cells. (F) The trajectories of CD8⁺ T cells from HPV⁻ and HPV⁺ TIL were largely overlapping.

Figure 3.. Dissection of transcriptional states and differentiation trajectories in CD4⁺ Tconv and Treg.

A total of 45,640 CD4⁺ T cells were recovered from all samples, with 41,889 CD4⁺ Tconv and 3,751 Treg. (A-B) Tconv cells were separated into 7 clusters by DRAGON. (C) Three-dimensional diffusion map embedding of all CD4⁺ Tconv samples reveals a branching between CD4⁺ Tconv from HPV⁺ and HPV⁻ TIL. (D) Two-dimensional regression planes show the differences in trajectories for CD4⁺ Tconv from HPV⁻ and HPV⁺ TIL. (E) Heatmap showing the top 50 differentially expressed genes from clusters 1 and 7, the two terminal branches by diffusion analysis. Cluster 1 is associated with at T follicular helper (TFH) phenotype, while cluster 7 has an effector memory phenotype. (F-G) We identified 6 clusters of CD4⁺ Treg cells from all patients. The majority (89%) of CD4⁺ Treg were derived from TIL. (H) Gene set enrichment revealed IFN-responsive clusters (2 and 4) and clusters enriched for TNF targets/signaling (3 and 6). (I) Diffusion map embedding of CD4⁺ Treg from HPV⁻ and HPV⁺ TIL revealed a consistent differentiation trajectory. (J) The density of CD4⁺ Treg along the DC1 axis was similar for HPV⁻ and HPV⁺ TIL, but a slightly higher frequency of CD4⁺ Treg from HPV⁻ TIL were earlier on the differentiation trajectory versus HPV⁺ TIL.

Figure 4.. Analysis of tonsil and TIL B cells reveals granular details of germinal center B cells, and a unique B cell population associated with HPV⁻ TIL.

A total of 16,736 B cells were recovered from all samples. (A-B) We identified a total of 11 clusters of B cell from tonsils, TIL and PBMC. (C) Gene set enrichment revealed a germinal center phenotype associated with clusters 1-4, enrichment of genes for plasma cells in cluster 5, and combinations of naïve, memory and switched B cells in other clusters. (D) Diffusion map embedding of all B cells colored by clusters as in (A). This three-dimensional embedding yielded axes related to germinal center formation (DC1), transition from naïve to memory B cells (DC4) and progression to plasma cells (DC3). Few HPV⁻ B cells progress along DC1 to become germinal center B cells. (E) Same diffusion map embedding of as in (D), but colored by sample types. (F) The majority of HPV⁻ B cells are concentrated on the right side of the DC1 axis, while HPV⁺ cells have a bimodal distribution along the DC1 axis (note log scale on the y axis in [F]).

Figure 5.. Unique states and potential plasticity of myeloid cells in the TME.

We identified a total of 26,599 myeloid cells consisting of 21,737 CD14+ monocytes/macrophages, 3,946 CD16+ monocytes, and 916 dendritic cells (DCs). (A) We identified 8 clusters of myeloid cells in PBMC, TIL and tonsil tissue by DRAGON. (B) Clusters 1 and 5-8 were enriched in tissue sites, while clusters 2-4 were largely present in PBMC. (C) Heatmap of the top 20 genes in each cluster was used to identify states of myeloid cells across clusters. Cluster 1 cells expressed characteristics of DCs, but also expressed IDO1. Cluster 6 was associated with a dendritic cell phenotype, while cluster 7 was strongly associated with cytokine and chemokine secretion. Cluster 8 was characterized by expression of complement and MRC1 (typically associated with type 2 macrophages). Cluster 5 appears to express genes associated with clusters 6-8, suggestive of a possible precursor relationship. (D-E) Three-dimensional diffusion map embedding of myeloid cells (colored by sample type [D] and clusters [E]) reveals separation between PBMC and TIL myeloid cells, but also differentiation from cluster 4 to 5 to 6/7/8. (F) The DC1 axis is associated with differentiation from peripheral to tissue myeloid cells, with overlapping portions from peripheral blood and tissue myeloid cells. (G) Two-dimensional regression planes show differentiation from cluster 5 to 6 and 7/8. It appears that cluster 6 can either differentiate from a spectrum of points along the 2D regression plane, or conversely that myeloid cells can transition from cluster 6 to cluster 7/8 states.

Figure 6.. Expression patterns of receptors and ligands differ between HPV⁻ and HPV⁺ TIL, and extensive putative cell-cell communication occurs in TIL.

Expression of genes for receptors and both cell-surface and secreted ligands was evaluated, as well as the cell-cell communication within and between cell types. (A) Log₂ fold change in ligand expression between HPV⁻ and HPV⁺ TIL reveals differences associated with chemokines/cytokines from HPV⁻ TIL and CD40LG and other B cell activating stimuli from HPV⁺ TIL. (B) Myeloid-associated receptors are highly expressed in HPV⁻ TIL, while receptors associated with B cells and germinal center activity are upregulated in HPV⁺ TIL. (C-F) Circos plots derived from the CellTalker algorithm, showing networks of interactions between ligands and receptors from different sample types. Extensive putative cell-cell communication is observable in healthy donor PBMC and tonsil tissue, and there is a massive expansion of cell-cell communication in TIL.

Figure 7.. Immunofluorescence analysis of sections from HNSCC provides insight into putative cell-cell communication based on spatial localization.

We performed immunofluorescence (IF) staining on tissue sections obtained from patients in our study. (A) Slides were clustered based on the frequency of cell types present, yielding a total of 5 clusters. (B-F) Representative IF images from each of the clusters identified in (A) shown on the left side of each panel, and the log odds of cell-cell proximity are shown for each image on the right side of each panel. A higher log odds of proximity implies that the cells tend to be in close in 2 dimensions. (G-H) CellTalker was used to identify putative interactions between the 7 clusters of CD4⁺ Tconv (Figure 3; CD4-1 to CD4-7) and 11 clusters of B cells (Figure 4; B-1 to B-11) identified by DRAGON. Ligand/receptor interactions between B cells and CD4⁺ Tconv were present in both HPV⁻ and HPV⁺ TIL, but interactions between germinal center B cells (B-1 through B-4) and effector TFH cells (CD4-1) were exclusively present in HPV⁺ TIL. (I) Progression free survival analysis based on gene set enrichment for TFH signature (defined in Figure 3E) from HNSCC patients from the TCGA. Patients with higher TFH enrichment had significantly longer progression free survival.