Bioinformatics analysis of the role of CXC ligands in the microenvironment of head and neck tumor

Fengyang Jing; Jianxiong Wang; Liming Zhou; Yujie Ning; Shengqian Xu; Youming Zhu

doi:10.18632/aging.203269

Research Paper Volume 13, Issue 13 pp 17789—17817

Bioinformatics analysis of the role of CXC ligands in the microenvironment of head and neck tumor

Fengyang Jing^{1,
,} , Jianxiong Wang^{2,
,} , Liming Zhou^1, , Yujie Ning^1, , Shengqian Xu^2, , Youming Zhu^1, ,

¹ Department of Dental Implant Center, Stomatologic Hospital and College, Anhui Medical University, Key Laboratory of Oral Diseases Research of Anhui Province, Hefei 230032, China
² Chief Physician, Department of Rheumatology and Immunology, The First Affiliated Hospital of Anhui Medical University, Hefei 230022, China

* Equal contribution

Received: February 10, 2021 Accepted: May 18, 2021 Published: July 11, 2021

https://doi.org/10.18632/aging.203269
How to Cite

Copyright: © 2021 Jing et al. This is an open access article distributed under the terms of the Creative Commons Attribution License (CC BY 3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Abstract

Chemokines play a significant role in cancer. CXC-motif chemokine ligands (CXCLs) are associated with the tumorigenesis and progression of head and neck squamous cell carcinoma (HNSC); however, their specific functions in the tumor microenvironment remain unclear. Here, we analyzed the molecular networks and transcriptional data of HNSC patients from the Oncomine, GEPIA, String, cBioPortal, Metascape, TISCH, and TIMER databases. To verify immune functions of CXCLs, their expression was analyzed in different immune cell types. To our knowledge, this is the first report on the correlation between CXCL9–12 and 14 expression and advanced tumor stage. CXCL2, 3, 8, 10, 13, and 16 were remarkably related to tumor immunity. Kaplan–Meier and TIMER survival analyses revealed that high expression of CXCL1, 2, 4, and 6–8 is correlated with low survival in HNSC patients, whereas high expression of CXCL9, 10, 13, 14, and 17 predicts high survival. Only CXCL13 and 14 were associated with overall survival in human papilloma virus (HPV)-negative patients. Single-cell datasets confirmed that CXCLs are associated with HNSC-related immune cells. Thus, CXCL1–6, 8–10, 12–14, and 17 could be prognostic targets for HNSC, and CXCL13 and 14 could be novel biomarkers of HPV-negative HNSC.

Introduction

Head and neck squamous cell carcinoma (HNSC) constitutes 5.7% of cancer-related mortality worldwide [1]. Lymph node metastasis is known to reduce the cancer-specific survival of HNSC patients by approximately 50% [2], while extranodal extension (ENE) increases regional and distant metastatic failure [3]. However, the pathogenicity of local tumor cell implantation during primary surgical intervention and the effects of wound healing growth factors on HNSC proliferation remain unclear [4]. Although human papilloma virus (HPV) infection is currently the only molecular biomarker available to individualize treatment options for oropharyngeal cancer, recent studies have elucidated biological interactions within the tumor microenvironment [5, 6].

Chemokines are a small class of cytokine-like molecules with four subtypes: C-X-C (CXC), C-C (CC), C-X3-C (CX3C), and X-C (XC), where X represents a conserved terminal cysteine residue [7]. Various CXC-motif chemokine ligands (CXCLs), termed CXCL1–17, have been identified to date [8]; all of these, except CXCL15, are found in humans (Supplementary Table 1). Chemokines indirectly participate in tumor development by affecting angiogenesis and the interaction between tumors and leukocytes, as well as by directly affecting tumor transformation, survival, growth, invasion, and metastasis. Although chemokines can kill tumors by regulating white blood cell infiltration and activating an immune response [9, 10], they can also promote tumor tissue angiogenesis, accelerate tumor cell proliferation, and promote basement membrane invasion, thereby promoting tumor growth and metastasis [11–13].

Several studies have shown that CXCLs may individually or synergistically play complex and distinct roles in HNSC by affecting tumor cell proliferation, migration, invasion, growth, and survival via various pathways to modulate tumor growth and metastasis [14–28]. Although these studies have elucidated the roles of different CXCLs in various oral cancers, the potential activation or inhibition mechanisms and functions of CXCLs in the HNSC tumor microenvironment have not yet been fully elucidated.

DNA and RNA sequencing techniques are important aspects of biological and biomedical research that have been revolutionized by the development of microarray technology; however, the roles of CXCLs in HNSC have not yet been analyzed using bioinformatic methods, and the immune functions of CXCLs in HNSC remain unclear. Based on published gene expression and copy number variation datasets, we analyzed the expression of CXCLs in HNSC patients to determine the potential functions and prognostic value of CXCLs. In addition, we verified the immune-related functions of CXCLs in HNSC using single-cell datasets and analyzed the effect of CXCL expression on the survival of HPV-positive and -negative patients.

Results

CXCL transcript expression in HNSC patients

First, we used Oncomine to compare the expression of the 16 CXCLs identified in human cells between HNSC and normal samples (Figure 1). Although the mRNA expression of CXCL12 and 17 was decreased in two datasets, that of all other CXCLs was significantly upregulated in HNSC patients from seven datasets [29–35] (Supplementary Table 2). Therefore, we compared the mRNA expression of CXCLs between HNSC and normal tissues using GEPIA. The expression of CXCL1, 8–11, and 13 was higher, while that of CXCL12 and 17 was lower, in HNSC than in normal tissues (Figure 2).

CXC chemokine mRNA expression in different types of cancer. Number of datasets with statistically significant CXC chemokine mRNA expression: red, upregulated; blue, downregulated. Fold changes and p values are shown in Supplementary Table 2.

Figure 1. CXC chemokine mRNA expression in different types of cancer. Number of datasets with statistically significant CXC chemokine mRNA expression: red, upregulated; blue, downregulated. Fold changes and p values are shown in Supplementary Table 2.

CXC-motif chemokine ligand (CXCL) expression in head and neck squamous cell carcinoma (HNSC). (A) Differential scatter diagram of CXCL expression in HNSC: red, upregulated; green, downregulated. (B) Differential box plot of CXCL expression in HNSC. *p

Figure 2. CXC-motif chemokine ligand (CXCL) expression in head and neck squamous cell carcinoma (HNSC). (A) Differential scatter diagram of CXCL expression in HNSC: red, upregulated; green, downregulated. (B) Differential box plot of CXCL expression in HNSC. *p < 0.05.

To verify the differences in CXCL expression between HNSC tumor and normal tissues as well as between HPV-positive and -negative HNSC patients, we used the Tumor IMmune Estimation Resource (TIMER) database. Consistently, various CXCLs were expressed differently among these samples, with CXCL5, 7–10, 13, 14, 16, and 17 displaying clear differences between HPV-positive and -negative patients (Supplementary Figures 1, 2). In addition, only CXCL9–12 and 14 expression varied significantly with HNSC tumor stage (Supplementary Figure 3). Together, these differences in the expression of members of the CXCL family in HNSC indicate that CXCLs play different but important roles in the HNSC tumor microenvironment.

CXCL gene expression and mutation in HNSC patients

To analyze the changes in CXCL expression and their correlations with HNSC, we used the cBioPortal online tool. CXCL expression was altered in samples collected from 488 HNSC patients, whereas the queried genes were altered in 218 (45%) of the 488 queried patients (Figure 3A). We also used cBioPortal to analyze CXCL mRNA expression (log RNA sequencing [RNA-Seq] version (v.)2 RSEM) and the correlations among CXCLs using Pearson’s correction. The following significant positive correlations were identified: CXCL1 correlated with CXCL2, 3, 6, and 8; CXCL2 with CXCL1, 3, and 8; CXCL3 with CXCL1, 2, 4, 5, and 8; CXCL4 with CXCL3 and 8; CXCL5 with CXCL3, 7, and 8; CXCL6 with CXCL1; CXCL7 with CXCL5; CXCL8 with CXCL1, 2, 3, 4, and 5; CXCL9 with CXCL10, 11, and 13; CXCL10 with CXCL9 and 11; CXCL11 with CXCL9 and 10; CXCL12 with CXCL13; and CXCL13 with CXCL9 and 12 (Figure 3B). These relationships were verified using the TIMER database, with similar outcomes (Supplementary Figure 4), suggesting that the functions of multiple members of the CXCL family in HNSC are related.

CXC-motif chemokine ligand (CXCL) gene expression and mutation analysis in head and neck squamous cell carcinoma (HNSC). (A) CXCL gene expression and mutations in HNSC were analyzed using cBioPortal. (B) The associations among CXCLs in HNSC were analyzed using the cBioPortal database. Darker colors indicate a stronger correlation.

Figure 3. CXC-motif chemokine ligand (CXCL) gene expression and mutation analysis in head and neck squamous cell carcinoma (HNSC). (A) CXCL gene expression and mutations in HNSC were analyzed using cBioPortal. (B) The associations among CXCLs in HNSC were analyzed using the cBioPortal database. Darker colors indicate a stronger correlation.

Predicted functions and pathways of CXCL-related genes in HNSC patients

The top 50 CXCL-related genes were detected using GEPIA2 (Supplementary Table 3) and then used to produce a network with CXCLs using String to visualize the relationships among them (Figure 4). Another network was produced using Metascape to visualize the function of these genes (Supplementary Figure 5), as follows: extracellular structure organization, blood vessel morphogenesis, cell–substrate adhesion, endothelium development, platelet alpha granule, calcium ion binding, calcium-dependent cell–cell adhesion via plasma membrane cell adhesion molecules, cytokine binding, protein homodimerization activity, protein localization to cell surface, regulation of bone mineralization, extracellular matrix glycoproteins, peptidyl-tyrosine phosphorylation, collagen trimer, cellular response to transforming growth factor beta stimulus, response to wounding, G protein-coupled receptor binding, cell surface interactions at the vascular wall, Rho GTPase cycle, and camera-type eye development.

Protein–protein interaction network. The relationships between CXCLs and the top 50 similar genes were visualized using the STRING database (see Supplementary Table 3 for a detailed gene list).

Figure 4. Protein–protein interaction network. The relationships between CXCLs and the top 50 similar genes were visualized using the STRING database (see Supplementary Table 3 for a detailed gene list).

The functions of the CXCLs were predicted using Gene Ontology (GO) analysis and Kyoto Encyclopedia of Genes and Genomes analysis (Figure 5) and visualized using the R package clusterProfiler [36]. In order to get more accurate results, we next performed Gene Set Enrichment Analysis (GSEA). We used the RNAseq datatype of TCGA-HNSC in the Linkedomics database. We set the minimum number of genes as 3 and the simulations were 500. Results were visualized using ggplot2 Rpackage (Figure 6). Taken together, the results of these functional enrichment analyses suggest that CXCLs play an important role in the HNSC tumor microenvironment via immune-related signaling pathways, which have been unreported previously.

Enrichment analysis of differentially expressed CXC-motif chemokine ligands (CXCLs) in head and neck squamous cell carcinoma (HNSC). GO enrichment analysis predicted the functional roles of target genes based on three aspects, including (A) biological processes (BP), (B) molecular functions (MF), and (C) cellular components (CC). (D) Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis.

Figure 5. Enrichment analysis of differentially expressed CXC-motif chemokine ligands (CXCLs) in head and neck squamous cell carcinoma (HNSC). GO enrichment analysis predicted the functional roles of target genes based on three aspects, including (A) biological processes (BP), (B) molecular functions (MF), and (C) cellular components (CC). (D) Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis.

Gene set enrichment analysis (GSEA) analyzed the Kyoto encyclopedia of genes and genomes (KEGG) enrichment of CXC-motif chemokine ligands (CXCLs) in head and neck squamous cell carcinoma (HNSC).

Figure 6. Gene set enrichment analysis (GSEA) analyzed the Kyoto encyclopedia of genes and genomes (KEGG) enrichment of CXC-motif chemokine ligands (CXCLs) in head and neck squamous cell carcinoma (HNSC).

CXCL expression in different immune cells in HNSC

Since the enriched pathways were related to immune function, we analyzed CXCL expression in two datasets from the Tumor Immune Single-cell Hub (TISCH) database: HNSC_GSE103322 (5,902 cells) and HNSC_GSE139324 (130,721 cells) (Figure 7). Detailed information about the distribution and proportion of immune cells in each dataset is shown in Figure 8. We found that CXCL2, 3, 8, 10, 12–14, and 16 expression differed in immune cells of the HNSC_GSE103322 dataset (Supplementary Table 4), whereas CXCL1–3, 8–10, 13, and 16 expression differed in the immune cells of HNSC_GSE139324 (Supplementary Table 5). This differential expression of CXCL1–3, 8–10, 12–14, and 16 in different immune cell types in HNSC suggests that they play an important role in the tumor microenvironment by affecting immune cells.

Figure 7. CXC-motif chemokine ligand (CXCL) expression in two single-cell datasets. CXCL expression in different immune cells in (A) HNSC_GSE103322 and (B) HNSC_GSE139324.

Figure 8. Infiltration site distribution (left) and proportion (right) of different cell types in HNSC. (A) HNSC_GSE103322. (B) HNSC_GSE139324.

CXCLs and immune cell infiltration in HNSC patients

Since TISCH analysis revealed that CXCLs are involved in inflammatory responses and immune cell infiltration in HNSC patients, we comprehensively explored the correlations between differentially expressed CXCLs and immune cell infiltration using the TIMER database (Supplementary Figure 6). Notably, CXCL2, 3, 5, 7–14, 16, and 17 were related to various types of immune cells (Supplementary Table 6). Therefore, these findings suggest that CXCL2, 3, 8, 10, 13, and 16 greatly affect immune functions and thereby influence HNSC tumorigenesis and tumor development.

Association between CXCL mRNA expression and HNSC prognosis

Next, we used the data from The Cancer Genome Atlas (TCGA) to obtain receiver operating characteristic curve (ROC) to evaluate the diagnostic efficacy of CXCLs in the development of HNSC (Figure 9). The area under the curve value of CXCL10, 11, and 13 was 0.8 fold higher than that of other CXCLs, indicating that the expression of these genes could distinguish tumors from non-tumors and used as a potential diagnostic biomarker.

Diagnostic efficacy of CXC-motif chemokine ligands (CXCLs) in head and neck squamous cell carcinoma (HNSC). (A) Receiver operating characteristic curve (ROC) of CXCL1, 2, 3, 4. (B) ROC of CXCL5, 6, 7, 8. (C) ROC of CXCL9, 10, 11, 12. (D) ROC of CXCL13, 14, 16, 17. TPR: True Positive Rate, FPR: False Positive Rate.

Figure 9. Diagnostic efficacy of CXC-motif chemokine ligands (CXCLs) in head and neck squamous cell carcinoma (HNSC). (A) Receiver operating characteristic curve (ROC) of CXCL1, 2, 3, 4. (B) ROC of CXCL5, 6, 7, 8. (C) ROC of CXCL9, 10, 11, 12. (D) ROC of CXCL13, 14, 16, 17. TPR: True Positive Rate, FPR: False Positive Rate.

Finally, we explored the effect of CXCL mRNA expression on the survival of HNSC patients using Kaplan–Meier Plotter curve and log rank test analyses. Decreased CXCL1, 2, 4, and 6–8 mRNA levels and increased CXCL9–14 and 17 mRNA levels were significantly associated with the overall survival (OS) of all HNSC patients (Figure 10A). In addition, patients with high CXCL5, 7, 14, or 17 mRNA levels or low CXCL1–4, 8, 11, or 12 mRNA levels were predicted to have high recurrence-free survival (RFS; Figure 10B). Therefore, we used TIMER to analyze the relationship between immune cell number and cumulative survival (Supplementary Figure 7), as well as the relationship between CXCL expression and the cumulative survival of HPV-positive and -negative HNSC patients. Decreased CXCL1, 2, 4–8, and 14 expression and increased CXCL9 and 13 expression were associated with high survival in HPV-positive HNSC patients; however, only high expression levels of CXCL13 and 14 were related to improved prognosis in HPV-negative HNSC patients (Figure 11). Together, our analysis of the relationship between CXCL expression and the prognosis of HNSC patients in large scale data from different databases verified that CXCL family members could exert significantly different effects on the long-term survival of HNSC patients.

Prognostic value of CXC-motif chemokine ligands (CXCLs) in head and neck squamous cell carcinoma (HNSC). (A) Prognostic value of CXCLs for the overall survival of HNSC patients. (B) Prognostic value of CXCLs for the recurrence-free survival of HNSC patients.

Figure 10. Prognostic value of CXC-motif chemokine ligands (CXCLs) in head and neck squamous cell carcinoma (HNSC). (A) Prognostic value of CXCLs for the overall survival of HNSC patients. (B) Prognostic value of CXCLs for the recurrence-free survival of HNSC patients.

Prognostic value of CXCLs for the overall survival of HPV-positive and -negative head and neck squamous cell carcinoma (HNSC) patients. Upper: HNSC-HPVpos, HPV-positive HNSC patients. Lower: HNSC-HPVneg, HPV-negative HNSC patients.

Figure 11. Prognostic value of CXCLs for the overall survival of HPV-positive and -negative head and neck squamous cell carcinoma (HNSC) patients. Upper: HNSC-HPVpos, HPV-positive HNSC patients. Lower: HNSC-HPVneg, HPV-negative HNSC patients.

Discussion

CXCLs play key roles in many cancers and are associated with diverse regulatory functions in HNSC [37]. Although some studies have reported the roles of various CXCLs in HNSC [38–40], no systematic analyses have been performed on this gene family. Moreover, few studies have investigated the relationship between CXCLs and the tumor microenvironment in HNSC. To our knowledge, this study is the first to analyze single-cell databases to explore the immune roles of CXCLs in the HNSC tumor microenvironment and identify key CXCLs that may play important roles in tumorigenesis and tumor progression.

Initially, we analyzed the relationship between CXCLs and HNSC stage and found that CXCL9–12 and 14 were differentially expressed in different tumor stages. Notably, CXCL9–11 and 14 were highly expressed in HNSC in the Oncomine, GEPIA, and TIMER databases and were associated with good prognosis. Therefore, changes in the expression of these molecules are likely associated with tumor progression and should be detected during the diagnosis, treatment, and prognostic evaluation of HNSC patients. Consistently, Yang et al. [39] reported that CXCL1 correlated negatively with the 5-year OS of oral squamous cell carcinoma patients (p < 0.05), whereas CXCL10 and 9 correlated positively with their 5-year OS (p < 0.05). Moreover, we found that high CXCL8 expression was associated with poor prognosis in HNSC, whereas high CXCL13, 14, and 17 expression may improve the OS of these patients; however, we found no difference between the expression of these molecules and RFS. Low CXCL2, 3, and 12 mRNA levels predicted high RFS, suggesting that they play important regulatory roles in promoting tumor recurrence. Although few studies have investigated CXCL14, high CXCL14 expression has been found to suppress tumor growth in HNSC [41]. We found that CXCL14 is associated with tumor stage, OS, and RFS, suggesting that it could be an important target for treating HNSC.

Previously, Szabo et al. [42] demonstrated that the upregulation of CXCL1 and 8 is not cancer specific, supporting the hypothesis that similar mechanisms exist in wound healing and oncogenesis. Here, we found a strong positive correlation between CXCL1 and 8 as well as many other CXCLs. Although most studies have focused on the roles of individual CXCLs in tumors, our results suggest that multiple molecules exert a combined effect on the tumor microenvironment. However, we found that CXCL4 is not associated with HNSC but is closely related to periodontitis and significantly and positively correlated with CXCL3 and 8. An in-depth understanding of the tumor microenvironment and tumorigenesis mechanism would clarify whether CXCL4, along with CXCL3 and 8, exerts potential synergistic effects on HNSC tumorigenesis.

In this study, our systematic enrichment analysis of CXCLs identified several signaling pathways that have not been previously reported to be involved in the function of CXCLs in HNSC. These findings suggest that HNSC progression is closely related to CXCL-mediated inflammation. Further analysis of single-cell datasets showed that the expression of many CXCLs was closely related to the role of immune cells in HNSC, confirming our hypothesis; however, CXCL16 was not closely associated with HNSC. Although no studies have yet reported the correlation between CXCL16 and HNSC, we observed a significant difference in CXCL16 expression in classically (M1) and alternatively (M2) activated macrophages in HNSC, suggesting that CXCL16 has immunomodulatory roles in HNSC. Further studies with more clinical samples are required to elucidate the specific molecules and mechanisms involved.

In conclusion, this study provides a thorough understanding of the heterogeneity and complexity of the molecular biological properties of HNSC by analyzing the expression and prognostic value of CXCLs in HNSC. In addition, we verified the function of CXCLs as biomarkers in HNSC and found, for the first time, that CXCL16 may play a significant role in the occurrence and development of HNSC by modulating immunity. Furthermore, we revealed that CXCL13 and 14 may be exclusive biomarkers for HPV-negative HNSC and that the top 50 CXCL-related genes enriched for signaling pathways may be closely related to HNSC development. Because of the lack of detailed single-cell data for HNSC, we were unable to analyze the role of CXCLs in immune infiltration in greater detail; therefore, we aim to improve these data in the future and conduct in-depth analyses to identify more significant and meaningful biomarkers, particularly for HPV-negative HNSC. Taken together, we believe that our findings improve the current knowledge of HNSC and could improve the diagnostic accuracy, treatment, and prognosis of HNSC patients.

Materials and Methods

Oncomine analysis

Oncomine is an online cancer gene expression microarray database that can be used to analyze CXCL transcription levels in different cancers (https://www.oncomine.org/resource/login.html) [43]. In this study, we compared CXCL mRNA expression in clinical cancer tissues with that in normal control tissues, with p value and fold change cutoffs of 0.05 and 2, respectively.

GEPIA dataset

GEPIA (http://gepia.cancer-pku.cn/) is a newly developed interactive web server that can analyze the RNA sequencing expression data of 9,736 tumors and 8,587 normal samples from TCGA and the Genotype Tissue Expression (GTEx) project using a standard processing pipeline [44]. Using this database, we verified differences in CXCL expression between tumor and normal tissues and identified correlations among CXCLs. The top 50 CXCL-related genes were detected using GEPIA2, a python package that provides rapid analysis and result retrieval.

cBioPortal

The cBio Cancer Genomics Portal (cBioPortal, http://www.cbioportal.org/) provides data for more than 5,000 tumor samples from 20 cancer studies [45, 46]. The genomic profiles (523 samples) include mutations, putative copy number alterations from the genomic identification of significant targets in cancer, and mRNA expression z scores relative to all samples (log RNA Seq V2 RSEM) with a threshold of ±2.0.

String

To visualize the relationships among the top 50 CXCL-related genes, we used the String database (https://string-db.org/), which uses a spring model to generate network images in which nodes are modeled as masses and edges are modeled as springs. Final node position is calculated by minimizing the “energy” of the system. High confidence edges are given a higher “spring strength” so that they reach the best position before lower confidence edges. By default, the database sets high confidence edge length to 80% of the normal length [47–57].

Metascape

Metascape (http://metascape.org/) is a web-based portal that allows experimental biologists to analyze and annotate using a comprehensive gene list. In particular, Metascape combines functional enrichment, interactive group analysis, gene annotation, and member search by combining more than 40 independent knowledge bases in one integrated portal. In addition, Metascape can conveniently compare and analyze datasets from several independent and orthogonal experiments [58]. In this study, we used Metascape to visualize the GO enrichment results.

TISCH

The TISCH (http://tisch.comp-genomics.org/) is a scRNA-seq database that provides detailed cell-type annotation at the single-cell level to explore the tumor microenvironment across different cancer types. TISCH integrates the single-cell transcriptomic profiles of nearly 2 million cells from 76 high-quality tumor datasets across 27 cancer types [59]. We analyzed CXCL expression in various immune cells at the single-cell level using two HNSC datasets in TISCH.

Kaplan–Meier plotter

The Kaplan–Meier plotter can assess the effect of 54,000 genes (mRNA, miRNA, and protein) on survival in 21 cancer types using source databases such as the intergovernmental Group on Earth Observations, European Genome-phenome Archive, and TCGA [60]. We analyzed OS and RFS using “autoselect best cutoff,” wherein all possible cutoff values between the lower and upper quartiles were computed, and the best performing threshold was used as a cutoff. Data from patients that survived over the selected follow-up threshold were not used to generate the plot.

TIMER

The TIMER database is a network resource that systematically evaluates the clinical effects of different immune cells in different types of cancer (http://cistrome.dfci.harvard.edu/TIMER/). A new statistical method produced by the developers was used to estimate the abundance of six immune cell types in the tumor microenvironment: B cells, CD4 T cells, CD8 T cells, neutrophils, macrophages, and dendritic cells [61–63]. We used TIMER to analyze the correlation among CXCLs and between CXCL expression and the infiltration of various immune cell types. In addition, we analyzed the difference in CXCL expression between tumor and normal tissues as well as between HPV-positive and -negative tumor tissues and analyzed the effect of CXCLs on OS in HPV-positive and -negative HNSC.

Data availability

The differences in gene expression were analyzed using Oncomine, GEPIA, TIMER, and TCGA (https://portal.gdc.cancer.gov/), while the relationships among CXCLs were identified using cBioPortal. The top 50 most similar genes were obtained using GEPIA, and the resulting network was assembled and analyzed using String, Metascape, LinkedOmics [64], and DAVID. We confirmed our results using two single-cell datasets (HNSC_GSE103322 and HNSC_GSE139324) from TISCH and analyzed immune infiltration using TIMER. Finally, the survival curves were analyzed using the Kaplan–Meier Plotter and TIMER.

Statistical analysis

All statistical analyses were carried out using SPSS Statistics 20.0 software (IBM SPSS, Chicago, USA). The statistical differences between the experimental groups were analyzed by Student’s t test. The relationship between CXCL expression and survival rates was analyzed using Kaplan–Meier (Km) curve. The p values were obtained using Student’s t test.

Supplementary Materials

Supplementary Figures

Supplementary Tables

Author Contributions

F.J. and J.W. designed the project. F.J., J.W., L.Z. and Y.N. performed the study, analyzed data and aided in writing the manuscript. S.X. and Y.Z. edited the manuscript. All authors read and approved the final manuscript.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Funding

This work was funded by the National Natural Science Foundation of China (31970677, 31501103 and 81571572), and the Research Level Improvement Plan of Anhui Medical University (2019xkjT010).

References

1. Patterson RH, Fischman VG, Wasserman I, Siu J, Shrime MG, Fagan JJ, Koch W, Alkire BC. Global Burden of Head and Neck Cancer: Economic Consequences, Health, and the Role of Surgery. Otolaryngol Head Neck Surg. 2020; 162:296–303. https://doi.org/10.1177/0194599819897265 [PubMed]
2. Chinn SB, Myers JN. Oral Cavity Carcinoma: Current Management, Controversies, and Future Directions. J Clin Oncol. 2015; 33:3269–76. https://doi.org/10.1200/JCO.2015.61.2929 [PubMed]
3. Reddy JP, Pettaway CA, Levy LB, Pagliaro LC, Tamboli P, Rao P, Jayaratna I, Hoffman KE. Factors associated with regional recurrence after lymph node dissection for penile squamous cell carcinoma. BJU Int. 2017; 119:591–97. https://doi.org/10.1111/bju.13686 [PubMed]
4. Lindgren G, Kjellén E, Wennerberg J, Ekblad L. Wound-healing factors can prime head and neck cancer cells to increase their tumor-forming capacity. Laryngoscope. 2016; 126:E213–17. https://doi.org/10.1002/lary.25907 [PubMed]
5. Leemans CR, Snijders PJ, Brakenhoff RH. The molecular landscape of head and neck cancer. Nat Rev Cancer. 2018; 18:269–82. https://doi.org/10.1038/nrc.2018.11 [PubMed]
6. Gualtieri T, Bonomo P, Sottili M, Mangoni M, Lavarone A, Russo ML, Desideri I, Livi L, Deganello A. Biomarkers in wound drainage fluids of head and neck squamous cell carcinoma patients receiving neck dissection: A pilot study. Clin Transl Radiat Oncol. 2020; 23:60–64. https://doi.org/10.1016/j.ctro.2020.04.008 [PubMed]
7. Otsubo Y, Hashimoto K, Kanbe T, Sumi M, Moriuchi H. Association of cord blood chemokines and other biomarkers with neonatal complications following intrauterine inflammation. PLoS One. 2017; 12:e0175082. https://doi.org/10.1371/journal.pone.0175082 [PubMed]
8. Zlotnik A, Yoshie O. Chemokines: a new classification system and their role in immunity. Immunity. 2000; 12:121–27. https://doi.org/10.1016/s1074-7613(00)80165-x [PubMed]
9. Hu DD, Li PC, He YF, Jia W, Hu B. Overexpression of Coiled-Coil Domain-Containing Protein 34 (CCDC34) and its Correlation with Angiogenesis in Esophageal Squamous Cell Carcinoma. Med Sci Monit. 2018; 24:698–705. https://doi.org/10.12659/msm.908335 [PubMed]
10. Li PC, Hu DD, Jia W, Hu B. Expression and Association of Tumor Necrosis Factor Receptor Associated Factor 4 (TRAF4) in Esophageal Squamous Cell Carcinoma. Med Sci Monit. 2019; 25:2368–76. https://doi.org/10.12659/MSM.915474 [PubMed]
11. D’Agostino G, Cecchinato V, Uguccioni M. Chemokine Heterocomplexes and Cancer: A Novel Chapter to Be Written in Tumor Immunity. Front Immunol. 2018; 9:2185. https://doi.org/10.3389/fimmu.2018.02185 [PubMed]
12. Seo EH, Namgung JH, Oh CS, Kim SH, Lee SH. Association of Chemokines and Chemokine Receptor Expression with Monocytic-Myeloid-Derived Suppressor Cells during Tumor Progression. Immune Netw. 2018; 18:e23. https://doi.org/10.4110/in.2018.18.e23 [PubMed]
13. Susek KH, Karvouni M, Alici E, Lundqvist A. The Role of CXC Chemokine Receptors 1-4 on Immune Cells in the Tumor Microenvironment. Front Immunol. 2018; 9:2159. https://doi.org/10.3389/fimmu.2018.02159 [PubMed]
14. Chiba Y, Mizoguchi I, Hasegawa H, Ohashi M, Orii N, Nagai T, Sugahara M, Miyamoto Y, Xu M, Owaki T, Yoshimoto T. Regulation of myelopoiesis by proinflammatory cytokines in infectious diseases. Cell Mol Life Sci. 2018; 75:1363–76. https://doi.org/10.1007/s00018-017-2724-5 [PubMed]
15. Wei LY, Lee JJ, Yeh CY, Yang CJ, Kok SH, Ko JY, Tsai FC, Chia JS. Reciprocal activation of cancer-associated fibroblasts and oral squamous carcinoma cells through CXCL1. Oral Oncol. 2019; 88:115–23. https://doi.org/10.1016/j.oraloncology.2018.11.002 [PubMed]
16. Sun W, Qiu Z, Huang W, Cao M. Gene expression profiles and protein-protein interaction networks during tongue carcinogenesis in the tumor microenvironment. Mol Med Rep. 2018; 17:165–71. https://doi.org/10.3892/mmr.2017.7843 [PubMed]
17. Dang H, Wu W, Wang B, Cui C, Niu J, Chen J, Chen Z, Liu Y. CXCL5 Plays a Promoting Role in Osteosarcoma Cell Migration and Invasion in Autocrine- and Paracrine-Dependent Manners. Oncol Res. 2017; 25:177–86. https://doi.org/10.3727/096504016X14732772150343 [PubMed]
18. Dufies M, Grytsai O, Ronco C, Camara O, Ambrosetti D, Hagege A, Parola J, Mateo L, Ayrault M, Giuliano S, Grépin R, Lagarde N, Montes M, et al. New CXCR1/CXCR2 inhibitors represent an effective treatment for kidney or head and neck cancers sensitive or refractory to reference treatments. Theranostics. 2019; 9:5332–46. https://doi.org/10.7150/thno.34681 [PubMed]
19. Wang N, Feng Y, Wang Q, Liu S, Xiang L, Sun M, Zhang X, Liu G, Qu X, Wei F. Neutrophils infiltration in the tongue squamous cell carcinoma and its correlation with CEACAM1 expression on tumor cells. PLoS One. 2014; 9:e89991. https://doi.org/10.1371/journal.pone.0089991 [PubMed]
20. Łukaszewicz-Zając M, Pączek S, Muszyński P, Kozłowski M, Mroczko B. Comparison between clinical significance of serum CXCL-8 and classical tumor markers in oesophageal cancer (OC) patients. Clin Exp Med. 2019; 19:191–99. https://doi.org/10.1007/s10238-019-00548-9 [PubMed]
21. Griffith JW, Sokol CL, Luster AD. Chemokines and chemokine receptors: positioning cells for host defense and immunity. Annu Rev Immunol. 2014; 32:659–702. https://doi.org/10.1146/annurev-immunol-032713-120145 [PubMed]
22. O’Brien JC, Rainwater YB, Malviya N, Cyrus N, Auer-Hackenberg L, Hynan LS, Hosler GA, Jacobe HT. Transcriptional and Cytokine Profiles Identify CXCL9 as a Biomarker of Disease Activity in Morphea. J Invest Dermatol. 2017; 137:1663–70. https://doi.org/10.1016/j.jid.2017.04.008 [PubMed]
23. Tan S, Wang K, Sun F, Li Y, Gao Y. CXCL9 promotes prostate cancer progression through inhibition of cytokines from T cells. Mol Med Rep. 2018; 18:1305–10. https://doi.org/10.3892/mmr.2018.9152 [PubMed]
24. Li Z, Qian J, Li J, Zhu C. Clinical Significance of Serum Chemokines in Esophageal Cancer. Med Sci Monit. 2019; 25:5850–55. https://doi.org/10.12659/MSM.916846 [PubMed]
25. Samarendra H, Jones K, Petrinic T, Silva MA, Reddy S, Soonawalla Z, Gordon-Weeks A. A meta-analysis of CXCL12 expression for cancer prognosis. Br J Cancer. 2017; 117:124–35. https://doi.org/10.1038/bjc.2017.134 [PubMed]
26. Garg R, Blando JM, Perez CJ, Abba MC, Benavides F, Kazanietz MG. Protein Kinase C Epsilon Cooperates with PTEN Loss for Prostate Tumorigenesis through the CXCL13-CXCR5 Pathway. Cell Rep. 2017; 19:375–88. https://doi.org/10.1016/j.celrep.2017.03.042 [PubMed]
27. Parikh A, Shin J, Faquin W, Lin DT, Tirosh I, Sunwoo JB, Puram SV. Malignant cell-specific CXCL14 promotes tumor lymphocyte infiltration in oral cavity squamous cell carcinoma. J Immunother Cancer. 2020; 8:e001048. https://doi.org/10.1136/jitc-2020-001048 [PubMed]
28. Xu YX, Sun J, Xiao WL, Liu YS, Yue J, Xue LF, Deng J, Zhi KQ, Wang YL. MiR-4513 mediates the proliferation and apoptosis of oral squamous cell carcinoma cells via targeting CXCL17. Eur Rev Med Pharmacol Sci. 2019; 23:3821–28. https://doi.org/10.26355/eurrev_201905_17809 [PubMed]
29. Ye H, Yu T, Temam S, Ziober BL, Wang J, Schwartz JL, Mao L, Wong DT, Zhou X. Transcriptomic dissection of tongue squamous cell carcinoma. BMC Genomics. 2008; 9:69. https://doi.org/10.1186/1471-2164-9-69 [PubMed]
30. Pyeon D, Newton MA, Lambert PF, den Boon JA, Sengupta S, Marsit CJ, Woodworth CD, Connor JP, Haugen TH, Smith EM, Kelsey KT, Turek LP, Ahlquist P. Fundamental differences in cell cycle deregulation in human papillomavirus-positive and human papillomavirus-negative head/neck and cervical cancers. Cancer Res. 2007; 67:4605–19. https://doi.org/10.1158/0008-5472.CAN-06-3619 [PubMed]
31. Ginos MA, Page GP, Michalowicz BS, Patel KJ, Volker SE, Pambuccian SE, Ondrey FG, Adams GL, Gaffney PM. Identification of a gene expression signature associated with recurrent disease in squamous cell carcinoma of the head and neck. Cancer Res. 2004; 64:55–63. https://doi.org/10.1158/0008-5472.can-03-2144 [PubMed]
32. Peng CH, Liao CT, Peng SC, Chen YJ, Cheng AJ, Juang JL, Tsai CY, Chen TC, Chuang YJ, Tang CY, Hsieh WP, Yen TC. A novel molecular signature identified by systems genetics approach predicts prognosis in oral squamous cell carcinoma. PLoS One. 2011; 6:e23452. https://doi.org/10.1371/journal.pone.0023452 [PubMed]
33. Estilo CL, O-charoenrat P, Talbot S, Socci ND, Carlson DL, Ghossein R, Williams T, Yonekawa Y, Ramanathan Y, Boyle JO, Kraus DH, Patel S, Shaha AR, et al. Oral tongue cancer gene expression profiling: Identification of novel potential prognosticators by oligonucleotide microarray analysis. BMC Cancer. 2009; 9:11. https://doi.org/10.1186/1471-2407-9-11 [PubMed]
34. Talbot SG, Estilo C, Maghami E, Sarkaria IS, Pham DK, O-charoenrat P, Socci ND, Ngai I, Carlson D, Ghossein R, Viale A, Park BJ, Rusch VW, Singh B. Gene expression profiling allows distinction between primary and metastatic squamous cell carcinomas in the lung. Cancer Res. 2005; 65:3063–71. https://doi.org/10.1158/0008-5472.CAN-04-1985 [PubMed]
35. Sengupta S, den Boon JA, Chen IH, Newton MA, Dahl DB, Chen M, Cheng YJ, Westra WH, Chen CJ, Hildesheim A, Sugden B, Ahlquist P. Genome-wide expression profiling reveals EBV-associated inhibition of MHC class I expression in nasopharyngeal carcinoma. Cancer Res. 2006; 66:7999–8006. https://doi.org/10.1158/0008-5472.CAN-05-4399 [PubMed]
36. Yu G, Wang LG, Han Y, He QY. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS. 2012; 16:284–87. https://doi.org/10.1089/omi.2011.0118 [PubMed]
37. Wolff HA, Rolke D, Rave-Fränk M, Schirmer M, Eicheler W, Doerfler A, Hille A, Hess CF, Matthias C, Rödel RM, Christiansen H. Analysis of chemokine and chemokine receptor expression in squamous cell carcinoma of the head and neck (SCCHN) cell lines. Radiat Environ Biophys. 2011; 50:145–54. https://doi.org/10.1007/s00411-010-0341-x [PubMed]
38. Koringa PG, Jakhesara SJ, Bhatt VD, Meshram CP, Patel AK, Fefar DT, Joshi CG. Comprehensive transcriptome profiling of squamous cell carcinoma of horn in Bos indicus. Vet Comp Oncol. 2016; 14:122–36. https://doi.org/10.1111/vco.12079 [PubMed]
39. Yang B, Dong K, Guo P, Guo P, Jie G, Zhang G, Li T. Identification of Key Biomarkers and Potential Molecular Mechanisms in Oral Squamous Cell Carcinoma by Bioinformatics Analysis. J Comput Biol. 2020; 27:40–54. https://doi.org/10.1089/cmb.2019.0211 [PubMed]
40. Li C, Zhou Y, Liu J, Su X, Qin H, Huang S, Huang X, Zhou N. Potential Markers from Serum-Purified Exosomes for Detecting Oral Squamous Cell Carcinoma Metastasis. Cancer Epidemiol Biomarkers Prev. 2019; 28:1668–81. https://doi.org/10.1158/1055-9965.EPI-18-1122 [PubMed]
41. Miyamoto C, Maehata Y, Motohashi K, Ozawa S, Ikoma T, Hidaka K, Wada-Takahashi S, Takahashi SS, Yoshino F, Yoshida A, Kubota E, Hata R, Lee MC. Fasudil, a Rho kinase inhibitor, suppresses tumor growth by inducing CXCL14/BRAK in head and neck squamous cell carcinoma. Biomed Res. 2014; 35:381–88. https://doi.org/10.2220/biomedres.35.381 [PubMed]
42. Szabo P, Valach J, Smetana K Jr, Dvořánková B. Comparative analysis of IL-8 and CXCL-1 production by normal and cancer stromal fibroblasts. Folia Biol (Praha). 2013; 59:134–37. [PubMed]
43. Rhodes DR, Yu J, Shanker K, Deshpande N, Varambally R, Ghosh D, Barrette T, Pandey A, Chinnaiyan AM. ONCOMINE: a cancer microarray database and integrated data-mining platform. Neoplasia. 2004; 6:1–6. https://doi.org/10.1016/s1476-5586(04)80047-2 [PubMed]
44. Tang Z, Li C, Kang B, Gao G, Li C, Zhang Z. GEPIA: a web server for cancer and normal gene expression profiling and interactive analyses. Nucleic Acids Res. 2017; 45:W98–102. https://doi.org/10.1093/nar/gkx247 [PubMed]
45. Gao J, Aksoy BA, Dogrusoz U, Dresdner G, Gross B, Sumer SO, Sun Y, Jacobsen A, Sinha R, Larsson E, Cerami E, Sander C, Schultz N. Integrative analysis of complex cancer genomics and clinical profiles using the cBioPortal. Sci Signal. 2013; 6:pl1. https://doi.org/10.1126/scisignal.2004088 [PubMed]
46. Cerami E, Gao J, Dogrusoz U, Gross BE, Sumer SO, Aksoy BA, Jacobsen A, Byrne CJ, Heuer ML, Larsson E, Antipin Y, Reva B, Goldberg AP, et al. The cBio cancer genomics portal: an open platform for exploring multidimensional cancer genomics data. Cancer Discov. 2012; 2:401–04. https://doi.org/10.1158/2159-8290.CD-12-0095 [PubMed]
47. Szklarczyk D, Gable AL, Lyon D, Junge A, Wyder S, Huerta-Cepas J, Simonovic M, Doncheva NT, Morris JH, Bork P, Jensen LJ, Mering CV. STRING v11: protein-protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets. Nucleic Acids Res. 2019; 47:D607–13. https://doi.org/10.1093/nar/gky1131 [PubMed]
48. Szklarczyk D, Morris JH, Cook H, Kuhn M, Wyder S, Simonovic M, Santos A, Doncheva NT, Roth A, Bork P, Jensen LJ, von Mering C. The STRING database in 2017: quality-controlled protein-protein association networks, made broadly accessible. Nucleic Acids Res. 2017; 45:D362–68. https://doi.org/10.1093/nar/gkw937 [PubMed]
49. Szklarczyk D, Franceschini A, Wyder S, Forslund K, Heller D, Huerta-Cepas J, Simonovic M, Roth A, Santos A, Tsafou KP, Kuhn M, Bork P, Jensen LJ, von Mering C. STRING v10: protein-protein interaction networks, integrated over the tree of life. Nucleic Acids Res. 2015; 43:D447–52. https://doi.org/10.1093/nar/gku1003 [PubMed]
50. Franceschini A, Lin J, von Mering C, Jensen LJ. SVD-phy: improved prediction of protein functional associations through singular value decomposition of phylogenetic profiles. Bioinformatics. 2016; 32:1085–87. https://doi.org/10.1093/bioinformatics/btv696 [PubMed]
51. Franceschini A, Szklarczyk D, Frankild S, Kuhn M, Simonovic M, Roth A, Lin J, Minguez P, Bork P, von Mering C, Jensen LJ. STRING v9.1: protein-protein interaction networks, with increased coverage and integration. Nucleic Acids Res. 2013; 41:D808–15. https://doi.org/10.1093/nar/gks1094 [PubMed]
52. Szklarczyk D, Franceschini A, Kuhn M, Simonovic M, Roth A, Minguez P, Doerks T, Stark M, Muller J, Bork P, Jensen LJ, von Mering C. The STRING database in 2011: functional interaction networks of proteins, globally integrated and scored. Nucleic Acids Res. 2011; 39:D561–68. https://doi.org/10.1093/nar/gkq973 [PubMed]
53. Jensen LJ, Kuhn M, Stark M, Chaffron S, Creevey C, Muller J, Doerks T, Julien P, Roth A, Simonovic M, Bork P, von Mering C. STRING 8--a global view on proteins and their functional interactions in 630 organisms. Nucleic Acids Res. 2009; 37:D412–16. https://doi.org/10.1093/nar/gkn760 [PubMed]
54. von Mering C, Jensen LJ, Kuhn M, Chaffron S, Doerks T, Krüger B, Snel B, Bork P. STRING 7--recent developments in the integration and prediction of protein interactions. Nucleic Acids Res. 2007; 35:D358–62. https://doi.org/10.1093/nar/gkl825 [PubMed]
55. von Mering C, Jensen LJ, Snel B, Hooper SD, Krupp M, Foglierini M, Jouffre N, Huynen MA, Bork P. STRING: known and predicted protein-protein associations, integrated and transferred across organisms. Nucleic Acids Res. 2005; 33:D433–37. https://doi.org/10.1093/nar/gki005 [PubMed]
56. von Mering C, Huynen M, Jaeggi D, Schmidt S, Bork P, Snel B. STRING: a database of predicted functional associations between proteins. Nucleic Acids Res. 2003; 31:258–61. https://doi.org/10.1093/nar/gkg034 [PubMed]
57. Snel B, Lehmann G, Bork P, Huynen MA. STRING: a web-server to retrieve and display the repeatedly occurring neighbourhood of a gene. Nucleic Acids Res. 2000; 28:3442–44. https://doi.org/10.1093/nar/28.18.3442 [PubMed]
58. Zhou Y, Zhou B, Pache L, Chang M, Khodabakhshi AH, Tanaseichuk O, Benner C, Chanda SK. Metascape provides a biologist-oriented resource for the analysis of systems-level datasets. Nat Commun. 2019; 10:1523. https://doi.org/10.1038/s41467-019-09234-6 [PubMed]
59. Sun D, Wang J, Han Y, Dong X, Ge J, Zheng R, Shi X, Wang B, Li Z, Ren P, Sun L, Yan Y, Zhang P, et al. TISCH: a comprehensive web resource enabling interactive single-cell transcriptome visualization of tumor microenvironment. Nucleic Acids Res. 2021; 49:D1420–30. https://doi.org/10.1093/nar/gkaa1020 [PubMed]
60. Nagy Á, Lánczky A, Menyhárt O, Győrffy B. Author Correction: Validation of miRNA prognostic power in hepatocellular carcinoma using expression data of independent datasets. Sci Rep. 2018; 8:11515. https://doi.org/10.1038/s41598-018-29514-3 [PubMed]
61. Li T, Fu J, Zeng Z, Cohen D, Li J, Chen Q, Li B, Liu XS. TIMER2.0 for analysis of tumor-infiltrating immune cells. Nucleic Acids Res. 2020; 48:W509–14. https://doi.org/10.1093/nar/gkaa407 [PubMed]
62. Li T, Fan J, Wang B, Traugh N, Chen Q, Liu JS, Li B, Liu XS. TIMER: A Web Server for Comprehensive Analysis of Tumor-Infiltrating Immune Cells. Cancer Res. 2017; 77:e108–10. https://doi.org/10.1158/0008-5472.CAN-17-0307 [PubMed]
63. Li B, Severson E, Pignon JC, Zhao H, Li T, Novak J, Jiang P, Shen H, Aster JC, Rodig S, Signoretti S, Liu JS, Liu XS. Comprehensive analyses of tumor immunity: implications for cancer immunotherapy. Genome Biol. 2016; 17:174. https://doi.org/10.1186/s13059-016-1028-7 [PubMed]
64. Vasaikar SV, Straub P, Wang J, Zhang B. LinkedOmics: analyzing multi-omics data within and across 32 cancer types. Nucleic Acids Res. 2018; 46:D956–63. https://doi.org/10.1093/nar/gkx1090 [PubMed]

Research Paper Volume 13, Issue 13 pp 17789—17817