Head and neck squamous cell carcinoma (HNSCC) is an aggressive malignancy with high morbidity and mortality rates and ranks as the sixth most common cancer all over the world. Despite numerous advancements in therapeutic methods, the prognosis of HNSCC patients still remains poor. Therefore, there is an urgent need to have a better understanding of the molecular mechanisms underlying HNSCC progression and to identify essential genes that could serve as effective biomarkers and potential treatment targets. In the present study, original data of three independent datasets were downloaded from the Gene Expression Omnibus database (GEO) and R language was applied to screen out the differentially expressed genes (DEGs). PYGM and TNNC2 were finally selected from the overlapping DEGs of three datasets for further analyses. Transcriptional and survival data related to PYGM and TNNC2 was detected through multiple online databases such as Oncomine, Gene Expression Profiling Interactive Analysis (GEPIA), cBioportal, and UALCAN. Quantitative real-time polymerase chain reaction (qPCR) analysis was adopted for the validation of PYGM and TNNC2 mRNA level in HNSCC tissues and cell lines. Survival curves were plotted to evaluate the association of these two genes with HNSCC prognosis. It was demonstrated that PYGM and TNNC2 were significantly down-regulated in HNSCC and the aberrant expression of PYGM and TNNC2 were correlated with HNSCC prognosis, implying the potential of exploiting them as therapeutic targets for HNSCC treatment or potential biomarkers for diagnosis and prognosis.

Head and neck squamous cell carcinoma (HNSCC) ranks as the sixth most common cancer all over the world and approximately 600000 new cases occur every year [1,2]. Despite the great efforts into the therapeutic advancements, the overall survival (OS) rate of HNSCC has not markedly improved and the overall prognosis still remains poor [3]. One of the possible reasons for the low survival rate and bad prognosis may be delayed diagnosis. It was reported that a majority of HNSCC patients were diagnosed at an advanced stage due to the asymptomatic nature of early stage disease and lack of effective screening modalities [4]. Therefore, reliable and effective means for the early detection of HNSCC and innovative treatment methods for this malignant tumor should be investigated.

For many years, microarrays based on high-throughput platforms served as invaluable and efficient tools for biomedical research by screening significant genetic alterations. Numerous studies have applied microarrays on different cancers to explore the ectopic gene expression to find out essential genes that may play essential roles in the pathological and prognosis of these cancers. However, due to the relatively high expense of microarrays and the limited sample size, the findings sometimes may not be comprehensive and could not be put into clinical application. Since it has increasingly been recognized that through the combinations of the results from different microarrays based on diverse sample origins, the statistical power is increased, the predictive power is more accurate, and the bias of individual studies can be overcome, we decided to make use of a variety of public database to detect relative reliable target genes for HNSCC [5].

The Cancer Genome Atlas (TCGA) and Gene Expression Omnibus (GEO) are two of the most extensive databases where large amounts of data generated from microarrays or high-sequencing technique are collected [6,7]. The original datasets can be retrieved from the database and therefore researchers could select the specific data that meets their study purpose. More importantly, there emerged a variety of online database, such as Oncomine, Gene Expression Profiling Interactive Analysis (GEPIA), UALCAN and c-BioPortal in which researchers have free access to the wide range of datasets and can even perform analysis directly. Therefore, dependent on the wide range of these databases, new studies that made use of the publicly archived gene expression data in various ways and the reuse of these public data can be very powerful [8]. In particular, the rational combination of the datasets has the potential to predict novel gene biomarkers and investigate important genes involved in disease progression so as to develop precision therapies [9]. For example, a study found several significant genes related with ovarian cancer progression and drug resistance by retrieving the data from TCGA, GEO and Oncomine [10]. Also, immune checkpoint HHLA2 was observed to be up-regulated in clear cell renal cell carcinoma and HHLA2 overexpression leads to a remarkable shorter OS and poorer prognosis, implying it could be a potential prognostic biomarker [11].

So far, a large number of genes have been illustrated to be associated with the process and prognosis of HNSCC. The emerging molecular signatures pioneer a promising field where ectopic gene expression may serve as effective biomarker for HNSCC initiation and treatment methods focused on gene regulation may improve the existing therapeutic strategies. Previous studies have suggested some growth factors, tumor suppressor proteins, cell cycle proteins and cytoskeletal proteins could function as potential molecular biomarkers for HNSCC [12]. In the present study, we utilized bioinformatics approaches to identify differentially expressed genes (DEGs) between HNSCC tissues and normal controls from three independent GEO datasets. Then, we selected PYGM and TNCC2, both were not investigated in HNSCC before, to evaluate their relative expression and the potential roles in HNSCC by several online databases. The study aimed to discover novel biomarkers and potential gene therapy targets for HNSCC.

GEO datasets and data processing

The GEO (http://www.ncbi.nlm.nih.gov/geo/) is a public repository for high-throughput microarray and next-generation sequence functional genomic datasets and the raw data, processed data and metadata can all be obtained from the resources [6]. In this study, we downloaded the original CEL data from three datasets (GSE13601, GSE31056, GSE30784) which performed microarray on HNSCC samples and matched normal tissues to identify DEGs. R language was applied to identify the DEGs and the main procedures were illustrated concisely as follows. We utilized the affy package to perform the background correction and data normalization, including conversion of raw data formats, imputation of missing values and background correction. Then, the samples were subjected to differential expression analysis by using the Limma package. P<0.01 and |log(FC)| > 2 were set as the threshold and the genes that met the criteria were screened out as DGEs. The intersection of DEGs from three datasets were performed by the VennDiagram package in R language at last.

Gene Ontology and KEGG pathway analysis

The Database for Annotation, Visualization and Integrated Discovery (DAVID, https://david.ncifcrf.gov/) is an online biological information database for gene functional analysis [13]. Therefore, to predict the possible functions of the overlapping DEGs and their neighboring genes and the potential pathways they may participate in, these genes underwent Gene Ontology (GO) functional and KEGG pathway enrichment analysis in DAVID. P<0.05 was considered as statistically significant.

Oncomine analysis

Oncomine (https://www.oncomine.org/), an online platform providing cancer microarray datasets and data-mining functions, was aimed at facilitating the discovery of significant genes whose expression was closely related with tumor development and progression [14]. Apart from presenting differential expression analyses comparing a variety of cancers with respective normal tissues, it also provides gene expression in different cancer subtypes and an overview of specific gene expression in a wide range of major cancers. Therefore, we explored the mRNA expression of PYGM and TNNC2 in various types of cancer and examined their levels in different HNSCC datasets by using Oncomine database. P-value <0.05, fold change of 1.5, and gene ranking of All were defined as the threshold and all the results are summarized in Figure 5.

GEPIA

GEPIA (http://gepia.cancer-pku.cn/), a web-based tool based on TCGA and GTEx data, could provide fast and interactive functionalities including differential expression analysis, profiling plotting, correlation analysis, patient survival analysis, similar gene detection and dimensionality reduction analysis [15]. In the current study, we overviewed the relative expression of PYGM and TNNC2 in a good number of cancers and |log2FC| = 1.5 and P-value =0.05 were set as the cutoff. UALCAN analysis, UALCAN (http://ualcan.path.uab.edu/), an easy to use, user-friendly and interactive web-portal for in-depth analyses of gene expression data on cancers from TCGA [16], could be utilized to analyze the relative expression of target genes across tumor and normal samples and relative clinicopathologic parameters. Therefore, we assessed the relationship of PYGM and TNCC2 expression with the grade and stage of the HNSCC samples by UALCAN analysis.

UCSC Xena browser

UCSC Xena (http://xenabrowser.net/) is an online database that allows users to explore functional genomic datasets for correlations between genomic and/or phenotypic variables. In the present study, we utilized this free online tool to detect whether the level of PYGM and TNNC2 was correlated with the survival of HNSCC patients from TCGA samples. Patients were grouped into relative high expression group and low- or medium-expression group accordingly and P<0.05 was considered statistically significant.

TCGA and c-BioPortal databases

TCGA is a comprehensive and coordinated project designed to discover essential cancer-causing genomic alterations. Based on these information and datasets, it is likely to improve diagnosis methods, treatment standards and prevent cancer ultimately. So far, TCGA has collected large-scale genome sequencing and pathological data analysis of over 30 human tumors, suggesting that it serves as a well-rounded and reliable cancer database [7]. c-BioPortal (www.cbioportal.org) is an online open-access resource where researches have access to explore, visualize and analyze multidimensional cancer genomics data [17]. In this study, c-BioPortal was applied to access HNSCC (TCGA, Provisional) data. The genomic profiles included mutations, putative copy-number alterations from GISTIC, mRNA expression Z-scores and protein expression Z-scores. OS and disease-free survival (DFS) plotter were drawn according to the instructions on c-BioPortal. Furthermore, to achieve a comprehensive understanding of the DEGs selected in the present study, the gene network of DEGs and the neighboring genes was generated for the analysis of the interaction between these genes.

Tissue specimens

A total of 60 pairs of HNSCC tissues and adjacent non-tumor tissues were collected from patients who received no chemotherapy or radiotherapy between October 2016 and September 2018 at the Department of Oral Maxillofacial-Head and Neck Oncology, Ninth People’s Hospital, Shanghai Jiao Tong University School of Medicine (Shanghai, China). The protocol of the study was approved by the Review Board of the Medical Ethics Committee of the Ninth People’s Hospital, Shanghai Jiao Tong University School of Medicine and all patients provided written informed consent prior to the present study. All the tissue samples were pathologically confirmed as HNSCC and tumor pathological differentiation and clinical stage were defined according to World Health Organization Classification of Tumors and the TNM classification system of the International Union Against Cancer (1988), respectively. OS was defined as the time from start of treatment until death. Follow-up ended in December 2017 or at death.

Cell culture

SCC-9, SCC-25 and CAL 27 cells were purchased from the American Type Culture Collection (ATCC, U.S.A.), and the human HNSCC cell lines HN4, HN6 and HN30 were kindly provided by the University of Maryland Dental School, U.S.A. SCC-4, SCC-9 and SCC-25 cells were cultured in Dulbecco’s modified Eagle’s medium/nutrient mixture F-12 (GIBCO-BRL, U.S.A.) supplemented with 10% fetal bovine serum (GIBCO-BRL, U.S.A.) and 1% penicillin and streptomycin at 37°C, 5% CO2 in a humidified atmosphere while other cells were maintained in DMEM with the same additives. Furthermore, normal oral epithelial cells which were obtained from primary culture were cultured in keratinocyte serum-free medium (GIBCO-BRL, U.S.A.) containing 0.2 ng/ml recombinant epidermal growth factor (Invitrogen, U.S.A.).

RNA extraction and quantitative real-time polymerase chain reaction

Total RNA of cells and tissues was extracted using TRIzol according to the manufacturer’s protocol (Takara, Japan). Afterward, total RNA (1 μg) was reverse transcribed into cDNA using a PrimeScript™ RT reagent kit (Takara, Japan). The quantitative real-time polymerase chain reaction (qPCR) experiment was carried out on an ABI StepOne real-time PCR system (Life Technologies, U.S.A.) by using a SYBR Premix ExTaq Reagent Kit (Takara, Japan). The detailed reaction conditions were as follows: 95°C for 5 min, 40 cycles of 5 s at 95°C and 30 s at 60°C. All the primers were designed and synthesized by Sangon Biotech (Shanghai) and the detailed sequences for the primers are as follows. PYGM (forward primer: 5′-CCATGCCCTACGATACGCC-3′ and reverse primer: 5′-TAGCCACCGACATTGAAGTCC-3′), TNNC2 (forward primer: 5′-TGGGGACATCAGCGTCAAG-3′ and reverse primer: 5′-CCAAGAACTCCTCGAAGTCGAT-3′), and GAPDH (forward primer: 5′-ACAACTTTGGTATCGTGGAAGG-3′ and reverse primer: 5′-GCCATCACGCCACAGTTTC-3′). GAPDH was used as the internal reference in the present study. The relative expression of target genes was normalized to GAPDH and 2−ΔΔCt method was applied to evaluate relative mRNA levels. All the experiments were repeated for three times.

Statistical analyses

Statistical analysis in the present study was performed using SPSS 21.0 statistical software package. All numerical data were expressed as mean ± standard deviation (SD) from triplicate experiments and the significant differences between two groups were calculated by Student’s two-tailed t test. Kaplan–Meier method was applied to generate survival curves with a log-rank test. The associations between gene expression and clinicopathological characteristics of HNSCC patients were analyzed by Mann–Whitney U-test and Kruskal–Wallis test. P<0.05 was defined as statistically significant.

Identification of overlapping DEGs in the GEO datasets

A total of three microarray datasets (GSE31056, GSE30784, GSE13601) were obtained from GEO to identify the DEGs between HNSCC tissues and matched normal tissues. The detailed information about GEO datasets are summarized in Table 1. |Log(FC)|>2 and adj.P-value <0.01 were defined as threshold. Volcano plots were generated by R language to visualize the distribution of expressed genes between cancer and normal controls from three independent studies (Figure 1A–C). Significantly up-regulated or down-regulated genes were represented by red or green dots, respectively. As suggested by Figure 2A, analysis results of R language manifested 11 genes commonly differentially expressed in these three datasets (PTHLH, POSTN, SPP1, KRT17, PLAU, MMP13, ISG15, STAT1, TNC, PYGM, TNNC2). Afterward, we utilized the heatmap2 package in R to generate heatmaps based on the expression levels of common DEGs in GSE31056, GSE30784 and GSE13601 (Figure 2B–D). Each column represented an individual sample and each row represented a specific gene. The color and the intensity indicated the relative expression levels of these DEGs between cancer tissues and normal controls.

The network of DEGs and their neighboring genes in HNSCC

To shed light on the potential mechanism of DEGs involved in HNSCC, a gene regulation network containing DEGs and the 50 most frequently altered neighboring genes was constructed in c-Bioportal. As illustrated by Figure 3, there existed two major modules in the network where STAT1, PYGM and TNNC2 from DEGs played significant roles. A great number of genes were in relationship with these three genes, suggesting the comprehensive and complex functions they played in the biological process of HNSCC. Furthermore, part of the 50 most frequently altered neighboring genes such as EGFR, PIK3CA and TP53, are well established as essential genes in cancer development and progression. No matter whether DEGs were associated with these important cancer-related genes directly or not, the network implied the possible regulatory role of DEGs in HNSCC. What is more, it was shown that PYGM was related with flavopiridol, an antineoplastic agent by inhibiting most cyclin-dependent kinases (cdks). Previous study has demonstrated that flavopiridol could potently inhibit cell proliferation and induce apoptosis in HNSCC cells [18], suggesting that regulation of PYGM may serve as a promising anti-cancer strategy in HNSCC.

GO analysis and KEGG analysis of the DEGs and neighboring genes

For a more in-depth understanding of the DEGs, GO analyses and KEGG pathway enrichment analyses were performed on the DEGs and the 50 most frequently altered neighboring genes. As illustrated in Figure 4A, for biological process, these genes were mostly associated with homeostatic process, response to virus, hypoxia or organic substance, protein amino acid phosphorylation, regulation of apoptosis and cell death and intracellular signaling cascade. It can be easily concluded that the biological process these genes participated in were closely correlated with cancer development and progression. Regarding cellular component, GO analysis results showed that the DEGs and their neighboring genes were mainly enriched in extracellular region, sarcomere, myofibril and contractile fiber (Figure 4B). For molecular function classification, the genes were significantly enriched in the following functions: calcium ion binding, cytokine and channel activity, transmembrane transporter activity, actin and ATP binding and peptidase activity (Figure 4C). The results from KEGG analysis showed that the among the pathways these genes particularly enriched, many were cancer-related pathways such as Toll-like receptor signaling pathway, JAK-STAT signaling pathway, focal adhesion and MAPK signaling pathway (Figure 4D). The above results indicated that these genes may modulate HNSCC proliferation and metastasis through multiple signaling pathways.

Overview of PYGM and TNNC2 expression in various cancers and different HNSCC datasets

Since we have screened out eleven DEGs from three GEO datasets, we then searched them in the Pubmed to see if they were explored by previous studies. Results manifested that POSTN, SPP1, KRT17, ISG15, STAT1 and TNC were all identified as potential biomarkers in HNSCC by former studies and the relative expression in HNSCC were also further validated by qPCR or IHC experiment [19–23]. Moreover, PTHLH, PLAU and MMP13 were demonstrated to be associated with the survival and prognosis of HNSCC patients and some biological functions of these genes were also validated in in vitro experiments [24–26]. However, there are no researches about the relationship of PYGM and TNNC2 with HNSCC and we decided to determine the level of PYGM and TNNC2 in HNSCC and to investigate whether they can act as effective biomarkers for HNSCC. Therefore, we overviewed the transcriptional expression differences of these two genes between cancer tissues and normal tissues in several cancer types and different HNSCC datasets in Oncomine and GEPIA. As shown in Figure 5A–D, mRNA levels of PYGM and TNNC2 were significantly down-regulated in a majority of cancer types and a multiple of independent analyses on HNSCC illustrated the down-regulation of PYGM and TNNC2 in tumor samples. The detailed data about several independent HNSCC datasets in Oncomine are listed in Table 2. From our perspective, the overall underexpression of PYGM and TNNC2 in various cancers may imply the role of them in cancer development and progression and further investigation was of great value.

Relationship of PYGM and TNNC2 expression with the clinicopathological parameters of patients with HNSCC

Since we have illuminated that the expression level of PYGM and TNNC2 was down-regulated in HNSCC samples, then we focused on whether mRNA expression of these genes was related to cancer grade or stage in individual patients. As shown in Figure 6A, the results indicated that there are significant differences in PYGM expression between patients with Grade 1 and Grade 2. While comparing with HNSCC patients with Grade 1, patients with a more advanced grade (Grade 2 and Grade 3) had a relative lower level of TNNC2 expression (Figure 6B). However, no statistical significance was observed among patients with different stages in both PYGM and TNNC2 expression, implying they may not influence the stage of HNSCC (Figure 6C,D). Taken together, these results indicated that PYGM and TNNC2 may serve as a biomarker for pathological grade in HNSCC and may somehow participate in the progression of HNSCC.

Prognostic value of PYGM and TNNC2 in HNSCC patients

To further explore the prognostic values of the mRNA expression of PYGM and TNNC2 in HNSCC patients, we conducted survival assay in UCSC Xena database. According to the survival curves in Figure 6E,F, patients with a higher level of PYGM had a better prognosis while the level of TNNC2 seemed to have no statistical influence on the survival time of HNSCC patients. Therefore, the present study suggested that PYGM may be a promising novel predictor for the prognosis of HNSCC.

Genetic alteration of PYGM and TNNC2 in HNSCC

To investigate whether genetic alterations of PYGM and TNNC2 play significant roles in HNSCC, we analyzed the genetic alterations of these two genes by using the cBioPortal online tool. From the OncoPrint schematic, gene alteration of PYGM and TNNC2 was 4 and 0.8% in HNSCC samples, respectively (Figure 7A). These included gene amplification, missense mutation and gene deep deletion. We further evaluated the relationship between gene alteration of these two genes and the survival of HNSCC patients. However, there are no significant differences in OS and DFS in cases with or without alterations in PYGM or TNNC2 (P-values, 0.423 and 0.367, respectively) (Figure 7B,C). All the above data suggested that genetic alteration of PYGM and TNNC2 did not play essential roles in HNSCC progression.

Validation of PYGM and TNNC2 expression in HNSCC cell lines and clinical samples

To further validate the expression level of PYGM and TNNC2 in HNSCC, we performed qPCR experiment in HNSCC cell lines and HNSCC tissues. As suggested by Figure 8A, compared with normal oral cell lines, HNSCC cell lines presented a remarkable lower level of PYGM and TNNC2. Moreover, the results of qPCR conducted on 62 pairs of HNSCC tissues and adjacent matched tissues illustrated that PYGM and TNNC2 were significantly down-regulated in HNSCC tissues (Figure 8B). Interestingly, statistical analysis suggested that PYGM and TNNC2 were both closely associated with smoking and drinking behaviors, with lower levels in patients with smoking or drinking history (Figure 8C,D). We suspected that the components in tobacco or alcohol may influence the expression of PYGM and TNNC2 and further investigation needs to be implemented. What is more, we evaluated the relationship of PYGM and TNNC2 with clinicopathological parameters of patients with HNSCC. It was found that the down-regulation of PYGM and TNNC2 in HNSCC was positively associated with lymph node metastasis and advanced tumor stage (Figure 8E,F; Tables 3 and 4), implying their potential role in HNSCC metastasis and progression. As concluded in Figure 8G,H, most of HNSCC cases came from tongue and gingival and it was suggested that the relative level of PYGM and TNNC2 in patients with gingival cancer were higher than those with tongue tumor. We supposed that the expression of PYGM and TNNC2 may be somehow related with the disease site of HNSCC.
Figure 8
Validation of mRNA level of PYGM and TNNC2 in HNNSCC cell lines and tissues

(A) The relative expression of PYGM and TNNC2 in HNSCC cell lines and compared normal cells. (B) PYGM and TNNC2 mRNA levels were determined by qPCR in 62 pairs of HNSCC tumor samples and adjacent normal tissues. (C,D) A significant lower level of PYGM and TNNC2 were observed in HNSCC patients with smoking or drinking history. (E,F) The relative lower expression level of PYGM and TNNC2 was correlated with lymph node metastasis and advanced TNM stage in HNSCC patients. (G,H) The expression of PYGM and TNNC2 in HNSCC patients with different disease sites.

Figure 8
Validation of mRNA level of PYGM and TNNC2 in HNNSCC cell lines and tissues

(A) The relative expression of PYGM and TNNC2 in HNSCC cell lines and compared normal cells. (B) PYGM and TNNC2 mRNA levels were determined by qPCR in 62 pairs of HNSCC tumor samples and adjacent normal tissues. (C,D) A significant lower level of PYGM and TNNC2 were observed in HNSCC patients with smoking or drinking history. (E,F) The relative lower expression level of PYGM and TNNC2 was correlated with lymph node metastasis and advanced TNM stage in HNSCC patients. (G,H) The expression of PYGM and TNNC2 in HNSCC patients with different disease sites.

Close modal

HNSCC ranks as the sixth most prevalent cancer and constitutes approximately 5% of all malignancies worldwide [27]. The recurrence, metastasis or therapeutic resistance may result in the poor prognosis of HNSCC patients in spite of great improvement on therapy methods. Although there have been numerous studies on a variety of potential biomarkers, the molecular mechanisms underlying the pathogenesis of HNSCC remains elusive [28]. Hence, identifying novel and reliable diagnostic markers and treatment targets is challenging and urgently needed.

HNSCC is characterized by aberrant molecular signatures with various activated oncogenes and inactivated tumor suppressor genes [29]. Recently, a good number of studies have discovered novel genes associated with HNSCC in different ways. It was verified that increased expression of BAG-1 might be a biomarker for cisplatin resistance in patients with primary or recurrent HNSCC [30]. Also, high PD-1/PD-L1 expression was illustrated to be strongly related with radiosensitivity, suggesting that the combination of radiotherapy and anti-PD-1/PD-L1 therapy may improve the prognosis of HNSCC patients [31]. Therefore, the clarification of molecular mechanisms underlying HNSCC development and progression would assist us to discover novel treatment strategy.

Gene expression profile microarrays provide data from different samples and can be utilized for the identification of DEGs by bioinformatics analysis [32]. A multiple of public online databases allowed for data mining and the comprehensive analysis of various databases enhanced the reliability and accuracy of the results. In the present study, we made use of three independent GEO datasets to detect DEGs in HNSCC. Given the GPL platform of these datasets were different (summarized in Table 1), we should analyze original CEL data in R language, respectively. We utilized the affy package to perform the background correction and data normalization, including conversion of raw data formats, imputation of missing values and background correction. Then, the samples were subjected to differential expression analysis by using the Limma package. Through the combinations of the results from different microarrays based on diverse sample origins, the statistical power is increased, the predictive power is more accurate and the bias of individual studies can be overcome. Finally, we screened out eleven DEGs commonly abnormally expressed in three GEO datasets. To trace the origins of gene dysregulation, we constructed an interactive network of eleven DEGs along with the top 50 frequently altered neighboring genes on HNSCC. From the gene network, some important cancer-related genes such as EGFR, TP53 and MYC, could be observed to be correlated with DEGs directly or indirectly. Since these genes have been demonstrated by numerous studies for their significant functions in mediating HNSCC progression and even exploited as potential targets for gene therapy [33–35]. We hypothesized that the gene dysregulation of DEGs might be explained by the abnormal activities of these genes and further studies should be implemented to validate the interaction between DEGs and these neighboring genes in the tumorigenesis of HNSCC. Subsequently, GO analysis conducted on the genes of the network showed that they mainly participate in the cellular homeostasis, regulation of apoptosis and cell death and response to virus or hypoxia, which are all important biological process involved in cancer development [36–38]. Meanwhile, most of the genes were enriched in extracellular region and it was established that the complex interactions between tumor cells and the ECM may play significant roles in tumor metastasis [39]. What is more, KEGG analysis results manifested genes in the network are closely associated with Toll-like receptor signaling pathway, Jak-STAT pathway and MAPK signaling, which were reported to be intimately involved in HNSCC progression [40–42]. KEGG pathways can help us have a better understanding of the HNSCC pathology and provide us with a guide to the potential mechanisms underlying DEGs.

Due to a majority of the DEGs investigated by previous studies, we then selected two genes, PYGM and TNNC2, for further analyses. Through data mining in a variety of databases, we illustrated that PYGM and TNNC2 were remarkably underexpressed in HNSCC, implying their possible role in cancer progression. Furthermore, the results of the research conducted by Dickinson et al. [43] showed that there was statistically significant difference regarding the protein level of TNNC2 in oral squamous cell carcinoma and healthy tissues. Specifically, the protein level of TNNC2 in oral squamous cell carcinoma samples was 11.5-fold lower than that in normal tissues. It showed that our findings were consistent with this study and suggested that TNNC2 may be exploited as efficient biomarkers for oral squamous cell carcinoma both in mRNA and protein levels. Subsequently, we further validated the relative lower mRNA levels of PYGM and TNNC2 in HNSCC tissues compared with adjacent normal tissues and demonstrated that the down-regulation of PYGM and TNNC2 was positively associated with lymph node metastasis and advanced tumor stage. PYGM, a gene involved in glycogen metabolism [44], has been illuminated to be down-regulated in breast cancer and its expression was related with patients survival time [45]. More importantly, previous study suggested that glycogen metabolism may play a central role in cancer progression. It was demonstrated that compared with normal tissues, solid tumor had a higher glycogen content and the correct storage and management of glycogen may be relevant to cancer cell survival [46–49], implying the possible modulatory effect of PYGM in cancer. Also, another study speculated that PYGM may be implicated in gastric cancer via the insulin resistance pathway [50]. However, there existed no findings about TNNC2 with any tumor type and we supposed that further studies are warranted to elucidate its potential functions in HNSCC carcinogenesis.

In conclusion, the present study identified two dysregulated genes, PYGM and TNNC2, in HNSCC and illustrated that they may serve as effective therapeutic targets and efficient biomarkers for HNSCC. However, there exists some limitations in the present study. The size of HNSCC samples used in our study for the validation of PYGM and TNNC2 is relatively small. Moreover, the biological effect of PYGM and TNNC2 on HNSCC progression and the underlying mechanism have not been investigated. Therefore, more comprehensive studies still need to be conducted in the future.

Y.Y. and Y.J. conceived and designed the study. Y.J. analyzed and interpreted the data. Y.J. wrote the manuscript and Y.Y. made a revised version.

This work was supported by the National Key R&D Program of China [grant numbers 2017YFC 0840100, 2017YFC 0840110].

The authors declare that there are no competing interests associated with the manuscript.

DAVID

Database for Annotation, Visualization and Integrated Discovery

DEG

differentially expressed gene

DFS

disease-free survival

FC

fold change

GEO

Gene Expression Omnibus

GEPIA

Gene Expression Profiling Interactive Analysis

GO

Gene Ontology

HNSCC

head and neck squamous cell carcinoma

KEGG

Kyoto Encyclopedia of Genes and Genomes

OS

overall survival

qPCR

quantitative real-time polymerase chain reaction

SD

standard deviation

TCGA

The Cancer Genome Atlas

1.
Siegel
R.
,
Ma
J.
,
Zou
Z.
and
Jemal
A.
(
2014
)
Cancer statistics, 2014
.
CA Cancer J. Clin.
64
,
9
29
[PubMed]
2.
Shield
K.D.
,
Ferlay
J.
,
Jemal
A.
,
Sankaranarayanan
R.
,
Chaturvedi
A.K.
,
Bray
F.
et al.
(
2017
)
The global incidence of lip, oral cavity, and pharyngeal cancers by subsite in 2012
.
CA Cancer J. Clin.
67
,
51
64
[PubMed]
3.
Leemans
C.R.
,
Braakhuis
B.J.
and
Brakenhoff
R.H.
(
2011
)
The molecular biology of head and neck cancer
.
Nat. Rev. Cancer
11
,
9
22
[PubMed]
4.
Matar
N.
and
Haddad
A.
(
2011
)
New trends in the management of head and neck cancers
.
J. Med. Liban.
59
,
220
226
[PubMed]
5.
Fei
H.J.
,
Chen
S.C.
,
Zhang
J.Y.
,
Li
S.Y.
,
Zhang
L.L.
,
Chen
Y.Y.
et al.
(
2018
)
Identification of significant biomarkers and pathways associated with gastric carcinogenesis by whole genome-wide expression profiling analysis
.
Int. J. Oncol.
52
,
955
966
[PubMed]
6.
Barrett
T.
,
Wilhite
S.E.
,
Ledoux
P.
,
Evangelista
C.
,
Kim
I.F.
,
Tomashevsky
M.
et al.
(
2013
)
NCBI GEO: archive for functional genomics data sets–update
.
Nucleic Acids Res.
41
,
D991
D995
[PubMed]
7.
Tomczak
K.
,
Czerwinska
P.
and
Wiznerowicz
M.
(
2015
)
The Cancer Genome Atlas (TCGA): an immeasurable source of knowledge
.
Contemp. Oncol. (Pozn)
19
,
A68
A77
[PubMed]
8.
Rung
J.
and
Brazma
A.
(
2013
)
Reuse of public genome-wide gene expression data
.
Nat. Rev. Genet.
14
,
89
99
[PubMed]
9.
Kannan
L.
,
Ramos
M.
,
Re
A.
,
El-Hachem
N.
,
Safikhani
Z.
,
Gendoo
D.M.
et al.
(
2016
)
Public data and open source tools for multi-assay genomic investigation of disease
.
Brief. Bioinform.
17
,
603
615
[PubMed]
10.
Liu
X.
,
Gao
Y.
,
Zhao
B.
,
Li
X.
,
Lu
Y.
,
Zhang
J.
et al.
(
2015
)
Discovery of microarray-identified genes associated with ovarian cancer progression
.
Int. J. Oncol.
46
,
2467
2478
[PubMed]
11.
Chen
D.
,
Chen
W.
,
Xu
Y.
,
Zhu
M.
,
Xiao
Y.
,
Shen
Y.
et al.
(
2019
)
Upregulated immune checkpoint HHLA2 in clear cell renal cell carcinoma: a novel prognostic biomarker and potential therapeutic target
.
J. Med. Genet.
56
,
43
49
[PubMed]
12.
Haddad
R.I.
and
Shin
D.M.
(
2008
)
Recent advances in head and neck cancer
.
N. Engl. J. Med.
359
,
1143
1154
[PubMed]
13.
Huang
D.W.
,
Sherman
B.T.
,
Tan
Q.
,
Kir
J.
,
Liu
D.
,
Bryant
D.
et al.
(
2007
)
DAVID Bioinformatics Resources: expanded annotation database and novel algorithms to better extract biology from large gene lists
.
Nucleic Acids Res.
35
,
W169
W175
,
[PubMed]
14.
Rhodes
D.R.
,
Yu
J.
,
Shanker
K.
,
Deshpande
N.
,
Varambally
R.
,
Ghosh
D.
et al.
(
2004
)
ONCOMINE: a cancer microarray database and integrated data-mining platform
.
Neoplasia
6
,
1
6
[PubMed]
15.
Tang
Z.
,
Li
C.
,
Kang
B.
,
Gao
G.
,
Li
C.
and
Zhang
Z.
(
2017
)
GEPIA: a web server for cancer and normal gene expression profiling and interactive analyses
.
Nucleic Acids Res.
45
,
W98
W102
[PubMed]
16.
Chandrashekar
D.S.
,
Bashel
B.
,
Balasubramanya
S.A.H.
,
Creighton
C.J.
,
Ponce-Rodriguez
I.
,
Chakravarthi
B.
et al.
(
2017
)
UALCAN: a portal for facilitating tumor subgroup gene expression and survival analyses
.
Neoplasia
19
,
649
658
[PubMed]
17.
Gao
J.
,
Aksoy
B.A.
,
Dogrusoz
U.
,
Dresdner
G.
,
Gross
B.
,
Sumer
S.O.
et al.
(
2013
)
Integrative analysis of complex cancer genomics and clinical profiles using the cBioPortal
.
Sci. Signal.
6
,
pl1
,
18.
Patel
V.
,
Senderowicz
A.M.
,
Pinto
D.
Jr
,
Igishi
T.
,
Raffeld
M.
,
Quintanilla- Martinez
L.
et al.
(
1998
)
Flavopiridol, a novel cyclin- dependent kinase inhibitor, suppresses the growth of head and neck squamous cell carcinomas by inducing apoptosis
.
J. Clin. Invest.
102
,
1674
1681
[PubMed]
19.
Choi
P.
,
Jordan
C.D.
,
Mendez
E.
,
Houck
J.
,
Yueh
B.
,
Farwell
D.G.
et al.
(
2008
)
Examination of oral cancer biomarkers by tissue microarray analysis
.
Arch. Otolaryngol. Head Neck Surg.
134
,
539
546
[PubMed]
20.
Zhang
X.
,
Zhang
L.
,
Tan
X.
,
Lin
Y.
,
Han
X.
,
Wang
H.
et al.
(
2018
)
Systematic analysis of genes involved in oral cancer metastasis to lymph nodes
.
Cell Mol. Biol. Lett.
23
,
53
[PubMed]
21.
Zhong
L.P.
,
Zhang
L.
,
Yang
X.
,
Pan
H.Y.
,
Zhou
X.J.
,
Wei
K.J.
et al.
(
2009
)
Comparative proteomic analysis of differentially expressed proteins in an in vitro cellular carcinogenesis model of oral squamous cell carcinoma
.
Proteomics Clin. Appl.
3
,
322
337
[PubMed]
22.
Sumino
J.
,
Uzawa
N.
,
Okada
N.
,
Miyaguchi
K.
,
Mogushi
K.
,
Takahashi
K.
et al.
(
2013
)
Gene expression changes in initiation and progression of oral squamous cell carcinomas revealed by laser microdissection and oligonucleotide microarray analysis
.
Int. J. Cancer
132
,
540
548
[PubMed]
23.
Pappa
E.
,
Nikitakis
N.
,
Vlachodimitropoulos
D.
,
Avgoustidis
D.
,
Oktseloglou
V.
and
Papadogeorgakis
N.
(
2014
)
Phosphorylated signal transducer and activator of transcription-1 immunohistochemical expression is associated with improved survival in patients with oral squamous cell carcinoma
.
J. Oral Maxillofac. Surg.
72
,
211
221
[PubMed]
24.
Lv
Z.
,
Wu
X.
,
Cao
W.
,
Shen
Z.
,
Wang
L.
,
Xie
F.
et al.
(
2014
)
Parathyroid hormone-related protein serves as a prognostic indicator in oral squamous cell carcinoma
.
J. Exp. Clin. Cancer Res.
33
,
100
[PubMed]
25.
Sepiashvili
L.
,
Hui
A.
,
Ignatchenko
V.
,
Shi
W.
,
Su
S.
,
Xu
W.
et al.
(
2012
)
Potentially novel candidate biomarkers for head and neck squamous cell carcinoma identified using an integrated cell line-based discovery strategy
.
Mol. Cell. Proteomics
11
,
1404
1415
[PubMed]
26.
Makinen
L.K.
,
Hayry
V.
,
Atula
T.
,
Haglund
C.
,
Keski-Santti
H.
,
Leivo
I.
et al.
(
2012
)
Prognostic significance of matrix metalloproteinase-2, -8, -9, and -13 in oral tongue cancer
.
J. Oral Pathol. Med.
41
,
394
399
[PubMed]
27.
Yang
B.
,
Chen
Z.
,
Huang
Y.
,
Han
G.
and
Li
W.
(
2017
)
Identification of potential biomarkers and analysis of prognostic values in head and neck squamous cell carcinoma by bioinformatics analysis
.
Onco Targets Ther.
10
,
2315
2321
[PubMed]
28.
Gunther
K.
,
Foraita
R.
,
Friemel
J.
,
Gunther
F.
,
Bullerdiek
J.
,
Nimzyk
R.
et al.
(
2017
)
The stem cell factor HMGA2 is expressed in non-HPV-associated head and neck squamous cell carcinoma and predicts patient survival of distinct subsites
.
Cancer Epidemiol. Biomark. Prev.
26
,
197
205
[PubMed]
29.
Califano
J.
,
van der Riet
P.
,
Westra
W.
,
Nawroz
H.
,
Clayman
G.
,
Piantadosi
S.
et al.
(
1996
)
Genetic progression model for head and neck cancer: implications for field cancerization
.
Cancer Res.
56
,
2488
2492
[PubMed]
30.
Liu
S.
,
Ren
B.
,
Gao
H.
,
Liao
S.
,
Zhai
Y.X.
,
Li
S.
et al.
(
2017
)
Over-expression of BAG-1 in head and neck squamous cell carcinomas (HNSCC) is associated with cisplatin-resistance
.
J. Transl. Med.
15
,
189
[PubMed]
31.
Lyu
X.
,
Zhang
M.
,
Li
G.
,
Jiang
Y.
and
Qiao
Q.
(
2019
)
PD-1 and PD-L1 expression predicts radiosensitivity and clinical outcomes in head and neck cancer and is associated with HPV infection
.
J. Cancer
10
,
937
948
[PubMed]
32.
Zou
A.E.
,
Ku
J.
,
Honda
T.K.
,
Yu
V.
,
Kuo
S.Z.
,
Zheng
H.
et al.
(
2015
)
Transcriptome sequencing uncovers novel long noncoding and small nucleolar RNAs dysregulated in head and neck squamous cell carcinoma
.
RNA
21
,
1122
1134
[PubMed]
33.
Tian
Y.
,
Lin
J.
,
Tian
Y.
,
Zhang
G.
,
Zeng
X.
,
Zheng
R.
et al.
(
2018
)
Efficacy and safety of anti-EGFR agents administered concurrently with standard therapies for patients with head and neck squamous cell carcinoma: a systematic review and meta-analysis of randomized controlled trials
.
Int. J. Cancer
142
,
2198
2206
[PubMed]
34.
Zhou
G.
,
Liu
Z.
and
Myers
J.N.
(
2016
)
TP53 mutations in head and neck squamous cell carcinoma and their impact on disease progression and treatment response
.
J. Cell. Biochem.
117
,
2682
2692
[PubMed]
35.
Strindlund
K.
,
Troiano
G.
,
Sgaramella
N.
,
Coates
P.J.
,
Gu
X.
,
Boldrup
L.
et al.
(
2017
)
Patients with high c-MYC-expressing squamous cell carcinomas of the tongue show better survival than those with low- and medium-expressing tumours
.
J. Oral Pathol. Med.
46
,
967
971
[PubMed]
36.
Wong
R.S.
(
2011
)
Apoptosis in cancer: from pathogenesis to treatment
.
J. Exp. Clin. Cancer Res.
30
,
87
[PubMed]
37.
Kujan
O.
,
Shearston
K.
and
Farah
C.S.
(
2017
)
The role of hypoxia in oral cancer and potentially malignant disorders: a review
.
J. Oral Pathol. Med.
46
,
246
252
[PubMed]
38.
Basanta
D.
and
Anderson
A.R.A.
(
2017
)
Homeostasis back and forth: an ecoevolutionary perspective of cancer
.
Cold Spring Harb. Perspect. Med.
7
,
[PubMed]
39.
Oudin
M.J.
,
Jonas
O.
,
Kosciuk
T.
,
Broye
L.C.
,
Guido
B.C.
,
Wyckoff
J.
et al.
(
2016
)
Tumor cell-driven extracellular matrix remodeling drives haptotaxis during metastatic progression
.
Cancer Discov.
6
,
516
531
[PubMed]
40.
Szczepanski
M.J.
,
Czystowska
M.
,
Szajnik
M.
,
Harasymczuk
M.
,
Boyiadzis
M.
,
Kruk-Zagajewska
A.
et al.
(
2009
)
Triggering of Toll-like receptor 4 expressed on human head and neck squamous cell carcinoma promotes tumor development and protects the tumor from immune attack
.
Cancer Res.
69
,
3105
3113
[PubMed]
41.
Lai
S.Y.
and
Johnson
F.M.
(
2010
)
Defining the role of the JAK-STAT pathway in head and neck and thoracic malignancies: implications for future therapeutic approaches
.
Drug Resist. Updat.
13
,
67
78
[PubMed]
42.
Ozawa
H.
,
Ranaweera
R.S.
,
Izumchenko
E.
,
Makarev
E.
,
Zhavoronkov
A.
,
Fertig
E.J.
et al.
(
2017
)
SMAD4 loss is associated with cetuximab resistance and induction of MAPK/JNK activation in head and neck cancer cells
.
Clin. Cancer Res.
23
,
5162
5175
[PubMed]
43.
Dickinson
A.
,
Saraswat
M.
,
Makitie
A.
,
Silen
R.
,
Hagstrom
J.
,
Haglund
C.
et al.
(
2018
)
Label-free tissue proteomics can classify oral squamous cell carcinoma from healthy tissue in a stage-specific manner
.
Oral Oncol.
86
,
206
215
[PubMed]
44.
Newgard
C.B.
,
Hwang
P.K.
and
Fletterick
R.J.
(
1989
)
The family of glycogen phosphorylases: structure and function
.
Crit. Rev. Biochem. Mol. Biol.
24
,
69
99
[PubMed]
45.
Dieci
M.V.
,
Smutna
V.
,
Scott
V.
,
Yin
G.
,
Xu
R.
,
Vielh
P.
et al.
(
2016
)
Whole exome sequencing of rare aggressive breast cancer histologies
.
Breast Cancer Res. Treat.
156
,
21
32
[PubMed]
46.
Skwarski
L.
,
Namiot
Z.
,
Stasiewicz
J.
,
Kemona
A.
,
Kralisz
M.
and
Gorski
J.
(
1998
)
Glycogen content in the gastric mucosa of partially resected stomach; a possible relationship with the development of cancer
.
Cancer Lett.
127
,
123
128
[PubMed]
47.
Yano
K.
,
Ohoshima
S.
,
Shimizu
Y.
,
Moriguchi
T.
and
Katayama
H.
(
1996
)
Evaluation of glycogen level in human lung carcinoma tissues by an infrared spectroscopic method
.
Cancer Lett.
110
,
29
34
[PubMed]
48.
Pelletier
J.
,
Bellot
G.
,
Gounon
P.
,
Lacas-Gervais
S.
,
Pouyssegur
J.
and
Mazure
N.M.
(
2012
)
Glycogen synthesis is induced in hypoxia by the hypoxia-inducible factor and promotes cancer cell survival
.
Front. Oncol.
2
,
18
[PubMed]
49.
Favaro
E.
,
Bensaad
K.
,
Chong
M.G.
,
Tennant
D.A.
,
Ferguson
D.J.
,
Snell
C.
et al.
(
2012
)
Glucose utilization via glycogen phosphorylase sustains proliferation and prevents premature senescence in cancer cells
.
Cell Metab.
16
,
751
764
[PubMed]
50.
Zhang
W.
,
Liu
S.
,
Zhan
H.
,
Yan
Z.
and
Zhang
G.
(
2018
)
Transcriptome sequencing identifies key pathways and genes involved in gastric adenocarcinoma
.
Mol. Med. Rep.
18
,
3673
3682
[PubMed]
This is an open access article published by Portland Press Limited on behalf of the Biochemical Society and distributed under the Creative Commons Attribution License 4.0 (CC BY).