Identification of Four Key Prognostic Genes and Three Potential Drugs in Human Papillomavirus Negative Head and Neck Squamous Cell Carcinoma

Abstract

Background

Head and neck squamous cell carcinoma (HNSCC) is a common tumor worldwide with poor prognosis. The pathogenesis of human papillomavirus (HPV)-positive and HPV-negative HNSCCs differs. However, few studies have considered the HPV status when identifying biomarkers for HNSCC. Thus, the identification of biomarkers for HPV-positive and HPV-negative HNSCCs is urgently needed.

Methods

Three microarray datasets from Gene Expression Omnibus (GEO) were analyzed, and the differentially expressed genes (DEGs) were obtained. Then, functional enrichment pathway analysis was performed and protein-protein interaction (PPI) networks were constructed. The expression of hub genes at both the mRNA and protein level was determined in Oncomine, The Cancer Genome Atlas (TCGA) and the Human Protein Atlas (HPA). In addition, survival analysis of the patient stratified by HPV status and the expression levels of key genes were performed based on TCGA data. The role of AREG, STAG3, CAV1 and C19orf57 in cancer were analyzed through Gene set enrichment analysis (GSEA). Finally, the top ten small molecule drugs were identified and the therapeutic value of zonisamide, NVP-AUY922, PP-2 and fostamatinib was further evaluated in six HPV-negative HNSCC cell lines.

Results

In total, 47 DEGs were obtained, 11 of which were identified as hub genes. Biological process analysis indicated that the hub genes were associated with the G1/S transition of the mitotic cell cycle. Survival analysis uncovered that the prognostic value of AREG, STAG3, C19orf57 and CAV1 differed between HPV-positive and HPV-negative patients. Gene set enrichment analysis (GSEA) showed the role of AREG, STAG3 and CAV1 in dysregulated pathways of tumor. Ten small molecules were identified as potential drugs specifically for HPV-positive or HPV-negative patients; three—NVP-AUY922, fostamatinib and PP-2—greatly inhibited the proliferation of six HPV-negative HNSCC cell lines in vitro, and NVP-AUY922 inhibited three HPV-negative HNSCC xenografts in vivo.

Conclusions

In conclusion, AREG, STAG3, C19orf57 and CAV1 are key prognostic factors and potential therapeutic targets in HPV-negative HNSCC. NVP-AUY922, fostamatinib and PP-2 may be effective drugs for HPV-negative HNSCC.

Background

Head and neck cancer is the ninth most common among all malignancies worldwide [1]. Head and neck squamous cell carcinomas (HNSCCs) account for approximately 90% of all head and neck cancers [2]. Despite considerable efforts, the five-year overall survival (OS) rate remains at 33%~71%. Local recurrence, distant metastasis, and drug resistance are the main causes of death [3]. Approximately 75% of HNSCCs are associated with the consumption of tobacco and alcohol; however, some HNSCCs, particularly oropharyngeal tumors, are caused by human papillomavirus (HPV) infection [4].

HPV is a DNA virus that can infect skin and mucous membranes [5]. HNSCCs are divided into HPV-related HNSCCs and HPV-unrelated HNSCCs according to the infection status of HPV. HPV-unrelated HNSCCs have a worse therapeutic response than HPV-related HNSCCs [6]. To determine the pathogenesis of HNSCC, some biomarkers have been identified. However, few studies have considered the HPV status, which may cause unreliable results. Thus, the identification of reliable markers for HPV-positive and HPV-negative tumors is urgently needed to establish effective diagnostic, prognostic, and therapeutic strategies.

During the past years, integrated bioinformatic methods have been established to analyze high-throughput data across platforms, which helps us to identify differentially expressed genes (DEGs) and their related functional pathways involved in HNSCC carcinogenesis and progression. Gene expression omnibus (GEO) contains a large of high-throughput functional genomic datasets [7]. Some key genes in HNSCC, such as FN1, APP, SERPINE1, PLAU and ACTA1, were identified through GEO datasets [8, 9].

This study aimed to identify biomarkers for HPV-positive as well as HPV-negative HNSCCs via integrated bioinformatic methods. To avoid bias, DEGs were screened in three GEO datasets. Then Gene Ontology (GO) term and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were performed and protein-protein interaction (PPI) networks were constructed to analyze the difference in the molecular mechanisms active in HPV-positive and HPV-negative HNSCC tissues. The hub genes were determined, and survival analysis suggested that amphiregulin (AREG), stromal antigen 3 (STAG3), chromosome 19 open reading frame 57 (C19orf57) and caveolin-1 (CAV1) are key prognostic biomarkers for HPV-negative HNSCCs. Gene set enrichment analysis showed the role of AREG, STAG3 and CAV1 in dysregulated pathways of tumor. Finally, 10 small molecules were screened as potential drugs. The therapeutic value of zonisamide, NVP-AUY922, PP-2 and fostamatinib was further evaluated in six HPV-negative HNSCC cell lines, and NVP-AUY922 was tested therapeutic value in three HPV-negative HNSCC xenografts in vivo.

Methods

Acquisition of the datasets and identification of the DEGs

Datasets were searched in GEO (http://www.ncbi.nlm.nih.gov/geo) and the three datasets with the largest sample size (GSE39366,GSE40774,GSE55550 [10-12]) were sorted for further analysis. GEO2R is a web analysis tool for identifying DEGs across different experimental conditions. HPV 16-positive and HPV 16-negative samples (the HPV-inactive [DNA+RNA−] samples in GSE55550 were excluded to obtain the most significant results) were separated into two groups and analyzed with GEO2R. The Benjamini-Hochberg (false discovery rate, FDR) method was applied to adjust the P-values (adjusted P-value, adj.P-value). The probe identifiers were converted to gene symbols, and gene symbols with duplication or loss were deleted. The DEGs meeting the criteria of an adj.P-value of <0.05 and a |log (FC) | of ≥1 were considered statistically significant DEGs. Metascape (https://metascape.org) is a web-based portal designed to provide a comprehensive gene list annotation and analysis resource [13]. Three groups of DEGs were uploaded to Metascape for meta-analysis to explore the most highly enriched pathways and the overlapping genes of the three datasets were visualized.

Go term and KEGG pathway enrichment analyses

The GO database is the largest database worldwide that provides information on the functions of genes, including their biological processes (BPs), cellular components (CCs) and molecular functions (MFs) [14]. KEGG is a web database for exploring advanced functions of biological systems [15]. The Database for Annotation, Visualization and Integrated Discovery (DAVID, https://david.ncifcrf.gov) tool is an online tool comprising integrated biological knowledge bases and analytical tools [16]. The overlapping DEGs were submitted to GO and KEGG pathway analyses by DAVID (version 6.8). To identify the direct interactions among the DEGs, “GOTERM_DIRECT” BP, CC and MF categories were selected, and a P-value of <0.05 was considered statistically significant.

Construction of the PPI network and identification of hub genes

The Search Tool for the Retrieval of Interacting Genes/Proteins (STRING, https://string-db.org/) database is widely used to predict PPI networks, as a biological database, it integrates PPI information from publicly available sources [17]. To examine the functional connections among the DEGs, STRING (version 11.0) was used to determine the direct interactions among the DEGs. An interaction score of > 0.4 was considered statistically significant. The PPI network was visualized by platform Cytoscape (version 3.8.0) [18]. The Cytoscape plugin Molecular Complex Detection (MCODE) [19] was used to identify the top two modules in the PPI network with the following parameters: degree cutoff=2, node density cutoff=0.2, node score cutoff=0.1, K-core=2, and max.depth=100. In addition, the cytoHubba plugin was used to calculate the degree of each protein node, and the top 11 genes were identified as hub genes. The expression of these hub genes in the three GSE datasets was visualized with R software (version 3.6.3).

Verification of hub gene mRNA expression levels

Oncomine (https://www.oncomine.org/) is a cancer microarray database and integrated data mining platform [20]. The mRNA expression levels of the 11 hub genes were verified in Oncomine, and only two HNSCC datasets (Zhai Cervix [21] and Slebos Head-Neck [22]) met the criteria of including the HPV infection status. Differential gene expression was visualized with the R software.

The Cancer Genome Atlas (TCGA, https://www.cancer.gov/tcga.) contains many patient samples across 33 cancer types, with clinical information as well as genome sequencing data. In addition, University of California Santa Cruz (UCSC) Xena (https://xenabrowser.net/) is an online database that could be used to explore correlations between genomic and phenotypic variables [23], and TCGA data are included. To verify mRNA expression, HNSCC samples and mRNA expression data of the hub genes were downloaded from TCGA by UCSC Xena. Data for 604 samples, including HNSCC and solid normal tissues were obtained. Samples without mRNA expression data were removed, and the differential expression of the key genes was analyzed using GraphPad Prism software (version 8.0.0). The Mann-Whitney U test was performed, and a two-tailed exact P-value of <0.05 was considered statistically significant.

Immunohistochemical analysis

The Human Protein Atlas (HPA, http://www.proteinatlas.org,version 19.3) is the largest database mapping human proteins [24]. Since HPV status data of HNSCC samples were not available in the HPA, immunohistochemical data for the key genes in HPV-positive and HPV-negative samples were not provided. Thus, immunohistochemical analysis of HNSCC and normal tissues (all normal tissues were oral mucosal tissues) was performed.

Survival analysis based on TCGA database data

Clinical information for the patients was downloaded from UCSC Xena, and the OS of patients stratified by the mRNA levels of the 11 key genes was analyzed. High expression or low expression was defined as an expression level of greater than or less than the median value, respectively. A log-rank (Mantel-Cox) test was performed, and a P-value of <0.05 was considered to indicate a significant difference. For the comparison between the HPV-positive (+) and the HPV-negative (-) groups, the hazard ratio (HR) was calculated by the log-rank test as follows: HR = expression in the HPV (+) group /expression in the HPV (-) group. For comparisons between other groups, the HR was calculated as follows: HR=high expression/low expression. The disease-free interval (DFI), disease-specific survival (DSS) and progression-free interval (PFI) of patients stratified by the expression levels of genes with statistical significance in the OS analysis were analyzed.

GSEA of AREG, CAV1, STAG3 and C19orf57

GSEA software (version 4.0.3, https://www.gsea-msigdb.org/gsea/downloads.jsp) was used to explore the association between gene expression and dysregulated pathways in cancer. The 73 samples from the TCGA database (HPV-negative) were divided into high and low expression groups based on the levels of AREG, CAV1, STAG3 and C19orf57, with the median value of each gene signifying the cutoff value. The predefined gene sets “c2.all.v7.1.symbols.gmt”, “c4.all.v7.1.symbols.gmt” and “c6.all.v7.1symbols.gmt” from the Molecular Signatures Database (MSigDB) were used for analyses. The normalized enrichment score (NES) was calculated, the cutoff for significance was defined as follows: a nominal p-value of <0.05 and FDR q-value of <0.25.

Identification of potential small molecule drugs

The Connectivity Map (https://clue.io/cmap) is a resource that uses cellular responses to perturbations to find relationships between diseases, genes, and therapeutics. The CLUE platform (https://clue.io) was used to identify potential small molecules based on data from 9 tumor cell lines [25]. The DEGs co-expressed in the three GEO datasets were divided into up- and downregulated genes and were then uploaded to run the query. The parameters were defined as follows: Gene expression (L1000) and Touchstone. After calculation, each small molecule was assigned an enrichment score ranging from -100 to 100. A positive score indicates similarity between the signature of a given perturbagen and that of the query, while a negative score indicates a disparity between the two signatures. Scores of +90 or higher and of -90 or lower were considered strong scores for developing hypotheses for further study. Here, the molecules with the top five positive and negative scores were identified for further study. The 2D structures of these molecules are available in PubChem (https://pubchem.ncbi.nlm.nih.gov/), a public database of small molecules and their biological properties [26].

Evaluation of the therapeutic value of small molecule drugs in vitro

HN4, HN6, HN30, SCC9, SCC25 and CAL27 cell lines were obtained from the Shanghai Key Laboratory of Stomatology; zonisamide, NVP-AUY922, PP-2 and fostamatinib were purchased from Absin Bioscience Inc. Zonisamide and PP-2 were dissolved in Dimethyl sulfoxide (DMSO) and diluted with medium with a gradient of 5µM. Similarly, NVP-AUY922 and fostamatinib were diluted with medium with a gradient of 5nM. Approximately 3×10³ cells per well were seeded in 96-well culture plates. After culture for 24 h, the medium was removed. The four drugs at different concentrations were added with the volume of 100 µl. For the negative control wells, drug-free medium was added. In addition, the background control wells were added only drug-free medium but without cells. After culture for 72 h, the medium was removed, 100 µl of fresh medium containing 10% Cell Counting Kit-8 (CCK8, Dojindo, Japan) solution was added. The medium in the control wells was replaced by fresh medium. After incubation for 3 h, the absorbance at 450 nm was measured in a multimode microplate reader (SpectraMax i3, Molecular Devices). The experiments were repeated in triplicate. Finally, the cell inhibition rate was calculated with the following formula: %cell inhibition rate = 100- (mean optical density (OD) of the treatment wells - mean OD of the background control wells)/(mean OD of the negative control wells – mean OD of the background control wells)×100, and the half-maximal inhibitory concentrations (IC50) of the drugs were calculated with GraphPad Prism software (version 8.0.0).

In vivo studies

Female BALB/c nude mice (6 weeks of age, weighing approximately 20g) were bred in SPF facilities of Shanghai Ninth People’s hospital. The experiment procedures were approved by the Laboratory Animal Care and Use Committees of the hospital. In the present study, HNSCC tumor models were established by inoculation 2×10⁶ HN6, SCC9, and SCC25 cells into the dorsal flank of the mice. When the tumor size reached about 5 mm, the mice were randomly assigned into control and NVP-AUY922 treatment group (n = 3 per group). Vehicle (10% DMSO, 5% Tween 20, 85% saline) and 50 mg/kg NVP-AUY922 were given intraperitoneal to each group three times a week for two weeks, and tumor size were measured by vernier calipers three times a week. The tumor volumes were calculated with the following formula: (length×width²)/2. Mice were sacrificed after 2 weeks.

Results

Identification of DEGs and meta-analysis of three datasets

Three GEO datasets (GSE39366, GSE40774, GSE55550) containing a total of 345 samples were analyzed by GEO2R, with the data divided into the HPV-positive and HPV-negative groups. The details of the three datasets was summarized in Table 1. DEGs were identified as genes meeting the following criteria: adj.P-value <0.05 and |log(FC)| ≥1. The screen identified 259 DEGs in GSE39366, 428 in GSE40774 and 726 in GSE55550, and 47 overlapping genes were identified across the three datasets (Figs. 1a-d). Metascape was used to perform a meta-analysis on the DEGs from three datasets. The overlapping DEGs and the top 20 enriched GO terms were visualized in Circos plot and heatmap, respectively, and the most enriched GO term was leukocyte migration (Figs. 1e and f).

GO and KEGG enrichment analysis of the 47 overlapping DEGs

DAVID was used to further analyze the BP and functional pathways enriched with the DEGs. The main enriched BPs were mammary gland development and leukocyte migration, similar to the results of the meta-analysis of all DEGs described above. In addition, positive regulation of the Toll-like receptor 3 signaling pathway was associated with pathogen recognition and activation of innate immunity, while the G1/S transition of the mitotic cell cycle was associated with carcinogenesis. These results showed that the DEGs are essential in tumor initiation and development (Fig. 2a). CC analysis indicated that the DEGs were most enriched in the extracellular space. Regarding the MF classification, the DEGs were enriched mainly in cytokine activity and protein kinase binding. KEGG analysis indicated that the DEGs were enriched mainly in cell cycle and Hippo signaling pathways, which regulate cell proliferation and apoptosis (Fig. 2b).

Construction of the PPI network and identification of hub genes

STRING database was used to construct the PPI network between DEGs (Fig. 2c). The network included 47 nodes and 23 edges, and the PPI enrichment P-value was 4.91e-05. The top 2 significant modules were identified (Figs. 2d and e). Furthermore, the hub genes were screened with cytoHubba, and since 4 genes occupied the same rank, 11 genes were finally identified as hub genes (Fig. 2f). CDKN2A was the top-ranked gene, followed by C19orf57, ZCWPW1 and KRT19. The expression levels of these 11 hub genes in the GEO datasets were visualized with heatmaps (Figs. 3a-c). NRG1, AREG, CAV1 and PTHLH were downregulated in HPV-positive HNSCCs, while the remaining genes were upregulated. The opposite pattern was observed in HPV-negative HNSCC samples.

Verification of hub gene expression patterns at the mRNA level based on Oncomine and TCGA data

To verify the expression of the hub genes, two datasets from the Oncomine database were analyzed (Fig. 3d). The four genes identified as downregulated in the three GEO datasets also had lower expression levels in Oncomine datasets. In addition, TCGA data were downloaded to confirm the mRNA expression patterns of the key genes between HPV-positive and HPV-negative HNSCCs as well as between HNSCC and normal tissues. The results confirmed that all 11 hub genes were significantly differentially expressed between HPV-positive and -negative HNSCC tissues. In addition, similarly, CAV1, PTHLH, NRG1 and AREG were downregulated but the other genes were upregulated in HPV-positive HNSCCs (Fig. 4).

The expression patterns of hub genes at the protein level

The expression of hub genes at the protein level was investigated in the HPA. Since the HPV status was not available in the HPA, the expression patterns of ten key genes between HNSCCs and normal samples were analyzed by immunohistochemistry. Except for the result of ZCWPW1 staining, which was not provided, as shown in figure 5, the protein expression levels of CDKN2A, KRT19, VCAM1 and CAV1 in HNSCCs were higher than those in normal tissues, although staining for C19orf57 was enhanced in normal tissues. In addition, no differences were observed for STAG3, NRG1, PTHLH and AREG staining. These results indicated that the protein expression levels of the hub genes between HNSCCs and normal tissues were consistent with their mRNA expression levels.

AREG, STAG3, C19orf57 and CAV1 were associated with the prognosis of HNSCC

To further explore the prognostic value of the hub genes, survival analysis was performed. First, the OS, DFI, DSS and PFI of HPV-positive and HPV-negative patients were compared. HPV-positive patients exhibited higher OS rates, longer PFI times and higher DSS rates (Fig. 7 and Additional file 1: Figure S1). Then, OS analysis of patients stratified by the expression levels of the 11 key genes revealed that AREG and CAV1 were associated with poor prognosis in HNSCC patients, C19orf57 and STAG3 were associated with favorable prognosis, and the other 7 genes were not significantly associated with OS (Additional file 2: Figure S2). Interestingly, AREG was also a prognostic factor in HPV-negative patients, and the OS rate was higher for patients with lower AREG expression levels (Fig. 6a and Fig. 7). Finally, DFI, DSS and PFI were further analyzed for patients stratified by the expression levels of the four genes identified as statistically significant in the OS analysis (Fig. 7 and Additional file 3: Figure S3). These results indicated that the AREG gene had the highest prognostic value for DFI, DSS and PFI in patients with HNSCC and that it was a prognostic factor for favorable PFI in HPV-positive patients. Additionally, for HPV-negative patients, AREG could be a prognostic factor for unfavorable DSS and PFI. Moreover, STAG3 and C19orf57 were prognostic factors for favorable DSS and PFI in HNSCCs. STAG3 was the only prognostic factor for unfavorable DFI in HPV-negative patients, while C19orf57 could be used as an indicator of poor PFI in HPV-positive patients. Similar results were found for CAV1, which was a factor for unfavorable DSS in HNSCCs as well as for unfavorable DSS and PFI in HPV-negative patients (Fig. 6b, 7 and Additional file 1: Figure S1). Overall, AREG, STAG3, C19orf57 and CAV1 were the four key factors with high prognostic value, and an overview of the P-values is shown in Figure 6.

AREG, CAV1 and STAG3 were associated with pathways dysregulated in cancer

AREG, CAV1 and STAG3 were associated with module 159 in our study (Figs. 8a, e and i). Module 159 is a translation regulation module in cancer and is defined by mining a large collection of tumor-oriented microarray data. AREG and CAV1 upregulation was related to the KRAS dependency signature in the collection of oncogenic signatures (Figs. 8b and f). In addition, AREG upregulated samples were enriched in the pathways of protein G-2 and S-phase expressed 1 (GTSE1) in G2/M progression after the G2 checkpoint (Fig. 8c), and associated with genes up-regulated in SCC12B2 cells (squamous cell carcinoma) by ultraviolet B (UV-B) irradiation (Fig. 8d). Increased CAV1 levels were associated with the regulation of the guanosine triphosphatase (GTPase) activity of RAS protein, which is stimulated by GTPase-activating proteins (GAPs) (Fig. 8g). In addition, CAV1 upregulation was related to oncogenesis regulated by Met (hepatocyte growth factor receptor) overexpression (Fig. 8h). STAG3 downregulation was linked to the retinoblastoma gene (RB) and P107 down regulation was associated with metastasis and poor differentiation of HNSCC (Figs. 8j-l). However, no enriched pathways met the cutoff criteria when C19orf57 was analyzed.

Ten Small molecules maybe candidate drugs for HNSCC

To further screen potential drugs for HNSCC, the up- and downregulated DEGs were uploaded to run the query. The molecules with the top five positive and negative scores were selected for further study. The details of these molecules were provided in Table 2. The molecules comprised 8 inhibitors, a blocker and an activator. A positive score indicated a similarity between a given perturbagen's signature and that of the query, while a negative score indicated a disparity between the two signatures. In other words, molecules with a positive score may suppress HPV-negative HNSCCs, while molecules with a negative score may suppress HPV-positive HNSCCs. Further studies may confirm the value of these molecules. In addition, the 2D structures of the molecules are shown in Fig. 9A.

NVP-AUY922, fostamatinib, PP-2 inhibited cell proliferation in six HPV-negative HNSCC cell lines

Four candidate small molecule drugs (NVP-AUY922, fostamatinib, PP-2 and zonisamide) were evaluated the therapeutic value in vitro in six HPV-negative HNSCC cell lines. The proliferation of all cell lines was inhibited by each of the four drugs, and the proliferation inhibition rate was negatively associated with the drug concentration (Figs. 9b-d). NVP-AUY922 was the most effective drug, followed by fostamatinib and PP-2. The IC50 values of NVP-AUY922, fostamatinib and PP-2 were range from 0.012 to 0.072 μM, 0.811 to 3.470 μM, 5.32 to 16.41 μM, respectively (Fig. 9e). However, zonisamide may not be therapeutically effective, because its IC50 ranged from 86.46 to 150.5 μM (Additional file 4: Figure S4).

NVP-AUY922 inhibited HPV-negative HNSCC xenografts

To test the therapeutic effect in vivo, the most effective drug (NVP-AUY922) and three cell lines (HN6, SCC9, SCC25) were selected for further research according to the results in vitro (Fig. 9e). The HNSCC xenografts were established in BALB/c nude mice, and treated with 50 mg/kg NVP-AUY922 three times a week for two weeks. As a result, NVP-AUY922 significantly inhibited the growth of all the three HNSCC xenografts in vivo compared with the control group (Fig. 10a-f).

Discussion

In this study, AREG, STAG3, C19orf57 and CAV1 were identified and showed different prognostic values for OS, DFI, DSS and PFI in HPV-positive and -negative patients. And GSEA revealed AREG, CAV1 and STAG3 were associated with dysregulated pathways in cancer. In addition, ten small molecules were screened as potential drugs for HNSCC. Three of them, NVP-AUY922, PP-2 and fostamatinib showed therapeutic value in six HPV-negative HNSCC cell lines, besides, NVP-AUY922 inhibited three HPV-negative HNSCC xenografts in vivo.

AREG is a ligand of the epidermal growth factor receptor (EGFR), AREG binds to EGFR can induce key intracellular signaling cascades controlling cell survival, proliferation and motility [27]. AREG is associated with several tumors, such as lung, breast, colorectal, ovary and prostate carcinomas, due to its role in tumorigenesis [28]. In a clinical trial of patients with HNSCC, patients with higher expression of AREG had shorter OS and PFI than patients with lower expression of AREG [29]. In addition, a recent study indicated vincristine resistance is promoted by AREG in oral squamous carcinoma [30]. Notably, in our study, AREG was a prognostic factor for unfavorable PFI and DSS in HPV-negative patients, but a prognostic factor for favorable PFI in HPV-positive patients. These results indicated that AREG may be strongly associated with the HPV status and is a potential therapeutic target in HNSCC.

STAG3 is a subunit of the cohesin complex that regulates the cohesion of sister chromatids during cell division. STAG3 is also required for chromosome pairing and synapsis, DNA repair and meiosis progression, and mutation of STAG3 may induce DNA repair process abnormalities [31]. A previous study demonstrated that STAG3 may be a tumor suppressor gene, and that loss of STAG3 may cause drug resistance in melanoma [32]. Here, a higher mRNA level of STAG3 was observed in HPV-positive HNSCC samples than in HPV-negative HNSCC samples, and STAG3 was suggested to be a biomarker for poor prognosis in HPV-negative patients. These outcomes were consistent with those of another study [33]. However, because of the small sample size in this study, more studies are needed to clarify the role of STAG3 in HPV-negative HNSCC patients.

Recently, C19orf57 was found to be significantly upregulated in HPV-active oropharyngeal squamous cell carcinoma (OPSCC) patients [34], indicating that C19orf57 may be a key gene associated with HPV. Similarly, our data support C19orf57 is a biomarker for unfavorable PFI times in patients with HPV-positive HNSCC. In addition, CAV1 is a major structural protein in caveolae and reported as an integral membrane protein associated with the progression of carcinoma. However, whether CAV1 act as an oncogene or a tumor suppressor gene in cancer progression is still unclear [35]. Several studies have indicated CAV1 is overexpressed in HNSCC and mediates tumor migration and invasion [36-38]. Similarly, we observed that CAV1 mRNA and protein levels were higher in HNSCC than in normal tissues, suggesting that CAV1 is a factor indicating poor OS. Thus, our results support the function of CAV1 as an oncogene.

Treatment of HNSCC remains challenging because approximately 67% of HNSCC patients present with advanced disease[4], and no effective drugs or suitable surgical methods are available. Thus, ten small molecules were screened as potential drugs in this study. SB-216763, ABT-751, vinorelbine, etacrynic acid and entinostat were suggested as candidate drugs for HPV-positive tumors, while PTB1, zonisamide, NVP-AUY922, PP-2 and fostamatinib were suggested as candidate drugs for HPV-negative tumors. NVP-AUY922, PP-2, and fostamatinib inhibited cell proliferation in vitro in six HPV-negative cell lines, with low IC50 values.

Heat shock proteins (HSPs) are chaperone proteins that can assist the folding of proteins associated with tumor growth. NVP-AUY922 is an HSP inhibitor that suppresses the growth, angiogenesis and metastasis of several kinds of tumors, including oral squamous cell carcinoma [39, 40]. These results support our finding that NVP-AUY922 was the most effective among the investigated drugs. Spleen tyrosine kinase (SYK) controls B-cell receptor (BCR) signal initiation and amplification, which promotes B-cell lymphoma progression [41]. A previous study reported fostamatinib, a SYK inhibitor, suppresses chronic lymphocytic leukemia [42] and prevents metastatic recurrence in breast cancer in vivo [43]. Fostamatinib has been approved for the treatment of immune thrombocytopenia autoimmune hemolytic anemia and IgA nephropathy [44], and the present study indicated that it may be useful in HPV-negative HNSCC patients. SRC is a non-receptor tyrosine kinase family member with a crucial role in tumor progression. PP-2, an inhibitor of SRC, may significantly ameliorate the invasiveness of breast cancer cells and enhance the radiosensitivity of glioma cells [45, 46]. Moreover, in our study, PP-2 showed therapeutic value in six HPV-negative HNSCC cell lines in vitro.

In summary, our study identified four key genes (AREG, STAG3, C19orf57 and CAV1) and ten small molecules for the treatment of HNSCC. Three of the identified small molecule drugs (NVP-AUY922, PP-2 and fostamatinib) showed inhibitory effects on six HPV-negative HNSCC cell lines, and NVP-AUY922 inhibited three HPV-negative HNSCC xenografts in vivo. These genes may be therapeutic targets, and the small molecule drugs may be used alone or in combination with other drugs in further research aiming to improve the prognosis of patients.

Conclusions

In conclusion, our study analyzed DEGs between HPV-positive and -negative HNSCCs. A total of 47 DEGs were screened, and 11 hub genes were identified as hub genes, among which AREG, STAG3, C19orf57 and CAV1 may be considered key prognostic biomarkers and therapeutic targets for HNSCC. In addition, NVP-AUY922, fostamatinib and PP-2 may be used as drugs to treat HPV-negative head and neck cancers.

Abbreviations

BCR: B-cell receptor

BP: Biological process

CC: Cellular component

CCK8: Cell Counting Kit-8

DAVID: Database for Annotation, Visualization and Integrated Discovery

DEG: Differentially expressed genes

DFI: Disease-free interval

DMSO: Dimethyl sulfoxide

DSS: Disease-specific survival

EGFR: Epidermal growth factor receptor

FDR: False discovery rate

GAPs: Guanosine triphosphatase activating proteins

GEO: Gene Expression Omnibus

GO: Gene Ontology

GSEA: Gene set enrichment analysis

GTPase: Guanosine triphosphatase

GTSE1: G-2 and S-phase expressed 1

HNSCC: Head and neck squamous cell carcinoma

HPA: Human Protein Atlas

HPV: Human papillomavirus

HR: Hazard ratio

HSPs: Heat shock proteins

IC50: Half-maximal inhibitory concentrations

KEGG: Kyoto Encyclopedia of Genes and Genomes

MCODE: Molecular Complex Detection

MF: Molecular function

MSigDB: Molecular Signatures Database

OD: Optical density

OPSCC: Oropharyngeal squamous cell carcinoma

OS: Overall survival

PFI: Progression-free interval

PPI: Protein-protein interaction

STRING: Search Tool for the Retrieval of Interacting Genes/Proteins

SYK: Spleen tyrosine kinase

TCGA: The Cancer Genome Atlas

UCSC: University of California Santa Cruz

UV-B: Ultraviolet B

Declarations

Ethics approval and consent to participate

All animal experiments were approved by and performed in accordance with the guidelines of the Shanghai Jiao Tong University School of Medicine.

Consent for publication

Not applicable.

Availability of data and materials

The datasets generated and/or analyzed during the current study are available in GEO (http://www.ncbi.nlm.nih.gov/geo); University of California Santa Cruz Xena (https://xenabrowser.net/); Oncomine (https://www.oncomine.org/); The Human Protein Atlas (http://www.proteinatlas.org) and Connectivity Map (https://clue.io/cmap).

Competing interests

The authors declare that they have no competing interests.

Funding

This work was supported by Shanghai Municipal Key Clinical Specialty (shslczdzk01601).

Authors' contributions

GCT, YF, ZYZ and XY designed the study, GCT, YF and DHZ were the major contributors in writing the manuscript. DHZ, JL, ZYZ and XY helped analyze the data and revise the manuscript with constructive suggestion. All authors read and approved the final manuscript.

Acknowledgments

The authors acknowledge The Ninth People's Hospital affiliated to the Shanghai Jiao Tong University School of Medicine for financial support.

Author details

¹Department of Oral and Maxillofacial-Head Neck Oncology, Shanghai Ninth People's Hospital, College of Stomatology, Shanghai Jiao Tong University School of Medicine, Shanghai, P.R. China

²National Clinical Research Center for Oral Diseases, Shanghai, P.R. China

³Shanghai Key Laboratory of Stomatology & Shanghai Research Institute of Stomatology, Shanghai, P.R. China

⁴Research Unit of Oral and Maxillofacial Regenerative Medicine, Chinese Academy of Medical Sciences, Shanghai, P.R. China

⁵Department of Oral Pathology, Ninth People's Hospital, Shanghai Jiao Tong University School of Medicine, Shanghai, P.R. China

References

1. Bray F, Ferlay J, Soerjomataram I, Siegel RL, Torre LA, Jemal A: Global cancer statistics 2018: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin 2018, 68(6):394-424.

2. Lydiatt WM, Patel SG, O'Sullivan B, Brandwein MS, Ridge JA, Migliacci JC, Loomis AM, Shah JP: Head and Neck cancers-major changes in the American Joint Committee on cancer eighth edition cancer staging manual. CA Cancer J Clin 2017, 67(2):122-137.

3. Zanoni DK, Montero PH, Migliacci JC, Shah JP, Wong RJ, Ganly I, Patel SG: Survival outcomes after treatment of cancer of the oral cavity (1985-2015). Oral oncology 2019, 90:115-121.

4. Cramer JD, Burtness B, Le QT, Ferris RL: The changing therapeutic landscape of head and neck cancer. Nature reviews Clinical oncology 2019, 16(11):669-683.

5. Bratman SV, Bruce JP, O'Sullivan B, Pugh TJ, Xu W, Yip KW, Liu FF: Human Papillomavirus Genotype Association With Survival in Head and Neck Squamous Cell Carcinoma. JAMA Oncol 2016, 2(6):823-826.

6. Kobayashi K, Hisamatsu K, Suzui N, Hara A, Tomita H, Miyazaki T: A Review of HPV-Related Head and Neck Cancer. J Clin Med 2018, 7(9).

7. Barrett T, Wilhite SE, Ledoux P, Evangelista C, Kim IF, Tomashevsky M, Marshall KA, Phillippy KH, Sherman PM, Holko M et al: NCBI GEO: archive for functional genomics data sets--update. Nucleic Acids Res 2013, 41(Database issue):D991-995.

8. Yang K, Zhang S, Zhang D, Tao Q, Zhang T, Liu G, Liu X, Zhao T: Identification of SERPINE1, PLAU and ACTA1 as biomarkers of head and neck squamous cell carcinoma based on integrated bioinformatics analysis. Int J Clin Oncol 2019, 24(9):1030-1041.

9. Jin Y, Yang Y: Identification and analysis of genes associated with head and neck squamous cell carcinoma by integrated bioinformatics methods. Mol Genet Genomic Med 2019, 7(8):e857.

10. Walter V, Yin X, Wilkerson MD, Cabanski CR, Zhao N, Du Y, Ang MK, Hayward MC, Salazar AH, Hoadley KA et al: Molecular subtypes in head and neck cancer exhibit distinct patterns of chromosomal gain and loss of canonical cancer genes. PLoS One 2013, 8(2):e56823.

11. Keck MK, Zuo Z, Khattri A, Stricker TP, Brown CD, Imanguli M, Rieke D, Endhardt K, Fang P, Bragelmann J et al: Integrative analysis of head and neck cancer identifies two biologically distinct HPV and three non-HPV subtypes. Clin Cancer Res 2015, 21(4):870-881.

12. Tomar S, Graves CA, Altomare D, Kowli S, Kassler S, Sutkowski N, Gillespie MB, Creek KE, Pirisi L: Human papillomavirus status and gene expression profiles of oropharyngeal and oral cancers from European American and African American patients. Head Neck 2016, 38 Suppl 1:E694-704.

13. 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. Nature communications 2019, 10(1):1523.

14. Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, Davis AP, Dolinski K, Dwight SS, Eppig JT et al: Gene ontology: tool for the unification of biology. The Gene Ontology Consortium. Nat Genet 2000, 25(1):25-29.

15. Kanehisa M, Sato Y, Furumichi M, Morishima K, Tanabe M: New approach for understanding genome variations in KEGG. Nucleic Acids Res 2019, 47(D1):D590-D595.

16. Huang da W, Sherman BT, Lempicki RA: Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protoc 2009, 4(1):44-57.

17. Szklarczyk D, Gable AL, Lyon D, Junge A, Wyder S, Huerta-Cepas J, Simonovic M, Doncheva NT, Morris JH, Bork P et al: STRING v11: protein-protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets. Nucleic Acids Res 2019, 47(D1):D607-D613.

18. Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, Amin N, Schwikowski B, Ideker T: Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res 2003, 13(11):2498-2504.

19. Bandettini WP, Kellman P, Mancini C, Booker OJ, Vasu S, Leung SW, Wilson JR, Shanbhag SM, Chen MY, Arai AE: MultiContrast Delayed Enhancement (MCODE) improves detection of subendocardial myocardial infarction by late gadolinium enhancement cardiovascular magnetic resonance: a clinical validation study. J Cardiovasc Magn Reson 2012, 14:83.

20. 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 (New York, NY) 2004, 6(1):1-6.

21. Zhai Y, Kuick R, Nan B, Ota I, Weiss SJ, Trimble CL, Fearon ER, Cho KR: Gene expression analysis of preinvasive and invasive cervical squamous cell carcinomas identifies HOXC10 as a key mediator of invasion. Cancer Res 2007, 67(21):10163-10172.

22. Slebos RJ, Yi Y, Ely K, Carter J, Evjen A, Zhang X, Shyr Y, Murphy BM, Cmelak AJ, Burkey BB et al: Gene expression differences associated with human papillomavirus status in head and neck squamous cell carcinoma. Clin Cancer Res 2006, 12(3 Pt 1):701-709.

23. Goldman M, Craft B, Hastie M, Repečka K, McDade F, Kamath A, Banerjee A, Luo Y, Rogers D, Brooks AN et al: The UCSC Xena platform for public and private cancer genomics data visualization and interpretation. bioRxiv 2019:326470.

24. Uhlen M, Fagerberg L, Hallstrom BM, Lindskog C, Oksvold P, Mardinoglu A, Sivertsson A, Kampf C, Sjostedt E, Asplund A et al: Proteomics. Tissue-based map of the human proteome. Science 2015, 347(6220):1260419.

25. Lamb J, Crawford ED, Peck D, Modell JW, Blat IC, Wrobel MJ, Lerner J, Brunet JP, Subramanian A, Ross KN et al: The Connectivity Map: using gene-expression signatures to connect small molecules, genes, and disease. Science 2006, 313(5795):1929-1935.

26. Li Q, Cheng T, Wang Y, Bryant SH: PubChem as a public resource for drug discovery. Drug Discov Today 2010, 15(23-24):1052-1057.

27. Berasain C, Avila MA: Amphiregulin. Semin Cell Dev Biol 2014, 28:31-41.

28. Busser B, Sancey L, Brambilla E, Coll JL, Hurbin A: The multiple roles of amphiregulin in human cancer. Biochim Biophys Acta 2011, 1816(2):119-131.

29. Tinhofer I, Klinghammer K, Weichert W, Knodler M, Stenzinger A, Gauler T, Budach V, Keilholz U: Expression of amphiregulin and EGFRvIII affect outcome of patients with squamous cell carcinoma of the head and neck receiving cetuximab-docetaxel treatment. Clin Cancer Res 2011, 17(15):5197-5204.

30. Hsieh MJ, Chen YH, Lee IN, Huang C, Ku YJ, Chen JC: Secreted amphiregulin promotes vincristine resistance in oral squamous cell carcinoma. Int J Oncol 2019, 55(4):949-959.

31. Hopkins J, Hwang G, Jacob J, Sapp N, Bedigian R, Oka K, Overbeek P, Murray S, Jordan PW: Meiosis-specific cohesin component, Stag3 is essential for maintaining centromere chromatid cohesion, and required for DNA repair and synapsis between homologous chromosomes. PLoS Genet 2014, 10(7):e1004413.

32. Shen CH, Kim SH, Trousil S, Frederick DT, Piris A, Yuan P, Cai L, Gu L, Li M, Lee JH et al: Loss of cohesin complex components STAG2 or STAG3 confers resistance to BRAF inhibition in melanoma. Nat Med 2016, 22(9):1056-1061.

33. Wood O, Woo J, Seumois G, Savelyeva N, McCann KJ, Singh D, Jones T, Peel L, Breen MS, Ward M et al: Gene expression analysis of TIL rich HPV-driven head and neck tumors reveals a distinct B-cell signature when compared to HPV independent tumors. Oncotarget 2016, 7(35):56781-56797.

34. Iman M, Narimani Z, Hamraz I, Ansari E: Network based identification of different mechanisms underlying pathogenesis of human papilloma virus-active and human papilloma virus-negative oropharyngeal squamous cell carcinoma. Journal of the Chinese Chemical Society 2018, 65(11):1307-1316.

35. Fu P, Chen F, Pan Q, Zhao X, Zhao C, Cho WC, Chen H: The different functions and clinical significances of caveolin-1 in human adenocarcinoma and squamous cell carcinoma. OncoTargets and therapy 2017, 10:819-835.

36. Nohata N, Hanazawa T, Kikkawa N, Mutallip M, Fujimura L, Yoshino H, Kawakami K, Chiyomaru T, Enokida H, Nakagawa M et al: Caveolin-1 mediates tumor cell migration and invasion and its regulation by miR-133a in head and neck squamous cell carcinoma. Int J Oncol 2011, 38(1):209-217.

37. Sun J, Gao J, Hu JB, Fan LF, Zhu XB, Subahan R, Chen HL: Expression of Cav-1 in tumour cells, rather than in stromal tissue, may promote cervical squamous cell carcinoma proliferation, and correlates with high-risk HPV infection. Oncol Rep 2012, 27(6):1733-1740.

38. Xue J, Chen H, Diao L, Chen X, Xia D: Expression of caveolin-1 in tongue squamous cell carcinoma by quantum dots. Eur J Histochem 2010, 54(2):e20.

39. Okui T, Shimo T, Fukazawa T, Mohammad Monsur Hassan N, Honami T, Ibaragi S, Takaoka M, Naomoto Y, Sasaki A: Novel HSP90 inhibitor NVP-AUY922 enhances the anti-tumor effect of temsirolimus against oral squamous cell carcinoma. Curr Cancer Drug Targets 2013, 13(3):289-299.

40. Eccles SA, Massey A, Raynaud FI, Sharp SY, Box G, Valenti M, Patterson L, de Haven Brandon A, Gowan S, Boxall F et al: NVP-AUY922: a novel heat shock protein 90 inhibitor active against xenograft tumor growth, angiogenesis, and metastasis. Cancer Res 2008, 68(8):2850-2860.

41. Friedberg JW, Sharman J, Sweetenham J, Johnston PB, Vose JM, Lacasce A, Schaefer-Cutillo J, De Vos S, Sinha R, Leonard JP et al: Inhibition of Syk with fostamatinib disodium has significant clinical activity in non-Hodgkin lymphoma and chronic lymphocytic leukemia. Blood 2010, 115(13):2578-2585.

42. Suljagic M, Longo PG, Bennardo S, Perlas E, Leone G, Laurenti L, Efremov DG: The Syk inhibitor fostamatinib disodium (R788) inhibits tumor growth in the Emu- TCL1 transgenic mouse model of CLL by blocking antigen-dependent B-cell receptor signaling. Blood 2010, 116(23):4894-4905.

43. Shinde A, Hardy SD, Kim D, Akhand SS, Jolly MK, Wang WH, Anderson JC, Khodadadi RB, Brown WS, George JT et al: Spleen Tyrosine Kinase-Mediated Autophagy Is Required for Epithelial-Mesenchymal Plasticity and Metastasis in Breast Cancer. Cancer Res 2019, 79(8):1831-1843.

44. Markham A: Fostamatinib: First Global Approval. Drugs 2018, 78(9):959-963.

45. Eom KY, Cho BJ, Choi EJ, Kim JH, Chie EK, Wu HG, Kim IH, Paek SH, Kim JS, Kim IA: The Effect of Chemoradiotherapy with SRC Tyrosine Kinase Inhibitor, PP2 and Temozolomide on Malignant Glioma Cells In Vitro and In Vivo. Cancer Res Treat 2016, 48(2):687-697.

46. Kwak SM, Seo J, Hwang JT, Sung GJ, Song JH, Jeong JH, Lee SH, Yoon HG, Choi HK, Choi KC: EGFR-c-Src-Mediated HDAC3 Phosphorylation Exacerbates Invasion of Breast Cancer Cells. Cells 2019, 8(8).

Tables

Table 1 Details of three GEO datasets.

Table 2 Ten most significant small molecule drugs.