Identification of key target genes and pathways in laryngeal carcinoma.

The purpose of the present study was to screen the key genes associated with laryngeal carcinoma and to investigate the molecular mechanism of laryngeal carcinoma progression. The gene expression profile of GSE10935 [Gene Expression Omnibus (GEO) accession number], including 12 specimens from laryngeal papillomas and 12 specimens from normal laryngeal epithelia controls, was downloaded from the GEO database. Differentially expressed genes (DEGs) were screened in laryngeal papillomas compared with normal controls using Limma package in R language, followed by Gene Ontology (GO) enrichment analysis and pathway enrichment analysis. Furthermore, the protein-protein interaction (PPI) network of DEGs was constructed using Cytoscape software and modules were analyzed using MCODE plugin from the PPI network. Furthermore, significant biological pathway regions (sub-pathway) were identified by using iSubpathwayMiner analysis. A total of 67 DEGs were identified, including 27 up-regulated genes and 40 down-regulated genes and they were involved in different GO terms and pathways. PPI network analysis revealed that Ras association (RalGDS/AF-6) domain family member 1 (RASSF1) was a hub protein. The sub-pathway analysis identified 9 significantly enriched sub-pathways, including glycolysis/gluconeogenesis and nitrogen metabolism. Genes such as phosphoglycerate kinase 1 (PGK1), carbonic anhydrase II (CA2), and carbonic anhydrase XII (CA12) whose node degrees were >10 were identified in the disease risk sub-pathway. Genes in the sub-pathway, such as RASSF1, PGK1, CA2 and CA12 were presumed to serve critical roles in laryngeal carcinoma. The present study identified DEGs and their sub-pathways in the disease, which may serve as potential targets for treatment of laryngeal carcinoma.

Introduction

Laryngeal cancer is one of the most common types of head and neck cancer, and the majority of laryngeal cancers are squamous cell carcinomas (1). The incidence rate of laryngeal cancer accounts for 1–5% of all malignancies, and is third most common head and neck cancer tumor (2). Laryngeal cancer has enormous functional and psychological consequences for the patients, in particular with regard to eating and communication (3). Laryngeal papillomas which are caused primarily by human papillomavirus (HPV) types 6 and 11 (4) are the most common benign tumors in the larynx. They are associated with a 3–7% risk of malignant transformation, in which irradiation and smoking appear to be risk factors (5). In order to further improve survival rates, diagnostic and preventive approaches, it is necessary to understand the mechanisms of laryngeal carcinogenesis.

The past decade has witnessed major progress in the understanding of the molecular alteration underlying the development of head and neck cancers (6,7). For laryngeal cancer, previous studies have revealed that the use of DNA microarray is helpful for the elucidation of laryngeal carcinogenesis (8,9). Colombo et al (10) identified that the differentially expressed genes (DEGs) in larynx tumors were involved in cellular processes associated with the cancer phenotype, including the cell cycle, DNA repair, and signal transduction among which cystatin B (CSTB) encoded the protein, which has previously been associated with antimetastatic function. Ma et al (6)reported that numerous molecular abnormalities encoding for cell cycle control and integrin-mediated cell adhesion, such as matrix metallopeptidase 12 (MMP12), keratin 19 (KRT19) and proline-rich protein BstNI subfamily 1 (PRB1) were identified in laryngeal carcinoma cells. In addition, the findings of Oblak et al (11) demonstrated that overexpression of suppressor of cytokine signaling 1 (SOCS1) was linked with reducing the aggressive and metastatic potential in laryngeal carcinoma due to the inhibition of cancer cell proliferation. Despite the fact that considerable efforts have been made over the years, the molecular mechanisms involved in the development and progression of laryngeal carcinoma have not been completely demonstrated and require further investigation.

Using the same gene expression profile, DeVoti et al (12) had demonstrated that altered expression of numerous genes was associated with cellular growth and differentiation and that there was a role for altered innate immunity in recurrent respiratory papillomatosis. In the present study, microarray analysis was used to identify DEGs in growing papillomas compared with adjacent laryngeal epithelia as controls. Comprehensive bioinformatics was used to analyze the significant pathways and functions and to construct the protein-protein interaction (PPI) network in addition to significant modules to identify the critical DEGs. Furthermore, significant biological pathway regions (sub-pathways) were predicted. The present study aimed to identify involvement of genes critical in laryngeal carcinoma and to get an improved understanding of the molecular circuitry in laryngeal carcinoma.

Materials and methods

Microarray data and data preprocessing

The mRNA expression profile of GSE10935 was downloaded from the National Center of Biotechnology Information (NCBI) Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo/), which was deposited by DeVoti et al (12). The platform is GPL96 [HG-U133A] Affymetrix Human Genome U133A Array and GPL571 [HG-U133A_2] Affymetrix Human Genome U133A 2.0 Array. In this dataset, there were 12 specimens from papillomas and 12 specimens from adjacent clinically normal, laryngeal epithelia from patients with recurrent respiratory papillomatosis.

The probes were annotated with gene symbols according to the annotation information in different platforms. Gene symbols which were common across the GPL96 and GPL571 platforms were retained to the subsequent analysis. Furthermore, cross slide normalization was performed by using the surrogate variable analysis (SVA) package (13) in R (14), and within slide normalization was performed by using preprocessCore package (15) in R. The gene expression matrix of specimens was received.

Screening of DEGs

A t-test (16) in the limma (17) package in R (14) was used to identify DEGs. Up-regulated and down-regulated genes were identified between papillomas and normal controls. The Bonferroni method (18) was applied to perform a multiple testing correction. The threshold for the DEGs was set as corrected P-value <0.05 and | log2 fold change (FC) | ≥1. In addition, the clustering analysis of DEGs was represented by a heat-map using gplots (19) in R (14).

Functional and pathway enrichment analysis of DEGs

In order to facilitate the functional annotation and pathway analysis, all the DEGs were analyzed using the Database for Annotation Visualization and Integrated Discovery (DAVID) (20) to perform the Gene Ontology (GO) (21) and Kyoto Encyclopedia of Genes and Genomes (KEGG) (22) analysis. The human genome was used as the background list and human was chosen as the species. The other default parameters of the tools were kept unchanged. The GO terms and KEGG pathways that were enriched by >5 genes and P<0.05 was considered to indicate a statistically significant difference.

PPI network construction and modules selection

The PPI network is represented by an undirected graph with nodes indicating the genes and edges indicating the mapped interactions of the proteins encoded by the genes (23). In the present study, all DEGs were imported into Cytoscape plugin to create network visualizations. The source of the interaction network databases were the Human Protein Reference Database (HPRD) (24) (http://www.hprd.org/) and the Biological General Repository for Interaction Datasets (BioGRID) (25) (http://thebiogrid.org/) database. Then the resulting PPI network was subjected to module analysis with the Plugin MCODE (26) with the default parameters (Degree cutoff ≥2, Node score cutoff ≥2, K-core ≥2, and Max depth=100).

Risk sub-pathway analysis

iSubpathwayMiner (http://cran.r-project.org/web/packages/iSubpathwayMiner/) is an R platform for graph-based construction and analysis of KEGG pathways (27). For a given KEGG pathway, the sub-pathways were obtained by searching all possible paths between start-points (membrane receptors or their ligands) and end-points (transcriptional factors or their immediate targets) in the adjacency matrix generated by node relationships (28).

Disease associated risk sub-pathway creation and presentation were based on iSubpathway Miner (28). The genes which had higher degrees in the PPI network enriched in the sub-pathways were identified as the critical genes in papilloma.

Results

Data preprocessing and DEGs screening

After normalization, 22,283 probe sets mapped to 12,496 gene symbols under GPL96 platform. Total 22,277 probe sets mapped to 12,495 gene symbols under GPL571 platform. Finally, 12,495 common gene symbols were identified and the gene expression matrix of specimens was received. The gene expression profile following normalization is presented in Fig. 1. The black line in the box tends towards the same straight line, indicating a good degree of standardization.

Figure 1.

Normalized expressed value data. The box in the black line means the median of each set of data, which determine the degree of standardization of data through its distribution. The green box indicates the normal sample and the red box represents the laryngeal carcinoma sample.

A total of 67 DEGs were identified in papilloma compared with the normal tissues, including 27 up-regulated genes and 40 down-regulated genes. A heat-map of the DEGs is presented in Fig. 2. The results demonstrated that the DEGs expression pattern could be used to distinguish the papillomas samples from normal controls.

Figure 2.

Heat-map overview of the differentially expressed genes. Brightness level indicates the gene expression value. When the gene expression value is higher, the color is brighter.

Functional analysis of DEGs

GO and pathway analysis indicated that up-regulated DEGs and down-regulated DEGs in papillomas samples were significantly enriched in different GO terms and KEGG pathways. A total of 17 GO terms were enriched, for example, ectoderm development, epithelium development, epithelial cell differentiation, extracellular space, and intermediate filament cytoskeleton (Table I). In addition, 2 KEGG pathways were significantly enriched, including bladder cancer and glycolysis/gluconeogenesis.

Table I.

The enriched gene ontology terms of differentially expressed genes.

Category	Term	P-value	Genes
BP	Ectoderm development	1.21E-04	PTHLH, KRT6A, KRT6B, KRT16, FLG, KRT14, AHNAK2
BP	Epithelium development	0.001897	PTHLH, FLG, VEGFA, KRT14, AHNAK2, CA2
BP	Epithelial cell differentiation	0.002007	PTHLH, FLG, VEGFA, KRT14, AHNAK2
BP	Epidermis development	0.005772	PTHLH, KRT16, FLG, KRT14, AHNAK2
BP	Regulation of cell proliferation	0.010959	PTHLH, CDKN1A, KRT6A, ID2, VEGFA, CLU, TGFBR3, FGFBP1, DPT
BP	Positive regulation of cell proliferation	0.037925	LDHA, ADH1C, FADS2, NQO1, ALDH3A1, CYP4B1, FTL
BP	Oxidation reduction	0.037925	LDHA, ADH1C, FADS2, NQO1, ALDH3A1, CYP4B1, FTL
BP	Negative regulation of molecular function	0.041489	PTHLH, CDKN1A, ID2, ID1, NQO1
CC	Extracellular region part	2.57E-09	PLAT, PZP, COL4A1, SPARCL1, CLU, LYZ, FXYD6, MMP1, PTHLH, METRN, APOL1, APOD, SERPINB5, VEGFA, NUCB2, TGFBR3, CA2, CFD, FGFBP1, DPT
CC	Extracellular space	5.20E-08	PLAT, PZP, CLU, LYZ, FXYD6, PTHLH, METRN, APOL1, APOD, SERPINB5, VEGFA, NUCB2, TGFBR3, CA2, CFD, FGFBP1
CC	Extracellular region	2.02E-04	PLAT, PZP, COL4A1, SPARCL1, CLU, LYZ, FXYD6, MMP1, PTHLH, METRN, APOL1, APOD, SERPINB5, VEGFA, NUCB2, TGFBR3, CA2, CFD, FGFBP1, DPT
CC	Intermediate filament	8.58E-04	KRT6A, KRT6B, KRT16, FLG, KRT14, BFSP1
CC	Intermediate filament cytoskeleton	9.46E-04	KRT6A, KRT6B, KRT16, FLG, KRT14, BFSP1
CC	Endoplasmic reticulum	0.038809	AADAC, EBP, NUCB2, TGFBR3, FADS2, SLC27A2, HLA-DOB, ALDH3A1, CYP4B1
CC	Proteinaceous extracellular matrix	0.040784	COL4A1, SPARCL1, VEGFA, MMP1, DPT
MF	Structural constituent of cytoskeleton	2.31E-04	KRT6A, KRT6B, KRT16, KRT14, BFSP1
MF	Structural molecule activity	0.043459	KRT6A, KRT6B, COL4A1, KRT16, FLG, KRT14, BFSP1

BP, biological process, CC, cellular component, MF, molecular function.

The PPI network included 984 nodes and 1,314 interactions (Fig. 3). The results demonstrated that a total of 34 genes had a node degree of >10, such as cyclin-dependent kinase inhibitor 1A (CDKN1A) (degree=208), N-myc downstream regulated 1 (NDRG1) (degree=87), lactate dehydrogenase A (LDHA) (degree=66), and Ras association (RalGDS/AF-6) domain family member 1 (RASSF1) (degree=48). In addition, 2 modules were identified from the PPI network (Fig. 4). Module 1 had 2 DEGs, RASSF1 and inhibitor of DNA binding 1 (ID1), and another 2 genes, suppressor of variegation 3–9 homolog 1 (SUV39H1) and cell division cycle 20 (CDC20). Besides, module 2 had 2 DEGs, NDRG1 and LDHA, and another gene, ubiquitin C (UBC).

Figure 3.

Protein-protein interaction network of DEGs. The red nodes indicate DEGs and green nodes represent genes not differentially expressed. DEGs, differentially expressed genes.

Figure 4.

Two highly connected modules extracted from the protein-protein interaction network.

Disease risk sub-pathway analysis

A total of 9 disease risk sub-pathways were identified as presented in Table II, such as glycolysis/gluconeogenesis, nitrogen metabolism, tyrosine metabolism, and phenylalanine metabolism. It should be noted that phosphoglycerate kinase 1 (PGK1), carbonic anhydrase II (CA2), and carbonic anhydrase XII (CA12) which were enriched in different sub-pathways including glycolysis/gluconeogenesis and nitrogen metabolism also had a higher node degree of >10. Moreover, aldehyde dehydrogenase 3 family, member A1 (ALDH3A1) was enriched in different pathways, including tyrosine metabolism and phenylalanine metabolism.

Table II.

The statistically significant sub-pathways identified of differentially expressed genes.

Pathway	Genes	FDR
Glycolysis/Gluconeogenesis	PGK1, PKM, ALDH3A1	0.000937
Nitrogen metabolism	CA2, CA12	0.005907
Metabolism of xenobiotics by cytochrome P450	CYP4B1, ALDH3A1	0.005907
Steroid biosynthesis	EBP	0.027276
Tyrosine metabolism	ALDH3A1	0.027276
Phenylalanine metabolism	ALDH3A1	0.040641
Drug metabolism-cytochrome P450	ALDH3A1, CYP4B1	0.040641
Fatty acid metabolism	SLC27A2	0.040641
Histidine metabolism	ALDH3A1	0.042535

FDR, represents false discovery rate.

Discussion

Advanced laryngeal cancer is an aggressive disease, which has a low rate of treatment efficacy and a high rate of recurrence (29). An improved understanding of the pathogenesis of laryngeal cancer should open therapeutic possibilities. In the present study, DEGs in laryngeal papillomas compared with normal controls were analyzed. The analysis of DEGs revealed that the DEGs were enriched in different GO terms and pathways, such as epithelial cell differentiation and the glycolysis/gluconeogenesis pathway. Besides, RASSF1 and ID1 were enriched in module 1 that extracted from the PPI network of DEGs. Furthermore, PGK1, CA2 and CA12 which had higher node degrees in the PPI network were identified to be involved in different disease risk sub-pathways, including glycolysis/gluconeogenesis and nitrogen metabolism.

RASSF1 encodes a protein similar to the RAS effector proteins and may be a potential tumor suppressor (30). Agathanggelou et al (31) have previously demonstrated that RASSF1 serves an important role in cell cycle regulation and apoptosis as well as microtubule stability and RAS-association domain family 1, isoform A gene (RASSF1A) inhibits tumor growth in both in vitro and in vivo systems. Furthermore, the findings of Sun et al (32) revealed that gene silencing through promoter hypermethylation of CpG islands is linked to the inactivation of tumor suppressor genes, such as RASSF1 in head and neck cancers. By contrast, ID1 is a helix-loop-helix (HLH) protein which inhibits the DNA binding and transcriptional activation ability of HLH family proteins with which it interacts (33). Coskunpinar et al (34) quantified the expression of tumorigenesis associated genes in laryngeal cancer and observed that a set of 16 genes had altered expression including ID1. Furthermore, Wang et al (35) reported that ID1 could interact with CDC20 and RASSF1A during early mitosis and then led to enhanced CDC20 activity. In line with the previous study, the present study identified that RASSF1A which was a hub protein in the PPI network could interact with ID1 in module 1 (Fig. 4), indicating that RASSF1A may interact with ID1 influencing cell cycle regulation and inhibiting tumor growth in laryngeal cancer.

In addition, the results of the present study showed that PGK1 were involved in the risk sub-pathway of glycolysis/gluconeogenesis. Glycolysis has been shown to be enhanced in almost all cancers, termed the ‘Warburg effect’ (36). Gatenby et al (37) demonstrated that tumor cells might sustain high glycolysis not only for biosynthetic precursor production and energy, but also for the eradication of adjacent normal cells within organ parenchyma. Besides, PGK1 is a glycolytic enzyme and may also act as a polymerase alpha cofactor (primer recognition protein) (38). Ning et al (39) reported that PGK1 in the glycolysis/gluconeogenesis pathway was significantly differentially expressed between laser micro-dissected malignant compared with benign clinical samples of prostate tissue. Moreover, Sun et al (32) observed that PGK1 were up-regulated in the tumor xenograft that originated from transketolase-like 1 (TKTL1) -expressing O11 cells which was from head and neck squamous cell carcinoma. Thus, PGK1 may serve an essential role in the glycolysis/gluconeogenesis pathway in promoting laryngeal cancer progression.

Furthermore, these findings revealed that CA2 and CA12 which had higher node degrees in PPI network were associated with nitrogen metabolism sub-pathway. Glutamine is the major source of nitrogen metabolism for nucleotide and amino acid synthesis, although many cells can metabolize glutamine in excess of their nitrogen requirement (40). Tennant et al (41)demonstrated that glutamine was an essential nutrient for many tumor cells and glutamine depletion may contribute to the effectiveness of the drug in acute lymphoblastic leukaemia. By contrast, CA2 and CA12 are two of several isozymes of carbonic anhydrase which can catalyze reversible hydration of carbon dioxide (42). Dasgupta et al (43) discovered that CA2 overexpressed in head and neck cancer cell line. In addition, Parkkila et al (44) showed that transmembrane CA12 were highly expressed in certain tumors, and may be implicated in the acidification of the extracellular enviroment surrounding cancer cells, thus creating a microenvironment beneficial to tumor growth and spread. Besides, Lascorz et al (45) found that CA2 and CA12 were enriched in the KEGG pathway nitrogen metabolism and the GO category carbonate dehydratase activity in microarray studies on colorectal carcinogenesis. In this context, CA2 and CA12 may serve crucial roles in nitrogen metabolism beneficial to laryngeal cancer progression.

In summary, the present study identified several key genes (RASSF1, PGK1, CA2 and CA12) that participated in different pathways (glycolysis/gluconeogenesis, and nitrogen metabolism) were involved in the mechanism of laryngeal cancer. RASSF1A may inhibit tumor growth in laryngeal cancer via interacting with ID1 involving in the cell cycle regulation. Besides, PGK1 may serve an essential role in the glycolysis/gluconeogenesis pathway in promoting laryngeal cancer progression. In addition, CA2 and CA12 may be crucial for nitrogen metabolism beneficial to laryngeal cancer progression. These efforts may provide an improved understanding of molecular mechanisms underlying the carcinogenesis of laryngeal cancer. Further investigation of the genes identified in this study and prospective studies to determine their potential clinical value are needed.