The Tumor Immune Microenvironment Architecture Correlates with Risk of Recurrence in Head and Neck Squamous Cell Carcinoma

The tumor immune microenvironment plays an important role in disease progression and response to therapy (1). For example, a higher concentration of regulatory T cells (Treg; ref. 2) and CD33⁺IL4α⁺ myeloid-derived suppressor cells (MDSC) correlate with the risk of recurrence. In contrast, a higher number of CD8⁺ tumor-infiltrating lymphocytes (TIL) increases recurrence-free survival in head and neck squamous cell carcinoma (HNSCC), breast cancer, lung cancer, and other human malignancies (3–7). Clinical staining of leukocytes is routinely performed with conventional one- or two-color IHC. However, studying more than two markers requires a careful selection of primary antibodies or the use of consecutive tissue sections. This is problematic for studying samples with low tissue availability and makes it extremely difficult to colocalize markers at the single-cell level. Alternatively, the tumor microenvironment is assessed by flow cytometry or single-cell (sc) RNA sequencing on single-cell suspensions from the tumors. These techniques allow for the simultaneous analysis of multiple markers and genes but require tissue dissociation. Thus, they do not reveal the topographic location of the leukocyte subsets and raise concerns about the possible underrepresentation of some leukocyte subsets that may be lost during tissue processing (8, 9). The emergence of multiplexed tissue imaging overcame these challenges. In particular, co-detection by indexing (CODEX) generates detailed information on the distribution of different cellular phenotypes while maintaining the morphologic context of healthy and diseased tissues (10–12). Although this technology has been used primarily on tumor tissues microarrays (11–15) that may exclude functional structures distant from the neoplastic edge, it revealed a complex organization of the tumor microenvironment into spatially defined cellular structures whose composition and frequency could be used to discriminate patients with HNSCC with lymph nodal metastases from those with N0 tumors (10–15).

Here we performed a prospective clinical study of patients with stage II to IV surgically resectable primary or recurrent HNSCC to (i) characterize the tumor microenvironment by multicolor flow cytometry and CODEX on large HNSCC sections; (ii) identify structures distant from the tumor edge; and (iii) evaluate whether the tumor architecture can be associated with prognosis. Compared with CODEX, we found that flow cytometry underestimated the number of neutrophils and polymorphonuclear (PMN) MDSCs and overestimated TILs frequency. Using CODEX data, we identified 20 cell types (CT) and 11 cellular neighborhoods (CN) that portray the HNSCC microenvironment. Moreover, we provided evidence that a “type 1” tertiary lymphoid structure (TLS1) and CD38⁺CD31⁺plasma cells are associated with fewer recurrences after salvage surgery in HNSCC.

Nine patients (median age 71, range 61–78) with stage II to IV HNSCC undergoing surgical resection were enrolled in the study under the University of Miami Institutional Review Board (IRB)–approved protocol no. 20200268 after signing an IRB-approved written informed consent. Seven patients had HNSCC in the oral cavity, one in the hypopharynx, and one in the oropharynx. Two patients were female, 5 had a smoking history, 6 had a recurrent tumor, and 5 had prior radiation to the field. No tumor was HPV related. Complete clinical and pathologic parameters are reported in Table 1.

Eight weeks old Balb/c female mice (Taconic) were injected orthotopically with 10⁵ 4T1 (ATCC) cells in the fat pad of the right abdominal mammary gland. Mice were euthanized when tumors reached 0.7 cm in diameter. All experiments with mice were performed under the Institutional Animal Care and Use Committee (IACUC) protocol no. 22051, approved by the ethical committee of the University of Miami. Part of each tumor and spleen was processed for flow cytometry, and part was snap frozen in OCT for subsequent CODEX analysis.

Human tumors were harvested at surgery and processed within 1 hour. Mouse tissues were processed immediately after euthanasia. Part of the tumor was flash-frozen in OCT compound for subsequent CODEX staining, and part was mechanically dissociated into single-cell suspension for multicolor flow cytometry analysis. Briefly, tumor specimens were washed twice with PBS and incubated for 20 minutes at 37°C with five volumes of PBS containing Clostridium histolyticum collagenase type IV (10 mg/mL; Sigma), MgCl₂ (100 μmol/L), and CaCl₂ (100 μmol/L). The solution was passed every 10′ in a needleless 5 mL syringe. Cells were then passed through a 70 μm cell strainer, washed with PBS, and resuspended for staining.

Single-cell suspensions (10⁶ cells) from HNSCC tumors were stained for 20′ at 4°C with the Zombie Violet Fixable Viability Dye (BioLegend), washed with PBS, and stained with the following antibody panels for 15′ at 4°C. (i) Myeloid cell panel: Anti-human CD33-FITC (clone HIM3–4; BD), anti-human Lox1-APC (clone 15C4; BioLegend), anti-human CD124-PE (clone 25463; R&D Systems), anti-human CD14-APC-H7 (clone MϕP9; BD Biosciences), anti-human CD15-BV711 (clone W6D3; BD Biosciences), anti-human HLA-DR V500 (clone G46–6; BD Biosciences), and anti-human CD11b-BV605 (clone ICRF44; BD Biosciences). (ii) Lymphoid cell panel: Anti-human CD3-Alexa Fluor 700 (clone OKT3; eBioscience), anti-human CD247-PE (clone 6B10.2; eBioscience), anti-human CD4-BV711 (clone SK3; BD Biosciences), anti-human CD8-BV605 (clone SK1; BD Biosciences), anti-human CD69-APC-Cy7 (clone FN50; BD Biosciences), and anti-human CD154-PE/Dazzle 594 (clone24–31; BioLegend). Cells in the lymphoid panel were then fixed and permeabilized using the Foxp3/Transcription Factor Staining Buffer Set (eBioscience) as per the manufacturer's instructions and with the anti-human Foxp3-APC (clone 236A/E7; eBioscience). Mouse tumors and spleens were stained with NIR live/dead dye (Thermo Fisher Scientific) and with rat anti-mouse antibodies specific for CD45 (AF700, clone AB_2572116; BioLegend), CD11b (bv711, clone AB_394002; BD Biosciences), Ly6 g (Pac.Blue, clone AB_1727562; BD Biosciences), Ly6c (FITC, clone AB_394628; BD Biosciences), CD3 (PE clone AB_1834427; Thermo Fisher Scientific), CD4 (FITC, clone AB_394585; BD Biosciences), and CD8 (Pac.Blue, clone AB_394571; BD biosciences). Data were acquired within 2 hours on a BD LSRII flow cytometer analyzer equipped with the following wavelengths lasers: 405 nm (50 mW), 488 nm (50 mW), 532 nm (150 mW), and 640 nm (40 mW) or with the 5 lasers Cytoflex LX-19C (Beckman Coulter). Analyses were performed in FCS Express 7 (De Novo software).

The following oligonucleotides conjugated antibodies against human epitopes were acquired from Akoya: CD45-BX001 (HI30), CD19-BX003 (HIB19), CD8-BX004 (SK1), CD38-BX007 (HB-7), CD279-BX014, (EH12.2H7), CD3-BX015 (UCHT1), pan-cytokeratin-BX019 (AE-1/AE-3), CD4-BX021 (SK3), CD90-BX022 (5E10), HLA-DR-BX026 (L243), CD31-BX032 (WM59), CD34-BX035 (561), and CD69-BX041 (FN50). In addition, we generated the following oligonucleotides-tagged antibodies against human epitopes using the CODEX Conjugation Kit (Akoya) as per manufacturer instructions: FOXP3-BX002 (AB_467556; Thermo Fisher Scientific), PDL1-BX005 (AB_2687808; BioXcell), CD14-BX006 (AB_830675; BioLegend), CD163-BX010 (AB_1088991; BioLegend), CD11b-BX013 (AB_314154; BioLegend), CD33-BX016 (AB_2562818; BioLegend), CD56-BX017 (AB_314444; BioLegend), IL4Ra-BX020 (AB_2126871; R&D Systems), CD45RA-BX024 (AB_314406; BioLegend), CD123-BX027 (RRID: AB_467453; Thermo Fisher Scientific), CCR7-BX028 (AB_2563726; BioLegend), LOX1-BX029 (AB_2562272; BioLegend), CD25-BX030 (AB_2561752; BioLegend), CD15-BX033 (AB_314194; BioLegend), CD66b-BX036 (AB_2728422; BioLegend), CD68-BX037 (AB_1089058; BioLegend), and CD16-BX042 (AB_314202; BioLegend). The following antibodies against mouse epitopes were purchased from Akoya: CD45-BX007 (clone AKYP0005), CD3-BX021 (clone AKYP0035), CD4-BX026 (clone AKYP0041), CD8a-BX029 (clone 53–6.7), CD11b-BX025 (cloneM1/17), and Ly6g-BX024 (clone AKYP0039). The anti-mouse Ly6c antibody (clone HK1.4; BioLegend) was conjugated with the BX027 barcode. Conjugation of anti-mouse and anti-human antibodies was verified by PAGE, and optimal titration was performed by single antibody staining on HNSCC tumor specimens of mouse spleen.

OCT embedded flash-frozen specimens were cut with a cryostat microtome in 5 μm sections and placed on a poly-l-lysine-coated coverslip. Sections were dried with drierite beads for 2′, fixed in acetone for 10′ at RT, rehydrated in hydration buffer (3 × 2′; Akoya), and fixed in prestain fix buffer (1.6% PFA in hydration buffer) for 10′. Sections were washed in hydration buffer and incubated for 20′ in staining buffer (Akoya). Tissues were then stained in a humidity chamber with the cocktail of oligonucleotides tagged antibodies diluted in blocking buffer (N, G, J, S, Akoya) for 3 hours at 4°C. Sections were then washed three times for 2′ in staining buffer and fixed by one incubation poststain fixative buffer (1.6% PFA in staining buffer) for 10′, three washes with PBS, one incubation in ice-cold methanol for 5′, three washed in PBS, and a 20′ incubation with the final fixative solution (Akoya) at RT. Sections were washed nine times in PBS and kept in storage buffer until image acquisition.

Fluorochrome conjugated reported oligonucleotides were hybridized on 12-cycle experiments with the Akoya PhenoCycler connected to a BZ-X700 microscope (Keyence), following manufacturer instructions. Images were automatically acquired on 121 adjacent 20× fields (5.6 mm × 4.2 mm, resolution 377.47 nm/pixel) with 7 z-slices (1.5 μm/z pitch) and 30% tile overlap, with 10, 250, 350, and 350 milliseconds acquisition time for DAPI, AF488, Cy3, and cy5 filters, respectively. The resulting 40,656 images per specimen were processed with Akoya's PhenoCycler Instrument Manager (version 1.3) and CODEX Processor (version 1.7.2) for segmentation using the default parameters.

Composite images were loaded in ImageJ (version 1.53) using the CODEX MAV extension (version 1.5.0.8) to check image quality and export single-cell fluorescence signals as CSV files.

Data have been preprocessed using Seurat 4.2.1 functions to analyze image-based spatial data (https://satijalab.org/seurat/articles/spatial_vignette_2.html) in R 4.2.2. Briefly, CSV data files from CODEX MAV have been loaded in R using Seurat LoadAkoya function, and protein signals have been normalized with the centered log-ratio based normalization. To detect the cell phenotypes in the HNSCC tissue (i.e., malignant, stromal, immune, or vascular cells), we first performed principal component analysis (PCA) for dimensionality reduction and then cluster analysis in the low-dimensional space. Before applying PCA, we scaled the data and determined the number of components for downstream analysis. Then, we ran dimensional reduction for each sample using the top 11 principal components and graph-based clustering at different resolutions (ranging from 0.2 to 0.8; Supplementary Table S1). Cell clusters have been visualized on a protein intensity-based uniform manifold approximation and projection (UMAP) and on their spatial location. To associate cell phenotypes to clusters, levels of protein markers have been displayed as dot plots for each cluster at any explored clustering resolution. For each sample, the optimal number of clusters (i.e., the clustering resolution) was determined based on visual inspection of cluster location and marker expression (Supplementary Table S1). Clusters with a similar morphological appearance in the tissue and similar marker expression profiles were merged, and artifacts were removed. Next, we manually assigned CT annotations to clusters at the selected resolution based on the average expression of protein markers in each cluster. After removing cells marked as artifacts, raw count and normalized data matrices of all samples have been merged using the merge function of the Seurat package. Supervised annotation and merging resulted in a final list of 20 CTs used to annotate all cells in every sample (Supplementary Table S1).

To identify CNs that are regions with a characteristic local composition of cell phenotypes, we implemented in R the neighborhood analysis presented in ref. 11. For each of the 792,349 cells from the nine samples, we identified a window consisting of W nearest neighboring cells (including the center cell) using the nn2 function of the RANN R package (version 2.6.1). The nn2 function uses a k-dimensional tree to find a given number of near neighbors (here, W) for each point identified by the X and Y coordinates in the input dataset. These windows were then clustered by their composition with respect to the 20 CTs previously identified by graph-based clustering and supervised annotation. Specifically, each window was converted to a vector of length 20 containing the frequency of each of the 20 CTs among the W neighbors. Subsequently, windows have been clustered using the MiniBatchKmeans function of the ClusterR package (version 1.2.9) implementing the Mini-batch K-means clustering algorithm with a given value of K. Each cell was then allocated to the CN of its surrounding window using the predict_MBatchKMeans function of the ClusterR package. We applied the entire procedure and identified CNs for different combinations of W (ranging from 7 to 200) and K (ranging from 7 to 14).

To quantify pairwise cell–cell contacts within CNs, we first determined the direct neighbors of each cell within each CN using the Delaunay triangulation implemented in the triangulate function of the RTriangle R package (version 1.6–0.11). First, we calculated the number of contacts between cells of types i and j as the set of edges N_ijk between cells of type i and j in CN k returned by the triangulation. Then, displayed pairwise contacts among different phenotype cells as circle plots using the chordDiagram function of the circlize R package (version 0.4.15).

To evaluate if cell–cell contacts of cells assigned to a given phenotype were influenced by the CN, we compared the relative cell–cell contact frequencies (RF) for all phenotypes in each pair of CNs (e.g., TLS1 and TLS2) across all patients containing both CNs. We calculated the relative cell–cell contact frequency between cells of types i and j in each CN k as N_ijk/N_ik where N_ijk is the number of edges of the triangulation between cells of type I and j and N_ik is the number of edges of type i in CN k. We compared the relative cell–cell contact frequency in pairs of CNs using the ttest function of the R stats package. We discarded cell–cell interactions between pairs of phenotypes present in less than three samples. P values have been adjusted for multiple comparisons with the p.adjust function of the R stats package using the method of Benjamini and Hochberg (16).

We used the procedure described in ref. 14 for tensor decomposition analysis, adapting the Python code available at https://github.com/nolanlab/NeighborhoodCoordination. For each patient, we constructed a joint composition matrix quantifying the frequency of each CT in each CN and collected these matrices in tridimensional tensors with dimensions given by patients by CTs by CNs. We generated distinct tensors for groups of patients with different clinical characteristics. Non-negative Tucker decomposition, as implemented in the Tensorly Python package, was applied to each tensor (17). We selected the ranks in each dimension equal to 3, 8, and 8 by visually inspecting the decomposition error for different combinations of the ranks in multiple decomposition runs. Modules for CT and CNs have been generated from the factor matrices of the CT and CN spaces, respectively. The interactions composing a tissue module were obtained from the 8 × 8 slices of the 3 × 8 × 8 core tensor linking the factor matrices. In the graphical representation of tensor decomposition analysis, the contribution of CTs and CNs to tissue modules was quantified by their values in the corresponding slice of the core tensor, and these values were used for color shading. Elements with a value lower than 0.1 in the corresponding slice of the core tensor have been greyed out.

The analysis of differential marker expression has been performed using the FindMarkers function of Seurat R package (version 4.2.1). To identify differentially expressed markers in the comparison of a given phenotype (e.g., plasma cells) in a CN (e.g., TLS1) versus the same phenotype in all other CNs, we set the fold change and FDR thresholds at 1.2 and 0.05, respectively.

Raw data for The Cancer Genome Atlas (TCGA) HNSC collection were downloaded as raw counts from the TCGA repository using the TCGAbiolinks R package (v. 2.29; ref. 18). Data normalization was performed using the edgeR R package (v. 3.30.2; ref. 18). Specifically, raw counts were normalized to counts per million mapped reads (CPM) and fragments per kilobase million (FPKM), and only genes with a CPM greater than 1 in at least three samples were retained. To identify two groups of tumors with either high or low TLS1 and CD31⁺plasma cells, we used the classifier described in ref. 19, which is a classification rule based on gene expression signature scores. Briefly, the signature scores have been obtained summarizing the standardized expression levels of CD38 and PECAM1 (CD31) genes into a combined score with zero mean. Tumors were classified as CD38/CD31 signature “Low” if the combined score was smaller than the median signature score and as CD38/CD31 signature “High” vice versa. This classification was applied to expression values of the TCGA HNSC samples with disease-specific survival information (n = 493).

Plasma cell fractions have been quantified using CIBERSORT (20) on the TCGA HNSC collection. Briefly, the FPMK expression matrix was uploaded to the CIBERSORT R script (version 1.04) as a mixture file, and CIBERSORT was run in absolute mode with the LM22 signature gene file, 100 permutations, and quantile normalization. In absolute mode, CIBERSORT scales relative cellular fractions into a score that reflects the absolute proportion of each CT in a mixture. Although not expressed as a fraction, the absolute score can be directly compared both between- and within-samples (21). Samples have been divided into two groups based on the median of the absolute scores for the plasma CT.

To evaluate the prognostic value of the CD38/CD31 signature and the plasma cell fraction, we estimated the probabilities of disease-specific survival using the Kaplan–Meier method. The Kaplan–Meier curves were compared using the log-rank (Mantel–Cox) test to confirm these findings. The P value was calculated according to the standard normal asymptotic distribution. Survival analysis was performed using functions of survival (version 3.4–0) and survminer (version 0.4.9) packages. Multivariate analysis of the association of age, sex pathological stage, histologic grade, plasma cell fraction, and CD38/CD31 signature classification factors with disease-free survival has been performed using the analyse_multivariate function of the survivalAnalysis R package (version 0.3.0). Forest plot representation of the multivariate analysis results has been obtained using the forest function of the forestplotter R package (version 0.2.3). To compare the multivariate models, we subtracted the residual deviance of the multivariate model with age, sex pathological stage, histologic grade, and plasma cell fraction (Model 1) from the one with age, sex pathologic stage, histological grade, plasma cell fraction, and CD38/CD31 signature classification (Model 2) and tested this difference against a chi-square distribution with one degree of freedom. All analyses have been done in R 4.2.2.

The single-cell data table of clustered, annotated CTs with metadata are publicly available in Mendeley (https://data.mendeley.com/datasets/t2yvtwnjx7/1). R scripts created for this study are available at https://github.com/bicciatolab/CODEX_HNSCC.

Additional data analyzed in this study were obtained from the TCGA repository at https://portal.gdc.cancer.gov/. All other raw data generated in this study are available upon request from the corresponding author.

To analyze the immunologic landscape of the HNSCC microenvironment, we performed highly multiplexed microscopy using CODEX on nine surgical specimens of patients with stage II to IV HPV-negative squamous cell carcinoma explants of the oral cavity and oropharynx. We selected and validated a panel of 30 antibodies that allows the discrimination of the main leukocyte subsets infiltrating the tumor and auxiliary markers for tumor, vascular, and stromal cells on OCT embedded flash-frozen sections. To minimize sampling errors, we processed an area of 42.6 mm² derived from 121 best-focused stitched 20× fields using the Akoya software and default parameters for quality control, image processing, and cell segmentation. We evaluated the image quality, segmentation, and fluorescence signals on CODEX MAV; we excluded staining artifacts and exported the fluorescence signals, metadata, and spatial coordinates of all cells for processing in R using an ad hoc computational framework. For each specimen, we clustered cells at different resolutions (ranging from 0.2. to 0.8) using the centered log-ratio normalized signals and the graph-based clustering algorithm of the Seurat package (22), and we visualized cell clusters on low-dimensional embeddings (t-SNE and UMAP) and on the spatial locations. We associated clusters to CTs at any explored resolution and determined the optimal cluster resolution by visual inspection of marker expression and spatial localization (Supplementary Table S1). We removed staining artifacts (n = 945 cells) and merged clusters with similar marker expression profiles and morphological appearance in the tissue. Finally, we generated a merged dataset accounting for 792,349 cells from all nine original specimens. This procedure allowed the definition of 20 phenotypes that, in a UMAP embedding, segregate cells according to the type and not on the specimen of origin (Fig. 1A; Supplementary Table S1). As expected, we identified cytokeratin-positive tumor cells, CD31⁺epithelial (blood vessel) cells, the main lymphocyte subsets (NK and NKT, CD8⁺cytotoxic and CD4⁺ T helper lymphocytes, Treg, B cells, and plasma cells), various populations of myeloid cells (monocytes, M1 and M2 macrophages, neutrophils, Lox⁺PMN-MDSCs; ref. 23), and hematopoietic stem and progenitor cells (HSPC; ref. 24). In addition, we observed the presence of CD11b⁺cells expressing the transcription factor FOXP3, IL4Rα, and PDL1 (Fig. 1B; Supplementary Fig. S1A), of a heterogeneous group of myeloid cells expressing PDL1 (PDL1⁺myeloid), and of CD90⁺ mesenchymal stem-like stromal cells (25) that surround the neoplastic islands (Fig. 1C; Supplementary Figs. S1B–S1K). Interestingly, the unsupervised clustering identified a population of cells that we named capsule, characterized by low or negative expression of most markers and, when present, contoring the neoplastic nests (Fig. 1C; Supplementary Figs. S1B–S1K). Finally, we compared the frequencies of each cell phenotype in all tumor specimens. Although some CTs were present in all or most patients (e.g., tumor, CD90⁺ cells, vessel, PMN-MDSC), others (as capsule, NKT, B cells, PDL1⁺myeloid) were present only in a limited number of samples (Fig. 1D; Supplementary Table S1).

In summary, the analysis of CODEX data revealed a complex tumor microenvironment characterized by multiple cell subsets whose composition differs across tumor specimens.

Flow cytometry is considered the gold standard for enumerating tumor-infiltrating leukocyte subsets. However, this technique requires extensive tissue processing to obtain a single-cell suspension, which may impact the proper representation of different cell subsets. In contrast, CODEX allows the evaluation of the tumor microenvironment using a limited amount of tissue and minimal specimen processing. Despite the limited sample processing, sampling errors may also affect CODEX analyses. We divided each surgical specimen into two parts to compare the performances of flow cytometry and CODEX in enumerating cell subsets. One part was processed and analyzed by multicolor flow cytometry; the other was frozen and evaluated by multiparametric CODEX immunofluorescence microscopy. We identified the major lymphoid and myeloid subsets in FCS Express using similar manual gating strategies on the data originated by the two techniques (Fig. 2A and B). Compared with CODEX, we observed that flow cytometry significantly overestimated the percentage of CD3⁺ lymphocyte subsets and dramatically underestimated the PMN-MDSC fraction (Fig. 2C and D). We reasoned that these differences might derive from tissue processing that, in the case of tumors preparation for flow cytometry, requires collagenase treatment, filtration, and multiple centrifugations, and we hypothesized that a tissue requiring a simpler processing (i.e., a mouse spleen that requires only mechanical dissociation, filtration, and RBC lysis) could minimize the differences between the two technologies. To test this hypothesis, we enumerated the leukocytes from the spleen and the tumor of mice bearing the 4T1 mammary carcinoma using both flow cytometry and CODEX. Briefly, spleens and tumors were harvested; half of each specimen was snap-frozen and analyzed by CODEX, whereas the other half was processed and analyzed by flow cytometry (Fig. 2E). As observed in human HNSCC, compared with CODEX, flow cytometry overestimated the number of T cells and mMDSC and underestimated the number of PMN-MDSC in the tumor. In contrast, no significant differences were observed when splenic leukocytes were enumerated by the two techniques (Fig. 2E). These results indicate that the low frequency of PMN-MDSC usually found by flow cytometry and scRNA sequencing in mouse and human tumors (26) and the lower PMN-MDSC/m-MDSC ratio observed in tumors (compared with the one observed in circulation or in the spleen) can be attributed to technical artifacts of the tissue preparation.

Recent data from Nolan's laboratory (10, 11, 14, 15, 27) show that host and neoplastic cells are not randomly distributed in the tumor microenvironment, but they are instead spatially organized in communities called CNs, which may provide the topological and phenotypic architecture for optimal cell–cell interaction and function. An example of a CN is the TLS, a spatial organization of infiltrating immune cells that can positively or negatively modulate antitumor immunity (28). To identify CNs in our samples, we quantified, for each cell of the merged dataset, the phenotype frequencies of its W nearest neighbors and clustered these topological and compositional vectors into K groups with distinctive enrichments of the 20 original CTs (11). To optimize the width of the capturing window W and the number K of CNs, we tested different combinations of neighboring cells (W range 7–200) and CNs (K range 7–14). Finally, we considered a window of W = 10 cells and opted for K = 14 CNs as this combination (i) maximized the number of CNs with a unique pattern of phenotype frequencies (enrichment score), (ii) minimized the number of CNs enriched in the most abundant CTs (i.e., tumor cells), and (iii) better recapitulated the cellular composition of the extensively characterized TLS neighborhood modeled around plasma cells, CD4⁺T cells, and stromal cells (29).

After evaluating the CT composition of each CN, we noticed that some clusters were very similar in leukocyte composition (Supplementary Fig. S2; Supplementary Table S1). We then merged the original 14 CNs into 11 CNs with a distinctive phenotype stoichiometry (Fig. 3A). We identified a “cold tumor” CN, mainly composed of neoplastic cells; a “hot tumor” CN, characterized by the presence of malignant cells and leukocytes; a “peritumoral” CN made of neoplastic and capsule cells; a “vascular” CN rich in vessel cells; and a “stromal CN” composed of CD90⁺MSC-like cells (Fig. 3A). In addition, the analysis revealed the presence of NK-rich, M1-rich, neutrophils-rich, and MDSC-rich CNs characterized by the high frequency of NK, TAM-M1, LOX1⁻neutrophil, and LOX1⁺PMN-PMN-MDSCs, respectively. Finally, different leukocyte subsets characterized two distinct CNs, that is, TLS1, rich in plasma cells and CD4⁺Foxp3⁻ lymphocytes, and TLS2, characterized by the presence of M2-like macrophages and Tregs. We assessed the CN frequencies in the nine tissue samples and noticed that cold tumor and TLS2 are the most abundant CNs, followed by peritumoral (when present), hot tumor, MDSC-rich, and stroma CNs (Fig. 3B; Supplementary Table S1). We observed substantial variation in the frequency of the other CNs: five samples lacked the peritumoral CN, and four were deficient in the neutrophil-rich and the M1-rich CNs (Fig. 3B; Supplementary Table S1). We then evaluated whether some correlation existed among the CN frequencies. We found that the stroma CNs negatively and positively correlated with the TLS1 (r = −0.80; P = 0.006) and the peritumoral CNs (r = 0.73; P = 0.020), respectively. Similarly, we observed a negative correlation (r = −0.67; P = 0.043) between the frequency of NK and cold tumor CNs (Supplementary Fig. S2B). As expected, the hot tumor and the peritumoral CNs separated the cold tumor CN from all other CNs (Fig. 3C; Supplementary Figs. S2C–S2T). Generally, the MDSC-rich CN and the stroma occupied the space between the neoplastic nests and contained the other CNs. In most samples, the NK-rich CN was located distantly from the tumoral CNs, whereas neutrophil-rich CN, when present, was in proximity and interacted with the neoplastic cells (Fig. 3C; Supplementary Figs. S2C–S2T).

Finally, we reasoned that the CNs could differ not only in the CT composition but also in the interaction between CTs. Thus, we quantified the pairwise cell–cell contacts for CTs within CNs (Fig. 4; Supplementary Fig. S3). We observed a higher degree of homotypic cell–cell contacts in some CNs (e.g., neoplastic cells with neoplastic cells, CD90⁺cells with CD90⁺cells, etc.) whereas, in other neighborhoods, cell–cell contacts were more heterotypic. For instance, in the peritumoral CN, capsule cells interacted with other capsule cells but also with tumor and leukocyte cells, forming a barrier between these latter CTs (Supplementary Fig. S3). We found the largest number of heterotypic interactions in the TLS1, TLS2, and hot tumor CNs (Fig. 4; Supplementary Fig. S3). Plasma cells dominated the cellular interactions in the TLS1 neighborhood. These cells interacted with other plasma cells but also with CD90⁺MSC-like cells, vessel cells, and both CD4⁺ and CD8⁺ T cells (Fig. 4). This finding is consistent with the type of TLS described in many cancer types and is generally associated with a good prognosis (30). M2 macrophages dominated the interactions within the TLS2. These macrophages interacted with Treg, PMN-MDSC, CD4, and CD8 T cells, suggesting a phenotype that may inhibit antitumor immunity (Fig. 4).

Because the same CTs can be found (although at different frequencies) in both TLS1 and TLS2, we evaluated whether cell–cell interactions for a given CT were influenced by the cell neighborhood it belonged to (i.e., did a plasma cell interact with different CTs if it belongs to TLS1 or TLS2?). Thus, we compared the relative cell–cell contact frequencies (11) for all phenotypes in different pairs of CNs (e.g., TLS1 and TLS2; hot and cold tumor) across all patients.

We found that, in TLS1, plasma cells interacted mostly with other plasma cells and were contacted by NK, CD90⁺stromal, and PNM-MDSC cells (Fig. 5A; Supplementary Table S1). Instead, in TLS2, M2 macrophages, PMN-MDSC, CD8 T lymphocytes, Treg, vessel cells, and CD90 cells interacted significantly more with M2 macrophages. In TLS2, additional significant contacts were those between Tregs with other Tregs and capsule cells with Tregs (Fig. 5A; Supplementary Table S1).

In the hot tumor CN, approximately one-third of the neoplastic cell contacts were with leukocyte subsets such as neutrophils, CD4 and CD8 T cells, but also with M2 macrophages and PMN-MDSC (Fig. 4; Supplementary Fig. S3). In contrast, in the cold tumor CN, most neoplastic cell interactions were homotypic (Fig. 4; Supplementary Fig. S3). The differential analysis of relative contact frequencies confirmed these observations and indicated that, in the cold tumor CN, neoplastic cells form a significantly higher number of homotypic interactions and that the few leukocytes interact primarily with tumor cells (Fig. 5B; Supplementary Table S1). Conversely, in the hot tumor neighborhood, CD8⁺ Tcells, vessels, Treg, M2 macrophages, and stromal cells interact significantly more with M2 macrophages (Fig. 5B; Supplementary Table S1).

These analyses revealed a complex tumor microenvironment in which CTs are organized in spatially distinct units that are different in cell composition, cell–cell interactions, and colocalization. These CNs may be functionally important for tumor growth and immune escape or the generation of effective antitumor immunity.

The composition of the tumor microenvironment is a major determinant of patients’ prognosis and tumor clinical and pathological features. Hence, we would expect that patterns of CTs and CNs in spatial regions are differently coupled and correlated, in patients, with different clinical phenotypes (e.g., recurring or not after surgery, presenting or not perineuronal invasion, or with and without primary tumor at surgery). As demonstrated by Schurch and colleagues, tensor decomposition can reveal how cell phenotypes and CNs are combined into distinct tissue modules in groups of patients with different clinical characteristics (11). We, therefore, constructed a joint composition matrix quantifying the frequency of each CT in each CN and collected these matrices in tridimensional tensors. Of the 9 analyzed patients, 4 recurred within a year from surgery (median recurrence time 0.65 years, range 0.35–0.92 years), 4 did not recur during the observation time (median observation time 2 years, range 1.2–2.1 years), and 1 patient was uncontactable after surgery. To determine whether any characteristic of the tumor microenvironment was associated with early (<1 year) tumor recurrence after salvage surgery, we constructed a 3D tensor for those patients that recurred after surgery and one for those that did not recur (Fig. 6A). Then, we applied non-negative Tucker tensor decomposition to the tensor of patients’ joint CT–CN matrices in recurrent and non-recurrent patient groups separately (11, 31). Tensor decomposition analysis returned a set of eight CN and eight CT modules (dashed inner rectangles in Fig. 6B) that shaped three different tissue modules in patients recurring and not recurring after surgery (solid outer rectangles in Fig. 6B). Graphically, we represented the contribution of each CT and CN to each tissue module based on its value in the corresponding slice of the core tensor and used this value for color shading. Elements with a value lower than 0.1 in the corresponding slice of the core tensor were not deemed to contribute to a given tissue module and were greyed out (Fig. 6B). In both recurrent and non-recurrent patients, the first tissue module was characterized by neoplastic, protumoral leukocytes, and stroma CN and CT components. Similarly, factors belonging to the stroma and stroma–tumor interface contributed to defining the CN and CT components in the second modules for both groups of patients, whereas the vasculature CN weighed in the definition of the second tissue module in those patients that recurred but did not contribute to the same module in non-recurrent patients. Strikingly, the third tissue module was defined by different CN and CT components in recurrent and non-recurrent patients: while in recurrent patients, NK CT weighed the most in the CT and CN components, in non-recurrent patients, the NK, M1 macrophages, and TLS1 weighed in the definition of the CN components whereas M1 macrophages, NK, CD4 T cells and plasma cells defined the CT component (Fig. 6B). In accordance with this finding, the TLS1 CN was more frequent in patients whose tumor did not recur (Fig. 7A) and its presence significantly correlated with the absence of recurrence (Supplementary Fig. S4A) and was predictive of recurrence-free survival (Fig. 7B). Similarly, plasma cells weighed in the tissue module of no recurrent tumors (Fig. 6B), were associated with the absence of recurrence (Supplementary Fig. S4B) and predicted recurrence-free survival (Fig. 7C). Of note, besides perineuronal invasion (P = 0.0042), no other clinical or pathological parameters (e.g., smoking history, age, sex, history of previous HNSCC) predicted recurrence after surgery in univariate analysis. Because of the limited amounts of patients, multivariate analysis could not be performed.

Interestingly, we also found that the frequency of cells from the cold tumor CN positively correlated with the perineural invasion (PNI) and the depth of invasion (DOI). In contrast, the peritumoral and MDSCs-rich CNs correlated with the tumor size and previously irradiated fields, respectively. We also observed a significant negative correlation between patient smoking history and the frequency of cells from the NK-rich CN (Supplementary Fig. S4A).

Given the small number of patients evaluated in this study, we investigated whether we could identify distinctive markers of plasma cells in TLS1 and used them to interrogate larger collections of patients with HNSCC. Hence, we performed a differential protein expression analysis between plasma cells in the TLS1 and those belonging to other CNs. We found that plasma cells in TLS1 expressed significantly higher levels of the CD38 and CD31 proteins (Fig. 7D and E). These markers define a subset of long-lived IgG-secreting plasma cells (32) usually found in primary and secondary lymphoid organs but not in circulation (33–35).

To evaluate whether the expressions of CD38 and CD31 in the tumor microenvironment had a prognostic value, we classified HNSCC tumors from the TCGA dataset based on the single and combined expression of these two genes. We found that higher transcriptional levels of CD38 and CD31 significantly correlated with better disease-specific survival (Fig. 7F), whereas the expression of either CD38 or CD31 alone was not or only slightly predictive (P > 0.01). Of note, the expression level of the CD38/CD31 signature was the only feature significantly (P = 0.0006) associated with a better prognosis in multivariate analysis with age, sex pathologic stage, and histological grade (Supplementary Fig. S4C). A similar significant (P = 0.0055) association was found when using as covariate the plasma cell fraction inferred using CIBERSORT (20) deconvolution algorithm (Supplementary Figs. S4D and S4E). Adding CD38/CD31 signature to this multivariate analysis conferred additional predictive value (Supplementary Fig. S4E). Similar correlations between the CD38/CD31 signature and a better prognosis were found using the Kaplan–Meier Plotter resource (www.kmplot.com; ref. 36) in renal clear cell carcinoma (HR, 0.42; P = 5×10⁻⁹), lung adenocarcinoma (HR, 0.65; P = 0.0061), and thymoma (HR, 0; P = 0.0001).

Taken together, these data support the notion that the spatial organization of the tumor microenvironment is essential for cell function and prognosis and points to TLS1 and CD31^highCD38^highplasma cells as important predictors of tumor recurrence.

Previous studies analyzed the interaction between the tumor and the immune system well and showed the pivotal role of different leukocyte subsets in neoplastic cell growth, recurrence, and metastases (37). However, only a few retrospective studies evaluated the tumor microenvironment globally, considering cell-to-cell interactions, spatial context, and relationships between CT phenotype and surrounding cells (11, 12). These studies revealed that the tumor is a highly organized tissue in which leukocytes interact within structures or CNs (11, 12). The presence, frequency, and composition (in terms of cell phenotypes) of these micro-niches might be relevant in shaping the optimal environment for leukocyte interaction and function and, ultimately, in determining tumor characteristics.

We performed the first prospective study in patients with HNSCCs to evaluate whether features of the fine architecture of the tumor microenvironment are associated with the risk of recurrence after salvage surgery. We analyzed the HNSCC samples using the CODEX-based multiparametric immunofluorescence method with an optimized 30 markers antibody panel and an ad hoc computational framework. Using a manual gating strategy, we compared the performances of CODEX in quantifying the frequency of the main leukocyte subsets with those of multicolor flow cytometry. We found that, compared with CODEX, flow cytometry overestimated the T-cell frequency and dramatically underestimated the polymorphonucleate populations and, in particular, the PMN-MDSCs. MDSCs emerged as a crucial population that inhibits tumor immunity and promotes tumor growth and metastases (38). Interestingly, different studies showed that the PMN-MDSC/M-MDSC ratio differs in tumors and in the periphery (26). These observations could indicate a real biological property (i.e., a preferential migration of M-MDSC to the tumor or their enhanced survival) or could derive from a selective loss of PMN-MDSC during isolation. Because PMN-MDSC and M-MDSC are functionally distinct (38), this conundrum is clinically essential for developing effective immunotherapies. Our data in mice and humans support the notion that the conventional techniques significantly underestimate PMN-MDSC, suggesting caution in determining the frequency of leukocyte subsets in the tumor microenvironment following enzymatic and mechanical processing of the tissue (Fig. 2).

Because cell phenotyping by hand-gating is subjective, sensitive to segmentation noise, and cumbersome with high-dimensional data (39, 40), we analyzed the normalized CODEX data with a clustering-based approach followed by manual annotation and identified 20 distinct CTs (Fig. 1). Besides the expected lymphoid and myeloid subsets present at various ratios across the analyzed tumors, this analysis revealed the presence of CD90⁺stromal cells, a population of cells lining around the neoplastic nests that we called capsule cells, and a rare population of CD11b⁺Foxp3⁺ cells, positive for PDL1 and IL4Rα, whose existence and biological relevance is however still controversial (41–43). The abundance of CD90⁺cells between the neoplastic island confirmed previous studies indicating CD90⁺mesenchymal stem cells as a major contributor to HNSCC stroma (25). CD90⁺mesenchymal stem cells drove glioma, colon carcinoma, breast cancer, and HNSCC progression by increasing neoplastic cell proliferation, migration, and adhesion (44).

The identification of the capsule cells was an unexpected finding of this study. These cells were negative for all tested lineage markers but delineated the neoplastic islets in 4 of the 9 patients and were positive for PDL1, suggesting a possible protective role from the immune response. Similar structures composed of PDL1⁺cancer associated fibroblasts (CAF) have been reported in breast cancer and shown to correlate with a better prognosis (45). Other studies, however, indicate CAFs as promoters of tumor-induced tolerance, neo-angiogenesis, and cancer stemness (46). Further studies are needed for a finer characterization of the capsule cells, their attribution to the CAF family, and the elucidation of their role in HNSCC progression.

As mentioned above, one of the most important findings of this paper is the identification of CNs as spatial entities accounting for the complexity of the HNSCC tumor microenvironment. Although some CNs were characterized by the homotypic interaction of a particular CT (e.g., MDSC-rich, M1-rich, NK-rich CNs), others were defined by the presence and interaction of different CTs. We focused on two CNs that we named TLS1 and TLS2. The TLS2 CN was characterized by M2 macrophages interacting with Treg, PMN-MDSC, stromal cells, vasculature, and different effector leukocytes such as CD4 and CD8 T cells, NK, and neutrophils. To our knowledge, a similar structure has been described only by a few reports and was shown to correlate with a worse prognosis in patients with breast and pancreatic cancers (47, 48). In contrast, the TLS1 neighborhood was dominated by plasma cells that interact with CD4⁺ and CD8⁺ effector T cells, monocytes, neutrophils, and M1 macrophages. This structure is consistent with a tertiary lymphoid organ mimicking the lymph nodal germinal centers and correlating with a good prognosis in many human malignancies (49). In agreement with previous studies in melanoma, colon carcinoma, and breast cancer (49), our data indicate that frequency TLS1 and plasma cells are good prognostic factors in HNSCC (Figs. 6 and 7). We hypothesize that phenotype and possibly the function of a particular CT depend on the spatial cellular context in which it resides. Indeed, differential protein expression analysis on plasma cells from TLS1 compared with plasma cells from other neighborhoods revealed the high expression of the CD38 and CD31 markers (Fig. 7). Interestingly, the co-expression of both CD38 and CD31 discriminates patients with a better prognosis in the TCGA HNSC dataset and adds prognostic value to deconvolution models using the plasma cells signature to classify the patients. This finding supports the notion that cell phenotype and function depend on the spatial context and highlights the power of spatial proteomic analysis. Of note, the expression of CD31 and CD38 is associated with a population of long-lived plasma cells usually found in the bone marrow and the lymph nodes (50). Long-lived plasma cells provide a lifelong production of antibodies in mice and humans (51) and are suspected of driving persistent inflammation in autoimmune diseases (51). Although these cells have not yet been described in the tumor microenvironment, their presence might indicate the secretion of high-affinity antibodies that facilitate antitumor immunity by antibody-mediated cell cytotoxicity or by favoring antigen uptake and cross-presentation by antigen-presenting cells (51).

In summary, our study highlights the power of CODEX immune fluorescence in immuno-monitoring, reveals the extraordinary complexity of the tumor microenvironment of HNSCC patients, and supports the notion that TLS1 and long-lived plasma cells are essential determinants of antitumor immunity and tumor progression. However, it is important to note that despite the large number of images (∼100,000) we evaluated, this explorative prospective study is limited by the small number of patients. Thus, the key finding needs to be further assessed in a larger group of patients, possibly using a smaller number of markers or technologies that reduce the data size to increase the speed of computational analysis.

No disclosures were reported.

D.T. Weed: Conceptualization, resources, writing–review and editing. S. Zilio: Investigation. C. McGee: Investigation. B. Marnissi: Resources, investigation. Z. Sargi: Resources. E. Franzmann: Resources. G. Thomas: Resources. J. Leibowitz: Resources. E. Nicolli: Resources. D. Arnold: Resources. S. Bicciato: Conceptualization, resources, data curation, software, formal analysis, supervision, funding acquisition, methodology, writing–original draft, project administration, writing–review and editing. P. Serafini: Conceptualization, resources, data curation, formal analysis, supervision, funding acquisition, investigation, methodology, writing–original draft, project administration, writing–review and editing.

Research reported in this publication was supported by: (i) Fondazione AIRC under 5 per mille 2019 - ID. 22759 program to S. Bicciato; (ii) SCCC-Translational Science award to P. Serafini and D.T. Weed; and (iii) the NCI of the NIH under award number P30CA240139. The content is solely the responsibility of the authors and does not necessarily represent the official views of the NIH or other funding agencies.

The publication costs of this article were defrayed in part by the payment of publication fees. Therefore, and solely to indicate this fact, this article is hereby marked “advertisement” in accordance with 18 USC section 1734.