Introduction
Head and neck squamous cell carcinomas (HNSCC), is the fifth most common type of cancer worldwide, with the most frequent type being laryngeal squamous cell carcinoma (LSCC) (1). As an aggressive neoplasm, LSCC is the second most predominant upper respiratory tract tumor and has a high mortality (2). In addition, it was reported that LSCC had a more frequent incidence among middle-aged and elderly males (3). Although advanced surgical treatments such as laryngectomy, radiotherapy and chemotherapy have been developed in the diagnosis and prevention of early LSCC, the only option for advanced LSCC has been radio-chemotherapy, which often did not achieve a positive clinical response, resulting in a poor prognosis of advanced LSCC (4,5). The 5-year survival rate for LSCC was 64.2%, as indicated in a meta-analysis study (6). Therefore, investigations have been conducted to identify new approaches such as biological therapeutic strategy, in order to provide more effective and accurate predictions and preventions of LSCC.
Previously, a cohort of gene markers was identified that was associated with LSCC, including matrix metalloproteinase 2 (MMP2), cluster of differentiation 24 (CD24) and U3 small nucleolar ribonucleoprotein protein (IMP3) (7–9). Notably, the mRNA expression of dopa decarboxylase (DDC) gene, which encodes the enzyme L-DOPA decarboxylase, was associated with tumor stage in HNSCC and thus this gene was proposed as a potential prognostic biomarker of LSCC (10). Furthermore, several genes were significantly distinct between early and late tumor stages of LSCC, including ERBB4, BCL7A and RECQL4 (11). A recent study included patients with LSCC in different tumor stages to clarify new biomarkers that had driven the progression of LSCC. Based on their finding that HOXC8 and HOXD11 were essential for FADU cell expansion, the HOX family members were suggested to be correlated with the development of LSCC (12). However, alterations of the specific genes expressed in the early and late stage of LSCC were rarely involved. Thus, we employed the microarray data GSE59102 and re-analyzed it, to identify crucial genes in both stages using differential analysis. Additionally, we emphasized the overlapped DEGs (differentially expressed genes) in the two stages and examined the potential associations between these DEGs using co-expression network and protein-protein interaction (PPI) network analyses. Transcription factors (TFs) and tumor-associated genes (TSGs) were also predicted following a comparison of the relevant databases. Using these bioinformatical methods, we aimed to identify novel gene markers exhibiting a tumor stage-dependent expression and explore the underlying regulatory mechanisms of LSCC progression.
Materials and methods
Microarray data
The expression profile with the accession number GSE59102 (12) was downloaded from the public Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo) database. The data set comprised 29 cancer and 13 margin samples (control) from patients undergoing surgical ablation of LSCC. The cancer samples were obtained from two different tumor stages: 14 in the early and 15 in the advanced stage. The platform for the gene expression calculation was Agilent-014850 Whole Human Genome Microarray 4×44K G4112F (Agilent Technologies, Palo Alto, CA, USA) and the annotation files in this platform were also downloaded for the analysis in the present study.
Data preprocessing and differential analysis
Following normalization, the data were presented in the form of a gene expression matrix. Uniform distribution of the gene mean value in different samples was obtained and the probe value was mapped into the gene symbol based on the annotation files. If ≥1 probe was mapped into one gene, the average probe value was calculated as the final expression value of this gene.
The differential analysis was performed between early and control samples and between advanced and control samples, respectively, recruiting the limma (Linear Models for Microarray Analysis) package of Bioconductor R (http://www.bioconductor.org/packages/release/bioc/html/limma.html) (13). The t-test method in limma was used to calculate the significance for the DEGs selection and the P-value was adjusted to FDR (false discover rate). The threshold were FDR <0.01 and |log2 (fold change)| ≥4.
Functional enrichment analysis of the DEGs
The gene ontology (GO, http://www.geneontology.org/) (14), biological process (BP) and Kyoto Encyclopedia of Genes and Genomes (KEGG, http://www.genome.jp/kegg/pathway.html) (15) pathway enrichment analyses were conducted for two sets of DEGs, respectively, using the Database for Annotation, Visualization and Integration Discovery (DAVID, http://david.abcc.Ncifcrf.gov/) database (16), which revealed the enriched GO terms and KEGG pathways based on the hypergeometric distribution. The cut-off value for significant terms was P<0.05.
Identification of TFs and TSGs
The two sets of DEGs were mapped into the TRASFAC database (17) to screen the potential TFs among these DEGs. Additionally, following the combination with the information from the tumor-suppressor gene (TSG) database (18) and TAG (http://www.binfo.ncku.edu.tw/TAG/) database, all known oncogenes and TEGs of the DEGs were extracted.
Co-expressed analysis of DEGs identified in both early- and advanced-stage samples
Considering that there may be several genes that are important in the early- and advanced-stage samples and that there may be certain correlations between these genes, the two sets of identified DEGs were merged to establish a co-expressed network, using criteria of Pearson's correlation coefficients >0.9.
Association analysis between the DEGs and different tumor stages
The overlapped genes in the two sets of DEGs were screened out, and any up- or downregulated genes showing a tumor stage-dependent increasing (the gene expression in different samples had the following trend: advanced > early > control) or decreasing expression (the gene expression in different samples had the following trend: advanced < early < control) were selected.
Subsequently, a KEGG pathway enrichment analysis was carried out for the DEGs that were positively and negatively associated with the tumor stage. As proteins often co-operate to exert their biological functions, examination of the PPI between them significantly contributes to identifying the underlying molecular mechanisms. Therefore, the above DEGs were mapped into the STRING (Search Tool for the Retrieval of Interacting Genes) (http://string-db.org/) database (19), to search the potential relationships between them from the protein level. The pairwise interactions of the DEGs (the protein products) with a combined score of >0.4 were selected to establish the PPI network, as visualized by the Cytoscape software (http://cytoscape.org/) (20).
Prediction of other disease-related genes
To examine the association between LSCC stage and DEGs, the screened overlapped tumor stage-dependent DEGs were mapped into the Comparative Toxicogenomics Database (CTD, http://ctdbase.org/) (21), which collected disease-related genes that were validated by experiment or reported in the literature, to search for the exact LSCC-associated DEGs.
Results
Identification of DEGs in different stages
Based on the aforementioned selection criteria, a set of 258 upregulated and 438 downregulated DEGs between early-stage and control samples, and a cohort of 294 upregulated and 328 downregulated DEGs between advanced-stage and control samples were identified. The heat map of the cluster analysis of the DEGs is shown in Fig. 1.
Dysregulated biological functions and pathways of the DEGs
To characterize the biological functions altered by the DEGs, the GO functional and KEGG pathway enrichment analyses were carried out. As shown in Table I, of the DEGs between early-stage and control samples, the upregulated genes were significantly enriched in the cell cycle-related process such as M phase (GO: 0000279), cell cycle phase (GO: 0022403), cell cycle (GO: 0007049) and cell cycle process (GO: 0022402) and pathways including cell cycle (Hsa04110), cytokine-cytokine receptor interaction (Hsa04060) and homologous recombination (Hsa03440), while the over-represented BP terms and pathways for the downregulated genes were chemical homeostasis (GO:0048878), cation homeostasis (GO:0055080), response to wounding (GO:0009611), cellular cation homeostasis (GO:0030003), salivary secretion (Hsa04970), drug metabolism - cytochrome P450 (Hsa00982), metabolism of xenobiotics by cytochrome P450 (Hsa00980) and retinol metabolism (Hsa00830).
The result of the functional enrichment analysis of the DEGs between the advanced-stage and control samples is presented in Table II. The upregulated genes were predominantly enriched in BP terms including ectoderm development (GO: 0007398), epidermis development (GO: 0008544), skeletal system development (GO: 0001501) and epidermal cell differentiation (GO: 0009913), as well as the pathways involved in cytokine-cytokine receptor interaction (Hsa04060), ECM-receptor interaction (Hsa04512) and chemokine signaling pathway (Hsa04062). For the downregulated genes, the notable BP terms were oxidation reduction (GO: 0055114), fatty acid metabolic process (GO: 0006631) and regulation of epithelial cell proliferation (GO: 0050678), while the predominant pathways included were drug metabolism-cytochrome P450 (Hsa00982), metabolism of xenobiotics by cytochrome P450 (Hsa00980) and retinol metabolism (Hsa00830).
TFs, oncogenes and TSGs identification among the DEGs
Based on the information in the relevant databases, the potential TFs, oncogenes and TSG among the DEGs were screened out. As shown in Table III, a set of upregulated HOX family members as well as the downregulated SMAD9, PPARG, LMO3, NKX3-1 and KLF15 were identified as the TFs. The upregulated HMGA2, HOXA10, LAMC2, CEP55 (centrosomal protein 55 kDa) as well as LMO2, CD24 and LCN2 were predicted as the oncogenes. Additionally, the upregulated genes including MMP11, E2F1 and BRCA2 as well as the downregulated MUC1, SERPINI2 and CEACAM1 were identified as TSGs, amongst the DEGs in the early-stage samples.
By contrast, in the advanced-stage samples, a number of upregulated HOX family members were identified as TFs, together with the downregulated DEGs such as NKX6-2, ELF5 and NKX3-1. The potential oncogenes were the upregulated LAMC2, HOXA10 and CEP55 as well as the downregulated LCN2 and WISP2. The upregulated MMP11, E2F1 and DEC1 as well as the downregulated DEGs including DEFB1, MUC1 and ACCN1 were deemed as the TSGs.
Co-expressed network of all the DEGs and the functional enrichment analysis
According to the Pearson's correlation coefficients, the co-expressed network was constructed, comprising 494 pairwise interactions and 205 DEGs (Fig. 2). A gene in the network served as a 'node' and the 'degree' of a node referred to the interaction number of that specific gene. In the network, the nodes with high degrees were prominent such as BUB1 (degree =35), BIRC5 (degree =31), TPX2 (degree =28), TTK (degree =25), FAM72D (degree =25), CEP55 (degree =24), CDC45 (degree =23), CDC25C (degree =23), E2F1 (degree =22) and DTL (degree =22).
The functional enrichment analysis revealed that, the BUB1, TTK and CDC45 genes in the network were mainly enriched in the cell cycle pathway, whereas BIRC5 was associated with the cancer pathway. Notably, E2F1 was involved in the two pathways.
PPI network of the tumor stage-dependent expressed DEGs and the enrichment analysis
Among the overlapped DEGs identified in the two samples, a set of 113 upregulated genes indicated an increasing expression in a tumor stage-dependent manner, while 115 downregulated genes exhibited a tumor stage-dependent decreasing expression. Combining the information from the STRING database, the PPI network of the tumor stage-dependent DEGs was constructed, comprising 125 interactions and 88 nodes (Fig. 3). As shown in Fig. 3, the highlighted nodes with high degrees were IL1B (degree =13), MMP9 (degree =10), MMP1 (degree =9), PLAU (degree =8), SERPINE1 (degree =7), MMP3 (degree =7) and COL1A1 (degree =7).
LSCC-related genes among the tumor stage-dependent genes
By combining the CTD database and the tumor stage-dependent genes, 65 upregulated and 67 downregulated DEGs were detected to be closely associated with LSCC via the curated literature. Notably, the ADH family members of ADH7 and ADH1B were considered as two markers of LSCC. The expression of ADH family members including ADH1A, ADH1C, ADH4 and ADH7 during different stages is shown in Fig. 4. Analysis of variance (ANOVA) of the expression of four genes among different samples are presented in Table IV, which indicated pronounced expression differences (FDR<0.05) among these samples.
Discussion
In the present study, we identified 696 DEGs (258 upregulated and 438 downregulated) between early-stage and control samples, which were enriched in the cell cycle-related pathways, as well as 622 DEGs (294 upregulated and 328 downregulated) between advanced-stage and control samples with the over-represented pathways including the cytokine-cytokine receptor interaction and ECM-receptor interaction pathways. In the co-expressed network, BUB1, TTK and E2F1, which were also predicted as potential TSGs, were prominent and mainly enriched in the cell cycle pathway. CEP55, also evident in the co-expressed network, was identified as an oncogene. The HOX family members were predicted as TFs in the two samples. The upregulated MMP1, MMP9 and MMP3 were the predominant nodes in the PPI network and had a tumor stage-dependent increasing expression pattern. Additionally, MMP3 was connected with MMP1 in the PPI network. ADH family (ADH1A, ADH1C, ADH4 and ADH7) were predicted to be correlated with LSCC.
The gene BUB1 encodes a serine/threonine-protein kinase that has a critical role in mitosis. It is essential for the localization of a mitotic checkpoint protein, BUBR1 (22), which was suggested to be associated with STK15 and this interaction was thought having significant roles in the tumorigenesis and development of LSCC (23). A low expression of the BUB1 protein was found to be highly correlated with the tumorigenesis and progression of LSCC and was proposed as a biomarker for the prognosis of LSCC (24). The present results indicated that BUB1 was markedly enriched in the cell cycle pathways, consistent with previous functional enrichment analysis in LSCC (25). Additionally, BUB1 was prominent in the co-expressed network, suggesting this gene has a significant role in the development of LSCC.
TTK is another mitotic checkpoint protein that is correlated with cell proliferation and is responsible for the accurate segregation of chromosomes during mitosis (26). Dysregulation of this gene was verified to be associated with human cancers such as bladder and ovarian cancer (27). Considering that this gene was a vital node in the co-expressed network and also involved in the cell cycle pathway in the present study, it may play a significant role in the progression of LSCC from the early to the advanced stage via participation of the cell cycle-mediated pathway.
The TF E2F1 belongs to the E2F family, which plays a critical role in the control of cell cycle and the action of tumor suppressor. A p53 independent target of E2F1, the death inducing protein KIAA0767 was differentially expressed in human lung SCC (27). Additionally, the overexpression of UHRF1, another target of E2F1, which played a pivotal function in promoting G1/S transition, was involved in the progression of LSCC and thus was considered as a useful biomarker for the detection of LSCC (29). Thus, the targets of E2F1 may contribute to the progression of LSCC. Combining our results that E2F1 was exhibited as a vital node in the co-expressed network, which was significantly enriched in the cell cycle pathway and predicted as a TSG, it may be hypothesized that the E2F1-mediated target expression may affect cell cycle and thus account for the progression of LSCC. However, more specific targets of the E2F1 should be established.
CEP55 is known to have a significant role in mitotic exit and cytokinesis. During interphase, CEP55 is localized to the mother centriole and the phosphorylation of this protein was required for the recruitment to midbody and cytokinesis (30). The oncogene FOXM1 (forkhead box M1) encoded protein is a transcriptional activator involved in cell proliferation. The evidence that two downstream targets of FOXM1, CEP55 and HELLS were positively associated with pre-malignancy and HNSCC development conferred the possibility that the two genes may be used as a biomarker in conjunction with FOXM1 for the early detection of LSCC (31). Of note, our results predicted CEP55 as an oncogene based on the information of TAG database, suggesting that CEP55 may exert its role in the development of LSCC as an oncogene.
Proteins encoded by homeobox genes are highly conserved TFs that have vital roles in morphogenesis in multicellular organisms. There are four similar homeobox gene clusters including HOXA, HOXB, HOXC and HOXD in mammals. It was found that HOX9 was highly expressed in LSCC and acted as an oncogene (32). HOTAIR is a long intergenic non-coding RNA, which is transcribed by HOX locus. Overexpression of HOTAIR was observed in primary LSCC and associated with advanced stage of LSCC and suggested as a novel biomarker in the treatment of LSCC (4). The aforementioned findings indicate that the HOX family has a significant role in the regulation of LSCC. Notably, HOX genes were identified as TFs in the present results. Although no relevant evidence was provided for the role of HOX genes as TFs, it was suggested that these genes are important in the progression of LSCC.
The proteolytic activity of MMPs against extracellular matrix confers the significant role of MMPs in the invasion and metastasis in tumors, which contributes to the progression into the advanced stage. The immunohistochemical experiment demonstrated the correlation of MMP1 and TIMP1 (a tissue inhibitor of MMPs), which was important in the progression of LSCC (33). Additionally, the MMP9 expression was significantly associated with the metastasis and prognosis of LSCC and was suggested as an indicator for lymph node metastasis and tumor invasion in LSCC (34,35). A further study corroborated that overexpression of the MMP9 protein is a characteristic in LSCC and closely correlated with the increase of the tumor-node-metastasis stage (36), consistent with our present observation that MMP9 was a tumor stage-dependent increasing DEG, suggesting the critical role of MMP9 in the development of LSCC.
MMP3 is also associated with metastatic progression. The upregulation of MMP3 was validated and identified as a signature in oral squamous cell carcinoma progression (37). Immunohistochemistry revealed the involvement of MMP3 in the events correlated with lip SCC (38). Although not validated, MMP3 was found to be upregulated in LSCC and associated with the development of invasion (39). Notably, MMP3 was confirmed to activate MMP1 and the high expression of MMP1 and MMP3 was connected with tumor invasion and metastasis in LSCC (40). Our PPI network indicated the relationships between MMP1and MMP3, which gave potent evidence that regulations may exist between the two MMPs in the progression of LSCC.
The ADH genes encode the alcohol dehydrogenase family protein. There was substantial evidence that ADH genes (such as ADH1B, ADH1C and ADH7) were involved in the risk of various types of cancer including gastric and upper aero-digestive cancers (41,42). With regard to their correlation with LSCC, a meta-analysis demonstrated that the polymorphisms of alcohol dehydrogenase genes (ADH1B and ADH1C) were closely associated with head and neck cancer risk (43). In the present study, four ADH genes such as ADH1A, ADH1C, ADH4 and ADH7 were downregulated in the two samples and showed a tumor stage-dependent decreasing pattern, suggesting that these genes could be used as signatures in the progression of LSCC.
PLAU encodes a serine protease responsible for the extracellular matrix degradation and tumor cell migration and proliferation (44), which may account for its positive tumor stage being correlated with the expression of LSCC, as revealed in our study. In breast cancer cells, PLAU was indicated to have a marked role in transformation by affecting invasiveness (45). Gene amplification has been extensively demonstrated to contribute to tumor cell growth. In prostatic cancer (PC), it was verified that PLAU amplification was preferentially found in the advanced stage, but not detected in the benign prostatic hyperplasia (46), suggesting that this gene may have a tumor stage-dependent expression pattern. Notably, MMP3 and PLAU were involved in the plasminogen activator pathway and induced by ERG (v-ets avian erythroblastosis virus E26 oncogene homolog) in PC (47). Given that MMP and PLAU were directly connected in our PPI network, it was hypothesized that a potential regulation may exist between these genes in the development of LSCC.
In conclusion, a set of critical genes in the progression of LSCC were identified including BUB1, TTK, E2F1, CEP55, MMPs (MMP1, MMP9 and MMP3), PLAU, HOX genes and ADH genes. Of these, MMPs, PLAU and ADH genes were used as tumor stage-dependent biomarkers for the prognosis of LSCC. During the progression of LSCC, BUB1, TTK, E2F1 and CEP55 may exert their roles via the involvement of the cell cycle pathway and HOX genes may serve as TFs, while MMP3 may be mediated by MMP1. Additionally, ADH genes may be used as indicators of the prognosis of LSCC. However, the findings remain to be verified.