Introduction

Head and neck squamous cell carcinoma (HNSCC), which arises from the oral cavity, larynx and pharynx, ranks as the sixth most common malignancy with an estimated 835,000 new cases and 43,000 associated deaths worldwide in 2018 [1, 2]. Unfortunately, diagnosis of HNSCC is usually made at advanced stages due to lack of symptoms in the early stage of head and neck tumorigenesis, and the 5-year survival rate is still less than 65% now [3]. It is widely believed that accumulation of numerous genetic alterations in epithelial cells is the essential process driven by the initiation and progression of HNSCC [4]. Therefore, investigation of the potential key biomarkers may help to further uncover the biological basis of HNSCC and improve clinical therapy.

Recently, microarrays based on high-throughput platforms for analysis of gene expression are increasingly valued as a promising and efficient tool to screen significant genetic alternations in carcinogenesis and identify biomarkers for diagnosis and prognosis of cancer [5]. A number of gene expression profiling microarrays have been conducted to find various differentially expressed genes (DEGs) in HNSCC [6]; however, considerable efforts in the identification of biomarker have met with limited success, primarily because of independent numbers of gene profiling. Now, by the means of integrated bioinformatics analysis for available microarray data, it is possible to make more reliable and precise screening results via overlapping relevant datasets.

In the current study, microarray data of gene expression profiles (GSE6631 [7], GSE58911 [8]) and the Cancer Genome Atlas (TCGA) HNSCC data [9] were integrated and analyzed by a series of biological informatics approaches, aberrantly DEGs and pathways were identified in HNSCC, protein–protein interaction (PPI) network was also constructed and hub genes were revealed. Subsequently, we investigated the relationship between the hub genes of subnetwork and overall survival (OS) in TCGA database, and tested the expression status of these hub genes in clinical samples at different stages of tumorigenesis. By this mean, we may bring to light the underlying mechanisms and identify the potential candidate biomarkers for diagnosis and prognosis of HNSCC.

Materials and methods

Microarray data

In the present study, the gene expression profiles (GSE6631, GSE58911) were obtained from Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/). Totally 22 paired samples of HNSCC and normal tissues were consisted in GSE6631, which based on GPL8300 platform (Affymetrix Human Genome U95 chips). GSE58911 dataset was already deposited in GPL6244 (Affymetrix Human Gene 1.0 ST Array), including 15 paired normal and HNSCC samples. Moreover, the TCGA HNSCC data (https://cancergenome.nih.gov/) were also downloaded, including 44 normal and 502 HNSCC tissues. We chose these 3 datasets for integrated analysis to identify commonly DEGs.

Data processing and identification of DEGs

The original raw array data were subjected to background correction, quartile data normalization, and converted into gene expression values. Data were normalized using the Bioconductor R package (https://cran.r-project.org/mirrors.html). Then, the DEGs between HNSCC samples and normal controls were identified using the empirical Bayes approach in linear models for the microarray data (limma) package. |logFC| > 1 and p < 0.05 were selected as the cutoff criterion.

Functional and pathway enrichment analysis of DEGs

To analyze the identified DEGs at the functional level, the significant gene ontology (GO) biological process terms [10] and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis [11] were performed using the Database for Annotation, Visualization and Integrated Discovery (DAVID, https://david.ncifcrf.gov/) with the thresholds of p < 0.05 and false discovery rate (FDR) < 0.01 [12].

Modules from the PPI network

To evaluate the interactive relationships among DEGs, we mapped the DEGs to the Search Tool for the Retrieval of Interacting Genes (STRING) database (https://string-db.org) [13]. Then, the interactive DEGs were selected to construct the PPI network (combined score ≥ 0.4) and visualized using Cytoscape [14]. The Molecular Complex Detection (MCODE) plugin in Cytoscape was used to screen the modules of PPI network with MCODE score > 3 and number of nodes > 4.

Survival analysis of the hub gene in TCGA database

The association between the corresponding genes in the top modules and patient OS for 5 years was analyzed using HNSCC samples from the TCGA data. All HNSCC patients were classified into high or low expression based on whether Z-score expression was > median (high) or < median (low), log-rank analysis and Kaplan–Meier plots were produced using Bioconductor R package.

Clinical samples and clinical staging system

A total of 52 paraffin-embedded HNSCC (39) and oral epithelial dysplasia (OED, 13) samples were obtained from the archives of the Department of Pathology of Shandong Provincial Hospital, Jinan, China. Of the HNSCC samples, there were 23 (59%) well, 10 (26%) moderately and 6 (15%) poorly differentiated HNSCC tissues; 10 matched adjacent non-cancerous oral mucosa (NOM) tissues were selected from the above-mentioned patients and detailed sample information is presented in Supplementary Table 1. For the use of these clinical materials for research purposes, prior patient consent and approval from the Institute Research Ethics Committee were obtained. The approval number was no. 2018-037. The stages of all HNSCC patients were classified according to the Union for International Cancer Control (UICC 2017).

Immunohistochemical analysis

Protein expression was evaluated on paraffin-embedded sections using microwave antigen retrieval with 0.01 M citrate buffer (pH 6.0). Rabbit polyclonal anti-SERPINE1 (1:200, sc-5297, Santa Cruz), rabbit polyclonal anti-PLAU (1:150, ab133563, Abcam) and rabbit polyclonal anti-ACTA1 (1:100, sc-5867, Santa Cruz) antibodies were utilized. The immunohistochemical procedure was as previously described [15, 16]. Interpretation of immunohistochemical staining was made independently by two specialists and the mean protein staining intensity (SI; 0, no staining; 1, weak; 2, moderate; 3, strong), labeling indices (LIs, defined as the percentage of positive cells in total cells), and mean labeling scores (LS, defined as LI × SI) in HNSCC, OED and NOM samples were calculated and compared among groups.

Statistical analysis

Comparison of two Kaplan–Meier curves was performed using the log-rank test of the R-package survival. The mean LS for HNSCC, OED and NOM samples was compared among the three groups by analysis of variance (ANOVA) using the SPSS 10.0 software package. p < 0.05 was considered statistically significant.

Results

Identification of aberrantly DEGs in HNSCC

Data from each microarray were separately analyzed to screen DEGs. As represented in Fig. 1, our integrated bioinformatics analysis indicated that a total of 132 genes were consistently and significantly deregulated in the same direction in these datasets, including 52 overlapping up-regulated (163 in GSE6631, 172 in GSE58911 and 5679 in TCGA) and 80 down-regulated genes (128 in GSE6631, 427 in GSE58911 and 3436 in TCGA) in HNSCC tissues, as compared to normal epithelial tissues (Table 1).

Fig. 1
figure 1

Identification of 132 common DEGs from the three cohort profile datasets (GSE6631, GSE58911 and TCGA) using Morpheus Website (https://software.broadinstitute.org). Different color areas represented different datasets. The cross areas meant the commonly changed DEGs. DEGs were identified with classical t test; statistically significant DEGs were defined with p < 0.05 and |logFC| > 1 as the cutoff criterion

Full size image
Table 1. 132 DEGs were identified from three profile datasets, including 52 up-regulated genes and 80 down-regulated genes in the HNSCC tissues, compared to normal tissue
Full size table

DEGs functional and pathway enrichment analysis

The top 5 significant terms of GO analysis in DAVID were illustrated in Table 2. In the biological process (BP) group, GO analysis results showed that up-regulated DEGs were significantly enriched in extracellular matrix organization, collagen catabolic process, extracellular matrix disassembly, cell adhesion and collagen fibril organization; the down-regulated DEGs were mainly enriched in muscle contraction and muscle filament sliding. For molecular function (MF), the enrichment of up-regulated DEGs was in metalloendopeptidase activity, extracellular matrix structural constituent, serine-type endopeptidase activity, collagen binding and endopeptidase activity, and down-regulated genes were enriched in structural constituent of muscle. Besides, GO cell component (CC) analysis indicated that the enrichment of up-regulated DEGs was predominantly in extracellular matrix, proteinaceous extracellular matrix, extracellular region, extracellular space and basement membrane, and down-regulated DEGs were enriched in extracellular exosome, muscle myosin complex and Z disc.

Table 2 Go analysis of DEGs associated with HNSCC
Full size table

We also determined the canonical signaling pathways associated with the commonly DEGs in the carcinogenesis of HNSCC by performing KEGG analysis. The activated pathways were enriched in ECM–receptor interaction, focal adhesion, PI3K–Akt signaling pathway, while the suppressed pathways were mainly involved in drug metabolism–cytochrome P450, tyrosine metabolism and tight junction (Table 3).

Table 3 KEGG pathway enrichment analysis of up-regulated and down-regulated DEGs
Full size table

PPI network construction and module analysis

Using the STRING online database and Cytoscape software, 90 DEGs (37 up-regulated and 53 down-regulated genes) of the 132 commonly DEGs were filtered into the PPI network complex, containing 131 nodes and 289 edges (Fig. 2), and the top 10 node degree genes were FN1, MMP9, SPP1, COL3A1, MMP13, POSTN, SPARC, COL4A1, ACTA1 and SERPINE1. According to the importance of degree, we chose two most significant modules from the PPI network complex for further analysis using Cytoscape MCODE. Pathway enrichment analysis showed that the module 1 consisted of 12 nodes and 40 edges, which were mainly associated with ECM–receptor interaction, focal adhesion and PI3K–Akt signaling pathway, while the module 2 included 12 nodes and 31 edges, which were also enriched in focal adhesion and ECM–receptor interaction (Fig. 3), suggesting that ECM–receptor interaction and focal adhesion signaling pathways were essential in the carcinogenesis of HNSCC.

Fig. 2
figure 2

DEGs PPI network complex and module analysis. Using the STRING online database, a total of 90 DEGs (37 up-regulated in red standing and 53 down-regulated genes in green standing) were filtered into the DEGs PPI network complex

Full size image
Fig. 3
figure 3

Top two modules from PPI network. a Module 1 and its enriched pathways (upper); b module 2 and its enriched pathways (bottom)

Full size image

The validation of hub genes as independent predictors for OS in the TCGA database

We subsequently sought to assess the significance of hub genes for HNSCC; the relationships between expression of hub genes and OS were verified in the TCGA HNSCC cohort. For most of the hub genes, our results showed that poor OS was associated only in those patients with high expression of SERPINE1 (p = 0.00054) or PLAU (p = 0.00289), as well as the low expression of ACTA1 (p = 0.04147), MYL1 (p = 0.01405), MYH2 (p = 0.04987) or MYLPF (p = 0.02122) (Fig. 4), suggesting that these candidate genes are associated with clinical outcome of HNSCC patients.

Fig. 4
figure 4

Kaplan–Meier curves exhibit the OS in the TCGA HNSCC cohort with high and low expression of SERPINE1, PLAU, ACTA1, MYL1, MYH2 and MYLPF

Full size image

SERPINE1, PLAU and ACTA1 are aberrantly expressed in the carcinogenesis of HNSCC.

To further clarify the potential biological role of these prognosis-associated genes in HNSCC transformation, we next characterized the expression changes of signature genes by microarray analysis (Supplementary Fig. 1). Among these hub genes that had large change levels in HNSCC samples, we compared their degrees in the highest ranked modules, two up-(SERPINE1, PLAU) and one down-regulated (ACTA1) genes were particularly selected to further test the protein expression in NOM, OED and HNSCC tissues. As expected, we found that the expression of SERPINE1 and PLAU increased from NOM to OED and HNSCC. The mean LS of SERPINE1 increased in OED (230.38 ± 53.046%) and HNSCC samples (205.85 ± 51.427%), compared with that of NOM (58.50 ± 15.072%). Likewise, the level of PLAU also increased significantly from NOM (51.80 ± 10.727%) through OED (194.54 ± 46.321%) to HNSCC samples (199.61 ± 47.895%), there was a significant difference in the mean LS of SERPINE1 or PLAU between NOM and OED (p = 0.000) or HNSCC samples (p = 0.000), respectively; however, the LS of SERPINE1 or PLAU between OED and HNSCC had no significant difference. Representative microphotographs of SERPINE1 and PLAU staining for NOM, OED and HNSCC are shown in Fig. 5. Instead, the level of ACTA1 showed an opposite trend; the LS was reduced from NOM (136.10 ± 49.249%) through OED (81.77 ± 5.403%) to HNSCC samples (60.66 ± 9.089%). There was a significant difference in the expression of ACTA1 among the 3 groups (Table 4). Thus, combined with the TCGA data analysis, these results suggested that SERPINE1, PLAU and ACTA1 are required for the initiation of head and neck tumorigenesis.

Fig. 5
figure 5

Immunohistochemical staining of SERPINE1, PLAU and ACTA1 in oral epithelium at different stages of head and neck carcinogenesis (× 40). The SERPINE1 was expressed in nucleus and cytoplasm, while PLAU was mainly expressed in cytoplasm; the cytoplasmic and/or nuclear staining intensity of basal layers were much denser than that of upper epithelial cells for SERPINE1 and PLAU. The expression of ACTA1 was mainly distributed at the cytoplasm of epithelial cells in prickle layer (red arrow, × 200). NOM, non-cancerous oral mucosa; OED, oral epithelial dysplasia

Full size image
Table 4 The mean LS of SERPINE1, PLAU and ACTA1 in NOM, OED and HNSCC tissues
Full size table

SERPINE1, PLAU and ACTA1 are correlated with clinical aggressiveness of HNSCC patients

As the expression levels of SERPINE1, PLAU and ACTA1 were validated as independently predicted factors for OS of HNSCC patients, we continued to define the association of these genes and clinical histology classification in HNSCC samples. As shown in Fig. 6, there was a significant difference in these three hub genes expression between well and moderately or poorly differentiated HNSCC. The results showed that the LS of SERPINE1 increased significantly from well (166.78 ± 13.426%) through moderately (234.60 ± 36.439%) to poorly differentiated HNSCC samples (282.00 ± 7.589%) (Table 5). A significant difference in the LS of PLAU was also found between poorly and well (p = 0.000) or moderately differentiated HNSCC tissues (p = 0.005). In contrast, the expression of ACTA1 showed an obviously downward trend between well and moderately or poorly HNSCC samples; the mean LS of ACTA1 in well-differentiated HNSCC was 64.78 ± 9.400%; however, the results of LS between moderately and poorly differentiated HNSCC were all zero. Thus, our findings indicated that SERPINE1, PLAU and ACTA1 were correlated with clinical malignancy of HNSCC patients.

Fig. 6
figure 6

Immunohistochemical staining of SERPINE1, PLAU and ACTA1 in tumor nests of HNSCC tissues, and the expression levels of three genes were associated with tumor differentiation (× 100)

Full size image
Table 5 The mean LS of SERPINE1, PLAU and ACTA1 in well, moderately and poorly differentiation of HNSCC tissues
Full size table

Discussion

Identifying oncogenic biomarkers and elucidating the underlying mechanism of the initiation and development of HNSCC would greatly benefit the early diagnosis and effective treatment for patients with high malignancy [17]. Emergency bioinformatics analysis has provided a powerful tool for the identification of biomarkers and therapeutic targets relevant to tumor progression and treatment response [18]. In the present study, we identified 52 up-regulated and 80 down-regulated DEGs through analyzing available data of gene expression profile datasets (GSE6631, GSE58911 and TCGA) in HNSCC by multiple bioinformatics tools. Functional analysis demonstrated that these DEGs are mainly associated with activation of ECM–receptor interaction and focal adhesion and suppression of drug metabolism–cytochrome P450 pathways. More importantly, based on TCGA dataset, our clinical experiments proved that a set of prognostic signatures including SERPINE1, PLAU and ACTA1 were identified as biomarkers for diagnosis and prognosis of HNSCC, which may provide novel insights for unraveling pathogenesis of HNSCC.

Recently, some basic studies have been conducted to identify the DEGs in HNSCC [19, 20]. For example, Yang et al. analyzed the gene expression profile of GSE6791 and identified 550 up-regulated and 261 down-regulated genes [21]. Similarly, Zhao found that PLAU, CLDN8 and CDKN2A could predict OS using gene expression profiles of GSE13601, GSE30784, GSE37991 and TCGA in oral squamous cell carcinoma [22]. Our integrated bioinformatics analysis indicated that 132 genes were consistently and significantly deregulated in GSE6631, GSE58911 and TCGA. Interestingly, our results revealed that there were also examples of genes that did not overlap compared with these reports; the main reason of this discrepancy may be because we used 3 different multiple profiles, which could greatly minimize the intra-tumoral heterogeneity and diversity of anatomical sites of HNSCC.

As was suggested by DAVID analysis, the up-regulated DEGs were mainly involved in extracellular matrix organization, collagen catabolic process, extracellular matrix disassembly, cell adhesion and collagen fibril organization at the level of BP. Extracellular matrix (ECM), as a crucial component of the cancer cell niche, provides the mechanical support for the tissue and mediates the cell–microenvironment interactions [23]. Significantly, collagens are one of the major proteins found within the ECM, and have themselves been implicated in many aspects of neoplastic transformation. Therefore, it is consistent with the findings that active functions of these cellular processes through ECM were the main cause for tumor development, progression and metastasis [24], whereas the down-regulated DEGs in HNSCC were mainly enriched in actin-mediated cell contraction and filament sliding, which were associated with decreased muscle function-mediated cytoskeleton remodeling in cancer development and progression [25]. Furthermore, the enriched KEGG pathway of up-regulated DEGs mainly induced ECM–receptor interaction, focal adhesion and PI3K–Akt signaling pathway. Significantly, 12 overlapping genes, including ITGA6, SPP1 and FN1, were identified to functionally involve in interactions between ECM and cells by activating these three signaling pathways, which lead to a direct or indirect control of cellular activities such as cell migration, differentiation and proliferation [26,26,28]. As a contrast, down-regulated DEGs were related to drug metabolism–cytochrome P450. The recent study has reported that the cytochrome P450 slowed metabolizers CYP2C9*2 and CYP2C9*3, which could directly regulate tumorigenesis via reduced epoxyeicosatrienoic acid production [29]. Together, these data suggested that deregulated pathways may be a major factor of HNSCC tumorigenesis, detecting these aberrant signaling pathways could precisely predict tumor progression [30].

After constructing PPI network with DEGs and listing the top degree of hub genes, the most significant two modules were filtered from the PPI network complex; consequent functional analysis showed that most of corresponding genes were associated with ECM–receptor interaction and focal adhesion. Furthermore, survival analysis of hub genes in these two modules revealed that SERPINE1, PLAU, ACTA1, MYL1, MYH2 and MYLPF were identified as prognostic markers for clinical outcome in the TCGA cohort. Among the up-regulated hub genes, PLAU, one of the major proteolytic enzymes involved in degradation of extracellular matrix, has been demonstrated to play critical roles in tissue remodeling and migration in the developmental as well as tumorigenesis process, whereas SERPINE1, as the most important physiological inhibitor of the PLAU, could in turn reverse this process and regulate the adhesion/ deadhesion balance of cells to the ECM [31]. However, it has been reported that SERPINE1 could induce the EMT process and promote tumor cell survival in breast and ovarian cancers [32, 33]. In our study, the bioinformatics analysis revealed significantly increased expressions of PLAU and SERPINE1 in HNSCC tissues, which were associated with poor clinical outcome. In contrast, for the down-regulated actin-family genes, ACTA1 gene encodes a protein exerting functions in cell motility, structure and integrity. Consistent with our observation, ACTA1 is also down-regulated in colorectal cancer [34]. In addition, our results showed that the other three specific down-regulation genes (MYL1, MYH2 and MYLPF) were involved in muscle contraction process, which might play a regulatory role in remodeling of muscle function in HNSCC tissues; however, the specific roles of these genes in cancers still need to be elucidated.

Of note, in view of the prognostic potency of these hub genes for HNSCC in TCGA database, by the validation of their top degree of genes and change levels of mRNA in microarrays, we selected SERPINE1, PLAU and ACTA1 to further detect their protein level by immunostaining. Our clinical analysis showed that SERPINE1, PLAU and ACTA1 were significantly changed in the progression of HNSCC. They were aberrantly expressed in the epithelium of OED and HNSCC and correlated with aggressiveness of HNSCC patients, which implied that these signature genes are possibly not only involved in the initiation of tumorigenesis but also late stages of cancer. Therefore, SERPINE1, PLAU and ACTA1 could be potentially utilized as diagnostic and prognostic biomarkers for HNSCC. More importantly, by comparing the extent of protein changes, the overexpressed SERPINE1 and PLAU are the most promising markers, and its detection could help to identify tumor cells in tissues.

In conclusion, the current study was intended to identify DEGs with comprehensive bioinformatics analysis to find the potential biomarkers and predict progression of HNSCC. We found that SERPINE1, PLAU and ACTA1 might be exploited as diagnostic and prognostic indicators for HNSCC. Finally, the results also suggested that the function of ECM–receptor interaction and focal adhesion may be essential signaling pathways in the development of HNSCC. Hence, our findings could significantly improve our understanding of the cause and underling molecular events of HNSCC, and provide potential targets for anticancer therapies.