Identification of genes involved in Epstein-Barr virus-associated nasopharyngeal carcinoma.

Nasopharyngeal carcinoma (NPC) is the most common cancer originating from the nasopharynx, and can be induced by infection with Epstein-Barr virus (EBV). To study the mechanisms of EBV-associated NPC, a microarray of the GSE12452 dataset was analyzed. GSE12452 was downloaded from Gene Expression Omnibus and consisted of 31 NPC samples and 10 normal healthy nasopharyngeal tissue samples. The differentially-expressed genes (DEGs) were screened using the linear models for microarray data package in R. Using Database for Annotation, Visualization and Integrated Discovery software, potential functions of the DEGs were predicted by Gene Ontology and pathway enrichment analyses. With the information from the Search Tool for the Retrieval of Interacting Genes/Proteins database, the protein-protein interaction (PPI) network was visualized by Cytoscape. Furthermore, modules of the PPI network were searched using ClusterONE in Cytoscape. A total of 951 DEGs were screened in the NPC samples compared with the normal healthy nasopharyngeal tissue samples. Function enrichment indicated that the upregulated genes were associated with the cell cycle, cytoskeleton organization and DNA metabolism. Meanwhile, the downregulated genes were mainly associated with cell differentiation, hormone metabolism, inflammatory response and immune response. PPI networks for the DEGs suggested that upregulated mitotic arrest deficient 2-like 1 (MAD2L1; degree=133), proliferating cell nuclear antigen (PCNA; degree=125) and cyclin B1 (CCNB1; degree=115), and downregulated member A1 of aldehyde dehydrogenase 1 (ALDH1A1; degree=15) may be of great importance as they exhibited higher degrees on interaction. Mucin 1 (MUC1) was a key node of module 4. Overall, the study indicated that MAD2L1, CCNB1, PCNA, ALDH1A1 and MUC1 may have a correlation with EBV-associated NPC.

Introduction

As the most common cancer originating from the nasopharynx (1), nasopharyngeal carcinoma (NPC) can be caused by viral influence, heredity and environmental factors (2). The viral influence is correlated with infection by Epstein-Barr virus (EBV), which is a B-lymphotropic herpes virus possessing growth-transforming properties (3). It has been reported that 95% of Americans are exposed to this virus in their thirties (4). In 2010, NPC led to 65,000 mortalities globally (5). Thus, there is an urgent requirement to study the mechanisms of EBV-associated NPC.

Recently, a number of studies have been performed to investigate the mechanisms of EBV-associated NPC. For example, as an early EBV antigen (6), BamH1-A Reading Frame-1 may be associated with the pathogenesis of NPC, such as the malignant transformation of human NPC epithelial cells (7,8). Aberrant hypermethylation of Ras association domain family 1 isoform A and the high viral load of EBV DNA may play an important role in NPC pathogenesis, thus, they may function as promising diagnostic markers for NPC (9,10). Serological results have shown that specifically expressed BRLF1 may be used in the diagnosis of NPC (11). Via the activation of Ets-1, c-Met can be induced by latent membrane protein-1 (LMP-1) and enhance the highly metastatic potential of NPC (12). By providing epitopes recognized by cytotoxic T-cells, EVB-encoded LMP2A may be a potential target and may be used in the immunotherapy of NPC (13).

In 2006, Sengupta et al (14) analyzed the expression of all latent EBV genes between NPC samples and normal healthy nasopharyngeal epithelium samples, and obtained a panel of differentially-expressed genes (DEGs). Using the same data by Sengupta et al (14), the present study aimed to further screen the DEGs and predict their underlying function by functional and pathway enrichment analyses. Furthermore, protein-protein interaction network (PPI) networks were constructed and modules of PPI network were searched to investigate the interaction associations between these DEGs.

Materials and methods

Microarray data

The expression profile of the GSE12452 dataset deposited by Sengupta et al (14) was downloaded from Gene Expression Omnibus (http://www.ncbi.nlm.nih.gov/geo/), which was based on the platform of the GPL570 [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array. GSE12452 consisted of a collection of 31 NPC samples and 10 normal healthy nasopharyngeal tissue samples. These samples were collected from Taiwanese patients who provided informed consent. After the tissues were resected and immediately flash frozen, samples were finally stored in liquid nitrogen.

DEG screening

Once GSE12452 had been downloaded, the microarray data was read using Affy package (www.bioconductor.org) (15) and Refseq annotation files, and then it was normalized by the Robust MultiArray Averaging method (16). The linear models for microarray data package (http://www.bioconductor.org) (17) in R was used to identify the DEGs between NPC samples and normal healthy nasopharyngeal tissue samples. The adjusted P-value of <0.01 and |logfold-change (FC)|>1 were used as the cut-off criteria.

Functional and pathway enrichment analysis

Gene Ontology (GO) terms can describe three wild-type gene products, including molecular function, biological process and subcellular location (18). The Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway database is used for the systematic analysis of gene functions, which can connect genomic information with functional information (19). The Database for Annotation, Visualization, and Integrated Discovery (DAVID) is a program that can integrate functional genomic annotations with intuitive summaries (20). Using DAVID software, GO and KEGG pathway enrichment analyses were performed for DEGs between NPC samples and normal healthy nasopharyngeal tissue samples. A P-value of <0.05 was used as the cut-off criterion.

PPI network and module construction

Interaction associations of the proteins encoded by the DEGs were searched by Search Tool for the Retrieval of Interacting Genes/Proteins online software (http://string-db.org/) (21), and then PPI networks were visualized by Cytoscape software (http://www.cytoscape.org) (22). The proteins in the PPI network were termed nodes and the degree of a node correlated with the number of its interactions. Using network statistics, a connectivity degree analysis was conducted for each node in the PPI network. ClusterONE (http://www.paccanarolab.org/clusterone/) (23) in Cytoscape was used to screen modules of the PPI network. A P-value of <0.01 was used as the cut-off criterion.

Results

DEG analysis

Compared with normal healthy nasopharyngeal tissue samples, a total of 951 DEGs were screened from the NPC samples, including 376 upregulated and 575 downregulated genes. There were more downregulated genes than upregulated genes.

The enriched GO functions for upregulated genes were divided into 28 clusters. Functions in cluster 1 were associated with the cell cycle, such as the M phase (P=5.27×10−30) and the cell cycle process (P=3.34×10−29). Functions in cluster 2 were associated with cytoskeleton organization, such as the microtubule-based process (P=1.09×10−11) and microtubule cytoskeleton organization (P=7.27×10−11). Functions in cluster 3 were implicated in the regulation of the cell cycle, such as the regulation of the mitotic cell cycle (P=6.43×10−9) and the regulation of cell cycle process (P=7.79×10−6). Functions in cluster 4 were involved in DNA metabolism, such as the DNA metabolic process (P=1.80×10−12) and the response to a DNA damage stimulus (P=5.83×10−8) (Table I).

Table I.

Enriched GO functions and KEGG pathways for the upregulated and downregulated genes[a].

A, Enriched GO functions for the upregulated genes
Category	Term		Description	Gene no.	Gene symbol[b]	P-value
BP	GO:0000279		M phase	53	DBF4, KNTC1, TTK	5.27×10−30
BP	GO:0022402		Cell cycle process	66	DBF4, KNTC1, TTK	3.34×10−29
BP	GO:0007017		Microtubule-based process	28	KIF23, CAV1, KIF4A	1.09×10−11
BP	GO:0000226		Microtubule cytoskeleton organization	21	KIF23, CAV1, KIF11	7.27×10−11
BP	GO:0007346		Regulation of mitotic cell cycle	19	CDC7, DLGAP5, TIPIN	6.43×10−9
BP	GO:0010564		Regulation of cell cycle process	13	CDC7, DLGAP5, TIPIN	7.79×10−6
BP	GO:0006259		DNA metabolic process	41	UNG, DBF4, TIPIN	1.80×10−12
BP	GO:0006974		Response to DNA damage stimulus	28	UNG, TIPIN, PRKDC	5.83×10−8
B, Enriched GO functions for the downregulated genes
Category	Term		Description	Gene no.	Gene symbol[b]	P-value
BP	GO:0030855		Epithelial cell differentiation	16	ELF3, FOXA1, ANXA1	4.04×10−7
BP	GO:0030216		Keratinocyte differentiation	7	SPRR1A, PPL, CNFN	3.48×10−3
BP	GO:0010817		Regulation of hormone levels	12	FAM3B, FOXA1, DUOX2	5.75×10−4
BP	GO:0034754		Cellular hormone metabolic process	5	DHRS9, UGT1A6, UGT1A10	4.25×10−2
BP	GO:0002526		Acute inflammatory response	9	F3, CLU, SAA4	1.55×10−3
BP	GO:0009611		Response to wounding	23	S100A9, ANXA1, PRDX5	3.86×10−3
BP	GO:0006959		Humoral immune response	8	CFB, FOXJ1, CLU	1.89×10−3
BP	GO:0002455		Humoral immune response mediated by circulating immunoglobulin	4	C7, CD55, CR2	3.12×10−2
C, Enriched KEGG pathways for the upregulated and downregulated genes
Regulation	Category	Term	Description	Gene no.	Gene symbol[b]	P-value
Up	KEGG	04110	Cell cycle	21	COL4A2, COL4A1, COL3A1	3.32×10−11
	KEGG	04512	ECM-receptor interaction	14	COL4A2, COL4A1, COL3A1	1.75×10−7
	KEGG	03030	DNA replication	9	RFC5, DNA2, RFC3	2.74×10−6
	KEGG	03430	Mismatch repair	7	RFC5, EXO1, MSH6	1.92×10−5
	KEGG	04115	p53 signaling pathway	9	CCNE2, CCNB1, CDK1	3.32×10−4
Down	KEGG	00980	Metabolism of xenobiotics by cytochrome P450	8	GSTA1, GSTA3, ADH1C	1.82×10−4
	KEGG	00982	Drug metabolism	8	GSTA3, ADH1C, ADH7	2.25×10−4
	KEGG	00830	Retinol metabolism	6	UGT1A1, UGT1A7, ALDH1A1	4.45×10−3
	KEGG	00590	Arachidonic acid metabolism	5	AKR1C3, GGT6, ALOX15	2.61×10−2
	KEGG	04640	Hematopoietic cell lineage	6	MS4A1, CD1C, CD1D	2.96×10−2

The enriched GO terms listed in the table are the functions for the top four clusters with the highest enriched score.

Only the most relevant genes in each category are listed. GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; ECM, extracellular matrix; BP, biological process. KNTC, kinetochore associated; KIF, kinesin family; CAV, caveolin; CDK, cyclin-dependent kinase; CDC, cell division cycle; DLGAP, discs large homolog associated protein; TIPIN timeless-interacting protein; UNG, uracil DNA glycosylase; PRKDC, protein kinase, DNA-activated, catalytic polypeptide; SPRR, small proline rich protein; PPL, periplakin; CNFN, cornifelin; DHRS, dehydrogenase/reductase (SDR family); UGT, UDP glucuronosyltransferase; ELF3, E74-like ETS transcription factor 3; FOX, forkhead box; ANX, annexin; FAM, family with sequence similarity; DUOX, dual oxidase; F3, coagulation factor III, tissue factor; CLU, clusterin; SAA4, serum amyloid A4, constitutive; S100A9, S100 calcium binding protein A9; PRDX, peroxiredoxin; CFB, complement factor B; C7, complement component 7; CD, cluster of differentiation; CR2, complement component 3d receptor 2; COL, collagen; RFC, replication factor C; DNA2, DNA replication helicase/nuclease 2; EXO, exonuclease; MSH, mutS homolog; CCN, cyclin; GSTA, glutathione S-transferase alpha; ADH, alcohol dehydrogenase; AKR, aldo-keto reductase; GGT, gamma-glutamyltransferase; ALOX15, arachidonate 15-lipoxygenase; MS4A1, membrane spanning 4-domains A1.

The enriched GO functions for downregulated genes were divided into 9 clusters. For instance, functions in cluster 1 were associated with cell differentiation, such as epithelial cell differentiation (P=4.04×10−7) and keratinocyte differentiation (P=3.48×10−3). Functions in cluster 2 were involved in hormone metabolism, such as the regulation of hormone levels (P=5.75×10−4) and the cellular hormone metabolic process (P=4.25×10−2). Functions in cluster 3 were associated with the inflammatory response, such as the acute inflammatory response (P=1.55×10−3) and the response to wounding (P=3.86×10−3). Functions in cluster 4 were implicated in the immune response, such as the humoral immune response (P=1.89×10−3) and the humoral immune response mediated by circulating immunoglobulin (P=3.12×10−2) (Table I).

The enriched KEGG pathways for upregulated genes are also listed in Table I, including the cell cycle (P=3.32×10−11), extracellular matrix-receptor interaction (P=1.75×10−7) and DNA replication (P=2.74×10−6). Additionally, pathway enrichment analysis was conducted for downregulated genes. As shown in Table I, pathways such as those for the metabolism of xenobiotics by cytochrome P450 (P=1.82×10−4), drug metabolism (P=2.25×10−4) and retinol metabolism (P=4.45×10−3) were enriched (Table I).

PPI network and module analysis

The PPI network for upregulated genes had 279 nodes and 5,690 interactions. In particular, upregulated mitotic arrest deficient 2-like 1 (MA D2L1; degree=133), replication factor C4 (degree=130), proliferating cell nuclear antigen (PCNA; degree=125), cyclin-dependent kinase 1 (degree=124) and cyclin B1 (CCNB1; degree=115) exhibited higher degrees in the PPI network. Modules 1 (Fig. 1) and 2 (Fig. 2) were obtained from the PPI network of upregulated genes. Module 1 had 144 nodes and 5,091 interactions. Module 2 had 23 nodes and 133 interactions. The cell cycle was an enriched pathway for the DEGs in modules 1 and 2.

Figure 1.

Module 1 obtained from the protein-protein interaction network of upregulated genes.

Figure 2.

Module 2 obtained from the protein-protein interaction network of upregulated genes.

The PPI network for downregulated genes had 260 nodes and 517 interactions. Importantly, downregulated glutathione S-transferase α1 (GSTA1; degree=19), GSTA3 (degree=19), dynein, light chain, roadblock-type 2 (degree=17), calmodulin 1 (degree=16), member A1 of aldehyde dehydrogenase 1 (ALDH1A1; degree=15) and sperm-associated antigen 6 (degree=15) had higher degrees in the PPI network. Modules 3 (Fig. 3) and 4 (Fig. 4) were obtained from the PPI network of downregulated genes. Module 3 had 13 nodes and 48 interactions. The enriched pathway for DEGs in module 3 included the mitogen-activated protein kinase signaling pathway. Module 4 had 11 nodes, such as mucin 1 (MUC1), and 43 interactions. O-glycan biosynthesis was an enriched pathway for DEGs in module 4.

Figure 3.

Module 3 obtained from the protein-protein interaction network of downregulated genes.

Figure 4.

Module 4 obtained from the protein-protein interaction network of downregulated genes.

Discussion

In the present study, 951 DEGs were screened in the NPC samples compared with normal healthy nasopharyngeal tissue samples. Function enrichment indicated that the upregulated genes were mainly associated with the cell cycle, cytoskeleton organization and DNA metabolism. The downregulated genes were associated with cell differentiation, hormone metabolism, the inflammatory response and the immune response. Meanwhile, upregulated MAD2L1 (degree=133), PCNA (degree=125) and CCNB1 (degree=115), and downregulated ALDH1A1 (degree=15) exhibited higher degrees of interaction in the PPI network.

Via the induction of mitotic arrest, MAD2 can cause chemosensitization to cisplatin in NPC cells and activate the apoptosis pathway (24). The aberrantly reduced expression of MAD2 can lead to a defective mitotic checkpoint and promote chromosomal instability in the disease (25). Sensitization to vincristine induced by MAD2 is associated with Raf/Bcl-2 phosphorylation and mitotic arrest in NPC cells (26). The M-phase events of MAD2 and CCNB1/CDC2 activation are essential to the paclitaxel-induced apoptosis of human NPC cells (27). Thus, the expression levels of MAD2L1 and CCNB1 may be associated with NPC. It has been reported that upregulated BCL-2 and a high PCNA labeling index may be implicated in local recurrence in NPC patients receiving the primary treatment of radiation therapy (28). Through inhibiting the expression of PCNA, a PCNA-small interfering RNA compound can effectively interfere with NPC cells and may be used in the gene therapy of NPC (29). These data indicate that PCNA may also play a role in NPC.

The expression of ALDH1 in the invasive front links, which is associated with tumor aggressiveness and epithelial-mesenchymal transition characteristics, may serve as a promising predictive factor in NPC (30,31). Budding cells characterized by a high level of ALDH1 may have invasive and metastatic properties, therefore, the degree of expression can be a useful prognostic marker in NPC patients (32). It has been reported that the overexpression of ALDH1A1 in NPC is associated with enhanced invasiveness (33). These data indicate that ALDH1A1 may have a close correlation with NPC. MUC1, a mucinous glycoprotein required for the detachment and release of tumor cells, combined with other invasiveness and angiogenic factors may function in a complex sequential process that ends in the metastasis of EBV-infected NPC cells (34). This may indicate that the expression level of MUC1 is associated with EBV-associated NPC.

In conclusion, in the present study, an integrated bioinformatics analysis of genes that may be involved in EBV-associated NPC was performed. A total of 951 DEGs were screened in the NPC samples compared with the normal healthy nasopharyngeal tissue samples. Certain DEGs, such as MAD2L1, CCNB1, PCNA, ALDH1A1 and MUC1, may have a correlation with NPC. However, future studies are required to advance the understanding of their mechanisms of action in NPC.