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

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

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

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

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.

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.

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.

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.

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

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.

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

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

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.

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.

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

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.

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

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.

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

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.