Single-cell analysis of myeloid cells in HPV+ tonsillar cancer

The incidence of human papillomavirus-positive (HPV⁺) tonsillar cancer has been sharply rising during the last decades. Myeloid cells represent an appropriate therapeutic target due to their proximity to virus-infected tumor cells, and their ability to orchestrate antigen-specific immunity, within the tonsil. However, the interrelationship of steady-state and inflammatory myeloid cell subsets, and their impact on patient survival remains unexplored. Here, we used single-cell RNA-sequencing to map the myeloid compartment in HPV⁺ tonsillar cancer. We observed an expansion of the myeloid compartment in HPV⁺ tonsillar cancer, accompanied by interferon-induced cellular responses both in dendritic cells (DCs) and monocyte-macrophages. Our analysis unveiled the existence of four DC lineages, two macrophage polarization processes, and their sequential maturation profiles. Within the DC lineages, we described a balance shift in the frequency of progenitor and mature cDC favoring the cDC1 lineage in detriment of cDC2s. Furthermore, we observed that all DC lineages apart from DC5s matured into a common activated DC transcriptional program involving upregulation of interferon-inducible genes. In turn, the monocyte-macrophage lineage was subjected to early monocyte polarization events, which give rise to either interferon-activated or CXCL-producing macrophages, the latter enriched in advanced tumor stages. We validated the existence of most of the single-cell RNA-seq clusters using 26-plex flow cytometry, and described a positive impact of cDC1 and interferon-activated DCs and macrophages on patient survival using gene signature scoring. The current study contributes to the understanding of myeloid ontogeny and dynamics in HPV-driven tonsillar cancer, and highlights myeloid biomarkers that can be used to assess patient prognosis.

1 Introduction

Tonsillar squamous cell carcinoma, or short tonsillar cancer (TC), is a head and neck cancer (HNC) caused by abnormal growth of epithelial cells of the tonsillar mucosa. Risk factors include tobacco smoking and alcohol abuse as well as high-risk human papillomavirus (HPV) infection (1). The incidence of HPV⁺ TC is sharply rising (2) and currently accounts for > 70% of TC cases (3, 4). Biomarkers associated with the immune compartment of the tumor microenvironment (TME) have shown prognostic value in HPV⁺ TC. For instance, several studies have highlighted a correlation between tumor-infiltrating CD8⁺ T-cells, T regulatory cells (Tregs) (5, 6), and PD-L1 expression (7, 8), respectively, to clinical outcome. However, immune-checkpoint blockade strategies, aiming at inducing T-cell responses, have shown limited clinical efficacy in HNCs, irrespective of HPV status (9, 10). A different approach to treat cancer concerns harnessing anti-tumor immunity by targeting myeloid cells due to their ability to modulate lymphoid cell activity (11). Increased knowledge about myeloid cell heterogeneity, their subset-specific features and functions within the TME of HPV⁺ TC, and their potential impact on survival is warranted in order to develop new treatment strategies.

Myeloid cells are divided into conventional dendritic cells (cDCs), monocytes, macrophages, and granulocytes, each characterized by cell-type specific functions (12, 13). DCs act as sentinel cells in peripheral tissue by taking up antigens, transferring them into lymphoid organs, and triggering T-cell activation (12). In the context of HNC, both cDC1 and cDC2 express CCR7 at protein level, highlighting their migratory potential (14). Furthermore, cDC2 subsets were specifically shown to induce and correlate with Th1 CD4⁺ T-cell abundance in the TME (14). In contrast, cDC1 migrate via XCL1 and CCL5, presumably produced by tumor infiltrating NK-cells (15). Using a combination of gene-signature scoring in bulk HNC RNA-seq, these chemokines were shown to correlate with both cDC1 and NK-cell signatures (15), stressing the relevance of the cDC1-NK-cell axis in HNC. Human cDC1s are known for their superior capacity to mount anti-tumoral CD8⁺ T-cell responses through antigen cross-presentation (16). However, in vitro studies have shown that human cDC2 from blood and tonsils also can cross-present antigens upon appropriate stimulation (16, 17). Compared to DCs, macrophages are sensors of tissue damage and infection, which help in clearing apoptotic cells and activating T-cells in situ. The functional role of macrophages in HNC is not clarified, and several retrospective IHC studies report contradictory results regarding prognostic impact (18–21). Studies of other cancers have presented macrophages as highly plastic cells involved in extracellular matrix remodeling through MMPs (22), angiogenesis via VEGF-α production (23), and attenuation of immune responses through TGF-β (24). Finally, monocytes can have both pro- and anti-tumoral functions, and differentiate into monocyte-derived DCs (25, 26) and tumor-associated macrophages (11, 27–29).

Recently, myeloid cell diversity has been revisited through the use of single-cell RNA-sequencing (scRNA-seq) in blood (30, 31), tonsil (32), and spleen (33) as well as in a variety of tumors (33–35), tumor draining lymph nodes (14, 36), and peritoneal ascites from cancer patients (26, 28). Although these studies report highly overlapping DC subsets, they also highlight differences in DC biology depending on the tissue of origin. In blood, the DC compartment, in addition to cDC1 and cDC2, also consists of the CD14⁺ CD163⁺ DC3 and Axl⁺ Siglec6⁺ DC5 subtypes (30). Analyses of the DC lineage in inflamed tonsil and spleen have demonstrated three additional populations closely related to cDC2s: a small proliferating DC cluster, a CLEC10A^-CLEC4A⁺ cDC2 population, and an activated CCR7⁺ cDC2 subset (26, 33). The heterogeneity of the monocyte-macrophage lineage (Mono-Macs) is even greater, including up to 15 different communities in a recent cross-tissue meta-analysis (35). The diversity of the Mono-Mac lineage is partly related to the disparity of markers and clusters reported, but also due to the low resemblance of in vivo Mono-Macs with the well characterized in vitro M1/M2 macrophage models (37). Additionally, macrophage ontogeny has been postulated to condition Mono-Mac function, and tissue resident macrophage populations such as Langerhans cells (LC) have recently gained attention (36, 38, 39). Collectively, the divergences observed in scRNA-seq studies possibly relate to the differentiation and maturation of myeloid cells, which in tissue are greatly influenced by the local microenvironment.

Myeloid cells represent potential therapeutic targets in TC since they can orchestrate anti- and pro-tumor T-cell responses. Given their divergent T-cell polarization capacity (26, 31, 32, 40), there is a need to profile their diversity within the TME and at healthy steady state settings. Several groups have assessed the transcriptional profiles of the whole immune compartment of HNC biopsies using scRNA-seq (41–43). However, these studies are limited in their capacity to resolve myeloid cell heterogeneity due to low frequency of these cells in the TME. Furthermore, HNCs arise in different anatomical sites and are heterogeneous in terms of immune-infiltration, and thus, it is important to delineate tumor subsite-specific nuances. Lastly, published studies most often use blood or tonsil specimens from patients with local inflammation as control material, which are informative, but may not represent appropriate steady state lymphoid tissue controls.

In this study, we evaluate myeloid cell heterogeneity and maturation status in HPV⁺ TC (cf. paired healthy tonsils (HT)), and the impact of these features on survival of HNC patients. We implement droplet-based 10X Genomics scRNA-seq to characterize FACS-sorted CD45⁺CD13⁺HLA-DR⁺ myeloid cells, describing transcriptomically unique clusters, which include novel cDC2, LC, and Mono-Mac populations. Notably, we describe a preferential recruitment of pre-cDC1s in the TME, giving rise to cDC1s that mature into activated DCs. We further show that patients enriched in gene-signatures from cDC1s and interferon-activated DCs and macrophages have a higher five-year overall survival. Better understanding of the myeloid complexity and functionality in TC can improve the design of myeloid-targeted therapies, which potentially can overcome the limitations of T-cell oriented immunotherapies in TC.

2 Materials and methods

Collection of tumor and contralateral tonsillar tissue was approved by the Swedish Ethical Review Authority (ref. no. 2017/580), and all participating patients granted written informed consent. Paired biopsies from TC (n=15) and contralateral HT (n=8) were obtained at Lund University Hospital at the time of diagnosis (Table S1).

2.1 Cell isolation

TC and HT samples were cut into small fragments in RPMI 1640 medium (Thermo Fisher Scientific, Bremen, Germany) supplemented with 0.1 mg/mL gentamycin (Sigma-Aldrich, St Louis, MO). The tissue fragments were enzymatically digested with Collagenase IV (Sigma-Aldrich) (2.0 mg/mL) and DNase I (Sigma Aldrich) (200 Kunits/mL) for 20 minutes at 37°C. Cells were filtered using a 70 μm cell strainer (BD Biosciences, San Jose, CA). For scRNA-seq, CD45⁺ CD13⁺ HLA-DR⁺ myeloid cells were isolated using FACSAria IIu (BD Biosciences) by cell sorting into FACS tubes containing 3 mL of fetal calf serum (FCS) (Life Technologies, Carlsbad, CA). Sorted myeloid cells were washed twice and resuspended in 45 μL of PBS supplemented with 0.04% w/v bovine serum albumin (BSA).

2.2 Flow cytometry

Cells were stained with Fixable viability stain 620 (BD Biosciences) to assess cell viability. Cells were washed and non-specific binding was blocked with ChromPure mouse IgG whole molecule (Jackson ImmunoResearch, West Grove, PA) for 15 minutes at room temperature. Cells were immediately stained with the appropriate antibody panel (Tables S2, S3) for 20 min at 4°C, washed, and analyzed on a FACSAria IIu instrument (BD Biosciences) or Cytek Aurora (Cytek Biosciences, Fremont, CA).

2.3 scRNA-seq library preparation and sequencing

The scRNA-seq libraries were prepared according to the user guide manual (CG000204) provided by 10X Genomics’ Chromium Single Cell 3’Reagent Kit (v3.1) (10X Genomics, Pleasanton, CA). Briefly, sorted cellular suspensions of myeloid cells were loaded on a Single Cell chip to encapsulate single cells using the 10X Chromium controller. Encapsulated cells were then subjected to in-drop lysis and reverse transcription reactions. The emulsion droplets were disrupted, and barcoded cDNA was purified using silene magnetic dynabeads (Thermo Fisher Scientific, Waltham, MA), followed by 12-14 cycles of PCR amplification using a C1000 Touch™ Thermal Cycler with 96–Deep Well Reaction Module (#1851197, Bio-Rad. Hercules, CA). The amplified cDNA was purified with the SPRIselect reagent kit (Beckman Coulter, Brea, CA) and subjected to enzymatic fragmentation, end repair, A-tailing, size selection with SPRIselect, adaptor ligation, post-ligation cleanup with SPRIselect, and sample index PCR followed by cleanup with SPRIselect. Library size and quality were assessed with High Sensitivity D1000 ScreenTape using a 4200 TapeStation System (G2991BA, Agilent, Santa Clara, CA). The cDNA libraries were sequenced on NovaSeq 6000 System (R1 – 28, i7 8, R2 – 91 cycles, Illumina, San Diego, CA).

2.4 scRNA-seq data analysis

Pre-processing of scRNA-seq was performed using Cell Ranger (v6.0). This pipeline included sample de-multiplexing, barcode processing, and 3’ gene counting. Reads were aligned to the GRCh38 transcriptome (v. 2020-A) and single cells within each sample were merged using the aggregate function in Cell Ranger. A total of 11,663 single cells with an average of 92,500 raw reads per cell were further analyzed in R (v4.0.3) and R studio (v1.4.1103) using the Seurat package (v4.0.1) (44). The datasets are available at GEO (GSE219210). A full collection of the pipeline that replicates the data analysis outlined in this manuscript is available as an R script repository on Github (https://github.com/dgomjim/Code_scRNAseq_myeloid_cells_in_HPV_TC).

We used Seurat to perform thresholding, normalization, integration, dimensionality reduction, cluster identification, and visualization as well as differential gene expression analysis. The gene-barcode matrix was filtered selecting cells with more than 700 genes and less than 6,000 genes detected per cell. Cells with more than 10% mitochondrial and more than 50% ribosomal UMI-counts were filtered out. Successful removal of doublets was ensured by comparing filtered cells to the barcodes selected using the DoubletFinder package (v2.0.3). Contaminating cells were identified and filtered out based on their expression of T-cell, B-cell, and NK-cell canonical genes (e.g., CD3D < 0.1 & CD3E < 0.1 & CD8A < 0.1 & CD8B < 0.1 & TIGIT < 0.1 & CD19 < 0.1 & MS4A1 < 0.1 & KLRB1 < 0.1 & CD7 < 0.1 & GNLY < 0.1 & CTLA4 < 0.1). Additionally, mitochondrial, and ribosomal genes were excluded before proceeding to data merging and normalization by log-transformation. Cell cycle scores were calculated using the “CellCycleScoring” function within Seurat, which generates a module score based on pre-set gene-sets characteristic of S, G2, and M phases of the cell cycle. Datasets from different samples were merged, and sample specific IDs were added to cellular barcodes allowing sample-based discrimination of cells during downstream steps. A total of 9,505 single cells were further selected for sample integration, which was performed via canonical correlation analysis (CCA) using the top 2,000 highly variable features across the datasets. The gene-barcode matrix was then scaled, and regression of gene expression was conducted based on UMI-counts and total number of unique UMIs as well as the percentage of mitochondrial and ribosomal genes. Principal component Analysis (PCA) was then performed on the normalized and scaled gene-barcode matrix. The first 20 principal components were selected for non-linear dimensionality reduction (tSNE and UMAP) based on their variance using the Elbow Plot approach. The shared nearest neighbors (SNN) algorithm was used to calculate distances between cells. Cell clusters were calculated using the community detection algorithm smart local moving (SLM). We inspected the granularity of the clusters, using step-wise increments of 0.1 in the resolution parameter (0.1-1) implemented in the “FindClusters” function, and plotted them using the ClustTree package (v0.4.3) (45). Differential gene expression analysis across cell clusters was performed using the “FindAllMarkers” function in Seurat. We applied Wilcoxon rank sum test to perform the analysis, only accounting for genes overexpressed in at least 10% of the cells in a cluster (p-value < 0.01, log₂FC > 1).

2.5 Gene signature analysis

We computed signature scores using the “AddModuleScore” function in Seurat to assess cluster identity. Scores were calculated by binning signature genes in 25 bins according to average expression. The signature’s average expression was corrected by subtracting the aggregated expression of 100 randomly selected genes from the same bin. We used published gene signatures for blood cDC1, cDC2, DC3, DC4, DC5, DC6, and monocytes 1-4 (30), in vitro activated cDC2 and moDC (26), M1/M2 in vitro polarized macrophages (37), skin LC (39), monocytes and monocyte-derived DC and macrophages (11, 35), myeloid derived suppressor cells (MDSC) (46), and tissue mast cells (35). Additionally, enrichment scores from a pre-annotated bone marrow single cell dataset (47) were computed across our myeloid clusters in order to identify progenitor DCs. The CITE-seq bone marrow dataset was accessed through the human cell atlas (48), and a supervised PCA was computed. Supervised PCA was used to map our query myeloid dataset onto the multimodal bone marrow reference, obtaining enrichment scores for each annotation.

2.6 Single-cell regulatory network inference and clustering analysis

We used single-cell regulatory network inference and clustering (SCENIC) to detect and score regulons in single cells (49). The pipeline consisted of three steps: generation of a co-expression modules using Genie3, module filtering based on enrichment of DNA binding sites using RcisTarget databases, and estimation of regulon activity with AUCell. The relative activity of the top ten enriched regulons per cluster was used to generate a heatmap. The regulon heatmap was complemented with a dot plot displaying the relative gene expression of the corresponding TF.

2.7 Developmental pathway analysis

We used the Velocyto python package to recount the spliced and unspliced reads. Selection of highly variable genes was performed on total reads, followed by a filtering step to select genes with more than 20 counts on both spliced and unspliced assays. After filtering, top 400-500 highly variable genes were used to calculate RNA velocity applying the scVelo “stochastic” model in the R package Velociraptor (50). Finally, we embedded the velocity vector in either a diffusion map (51) or UMAP to visualize developmental trajectories. The frequency of maturing DCs was calculated as the percentage of DCs with high RNA velocity values transitioning into LAMP3+ expressing actDC, out of the total number of DC in each lineage.

2.8 Survival analysis

TCGA-HNSC normalized gene expression data as well as clinical metadata were downloaded from the GDC portal repository of the cancer genome atlas (TCGA) (52). Individual samples were normalized and z-scored for cluster-specific gene-sets using the GSVA package (53). Then, samples were divided into quartiles, and we defined high (Q1) and low (Q4) enriched samples for each cluster-specific gene set. We performed five-year overall survival analysis using the survival and survminer packages (54). Difference between survival curves were assessed by log-rank test (p-value < 0.05). We used uni- and multivariate Cox proportional-hazards models to estimate the relation between the gene-signature scores and the mortality rate. The assumptions of the Cox proportional-hazards model were tested assessing the Schoenfeld residuals.

2.9 Pathway analysis

Gene set enrichment analysis (GSEA) was performed using EnrichR (v3.0) (55) by querying differentially expressed genes among myeloid clusters onto the gene ontology database (GO, GO_bioprocess_2018) (56). We selected reoccurring and highly ranked pathways (adj.p-value < 0.01) related to APC function. We harmonized these pathways throughout clusters using gene set variation transformation (53). The resulting GSVA scores were plotted as a clustered heatmap using pheatmap (v1.0.12).

2.10 RL interaction analysis

To elucidate potential RL interactions of myeloid cells, we built a “meta-dataset” containing our in-house dataset and tumor infiltrating leukocytes from TC and HT samples of a publicly available HNC scRNA-seq dataset (41). In brief, we subsetted single cells from three TCs and five HT samples and performed normalization, sample integration via CCA, data scaling, dimensionality reduction, cluster detection, and annotation based on canonical markers. Next, we integrated leukocytes of the external and in-house datasets, normalized them via “SCTransform”, and re-calculated PCA and UMAP. Once we obtained the harmonized “meta-dataset”, we proceeded to RL analysis using the R package CellChat (57). First, we identified RL interactions enriched in specific clusters in TC and HT separately. To do so we split the “meta-dataset” into cancer and healthy, and selected overexpressed genes detected in at least 25 of the cells per cluster (p-value < 0.05, FC > 2) before calculating the communication probability. Next, we merged the cancer and healthy CellChat objects and filtered cluster specific RL interactions that were enriched in cancer. Finally, we visualized cluster specific RL interactions of cDC1s, actDCs, and act Macros with subsets of T/NK-cells using the “netVisual_bubble” function in the CellChat package.

3 Results

3.1 Expansion of the myeloid compartment in HPV⁺ TC include subsets of DCs, monocytes, macrophages, and a small population of HLA-DR⁺ granulocytes

To evaluate the heterogeneity of the myeloid compartment in HPV⁺ TC, we dissociated fresh biopsies into single cells and used flow cytometry to characterize and quantify CD45⁺CD13⁺HLA-DR⁺ myeloid cell populations (Figure S1A). We observed a significant increase in the frequency of CD13⁺HLA-DR⁺ cells in TC (n=8) as compared to contralateral HT (n=5), indicating that the myeloid compartment expands in TC (Figure S1B). Within the myeloid gate, the greatest frequency increase in TC as compared to HT was observed in the Mono-Mac and LC populations, followed by cDC1s. Next, we addressed the diversity of myeloid cells in treatment-naïve TC patients at transcriptomic level. We sorted CD45⁺CD13⁺HLA-DR⁺ cells from five HPV⁺ TC biopsies and one paired contralateral HT, followed by scRNA-seq using the droplet-based 10X Genomics platform (Figure 1A). After QC, normalization, and sample integration, we used UMAP to visualize single cells in a low dimensional space. Unsupervised clustering of 9,502 single cells revealed 12 clusters (Figure 1B; Figure S1C), characterized by expression of MHC class II coding transcripts (e.g., HLA-DRA, HLA-DRB1) and heterogeneous expression of CD14 (Figure 1C), with different distribution in TC and HT (Figure 1D).

FIGURE 1

Figure 1 Heterogeneity of myeloid populations in TC and HT assessed by scRNAseq. (A) Viable CD45⁺ CD13⁺ HLA-DR⁺ myeloid cells were sorted from single cell suspensions of five TC samples and one contralateral HT and subjected to scRNA-seq using 10X genomics. The frequencies in the dot plots represent the average in TC calculated from viable CD45⁺ cells. (B) UMAP visualization of single cell transcriptomes displaying twelve clusters grouped into six major myeloid subtypes. (C) UMAP visualization, as in B, showing expression levels of HLA-DRA, HLA-DRB1, and CD14. (D) Frequency of myeloid populations across tissue of origin. (E) Dot plot displaying three canonical markers among the top ten differentially expressed genes across clusters. (F) Violin plots showing scores of indicated signatures. The score in each cell is shown in relative units, along with the density distribution shown by the shape of the plot. G2M and S scores represent the relative expression of cell cycle phase genes in a cell. Scale bars represent normalized counts in (C) and (E).

Next, we annotated the twelve clusters using a combination of canonical marker gene expression, after differential gene expression analysis (DGEA) across clusters and gene signature scoring (11, 26, 30, 47) (Figures 1E, F). The cDC1 population (cluster 3) expressed high levels of CLEC9A, IRF8, and XCR1, along with a high score of a described blood cDC1 gene-set. Similar to that observed at protein level (Figure S1B), the cDC1 frequency in the scRNA-seq dataset increased in TC as compared to paired HT. In turn, cDC2-like cells (clusters 0 and 6) displayed CLEC10A, CD1C, and FCER1A expression and enrichment in a blood cDC2 gene set. Two additional DC clusters were predicted, representing different cell states of DC development and maturation. Progenitor DCs (cluster 9) exhibited markers related to cell cycle and progenitor cells (MKI67, TOP2A, and LTB), as illustrated by their high S and G2/M cell cycle phase scores. Additionally, this cluster featured high prediction scores of bone-marrow DC progenitor cells as well as low scores of hematopoietic stem cells and granulocyte-monocyte progenitors (Figure S1D). The 4^th DC population, activated DC (actDC, cluster 7), was characterized by high expression of genes related to DC maturation, such as CCR7, LAMP3, and CCL19, and high scores for a gene-set from activated DCs. Six clusters of CD14⁺ cells (clusters 1, 2, 4, 5, 8, and 10) conformed to the Mono-Mac lineage (Figures 1E, F). The Mono-Mac population was characterized by expression of canonical genes S100A8/9 and C5AR1 as well as a high CD14⁺ monocyte signature score. As observed in the flow cytometry analysis (Figure S1B), Mono-Macs were more frequent in TC as compared to HT (Figure 1D). Lastly, cluster 11 showed high expression of KIT, GATA2, and TPSB2 as well as an exclusive mast cell signature score, which only was detected in TC samples. Mast cells expressed lower levels of MHC class II related transcripts compared to other myeloid cells, and expression of co-stimulatory molecule coding genes was not detected (Figure 1C; Figures S1E).

3.2 DC lineage-specific subpopulations and analysis of cell identity dynamics.

To detect rare populations, we divided the dataset into DCs and Mono-Macs and re-clustered these lineages separately. Re-clustering followed by DGEA revealed seven transcriptomically different communities in the DC lineages across TC and HT (Figures 2A-C; Table S4). Three clusters of cDC2-like cells were characterized by expression of CLEC10A, CD1C, and FCER1A (Figure 2C), high scoring of blood and tissue cDC2 gene-sets (Figure S2A), and heterogeneous distribution in TC and HT (Figure 2B). We identified two clusters of cDC2s with different expression of the transcriptional factor RUNX3. RUNX3^- cDC2s exhibited preferential expression of CLEC10A, JAML, and TXNIP and were only present in HT. In contrast, RUNX3⁺ cDC2s were characterized by expression of interferon-inducible genes (IFI30, IFITM3), and were more abundant in TC as compared to HT. Additionally, MAFB⁺ LCs uniquely expressed CD1A, TGFB, CD14, and C1QA/B/C and were only found in TC. This population was enriched in a LC and moDC gene-sets, but not in a previously described blood CD14⁺ DC3 gene-signature (Figure S2A). MAFB⁺ LC displayed high G2M scores, highlighting their self-renewal capacity and further confirming their identity as LC (38) (Figure S2C). A rare population of DC5s expressing SIGLEC6, LILRA4, and IL3RA was found closely associated with the cDC2-like cell cluster. DC5s featured a high blood DC5 gene-set score (Figure S2A) and an increased frequency in HT as compared to TC. Lastly, cDC1s, actDCs, and progenitor DCs corresponded to the three remaining DC populations, equally observed in the previous clustering (Figure 1B). The correspondence between the two clustering steps showed the homogeneity of cDC1, actDC, and progenitor DC communities at the present sequencing depth.

FIGURE 2

Figure 2 Dendritic cell dynamics and activation differ according to the tissue of origin. (A) tSNE plot showing seven distinct DC clusters DC clusters in TC and HT, or (B) separated by tissue of origin. (C) Heatmap displaying top ten DEG per cluster after DGEA (FC > 1.5, p-value < 0.05). (D) Heatmap of the top 10 regulons per cluster colored by regulon activity (upper) and dot plot displaying the gene expression levels of each regulon’s TF per cluster (lower). (E) Dendrogram of DC clusters based on the Euclidean distance. Diffusion map combined with RNA velocity analysis to infer the development trajectory of all DCs (F), cDC1s (G), and DC2-like cells (H). (I) Frequency of maturating cDC1s and DC2-like cells annotated by sample of origin, as inferred from RNA velocity analysis. DEG, differentially expressed genes; DGEA, differential gene expression analysis; DC, Dendritic cell; TF, Transcription factor.

We used SCENIC to assess different activity of gene regulatory networks in our clusters (Figure 2D; Table S5). We noted that cDC1s progenitor DCs and actDCs featured the most distinct regulons in the DC compartment, while cDC2-like cells and DC5s partly shared their regulatory networks. actDCs were characterized by high activity of RELB, HIVEP1, and IRF1/2 regulons. In turn, cDC1s featured high activity, and expression of STAT2, RUNX1, and IRF8 regulons. Progenitor DCs uniquely displayed activity and expression of E2F1/7/8 as well as BRCA1 and MYBL2. DC5s uniquely showed TCF4 and RUNX2 regulon activity and expression, and preferential RARA enrichment as compared to MAFB⁺ LC. On the other hand, MAFB⁺ LC exhibited exclusive MAFB activity and expression, and preferential enrichment in RXRA and DAB2 among other regulons and TFs. Finally, RUNX3^-, RUNX3⁺ cDC2s, and MAFB⁺ LC shared several active regulons including SNAI3, CEBPD, and TFEB. However, while RUNX3⁺ cDC2s exhibited higher activity of VDR, RUNX3^- cDC2s were characterized by higher NFYB and SMAD3 activity. Next, we performed hierarchical clustering based on Euclidean distance between the clusters to assess their transcriptomic similarity (Figure 2E). We observed that actDCs were the most divergent group, clustering separately in the dendrogram. cDC2-like cells and DC5s clustered together and were distant from the cDC1 and progenitor DC pair, in line with the results from our gene regulatory network analysis.

Next, we sought to study cell identity dynamics of DC communities. Thus, we projected single cell transcriptomes along with RNA velocity information onto diffusion map embeddings to infer future states of cells (Figures 2F–H). Overall, we observed directionality in RNA velocity of both cDC1 and cDC2-like cell branches towards actDCs. cDC1s and cDC2-like maturation into actDCs was confirmed by the sequential loss of CLEC9A and CLEC10A, respectively. These events were followed by a sequential increase in LAMP3 expression along the trajectory, as observed in a 3D UMAP embedding (movie S1). The maturation of the cDC1 lineage was investigated by subdividing cDC1s along with progenitors enriched in cDC1 gene-set and actDC. RNA velocity analysis displayed a continuum of states along the progenitors into the fully differentiated cDC1 cluster, progressing to actDCs (Figure 2E). We evaluated transcriptomic changes along the developmental trajectories by assessing the velocity of highly variable genes with cluster-specific splicing dynamics (Figure S2B). We observed that progenitor DCs had higher unspliced levels of DNAJB1, FUCA1 and DNMT1, on their development towards cDC1s. In turn, SGO2 and CTSL were downregulated on the progression of progenitor DCs towards cDC1s and actDCs. cDC1s matured into actDCs displaying an initial upregulation of CD40, followed by CCR7, PRDM1, NFKB1, LAMP3, IDO1, and MIR155HG. In parallel, we analyzed cDC2-like cells likewise and detected a strong directionality from RUNX3^- cDC2s and MAFB⁺ LC into RUNX3⁺ cDC2s, which later progressed into actDCs (Figure 2H). We observed that all DC2-like clusters presented active transcription of CD1C, IRF4, and HIVEP1, and the last two genes were markedly upregulated upon maturation into actDCs (Figure S2C). MAFB⁺ LC and RUNX3⁺ cDC2s exhibited upregulation of RUNX3, while only RUNX3^- cDC2s upregulated ISG20 upon progression into RUNX3⁺ cDC2s. MAFB⁺ LCs specifically downregulated CD1A and CD207 upon progression into RUNX3⁺ cDC2s. DC5s displayed active transcription of TCF4, while LILRA4 was downregulated. In turn RUNX3⁺ cDC2s sequentially upregulated CD274, LAMP3, and TCF4, followed by IRF1, upon maturation into actDCs. Contrary to progenitor cDC1s, we did not detect a directional flow in DC5s nor in CLEC10A⁺ progenitor DCs. Finally, we used RNA velocity information to estimate the frequency of cDC1 and DC2-like cells maturing into actDC expressing LAMP3 (Figure 2I). We observed that DC2-like cells matured more frequently than cDC1s. We also noted that cDC1 and DC2-like maturation in the same TC samples did not correlate to each other, indicating that these two processes are independent.

3.3 Mono-Mac subsets display tumor-stage specific patterns

Similar to DC lineages, we subsetted Mono-Macs and re-clustered this lineage to address its heterogeneity. Using the SLM clustering algorithm, followed by DGEA, we detected and characterized six transcriptomically unique clusters of Mono-Macs, with different distributions based on TNM classification of TCs (Figures 3A-C; Figures S3A, B and Table S6). Monocytes, expressing high levels of FCN1, CTSS, and JUNB, were highly enriched in blood CD14⁺ and CD16⁺ monocyte as well as MDSC signatures. In turn, ISG Monos featured high levels of IFN inducible genes (IFIT1, IFIT3, IFITM), and a high score of an ISG Mono gene-set. NLRP3 Macros expressed high levels of NLRP3, CD300E, and AREG, and were preferentially enriched in an NLRP3 and MDSC gene-sets. C1Q Macros expressed C1QA/B/C genes similarly to MAFB⁺ LC, but uniquely expressed FCGR3A and were enriched in a C1Q tumor-associated macrophage (TAM) gene-set. A cluster displaying expression of IL7R, CCR7, MIR155HG, and CD274 was annotated as activated macrophages (act Macro). Finally, CXCL Macros showed high expression of CXCL1/3/5, SPP1, and MMP1 as well as high score of a SPP1 TAM gene-set. Remarkably, signature scores corresponding to M1 and M2 activation states (37) did not distinguished the Mono-Mac clusters, indicating that the M1/M2 dichotomy did not resolve the differences between these (Figure S3A). Apart from CXCL Macros, all Mono-Mac clusters were found in HT and TC, regardless of tumor stage. Instead, CXCL Macros were only detected in advanced primary tumors (n=2, T3N1M0, Figure 3B).

FIGURE 3

Figure 3 Mono-Mac subsets display tumor-stage specific patterns. (A) tSNE plot showing the six transcriptomically unique Mono-Mac clusters across cancer and healthy tissues. (B) tSNE plot of Mono-Mac distribution according to TNM classification. (C) Heatmap of the top 10 DEG per cluster after DGEA (FC > 1.5, p-value < 0.05). (D) Heatmap of the top 10 regulons per cluster colored by regulon activity (upper) and dot plot displaying the gene expression levels of each regulon’s TF per cluster (lower). (E) Dendogram of Mono-Mac clusters based on Euclidean distance. (F) UMAP plot showing the inferred development trajectory of Mono-Macs using RNA velocity. (G) DEG between CXCL and act Macro clusters, and enriched GO pathways obtained after GSEA. DEG, differentially expressed genes; DGEA, differential gene expression analysis; GSEA, gene set enrichment analysis.

We evaluated gene regulatory network activity and transcription factor expression in Mono-Mac clusters (Figure 3D; Table S7). Briefly, the Mono-Mac lineage featured common expression of FOS and FOSB, where monocytes presented the higher regulon activity. In addition, Monocytes and NLRP3 Macros shared preferential expression and regulon activity of HES1, RARA, RXRA, and NFIL. In comparison, ISG Monocytes displayed expression and high regulon activity of IRF7/8 and STAT1, and shared RUNX3 and MEF2A regulon activity with C1Q Macros. Finally, CXCL and act Macros shared expression and activity of certain regulons including STAT4 and HIVP1, but while the former displayed preferential NR3C1 and TAF7 activity and expression, the latter featured higher STAT5A, NFKB1, and REL activity. We compared the transcriptomic similarity between the Mono-Mac clusters by performing hierarchical clustering based on Euclidean distance (Figure 3E). Monocytes and ISG monocytes, NLRP3 and C1Q Macros, and CXCL and act Macros, were paired together, respectively.

Next, we projected single cell transcriptomes along with RNA velocity information onto UMAP embeddings to decipher future state of cells (Figure 3F). We observed that Monocytes and ISG mono differentiated into NLRP3 and C1Q Macros, which in turn polarized into CXCL and act Macros, respectively. Furthermore, we evaluated transcriptomic changes across the Mono-Mac developmental trajectories by analyzing highly variable genes with dynamic changes in their velocity (Figure S3C). We noted that Monocytes, ISG Monos, NLRP3, and C1Q Macros upregulated the expression of IL1B and OLR1, while the expression of these genes was drastically downregulated upon activation into CXCL and act Macros. Monocytes and NLRP3 Macros featured active transcription of NLRP3, S100A9, AQP9, and FCN1. Upon progression into CXCL Macros these upregulated MMP10 and CXCL5. On the contrary, ISG Monos presented active transcription of ISG20 and IDO1, and sequentially upregulated HLA-DPB1, IL4I1, and IRF4 as they differentiated into C1Q and act Macros, respectively. Together these results highlighted that there was a core transcriptional program in different stages of the Mono-Mac lineage trajectory.

To gain insights into the two parallel Mono-Mac differentiation processes, we performed pairwise DGEA between act Macro and CXCL Macro clusters followed by GSEA (Figure 3G). act Macros featured a higher expression of antigen presentation related transcripts and enrichment of the corresponding pathways. In turn, CXCL Macros were characterized by high expression of chemotactic genes related to granulocyte and neutrophil chemotaxis (CXCL1/3/5/8), angiogenesis (VEGFA), and genes related to matrix remodeling pathways (MMP1, SPP1). These results highlighted distinct functions of the terminal macrophage states in the two differentiation pathways present in our dataset.

3.4 Validation of myeloid heterogeneity in TC and paired HT

Next, we sought to define a flow cytometry antibody panel to validate the myeloid communities detected by scRNA-seq at protein level. Protein targets were selected by performing DGEA analysis across myeloid cell clusters, followed by a filtering step to select cell membrane protein coding transcripts (Table S8). We evaluated TC and HT-derived single-cell suspensions using an antibody panel directed towards proteins encoded by genes identified in our scRNA-seq analysis (XCR1, CD1c, CD14, CD5, CCR7, CD1A, CD207, SDC2, HLA-DRA, CD300E, FCGR3A, and TNFSF10) (Figure 4A). The panel was complemented with antibodies specific for immune-checkpoints PD-L1, PD-L2, LAG3, and BTLA as well as the co-stimulatory protein CD40, given that their corresponding genes exhibited cluster/lineage specific expression in our scRNA-seq dataset (Figure 4B). Antibodies against CD19 and CD20 targeting B-cells, and CD123 targeting plasmacytoid DCs, were included in the panel to increase resolution and compare myeloid cells to other APC populations. We also included other markers in our panel that, despite not being detected in our analysis, have recently been reported to mark novel myeloid populations. CD163 has been reported to mark a CD14⁺ DC3 population possibly distinct from CD14⁺ LC and moDC (30, 31). CCR2 was included in the panel since it is downregulated by monocytes upon differentiation (36, 58), while CD34 was included to detect cDC progenitors (30).

FIGURE 4

Figure 4 Validation of scRNA-seq clusters using 26-plex flow cytometry. (A) Gene expression of cluster-specific markers mapping to cell membrane from scRNA-seq analysis. (B) Expression of genes related to APC function from scRNA-seq analysis. (C) tSNE plots showing individual marker expression on 45,712 viable CD45⁺HLA-DR⁺CD19^-CD20^- single cells from TC (n=7) and paired HT (n=3). (D) tSNE plots showing validated populations and their distribution in TC and HT tissues. (E) Gating strategy used to identify scRNA-seq clusters using flow cytometry. (F) Expression of genes related to APC function from flow cytometry. DC, Dendritic cell; TC, Tonsillar cancer; HT, Healthy tonsil; pDC, plasmacytoid DC.

We analyzed ten cryopreserved single cell suspensions, corresponding to biopsies from seven TCs and three paired HTs, by flow cytometry (Table S1). After singlet filtering, viable CD45⁺HLA-DR⁺CD19^-CD20^- cells were concatenated and visualized within a tSNE space. We observed a clear signal from all markers apart from CD16, LAG3, TRAIL, BTLA, CD34, and CD300E (Figure 4C). We established a gating strategy to define several myeloid populations in human TC and HT and compared their levels of PD-L1, PD-L2, HLA-DR, CCR7, and CD40 (Figures 4D-F). We obtained a clear separation of the dataset based on CD11c and CD123, with minimal overlap between the two markers. Careful inspection of the CD123⁺ gate revealed two small CD5⁺ populations corresponding to DC5s, which were distinguishable by CD11c expression, as previously reported in blood (30). Compared to pDCs, DC5s featured higher levels of HLA-DR, CD40, and CCR7. In turn, inspection of the CD11c⁺CD14⁺CD1c^- Mono-Mac gate revealed a parallel increase of HLA-DR along the SDC2 signal. SDC2⁺ HLA-DR^high corresponded to act Macros according to their co-expression of PD-L1, CD40, and CCR7 at protein level, as observed at RNA level. In addition, this population exclusively expressed CD163 compared to other CD14⁺ Mono-Mac populations. Next, we focused on resolving the heterogeneity of the CD1c⁺ gates. Similar to our scRNA-seq analysis, we observed co-expression of CD1c and CD207 in CD14^+/- cDC2s, suggesting that these corresponded to our newly identified RUNX3⁺ cDC2 cluster. In turn, CD1a⁺CD207⁺CD14^+/- was identified as LC. Interestingly, we observed a co-expression of CD163 and CD14 at protein level in CD1c⁺CD1A^-CD207^- cells. This indicates that DC3s are present in TC, as previously reported in blood (30, 31) and other HNCs (40), but are indistinguishable from other CD1c⁺ cells by scRNA-seq at our sequencing depth. In turn, CD5 expression inversely correlated with CD14 in CD1c⁺ cells, suggesting that these were bona fide cDC2s, as described in blood (31). While all CD1c⁺ populations had some expression of PD-L1, PD-L2, HLA-DR, CCR7, and CD40, CD207^-CD1A^-CD1c⁺ cDC2s showed the most immature profile according to CD40 and CCR7 levels. Finally, XCR1⁺ cDC1s were clearly distinguishable by their high CD40 and HLA-DR levels, CCR7 expression, and residual PD-L1 and PD-L2 expression, indicating a mature profile of these population.

To assess transcriptomic changes induced by the TME, we performed DGEA comparing TC and HT myeloid populations from our scRNA-seq dataset (Figure S4, Table S9). Progenitor DCs in TC expressed higher levels of cDC1 canonical genes (CLEC9A, BATF3, XCR1) and cell cycle genes (CCND1, CAMK2D), as illustrated by the enrichment in blood cDC1 and cell cycle scores, respectively. In contrast, progenitor DCs in HT featured higher levels of canonical cDC2 genes (CLEC10A, CD1C) and blood cDC2 gene-set score. On comparing DEG between TC and HT, we observed a common increase of IFN response-related genes (ISG15, ISG20, TXN) in all DC clusters. Conversely, all myeloid populations in TC showed a downregulation of inflammasome-regulation genes MTRNR2L12, DNASE1L3, and LGALS2 as well as MHC-II coding genes.

3.5 cDC1s and activated DCs and macrophages correlate with increased survival in HNC patients

We sought to assess the presence of myeloid populations described here in TC and HT, in other HNCs. To this end, we accessed publicly available scRNA-seq data (41), subset myeloid cells and integrated them using our in-house dataset as reference. We observed that all clusters identified in our in-house TC dataset were also present in other subtypes of HNC, although in varying frequencies (Figure S5A). RUNX3^- cDC2s were only marginally present in HNC subtypes (4 predicted cells), indicating that this cell state was specific to healthy tissue as described above (Figure 2A) thus validating our prediction. Surprisingly, the diversity of tumor-infiltrating myeloid cells was similar between HPV⁺ and HPV^- tumors (Figure S5B).

To investigate the prognostic impact of myeloid populations in HNC, we defined gene signatures that enabled us to infer the relative abundance of myeloid populations within bulk HNC RNA-seq. We performed DGEA across all myeloid clusters that were isolated from TC, obtaining a total of 998 DEG (FC ≥ 2, p-value ≤ 0.01, Table S8). Then, we filtered our candidate genes using a publicly available scRNA-seq HNC dataset (41), to select cluster-specific genes that were not expressed by other tumor infiltrating-leukocytes (Figures 5A, B). We obtained 3-gene transcriptomic signatures for all our clusters except progenitor DCs, monocytes, and C1Q Macros. We further investigated the impact of these 3-gene signatures within bulk RNA-seq transcriptomes of HNC from the TCGA (52). We applied Cox proportional-hazards regression to compare the survival rate of high (Q4) and low (Q1) scoring quartiles of each gene signature (Figure 5C, Table 1). We found that high scores of cDC1, actDC, and act Macro gene-sets correlated with increased five-year survival in HNC, independently of age and sex. Of note, high score of the cDC1 signature was also positively correlated to survival when TC patients were assessed separately (Figure S5C).

FIGURE 5

Figure 5 Impact of myeloid populations on patient survival, as well as functions and interactions thereof. (A) Dot plot displaying cluster-specific signature genes in the myeloid dataset and (B) in a publicly available scRNA-seq dataset (41). (C) Survival curves showing five-year overall survival rates for the indicated gene-sets, from a publicly available bulk RNA-seq dataset of HNC patients (n=529) (52). Comparison of survival probability was performed using log-rank test and Cox proportional-hazards model. Years: years since diagnosis. (D) Heatmap colored by GSVA score of enriched GO pathways in the myeloid clusters (adjusted p-value < 0.05). (E) Chord diagrams displaying the interaction between cell pairs. Interaction strength is presented as the weight of the edges, and the receiver cell cluster is indicated by the apex of the arrow. (F) Bubble plot showing statistically significant receptor-ligand (R-L) coding gene pairs in TC between myeloid cell subsets and T/NK- cells from a publicly available dataset (41). GSVA, gene set variation analysis; FC, fold change; p, p-value; CM, central memory; EM, effector memory; Exh, exhausted; IFN-act, IFN-activated; NK, natural killer; TFH, T follicular helper; Transi, transitional; TRM, tissue-resident memory.

TABLE 1

Table 1 Adjusted hazard ratios (HR) comparing the five-year overall survival rate between HNC patients with high (Q4) and low (Q1) cluster-specific gene-set scores.

Next, we aimed at elucidating mechanisms by which cDC1s, actDCs, and act Macros impacted survival. Pathway enrichment analysis (GO bioprocesses 2018 (56), adjusted p-value < 0.05) revealed subset specific functions within the myeloid compartment in TC (Table S10). Highly scoring pathways related to myeloid cell function, such as antigen-presentation and chemotaxis, were subjected to GSVA score transformation and compared across clusters (Figure 5D). Overall, we found that all DCs, but also ISG monos, featured higher scores for antigen presentation pathways compared to macrophages. Specifically, the cDC1 cluster displayed the highest score in pathways related to cross-presentation via MHC-I, followed by ISG monos, DC5s, and actDCs. In contrast, all DC clusters presented similar scores in presentation via MHC-II, which was higher compared to Mono-Mac clusters. In turn, macrophage clusters were characterized by higher expression of genes related to chemotaxis, phagocytosis, ECM disassembly, angiogenesis, and complement pathway. Notably, ISG monocytes were preferentially enriched in cellular response to IFNG and negative regulation of viral life cycle, suggesting a TME-specific polarization of this cluster.

To investigate which leukocytes that interacted with cDC1s, actDCs, and act Macros, we accessed a publicly available HNC scRNA-seq dataset (41). After annotating single cells from TC and HT based on canonical marker expression, we integrated all leukocytes of the external dataset and our in-house dataset (Figure S5D). Following normalization, we performed receptor ligand analysis to identify major cell-cell interactions in TC as well as receptor-ligand gene pairs overexpressed in TC as compared to HT (Table S11). We focused our analysis on the interaction of myeloid cells with T/NK-cell subsets, because of the high abundance of the latter in the TME (Figures 5E, F). cDC1s featured higher signaling strength with NK and tissue resident memory (TRM) CD8 T-cells, followed by Treg subsets, as compared to other T/NK-cells. Furthermore, cDC1s interacted with TRM CD8 T-cells, effector-memory (EM) CD8 T-cells, and NK-cells via XCL1/XCL2-XCR1 as well as CLEC2B-KLRB1 axes. Based on the predicted signaling strength, actDCs and act Macros interactions were more diverse than for cDC1s. Both actDCs and act Macros were predicted to interact with several T/NK-subsets via NECTIN2-TIGIT. However, while actDCs were uniquely predicted to interact with T/NK-subsets via CCL19-CCR7, act Macros displayed a high probability of interacting via ANXA1-FPR1. The interaction via CD40-CD40LG axis was predicted in several myeloid and T/NK-subset combinations, but it was highest between actDC and T-follicular helper (TFH), naïve, and transitional CD4 T-cells. Finally, the common predicted interactions of cDC1, actDC, and act Macro clusters in TC included: IFNG-IFNGR1/2 with NK and EM CD8 T-cells, as well as CD86-CTLA4 with Tregs, exhausted CD8 T-cells, and TFH cells.

4 Discussion

HPV-driven TC is a subset of HNC characterized by type 1 T-cell inflammatory responses (40). Effective cellular immune responses are complex phenomena involving an interplay of different cell types. Myeloid cells represent a particularly interesting entity due to their marked ability to prime T-cells and their wide diversity (30). Motivated by the expansion of the CD13⁺HLA-DR⁺ population in HPV⁺ TC compared to paired HT, we sought to unbiasedly assess their heterogeneity. Here, we present a high-resolution characterization of TC infiltrating and lymphoid tissue resident myeloid cells. We generated a scRNA-seq dataset of 9,505 CD45⁺CD13⁺HLA-DR⁺ myeloid cells from TC biopsies and paired HT from treatment naïve patients using 10X Genomics. This dataset represents a tool for understanding myeloid identity dynamics in the tumor as well as steady state lymphoid tissue settings. By combining transcriptome profiling, gene signature scoring, regulon activity assessment, and RNA velocity, we obtained a dynamic snapshot that illustrates the interrelationship of different myeloid communities. In addition, pathway enrichment, R-L, and survival analyses of in-house and publicly available HNC datasets yielded information on how myeloid communities impact survival and on the underlying molecular mechanisms and cell-cell interactions.

DC identity has recently been re-assessed in blood (30), lymphoid organs (32, 33), and several types of cancer (11, 28, 35). Here we identified seven DC populations corresponding to four different lineages and diverse cell states. We observed that the cDC1 cluster was consistently homogeneous through tissue types, and that its identity was dictated by the IRF8 and STAT2 regulons, as described previously (59, 60). Interestingly, cDC1s also featured high regulon activity and expression of SMARCA5. This TF has recently been characterized in viral DNA sensing and induction of IFN pathways, enabling DCs to prime CTL responses with a marked IFN-γ and IL-12 production (61, 62). In this context, our in-house generated cDC1 signature highlights a positive impact of cDC1s on HNC patient survival, and for TC patients when assessed separately. Pathway enrichment analysis evidenced a strong expression of genes related to MHC-I-mediated cross-presentation, which is a well characterized trait of the cDC1 lineage (16). Complementing the hypothesis that cDC1s migrate in the tumor through the XCL1/2 – XCR1 axis (15), we show that, besides for NK-cells, XCL1/XCL2 were also expressed by a subpopulation of ENTPD1⁺ ITGAE⁺ CD8⁺ T-cells. In fact, the interaction strength between these CD8⁺ T-cells and cDC1s was also mediated by the CLEC2B-KLRB1 axis. All in all, our data suggest a key interaction between cDC1 and a small subpopulation of TRM CD8⁺ T-cells that has been reported to exert a tumor-reactive function in the TME of several cancer types including HNC, but also ovarian, lung, and colon cancer (63). Interestingly, RNA velocity analysis indicated that the cDC1 population stems from a cycling progenitor DC-population to later mature into CCR7⁺ LAMP3⁺ actDC, suggesting that cDC1s acquire a migratory potential in the TME. By comparing the signature scores of the progenitor DC population in TC and HT, we observed a marked increase in cell cycle genes and a preferential blood cDC1 signature enrichment in TC’s progenitor DCs. This indicates that progenitor DCs infiltrating TC are mainly pre-committed to a cDC lineage as is the case in blood (64). In addition, it suggests that, in the context of virally driven TC, there is a high production rate of cDC, and that these are pre-committed to cDC1, to the detriment of cDC2. This hypothesis is supported by the overall increase in myeloid APC frequency in TC, including cDC1s, while cDC2 populations were similarly abundant, as estimated both in scRNA-seq analysis and by flow cytometry. Brown et al., previously reported the presence of mitotic pre-cDCs in the spleen (33). To our knowledge, our study is the first report of cycling pre-cDCs in cancer, shifting the balance towards the cDC1 lineage.

Similar to what has been reported for blood (30) and spleen (33), DC2-like cells displayed higher heterogeneity compared to cDC1s both at RNA and protein level. Our dataset recapitulates transcriptomically distinct clusters with common CEBPB regulon activity (65). In scRNA-seq, MAFB⁺ LC acted as a source of DC2-like cells in TC, while RUNX3^- cDC2s were the unique source in HT, as seen in the RNA velocity analysis. MAFB⁺ LC uniquely exhibited MAFB and RXRA regulon activity, as well as expression of CD207, CD14, CD1A and CD1C at gene and protein level. Furthermore, these featured high scores for skin LC and tissue cDC2 gene sets (35), and no enrichment in a DC3 gene-signature, indicating that these cells did not correspond to the recently characterized DC3 subset in blood (30, 66). In fact, our flow cytometry results showed that LC and DC3 are two distinct entities judging by the lack of CD163 co-expression with CD207 and CD1a. In addition, some MAFB⁺ LC exhibited high scores of G2 and M cell cycle phases, indicating that these cells had self-renewal capacity in tonsillar tissue, as they do in skin (38, 67). MAFB⁺ LCs featured high RUNX3 upregulation, and downregulated the activity of the MAFB regulon on their differentiation into RUNX3⁺ cDC2s. In line with our results, previous studies have reported expression of MAFB and RUNX3 during LC differentiation in mouse oral mucosa and skin (68–70). In contrast, RUNX3^- cDC2s most likely represented bona-fide lymphoid organ resident cDC2s, as judged by their exclusive presence in HT, absence in TC and other HNC subtypes, and the lack of monocyte lineage markers and regulons at RNA and protein level. We observed that cDC2s upregulated RUNX3 along with CD207 and interferon response genes ISG15, IFI30, and IFITM3 in TC. We hypothesize that, in the TME, cDC2s are subjected to IFN-I dependent maturation, potentially via TLR stimulation as it occurs upon viral or poly I:C stimulation (71, 72). Consistent with our observation at RNA level, CD207⁺CD1a^- cDC2s featured a more mature profile at protein level as compared to their CD207^- counterparts. Upon maturation into actDCs, DC2-like cells downregulated the expression of RUNX3, which is a known repressor of CCR7 expression and cDC2 maturation (73). Collectively, our RNA-velocity analysis showed that both cDC2 and LC matured into actDC gaining CCR7 expression, as reported previously by Bennewies et al. and Reynolds et al. (14, 39). Compared to other myeloid populations, MAFB⁺ LC expressed high levels of TGFB, which has been shown to inhibit CCR7 expression by DC in vitro and cause intra-tumoral retention of act DCs in mouse models (74, 75). Survival analysis based on signature scoring of the TCGA dataset did not suggest a potential contribution to HNC or TC patient survival of either of the three DC2-like populations. However, DC2-like cells, as well as cDC1s, converged into an actDC gene-program governed by RELB (76) and IRF1/2 (77), which was associated with an increase in survival, as assessed by gene-signature scoring. Pathway enrichment analysis suggested that actDCs perform antigen presentation via MHC-I and II in the TME and interact with diverse T-cell populations. Interaction through the CCL19-CCR7 axis is indicative of active T-cell recruitment into TC’s TME (78, 79), while the CD40-CD40LG axis suggests a positive feedback-loop between type-1 T-cell polarization and DC maturation (80). The synergy between sequential ligation of IFNAR1/2 and CD40 on the surface of DC is known to boost DC maturation (81). Thus, targeting myeloid cells through CD40 in an IFN-rich microenvironment, as the one of HPV⁺ TC, may represent a viable therapeutic alternative to immune-checkpoint blockade strategies.

We also detected a rare DC cluster expressing ILR3A, SIGLEC6, and AXL in close transcriptomic proximity to DC2-like cells, which has been recently termed DC5 (30), AS DC (33), and transitional DC (82). As reported previously, DC5 identity was linked to TCF4 and RUNX2 regulon activity (33, 82). As for cDC1 and DC2-like cells, DC5s expressed high levels of the transcription factor ZNF366, which is a known repressor of the pDC developmental program, thus suggesting a myeloid origin of this population, congruent with our myeloid focused sorting strategy (83). DC5s have been reported to differentiate into cDC2-like cells in vitro (30): however, despite detecting CD11c expression on some TC derived DC5s by FC, we did not observe directionality towards DC2-like cells in RNA-velocity analysis. Together these results suggest a limited polarization of DC5s into cDC2 in the TME. Nonetheless, DC5s featured higher HLA-DR, CD40, and CCR7 protein expression compared to plasmacytoid DCs. Hence, specific functions of this rare subpopulation remain to be determined. Additionally, we detected a granulocyte lineage within the CD13⁺HLA-DR⁺ population, corresponding to mast cells with moderate expression of MHC-II related transcripts. Some studies have reported that granulocytes gain MHC-II expression during inflammation, highlighting their plasticity within an appropriate niche (84). Specifically, mast cells gain MHC-II expression upon exposure to IFN-γ, IL-4, and GM-CSF in vitro, although their specific ability to present antigens and activate T-cells in vivo is under debate (85).

We have previously demonstrated the expansion of the Mono-Mac lineage in HPV⁺ TC compared to paired HT (86), which indicates an active recruitment of monocytes during carcinogenesis. In this study, we characterized a parallel differentiation process that stemmed from different monocyte polarization events, using RNA velocity analysis. Our results challenge the in vitro M1/M2 paradigm, similarly to other scRNA-seq studies (11, 35). Monocytes expressing high levels of FOS and FOSB differentiated sequentially into NLRP3 and CXCL Macros. Of note, monocytes and NLRP3 Macros displayed the highest level of MDSC gene signature, which have been reported to differentiate into tumor-associated macrophages in tumor mice models (27). However, we did not observe key immunosuppressive features of MDSCs such as expression of IL10 and downregulation of TNF in our analysis (87). Together these observations suggest that monocytes and NLRP3 macrophages have limited immunosuppressive capacity in HPV⁺ TC. Nonetheless, NLRP3 macrophages differentiated into CXCL Macros, which were only found in high tumor stages and featured the highest activity of the hypoxia-related regulons HIF1A and CREB and their target genes downstream including VEGFA and CXCL1/3/5/8 (88, 89). Our data suggest that in high tumor stages CXCL macrophages develop in the presence of hypoxic conditions, and trigger angiogenesis and recruitment of polymorphonuclear-MDSC possibly via CXCR2 and ACKR1 (90, 91). Indeed, several studies have reported accumulation of polymorphonuclear-MDSCs in high tumor stages (92) and an association with poor prognosis for HNC patients (93). The second branch of the Mono-Mac lineage sourced from ISG Monos, which uniquely featured activity and expression of the IFNα master regulator IRF7 (94), possibly through trimerization of STAT1/2 and IRF9 (95). ISG Monos expressed TRAIL, which may be involved in tumoricidal mechanisms, as described in vitro for IFNα-polarized monocytes (96). In addition, ISG Monos expressed the highest levels of CXCL9/10/11, to potentially recruit CD8 T-cell in the TME through CXCR3, and displayed enrichment of TAP-independent antigen presentation via the MHC-I pathway. Classically, these traits have been attributed to cDC1s (97), but they might also be governed by IRF8 in human monocytes (98, 99) as described in mice (100, 101). Interestingly, ISG Monos sequentially differentiated into C1Q and act Macros, the latter of which had a moderate positive impact on HNC patient survival. Like actDCs, act Macros TF expression profile was associated with RELB and NFKB1/2, but the expression levels and gene-network activity were not as pronounced. We further validated the presence of act Macros by FC, proposing SDC2, CD163, CCR7, PD-L1, and CD40 as phenotypic surface markers. Mulder et al. described similar macrophage populations in lung, liver, and colon cancer: however, they linked their activity to recruitment of Treg and attenuation of T-cell responses via PDL1 and PDL2, features that we did not predict in our analysis (35). The function and location of CCR7⁺ act Macros may deserve further attention to elucidate if this population can migrate to draining lymph nodes, or if it is rather involved in intra-tumoral gradients of CCL19/21 such as those derived from tertiary lymphoid structures in HPV⁺ tumors (79). act Macros show marker similarity with other reported populations such as human tonsillar macrophages and CD169⁺ subcapsular sinus macrophages (32, 102). Hence, act Macros could be involved in the induction of TFH polarization (32) and the transfer of antigens to cDC1, favoring cross-presentation to CD8⁺ T-cells as described in draining lymph nodes and spleen (102, 103).

In summary, we have provided a comprehensive characterization of myeloid cell identity and dynamics in HPV⁺ TC and HT and validated our findings using 26-plex flow cytometry. Consistent with the idea that HPV⁺ TC is enriched in type-1 immune responses, we have characterized distinct polarization states of DC and Mono-Mac lineages, highlighting the impact of a type-I/II IFN-rich TME on myeloid cells. Among all populations, cDC1s stand as an ideal therapeutic target due to the biased production of their precursors, their increased abundance and capacity to mature in the tumor lesion, and their positive impact on TC and HNC patient survival. Our data support the conceptualization of targeted-therapeutic strategies to modulate intra-tumoral cDC1 activity, and delineate myeloid cell profiles that potentially may be used for patient stratification and treatment selection.

Data availability statement

The original contributions presented in the study are publicly available. This data can be found here: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE219210, Accession number; GSE219210.

A full collection of the pipeline that replicates the data analysis outlined in this article will also be made available as an R script repository on Github upon request.

Ethics statement

The studies involving human participants were reviewed and approved by Swedish Ethical Review Authority (ref. no. 2017/580). The patients/participants provided their written informed consent to participate in this study.

Author contributions

Conceptualization, DJ, LG, and ML. Experiments, DJ, AS and CA. Data analysis, DJ and AA. Data visualization, DJ. Providing human samples, SS, DA, and LG. Supervision LG and ML. Writing—original draft, DJ, LG and ML. Writing—review and editing, all authors. All authors contributed to the article and approved the submitted version.

Funding

This work was supported by grants from EU Horizon 2020 Framework programme for research and innovation (EU-MSCACOFUND, 754299 and 847583, CanFaster), the Cancera Foundation, Laryngfonden, ACTA Oto-Laryngologica Foundation, Henning and Ida Persson’s Research Foundation, Mrs. Berta Kamprad’s Cancer Foundation (FBKS-2022-8-368), and the Faculty of Engineering, Lund University.

Acknowledgments

We thank the Clinical Genomics Lund, (SciLifeLab) and Center for Translational Genomics (Lund University) for providing expertise and service with sequencing; National Bioinformatics Infrastructure Sweden (SciLifeLab) for advice in bioinformatic analysis; Cytek’s Nordic specialist team for insights in flow cytometry panel optimization; Immunology Section, Lund University for access to the Cytek Aurora instrument; and Mrs. Charlotte Cervin Hoberg, Department of ORL, Head & Neck Surgery, Skåne University Hospital, for clinical assistance. The right to re-use and adapt figures previously published in Springer Nature was obtained from the Copyright Clearance Center.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2022.1087843/full#supplementary-material

References

1. Gillison ML, Chaturvedi AK, Anderson WF, Fakhry C. Epidemiology of human papillomavirus-positive head and neck squamous cell carcinoma. J Clin Oncol (2015) 33(29):3235–42. doi: 10.1200/JCO.2015.61.6995

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Hammarstedt L, Lindquist D, Dahlstrand H, Romanitan M, Dahlgren LO, Joneberg J, et al. Human papillomavirus as a risk factor for the increase in incidence of tonsillar cancer. Int J Cancer (2006) 119(11):2620–3. doi: 10.1002/ijc.22177

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Haeggblom L, Attoff T, Yu J, Holzhauser S, Vlastos A, Mirzae L, et al. Changes in incidence and prevalence of human papillomavirus in tonsillar and base of tongue cancer during 2000-2016 in the Stockholm region and Sweden. Head Neck (2019) 41(6):1583–90. doi: 10.1002/hed.25585

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Faust H, Eldenhed Alwan E, Roslin A, Wennerberg J, Forslund O. Prevalence of human papillomavirus types, viral load and physical status of Hpv16 in head and neck squamous cell carcinoma from the south Swedish health care region. J Gen Virol (2016) 97(11):2949–56. doi: 10.1099/jgv.0.000611

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Nasman A, Romanitan M, Nordfors C, Grun N, Johansson H, Hammarstedt L, et al. Tumor infiltrating Cd8+ and Foxp3+ lymphocytes correlate to clinical outcome and human papillomavirus (Hpv) status in tonsillar cancer. PloS One (2012) 7(6):e38711. doi: 10.1371/journal.pone.0038711

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Nordfors C, Grun N, Tertipis N, Ahrlund-Richter A, Haeggblom L, Sivars L, et al. Cd8+ and Cd4+ tumour infiltrating lymphocytes in relation to human papillomavirus status and clinical outcome in tonsillar and base of tongue squamous cell carcinoma. Eur J Cancer (2013) 49(11):2522–30. doi: 10.1016/j.ejca.2013.03.019

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Hong AM, Vilain RE, Romanes S, Yang J, Smith E, Jones D, et al. Pd-L1 expression in tonsillar cancer is associated with human papillomavirus positivity and improved survival: Implications for anti-Pd1 clinical trials. Oncotarget (2016) 7(47):77010–20. doi: 10.18632/oncotarget.12776

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Polesel J, Menegaldo A, Tirelli G, Giacomarra V, Guerrieri R, Baboci L, et al. Prognostic significance of pd-L1 expression in patients with primary oropharyngeal squamous cell carcinoma: A meta-analysis. Front Oncol (2021) 11:787864. doi: 10.3389/fonc.2021.787864

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Ferris RL, Blumenschein G Jr., Fayette J, Guigay J, Colevas AD, Licitra L, et al. Nivolumab for recurrent squamous-cell carcinoma of the head and neck. N Engl J Med (2016) 375(19):1856–67. doi: 10.1056/NEJMoa1602252

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Burtness B, Harrington KJ, Greil R, Soulieres D, Tahara M, de Castro G Jr., et al. Pembrolizumab alone or with chemotherapy versus cetuximab with chemotherapy for recurrent or metastatic squamous cell carcinoma of the head and neck (Keynote-048): A randomised, open-label, phase 3 study. Lancet (2019) 394(10212):1915–28. doi: 10.1016/S0140-6736(19)32591-7

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Zhang L, Li Z, Skrzypczynska KM, Fang Q, Zhang W, O'Brien SA, et al. Single-cell analyses inform mechanisms of myeloid-targeted therapies in colon cancer. Cell (2020) 181(2):442–59 e29. doi: 10.1016/j.cell.2020.03.048

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Cabeza-Cabrerizo M, Cardoso A, Minutti CM, Pereira da Costa M, Reis e Sousa C. Dendritic cells revisited. Annu Rev Immunol (2021) 39(1):131–66. doi: 10.1146/annurev-immunol-061020-053707

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Dutertre C-A, Becht E, Irac SE, Khalilnezhad A, Narang V, Khalilnezhad S, et al. Single-cell analysis of human mononuclear phagocytes reveals subset-defining markers and identifies circulating inflammatory dendritic cells. Immunity (2019) 51(3):573–89.e8. doi: 10.1016/j.immuni.2019.08.008

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Binnewies M, Mujal AM, Pollack JL, Combes AJ, Hardison EA, Barry KC, et al. Unleashing type-2 dendritic cells to drive protective antitumor Cd4+ T cell immunity. Cell (2019) 177(3):556–71.e16. doi: 10.1016/j.cell.2019.02.005

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Bottcher JP, Bonavita E, Chakravarty P, Blees H, Cabeza-Cabrerizo M, Sammicheli S, et al. Nk cells stimulate recruitment of Cdc1 into the tumor microenvironment promoting cancer immune control. Cell (2018) 172(5):1022–37 e14. doi: 10.1016/j.cell.2018.01.004

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Segura E, Durand M, Amigorena S. Similar antigen cross-presentation capacity and phagocytic functions in all freshly isolated human lymphoid organ-resident dendritic cells. J Exp Med (2013) 210(5):1035–47. doi: 10.1084/jem.20121103

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Nizzoli G, Krietsch J, Weick A, Steinfelder S, Facciotti F, Gruarin P, et al. Human Cd1c+ dendritic cells secrete high levels of il-12 and potently prime cytotoxic T-cell responses. Blood (2013) 122(6):932–42. doi: 10.1182/blood-2013-04-495424

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Schneider K, Marbaix E, Bouzin C, Hamoir M, Mahy P, Bol V, et al. Immune cell infiltration in head and neck squamous cell carcinoma and patient outcome: A retrospective study. Acta Oncol (2018) 57(9):1165–72. doi: 10.1080/0284186X.2018.1445287

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Seminerio I, Kindt N, Descamps G, Bellier J, Lechien JR, Mat Q, et al. High infiltration of Cd68+ macrophages is associated with poor prognoses of head and neck squamous cell carcinoma patients and is influenced by human papillomavirus. Oncotarget (2018) 9(13):11046–59. doi: 10.18632/oncotarget.24306

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Marcus B, Arenberg D, Lee J, Kleer C, Chepeha DB, Schmalbach CE, et al. Prognostic factors in oral cavity and oropharyngeal squamous cell carcinoma. Cancer (2004) 101(12):2779–87. doi: 10.1002/cncr.20701

PubMed Abstract | CrossRef Full Text | Google Scholar

21. He KF, Zhang L, Huang CF, Ma SR, Wang YF, Wang WM, et al. Cd163+ tumor-associated macrophages correlated with poor prognosis and cancer stem cells in oral squamous cell carcinoma. BioMed Res Int (2014) 2014:838632. doi: 10.1155/2014/838632

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Pettersen JS, Fuentes-Duculan J, Suarez-Farinas M, Pierson KC, Pitts-Kiefer A, Fan L, et al. Tumor-associated macrophages in the cutaneous scc microenvironment are heterogeneously activated. J Invest Dermatol (2011) 131(6):1322–30. doi: 10.103/jid.2011.9

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Lewis JS, Landers RJ, Underwood JC, Harris AL, Lewis CE. Expression of vascular endothelial growth factor by macrophages is up-regulated in poorly vascularized areas of breast carcinomas. J Pathol (2000) 192(2):150–8. doi: 10.1002/1096-9896(2000)9999:9999<::AID-PATH687>3.0.CO;2-G

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Kelly A, Gunaltay S, McEntee CP, Shuttleworth EE, Smedley C, Houston SA, et al. Human monocytes and macrophages regulate immune tolerance Via integrin Alphavbeta8-mediated tgfbeta activation. J Exp Med (2018) 215(11):2725–36. doi: 10.1084/jem.20171491

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Goudot C, Coillard A, Villani AC, Gueguen P, Cros A, Sarkizova S, et al. Aryl hydrocarbon receptor controls monocyte differentiation into dendritic cells versus macrophages. Immunity (2017) 47(3):582–96 e6. doi: 10.1016/j.immuni.2017.08.016

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Tang-Huau TL, Gueguen P, Goudot C, Durand M, Bohec M, Baulande S, et al. Human in vivo-generated monocyte-derived dendritic cells and macrophages cross-present antigens through a vacuolar pathway. Nat Commun (2018) 9(1):2570. doi: 10.1038/s41467-018-04985-0

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Corzo CA, Condamine T, Lu L, Cotter MJ, Youn JI, Cheng P, et al. Hif-1alpha regulates function and differentiation of myeloid-derived suppressor cells in the tumor microenvironment. J Exp Med (2010) 207(11):2439–53. doi: 10.1084/jem.20100587

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Zhang Q, He Y, Luo N, Patel SJ, Han Y, Gao R, et al. Landscape and dynamics of single immune cells in hepatocellular carcinoma. Cell (2019) 179(4):829–45 e20. doi: 10.1016/j.cell.2019.10.003

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Coillard A, Segura E. In vivo differentiation of human monocytes. Front Immunol (2019) 10:1907. doi: 10.3389/fimmu.2019.01907

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Villani AC, Satija R, Reynolds G, Sarkizova S, Shekhar K, Fletcher J, et al. Single-cell rna-seq reveals new types of human blood dendritic cells, monocytes, and progenitors. Science (2017) 356(6335):eaah4573. doi: 10.1126/science.aah4573

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Bourdely P, Anselmi G, Vaivode K, Ramos RN, Missolo-Koussou Y, Hidalgo S, et al. Transcriptional and functional analysis of Cd1c(+) human dendritic cells identifies a Cd163(+) subset priming Cd8(+)Cd103(+) T cells. Immunity (2020) 53(2):335–52 e8. doi: 10.1016/j.immuni.2020.06.002

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Durand M, Walter T, Pirnay T, Naessens T, Gueguen P, Goudot C, et al. Human lymphoid organ Cdc2 and macrophages play complementary roles in T follicular helper responses. J Exp Med (2019) 216(7):1561–81. doi: 10.1084/jem.20181994

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Brown CC, Gudjonson H, Pritykin Y, Deep D, Lavallée V-P, Mendoza A, et al. Transcriptional basis of mouse and human dendritic cell heterogeneity. Cell (2019) 179(4):846–63.e24. doi: 10.1016/j.cell.2019.09.035

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Ramos RN, Missolo-Koussou Y, Gerber-Ferder Y, Bromley C, Bugatti M, Núñez NG, et al. Tissue-resident Folr2+ macrophages associate with tumor-infiltrating Cd8+ T cells and with increased Survival of breast cancer patients. BioRxiv (2021). doi: 10.1101/2021.04.12.439412

CrossRef Full Text | Google Scholar

35. Mulder K, Patel AA, Kong WT, Piot C, Halitzki E, Dunsmore G, et al. Cross-tissue single-cell landscape of human monocytes and macrophages in health and disease. Immunity (2021) 54(8):1883–900.e5. doi: 10.1016/j.immuni.2021.07.007

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Nalio Ramos R, Missolo-Koussou Y, Gerber-Ferder Y, Bromley CP, Bugatti M, Nunez NG, et al. Tissue-resident Folr2(+) macrophages associate with Cd8(+) T cell infiltration in human breast cancer. Cell (2022) 185(7):1189–207 e25. doi: 10.1016/j.cell.2022.02.021

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Martinez FO, Gordon S, Locati M, Mantovani A. Transcriptional profiling of the human monocyte-to-Macrophage differentiation and polarization: New molecules and patterns of gene expression. J Immunol (2006) 177(10):7303–11. doi: 10.4049/jimmunol.177.10.7303

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Kanitakis J, Morelon E, Petruzzo P, Badet L, Dubernard JM. Self-renewal capacity of human epidermal langerhans cells: Observations made on a composite tissue allograft. Exp Dermatol (2011) 20(2):145–6. doi: 10.1111/j.1600-0625.2010.01146.x

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Reynolds G, Vegh P, Fletcher J, Poyner EFM, Stephenson E, Goh I, et al. Developmental cell programs are Co-opted in inflammatory skin disease. Science (2021) 371(6527):eaba6500. doi: 10.1126/science.aba6500

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Santegoets SJ, Duurland CL, Jordanova EJ, van Ham VJ, Ehsan I, Loof NM, et al. Cd163(+) cytokine-producing Cdc2 stimulate intratumoral type 1 T cell responses in Hpv16-induced oropharyngeal cancer. J Immunother Cancer (2020) 8(2):8:e001053. doi: 10.1136/jitc-2020-001053

CrossRef Full Text | Google Scholar

41. Cillo AR, Kurten CHL, Tabib T, Qi Z, Onkar S, Wang T, et al. Immune landscape of viral- and carcinogen-driven head and neck cancer. Immunity (2020) 52(1):183–99 e9. doi: 10.1016/j.immuni.2019.11.014

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Kurten CHL, Kulkarni A, Cillo AR, Santos PM, Roble AK, Onkar S, et al. Investigating immune and non-immune cell interactions in head and neck tumors by single-cell rna sequencing. Nat Commun (2021) 12(1):7338. doi: 10.1038/s41467-021-27619-4

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Puram SV, Tirosh I, Parikh AS, Patel AP, Yizhak K, Gillespie S, et al. Single-cell transcriptomic analysis of primary and metastatic tumor ecosystems in head and neck cancer. Cell (2017) 171(7):1611–24.e24. doi: 10.1016/j.cell.2017.10.044

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Hao Y, Hao S, Andersen-Nissen E, Mauck WM 3rd, Zheng S, Butler A, et al. Integrated analysis of multimodal single-cell data. Cell (2021) 184(13):3573–87 e29. doi: 10.1016/j.cell.2021.04.048

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Zappia L, Oshlack A. Clustering trees: A visualization for evaluating clusterings at multiple resolutions. Gigascience (2018) 7(7):giy083. doi: 10.1093/gigascience/giy083

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Alshetaiwi H, Pervolarakis N, McIntyre LL, Ma D, Nguyen Q, Rath JA, et al. Defining the emergence of myeloid-derived suppressor cells in breast cancer using single-cell transcriptomics. Sci Immunol (2020) 5(44):eaay6017. doi: 10.1126/sciimmunol.aay6017

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Mende N, Bastos HP, Santoro A, Sham K, Mahbubani KT, Curd A, et al. Quantitative and molecular differences distinguish adult human medullary and extramedullary haematopoietic stem and progenitor cell landscapes. BioRxiv (2020). doi: 10.1101/2020.01.26.919753

CrossRef Full Text | Google Scholar

48. Regev A, Teichmann SA, Lander ES, Amit I, Benoist C, Birney E, et al. The human cell atlas. Elife (2017) 6. doi: 10.7554/eLife.27041

CrossRef Full Text | Google Scholar

49. Aibar S, Gonzalez-Blas CB, Moerman T, Huynh-Thu VA, Imrichova H, Hulselmans G, et al. Scenic: Single-cell regulatory network inference and clustering. Nat Methods (2017) 14(11):1083–6. doi: 10.1038/nmeth.4463

PubMed Abstract | CrossRef Full Text | Google Scholar

50. La Manno G, Soldatov R, Zeisel A, Braun E, Hochgerner H, Petukhov V, et al. Rna velocity of single cells. Nature (2018) 560(7719):494–8. doi: 10.1038/s41586-018-0414-6

PubMed Abstract | CrossRef Full Text | Google Scholar

51. Angerer P, Haghverdi L, Buttner M, Theis FJ, Marr C, Buettner F. Destiny: Diffusion maps for Large-scale single-cell data in r. Bioinformatics (2016) 32(8):1241–3. doi: 10.1093/bioinformatics/btv715

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Cancer Genome Atlas N. Comprehensive genomic characterization of head and neck squamous cell carcinomas. Nature (2015) 517(7536):576–82. doi: 10.1038/nature14129

PubMed Abstract | CrossRef Full Text | Google Scholar

53. Hanzelmann S, Castelo R, Guinney J. Gsva: Gene set variation analysis for microarray and rna-seq data. BMC Bioinf (2013) 14:7. doi: 10.1186/1471-2105-14-7

CrossRef Full Text | Google Scholar

54. Terry M Therneau TL, Atkinson E, Crowson C. Package ‘Survival’, Vol. 3. (2022). pp. 3–1, CRAN.

Google Scholar

55. Kuleshov MV, Jones MR, Rouillard AD, Fernandez NF, Duan Q, Wang Z, et al. Enrichr: A comprehensive gene set enrichment analysis web server 2016 update. Nucleic Acids Res (2016) 44(W1):W90–7. doi: 10.1093/nar/gkw377

PubMed Abstract | CrossRef Full Text | Google Scholar

56. Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, et al. Gene ontology: Tool for the unification of biology. Gene Ontol. Consortium. Nat Genet (2000) 25(1):25–9. doi: 10.1038/75556

CrossRef Full Text | Google Scholar

57. Jin S, Guerrero-Juarez CF, Zhang L, Chang I, Ramos R, Kuan CH, et al. Inference and analysis of cell-cell communication using cellchat. Nat Commun (2021) 12(1):1088. doi: 10.1038/s41467-021-21246-9

PubMed Abstract | CrossRef Full Text | Google Scholar

58. Phillips RJ, Lutz M, Premack B. Differential signaling mechanisms regulate expression of cc chemokine receptor-2 during monocyte maturation. J Inflammation (Lond) (2005) 2:14. doi: 10.1186/1476-9255-2-14

CrossRef Full Text | Google Scholar

59. Hambleton S, Salem S, Bustamante J, Bigley V, Boisson-Dupuis S, Azevedo J, et al. Irf8 mutations and human dendritic-cell immunodeficiency. N Engl J Med (2011) 365(2):127–38. doi: 10.1056/NEJMoa1100066

PubMed Abstract | CrossRef Full Text | Google Scholar

60. Rosa FF, Pires CF, Kurochkin I, Halitzki E, Zahan T, Arh N, et al. Single-cell transcriptional profiling informs efficient reprogramming of human somatic cells to cross-presenting dendritic cells. Sci Immunol (2022) 7(69):eabg5539. doi: 10.1126/sciimmunol.abg5539

PubMed Abstract | CrossRef Full Text | Google Scholar

61. Cao L, Liu S, Li Y, Yang G, Luo Y, Li S, et al. The nuclear matrix protein safa surveils viral rna and facilitates immunity by activating antiviral enhancers and super-enhancers. Cell Host Microbe (2019) 26(3):369–84 e8. doi: 10.1016/j.chom.2019.08.010

PubMed Abstract | CrossRef Full Text | Google Scholar

62. Sun L, Kong B, Sheng X, Sheu JJ, Shih Ie M. Dendritic cells transduced with rsf-1/Hbxap gene generate specific cytotoxic T lymphocytes against ovarian cancer in vitro. Biochem Biophys Res Commun (2010) 394(3):633–8. doi: 10.1016/j.bbrc.2010.03.038

PubMed Abstract | CrossRef Full Text | Google Scholar

63. Duhen T, Duhen R, Montler R, Moses J, Moudgil T, de Miranda NF, et al. Co-Expression of Cd39 and Cd103 identifies tumor-reactive Cd8 T cells in human solid tumors. Nat Commun (2018) 9(1):2724. doi: 10.1038/s41467-018-05072-0

PubMed Abstract | CrossRef Full Text | Google Scholar

64. Breton G, Zheng S, Valieris R, Tojal da Silva I, Satija R, Nussenzweig MC. Human dendritic cells (Dcs) are derived from distinct circulating precursors that are precommitted to become Cd1c+ or Cd141+ dcs. J Exp Med (2016) 213(13):2861–70. doi: 10.1084/jem.20161135

PubMed Abstract | CrossRef Full Text | Google Scholar

65. Heidkamp GF, Sander J, Lehmann CHK, Heger L, Eissing N, Baranska A, et al. Human lymphoid organ dendritic cell identity is predominantly dictated by ontogeny, not tissue microenvironment. Sci Immunol (2016) 1(6):eaai7677. doi: 10.1126/sciimmunol.aai7677

PubMed Abstract | CrossRef Full Text | Google Scholar

66. Villar J, Segura E. Decoding the heterogeneity of human dendritic cell subsets. Trends Immunol (2020) 41(12):1062–71. doi: 10.1016/j.it.2020.10.002

PubMed Abstract | CrossRef Full Text | Google Scholar

67. Pellicioli ACA, Bingle L, Farthing P, Lopes MA, Martins MD, Vargas PA. Immunosurveillance profile of oral squamous cell carcinoma and oral epithelial dysplasia through dendritic and T-cell analysis. J Oral Pathol Med (2017) 46(10):928–33. doi: 10.1111/jop.12597

PubMed Abstract | CrossRef Full Text | Google Scholar

68. Wu X, Briseno CG, Durai V, Albring JC, Haldar M, Bagadia P, et al. Mafb lineage tracing to distinguish macrophages from other immune lineages reveals dual identity of langerhans cells. J Exp Med (2016) 213(12):2553–65. doi: 10.1084/jem.20160600

PubMed Abstract | CrossRef Full Text | Google Scholar

69. Fainaru O, Woolf E, Lotem J, Yarmus M, Brenner O, Goldenberg D, et al. Runx3 regulates mouse tgf-Beta-Mediated dendritic cell function and its absence results in airway inflammation. EMBO J (2004) 23(4):969–79. doi: 10.1038/sj.emboj.7600085

PubMed Abstract | CrossRef Full Text | Google Scholar

70. Saba Y, Aizenbud I, Matanes D, Koren N, Barel O, Zubeidat K, et al. Early antitumor activity of oral langerhans cells is compromised by a carcinogen. Proc Natl Acad Sci USA (2022) 119(3):e2118424119. doi: 10.1073/pnas.2118424119

PubMed Abstract | CrossRef Full Text | Google Scholar

71. Verdijk RM, Mutis T, Esendam B, Kamp J, Melief CJ, Brand A, et al. Polyriboinosinic polyribocytidylic acid (Poly(I:C)) induces stable maturation of functionally active human dendritic cells. J Immunol (1999) 163(1):57–61. doi: 10.4049/jimmunol.163.1.57

PubMed Abstract | CrossRef Full Text | Google Scholar

72. Honda K, Sakaguchi S, Nakajima C, Watanabe A, Yanai H, Matsumoto M, et al. Selective contribution of ifn-Alpha/Beta signaling to the maturation of dendritic cells induced by double-stranded rna or viral infection. Proc Natl Acad Sci USA (2003) 100(19):10872–7. doi: 10.1073/pnas.1934678100

PubMed Abstract | CrossRef Full Text | Google Scholar

73. Fainaru O, Shseyov D, Hantisteanu S, Groner Y. Accelerated chemokine receptor 7-mediated dendritic cell migration in Runx3 knockout mice and the spontaneous development of asthma-like disease. Proc Natl Acad Sci USA (2005) 102(30):10598–603. doi: 10.1073/pnas.0504787102

PubMed Abstract | CrossRef Full Text | Google Scholar

74. Sato K, Kawasaki H, Nagayama H, Enomoto M, Morimoto C, Tadokoro K, et al. Tgf-beta 1 reciprocally controls chemotaxis of human peripheral blood monocyte-derived dendritic cells Via chemokine receptors. J Immunol (2000) 164(5):2285–95. doi: 10.4049/jimmunol.164.5.2285

PubMed Abstract | CrossRef Full Text | Google Scholar

75. Imai K, Minamiya Y, Koyota S, Ito M, Saito H, Sato Y, et al. Inhibition of dendritic cell migration by transforming growth factor-Beta1 increases tumor-draining lymph node metastasis. J Exp Clin Cancer Res (2012) 31(1):3. doi: 10.1186/1756-9966-31-3

PubMed Abstract | CrossRef Full Text | Google Scholar

76. Hernandez A, Burger M, Blomberg BB, Ross WA, Gaynor JJ, Lindner I, et al. Inhibition of nf-kappa b during human dendritic cell differentiation generates anergy and regulatory T-cell activity for one but not two human leukocyte antigen Dr mismatches. Hum Immunol (2007) 68(9):715–29. doi: 10.1016/j.humimm.2007.05.010

PubMed Abstract | CrossRef Full Text | Google Scholar

77. Hu Y, Park-Min KH, Yarilina A, Ivashkiv LB. Regulation of stat pathways and Irf1 during human dendritic cell maturation by tnf-alpha and Pge2. J Leukoc Biol (2008) 84(5):1353–60. doi: 10.1189/jlb.0107040

PubMed Abstract | CrossRef Full Text | Google Scholar

78. Johansson-Percival A, Ganss R. Therapeutic induction of tertiary lymphoid structures in cancer through stromal remodeling. Front Immunol (2021) 12:674375. doi: 10.3389/fimmu.2021.674375

PubMed Abstract | CrossRef Full Text | Google Scholar

79. Ruffin AT, Cillo AR, Tabib T, Liu A, Onkar S, Kunning SR, et al. B cell signatures and tertiary lymphoid structures contribute to outcome in head and neck squamous cell carcinoma. Nat Commun (2021) 12(1):3349. doi: 10.1038/s41467-021-23355-x

PubMed Abstract | CrossRef Full Text | Google Scholar

80. Cella M, Scheidegger D, Palmer-Lehmann K, Lane P, Lanzavecchia A, Alber G. Ligation of Cd40 on dendritic cells triggers production of high levels of interleukin-12 and enhances T cell stimulatory capacity: T-T help Via apc activation. J Exp Med (1996) 184(2):747–52. doi: 10.1084/jem.184.2.747

PubMed Abstract | CrossRef Full Text | Google Scholar

81. Greyer M, Whitney PG, Stock AT, Davey GM, Tebartz C, Bachem A, et al. T Cell help amplifies innate signals in Cd8(+) dcs for optimal Cd8(+) T cell priming. Cell Rep (2016) 14(3):586–97. doi: 10.1016/j.celrep.2015.12.058

PubMed Abstract | CrossRef Full Text | Google Scholar

82. Leylek R, Alcantara-Hernandez M, Granja JM, Chavez M, Perez K, Diaz OR, et al. Chromatin landscape underpinning human dendritic cell heterogeneity. Cell Rep (2020) 32(12):108180. doi: 10.1016/j.celrep.2020.108180

PubMed Abstract | CrossRef Full Text | Google Scholar

83. Nutt SL, Chopin M. Transcriptional networks driving dendritic cell differentiation and function. Immunity (2020) 52(6):942–56. doi: 10.1016/j.immuni.2020.05.005

PubMed Abstract | CrossRef Full Text | Google Scholar

84. Singhal S, Bhojnagarwala PS, O'Brien S, Moon EK, Garfall AL, Rao AS, et al. Origin and role of a subset of tumor-associated neutrophils with antigen-presenting cell features in early-stage human lung cancer. Cancer Cell (2016) 30(1):120–35. doi: 10.1016/j.ccell.2016.06.001

PubMed Abstract | CrossRef Full Text | Google Scholar

85. Poncet P, Arock M, David B. Mhc class ii-dependent activation of Cd4+ T cell hybridomas by human mast cells through superantigen presentation. J Leukoc Biol (1999) 66(1):105–12. doi: 10.1002/jlb.66.1.105

PubMed Abstract | CrossRef Full Text | Google Scholar

86. Jimenez DG, Sobti A, Askmyr D, Sakellariou C, Santos SC, Swoboda S, et al. Tonsillar cancer with high Cd8(+) T-cell infiltration features increased levels of dendritic cells and transcriptional regulation associated with an inflamed tumor microenvironment. Cancers (Basel) (2021) 13(21):5341. doi: 10.3390/cancers13215341

PubMed Abstract | CrossRef Full Text | Google Scholar

87. Bergenfelz C, Leandersson K. The generation and identity of human myeloid-derived suppressor cells. Front Oncol (2020) 10:109. doi: 10.3389/fonc.2020.00109

PubMed Abstract | CrossRef Full Text | Google Scholar

88. Lucia K, Wu Y, Garcia JM, Barlier A, Buchfelder M, Saeger W, et al. Hypoxia and the hypoxia inducible factor 1alpha activate protein kinase a by repressing rii beta subunit transcription. Oncogene (2020) 39(16):3367–80. doi: 10.1038/s41388-020-1223-6

PubMed Abstract | CrossRef Full Text | Google Scholar

89. Fang HY, Hughes R, Murdoch C, Coffelt SB, Biswas SK, Harris AL, et al. Hypoxia-inducible factors 1 and 2 are important transcriptional effectors in primary macrophages experiencing hypoxia. Blood (2009) 114(4):844–59. doi: 10.1182/blood-2008-12-195941

PubMed Abstract | CrossRef Full Text | Google Scholar

90. Salazar N, Zabel BA. Support of tumor endothelial cells by chemokine receptors. Front Immunol (2019) 10:147. doi: 10.3389/fimmu.2019.00147

PubMed Abstract | CrossRef Full Text | Google Scholar

91. Greene S, Robbins Y, Mydlarz WK, Huynh AP, Schmitt NC, Friedman J, et al. Inhibition of mdsc trafficking with sx-682, a Cxcr1/2 inhibitor, enhances nk-cell immunotherapy in head and neck cancer models. Clin Cancer Res (2020) 26(6):1420–31. doi: 10.1158/1078-0432.CCR-19-2625

PubMed Abstract | CrossRef Full Text | Google Scholar

92. Trellakis S, Bruderek K, Hutte J, Elian M, Hoffmann TK, Lang S, et al. Granulocytic myeloid-derived suppressor cells are cryosensitive and their frequency does not correlate with serum concentrations of colony-stimulating factors in head and neck cancer. Innate Immun (2013) 19(3):328–36. doi: 10.1177/1753425912463618

PubMed Abstract | CrossRef Full Text | Google Scholar

93. Lang S, Bruderek K, Kaspar C, Hoing B, Kanaan O, Dominas N, et al. Clinical relevance and suppressive capacity of human myeloid-derived suppressor cell subsets. Clin Cancer Res (2018) 24(19):4834–44. doi: 10.1158/1078-0432.CCR-17-3726

PubMed Abstract | CrossRef Full Text | Google Scholar

94. Marie I, Durbin JE, Levy DE. Differential viral induction of distinct interferon-alpha genes by positive feedback through interferon regulatory factor-7. EMBO J (1998) 17(22):6660–9. doi: 10.1093/emboj/17.22.6660

PubMed Abstract | CrossRef Full Text | Google Scholar

95. Mogensen TH. Irf and stat transcription factors - from basic biology to roles in infection, protective immunity, and primary immunodeficiencies. Front Immunol (2018) 9:3047. doi: 10.3389/fimmu.2018.03047

PubMed Abstract | CrossRef Full Text | Google Scholar

96. Tecchio C, Huber V, Scapini P, Calzetti F, Margotto D, Todeschini G, et al. Ifnalpha-stimulated neutrophils and monocytes release a soluble form of tnf-related apoptosis-inducing ligand (Trail/Apo-2 ligand) displaying apoptotic activity on leukemic cells. Blood (2004) 103(10):3837–44. doi: 10.1182/blood-2003-08-2806

PubMed Abstract | CrossRef Full Text | Google Scholar

97. Spranger S, Dai D, Horton B, Gajewski TF. Tumor-residing Batf3 dendritic cells are required for effector T cell trafficking and adoptive T cell therapy. Cancer Cell (2017) 31(5):711–23 e4. doi: 10.1016/j.ccell.2017.04.003

PubMed Abstract | CrossRef Full Text | Google Scholar

98. Langlais D, Barreiro LB, Gros P. The macrophage Irf8/Irf1 regulome is required for protection against infections and is associated with chronic inflammation. J Exp Med (2016) 213(4):585–603. doi: 10.1084/jem.20151764

PubMed Abstract | CrossRef Full Text | Google Scholar

99. Oosenbrug T, van de Graaff MJ, Haks MC, van Kasteren S, Ressing ME. An alternative model for type I interferon induction downstream of human Tlr2. J Biol Chem (2020) 295(42):14325–42. doi: 10.1074/jbc.RA120.015283

PubMed Abstract | CrossRef Full Text | Google Scholar

100. Kim S, Bagadia P, Anderson DA 3rd, Liu TT, Huang X, Theisen DJ, et al. High amount of transcription factor Irf8 engages Ap1-irf composite elements in enhancers to direct type 1 conventional dendritic cell identity. Immunity (2020) 53(4):759–74 e9. doi: 10.1016/j.immuni.2020.07.018

PubMed Abstract | CrossRef Full Text | Google Scholar

101. Nixon BG, Kuo F, Liu M, Capistrano K, Do M, Franklin RA, et al. Irf8 governs tumor-associated macrophage control of T cell exhaustion. bioRxiv (2020). doi: 10.2139/ssrn.3554068

CrossRef Full Text | Google Scholar

102. van Dinther D, Veninga H, Iborra S, Borg EGF, Hoogterp L, Olesek K, et al. Functional Cd169 on macrophages mediates interaction with dendritic cells for Cd8(+) T cell cross-priming. Cell Rep (2018) 22(6):1484–95. doi: 10.1016/j.celrep.2018.01.021

PubMed Abstract | CrossRef Full Text | Google Scholar

103. Grabowska J, Lopez-Venegas MA, Affandi AJ, den Haan JMM. Cd169(+) macrophages capture and dendritic cells instruct: The interplay of the gatekeeper and the general of the immune system. Front Immunol (2018) 9:2472. doi: 10.3389/fimmu.2018.02472

PubMed Abstract | CrossRef Full Text | Google Scholar