Abstract
Introduction:
Head and neck cancers (HNCs) are aggressive malignancies, with recurrence contributing substantially to poor patient outcomes. We used spatial transcriptomics to identify tumour and immune features associated with recurrence.
Methods:
Pathologist-annotated tumour and stromal regions from recurrent and non-recurrent HNC samples were analysed using spatially resolved transcriptomics. Exploratory differential expression analysis was used to nominate tumour-associated candidate markers for recurrence. Selected candidates were assessed by immunohistochemistry and further evaluated using independent spatial transcriptomic and bulk RNA-sequencing cohorts. Immune-cell abundance and spatial co-localisation were also examined.
Results:
PAEP, MARCO, PNPLA3, and FMO2 emerged as leading candidate markers associated with recurrence. Increased expression of PAEP, MARCO, and PNPLA3 was supported by immunohistochemistry. Integrative analysis across independent spatial and bulk RNA cohorts demonstrated consistent PAEP overexpression in recurrent HNC, particularly within tumour regions, with this pattern reproduced across datasets and profiling platforms. Non-recurrent tumours showed a trend towards greater T- and B-cell infiltration, with these populations frequently co-localising within the tumour microenvironment. In recurrent tumour regions, PAEP expression showed an inverse trend with B-cell and plasma-cell abundance, which was higher in non-recurrent samples.
Discussion:
Spatial transcriptomics identified tumour-intrinsic and immune microenvironmental features associated with HNC recurrence. PAEP emerged as the most reproducible candidate biomarker across internal and external cohorts and may be associated with reduced B-cell and plasma-cell infiltration. These findings support further validation of PAEP and related immune signatures for recurrence-risk stratification and personalised treatment.
1 Introduction
Head and neck cancers (HNC) are responsible for over 900, 000 new cases and 450, 000 deaths annually, making them one of the deadliest and most debilitating cancers globally (). Despite advancements in diagnostic and therapeutic strategies, HNC patients’ overall 5-year survival remains low, at 50-60% (–), with recurrence being a key driver of these outcomes.
Recurrence is known to often arise from minimal residual disease that persists after curative surgery (–). Additionally, tumors consist of dynamically evolving subpopulations of tumor cells that frequently acquire genetic and epigenetic changes (, ), underpinning plasticity that may lead to therapy resistance and tumor recurrence. Similarly, the ability of tumor cells to disseminate into secondary sites contributes to cancer relapse (, ). Moreover, certain metastasized cells, particularly cancer stem cells, can remain dormant for extended periods before germinating at the new site, escaping treatment that targets active tumor cells ().
Tumor microenvironment (TME) plays a pivotal role in recurrence. Immune infiltration, such as by CD8+ T lymphocytes, is associated with improved outcomes (–), while cancer-associated fibroblasts have been shown to promote the invasion and migration of cancer cells and hinder immune cell activity (–). Recent single-cell RNA sequencing (scRNA-seq) studies have illuminated the dynamic evolution of tumor and stromal compartments during HNC progression. One investigation delineated malignant epithelial trajectories spanning precancerous to recurrent disease, uncovering a tumorigenic subcluster regulated by TFDP1, alongside POSTN+ fibroblasts and SPP1+ macrophages that remodel the TME to drive progression (). In addition, exhausted CXCL13+ CD8+ T cells were identified as facilitators of extranodal expansion during recurrence. Collectively, these findings underscore how intratumoral heterogeneity and immune dysfunction promote relapse. Technologies that can reveal interactions between cancer, immune, and stromal cells within TME, such as spatial transcriptomics, may lead to the discovery of predictive biomarkers of recurrence.
Current prognosis parameters include tumor stages, anatomical sites, HPV/p16 status, comorbidities, and smoking or alcohol history. For example, the presence of p16 as a marker for HPV infection indicates a better prognosis, longer survival, and response to HPV-driven carcinogenesis (). While these factors contribute to stratifying risk, they do not capture biological complexity underlying recurrence. Early or regional recurrence increases the risk of death, yet molecular markers reflecting the progression at early stages remain lacking, leading to underestimation of survival variability among patients.
In this study, we aimed to discover biomarkers and characterize the TMEs of recurrent HNC by applying Visium Spatial Transcriptomics (ST) to compare recurrent and non-recurrent HNC, unbiasedly exploring spatially resolved expression of 18, 000 genes within the preserved tissue context. By integrating spatial data with histopathological assessment of cancer cores, tumor edges, and stromal regions, we identified potential biomarkers of recurrence. We assessed these candidates at the protein level and in external bulk and spatial transcriptomic datasets to identify potential mechanisms associated with recurrence in HNC. This work offers a valuable resource reporting spatially resolved expressions for advancing our understanding of tumor biology and improving prognostic modelling.
2 Materials and methods
Figure 1 displays the overall workflow of the study.
Figure 1
2.1 Patient recruitment
HNC patients were recruited from Royal Brisbane and Women’s Hospital (RBWH), Brisbane, following the approved ethics by the Metro South Human Research Ethics Committee, Brisbane, QLD, Australia (HREC/12/QPAH/381). Written informed consent was obtained from all HNC patients before enrolment in the study.
2.2 Samples preparation and sequencing
Formalin-fixed paraffin-embedded (FFPE) tumor tissues from 10 HNC patients, comprising recurrent (N = 4) and non-recurrent (N = 6), were used as described in Supplementary Table 1. One patient (R3) had two distinct tumor biopsies embedded into separate FFPE blocks. Both blocks were sectioned individually onto the Visium array slides but treated as a single sample during subsequent computational analyses. Clinical details are outlined in Supplementary Table 1.
To ensure tissue suitability for spatial transcriptomic analysis, RNA quality was assessed using DV200 score, reflecting the percentage of total profiled RNA fragments longer than 200 bp. About 5–10 sections (5µm thick) were used for RNA quality assessment using the RNeasy FFPE Kit (#73504, Qiagen). The DV200 score was determined via the BioAnalyzer electrophoresis using the RNA 6000 Pico Kit (#5067-1513, Agilent), applying a threshold of DV200>30% for study inclusion ().
Library preparation followed the Visium Spatial Gene Expression for FFPE User Guide (CG0000407, CG000408, CG000409 – 10x Genomics), with modifications described by Vo et al. (). FFPE blocks were sectioned at 5µm by a rotary microtome (Leica RM2245 Rotary Microtome). Sections were transferred into a preheated (41 °C) water bath for 3 minutes before being placed onto the Visium Spatial Gene Expression slide. Multiplexing was employed to optimize capture area usage. The slide was then dried at 42 °C for 3 hours and stored overnight at room temperature in a parafilm-seal 50mL-Falcon tube containing silica bead desiccants.
The following day, the slide was incubated at 60 °C for 2 hours on a thermocycler (Bio-Rad C1000 Thermal Cycler) to melt paraffin wax, immediately deparaffinized in xylene (twice, for 10 minutes), and rehydrated through an ethanol gradient (100% thrice for 3 minutes; 96% twice for 3 minutes; 70% for 3 minutes; followed by a brief water immersion). Haematoxylin and eosin (H&E) staining and imaging was performed according to the Visium Spatial for FFPE, including Deparaffinization, H&E Staining, Imaging & Decrosslinking, following the Demonstrated Protocol (CG000409 - 10x Genomics). Slide images were acquired at 20× magnification using a Zeiss Axioscan Z1 Fluorescent Imager. Following imaging, tissue sections underwent decrosslinking in HCl-preconditioned TE buffer (pH 8.0) at 70 °C for 1 hour.
For expression profiling, 25 bp probe pairs targeting the whole human transcriptome (#2000449, #2000450, 10x Genomics) were hybridized to their complementary transcript targets in situ overnight. Following hybridization, adjacent probe pairs were ligated to generate 50 bp ligation products. Endogenous tissue mRNA was subsequently digested using an RNase mix, releasing the stable ligated probe products from the tissue. The synthetic poly(A) tails of these released probes were captured by the slide-bound, spatially barcoded oligo(dT) primers.
Slide-bound cDNA was synthesized from the captured, ligated probes and eluted into 0.2 mL PCR tubes. Each generated cDNA molecule contained a unique spatial barcode, a Unique Molecular Identifier (UMI), and the specific probe sequence. Short-read Illumina-compatible libraries were constructed using the Visium Spatial Gene Expression Reagent Kits for FFPE (CG000407, 10x Genomics). Quantitative PCR was performed to determine the optimal Sample Index PCR cycle number, and the libraries were amplified for 15 cycles. Final library molecules were purified using SPRIselect beads, and cDNA quality and concentration were verified on a Bioanalyzer. Sequencing was conducted by the IMB Sequencing Facility on an Illumina NovaSeq platform (100-cycle kit), targeting a sequencing depth of approximately 100 million reads per library.
The corresponding H&E images were subsequently annotated by an HNC-specialized anatomical pathologist.
2.3 Data analyses
For full reproducibility, the code detailing the following analyses is publicly available at https://github.com/BiomedicalMachineLearning/RecurrenceHN.
2.3.1 Pre-processing and quality control
Raw sequencing data and high-resolution H&E slide images were processed using the 10x Genomics SpaceRanger pipeline. Transcripts were mapped to the GRCh38 human reference transcriptome to generate spatially resolved gene expression matrices aligned with their corresponding tissue histology. The 10x Genomics Loupe Browser was utilized to manually overlay metadata, including patient recurrence status and anatomical region. To establish morphological context, discrete tissue spots were assigned a dominant cell type or histological domain based on manual annotations provided by an HNC-specialized anatomical pathologist on the matching H&E images.
Downstream quality control (QC), filtering, and integration were conducted using the Seurat (v4.1.0) R package (). To eliminate low-quality data and technical artifacts, spots capturing fewer than 100 unique genes were excluded, and genes detected in fewer than three spots across the entire dataset were removed. Additionally, any spots located outside the boundaries of the primary tissue section or lacking pathologist-guided histological annotations were discarded.
To account for differences in sequencing depth across capture areas, the filtered raw count matrices were normalized and variance-stabilized via the SCTransform framework (). This method applies a regularized negative binomial regression model to effectively remove technical variation while preserving biological heterogeneity. To mitigate technical batch effects and ensure cross-sample comparability, a unified integrated assay was generated for downstream clustering. This was achieved using Seurat’s anchor-based integration workflow, employing Canonical Correlation Analysis (CCA) to identify shared biological states (“anchors”) across individual sample matrices. The resulting integration vectors were applied to harmonize the datasets, ensuring that downstream clustering was driven by genuine biological features rather than inter-sample technical variation.
2.3.2 Unsupervised clustering
Unsupervised clustering was performed on the SCTransform normalized and integrated dataset using the Seurat FindClusters function. The top 50 principal components (PC) were used to build the shared nearest neighors (SNN). Sub-clustering was performed on cluster 1 using the FindSubCluster function at 0.1 resolution in Seurat to gain better annotation of that cluster. The FindMarker function with Wilcoxon rank sum test () was used to identify upregulated marker genes for each cluster. Genes with an adjusted P<0.05 were selected to perform GO term analysis to identify gene groups and/or functions associated with each cluster group were analyzed by clusterProfiler R package, with a significance threshold of Benjamini-Hochberg- (BH-) adjusted P<0.05. Results of cluster gene markers and enriched pathway were visualised using the dotplot, cnetplot and emapplot functions (clusterProfiler), producing dot plots, Gene-Concept Network and enrichment maps respectively (). To find enriched clusters in recurrent or non-recurrent samples, we performed a permutation test as in the scProportionTest R package (), using the permutation_test function with the default parameters (n=1, 000 permutation tests). The test calculates the relative difference in abundance of the spots annotated by a cluster between the two different conditions, providing the log2FC and false discovery rate (FDR). The permutation_plot function was used to visualize the FC between recurrence status for all clusters.
2.3.3 Cell type annotation at spot level by deconvolution
Visium spot deconvolution was performed using the Conditional Autoregressive Deconvolution (CARD) R package (). This software estimates cell-type composition within spatial transcriptomics spots by integrating cell-type-specific transcriptional profiles derived from a single-cell reference dataset with a conditional autoregressive model that accounts for spatial correlation between neighboring spots. Publicly available scRNA-seq data, generated from 37 human HNC specimens, was obtained from the National Center of Biotechnology Information’s Gene Expression Omnibus (GSE181919) () and used as reference. This scRNA-seq dataset consisted of samples with or without HPV, from different tissue origins and from primary and metastatic patients, matching the diversity of our dataset.
Raw count data from both spatial transcriptomics and single-cell data was used to run CARD. The cell type with the highest predicted proportion was assigned as the dominant identity for each spot. CARD was also applied to impute single-cell resolution mappings from spatial transcriptomics data. This was achieved by constructing a refined spaced grid and estimating new cell locations through calculating a conditional mean established between the original and imputed cell coordinates. This CARD function was used to generate single-cell mapped spatial data for each sample.
2.3.4 Differential analysis of mRNA expression data
Differential expression (DE) analysis was performed between recurrent and non-recurrent samples in regions labeled by a pathologist as tumor, stroma and stroma-tumor infiltrating lymphocytes (TIL). DE analysis was also conducted for tumor p16-positive (n=4) and p16-negative (n=4) FFPE sections separately. By focusing the DE analysis specifically on the spots in the malignant compartment, we distill the non-specific signal from other regions of the tissue. Biological samples were randomly partitioned into three distinct groups of spots to generate triplicate pseudo-samples per patient. The exceptions were samples NR1 and NR4, which were aggregated as single pseudo-samples due to their limited number of tumor spots. For each partition, transcriptomic data were averaged across all constituent spots. This approach yielded a total of 26 pseudo-bulk profiles for differential expression (DE) analysis. Differential expression analysis was then performed on these aggregated profiles using the edgeR (v3.36.0) R package (), which is highly robust to the sparsity and variable library sizes inherent to spatial pseudo-bulk data.
A limitation of this approach is that the three pseudo-samples derived from the same patient are not independent biological replicates. Because patient-level correlation was not modelled within a mixed-effects framework, the resulting p-values are affected by an inflated effective sample size. Strict patient-level pseudo-bulk aggregation (i.e., one pseudo-bulk per patient) was not used due to insufficient degrees of freedom, which arises from the small cohort size (n = 10), in the context of the high baseline clinical and anatomical heterogeneity of the samples. Consequently, this initial analysis was used strictly as an exploratory feature-selection tool to rank candidate genes for subsequent independent validation, rather than as a definitive statistical assessment of differential expression. We also implemented a DE analysis approach taking into account patient-level variation using linear mixed-effects model as described below.
2.3.5 Gene set enrichment and GO term analysis on DE genes specific to cell types
Gene set enrichment analysis (GSEA) and Gene Ontology (GO) analysis were applied to find enriched pathways among up and down regulated genes in the recurrent tumor samples. For GO analysis the list of differentially expressed genes was filtered based on a p-value (p<0.05) and log2fold change (log2FC) >1.0 to obtain biologically meaningful pathways significantly changing between recurrence and nonrecurrence. For GSEA, the log2FC-ranked genes were used to calculate the enrichment score, using all expressed genes in the DE analysis (N = 13, 266) as the background genes (). Significance was determined using the BH-adjusted p-value (p<0.05) (). The following R packages were used to perform GO analysis: org.Hs.eg.db (v3.14.0), DOSE (3.24.2) (), enrichplot (1.18.3), clusterProfiler (v4.2.2) () and ReactomePA (1.42.0) () and fgsea (1.24.0) ().
2.3.6 Identification of candidate genes associated with tumor recurrence using DE genes
To identify top candidate genes for tumor recurrence from a large pool of differentially expressed genes, we implemented a systematic prioritization strategy leveraging both statistical thresholds and spatial expression patterns.
First, genes were initially filtered based on an effect size threshold of |log2FC| > 1.5. From this subset, the top 10 genes were selected for further evaluation based on a combination of fold change magnitude and statistical significance. Selecting genes with a higher log2FC increases the likelihood of capturing robust biological changes that translate to the protein level, while prioritizing the lowest P values minimizes the false discovery rate.
Second, we leveraged the unique spatial resolution of our data to map gene expression directly back onto the corresponding H&E tissue images. Using the SpatialDimPlot function from Seurat v4.1.0, we assessed whether candidate gene expression was consistently and uniformly upregulated across malignant domains in all recurrent samples. This spatial filtering step was critical to prevent single sample expression outliers from biasing our pseudo bulk statistical significance, which represents an essential safeguard given our discovery cohort size. Genes that demonstrated consistent upregulation across multiple recurrent patients, rather than localized or sample specific spikes, were prioritized as candidate biomarkers for downstream protein validation.
To account for the intra-patient correlation inherent in partitioning spatial spots, the four candidate genes were re-evaluated using a linear mixed-effects framework. Spot-level gene expression was modelled using the fixed effect of the clinical recurrence status, with patient identity included as random intercepts. This approach controls for the possible inflation of effective sample sizes by treating individual patients as the primary unit of biological replication.
2.3.7 Validating transcriptionally-defined candidate recurrence genes through immunohistochemical staining of proteins
For IHC staining of candidate proteins, FFPE cancer tissue slides (for FMO2) or tissue microarray slides (for MARCO, PNPLA3 and PAEP) were deparaffinized in xylene and rehydrated in ethanol series, and heat induced antigen retrieval was performed with either Sodium Citrate Buffer (10mM Sodium Citrate, 0.05% Tween 20, pH 6.0, at 100C for 10 minutes, for FMO2), or Diva pH6 decloaker buffer (Cat# DV2004, Biocare medical, CA, USA, at 100C for 20 minutes, for PNPLA3, PAEP and MACO). Sections were then incubated with the following antibodies: Anti-MARCO polyclonal antibody (Cat#PA5-64134, Invitrogen, 1:1000, 90 mins at RT), anti-PNPLA3 polyclonal antibody (Cat#PA5-18901, Invitrogen, 1:500, 90 mins at RT), anti-PAEP polyclonal antibody (Cat#PA5-54152, Invitrogen, 1:500, 90 mins at RT) and anti-FMO2 polyclonal antibody (Cat#ab95977, Abcam, 1:500, 120 mins at RT). Then, the samples were rinsed twice with PBS for 5 minutes, and endogenous peroxidase activity was blocked in 2.5% hydrogen peroxide for 10 minutes at room temperature. The samples were rinsed again in PBS, and incubated for 1 hour at RT with secondary antibodies: MACH – 1 universal polymer (Biocare medical, 30 minutes at RT, for MARCO and PAEP), Goat Probe and polymer (Biocare Medical, 20 minutes at RT for PNPLA3) and biotinylated anti rabbit secondary antibody (Thermofisher, 2 hours at RT, for FMO2). Streptavidin-hourseradish peroxidase was added to samples stained for FMO2 (20 minutes at RT). Diaminobenzidine (Abcam) was used as substrate, and sections were counterstained with haematoxylin. Slides were imaged with a microscope (Olympus BX60) and processed by ImageJ. Cellular scoring was conducted on a semi-quantitative four-point scale by a qualified anatomical pathologist, assessing IHC intensity as follows: negative (0), mild positivity (1), moderate positivity (2), and strong positivity (3).
2.3.8 Validating candidate recurrence genes with external TCGA bulk RNA-seq data
Bulk RNA-seq data were obtained from TCGA-HNSC using the following criteria: (1) diagnosis of primary HNSC without prior treatment, (2) availability of RNA-seq data from the primary tumor sample, and (3) complete data for age, gender, race, alcohol use, and smoking exposure. Using additional patient follow-up information from the TCGA database, patients were classified into two groups: the recurrence group (patients who experienced tumor progression or recurrence during follow-up) and the non-recurrence group (patients with no evidence of progression or recurrence and who were tumor-free at last contact). Age was further categorized into four groups: 0–39, 40–59, 60–79, and ≥80 years. Details on the dataset are provided in Supplementary Table 2.
Because the TCGA data represents deeply sequenced, standard bulk RNA counts, gene expression analysis was performed using DESeq2 (), comparing recurrence and non-recurrence groups while adjusting for all relevant confounding factors mentioned above.
2.3.9 Validating candidate recurrence genes with independent spatial transcriptomics datasets
Marker gene expression was analyzed using two additional public Visium datasets from HNC. The first dataset comprised paired samples collected before and after cancer, and the second consisted of a sample from a patient with no evidence of recurrence after a two-year MRI follow-up [relapse, details at ()]. In addition, Xenium data from an HNC patient with unknown recurrence status were analyzed to validate the cell-type specificity of PAEP. Gene expression values were normalized using the SCT method (). Cell type identities for the Xenium dataset were determined via automated label transfer using a reference single-cell RNA-sequencing (scRNA-seq) dataset from Quah et al. (). Following log-normalization of the Xenium transcript counts using the NormalizeData function, we implemented the Seurat integration workflow. Specifically, ‘anchors’ were identified between the reference and query datasets using CCA, and cell type labels were projected onto the Xenium data using the TransferData function. This approach assigns a prediction score to each cell, ensuring that spatial cell type assignments are consistent with established transcriptomic profiles.
2.3.10 Public cohort survival analysis
To evaluate the prognostic relevance of PAEP expression in a larger clinical context, a formal survival analysis was performed using the KM-plotter web portal (Győrffy, 2024), which integrates public transcriptomic and clinical datasets for head and neck squamous cell carcinoma. Patient cohorts were stratified into high- and low-expression groups based on PAEP mRNA levels using the auto-select best cutoff feature. Overall survival (OS) was evaluated at two distinct clinical endpoints: a short-term, high-risk window restricted to 30 months follow-up, and an unconstrained long-term follow-up extending up to 200 months. Hazard ratios (HR) with 95% confidence intervals (CI) and log-rank p-values were calculated automatically by the portal interface.
3 Results
3.1 Cell type and cell-type specific gene signature changes between recurrent and non-recurrent patients
To ensure the reliability of the spatial analysis, we first performed data preprocessing and normalization. Given the high degree of tissue complexity and structural diversity in HNC, as documented by expert pathologist annotations (Supplementary Figures 1A–C), robust batch correction was essential to mitigate technical bias and ensure that the resulting signals reflected true biological variation across samples (Supplementary Figures 2A–D). Following these quality control measures, we combined unsupervised analyses with histological annotation to organize the intricate histological landscape into biologically meaningful molecular units. We performed unsupervised clustering of all Visium spots based on gene expression similarities, which identified five distinct clusters (Figures 2A, B).
Figure 2
Spatial mapping revealed that Clusters 0 and 2 were predominantly localized to malignant areas, showing high spatial concordance with pathologist-defined tumor regions (Figure 2C; Supplementary Figure 1). This alignment was quantitatively confirmed, with 70.7% of spots in Cluster 0 and 76.7% in Cluster 2 falling within annotated tumor boundaries (Supplementary Figure 3). In contrast, Cluster 1 exhibited a heterogeneous distribution, spanning both stromal and immune-enriched regions. This mixed identity was reflected in its enriched GO terms, which broadly combined stromal terms (e.g., collagen fibril organization) with general immune functions (e.g., complement activation) (Supplementary Figure 4A).
To resolve this heterogeneity, Cluster 1 was subclustered into two distinct spatial domains: Subcluster 1_0 and Subcluster 1_1 (Supplementary Figure 4B). This separation successfully deconvoluted the biological signals: Subcluster 1_0 localized to immune-rich areas and displayed distinct innate and adaptive immune signatures (e.g., regulation of B cell activation), whereas Subcluster 1_1 aligned with the stroma and was enriched for connective tissue functions (Supplementary Figures 3, 4C). Notably, the specific adaptive immune pathways detected in Subcluster 1_0 were obscured when analyzing Cluster 1 as a whole, highlighting the necessity of subclustering to capture the distinct immune microenvironment.
To identify spatial domains associated with clinical outcome, we applied permutation testing and bootstrap analysis (Figure 2D). We found that Subcluster 1_1, Cluster 3, and Cluster 4 were significantly enriched in non-recurrent samples, whereas Cluster 2 was enriched in recurrent samples. Subcluster 1_1 (log2FC < -0.58), identified above as a stromal compartment, showed significant enrichment in the non-recurrent group, consistent with a role for stromal architecture in tumor containment. Cluster 3 showed the strongest nominal enrichment (log2FC = –4.44) and was defined by muscle-related terms such as ‘striated muscle contraction’ (Supplementary Figures 4D, E). Pathologist annotations localized this cluster predominantly to samples NR4 and NR6, both of which are mandible tumors, suggesting this signature reflects adjacent muscle tissue captured in these specific resections rather than a tumor-intrinsic mechanism. Cluster 4 (log2FC = –1.62) also favored non-recurrent samples. Although this cluster contained relatively few spots, it mapped predominantly to immune-rich regions and was associated with the negative regulation of leukocyte migration. Conversely, Cluster 2, which showed enrichment in recurrent samples (log2FC > 0.58), mapped largely to tumor and normal mucosal regions. This cluster was characterized by nine enriched GO terms, including immune-related processes (e.g., acute inflammatory response, regulation of humoral response) as well as pathways linked to cell proliferation (Supplementary Figures 4D, E).
To further resolve the cellular composition within individual Visium spots and characterize the heterogeneity of the tumor microenvironment, we performed transcriptomic deconvolution, which suggested a divergent immune landscape between outcomes. Non-recurrent samples displayed a trend toward higher immune infiltration, particularly regarding T and B cell populations (Figure 3A). Quantitatively, both the average proportion of the two lymphocytes per spot (Supplementary Figure 5A) and the frequency of immune-dominated spots (Supplementary Figure 5B) were elevated in the non-recurrent group relative to recurrent samples. However, these differences did not reach statistical significance (Wilcoxon test, p > 0.05), likely reflecting the limited statistical power of the cohort size. Despite this, the spatial distribution of these estimated cell types aligned closely with pathologist-annotated regions, reinforcing the biological validity of the deconvolution signal (Figures 3B, C; Supplementary Figure 6A). For instance, the ‘myocyte’ signal observed within histologically benign salivary tissue accurately captures endogenous myoepithelial cells, which physiologically surround salivary acini and share strong transcriptional overlap with muscle lineages, a biological nuance consistent with the cluster definitions in our reference atlas. For sanity check, we confirmed that the tumor compartment demonstrates clear and robust enrichment for epithelial markers, including KRT5, KRT14, TP63, and EPCAM. Conversely, the stroma regions display expected enrichment for the canonical mesenchymal/stromal marker VIM and VWF (Supplementary Figure 6B). However, because these spatial resolution maps are computationally inferred by CARD rather than directly observed, they represent statistical cell-type proportions within each multicellular spot.
Figure 3
Spatial analysis highlighted a key distinction in immune topography: within non-recurrent tumors, T cells and B cells were frequently observed colocalizing in the same histological niches. In contrast, recurrent samples tend to lack both lineages (Supplementary Figures 7A–D), suggestive of an immune-excluded phenotype. This trend was further confirmed by quantitative density profiling: while recurrent tumors displayed a single density peak near zero for both cell types, non-recurrent samples revealed a distinct secondary peak, representing a specific subpopulation of tumor spots containing approximately 5% T cells and 1% B/plasma cells (Supplementary Figure 8). These results indicate that the coordinated infiltration of both T and B cells into the tumor core is a specific feature of non-recurrent disease, distinguishing it from the immunologically quiescent environment of recurrent tumors.
3.2 Spatial DE analysis of cell types identified genes linked to recurrence in HNC
To identify candidate biomarkers associated with HNC recurrence, we focused our comparative analysis specifically on the malignant compartment (Figure 4A). We isolated spatial gene expression data for spots located within pathologist-annotated tumor regions (Supplementary Figure 2A) and performed DE analysis to contrast the transcriptional profiles of recurrent versus non-recurrent tumors. We compared 13, 266 genes from 26 pseudo-samples (see Methods; also refer to the linear mixed-effects approach) and identified 271 and 265 genes significantly up- and downregulated in recurrent tumor spots compared to non-recurrent tumor spots (P<2.02×10-3), respectively (Figure 4B).
Figure 4
To explore functional pathways associated with tumor recurrence, we first performed GSEA using the Reactome database (). Non-recurrence samples showed upregulation of MAPK6/MAPK4 and PCP/CE pathways (Supplementary Figure 9), consistent with the context-dependent pro- and anti-oncogenic roles of these kinases () and the PCP/CE pathway (). Broad GO enrichment analysis of tumor DEGs did not identify significant terms, likely due to disease heterogeneity.
To address this, recurrence DEGs were stratified by p16 status. In p16-positive tumor regions, pseudo-bulk analysis identified 1, 486 upregulated and 1, 677 downregulated DEGs (filtering threshold P < 1.12×10-2), whereas p16-negative regions had 231 upregulated and 404 downregulated DEGs in recurrent samples (P < 2.14×10-3) (Supplementary Figures 10A, B). Stratification revealed significant GO term enrichment (P < 0.05), including: p16-positive upregulated DEGs were enriched for collagen fibril organization and xenobiotic metabolism, while downregulated DEGs were associated with metabolic changes. In p16-negative recurrent tumors, downregulated DEGs were enriched for viral regulation, IFN-β responses, and apoptosis, while no upregulated DEGs were found (Supplementary Figure 10C). Clustering of enriched GO terms suggested that p16-positive and p16-negative tumors may use distinct pathways to drive recurrence (Supplementary Figure 10D).
Analyses of stromal regions (Supplementary Figure 11) and stromal–TIL–tumor interfaces (Supplementary Figure 12) identified additional genes significantly upregulated in recurrent tumors. These recurrence-associated markers showed enrichment in proliferation, protein processing, and immune regulatory pathways. Although limited by the small cohort size, these findings in the stroma and interface compartments are broadly consistent with the immune-suppressed and proliferative profile observed within the recurrent tumor regions.
3.3 Identification and validation of candidate genes associated with tumor recurrence
To find top recurrence candidate biomarkers, we filtered over 200 upregulated DEGs (Supplementary Table 3) in recurrent tumor samples, narrowing it down to 104 DEGs with a fold change (FC) threshold (|log2FC|>1.5) (Figure 4B). The top 10 genes, ranked by fold-change and p-value, underwent systematic spatial evaluation to verify their distribution. Candidates were defined as ‘robust’ if they demonstrated consistent upregulation within the tumor compartments of at least 50% of the recurrent samples, ensuring that the prioritized markers represented the broader recurrent phenotype rather than isolated outliers.
For instance, CWH43 was excluded despite a high fold-change (log2FC = 3.97) because its expression was restricted almost exclusively to sample R2. In contrast, FMO2, PNPLA3, PAEP, and MARCO emerged from this curation as the most robust candidates for further validation, though they exhibited distinct patterns of heterogeneity (Figures 4C, D). Specifically, PAEP and PNPLA3 demonstrated strong, consistent upregulation across samples R1, R2, and R4. Conversely, the signals for FMO2 and MARCO were predominantly driven by sample R3, reflecting a distinct transcriptional profile in this tumor. Importantly, while the expression of FMO2 and MARCO appeared lower in R1, R2, and R4 compared to R3, it remained consistently elevated relative to the non-recurrent baseline. By retaining both sets of markers, we aimed to account for inter-tumor heterogeneity and ensure detection across diverse recurrence phenotypes.
Protein-level validation was performed by IHC on tumor regions from 4 patients with recurrence and 6 without recurrence. Representative images are shown in Figures 5A–C. PNPLA3 expression was higher in the recurrence group (scores 1–3, median 2) compared to the non-recurrence group (scores 0–3, median 1). MARCO expression was also elevated in patients with recurrence (scores 0–3, median 1.5) relative to those without recurrence (scores 0–2, median 0). PAEP expression was low overall, with scores of 0–1 (median 0) in the recurrence group, and uniformly negative in the non-recurrence group.
Figure 5
To further assess the association between candidate gene expression and tumor recurrence, we examined published spatial transcriptomics data comparing recurrent and non-recurrent HNC samples. Among the four markers, PAEP exhibited the most distinct association with recurrence status (Figure 5D; Supplementary Table 4). Notably, distinct PAEP-positive spots (highlighted in orange) were exclusively detected in the baseline primary tumors of patients who subsequently relapsed (‘Recurrent’ group). In contrast, PAEP expression was absent in both the matched post-treatment recurrences (‘Recurrent treated’) and the primary tumors of patients who remained disease-free (‘Non-recurrent’). This specific localization to the pre-treatment state of high-risk patients suggests that PAEP may serve as a candidate biomarker associated with recurrence risk. Similarly, PNPLA3 showed elevated expression specifically in the pre-relapse primary tumors compared to the other clinical groups.
For generalization, we also analyzed bulk RNA-seq data from 60 recurrence tumor samples and 116 non-recurrence tumor samples obtained from the TCGA database. After statistically adjusting for potential confounding factors, including primary tumor anatomical site, age, gender, race, and history of exposure to tobacco and alcohol, the expression of PAEP remained significantly elevated in the recurrence group compared to the non-recurrence group log2FC = 0.88, p = 0.007; Figure 5E), supporting its potential role as a biomarker for recurrence. Because this analysis was designed specifically to test the directional hypothesis of our four pre-selected candidate genes rather than to perform an exploratory genome-wide screen, unadjusted p-values were utilized. This association was not observed for the other three genes examined, which may be due to the lack of cell type specificity inherent in bulk RNA-seq data (i.e., expression signals from tumor cells may be diluted by non-tumor cells present in the sample). In contrast, PAEP expression in HNC is largely confined to cancer cells, as indicated by an external Xenium data which shows PAEP expression at single-cell resolution (Figures 5F, G). This aligns with Human Protein Atlas reports showing minimal PAEP expression in non-glandular and non-luminal normal cell types, underscoring the value of single-cell and spatial transcriptomics for identifying cell type–specific markers. Interestingly, plasma cells and T cells emerged as the next most prominent PAEP-expressing populations (Figure 5G), some of which can have higher PAEP expression than cancer cells (Supplementary Figure 13), suggesting a potential functional role for PAEP within immune compartments. Leveraging spatial transcriptomics, we observed a strong negative colocalization of PAEP with B cells and/or plasma cells in nine of ten HNC samples (Figure 6), whereas little spatial association was detected with T cells (Supplementary Figure 14).
Figure 6
To assess the relevance of these spatial findings to clinical outcomes, we analyzed PAEP expression within an integrated public HNC cohort via the KM-plotter analytical database. Given that a 30-month threshold effectively captures the majority of high-risk HNC recurrence events (Haring, 2023), we evaluated survival within this clinically relevant window. High PAEP expression was significantly associated with reduced survival within this timeframe (HR = 1.57 [1.14 to 2.15], logrank P = 0.0049; Supplementary Figure 15A), whereas the prognostic signal was attenuated over long-term follow-up (HR = 1.18, logrank P = 0.27; Supplementary Figure 15B). This association supports PAEP as a potential marker of early treatment failure, providing clinical context that aligns with our recurrence-focused spatial discovery cohort. Because comprehensive recurrence-free survival data were unavailable in this public dataset, overall survival was utilised as a proxy metric. This endpoint includes all-cause mortality, which may introduce confounding variables unrelated to disease recurrence. Additionally, the automated cut-off selection used to stratify expression groups introduces statistical optimism and may overestimate the significance of the survival difference. These retrospective results represent a baseline characterization of a candidate feature rather than a validated clinical prognostic tool.
To evaluate the statistical robustness of the four selected candidate genes against intra-patient correlation, we re-analysed the spatial spot data using a linear mixed-effects model with patient identity specified as a random effect. Under this framework, PAEP emerged as the only candidate to retain statistical significance (p-value threshold 0.05), demonstrating a positive association with the recurrent phenotype (effect size = 0.107, p = 0.048, Supplementary Table 5). The statistical signals for MARCO, PNPLA3, and FMO2 did not meet the threshold for significance after accounting for patient-level variance. This computational outcome aligns with our downstream validation experiments.
4 Discussion
This study produced a valuable data resource of spatially resolved transcriptomics to contrast baseline tumor characteristics between multiple recurrent and non-recurrent HNC samples. In contrast to conventional bulk RNA-seq, spatial transcriptomics enables the resolution of gene expression patterns within defined tumor compartments, thereby enhancing biomarker discovery specificity by reducing noise from surrounding non-tumor signals. By separately interrogating tumor and stromal regions annotated by the pathologist, we provide high-resolution insights into the cellular and molecular determinants of recurrence risk.
Spatial analysis suggested a trend toward a higher abundance of immune cells in non-recurrent tumors, with a clear enrichment of B cells and plasma cells. While the prognostic value of T-cell infiltration, especially CD8+ T cells, is well established across solid tumors, including HNC (–), the contribution of B cells and plasma cells has been less clear. Our exploratory findings align with emerging evidence that B cells and plasma cells may enhance antitumor immunity (, ). Patients without recurrence exhibited a trend toward greater immune infiltration within malignant regions, suggesting that coordinated T- and B-cell activity may underlie a more effective endogenous immune response.
Notably, reduced B-cell and plasma cell infiltration was not uniformly associated with recurrence. For example, sample R3 recurred despite high B-cell abundance, possibly reflecting its p16-positive status, as distinct pathway enrichments were observed in p16-positive tumors. This highlights the complexity of recurrence, which is unlikely to be explained by a single factor. Larger cohorts will be needed to identify additional drivers with statistical power, though this comes with the trade-off of balancing depth with breadth of analysis.
Through combined statistical testing and careful evaluation of robustness, we identified tumor-specific DEGs most strongly associated with recurrence. FMO2, PAEP, MARCO, and PNPLA3 emerged as top candidates. Of these, three exhibited an observational pattern of elevated protein expression in the recurrent group, with further evidence provided by external datasets and published literature.
Among these, PAEP was the most notable and statistically significant in the small discovery cohort. PAEP encodes glycodelin, a secreted glycoprotein with established roles in reproduction, including the prevention of premature capacitation () and the provision of immunoprotection during implantation and placentation (). Beyond these functions, glycodelin exerts potent immunomodulatory effects. It modulates T-cell responses by raising the activation threshold through interactions with surface molecules such as CD45 (and possibly CD7), thereby skewing cytokine production toward a Th2 profile and inducing apoptosis in activated lymphocytes via a caspase-dependent mitochondrial pathway (, ). Glycodelin also suppresses B-cell proliferation and alters monocyte cytokine secretion, collectively fostering an immunosuppressive environment critical for maternal tolerance during pregnancy (). Literature further suggests that glycodelin’s N-glycosylation–derived oligosaccharide chains may engage CD22 on B cells to mediate immunosuppressive activity (), though this interaction has not yet been validated in HNC.
In our analysis, PAEP expression in tumor regions showed an inversely correlated trend with B-cell and plasma cell abundance, suggesting an immunomodulatory role that could weaken humoral immune responses and contribute to recurrence. This finding is consistent with prior reports of PAEP upregulation in immune-resistant metastatic tumors and correlations with immune infiltration and checkpoint expression across cancer types (). To our knowledge, this spatial comparison suggests the new evidence implicating PAEP in shaping the immune microenvironment of HNC, providing a potential mechanistic link between PAEP expression, tumor immune evasion, and recurrence risk.
Other recurrence-associated DEGs also point to immune modulation. FMO2, previously linked to poor prognosis in ovarian and gastric cancers (–59), was enriched in activated CAFs and correlated with tumor-associated macrophage infiltration (58). MARCO, a macrophage scavenger receptor, was overexpressed in recurrent HNC and is known to promote tumor progression and immune evasion, with anti-MARCO antibodies showing preclinical promise (, 60–63). PNPLA3 remains poorly characterized in cancer; while certain variants predispose to liver cancers and metabolic liver disease (64, 65), its role in extrahepatic tumors remains unclear (66). For these genes, associations with HNC recurrence could not be consistently validated across datasets, likely reflecting the biological complexity of recurrence and the limitation of sample size.
A limitation of this study is the relatively restricted cohort size, consisting of 10 samples (4 recurrent and 6 non-recurrent) distributed across four slide captures. While this sample size limits the statistical power to robustly validate subtle transcriptomic differences or generalize findings across heterogeneous HNC populations, the study was designed to prioritize spatial resolution over cohort breadth. By generating high-dimensional data within defined histological compartments, we were able to identify robust molecular signals in the tumor microenvironment that might be obscured in larger, bulk RNA-seq cohorts.
External datasets were used to validate the findings. Furthermore, the small cohort size combined with high baseline clinical and anatomical heterogeneity precluded the use of strict patient-level pseudo-bulk aggregation due to insufficient degrees of freedom. Consequently, our discovery-phase differential expression analysis relied on intra-tumour pseudo-samples, which introduces intra-patient correlation and likely understates the initial p-values. Additionally, we performed linear mixed effects analysis, simultaneously taking into account inter-patient and intra-patient variation, to provide a comprehensive view of possible statistical approaches.
To mitigate the constraint by a small sample size, the discovery outputs were treated strictly as hypothesis-generating candidate features, with the statistical burden of proof shifted to our independent, patient-level external validation cohort (n = 176) where clinical confounders were controlled. While multivariate analysis of this large-scale external cohort demonstrates that the upregulation of molecular features like PAEP occurs independently of anatomical sub-locale, future site-matched spatial transcriptomic cohorts with larger sample sizes are required to decouple generalised recurrence-associated immune architectures from baseline tissue variations.
The linear mixed-effects framework was applied to take into account patient-level variance as applied to assess the four selected candidates. This re-evaluation confirmed PAEP as the top candidate marker maintaining statistical significance under a mixed-effects framework, whereas MARCO, PNPLA3, and FMO2 did not retain significance once patient-level variation was controlled. This computational result explains our downstream empirical findings, where PAEP was the top candidate feature to receive reproducible support in protein-level assessment and external cohorts.
Our findings demonstrate that spatially resolved transcriptomics can suggest cellular and molecular features associated with HNC recurrence, proposing the roles of both tumor-intrinsic markers and the immune microenvironment. By identifying PAEP as a candidate feature alongside descriptive observations regarding the contribution of B cells and plasma cells, this work demonstrates the utility of spatially informed approaches for candidate feature discovery. With continued advances in image-based ST and AI-driven analysis (67), such markers hold promise for the development of diagnostic tools that can guide recurrence risk stratification and personalized treatment in HNC.
Statements
Data availability statement
The raw human spatial transcriptomic data presented in this study have been deposited in the University of Queensland eSpace repository under DOI: https://doi.org/10.48610/5842caa. These data are not publicly accessible because they contain potentially identifiable human research data and are subject to ethical, institutional, and participant-consent restrictions. Requests for access should be submitted through the UQ eSpace repository and will be considered following approval by the relevant data custodian and institutional ethics requirements. The processed expression matrices, associated metadata, and other derived datasets supporting the conclusions of this article are publicly available in Zenodo at https://doi.org/10.5281/zenodo.18238293. The complete computational analysis pipeline, including preprocessing and spatial analysis scripts, is publicly available on GitHub at Available at: https://github.com/BiomedicalMachineLearning/RecurrenceHN.
Ethics statement
The studies involving humans were approved by Metro South Human Research Ethics Committee. The studies were conducted in accordance with the local legislation and institutional requirements. The participants provided their written informed consent to participate in this study.
Author contributions
QN: Investigation, Conceptualization, Funding acquisition, Resources, Methodology, Supervision, Writing – review & editing, Data curation, Writing – original draft. CW: Formal analysis, Data curation, Writing – original draft, Investigation, Writing – review & editing. YC: Writing – review & editing, Formal analysis, Methodology, Writing – original draft, Validation, Visualization, Investigation, Data curation. KN: Data curation, Investigation, Formal analysis, Writing – review & editing, Writing – original draft, Visualization. XZ: Formal analysis, Writing – original draft, Data curation, Methodology, Validation, Writing – review & editing. GH: Methodology, Data curation, Writing – original draft, Writing – review & editing. AC: Visualization, Writing – original draft, Methodology, Writing – review & editing, Formal analysis. OM: Methodology, Writing – review & editing, Writing – original draft, Validation, Visualization, Data curation, Formal analysis. PP: Methodology, Writing – review & editing, Formal analysis, Writing – original draft. PL: Writing – review & editing, Data curation, Writing – original draft, Methodology. JG: Resources, Writing – review & editing, Methodology, Writing – original draft. SV: Writing – review & editing, Writing – original draft, Resources, Methodology. LK: Resources, Writing – original draft, Methodology, Writing – review & editing. TT: Resources, Writing – original draft, Methodology, Writing – review & editing. OB: Writing – review & editing, Methodology, Writing – original draft, Resources. CP: Writing – review & editing, Investigation, Conceptualization, Writing – original draft, Supervision, Funding acquisition, Resources, Methodology.
Funding
The author(s) declared that financial support was received for this work and/or its publication. CP is currently receiving funds from the National Health and Medical Research Council (APP 2002576 and APP 2012560), the Garnett Passe and Rodney Williams Foundation, NIH R21 and the Metro North Collaborative Grant Scheme. QN is supported by the National Health and Medical Research Council Ideas Grant (APP 2001514) and Investigator Grant (GNT2008928).
Acknowledgments
The authors thank RBWH, cancer care services clinical coordinators for patient recruitment. The results published here are in part based upon data generated by the TCGA Research Network: https://www.cancer.gov/tcga. The work was also supported by the Australian Cancer Research Foundation Center for Optimized Cancer Therapy and the Queensland Institute of Medical Research Berghofer Spatial Tissue and Artificial Intelligence Research Center.
Conflict of interest
The author(s) declared that this work was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Generative AI statement
The author(s) declared that generative AI was not used in the creation of this manuscript.
Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.
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.2026.1713732/full#supplementary-material
Generative AI statement
The author(s) declared that generative AI was not used in the creation of this manuscript.
References
1
SungHFerlayJSiegelRLLaversanneMSoerjomataramIJemalAet al. Global cancer statistics 2020: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA: A Cancer J For Clin. (2021) 71:209–49. doi: 10.3322/caac.21660
2
ArgirisALiSSavvidesPOhrJPGilbertJLevineMAet al. Phase III randomized trial of chemotherapy with or without bevacizumab in patients with recurrent or metastatic head and neck cancer. J Clin Oncol Off J Am Soc Clin Oncol. (2019) 37:3266–74. doi: 10.1200/jco.19.00555
3
ModyMDRoccoJWYomSSHaddadRISabaNF. Head and neck cancer. Lancet. (2021) 398:2289–99. doi: 10.1016/s0140-6736(21)01550-6
4
KulasingheAKennyLPerryCThieryJPJovanovicLVelaIet al. Impact of label-free technologies in head and neck cancer circulating tumour cells. Oncotarget. (2016) 7:71223–34. doi: 10.18632/oncotarget.12086
5
SchmidtHKulasingheAPerryCNelsonCPunyadeeraC. A liquid biopsy for head and neck cancers. Expert Rev Mol Diagn. (2016) 16:165–72. doi: 10.1586/14737159.2016.1127758
6
HillerJGPerryNJPoulogiannisGRiedelBSloanEK. Perioperative events influence cancer recurrence risk after surgery. Nat Rev Clin Oncol. (2018) 15:205–18. doi: 10.1038/nrclinonc.2017.194
7
ZhangXEkanayake WeeramangeCHughesBGMVasaniSLiuZYWarkianiMEet al. Application of circulating tumour cells to predict response to treatment in head and neck cancer. Cell Oncol (Dordrecht). (2022) 45:543–55. doi: 10.1007/s13402-022-00681-w
8
KulasingheAKennyLPunyadeeraC. Circulating tumour cell PD-L1 test for head and neck cancers. Oral Oncol. (2017) 75:6–7. doi: 10.1016/j.oraloncology.2017.10.011
9
HuangXLeoPJonesLDuijfPHGHartelGKennyLet al. A comparison between mutational profiles in tumour tissue DNA and circulating tumour DNA in head and neck squamous cell carcinoma – a systematic review. Mutat Res Rev Mutat Res. (2024) 793:108477. doi: 10.1016/j.mrrev.2023.108477
10
HuangXDuijfPHGSriramSPereraGVasaniSKennyLet al. Circulating tumour DNA alterations: emerging biomarker in head and neck squamous cell carcinoma. J BioMed Sci. (2023) 30:65. doi: 10.1186/s12929-023-00953-z
11
TurajlicSSottorivaAGrahamTSwantonC. Resolving genetic heterogeneity in cancer. Nat Rev Genet. (2019) 20:404–16. doi: 10.1038/s41576-019-0114-6
12
GuoMPengYGaoADuCHermanJG. Epigenetic heterogeneity in cancer. biomark Res. (2019) 7:23. doi: 10.1186/s40364-019-0174-y
13
LambertAWPattabiramanDRWeinbergRA. Emerging biological principles of metastasis. Cell. (2017) 168:670–91. doi: 10.1016/j.cell.2016.11.037
14
FaresJFaresMYKhachfeHHSalhabHAFaresY. Molecular principles of metastasis: a hallmark of cancer revisited. Signal Transduction Tgt Ther. (2020) 5:28. doi: 10.1038/s41392-020-0134-x
15
GaoXLZhangMTangYLLiangXH. Cancer cell dormancy: mechanisms and implications of cancer recurrence and metastasis. OncoTargets Ther. (2017) 10:5219–28. doi: 10.2147/ott.s140854
16
FridmanWHPagèsFSautès-FridmanCGalonJ. The immune contexture in human tumours: impact on clinical outcome. Nat Rev Cancer. (2012) 12:298–306. doi: 10.1038/nrc3245
17
KansyBAConcha-BenaventeFSrivastavaRMJieH-BShayanGLeiYet al. PD-1 status in CD8+ T cells associates with survival and anti-PD-1 therapeutic outcomes in head and neck cancer. Cancer Res. (2017) 77:6353–64. doi: 10.1158/0008-5472.CAN-16-3167
18
SoYKByeonS-JKuBMKoYHAhnM-JSonY-Iet al. An increase of CD8+ T cell infiltration following recurrence is a good prognosticator in HNSCC. Sci Rep. (2020) 10:20059. doi: 10.1038/s41598-020-77036-8
19
ElmusratiAAPilboroughAEKhurramSALambertDW. Cancer-associated fibroblasts promote bone invasion in oral squamous cell carcinoma. Br J Cancer. (2017) 117:867–75. doi: 10.1038/bjc.2017.239
20
NewJArnoldLAnanthMAlviSThorntonMWernerLet al. Secretory autophagy in cancer-associated fibroblasts promotes head and neck cancer progression and offers a novel therapeutic target. Cancer Res. (2017) 77:6679–91. doi: 10.1158/0008-5472.can-17-1077
21
TakahashiHSakakuraKKawabata-IwakawaRRokudaiSToyodaMNishiyamaMet al. Immunosuppressive activity of cancer-associated fibroblasts in head and neck squamous cell carcinoma. Cancer Immunol Immunother. (2015) 64:1407–17. doi: 10.1007/s00262-015-1742-0
22
ChoiJ-HLeeB-SJangJYLeeYSKimHJRohJet al. Single-cell transcriptome profiling of the stepwise progression of head and neck cancer. Nat Commun. (2023) 14:1055. doi: 10.1038/s41467-023-36691-x
23
PanCIssaevaNYarbroughWG. HPV-driven oropharyngeal cancer: current knowledge of molecular biology and mechanisms of carcinogenesis. Cancers Head Neck. (2018) 3:12. doi: 10.1186/s41199-018-0039-3
24
Illumina. Evaluating RNA quality from FFPE samples (2016). Available online at: https://www.illumina.com/content/dam/illumina-marketing/documents/products/technotes/evaluating-rna-quality-from-ffpe-samples-technical-note-470-2014-001.pdf (Accessed July 25, 2026).
25
VoTPrakrithiPJonesKYoonSLamPYKaoY-Cet al. Benchmarking robust spatial transcriptomics approaches to capture the molecular landscape and pathological architecture of archived cancer tissues. J Pathol. (2025) 265:274–88. doi: 10.1002/path.6383
26
HaoYHaoSAndersen-NissenEMauckWMZhengSButlerAet al. Integrated analysis of multimodal single-cell data. Cell. (2021) 184:3573–3587.e29. doi: 10.1016/j.cell.2021.04.048
27
HafemeisterCSatijaR. Normalization and variance stabilization of single-cell RNA-seq data using regularized negative binomial regression. Genome Biol. (2019) 20:296. doi: 10.1186/s13059-019-1874-1
28
YuGWangLGHanYHeQY. clusterProfiler: an R package for comparing biological themes among gene clusters. Omics. (2012) 16:284–7. doi: 10.1089/omi.2011.0118
29
MillerSAPolicastroRASriramkumarSLaiTHuntingtonTDLadaikaCAet al. LSD1 and aberrant DNA methylation mediate persistence of enteroendocrine progenitors that support BRAF-mutant colorectal cancer. Cancer Res. (2021) 81:3791–805. doi: 10.1158/0008-5472.can-20-3562
30
MaYZhouX. Spatially informed cell-type deconvolution for spatial transcriptomics. Nat Biotechnol. (2022) 40:1349–59. doi: 10.1038/s41587-022-01273-7
31
RobinsonMDMcCarthyDJSmythGK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. (2010) 26:139–40. doi: 10.1093/bioinformatics/btp616
32
TimmonsJASzkopKJGallagherIJ. Multiple sources of bias confound functional enrichment analysis of global -omics data. Genome Biol. (2015) 16:186. doi: 10.1186/s13059-015-0761-7
33
BenjaminiYHochbergY. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Stat Soc: Ser B (Methodol). (1995) 57:289–300. doi: 10.1111/j.2517-6161.1995.tb02031.x
34
YuGWangL-GYanG-RHeQ-Y. DOSE: an R/Bioconductor package for disease ontology semantic and enrichment analysis. Bioinformatics. (2014) 31:608–9. doi: 10.1093/bioinformatics/btu684
35
YuG. Thirteen years of clusterProfiler. Innovation. (2024) 5(6):100722. doi: 10.1016/j.xinn.2024.100722
36
YuGHeQ-Y. ReactomePA: an R/Bioconductor package for reactome pathway analysis and visualization. Mol Biosyst. (2016) 12:477–9. doi: 10.1039/c5mb00663e
37
KorotkevichGSukhovVBudinNShpakBArtyomovMNSergushichevA. Fast gene set enrichment analysis. bioRxiv. (2021), 060012. doi: 10.1101/060012
38
LoveMIHuberWAndersS. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. (2014) 15:550. doi: 10.1186/s13059-014-0550-8
39
CauserATanXLuXMoseleyPTeohSMMolotkovNet al. Deep spatial-omics analysis of head & neck carcinomas provides alternative therapeutic targets and rationale for treatment failure. NPJ Precis Oncol. (2023) 7:89. doi: 10.1038/s41698-023-00444-2
40
QuahHSCaoEYSutejaLLiCHLeongHSChongFTet al. Single cell analysis in head and neck cancer reveals potential immune evasion mechanisms during early metastasis. Nat Commun. (2023) 14:1680. doi: 10.1038/s41467-023-37379-y
41
FabregatASidiropoulosKViteriGFornerOMarin-GarciaPArnauVet al. Reactome pathway analysis: a high-performance in-memory approach. BMC Bioinf. (2017) 18:142. doi: 10.1186/s12859-017-1559-2
42
KostenkoSDumitriuGMoensU. Tumour promoting and suppressing roles of the atypical MAP kinase signalling pathway ERK3/4-MK5. J Mol Signaling. (2012) 7:9. doi: 10.1186/1750-2187-7-9
43
GeorgoudakiAMProkopecKEBouraVFHellqvistESohnSÖstlingJet al. Reprogramming tumor-associated macrophages by antibody targeting inhibits cancer progression and metastasis. Cell Rep. (2016) 15:2000–11. doi: 10.1016/j.celrep.2016.04.084
44
van der LeunAMThommenDSSchumacherTN. CD8+ T cell states in human cancer: insights from single-cell analysis. Nat Rev Cancer. (2020) 20:218–32. doi: 10.1038/s41568-019-0235-4
45
NguyenNBellileEThomasDMcHughJRozekLViraniSet al. Tumor infiltrating lymphocytes and survival in patients with head and neck squamous cell carcinoma. Head Neck. (2016) 38:1074–84. doi: 10.1002/hed.24406
46
BalermpasPMichelYWagenblastJSeitzOWeissCRödelFet al. Tumour-infiltrating lymphocytes predict response to definitive chemoradiotherapy in head and neck cancer. Br J Cancer. (2014) 110:501–9. doi: 10.1038/bjc.2013.640
47
ShimizuSHiratsukaHKoikeKTsuchihashiKSonodaTOgiKet al. Tumor-infiltrating CD8+ T-cell density is an independent prognostic marker for oral squamous cell carcinoma. Cancer Med. (2019) 8:80–93. doi: 10.3802/jgo.2023.34.s2.is-020
48
LiuHLiZHanXLiZZhaoYLiuFet al. The prognostic impact of tumor-infiltrating B lymphocytes in patients with solid Malignancies: a systematic review and meta-analysis. Crit Rev Oncol/Hematol. (2023) 181:103893. doi: 10.1016/j.critrevonc.2022.103893
49
LaumontCMBanvilleACGilardiMHollernDPNelsonBH. Tumour-infiltrating B cells: immunological mechanisms, clinical impact and therapeutic opportunities. Nat Rev Cancer. (2022) 22:414–30. doi: 10.1038/s41568-022-00466-1
50
KoistinenHKoistinenRDellAMorrisHREastonRLPatankarMSet al. Glycodelin from seminal plasma is a differentially glycosylated form of contraceptive glycodelin-A. Mol Hum Reprod. (1996) 2:759–65. doi: 10.1093/molehr/2.10.759
51
JulkunenMRutanenE-MKoskimiesARantaTBohnHSeppäläM. Distribution of placental protein 14 in tissues and body fluids during pregnancy. BJOG: Int J Obstet Gynaecol. (1985) 92:1145–51. doi: 10.1111/j.1471-0528.1985.tb03027.x
52
LeeCLChiuPCLamKKSiuSOChuIKKoistinenRet al. Differential actions of glycodelin-A on Th-1 and Th-2 cells: a paracrine mechanism that could produce the Th-2 dominant environment during pregnancy. Hum Reprod. (2011) 26:517–26. doi: 10.1093/humrep/deq381
53
LeeC-LLamEYFLamKKWKoistinenHSeppäläMNgEHYet al. Glycodelin-A stimulates interleukin-6 secretion by human monocytes and macrophages through L-selectin and the extracellular signal-regulated kinase pathway*. J Biol Chem. (2012) 287:36999–7009. doi: 10.1074/jbc.m112.385336
54
AlokAKarandeAA. The role of glycodelin as an immune-modulating agent at the feto–maternal interface. J Reprod Immunol. (2009) 83:124–7. doi: 10.1016/j.jri.2009.06.261
55
DellAMorrisHREastonRLPanicoMPatankarMOehningerSet al. Structural analysis of the oligosaccharides derived from glycodelin, a human glycoprotein with potent immunosuppressive and contraceptive activities. J Biol Chem. (1995) 270:24116–26. doi: 10.1074/jbc.270.41.24116
56
SchneiderMAGranzowMWarthASchnabelPAThomasMHerthFJFet al. Glycodelin: a new biomarker with immunomodulatory functions in non–small cell lung cancer. Clin Cancer Res. (2015) 21:3529–40. doi: 10.1158/1078-0432.ccr-14-2464
57
ChenYChenZTangYXiaoQ. The involvement of noncanonical Wnt signaling in cancers. Biomed Pharmacotherapy. (2021) 133:110946. doi: 10.1016/j.biopha.2020.110946
58
DolphinCTBeckettDJJanmohamedACullingfordTESmithRLShephardEAet al. The Flavin-containing Monooxygenase 2 Gene (FMO2) of humans, but not of other primates, encodes a truncated, nonfunctional protein. J Biol Chem. (1998) 273:30599–607. doi: 10.1074/jbc.273.46.30599
59
YuSYangRXuTLiXWuSZhangJ. Cancer-associated fibroblasts-derived FMO2 as a biomarker of macrophage infiltration and prognosis in epithelial ovarian cancer. Gynecol Oncol. (2022) 167:342–53. doi: 10.1016/j.ygyno.2022.09.003
60
LenhardMHeubleinSKunert-KeilCVrekoussisTLombaIDitschNet al. Immunosuppressive Glycodelin A is an independent marker for poor prognosis in endometrial cancer. BMC Cancer. (2013) 13:616. doi: 10.1186/1471-2407-13-616
61
ShiBChuJHuangTWangXLiQGaoQet al. The scavenger receptor MARCO expressed by tumor-associated macrophages are highly associated with poor pancreatic cancer prognosis. Front Oncol. (2021) 11:771488. doi: 10.3389/fonc.2021.771488
62
EisingerSSarhanDBouraVFIbarlucea-BenitezITyystjärviSOliynykGet al. Targeting a scavenger receptor on tumor-associated macrophages activates tumor cell killing by natural killer cells. Proc Natl Acad Sci. (2020) 117:32005–16. doi: 10.1073/pnas.2015343117
63
La FleurLBotlingJHeFPelicanoCZhouCHeCet al. Targeting MARCO and IL37R on immunosuppressive macrophages in lung cancer blocks regulatory T cells and supports cytotoxic lymphocyte function. Cancer Res. (2021) 81:956–67. doi: 10.1158/0008-5472.can-20-1885
64
WijarnpreechaKScribaniMRaymondPHarnoisDMKeavenyAPAhmedAet al. PNPLA3 gene polymorphism and liver- and extrahepatic cancer-related mortality in the United States. Clin Gastroenterol Hepatol. (2021) 19:1064–6. doi: 10.1016/j.cgh.2020.04.058
65
PingitorePRomeoS. The role of PNPLA3 in health and disease. Biochim Biophys Acta (BBA) - Mol Cell Biol Lipids. (2019) 1864:900–6. doi: 10.1016/j.bbalip.2018.06.018
66
HassanMMKasebAEtzelCJEl-SeragHSpitzMRChangPet al. Genetic variation in the PNPLA3 gene and hepatocellular carcinoma in USA: risk and prognosis prediction. Mol Carcinog. (2013) 52:139–47. doi: 10.1002/mc.22057
67
TanXMulayOXieJMacDonaldSKimTYeNet al. Robust and interpretable prediction of gene markers and cell types from spatial transcriptomics data. Nat Commun. (2026) 17:1781. doi: 10.1038/s41467-026-68487-0
Summary
Keywords
biomarkers, head and neck cancer, spatial transcriptomics, tumor microenvironment, tumor recurrence
Citation
Weeramange CE, Chen Y, Nguyen K-M, Zhang X, Hartel G, Causer A, Mulay O, Pavithra P, Lam PY, Gonzalez Cruz J, Vasani S, Kenny L, Taheri T, Breik O, Nguyen Q and Punyadeera C (2026) Spatial transcriptomics and histopathological analyses reveal potential biomarkers and immune features linked to recurrence in head and neck cancer. Front. Immunol. 17:1713732. doi: 10.3389/fimmu.2026.1713732
Received
26 September 2025
Revised
09 July 2026
Accepted
20 July 2026
Published
11 August 2026
Volume
17 - 2026
Edited by
Qi Zhang, Yale University, United States
Reviewed by
Yue Lyu, University of Texas MD Anderson Cancer Center, United States
Cristina Conde-Lopez, German Cancer Research Center (DKFZ), Germany
Updates
Copyright
© 2026 Weeramange, Chen, Nguyen, Zhang, Hartel, Causer, Mulay, Pavithra, Lam, Gonzalez Cruz, Vasani, Kenny, Taheri, Breik, Nguyen and Punyadeera.
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Chamindie Punyadeera, c.punyadeera@griffith.edu.au; Quan Nguyen, quan.nguyen@qimrberghofer.edu.au
†These authors have contributed equally to this work
‡ORCID: Quan Nguyen, orcid.org/0000-0001-7870-5703; Chamindie Punyadeera, orcid.org/0000-0001-9039-8259
Disclaimer
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.