Dissecting Tissue Compartment-Specific Protein Signatures in Primary and Metastatic Oropharyngeal Squamous Cell Carcinomas.

Head and neck squamous cell carcinoma (HNSCC) often presents with locoregional or distant disease, despite multimodal therapeutic approaches, which include surgical resection, chemoradiotherapy, and more recently, immunotherapy for metastatic or recurrent HNSCC. Therapies often target the primary and nodal regional HNSCC sites, and their efficacy at controlling occult distant sites remains poor. While our understanding of the tumor microenvironment conducive to effective therapies is increasing, the biology underpinning locoregional sites remains unclear. Here, we applied targeted spatial proteomic approaches to primary and lymph node metastasis from an oropharyngeal SCC (OPSCC) cohort to understand the expression of proteins within tumors, and stromal compartments of the respective sites in samples of both matched and unmatched patients. In unmatched analyses of n = 43 primary and 11 nodal metastases, our data indicated that tumor cells in nodal metastases had higher levels of Ki-67, PARP, BAD, and cleaved caspase 9, suggesting a role for increased proliferation, DNA repair, and apoptosis within these metastatic cells. Conversely, in matched analyses (n = 7), pro-apoptotic markers BIM and BAD were enriched in the stroma of primary tumors. Univariate, overall survival (OS) analysis indicated CD25 in tumor regions of primary tumors to be associated with reduced survival (HR = 3.3, p = 0.003), while progesterone receptor (PR) was associated with an improved OS (HR = 0.33, p = 0.015). This study highlights the utility of spatial proteomics for delineating the tumor and stromal compartment composition, and utility toward understanding these properties in locoregional metastasis. These findings indicate unique biological properties of lymph node metastases that may elucidate further understanding of distant metastatic in OPSCC.

Introduction

Head and neck squamous cell carcinoma (HNSCC) is the 7th leading cause of cancer worldwide, with approximately 890,000 new cases and 450,000 deaths (1). HNSCC is considered a heterogeneous malignancy arising from the upper aerodigestive tract, particularly from the squamous mucosal line. The lip, oral and nasal cavity, paranasal sinuses, larynx, nasopharynx, oropharynx, and hypopharynx are the areas involved in HNSCC (1). Several risk factors contribute to the development of HNSCC, including both tobacco and alcohol consumption. Viral-driven HNSCC is also found, with Epstein–Barr Virus (EBV) and human papillomavirus (HPV) responsible for the nasopharynx and oropharynx malignancies, respectively (1).

Oropharyngeal SCC (OPSCC) is responsible for a quarter of HNSCC (2), with a tendency to occur more commonly in non-smokers and frequent nodal involvement. Treatment of OPSCC is aimed at curing and organ preservation using a multimodality approach. In locoregional disease, patients are often treated with chemoradiation. Most notably, HPV-induced OPSCC has a better prognosis compared to HPV-negative OPSCC, namely, better radiation sensitivity and overall survival. The 5-year survival rate for patients with HPV-negative and HPV-positive OPSCC is 46 and 57.4%, respectively (3).

In approximately two-thirds of OPSCC patients, locoregional metastasis has been reported. Advanced nodal status, especially extranodal extension, is a poor prognostic predictor, and defining the molecular phenotypes of these multiple sites is important to understand not only the primary tumor, but also the influence of lymph node metastasis to develop effective therapies (4). To gain a deeper understanding of the primary and metastatic tissues, studies of the tumor microenvironment (TME) are needed, ideally by comparing primary and lymph node metastasis, and where possible, where these are from matched patients. The identification and characterization of potential TME biomarkers could have significant predictive and prognostic value in the treatment selection for OPSCC patients (5). However, this can be challenging due to inter- or intra-tumoral heterogeneity (6). The cellular composition and molecular interactions between the TME and the host immune cells could play a role in disease progression and treatment resistance. It is thought that locoregional metastasis are immune-cold and therefore treatment-resistant to any immuno-modulatory treatment strategies (1, 7).

To delineate primary and nodal disease in OPSCC, we applied digital spatial profiling (DSP) of the TME using targeted-panel multiplex immunohistochemistry of tumor- and stromal-compartments. Our study found in unmatched analyses that the tumor compartments of the primary were enriched for ARG1, PD-L1, and nodal metastases were enriched for Ki-67, PARP, BAD, and cleaved caspase 9. Stromal compartments of the primary were enriched for VISTA and IDO1. In matched analyses, BIM and BAD were enriched in the stroma of primary tumors. Furthermore, protein signatures were identified to discriminate matched primary/nodal tissues, and survival associations were investigated compartmentally by tumor sample site.

Material and Methods

This study has the approval of the Queensland University of Technology Human Research Ethics Committee (UHREC #2000000494) and University of Queensland ratification. A tissue microarray (TMA) of OPSCC specimens was sourced from the Tristar Technologies Group (USA) and contained forty-three primary tumors, eleven nodal metastases, and seven matched pairs of primary tumor/nodal metastases with concordant clinicopathological annotations.

Nanostring GeoMx Digital Spatial Profiling (DSP)

The TMA slides from OPSCC samples were obtained and analyzed using the Nanostring GeoMX Digital Spatial Profiling (DSP) technology by the Systems Biology and Data Science Group at Griffith University (Gold Coast, Australia). Pan-cytokeratin and CD45 were the visualization markers used by the instrument to stain the tumor and lymphocytes, respectively. A protein panel of 68 antibodies was used, namely, immune activation, pan-tumor, immuno-oncology (IO) drug target, cell death, immune cell typing, human immune cell core panel, and PI3K/AKT panels. The slides were prepared according to the instructions of the manufacturer, and the tumor/stroma distinction was achieved by masking on PanCK+ or PanCK− regions, respectively. Using the Nanostring nCounter® platform, antibody barcodes were counted in accordance with the instructions of the manufacturer. In DSP analysis, external RNA Controls Consortium (ERCC) QC was employed to prepare the data for further bioinformatic analysis.

Bioinformatic Analysis

Data analysis was conducted in collaboration with the Queensland Cyber Infrastructure Foundation (QCIF, QLD, Australia). The quality of data was investigated using principal component analysis, and the suitability of the RUV-III normalization method was determined using coefficients of variation (8, 9). Differential analysis was carried out using Limma packages (10). Sparse partial least squares-discriminant analysis (sPLS-DA) within the mixOmics package was used to identify multivariate minimal protein signatures (11). The Kaplan–Meier survival analysis and Cox proportional hazards models were constructed within R studio (12) using Survival package (13) and plots generated by ggplot2 (14). Data shown are not adjusted for multiple testing. The false discovery rate (FDR) adjusted results were not significant within this cohort.

Results

OPSCC Patient Cohort

To investigate protein expression in primary tumors and nodal metastases, we evaluated tissue from a tissue microarray (TMA). The TMA included unmatched specimens from forty-three primary tumors and eleven nodal metastases, and seven matched primary and nodal metastases. All primary tumors were resected from oropharyngeal regions and had squamous cell carcinoma histology ( ). The TNM staging (8th edition) ranged from T2-4, N0-2, and M0 ( ). Only one patient within the cohort was HPV-16 positive.

Table 1

OPSCC cohort characteristics.

	Nodal Metastasis, n = 11	Primary, n = 43	Matched Nodal Met/Primary, n = 7
Gender
Female	0 (0%)	7 (16%)	1 (14%)
Male	11 (100%)	36 (84%)	6 (86%)
Age
25–50	5 (45%)	5 (12%)	3 (42%)
50–95	6 (55%)	38 (88%)	4 (58%)
Status
Alive	9 (82%)	28 (65%)	2 (29%)
Deceased	2 (18%)	15 (35%)	5 (71%)
Radiation
Yes	3 (27%)	21 (49%)	1 (14%)
No	0	22 (51%)	6 (86%)
N/A	8 (73%)	N/A	N/A
Adjuvant chemotherapy
Yes	1 (9.1%)	8 (19%)	3 (43%)
No	2 (18%)	35 (81%)	4 (57%)
N/A	8 (73%)	N/A	N/A

Identification of Differentially Expressed Proteins by the Nanostring GeoMX DSP Assay

Nanostring Digital Spatial Profiler (DSP) was applied to investigate the protein expression of 68 TMA cores. By masking on PanCK+/PanCK− regions, we compartmentalized protein expression within the tumor and stromal compartments ( ). Differential expression (DE) was performed on normalized protein counts to identify compartmental enrichment of proteins in unmatched primary tumors compared to nodal metastases and between patient-matched primary tumors and nodal metastases. We discovered significant differences in key deregulated proteins between the two comparison groups. Additionally, DE within the stroma between unmatched primary tumors and nodal metastasis revealed that V-domain IG suppressor of T cell activation (VISTA) and Indoleamine 2,3-dioxygenase 1 (IDO1) exhibited higher expression in primary tumors ( ). Parallel analysis of tumor compartments uncovered DE proteins. Programmed death-ligand 1 (PD-L1), Fibroblast activation protein-alpha (FAP-a), Poly (ADP-ribrose) polymerase (PARP), Ki-67, Progesterone receptor (PR), Arginase 1 (ARG1), CD56, BCL2-antagonist of cell death (BAD), and Cleaved Caspase 9. Accordingly, PD-L1, FAP-a, PR, and ARG1 had higher expression in primary tumors, while PARP, Ki-67, CD56, BAD, and Cleaved Caspase 9 were enriched in nodal metastases ( ). Furthermore, analysis of stromal compartments across matched primary and nodal metastasis specimens unveiled seven significant differentially expressed proteins, namely, smooth muscle actin (SMA), phosphatase and tensin homolog deleted on chromosome 10 (PTEN), CD163, Bcl-2-like protein 11 (BIM), BAD, PD-L1, and CD25. BIM, BAD, and CD25 were enriched in matched primary tumors, while SMA, PTEN, CD163, and PD-L1 were higher in matched nodal metastases ( ). Analysis of tumor compartments between matched primary and nodal metastasis specimens indicated no significant DE proteins.

Figure 1

Spatial profiling was performed on tumor microarray cores from (A) primary and (B) nodal metastasis. Tissues were stained for PanCK+ (Tumor) and PanCK− (Stroma) areas. Green, PanCK; Red, CD45. Tissue segmentation strategy to capture (C) Tumor mask in purple and (D) stromal regions in green. Masks were generated per PanCK+/− feature to liberate barcodes for digital counting by nCounter.

Figure 2

Differential protein expression compared to specimens from unmatched OPSCC primary (n = 43) and nodal metastasis (n = 11). (A) Upper panel. Volcano scatter plot showing stromal enrichment of proteins in primary (left) vs nodal metastases (right) ranked by significance (−log10 P-value). Lower panel. List of top two significant deregulated proteins ranked by P-value. (B) Boxplots indicating VISTA and IDO1 enrichment in primary tumors. (C) Upper panel. Volcano scatter plot showing tumor region enrichment of proteins from primary (left) vs nodal metastases (right) ranked by significance (−log10 P-value). Lower panel. List of top nine significant deregulated proteins ranked by P-value. (D) Boxplots indicating enrichment of PD-L1 and PARP in primary tumors and nodal metastases, respectively.

Figure 3

Differential protein expression comparing specimens from matched OPSCC primary tumors and nodal metastasis (n = 7). (A) Upper panel. Volcano scatter plot showing stromal enrichment of proteins from nodal metastases (left) vs primary tumors (right) ranked by significance (−log10 P-value). Lower panel. List of top seven significant deregulated proteins ranked by P-value. (B) Representative boxplots indicating SMA and BIM enrichment in matched nodal metastasis and matched primary specimens, respectively.

Survival Associations of OPSCC Primary Tumor and Nodal Metastases

To assess the association between our protein expression and overall survival (OS), we performed a Cox proportional hazards model on all proteins. Analysis of the stromal compartment in primary tumor specimens revealed that the expression of NF1 (HR = 0.748, p = 0.025), CD27 (HR = 0.279, p = 0.035), and CD80 (HR = 0.703, p = 0.02) was associated with a better OS ( ). Moreover, tumoral compartment analysis across primary tumor samples indicated that PR (HR = 0.332, p = 0.015) was associated with better OS, while CD25 (HR = 3.311, p = 0.003) was associated with worse OS ( ). Survival analysis within the stromal compartment from nodal metastasis specimens showed that NY.ESO.1 (HR = 0.309, p = 0.016), and B7.H3 (HR = 0.249, p = 0.026) were associated with better OS ( ). Evaluation of tumoral compartment of nodal metastasis samples found several proteins associated with better OS, namely, SMA (HR = 0.685, p = 0.034), CD45 (HR = 0.463, p = 0.045), CD8 (HR = 0.251, p = 0.028), Fibronectin (HR = 0.472, p = 0.035), and STING (HR = 0.239, p = 0.024), however, BIM (HR = 1.979, p = 0.048), GZMA (HR = 4.332, p = 0.036), FOXP3 (HR = 3.258, p = 0.011), and PR (HR = 9.34, p = 0.017) were associated with worse OS ( ).

Figure 4

Identification of proteins with overall survival associations. (A, B) Forest plot indicating hazard ratio with 95% confidence interval for proteins from primary specimens. (C, D) Forest plot indicating hazard ratio with 95% confidence interval for proteins from nodal metastases. HR >1 demonstrates association with poorer outcome.

Multivariate Discrimination of OPSCC Primary Tumors From Nodal Metastases

Multivariate analysis by sparse partial least-squares discriminant analysis (sPLSDA) was employed to identify minimal protein signatures that collectively distinguish primary tumors from nodal metastases. Signatures within the stroma of matched primary vs nodal metastases stratified samples effectively ( ). The first signature ( ) included levels of SMA, PTEN, cleaved caspase 9, and CD25 (AUC = 0.979) ( ). A second signature comprised of CD95, CD80, and CD27 distinguished matched sample types as well (AUC = 0.918) ( ).

Figure 5

A multi-protein signature differentiates OPSCC tumor progression. (A) sPLSDA distinguishes the groups (Matched Primary vs Matched Nodal Mets) by protein signatures in stroma. (B) Features of discriminating proteins per component in stroma. (C, D) ROC curve of each signature was used to differentiate the groups (Matched Primary vs Matched Nodal Mets). Color of component loadings indicates patient group in which feature was maximally expressed. Positive or negative values in bar chart indicate positive or negative loading to the discriminant signature.

Discussion

HNSCC has a high risk of locoregional nodal metastasis, which affects patient prognosis and treatment outcomes (15). Patients with nodal metastasis are considered to have locoregionally advanced disease with a lower chance of remission (16). Currently, the predictors of nodal metastasis include tumor thickness and size, which have been shown to be unreliable predictors (17). It has been difficult to manage clinically negative neck nodes (N0) due to a lack of reliable predictors of occult metastasis (16, 18). In a study conducted by Shah et al., the authors found that there was a 40% chance of nodal metastasis in clinically node-negative neck dissections (19). For the purpose of distinguishing patients with a high risk of nodal metastasis, various pathological and clinical factors, namely, lymphovascular invasion, tumor differentiation, depth of invasion (DOI), and pattern of invasion (POI), have been reported (20, 21). Moreover, modern imaging modalities such as MRI, CT imaging, and PET/CT scanning, have been used to aid in the detection of locoregional nodal metastases. However, some radiographic features, namely, “subclinical”, or “microscopic”, or “occult” disease, remain difficult to diagnose using any of these approaches (22, 23). Therefore, companion diagnostics tools are needed to improve the prediction of the likelihood of the development of locoregional metastasis. Ideally, this would be possible by interrogating the primary tissue to determine its aggressiveness and propensity for metastasis.

Tumor tissue analysis by bulk expression or single cell RNA sequencing offers an overview of the molecular features of HNSCC tumors and their TME. These methods are incapable of revealing the spatial cellular properties required for the anti-tumor immune responses (24). Spatial proteomic approaches can provide compartment-specific tumor information to aid in delineating tumor composition. To garner insight into these properties that distinguish primary OPSCC tumors from their metastatic nodal counterparts, we have employed Digital Spatial Profiling to address a targeted profile of proteins and present this data as a first step in profiling HNSCC nodal involvement in an OPSCC cohort.

Unmatched analyses between primary tumors (n = 43) and nodal metastases (n = 11) provided insight into the potential dysregulation of several proteins, despite the inherent limitations associated with such a sampling strategy. VISTA, IDO1, and PD-L1 are key immune checkpoints, with VISTA and IDO1 appearing more abundant in the stromal compartment of primary tumors, whereas PD-L1 indicated higher expression in their respective tumor compartments. The V-domain Ig suppressor of T cell activation (VISTA) is an inhibitory immune checkpoint protein that is typically expressed on naïve CD4+ and Foxp3+ Tregs and functions by inhibiting T-cell proliferation and promoting naive to Treg conversion (25). VISTA was associated with several immune cell regions in the stroma but not in HNSCC tumors (26). Blockade of VISTA was found to boost anti-tumor immunity in the tumor microenvironment by increasing the number of activated dendritic cells (DCs) and decreasing the number of myeloid-derived suppressor cells (MDSCs). IDO1 induces T-cell apoptosis through activation of caspase 8 and releases mitochondrial cytochrome C, functioning in an immunosuppressive capacity (27). Programmed death-ligand 1 (PD-L1) is a canonical inhibitory immune checkpoint that binds PD-1 on the surface of tumors and immune cells (28). PD-L1 expression on the surface of HNSCC tumor cells is associated with a more robust anti-tumor immune response (29, 30).

Apoptotic pathways play an important role in tumorigenesis, and our results indicate that pro-apoptotic BAD and cleaved caspase 9 were enriched in tumor cells of unmatched nodal metastases. Conversely, in the matched analysis, BIM and BAD were enriched in the stroma of primary tumors. Under various physiological and patho-physiological conditions, Bcl-2 interacting mediator of cell death (BIM) promotes the intrinsic apoptotic pathway (31). Bcl2-associated agonist of cell death (BAD) is a member of the BCL2 family of proteins that act as pro-apoptotic regulators (32). The expression of BAD has been linked to chemoresistance in cancer patients (33, 34). Caspase 9 functions as a pro-apoptotic regulator, allowing the activation of effector caspases 3 and 7 (35). Caspase 9-induced apoptosis has been linked to chemotherapy response. Studies have shown that HNSCC tumors may be resistant to cisplatin if they have a reduced expression of caspase 9 (36).

It is interesting to note that despite an imbalance in samples in the unmatched analysis, several proteins appear enriched within nodal metastases relative to primary tumors. In addition to the pro-apoptotic markers above, Ki-67 and PARP appear enriched in nodal tumor cells. Ki-67 is an established proliferation marker (37), while PARP responds to DNA damage by recruiting effector proteins to repair single-strand breaks (38, 39). PARP inhibitors have been studied as a promising drug to overcome the limitations of conventional therapies that cause DNA damage, such as chemotherapy or radiotherapy (40). This pro-apoptotic, proliferative, and DNA damage phenotype of unmatched nodal metastatic cells is a novel finding in our data that requires further validation, perhaps indicating tumor evolution or response to changes in the cellular ecosystem of the lymph node.

Several other notable features of our data include increased expression of CD25 in the stroma of matched primary tumors. CD25, also known as the IL-2 receptor alpha, is a protein found on activated T cells, specifically Tregs (41). Interestingly, we found that it was CD25 expression within primary tumor regions, not stroma, that was associated with poorer OS. Additionally, PR expression appeared associated with better survival within primary tumors and was also enriched within their tumor regions relative to nodal metastases.

In addition to the differential expression of each individual protein, we applied a multivariate statistical model (sPLSDA) to further discern features that collectively discriminated between our matched patient samples. Of note, this model only performed effectively in stratifying these sample types by their stroma. Expression of CD25 and cleaved caspase 9 in the primary samples and PTEN and SMA in nodal samples could discern samples. Similarly, levels of CD80 and CD27 within primary samples and CD95 in nodal samples could separate these samples. This method offers an alternative to traditional differential expression that may provide insight into contributing differences in observed phenotypes using a multivariate approach.

A protein association with overall survival (OS) was investigated using the Univariate Cox proportional hazards model. Stromal expression of CD27 in primary tumor specimens, and NY-ESO.1 and B7-H3 expression in the nodal metastasis samples, was associated with improved OS. In nodal metastases tumor regions, CD8 and STING were associated with improved OS; however, GZMA was associated with poorer OS. Interestingly, the expression of PR in the tumoral compartment of primary and nodal metastasis tumors demonstrated a different survival pattern. Although the PR expression was associated with improved OS in primary tumors, it was associated with worse OS in nodal metastasis specimens. Immune response protein markers, CD27, and NY.ESO.1 within the stromal compartment of primary and nodal metastasis specimens, were associated with improved OS in our study. Cluster of differentiation 27 (CD27) belongs to the tumor necrosis factor (TNF) receptor superfamily and is involved in T and B cell co-stimulation (42). New York esophageal squamous cell carcinoma 1 (NY.ESO.1) is a member of the cancer testis antigen (CTA) family, which regulates both humoral and cellular immune responses. NY-ESO.1 expression has been linked to higher tumor differentiation grade and stage, and lymph node metastasis (43). In our study, we found that B7.H3 and STING protein expression, NF-κB pathway markers, were linked to improved OS. B7.H3, also known as CD276, promotes anti-tumor immune response by activating T and NK cells (44). STING, on the other hand, contributes to the immune response to tumor cells through the upregulation of interferon gamma 1 (IFN1) (45). Progesterone receptor (PR) is a type of androgen receptor and a member of the nuclear receptor family of transcription factors that regulates target gene expression networks in response to its ligand progesterone (46).

Our study has identified tumor and stromal compartment-specific proteins and signatures that may have predictive and prognostic implications for HNSCC and the development of nodal metastasis. Nevertheless, the study is impacted by the number of samples for each cohort, in particular the matched group. We propose further investigation to profile primary, locoregional, and distant metastasis from matched patient samples to understand the molecular features driving the development of metastasis in OPSCC.

Data Availability Statement

The data presented in the study are deposited in the Geo repository, under accession number GEO Submission GSE200601.

Ethics Statement

This study has the approval from the Queensland University of Technology Human Research Ethics Committee (UHREC #2000000494). The patients/participants provided their written informed consent to participate in this study.

Author Contributions

Experimental design: HR, BH, and AK. Methodology: HR, JM, AM, BH, and AK. Data analysis and reporting: HR, JM, AM, KO, BGMH, and AK. Manuscript preparation and critical review: all authors. All authors listed have made a substantial, direct, and intellectual contribution to the work and approved it for publication.

Funding

This project is funded by the Garnett Passe and Rodney Williams Memorial Foundation (GPRWMF) conjoint grant for AK and BH. AK is supported by a fellowship from the NHMRC (1157741). KO’B and AK are supported by the PA Research Foundation.

Conflict of Interest

AM is employed by QCIF Bioinformatics.

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

Publisher’s Note

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