This article is mentioned in:

Introduction

Nasopharyngeal carcinoma (NPC) is a rare malignancy among most populations worldwide; based on previous epidemiological data, global morbidity and mortality of the disease are relatively low, both <2 per 100,000 persons (1). However, NPC is an endemic disease that is prevalent in certain regions, such as southern China, Southeast Asia, Middle East and North Africa (2). For example, in several regions of southern China, the prevalence is as high as 8–27 per 100,000 persons (3). Etiology of NPC is complicated, and family history, age and Epstein-Barr virus (EBV) infection are risk factors for disease development (4); plasma EBV DNA has been recognized as a biomarker for the prognosis of NPC (5). Genetic factors also contribute to NPC progression, such as overexpression of CD109 (6) and far upstream element-binding protein 1 (7). MicroRNAs (miRNAs) regulate genes at the transcriptional level and serve important roles during NPC progression. miRNA (miR)-9 was previously predicted to be a prognostic factor for NPC, and decreased expression has been associated with advanced tumor-node-metastasis stage (8).

However, these findings are insufficient to comprehensively reveal the pathogenesis of NPC. Epigenetic alterations also serve important roles in disease development. DNA hypermethylation is the most common epigenetic alteration; it involves the binding of a methyl group to CpG dinucleotides within a gene promoter, and is considered a hallmark in various cancer types (9). The silencing of a tumor suppressor gene may be caused by abnormal hypermethylation at the promoter. In NPC, numerous tumor suppressor genes were previously identified to be methylated (9); for example, expression of the tumor suppressor gene DLC1 Rho GTPase-activating protein was revealed to be repressed by promoter hypermethylation in NPC (10). Previous in vitro experiments indicated that overexpression of cytochrome b5 reductase 2 (CYB5R2) had an inhibitory effect on NPC; however, CYB5R2 is often inactive in human NPC owing to promoter hypermethylation (11). Other tumor suppressor genes, such as Ras association domain family member 1, cyclin dependent kinase inhibitor 2A, deleted in lung and esophageal cancer 1, death associated protein kinase 1 and in particular, ubiquitin C-terminal hydrolase L1, are also methylated in NPC patients with different incidences (9). These data suggested that there is a requirement for identifying methylation status as a prognostic biomarker in NPC.

A previous study using whole-genome methylation analysis comprised a panel of six hypermethylated genes that were identified as being correlated with poor survival in patients with NPC (12). Another study compared methylome data from methylation chip array analysis with data from The Cancer Genome Atlas Database and demonstrated that hypermethylation was a predominant pattern in NPC, and the methylation was associated with CpG islands in tumor tissues (13). Nevertheless, there is an abundance of useful data that remains unexploited in the aforementioned methylated data sets, which may contribute to our comprehensive understanding of the correlations between methylation and NPC etiology.

Therefore, the present study used the two methylation data sets to find the overlapping methylated regions in NPC. The identified genes were further investigated, as well as the potential correlations with each other or with miRNAs. Based on these comprehensive bioinformatics analyses, the results from the present study were expected to provide novel insights into the epigenetic alteration mechanisms in NPC and to identify additional biomarkers for disease prognosis.

Materials and methods

Data resource

Two data sets (GSE52068 and GSE62336) relating to NPC methylation were downloaded from the Gene Expression Omnibus database (http://www.ncbi.nlm.nih.gov/geo); both data sets were obtained from the Illumina HumanMethylation450 BeadChip platform (Illumina, Inc., San Diego, CA, USA). There were a total of 48 samples in the GSE52068 data set, which included 24 NPC cases and 24 controls (12), and 50 samples in the GSE62336 data set, which included 25 NPC cases and 25 controls (13).

Pretreatment of raw data

Signal intensity files within the two data sets were downloaded, which provided methylated and unmethylated signal values of the CpG site probes in each sample. Subsequently, the β-value that indicated the methylation status of a gene was calculated using the following formula:

The calculated value was pretreated with the β-mixture quantile normalization method as previously described (14).

Prediction of differentially methylated regions (DMRs)

The COHCAP package, which is a part of the R and Bioconductor online software suite (http://www.bioconductor.org/packages/3.0/bioc), was used to predict DMRs between NPC cases and controls in the two data sets (15). Hypermethylation and hypomethylation status were selected based on the recommendation of Sproul et al (16), which stated that if the CpG site β-value is greater than a certain value, it is defined as hypermethylation, and if the CpG site β-value is less than a certain value, it is defined as hypomethylation (16). Using COHCAP to compare NPC and control, the Δβ-values, P-values and false discovery rates (FDRs) were obtained according to the files containing β values of the CpG site probe. DMRs between the two samples were identified using the criteria of FDR <0.05 and |Δβ| >0.1 (Δβ value >0.1 indicated hypermethylation, and Δβ <-0.1 indicated hypomethylation). Subsequently, the corresponding gene symbols of the DMRs were established based on the annotation files in the methylation microarray data sets. Conversely, DEMs that did not correspond to a gene symbol and DMRs that corresponded to >1 gene symbol were filtered out.

Differentially methylated CpG islands (DMCs)

In addition to predicting the methylated site, COHCAP was also used to predict CpG island methylation, which was established by calculating the abnormal differentially methylated site counts in each gene and the average β-values. The ‘COHCAP.avg.by.island’ function implemented in COHCAP package was used to identify DMCs between NPC and normal samples in each data set, with the parameters: ‘num.sites=4, which reflected that the minimal number of differentially methylated sites was 4; |Δβ| >0.1; and FDR <0.05.

Common DMRs and DMCs in two data sets

Comparing the DMR and DMC information from the two databases, overlapping DMRs and DMCs were identified and extracted. If a DMR appeared in both data sets, it suggested that the corresponding genes may be highly associated with NPC. From these resulting data the corresponding genes of the overlapping DMRs were obtained.

Enrichment analysis of the overlapped genes

The Database for Annotation, Visualization and Integration Discovery (version 6.8; http://david.abcc.Ncifcrf.gov) software was used to perform Gene Ontology (GO; http://www.geneontology.org) term analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG; http://www.genome.jp/kegg/pathway.html) pathway enrichment analysis, under the conditions of gene count ≥2 and P<0.05 based on hypergeometric test, as previously described (17).

Construction of protein-protein interaction (PPI) network of overlapping DMR genes

The identified DMR genes were mapped using the Search Tool for the Retrieval of Interacting Genes (http://string-db.org) database, version 10.0 (18), to explore their potential interactions at protein level. A PPI network for these genes was constructed using the parameters: species, Homo sapiens (hsa); and PPI score=0.4. In addition, it was required that these PPIs were derived from validated experiments. The PPI network was visualized with Cytoscape (http://cytoscape.org) software, and the topological property was analyzed based on the degree of a node. For nodes in the PPI network, the random walk algorithm was used to analyze them to obtain marker genes of NPC, as previously described (19).

Identification of miRNAs associated with DMRs and their targets

Overlapping DMR genes in the two data sets were further examined to determine whether they encoded miRNAs. The putative miRNA-coding genes were combined with target gene information in the miRWalk database to predict their target genes using the ‘validated target module’ (20). The target genes were compared with DMRs to select the target genes connecting with other DMR genes.

NPC related gene analysis

The Comparative Toxicogenomics Database (CTD) is a comprehensive database that records disease-related genes that have been previously reported (21). Based on information in this database, NPC related genes were downloaded to determine whether they were methylated.

Results

DMRs between NPC samples and controls in two data sets

Using the thresholds of FDR <0.05 and |Δβ| >0.1, a total of 4,218 hypermethylated and 2,178 hypomethylated CpG sites were identified in the GSE52068 data set, which corresponded to 2,243 and 2,013 genes, respectively. In the GSE62336 data set, 11,347 hypermethylated and 39 hypomethylated CpG sites were identified, which corresponded to 4,372 and 37 genes, respectively.

Overlapped DMR related genes and their enriched functions and pathways

The identified genes from the two data sets were compared, and the overlapping DMRs were selected. A total of 3,306 shared CpG sites were identified as hypermethylated, and corresponded to 1,854 genes; and 18 shared CpG sites were identified as hypomethylated, corresponding to 18 genes.

According to KEGG pathway enrichment and GO term analyses (Fig. 1A and B, respectively), the selected genes were revealed to be highly associated with pathways including mitogen-activated protein kinase (MAPK) signaling pathway, phosphatidylinositol 3-kinase (PI3K)/AKT signaling pathway, calcium signaling pathway, pathways in cancer and neuroactive ligand receptor interaction. With regards to the functions, these genes were significantly enriched in biological process including positive regulation of transcription from RNA polymerase II promoter, transcription from RNA polymerase II promoter and multicellular organism development. Enriched cellular component annotations included integral component of plasma membrane, plasma membrane and cell junction. Molecular function annotations included the following: transcription factor activity, sequence-specific DNA binding;, transcription factor binding; sequence-specific DNA binding.

Figure 1. Enriched GO functions and KEGG pathways of overlapped genes with differentially methylated regions. (A) KEGG pathway analysis. (B) GO term function category. BP, biological process; CC, cellular component; GO, gene ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; MF, molecular function.

PPI network of the overlapped DMR corresponding genes

Based on the aforementioned criteria, including PPI score=0.4, a PPI network was established that contained 677 nodes and 991 interactions. The random walk algorithm identified 20 nodes that were predominant with high degree (Fig. 2), such as proto-oncogene tyrosine-protein kinase SRC (degree=50), SMAD family member 3 (SMAD3; degree=31), tyrosine 3-monooxygenase/tryptophan 5-monooxygenase activation protein ζ (YWHAZ; degree=28), β-actin (ACTB; degree=27) and Heat shock protein family A member 4 (HSPA4; degree=24). Notably, these genes were all hypomethylated.

Figure 2. Protein-protein interaction network of the overlapped genes with differentially methylated regions. Orange nodes indicate hypomethylation; blue nodes indicate hypermethylation. ‘*’ Indicates the following genes: SRC, SMAD3, YWHAZ, ACTB and HSPA4. SRC, proto-oncogene tyrosine-protein kinase; SMAD3, SMAD family member 3; YWHAZ, tyrosine 3-monooxygenase/tryptophan 5-monooxygenase activation protein ζ; ACTB, β-actin; HSPA4, Heat shock protein family A member 4.

Abnormally methylated miRNAs among the overlapped DMR genes

Among the overlapped DMR genes, 14 were identified as miRNA-encoding genes, which regulated 306 other DMR genes, such as hsa-miR-148a-3p, hsa-miR-129-5p, hsa-miR-150-5p, hsa-miR-124-3p and hsa-miR-196a-5p (Fig. 3). Based on the prediction of miRNA-target relationships, HSPA4 was predicted as a target of hsa-miR-148a-3p, and ACTB was the target of hsa-miR-124-3p and hsa-miR-196a-5p.

Figure 3. miRNA-target network of methylated miRNAs between the overlapped genes with differentially methylated regions. Circles represent target genes and triangle stands for miRNAs. miRNA, microRNA.

DMC analysis between two different samples

Using the aforementioned parameters of num.sites=4, |Δβ| >0.1 and FDR <0.05, a total of 195 DMCs were identified in the GSE52068 data set, and 798 DMCs were identified in the GSE62336 data set. Between them, 150 were overlapped and hypermethylated, and 47 of the overlapped genes were related to NPC, which were defined as NPC genes; a heat map of the 47 identified NPC genes is presented in Fig. 4. Notably, an miRNA-encoding gene was also identified, miR129-2 chr11:43602545-43603215.

Figure 4. Heat map of the differentially methylated CpG islands between NPC samples and normal samples in the (A) GSE52068 and (B) GSE62336 data sets. NPC, nasopharyngeal carcinoma.

Identification of NPC related genes

Based on information obtained from the CTD database, it was revealed that 583 NPC genes had been previously reported to associate with hypermethylated DMRs and three NPC genes were related to hypomethylated DMRs. Notably, 17 of them were highlighted in the PPI network, such as cell division cycle 42 (CDC42), α1 actin (ACTA1), paired box 6 (PAX6), SMAD3, matrix metallopeptidase 2 (MMP2), collagen type I α1 (COL1A1), protein kinase Cγ (PRKCG), SRC, phospholipase Cγ (PLCG1), ubiquitin-conjugating enzyme E2 G2 (UBE2G2), synuclein α (SNCA), HSPA4, sequestosome 1 (SQSTM1), ACTB, tyrosine 3-monooxygenase/tryptophan 5-monooxygenase activation protein γ (YWHAG), fibroblast growth factor receptor 1 (FGFR1) and YWHAZ. Among them, MMP2 was determined to be an NPC marker gene.

Discussion

The present study identified a number of crucial genes that were indicted to be hypomethylated, instead of hypermethylated, in the NPC samples, including SRC, SMAD3, YWHAZ and HSPA4. These genes were predominant in the PPI network. Among them, HSPA4 was predicted as a target of hsa-miR-148a-3p. SRC and SMAD3 were correlated with NPC based on the CTD database. Most of the identified DMCs were revealed to be hypermethylated, and 47 were associated with NPC. A miRNA-encoding gene was also highlighted, miR129-2 chr11:43602545-43603215.

Promoter methylation is implicated in the regulation of gene expression, and hypomethylation is commonly exhibited in several cancer types, such as melanoma and squamous cell carcinoma (SCC) (22). Hypomethylation of deleted in split hand/split foot 1 was demonstrated to result in the overexpression of this gene, which was associated with poor prognosis in patients with SCC (22). SRC encodes a tyrosine-protein kinase protein; it is an oncogene that has a similar function as the v-src gene of the Rous sarcoma virus. Activation of cellular SRC (c-SRC) was previously reported to promote metastasis of NPC by regulating the PI3K/AKT signaling pathway (23). The long non-coding RNA actin filament associated protein 1-antisense RNA 1 (AFAP1-AS1) was demonstrated to act as an adaptor that connects other proteins like SRC in the mediation of actin filament integrity (24). In NPC, upregulation of AFAP1-AS1 expression may lead to tumor metastasis and has been associated with poor prognosis (24). These data suggested that overexpression of SRC may be associated with metastasis and poor prognosis of NPC. In the present study, SRC was demonstrated to be hypomethylated and was implicated in NPC. On this basis, it was inferred that hypomethylation of SRC may result in the upregulation of this gene, which may contributes to the NPC metastasis and poor prognosis.

SMAD3 is a member of SMAD family of signal transducers and transcriptional regulators. Activated by transforming growth factor-β (TGF-β), SMAD3 serves a crucial role in promoting carcinogenesis (25). A previous study on TGF-β/Smad signaling in NPC reported that the expression of SMAD3 was significantly increased in the CNE2 NPC cell line (26). In addition, downregulation of SMAD3 expression by miR-145 was demonstrated to inhibit the invasion and metastasis of NPC (27), which suggested that the elevated SMAD3 expression level may account for metastasis of NPC. Other previous studies have reported the relationship between SMAD3 expression and methylation, and it was demonstrated that different SMAD3 expression levels had similar promoter methylation patterns in the maintenance of self-tolerance (28). However, unlike these autoimmune-associated lesions, results from the present study indicated that SMAD3 was hypomethylated in NPC, which suggested that SMAD3 hypomethylation may result in the upregulated expression of this gene and may account for NPC metastasis. These results indicated that the upregulation of SMAD3 by hypomethylation may be linked to poor prognosis of NPC.

HSPA4 belongs to the HSP70 protein family. Induced by Nibrin, an important DNA repair protein that maintains genomic integrity in the advanced stages NPC, HSPA4 was reported to facilitate the metastatic activity in tumor cells (29). In addition, polymorphism of another HSP70 family member, HSP70-2, was demonstrated to be involved with susceptibility to NPC (30). These data suggested that abnormal expression of HSP70-related genes may be highly associated with NPC progression. However, a connection between HSPA4 expression and its methylation status have not been reported. In the present study, HSPA4 was indicated to be hypomethylated in NPC samples, which may account for an upregulation in expression. These results indicated that HSPA4 may be important to NPC metastasis and hypomethylation may be a predictive factor for disease prognosis.

In terms of miRNAs, miR-148a was previously reported to be downregulated in NPC (31), and it may be involved in epithelial-mesenchymal transition (EMT) and metastasis in hepatocellular carcinoma (32); however, this function has not been reported in NPC. The present study predicted that HSPA4 was a target of hsa-miR-148a-3p. Combining this with the hypomethylation of HSPA4, it may be inferred that hsa-miR-148a-3p binds to the gene promoter region that controls methylation of this gene during NPC metastasis. However, this hypothesis requires further analysis through dual-luciferase reporter system.

By binding to the phosphoserine-containing proteins, the YWHAZ protein serves an important role in the regulation of signal transduction (33). In NPC cells, EBV-miR-BART10-3p was reported to induce the expression of β-catenin by inhibiting the expression of β-transducin repeat-containing E3 ubiquitin protein ligase (BTRC), which facilitates EMT and promotes metastasis (34). Formation of the β-catenin/YWHAZ complex was reported to inhibit the binding of β-catenin to BTRC, which ensures the stability of β-catenin and thus promotes EMT and metastasis in lung cancer (35). Neither the function nor the methylation status of YWHAZ in NPC has been reported; however, hypermethylation in the gene promoter was previously demonstrated to suppress the expression of this gene. The hypomethylation of YWHAZ reported in the present study may indicate an increase in gene expression. Therefore, it was speculated that, similar to lung cancer, the increase in YWHAZ expression may also lead to the EMT in NPC and promote metastasis. This hypothesis, however, needs to be validated in future studies.

miR129-2 was previously demonstrated to act as a tumor suppressor, and the methylation of miR129-2 was reported to be a frequent event in lymphoid malignancies (36). Differential methylation of this miRNA has been indicated in rectal cancer and colorectal cancer (37,38). In the present study, miR129-2 was revealed to be hypermethylated in NPC, which suggested that expression may be inhibited by hypermethylation in NPC development.

It should be noted that the inferred effects of the identified methylated genes and miRNAs, as well as their predictive correlations in NPC, require further experimental validation. Despite this limitation, the present study may be of great value for providing novel methylation biomarkers for NPC detection and prognosis.

In conclusion, several novel methylated genes and miRNA were identified that may serve as potential biomarkers for NPC prognosis. Hypomethylation of SRC, SMAD3, YWHAZ and HSPA4, and the hypermethylation of miR129-2, may be linked with poor prognosis of NPC. Nevertheless, all the predicted results need to be further validated by experiments.

Glossary

Abbreviations

Abbreviations:

References